Hi,
I have been trying with occasional success to invert either a cholesky or ldl decomposition.

Basically I require the A^-1 from a given symmetric positive-definite matrix A.

However my test code (A^-1*A = I) only occasionally results in the identity matrix as expected.

For an ldl I do the following

Code: Select all

        // get the inverse of the linear system
       Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> Li = ldlt.matrixL();
       const typename Eigen::LDLT<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>, Eigen::Upper>::TranspositionType& trans = ldlt.transpositionsP();
        Li = (trans.transpose()*Li).inverse();
       const Eigen::Diagonal<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> >& vectorD = ldlt.vectorD();
        T tolerance = std::max(vectorD.array().abs().maxCoeff()*Eigen::NumTraits<T>::epsilon(), T(1) / Eigen::NumTraits<T>::highest()); // motivated by LAPACK's xGELSS
       for(int i = 0; i < vectorD.size(); ++i) {
//         const T& v = vectorD(trans(i));
         const T& v = vectorD(i);
            std::cout << "D[" << trans(i) << "][" << v << "]" << std::endl;
           if(std::abs(v) > tolerance)
              Li.row(i) /= std::sqrt(v);
           else
              Li.row(i).setZero();
        }
        Ai_ = Li.transpose()*Li;


#if 1
        Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> test2 = A_*Ai_;
        for(int i = 0; i < vars_n_; ++i) {
           if(std::abs(test2(i, i)-static_cast<T>(1)) > T(1e-3)) {
               std::cout << "Inverse2 Failed[" << i << "] : " << std::abs(test2(i, i)-static_cast<T>(1)) << std::endl;
             std::cout << test2 << std::endl;            
             break;
         }
        }
#endif


Can anyone point out how I might be able to correct this code to obtain consistently correct solution for A^-1.

Similar erroneous results are obtained with a cholesky decomposition.

Cheers,
David...

I did not check your code but the following should do the job:

A_inv = A.ldlt().solve(Matrix<T,Dynamic,Dynamic>::Identity(n,n));

Note that it is not recommended to compute the inverse of a matrix, it is often preferable to use the solve() function.

Thank you for your suggestion, your code is a simpler way to achieve the same result.

Unfortunately I am still getting failure cases. Here is an example, performing these operations in Sage I get

A = Matrix([
[ 105264, 239778, -2179.68, -55000.9, 23846.2, 1840.87],
[ 239778, 587291, -5428.58, -57529.7, 26543.8, 4892.12],
[ -2179.68, -5428.58, 852466, -171106, 86453.6, 28457.1],
[ -55000.9, -57529.7, -171106, 2.47209e+06, -999315, 51543.9],
[ 23846.2, 26543.8, 86453.6, -999315, 413856, -2507.14],
[ 1840.87, 4892.12, 28457.1, 51543.9, -2507.14, 35032]])

A.inverse()
[ 0.000143904042637198 -0.0000585381756708675 4.68631395737500e-7 -3.31880291656818e-6 -0.0000126232888954470 4.21177630828247e-6]
[-0.0000585381756708675 0.0000255231624212554 -1.58642539153700e-7 1.38267220886758e-6 5.09545501630961e-6 -2.02900189178896e-6]
[ 4.68631395737500e-7 -1.58642539153700e-7 1.23605133964706e-6 1.32116935115263e-8 -2.49457718309198e-7 -1.04382950871681e-6]
[ -3.31880291643506e-6 1.38267220880213e-6 1.32116935137451e-8 0.246118539178948 0.592351880792025 -0.319730569418160]
[-0.0000126232888951266 5.09545501615208e-6 -2.49457718303858e-7 0.592351880792025 1.42566016951603 -0.769519041500798]
[ 4.21177630810953e-6 -2.02900189170393e-6 -1.04382950871970e-6 -0.319730569418160 -0.769519041500798 0.415388798940202]

