Solve for inverse square root

Question

I have a positive definite matrix A of which I already computed the cholesky decomposition: A=LDL^T. For some vector x, I would like to compute S^{-1}x, where S is a square root of A. For now, I do

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

Is this a stable way to do this computation? I though computing inverses was a bad thing in general. Is there a way to use the already performed LDLT decomposition? I feel it is possible, because that's what's actually happening behind the scenes in LDLT::solve()!

Answer 1

Here is complete code to solve the problem for a symmetric matrix A and a general right hand side b (a vector or a matrix). I was not able to find anything online that I can play with (or just copy-paste) so I wrote one.

The method stable_cholesky_solver does the job of solving for the square root using the stable decomposition lldt() that uses pivoting. The main verifies that it does whatever it is supposed to do and presents a way to achieve the same goal, using the less stable (but faster) llt() decomposition. See the first few lines of the documentation to understand my notation of 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;

}