Dense linear problems and decompositions » Inplace matrix decompositions module

Starting from Eigen 3.3, the LU, Cholesky, and QR decompositions can operate inplace, that is, directly within the given input matrix. This feature is especially useful when dealing with huge matrices, and or when the available memory is very limited (embedded systems).

To this end, the respective decomposition class must be instantiated with a Ref<> matrix type, and the decomposition object must be constructed with the input matrix as argument. As an example, let us consider an inplace LU decomposition with partial pivoting.

Let's start with the basic inclusions, and declaration of a 2x2 matrix A:

code	output
#include <iostream> #include <Eigen/Dense> using namespace std; using namespace Eigen; int main() { MatrixXd A(2,2); A << 2, -1, 1, 3; cout << "Here is the input matrix A before decomposition:\n" << A << endl;	Here is the input matrix A before decomposition: 2 -1 1 3

code

output

#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main()
{
  MatrixXd A(2,2);
  A << 2, -1, 1, 3;
  cout << "Here is the input matrix A before decomposition:\n" << A << endl;

Here is the input matrix A before decomposition:
 2 -1
 1  3

No surprise here! Then, let's declare our inplace LU object lu, and check the content of the matrix A:

  PartialPivLU<Ref<MatrixXd> > lu(A);
  cout << "Here is the input matrix A after decomposition:\n" << A << endl;

Here is the input matrix A after decomposition:
  2  -1
0.5 3.5

Here, the lu object computes and stores the L and U factors within the memory held by the matrix A. The coefficients of A have thus been destroyed during the factorization, and replaced by the L and U factors as one can verify:

  cout << "Here is the matrix storing the L and U factors:\n" << lu.matrixLU() << endl;

Here is the matrix storing the L and U factors:
  2  -1
0.5 3.5

Then, one can use the lu object as usual, for instance to solve the Ax=b problem:

  MatrixXd A0(2,2); A0 << 2, -1, 1, 3;
  VectorXd b(2);    b << 1, 2;
  VectorXd x = lu.solve(b);
  cout << "Residual: " << (A0 * x - b).norm() << endl;

Residual: 0

Here, since the content of the original matrix A has been lost, we had to declared a new matrix A0 to verify the result.

Since the memory is shared between A and lu, modifying the matrix A will make lu invalid. This can easily be verified by modifying the content of A and trying to solve the initial problem again:

  A << 3, 4, -2, 1;
  x = lu.solve(b);
  cout << "Residual: " << (A0 * x - b).norm() << endl;

Residual: 15.8114

Note that there is no shared pointer under the hood, it is the responsibility of the user to keep the input matrix A in life as long as lu is living.

If one wants to update the factorization with the modified A, one has to call the compute method as usual:

  A0 = A; // save A
  lu.compute(A);
  x = lu.solve(b);
  cout << "Residual: " << (A0 * x - b).norm() << endl;

Residual: 0

Note that calling compute does not change the memory which is referenced by the lu object. Therefore, if the compute method is called with another matrix A1 different than A, then the content of A1 won't be modified. This is still the content of A that will be used to store the L and U factors of the matrix A1. This can easily be verified as follows:

  MatrixXd A1(2,2);
  A1 << 5,-2,3,4;
  lu.compute(A1);
  cout << "Here is the input matrix A1 after decomposition:\n" << A1 << endl;

Here is the input matrix A1 after decomposition:
 5 -2
 3  4

The matrix A1 is unchanged, and one can thus solve A1*x=b, and directly check the residual without any copy of A1:

  x = lu.solve(b);
  cout << "Residual: " << (A1 * x - b).norm() << endl;

Residual: 2.48253e-16

Here is the list of matrix decompositions supporting this inplace mechanism: