Смущенный по поводу собственного разложения QR - PullRequest
Новая
       44

Смущенный по поводу собственного разложения QR

3 голосов
/ 08 мая 2019

Я запутался насчет 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;хотя один столбец правильный, а другой неправильный.

1 Ответ

0 голосов
/ 11 мая 2019

Как указала @geza, матрица R QR-разложения (в общем случае) не будет содержать собственных значений исходной матрицы, жизнь была бы слишком простой, если бы это было так:)

К вашей другой проблеме, если вы хотите восстановить A из Q и R, вам нужно смотреть только на верхнюю треугольную часть QR.matrixQR()

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

Кроме того, я бы посоветовал быть осторожным с использованием auto в сочетании с шаблонами выражений (в вашем случае ничего серьезного, кроме того, что Rdiag оценивается как минимум дважды).

Кроме того, использование long double - едва ли хорошая идея для современных процессоров. Сначала убедитесь, что используемые вами алгоритмы численно устойчивы, и, если точность действительно является проблемой, рассмотрите возможность использования произвольных значений точности (например, MPFR).

...