对特征QR分解感到困惑

Question

我对特征的QR分解感到困惑。我的理解是，矩阵Q被隐式存储为一系列居民转换，并且矩阵R被存储为上三角矩阵，并且R的对角线包含A的特征值我所关心的(。

但是，我编写了以下程序，该程序通过两种使用Eigen::EigenSolver的方法计算矩阵A的特征值，另一个使用QR。我知道我的QR方法正在返回错误的结果，并且EigenSolver方法正在返回正确的结果。

我在这里误会了什么？

#include <iostream>
#include <algorithm>
#include <Eigen/Dense>
int main()
{
    using Real = long double;
    long n = 2;
    Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> A(n,n);
    for(long i = 0; i < n; ++i) {
        for (long j = 0; j < n; ++j) {
            A(i,j) = Real(1)/Real(i+j+1);
        }
    }
    auto QR = A.householderQr();
    auto Rdiag = QR.matrixQR().diagonal().cwiseAbs();
    auto [min, max] = std::minmax_element(Rdiag.begin(), Rdiag.end());
    std::cout << "u03BAu2082(A) = " << (*max)/(*min) << "n";
    std::cout << "u2016Au2016u2082 via QR = " << (*max) << "n";
    std::cout << "Diagonal of R =n" << Rdiag << "n";

    // dblcheck:
    Eigen::SelfAdjointEigenSolver<decltype(A)> eigensolver(A);
    if (eigensolver.info() != Eigen::Success) {
        std::cout << "Something went wrong.n";
        return 1;
    }
    auto absolute_eigs = eigensolver.eigenvalues().cwiseAbs();
    auto [min1, max1] = std::minmax_element(absolute_eigs.begin(), absolute_eigs.end());
    std::cout << "u03BAu2082(A) via eigensolver = " << (*max1)/(*min1) << "n";
    std::cout << "u2016Au2016u2082 via eigensolver = " << (*max1) << "n";
    std::cout << "The absolute eigenvalues of A via eigensolver are:n" << absolute_eigs << "n";
}

输出：

κ₂(A) = 15
‖A‖₂ via QR = 1.11803
Diagonal of R =
  1.11803
0.0745356
κ₂(A) via eigensolver = 19.2815
‖A‖₂ via eigensolver = 1.26759
The absolute eigenvalues of A via eigensolver are:
0.0657415
  1.26759

其他信息：

• 从头开始克隆特征：

$ hg log | more
changeset:   11993:20cbc6576426
tag:         tip
date:        Tue May 07 16:44:55 2019 -0700
summary:     Fix AVX512 & GCC 6.3 compilation

• 使用有或没有-ffast-math的G -8，G -9和Apple Clang汇编时会发生。我获得了Eigen::FullPivHouseholderQR的错误结果。

• 我还研究了源HouseholderQR.h，他们通过m_qr.diagonal().prod()计算了行列式，这使我对正确使用API的使用感到更加自信。从eigensolver中获取特征值的乘积返回与QR.absDeterminant()相同的值。

• 以下代码片段不返回原始矩阵a ：

Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> R = QR.matrixQR();
Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> Q = QR.householderQ();
std::cout << "Q*R = n" << Q*R << "n";

• 我检查了Q具有所有必要的属性：q^-1 = q^t，q^tq = i，和| det(q(|= 1。

• 我还验证了 QR.householderQ().transpose()*QR.matrixQR()不相等的 A;尽管一列是正确的，而另一列是错误的。

Answer 1

正如@Geza指出的那样，QR分解的R矩阵将(通常(不包含原始矩阵的特征值，如果是这种情况，生活将太容易了

对于您的另一个问题，如果您想从Q重建A和R您只需要查看QR.matrixQR()的上部三角形部分

Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>
              R = QR.matrixQR().triangularView<Eigen::Upper>();

除此之外，我建议您谨慎使用auto与表达模式结合使用(在您的情况下没有严重错误，除了评估Rdiag至少两次(。

此外，在现代CPU上，使用long double几乎不是一个好主意。首先确保您使用的算法在数值上是稳定的，并且如果精度确实是一个问题，请考虑使用任意的精确浮子(如MPFR(。