求解平方根反比

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

该方法stable_cholesky_solver使用使用枢轴的稳定分解lldt()求解平方根。main验证它是否做了它应该做的任何事情，并提供了一种实现相同目标的方法，使用不太稳定（但更快）的llt()分解。请参阅文档的前几行以了解我的 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 << "nnUsing 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 << "nnUsing 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;
}