Confused about Eigen QR decomposition

Question

I am confused about Eigen's QR decomposition. My understanding is that the matrix Q is stored implicitly as a sequence of Householder transformations, and that the matrix R is stored as an upper triangular matrix, and that the diagonal of R contains the eigenvalues of A (at least up to phase, which is all I care about).

However, I wrote the following program which computes the eigenvalues of a matrix A via two different methods, one using the Eigen::EigenSolver, and the other using QR. I know that my QR method is returning the wrong results, and that the EigenSolver method is returning the correct results.

What am I misunderstanding here?

#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";
}

Output:

κ₂(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

Other info:

• Cloned eigen from the head:

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

• Occurs when compiled with g++-8, g++-9, and Apple Clang, with and without -ffast-math. I obtain the same wrong result with Eigen::FullPivHouseholderQR.

• I also looked into the source HouseholderQR.h, and they compute the determinant via m_qr.diagonal().prod(), which makes me feel somewhat more confident that I'm using the API correctly. Taking the product of the eigenvalues from the EigenSolver returns the same values as QR.absDeterminant().

• The following code snippet does not return the original matrix 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";

• I checked that Q has all the requisite properties: Q^-1 = Q^T, Q^TQ = I, and |det(Q)| = 1.

• I've also verified that QR.householderQ().transpose()*QR.matrixQR() is not equal A; although one column is correct and another is wrong.

Answer 1

As @geza pointed out, the R matrix of a QR decomposition will (in general) not contain the Eigenvalues of the original matrix, life would be too easy if that was the case :)

To your other problem, if you want to reconstruct A from Q and R you need to only look at the upper triangular part of QR.matrixQR()

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

Besides that, I'd suggest being careful with using auto in combination with expression-templates (nothing severely wrong in your case, except that Rdiag is evaluated at least twice).

Also, using long double is barely a good idea on modern CPUs. First make sure that the algorithms you use are numerically stable and if precision really is an issue, consider using arbitrary precision floats (like MPFR).