对称正定矩阵的Eigen高效逆矩阵

Question

在Eigen中，如果我们有对称正定矩阵 A，则可以通过以下方式计算其逆矩阵：

A.inverse();

或者

A.llt().solve(I);

I是与A大小相同的单位矩阵。但是有没有更有效的方法来计算对称正定矩阵的逆呢？例如，如果我们将A的Cholesky分解写为A = LL^{T}，那么L^{-T} L^{-1}就是A的逆，因为A L^{-T} L^{-1} = LL^{T} L^{-T} L^{-1} = I（其中L^{-T}表示L的转置的逆）。

因此，我们可以获得A的Cholesky分解，计算其逆，然后获得该逆的交叉乘积以找到A的逆。但我的直觉是，计算这些显式步骤会比上面的A.llt().solve(I)慢。

在任何人问之前，我确实需要一个显式的逆 - 这是Gibbs采样器的一部分计算。

- dpritch

---

尽管被接受的答案没有说明是否有更快的方法来计算Eigen中对称正定矩阵的逆，但据我所知，我在问题中提到的显式方法是最快的方法（其具有O((1/3)n^3+2n^2)的阶数） - 这显然就是Eigen在幕后所做的。 - dpritch

Answer 1

使用 A.llt().solve(I)，您假设矩阵A是对称正定矩阵，并应用Cholesky分解来解方程Ax=I。数学求解过程与您显式的方法完全相同，因此如果每个步骤都正确执行，则性能应该相同。

另一方面，使用A.inverse()，您正在进行通用矩阵求逆，针对大型矩阵使用LU分解。因此，性能应低于A.llt().solve(I);。

- kangshiyin

---

谢谢，这节省了我大量时间来分析Cholesky求解方法和显式方法，在一天结束时仍然不能确定它们是否相同。 - dpritch

Answer 2

你应该针对你的具体问题对代码进行性能分析，以获得最佳答案。我在使用googletest库和this repo时对代码进行了基准测试，以评估两种方法的可行性：

#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()作为基准的精度不足是明显的。

- joe.dinius