当前位置：首页 > news >正文

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

news 2025/7/24 13:15:38

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

一、定义与目标

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

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

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

$A_k = Q_k R_k$

然后重新组合得到：

$A_{k+1} = R_k Q_k$

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

二、原理与推导

假设 $A$ 可对角化：

$\Lambda V^{-1}$

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

$A_k = Q_k R_k \Rightarrow A_{k+1} = R_k Q_k$

则：

$A_{k+1} = Q_k^T A_k Q_k$

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

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

三、算法步骤

原始 QR 算法（无加速）

初始化：$A_0 = A$
对 $A_k$ 做 QR 分解：

$A_k = Q_k R_k$
计算 $A_{k+1} = R_k Q_k$
判断是否收敛（如 $|A_{k+1} - A_k| < \epsilon$）
重复步骤 2~4

改进：带位移的 QR 算法（Shifted QR）

通过引入位移（shift）$\mu_k$ 加快收敛速度：

$A_k - \mu_k I = Q_k R_k \Rightarrow A_{k+1} = R_k Q_k + \mu_k I$

位移常选为 $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());// 右侧向量 bVectorXd 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 = bVectorXd 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()` 解方程组