求解反平方根

Question

我有一个正定矩阵A，其中我已经计算了胆甾醇分解：A = LDL^T。 对于一些向量x，我想计算小号^ { - 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（矢量或矩阵）的问题。我无法在网上找到任何可以玩的东西（或只是复制粘贴），所以我写了一个。

方法stable_cholesky_solver做使用pivoting的稳定分解lldt()解决平方根的工作。 main验证它做它应该做的任何事情，并提出一种方法来实现相同的目标，使用不太稳定（但更快）的分解。 请参阅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; 

}