Confused about Eigen QR decomposition
Asked Answered
W

1

6

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.

Wadsworth answered 8/5, 2019 at 19:6 Comment(2)
I don't know Eigen. But QR decomposition doesn't give you the eigenvalues. It needs further processing to get eigenvalues from QR decomposition. It's just the case that the product of diagonal is the determinant (because R is triangular, and Q has a determinant of plus/minus one).Swop
@geza: You're right. Probably ought to delete this question.Wadsworth
P
4

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).

Phoney answered 11/5, 2019 at 18:44 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.