矩阵中QR算法分解简介和基于Eigen库使用示例

发布于:2025-07-25 ⋅ 阅读:(17) ⋅ 点赞:(0)

QR 算法是一种用于**特征值分解(Eigen Decomposition)**的迭代数值方法。广泛应用于求解实对称矩阵的特征值与特征向量,是现代数值线性代数中核心算法之一。


一、定义与目标

**QR算法(QR Algorithm)**是一种迭代方法,目标是对矩阵 $A$ 求其特征值:

  • 输入:A∈Rn×nA \in \mathbb{R}^{n \times n}ARn×n(一般为实对称矩阵)
  • 输出:AAA 的特征值和特征向量

该算法通过将矩阵分解为:

Ak=QkRk A_k = Q_k R_k Ak=QkRk

然后重新组合得到:

Ak+1=RkQk A_{k+1} = R_k Q_k Ak+1=RkQk

重复这个过程,多次迭代后 A_kA\_kA_k 收敛为一个上三角矩阵,其对角元素即为特征值。


二、原理与推导

假设 $A$ 可对角化:

A=VΛV−1 A = V \Lambda V^{-1} A=VΛV1

若我们不断对 $A_k$ 做 QR 分解:

Ak=QkRk⇒Ak+1=RkQk A_k = Q_k R_k \Rightarrow A_{k+1} = R_k Q_k Ak=QkRkAk+1=RkQk

则:

Ak+1=QkTAkQk A_{k+1} = Q_k^T A_k Q_k Ak+1=QkTAkQk

这是一个相似变换(Similarity Transformation),意味着 $A_{k+1}$ 与 $A_k$ 拥有相同特征值。

由于 $A_k$ 是逐步对角化的,最终趋于上三角矩阵,其对角线收敛到 $A$ 的特征值。


三、算法步骤

原始 QR 算法(无加速)

  1. 初始化:$A_0 = A$

  2. 对 $A_k$ 做 QR 分解:

    Ak=QkRk A_k = Q_k R_k Ak=QkRk

  3. 计算 $A_{k+1} = R_k Q_k$

  4. 判断是否收敛(如 $|A_{k+1} - A_k| < \epsilon$)

  5. 重复步骤 2~4

改进:带位移的 QR 算法(Shifted QR)

通过引入位移(shift)$\mu_k$ 加快收敛速度:

Ak−μkI=QkRk⇒Ak+1=RkQk+μkI A_k - \mu_k I = Q_k R_k \Rightarrow A_{k+1} = R_k Q_k + \mu_k I AkμkI=QkRkAk+1=RkQk+μkI

位移常选为 $A_k$ 的最后一行最后一列的元素(Wilkinson shift)。


四、应用场景

应用 描述
特征值计算 高效计算中小型矩阵的所有特征值
PCA QR 算法可用于求协方差矩阵的特征值(降维)
结构力学 解析振动系统特征频率
图像压缩 SVD(特征值分解的推广)相关算法
SLAM 优化 QR 用于线性系统求解,如 BA 优化中稀疏 QR

五、C++ 实现示例(使用 Eigen 库)

示例1

#include <Eigen/Dense>
#include <iostream>

using namespace Eigen;
using namespace std;

int main() {
    MatrixXd A(3, 3);
    A << 2, -1, 0,
         -1, 2, -1,
          0, -1, 2;

    MatrixXd Ak = A;
    int maxIter = 100;
    double tol = 1e-8;

    for (int k = 0; k < maxIter; ++k) {
        HouseholderQR<MatrixXd> qr(Ak);
        MatrixXd Q = qr.householderQ();
        MatrixXd R = qr.matrixQR().triangularView<Upper>();
        Ak = R * Q;

        // 检查是否近似对角化
        MatrixXd offDiag = Ak - Ak.diagonal().asDiagonal();
        if (offDiag.norm() < tol) break;
    }

    cout << "近似对角矩阵 Ak:\n" << Ak << endl;
    cout << "特征值(近似):\n" << Ak.diagonal().transpose() << endl;

    return 0;
}

示例2-Eigen::SparseQR 示例

依赖:

  • #include <Eigen/Sparse>
  • #include <Eigen/SparseQR>

示例代码:稀疏 QR 解线性系统 $Ax = b$

#include <iostream>
#include <Eigen/Sparse>
#include <Eigen/SparseQR>

int main() {
    using namespace Eigen;

    typedef SparseMatrix<double> SpMat;
    typedef Triplet<double> T;

    // 构造稀疏矩阵 A (5x5 三对角矩阵)
    std::vector<T> tripletList;
    tripletList.reserve(13);
    for (int i = 0; i < 5; ++i) {
        tripletList.emplace_back(i, i, 2.0);       // 对角
        if (i > 0) tripletList.emplace_back(i, i - 1, -1.0); // 下对角
        if (i < 4) tripletList.emplace_back(i, i + 1, -1.0); // 上对角
    }

    SpMat A(5, 5);
    A.setFromTriplets(tripletList.begin(), tripletList.end());

    // 右侧向量 b
    VectorXd b(5);
    b << 1, 2, 3, 4, 5;

    // 使用稀疏 QR 分解
    SparseQR<SpMat, COLAMDOrdering<int>> solver;
    solver.compute(A);
    if (solver.info() != Success) {
        std::cerr << "分解失败!" << std::endl;
        return -1;
    }

    // 解 Ax = b
    VectorXd x = solver.solve(b);
    if (solver.info() != Success) {
        std::cerr << "求解失败!" << std::endl;
        return -1;
    }

    std::cout << "解 x:\n" << x << std::endl;

    return 0;
}

相关注意事项

项目 说明
SparseQR 用于稀疏矩阵的列主元 QR 分解
COLAMDOrdering 使用最少填充的列重排序算法,提高分解效率
Triplet 构造 推荐方式,避免反复插入
矩阵类型 使用 Eigen::SparseMatrix<T>,不要使用 MatrixXd
求解过程 compute() 先分解,solve() 解方程组


六、注意事项与扩展

  • 对称矩阵 QR 收敛较快,非对称时需要更复杂处理(如 Hessenberg 预处理)
  • Eigen 库中已经内建 SelfAdjointEigenSolver 使用分块 + QR 实现
  • 稀疏矩阵场景中用 SparseQR 更合适
  • 在 SLAM 中,QR 分解(尤其是列主元排序 + 稀疏结构保留)常用于优化线性子问题


网站公告

今日签到

点亮在社区的每一天
去签到