1

I am decomposing a sparse SPD matrix A using Eigen. It will either be a LLt or a LDLt deomposition (Cholesky), so we can assume the matrix will be decomposed as A = P-1 LDLt P where P is a permutation matrix, L is triangular lower and D diagonal (possibly identity). If I do

SolverClassName<SparseMatrix<double> > solver;
solver.compute(A);

To solve Lx=b then is it efficient to do the following?

solver.matrixL().TriangularView<Lower>().solve(b)

Similarly, to solve Px=b then is it efficient to do the following?

solver.permutationPinv()*b

I would like to do this in order to compute bt A-1 b efficiently and stably.

yannick
  • 397
  • 3
  • 19

1 Answers1

0

Have a look how _solve_impl is implemented for SimplicialCholesky. Essentially, you can simply write:

Eigen::VectorXd x = solver.permutationP()*b; // P not Pinv!
solver.matrixL().solveInPlace(x); // matrixL is already a triangularView

// depending on LLt or LDLt use either:
double res_llt = x.squaredNorm();
double res_ldlt = x.dot(solver.vectorD().asDiagonal().inverse()*x);

Note that you need to multiply by P and not Pinv, since the inverse of A = P^-1 L D L^t P is 

P^-1 L^-t D^-1 L^-1 P

because the order of matrices reverses when taking the inverse of a product.

chtz
  • 17,329
  • 4
  • 26
  • 56