However using your suggested approach I get A^-1 =
0.00014391 -5.85409e-05 4.6865e-07 8.40998e-06 1.56051e-05 -1.1025e-05
-5.85409e-05 2.55243e-05 -1.58651e-07 -3.41675e-06 -6.45559e-06 4.20587e-06
4.6865e-07 -1.58651e-07 1.23605e-06 1.38723e-07 5.26197e-08 -1.20688e-06
8.40998e-06 -3.41675e-06 1.38723e-07 867773 2.08853e+06 -1.12732e+06
1.56051e-05 -6.45558e-06 5.26188e-08 2.08853e+06 5.02664e+06 -2.7132e+06
-1.1025e-05 4.20587e-06 -1.20688e-06 -1.12732e+06 -2.7132e+06 1.46448e+06

Note that the first 3x3 is about the same but the rest of the matrix differs. This kind of example is consistent with other larger matrices 6nx6n whereby the first 3nx3n portion appears to be calculated correctly but the rest is not.

Any ideas?

Sorry but I cannot reproduce, here my test program:

Code: Select all

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

int main()
{
  MatrixXd A(6,6), iA1, iA2;
  A << 105264, 239778, -2179.68, -55000.9, 23846.2, 1840.87,
  239778, 587291, -5428.58, -57529.7, 26543.8, 4892.12,
  -2179.68, -5428.58, 852466, -171106, 86453.6, 28457.1,
  -55000.9, -57529.7, -171106, 2.47209e+06, -999315, 51543.9,
  23846.2, 26543.8, 86453.6, -999315, 413856, -2507.14,
  1840.87, 4892.12, 28457.1, 51543.9, -2507.14, 35032;

  iA1 = A.llt().solve(MatrixXd::Identity(6,6));
  iA2 = A.ldlt().solve(MatrixXd::Identity(6,6));
  std::cout << "with LLT:\n" <<  iA1 << "\n\n";
  std::cout << "with LDLT:\n" <<  iA2 << "\n\n";

  return 0;
}



and the result:

Code: Select all

with LLT:
 0.000143904 -5.85382e-05  4.68631e-07  -3.3188e-06 -1.26233e-05  4.21178e-06
-5.85382e-05  2.55232e-05 -1.58643e-07  1.38267e-06  5.09546e-06   -2.029e-06
 4.68631e-07 -1.58643e-07  1.23605e-06  1.32117e-08 -2.49458e-07 -1.04383e-06
 -3.3188e-06  1.38267e-06  1.32117e-08     0.246119     0.592352    -0.319731
-1.26233e-05  5.09546e-06 -2.49458e-07     0.592352      1.42566    -0.769519
 4.21178e-06   -2.029e-06 -1.04383e-06    -0.319731    -0.769519     0.415389

with LDLT:
 0.000143904 -5.85382e-05  4.68631e-07  -3.3188e-06 -1.26233e-05  4.21178e-06
-5.85382e-05  2.55232e-05 -1.58643e-07  1.38267e-06  5.09546e-06   -2.029e-06
 4.68631e-07 -1.58643e-07  1.23605e-06  1.32117e-08 -2.49458e-07 -1.04383e-06
 -3.3188e-06  1.38267e-06  1.32117e-08     0.246119     0.592352    -0.319731
-1.26233e-05  5.09546e-06 -2.49458e-07     0.592352      1.42566    -0.769519
 4.21178e-06   -2.029e-06 -1.04383e-06    -0.319731    -0.769519     0.415389


Thanks for your further reply.

I had omitted some of the detail from my explanation of the problem. I was only filling the upper triangular portion of the matrix with non-zeros. To solve I was utilizing this fact as follows :

Code: Select all

        Ai_ = A_.template selfadjointView<Eigen::Upper>().ldlt().solve(Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>::Identity(vars_n_, vars_n_));


It seems that I was actually making a rather simple error in my test case I was trying to recover the identity by doing :

Code: Select all

Ai_*A_


instead of :

Code: Select all

Ai_*A_.template selfadjointView<Eigen::Upper>()


Your test case helped me focus in a little better on the problem, many thanks!