对 Eigen QR 分解感到困惑

Question

我对 Eigen 的 QR 分解感到困惑。我的理解是矩阵 Q 隐式存储为一系列 Householder 变换，矩阵 R 存储为上三角矩阵，并且 R 的对角线包含A 的特征值（至少到阶段，这是我所关心的）。

但是，我编写了以下程序，通过两种不同的方法计算矩阵 A 的特征值，一种使用 Eigen::EigenSolver，另一种使用 QR。我知道我的 QR 方法是 return 错误的结果，而 EigenSolver 方法是 return 正确的结果。

我误解了什么？

#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 << "\u03BA\u2082(A) = " << (*max)/(*min) << "\n";

    std::cout << "\u2016A\u2016\u2082 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 << "\u03BA\u2082(A) via eigensolver = " << (*max1)/(*min1) << "\n";
    std::cout << "\u2016A\u2016\u2082 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

在使用 g++-8、g++-9 和 Apple Clang 编译时发生，有和没有 -ffast-math。我用 Eigen::FullPivHouseholderQR.
我还查看了来源 HouseholderQR.h，他们通过 m_qr.diagonal().prod() 计算行列式，这让我对使用 [=71] 更有信心=] 正确。从 EigenSolver returns 中获取与 QR.absDeterminant().
下面的代码片段没有return原始矩阵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 和 R 重建 A 你只需要看 QR.matrixQR()[= 的上三角部分19=]

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

除此之外，我建议谨慎使用 auto 和 expression-templates （在你的情况下没有严重错误，除了 Rdiag 至少被评估两次） .

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

对 Eigen QR 分解感到困惑

Confused about Eigen QR decomposition

c++

linear-algebra

eigen

eigen3

eigenvalue