Собственная эффективная обратная симметричная положительно определенная матрица

В Эйгене, если у нас есть симметричная положительно определенная матрица A, мы можем вычислить обратную A по формуле

A.inverse();

or

A.llt().solve(I);

где I — единичная матрица того же размера, что и A. Но есть ли более эффективный способ вычисления обратной симметричной положительно определенной матрицы?

Например, если мы запишем разложение Холецкого A как A = LL^{T}, то L^{-T} L^{-1} будет инверсией A, начиная с A L^{-T} L^{-1} = LL^{T} L^{-T} L^{-1} = I (и где L^{-T} обозначает инверсию транспонирования L).

Таким образом, мы могли бы получить разложение Холецкого для A, вычислить его обратное, а затем получить векторное произведение этого обратного, чтобы найти обратное A. Но я чувствую, что вычисление этих явных шагов будет медленнее, чем использование A.llt().solve(I), как указано выше.

И прежде чем кто-нибудь спросит, мне действительно нужна явная инверсия — это вычисление части сэмплера Гиббса.


c++ eigen matrix-inverse
person dpritch    schedule 28.07.2016    source источник
comment
Хотя в принятом ответе не указано, существует ли более быстрый способ вычисления обратной симметричной положительно определенной матрицы с собственными числами, на мой взгляд, явный метод, который я упомянул в вопросе, является самым быстрым способом сделать это (и имеет порядок O (( 1/3)n^3 + 2n^2) ) - очевидно, что Эйген делает под капотом.   -  person dpritch    schedule 29.07.2016

Ответы (2)


arrow_upward
6
arrow_downward

С помощью A.llt().solve(I) вы предполагаете, что A является матрицей SPD, и применяете разложение Холецкого для решения уравнения Ax=I. Математическая процедура решения уравнения точно такая же, как и ваш явный способ. Таким образом, производительность должна быть такой же, если вы делаете каждый шаг правильно.

С другой стороны, с A.inverse() выполняется общая инверсия матрицы, в которой используется LU-разложение для большой матрицы. Таким образом, производительность должна быть ниже A.llt().solve(I);.

person kangshiyin    schedule 28.07.2016
comment
Спасибо, это сэкономило мне кучу времени на профилировании метода решения Холецкого по сравнению с явным методом, и в конце дня я все еще не знаю наверняка, что они одинаковы. - person dpritch; 29.07.2016

arrow_upward
0
arrow_downward

Вы должны профилировать код для вашей конкретной проблемы, чтобы получить лучший ответ. Я тестировал код, пытаясь оценить жизнеспособность обоих подходов, используя библиотеку googletest и этот репозиторий:

#include <gtest/gtest.h>

#define private public
#define protected public

#include <kalman/Matrix.hpp>
#include <Eigen/Cholesky>
#include <chrono>
#include <iostream>

using namespace Kalman;
using namespace std::chrono;

typedef float T;
typedef high_resolution_clock Clock;

TEST(Cholesky, inverseTiming) {
    Matrix<T, Dynamic, Dynamic> L;
    Matrix<T, Dynamic, Dynamic> S;
    Matrix<T, Dynamic, Dynamic> Sinv_method1;
    Matrix<T, Dynamic, Dynamic> Sinv_method2;
    int Nmin = 2;
    int Nmax = 128;
    int N(Nmin);

    while (N <= Nmax) {

        L.resize(N, N);
        L.setRandom();
        S.resize(N, N);
        // create a random NxN SPD matrix
        S = L*L.transpose();
        std::cout << "\n";
        std::cout << "+++++++++++++++++++++++++ N = " << N << " +++++++++++++++++++++++++++++++++++++++" << std::endl;
        auto t1 = Clock::now();
        Sinv_method1.resize(N, N);
        Sinv_method1 = S.inverse();
        auto dt1 = Clock::now() - t1;
        std::cout << "Method 1 took " << duration_cast<microseconds>(dt1).count() << " usec" << std::endl;
        auto t2 = Clock::now();
        Sinv_method2.resize(N, N);
        Sinv_method2 = S.llt().solve(Matrix<T, Dynamic, Dynamic>::Identity(N, N));
        auto dt2 = Clock::now() - t2;
        std::cout << "Method 2 took " << duration_cast<microseconds>(dt2).count() << " usec" << std::endl;
        for(int i = 0; i < N; i++)
        {
            for(int j = 0; j < N; j++)
            {
                EXPECT_NEAR( Sinv_method1(i, j), Sinv_method2(i, j), 1e-3 );
            }
        }

        N *= 2;
        std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
        std::cout << "\n";
    }
}

Приведенный выше пример показал мне, что для моей проблемы с размером ускорение было незначительным при использовании method2, тогда как отсутствие точности (использование вызова .inverse() в качестве эталона) было заметным.

person joe.dinius    schedule 18.10.2018