Я запутался насчет QR-разложения Эйгена.Насколько я понимаю, матрица Q
хранится неявно как последовательность преобразований Домохозяина, а матрица R
хранится в виде верхней треугольной матрицы и что диагональ R
содержит собственные значения A
(по крайней мере, до фазы, и это все, что меня волнует).
Однако я написал следующую программу, которая вычисляет собственные значения матрицы A
двумя разными методами, один из которых Eigen::EigenSolver
, идругой использует 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 << "\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()
, что заставляет меня чувствовать себя несколько увереннеечто я правильно использую 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
;хотя один столбец правильный, а другой неправильный.