Tutorial 3/4 - Advanced linear algebra

Overview | Core features | Geometry | Advanced linear algebra | Sparse matrix

This tutorial chapter explains how you can use Eigen to tackle various problems involving matrices: solving systems of linear equations, finding eigenvalues and eigenvectors, and so on.

Table of contents

Solving linear problems

This part of the tutorial focuses on solving systems of linear equations. Such systems can be written in the form $ A \mathbf{x} = \mathbf{b} $, where both $ A $ and $ \mathbf{b} $ are known, and $ \mathbf{x} $ is the unknown. Moreover, $ A $ is assumed to be a square matrix.

The equation $ A \mathbf{x} = \mathbf{b} $ has a unique solution if $ A $ is invertible. If the matrix is not invertible, then the equation may have no or infinitely many solutions. All solvers assume that $ A $ is invertible, unless noted otherwise.

Eigen offers various algorithms to this problem. The choice of algorithm mainly depends on the nature of the matrix $ A $, such as its shape, size and numerical properties.

The methods described here can be used whenever an expression involve the product of an inverse matrix with a vector or another matrix: $ A^{-1} \mathbf{v} $ or $ A^{-1} B $.

LU decomposition (with partial pivoting)

This is a general-purpose algorithm which performs well in most cases (provided the matrix $ A $ is invertible), so if you are unsure about which algorithm to pick, choose this. The method proceeds in two steps. First, the LU decomposition with partial pivoting is computed using the MatrixBase::partialPivLu() function. This yields an object of the class PartialPivLU. Then, the PartialPivLU::solve() method is called to compute a solution.

As an example, suppose we want to solve the following system of linear equations:

\[ \begin{aligned} x + 2y + 3z &= 3 \\ 4x + 5y + 6z &= 3 \\ 7x + 8y + 10z &= 4. \end{aligned} \]

The following program solves this system:

#include <Eigen/Core>
#include <Eigen/LU>
#include <iostream>

using namespace std;
using namespace Eigen;

int main(int, char *[])
{
   Matrix3f A;
   Vector3f b;
   A << 1,2,3,  4,5,6,  7,8,10;
   b << 3, 3, 4;
   cout << "Here is the matrix A:" << endl << A << endl;
   cout << "Here is the vector b:" << endl << b << endl;
   Vector3f x = A.lu().solve(b);
   cout << "The solution is:" << endl << x << endl;
}
output: 
Here is the matrix A:
 1  2  3
 4  5  6
 7  8 10
Here is the vector b:
3
3
4
The solution is:
-2
1
1

There are many situations in which we want to solve the same system of equations with different right-hand sides. One possibility is to put the right-hand sides as columns in a matrix $ B $ and then solve the equation $ A X = B $. For instance, suppose that we want to solve the same system as before, but now we also need the solution of the same equations with 1 on the right-hand side. The following code computes the required solutions:

Matrix3f A(3,3);
A << 1,2,3,  4,5,6,  7,8,10;
Matrix<float,3,2> B;
B << 3,1, 3,1, 4,1;
Matrix<float,3,2> X;
X = A.lu().solve(B);
cout << "The solution with right-hand side (3,3,4) is:" << endl;
cout << X.col(0) << endl;
cout << "The solution with right-hand side (1,1,1) is:" << endl;
cout << X.col(1) << endl;
output: 
The solution with right-hand side (3,3,4) is:
-2
1
1
The solution with right-hand side (1,1,1) is:
-1
1
-2.98e-07

However, this is not always possible. Often, you only know the right-hand side of the second problem, and whether you want to solve it at all, after you solved the first problem. In such a case, it's best to save the LU decomposition and reuse it to solve the second problem. This is worth the effort because computing the LU decomposition is much more expensive than using it to solve the equation. Here is some code to illustrate the procedure. It uses the constructor PartialPivLU::PartialPivLU(const MatrixType&) to compute the LU decomposition.

Matrix3f A(3,3);
A << 1,2,3,  4,5,6,  7,8,10;
PartialPivLU<Matrix3f> luOfA(A); // compute LU decomposition of A
Vector3f b;
b << 3,3,4;
Vector3f x;
x = luOfA.solve(b);
cout << "The solution with right-hand side (3,3,4) is:" << endl;
cout << x << endl;
b << 1,1,1;
x = luOfA.solve(b);
cout << "The solution with right-hand side (1,1,1) is:" << endl;
cout << x << endl;
output: 
The solution with right-hand side (3,3,4) is:
-2
1
1
The solution with right-hand side (1,1,1) is:
-1
1
-2.38e-07

