对称正定矩阵的特征有效逆

Eigen efficient inverse of symmetric positive definite matrix

本文关键字：特征 有效 对称      更新时间：2023-10-16

在Eigen中，如果我们有对称正定矩阵A，那么我们可以通过计算A的逆

A.inverse();

或

A.llt().solve(I);

其中CCD_ 3是与CCD_ 4大小相同的单位矩阵。但是，有没有一种更有效的方法来计算对称正定矩阵的逆？

例如，如果我们将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)慢。

在有人问之前，我确实需要一个明确的逆——这是吉布斯采样器的一部分的计算。

使用A.llt().solve(I)，假设A是SPD矩阵，并应用Cholesky分解来求解方程Ax=I。解方程的数学过程和你的显式方法完全一样。因此，如果你每一步都做得正确，性能应该是一样的。

另一方面，使用A.inverse()进行一般矩阵求逆，它对大矩阵使用LU分解。因此性能应该低于CCD_ 19。

您应该为您的特定问题分析代码，以获得最佳答案。我在测试代码的同时，试图使用谷歌测试库和这个回购来评估这两种方法的可行性：

#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()调用作为基准）是显而易见的。

相关文章：