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

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

---

一、定义与目标

**QR算法（QR Algorithm）**是一种迭代方法，目标是对矩阵 $A$ 求其特征值：

• 输入：A∈Rn×nA \\in \\mathbb{R}^{n \\times n}A∈Rn×n（一般为实对称矩阵）
• 输出：$A$ 的特征值和特征向量

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

Ak=QkRkA_k = Q_k R_kAk​=Qk​Rk​

然后重新组合得到：

Ak+1=RkQkA_{k+1} = R_k Q_kAk+1​=Rk​Qk​

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

---

二、原理与推导

假设 $A$ 可对角化：

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

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

Ak=QkRk⇒Ak+1=RkQkA_k = Q_k R_k \\Rightarrow A_{k+1} = R_k Q_kAk​=Qk​Rk​⇒Ak+1​=Rk​Qk​

则：

Ak+1=QkTAkQkA_{k+1} = Q_k^T A_k Q_kAk+1​=QkT​Ak​Qk​

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

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

---

三、算法步骤

原始 QR 算法（无加速）

1. 初始化：$A_0 = A$

2. 对 $A_k$ 做 QR 分解：

   Ak=QkRkA_k = Q_k R_kAk​=Qk​Rk​

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+μkIA_k - \\mu_k I = Q_k R_k \\Rightarrow A_{k+1} = R_k Q_k + \\mu_k IAk​−μk​I=Qk​Rk​⇒Ak+1​=Rk​Qk​+μk​I

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

---

四、应用场景

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

---

五、C++ 实现示例（使用 Eigen 库）

示例1

#include #include 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
• #include

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

#include #include #include 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，不要使用 MatrixXd 求解过程 compute() 先分解，solve() 解方程组

---

---

六、注意事项与扩展

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

---