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 withEigen::FullPivHouseholderQR
.I also looked into the source
HouseholderQR.h
, and they compute the determinant viam_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 asQR.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 equalA
; although one column is correct and another is wrong.