c++ - 求解平方根反比

Question

我有一个正定矩阵 A，我已经计算了它的 cholesky 分解:A=LDL^T。对于一些 vector x，我想计算 S^{-1}x，其中 S 是 A 的平方根。现在，我做

Eigen::SelfadjointEigenSolver<Eigen::MatrixXd> es(A);
Eigen::MatrixXd Si(es.operatorInverseSqrt());
return Si*get_x();

这是进行此计算的稳定方法吗？我虽然计算逆通常是一件坏事。有没有办法使用已经执行的 LDLT 分解？我觉得这是可能的，因为这就是 LDLT::solve() 幕后实际发生的事情!

Answer 1

这是解决对称矩阵 A 和一般右侧 b( vector 或矩阵)问题的完整代码。我无法在网上找到任何可以玩(或只是复制粘贴)的东西，所以我写了一个。

stable_cholesky_solver 方法使用旋转的稳定分解 lldt() 求解平方根。 main 验证它做了它应该做的任何事情，并提出了一种实现相同目标的方法，使用不太稳定(但更快)的 llt() 分解。查看 documentation 的前几行理解我对 L、P、D 的记法。

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;

Matrix<double, Dynamic, Dynamic> stable_cholesky_solver( 
                     LDLT<MatrixXd> ldltDecomp,
                     Matrix<double, Dynamic, Dynamic> A,
                         bool transpose = false )
{

    // Preparations:

    // For some reason if I sub it below I get error
    Matrix<double, Dynamic, Dynamic> L = ldltDecomp.matrixL();

    // Number of rows is all that matters, regardless if rhs is a
    // matrix or a vector
    int k = A.rows(); 

    // Manually inverting D. This procedure has the advantage that
    // D^{-1/2} can also be applied to matrices.
    VectorXd diag;
    diag.resize(k);
    for( int i = 0 ; i < k ; ++i ) 
      diag(i) = 1. / sqrt( ldltDecomp.vectorD()(i) ) ; // Manual inversion
    DiagonalMatrix<double, Dynamic > sqrtInvD = diag.asDiagonal();

    // The permutation "matrix" P
    Transpositions<Dynamic> P = ldltDecomp.transpositionsP(); 

    // Holds all the computations
    Matrix<double, Dynamic, Dynamic> x;

    // Now the actual computation

    if( !transpose ) {

        // x = PA
        x = P * A;

        // x = L^{-1}PA
        x = L.triangularView<Lower>().solve<OnTheLeft>(x);

        // x = D^{-1/2}L^{-1}PA
        x = sqrtInvD * x;

    } else {

        // x = D^{-1/2}A
        x = sqrtInvD * A;

        // x = L^{-t}D^{-1/2}A
        x = L.triangularView<Lower>().transpose().solve<OnTheLeft>(x);

        // x = P^tL^{-t}D^{-1/2}A
        x = P.transpose() * x; 
   }
   return x;

}



int main()
{

    int k = 3; // Dimensionality

    // Define, declare and enter the problem's data
    MatrixXd A;
    A.resize(k, k);
    MatrixXd b;
    b.resize(k, 2 );

    A <<
      13, 5, 7 ,
      5 , 9, 3 ,
      7 , 3, 11;
    b <<
      3, 3, 4,
      1,-2, 9;

    cout << "Here is the " << A.rows() << " by " << A.cols() << " matrix A:\n" << A << endl;
    cout << "Here is the " << b.rows() << " by " << b.cols() << " matrix b:\n" << b << endl;
    cout << "Let's solve Ax = b using different methods.\n" <<endl;

    // Two placeholders that will be used throughout
    MatrixXd L;
    MatrixXd x;

    // ldlt()
    cout << "\n\nUsing the stable Cholesky decompostion ldlt()" << endl;

    // The object that actually holds the entire decomposition
    LDLT<MatrixXd> ldltDecomp = A.ldlt();

    // Direct solution, using Eigen's routines
    x = ldltDecomp.solve(b);
    cout << "Direct x =\n" << x << endl;
    cout << "Direct b =\n" << A*x << endl;

    // Manual solution - implementing the process Eigen is taking, step
    // by step (in the function defined above). 
    x = stable_cholesky_solver( ldltDecomp, b );
    x = stable_cholesky_solver( ldltDecomp, x, true );

    cout << "Manual x =\n" << x << endl;
    cout << "Manual b =\n" << A*x << endl;


    // llt()
    cout << "\n\nUsing the less stable, but faster Cholesky decomposition " << "without pivoting llt()" << endl;

    // Here A.llt() is the object that actually holds the decomposition
    // (like ldltDecomp before) but we only need the matrix L.
    L = A.llt().matrixL();

    x = L.triangularView<Lower>().solve<OnTheLeft>(b);
    x = L.triangularView<Lower>().transpose().solve<OnTheLeft>(x);
    cout << "Manual x =\n" << x << endl;
    cout << "Manual b =\n" << A*x << endl;

}