减少本征的 QR 分解

Question

我正在尝试使用特征argmin (Ax - b)^2解决经典的最小二乘问题。我知道系统是过度确定的(即，m >= n其中 A 是m x n(，并且 A 具有完整的等级(n(。我选择使用本征的稀疏 QR 部分，因为 A 本身就是稀疏的。我得到了这个：

#include <Eigen/Sparse>
#include <Eigen/SparseQR>
void solve_least_squares(const Eigen::SparseMatrix<double>& matrix,
const Eigen::VectorXd& rhs)
{
int m = matrix.rows();
int n = matrix.cols();
assert(m >= n);
assert(matrix.isCompressed());
Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> decomposition(matrix);
Eigen::VectorXd t = (decomposition.matrixQ().transpose() * rhs).topRows(n);
std::cout << t << std::endl;
Eigen::VectorXd s = decomposition.matrixR().topRows(n).triangularView<Eigen::Upper>().solve(t);
return decomposition.colsPermutation().inverse()*s;
}

但我想知道这是否是最有效的实现：首先，Eigen 似乎确定了完整的QR 分解而不是简化的分解(Q是m x m而不是m x n(。这似乎很浪费，因为我只需要前n行。

在密集的情况下，HouseholderQR 类中有任何示例：

MatrixXf A(MatrixXf::Random(5,3)), thinQ(MatrixXf::Identity(5,3)), Q;
A.setRandom();
HouseholderQR<MatrixXf> qr(A);
Q = qr.householderQ();
thinQ = qr.householderQ() * thinQ; // <- here
std::cout << "The complete unitary matrix Q is:n" << Q << "nn";
std::cout << "The thin matrix Q is:n" << thinQ << "nn";

这里，使用矩阵乘法来获得约Q，这似乎比切片整个矩阵更有效率。

由于特征的(密集(SVD 分解提供了计算薄 SVD 的选项，因此 QR 分解是否有类似的选项？

Answer 1

在分解中，Q因子存储为一系列 Householder 向量，matrixQ()方法实质上返回对该矩阵的引用(它以"好像"乘以实际矩阵的方式重载乘法(。当存储为 Householder 矩阵时，Q表示Q的全部或薄部分没有区别(实际上，将向量乘以完整矩阵可以有效地就地完成(。

如果你想求解你的线性系统，你应该写

Eigen::VectorXd s = decomposition.solve(rhs);

---

顺便说一句：如果你可以失去一些精度(并且你的矩阵条件足够好(，那么通过求解正态方程(Matlab符号(可以(在大多数情况下(更快：

x = (A'*A)  (A'*b) 

并使用稀疏的乔列斯基分解(例如，Eigen::SimplicialLLT或Eigen::SimplicialLDLT(。但是用你的实际数据进行基准测试，并检查准确性是否足以满足你的用例(也许在细化步骤之后，重新使用Cholesky分解(。