3
$\begingroup$

My problem needs to solve dense symmetric linear systems something like:

A x = b, A y = x, A z = y+x,...in sequence.

In Eigen C++, if I take advantage of the symmetry of A by using:

x=A.ldlt().solve(b);
y=A.ldlt().solve(x);
z=A.ldlt().solve(y+x);

The same matrix A has to be Cholesky factorized many times.

Similar issue may exists for householderQR and LU algorithms.

How can I reuse the coefficient matrix decomposition result for other systems?

It seems this is equivalent to ask whether there is a simple back-substitution solver for upper triangular linear systems.

$\endgroup$
6
$\begingroup$

You need to declare a LDLT object like this:

LDLT<MatrixXd> ldlt(A);
x = ldlt.solve(b);
y = ldlt.solve(x);
...
$\endgroup$
  • $\begingroup$ Thank you. So if I have A changed in a loop, then just re-declare it? $\endgroup$ – LCFactorization Dec 24 '13 at 15:51
  • 1
    $\begingroup$ There is a compute method like this: ldlt.compute(B); $\endgroup$ – ggael Dec 25 '13 at 10:00
2
$\begingroup$

Given that the accepted answer has typos and does not give full code, I decided to post mine.

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
    MatrixXd A;
    A.resize(3,3);
    VectorXd b;
    b.resize(3);

    A <<
      13, 5, 7 ,
      5 , 9, 3 ,
      7 , 3, 11;
    b << 3, 3, 4;

    cout << "Here is the matrix A:\n" << A << endl;
    cout << "Here is the vector b:\n" << b << endl;
    LDLT<MatrixXd> L = A.ldlt();
    VectorXd x = L.solve(b);
    cout << "The solution is:\n" << x << endl;
    cout << "reconstructed b is:\n"  << A*x << endl; 

}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.