Warning: All this code presumes that the matrix $ A $ is invertible, so that the system $ A \mathbf{x} = \mathbf{b} $ has a unique solution. If the matrix $ A $ is not invertible, then the system $ A \mathbf{x} = \mathbf{b} $ has either zero or infinitely many solutions. In both cases, PartialPivLU::solve() will give nonsense results. For example, suppose that we want to solve the same system as above, but with the 10 in the last equation replaced by 9. Then the system of equations is inconsistent: adding the first and the third equation gives $ 8x + 10y + 12z = 7 $, which implies $ 4x + 5y + 6z = 3\frac12 $, in contradiction with the second equation. If we try to solve this inconsistent system with Eigen, we find:

Matrix3f A;
Vector3f b;
A << 1,2,3,  4,5,6,  7,8,9;
b << 3, 3, 4;
cout << "Here is the matrix A:" << endl << A << endl;
cout << "Here is the vector b:" << endl << b << endl;
Vector3f x;
x = A.lu().solve(b);
cout << "The solution is:" << endl << x << endl;
output: 
Here is the matrix A:
1 2 3
4 5 6
7 8 9
Here is the vector b:
3
3
4
The solution is:
nan
inf
-inf

The LU decomposition with full pivoting (class FullPivLU) and the singular value decomposition (class SVD) may be helpful in this case, as explained in the section Other solvers (for singular matrices and special cases) below.

See also:
LU module, MatrixBase::partialPivLu(), PartialPivLU::solve(), class PartialPivLU.

Cholesky decomposition

If the matrix $ A $ is symmetric positive definite, then the best method is to use a Cholesky decomposition: it is both faster and more accurate than the LU decomposition. Such positive definite matrices often arise when solving overdetermined problems. These are linear systems $ A \mathbf{x} = \mathbf{b} $ in which the matrix $ A $ has more rows than columns, so that there are more equations than unknowns. Typically, there is no vector $ \mathbf{x} $ which satisfies all the equation. Instead, we look for the least-square solution, that is, the vector $ \mathbf{x} $ for which $ \| A \mathbf{x} - \mathbf{b} \|_2 $ is minimal. You can find this vector by solving the equation $ A^T \! A \mathbf{x} = A^T \mathbf{b} $. If the matrix $ A $ has full rank, then $ A^T \! A $ is positive definite and thus you can use the Cholesky decomposition to solve this system and find the least-square solution to the original system $ A \mathbf{x} = \mathbf{b} $.

Eigen offers two different Cholesky decompositions: the LLT class provides a $ LL^T $ decomposition where L is a lower triangular matrix, and the LDLT class provides a $ LDL^T $ decomposition where L is lower triangular with unit diagonal and D is a diagonal matrix. The latter includes pivoting and avoids square roots; this makes the LDLT decomposition slightly more stable than the LLT decomposition. The LDLT class is able to handle both positive- and negative-definite matrices, but not indefinite matrices.

The API is the same as when using the LU decomposition.

#include <Eigen/Cholesky>
MatrixXf D = MatrixXf::Random(8,4);
MatrixXf A = D.transpose() * D;
VectorXf b = A * VectorXf::Random(4);
VectorXf x_llt = A.llt().solve(b);   // using a LLT factorization
VectorXf x_ldlt = A.ldlt().solve(b); // using a LDLT factorization

The LLT and LDLT classes also provide an in place API for the case where the value of the right hand-side $ b $ is not needed anymore.

A.llt().solveInPlace(b);

This code replaces the vector $ b $ by the result $ x $.

As before, you can reuse the factorization if you have to solve the same linear problem with different right-hand sides, e.g.:

// ...
LLT<MatrixXf> lltOfA(A);
lltOfA.solveInPlace(b0);
lltOfA.solveInPlace(b1);
// ...
See also:
Cholesky module, MatrixBase::llt(), MatrixBase::ldlt(), LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT.

Triangular solver

If the matrix $ A $ is triangular (upper or lower) and invertible (the coefficients of the diagonal are all not zero), then the problem can be solved directly using the TriangularView class. This is much faster than using an LU or Cholesky decomposition (in fact, the triangular solver is used when you solve a system using the LU or Cholesky decomposition). Here is an example:

Matrix3f A;
Vector3f b;
A << 1,2,3,  0,5,6,  0,0,10;
b << 3, 3, 4;
cout << "Here is the matrix A:" << endl << A << endl;
cout << "Here is the vector b:" << endl << b << endl;
Vector3f x = A.triangularView<Upper>().solve(b);
cout << "The solution is:" << endl << x << endl;
output: 
Here is the matrix A:
 1  2  3
 0  5  6
 0  0 10
Here is the vector b:
3
3
4
The solution is:
1.56
0.12
0.4

The MatrixBase::triangularView() function constructs an object of the class TriangularView, and TriangularView::solve() then solves the system. There is also an in place variant:

Matrix3f A;
Vector3f b;
A << 1,2,3,  0,5,6,  0,0,10;
b << 3, 3, 4;
A.triangularView<Upper>().solveInPlace(b);
cout << "The solution is:" << endl << b << endl;
output: 
The solution is:
1.56
0.12
0.4
See also:
MatrixBase::triangularView(), TriangularView::solve(), TriangularView::solveInPlace(), TriangularView class.

Direct inversion (for small matrices)

The solution of the system $ A \mathbf{x} = \mathbf{b} $ is given by $ \mathbf{x} = A^{-1} \mathbf{b} $. This suggests the following approach for solving the system: compute the matrix inverse and multiply that with the right-hand side. This is often not a good approach: using the LU decomposition with partial pivoting yields a more accurate algorithm that is usually just as fast or even faster. However, using the matrix inverse can be faster if the matrix $ A $ is small (4) and fixed size, though numerical stability problems may still remain. Here is an example of how you would write this in Eigen:

Matrix3f A;
Vector3f b;
A << 1,2,3,  4,5,6,  7,8,10;
b << 3, 3, 4;
Vector3f x = A.inverse() * b;
cout << "The solution is:" << endl << x << endl;
output: 
The solution is:
-2
1
1

Note that the function inverse() is defined in the LU module.

See also:
MatrixBase::inverse().

Other solvers (for singular matrices and special cases)

Finally, Eigen also offer solvers based on a singular value decomposition (SVD) or the LU decomposition with full pivoting. These have the same API as the solvers based on the LU decomposition with partial pivoting (PartialPivLU).

The solver based on the SVD uses the class SVD. It can handle singular matrices. Here is an example of its use:

#include <Eigen/SVD>
// ...
MatrixXf A = MatrixXf::Random(20,20);
VectorXf b = VectorXf::Random(20);
VectorXf x = A.svd().solve(b);
SVD<MatrixXf> svdOfA(A);
x = svdOfA.solve(b);

LU decomposition with full pivoting has better numerical stability than LU decomposition with partial pivoting. It is defined in the class FullPivLU. The solver can also handle singular matrices.

#include <Eigen/LU>
// ...
MatrixXf A = MatrixXf::Random(20,20);
VectorXf b = VectorXf::Random(20);
VectorXf x = A.lu().solve(b);
FullPivLU<MatrixXf> luOfA(A);
x = luOfA.solve(b);

See the section LU below.

See also:
class SVD, SVD::solve(), SVD module, class FullPivLU, LU::solve(), LU module.

top

LU

Eigen provides a rank-revealing LU decomposition with full pivoting, which has very good numerical stability.

You can obtain the LU decomposition of a matrix by calling lu() , which is the easiest way if you're going to use the LU decomposition only once, as in 

#include <Eigen/LU>
MatrixXf A = MatrixXf::Random(20,20);
VectorXf b = VectorXf::Random(20);
VectorXf x = A.lu().solve(b);

Alternatively, you can construct a named LU decomposition, which allows you to reuse it for more than one operation: 

#include <Eigen/LU>
MatrixXf A = MatrixXf::Random(20,20);
Eigen::FullPivLU<MatrixXf> lu(A);
cout << "The rank of A is" << lu.rank() << endl;
if(lu.isInvertible()) {
  cout << "A is invertible, its inverse is:" << endl << lu.inverse() << endl;
}
else {
  cout << "Here's a matrix whose columns form a basis of the kernel a.k.a. nullspace of A:"
       << endl << lu.kernel() << endl;
}
See also:
LU module, LU::solve(), class FullPivLU

top

Cholesky

todo

See also:
Cholesky module, LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT

top

QR

todo

See also:
QR module, class QR

top

Eigen value problems

todo

See also:
class SelfAdjointEigenSolver, class EigenSolver
Generated on Fri Mar 19 06:59:32 2010 for Eigen by  doxygen 1.6.3