Tutorial 3/3 - Advanced linear algebra

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

Table of contents

Solving linear problems

This part of the tutorial focuses on solving linear problem of the form $ A \mathbf{x} = b $, where both $ A $ and $ b $ are known, and $ x $ is the unknown. Moreover, $ A $ assumed to be a squared matrix. Of course, 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} B $. Eigen offers various algorithms to this problem, and its choice mainly depends on the nature of the matrix $ A $, such as its shape, size and numerical properties.

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 MatrixBase::solveTriangular(), or better, MatrixBase::solveTriangularInPlace(). Here is an example:

Matrix3d m = Matrix3d::Zero();
m.part<Eigen::UpperTriangular>().setOnes();
cout << "Here is the matrix m:" << endl << m << endl;
Matrix3d n = Matrix3d::Ones();
n.part<Eigen::LowerTriangular>() *= 2;
cout << "Here is the matrix n:" << endl << n << endl;
cout << "And now here is m.inverse()*n, taking advantage of the fact that"
        " m is upper-triangular:" << endl
     << m.marked<Eigen::UpperTriangular>().solveTriangular(n);
output:
Here is the matrix m:
1 1 1
0 1 1
0 0 1
Here is the matrix n:
2 1 1
2 2 1
2 2 2
And now here is m.inverse()*n, taking advantage of the fact that m is upper-triangular:
 0 -1  0
 0  0 -1
 2  2  2

See MatrixBase::solveTriangular() for more details.

Direct inversion (for small matrices)

If the matrix $ A $ is small ( $ \leq 4 $) and invertible, then a good approach is to directly compute the inverse of the matrix $ A $, and then obtain the solution $ x $ by $ \mathbf{x} = A^{-1} b $. With Eigen, this can be implemented like this:

#include <Eigen/LU>
Matrix4f A = Matrix4f::Random();
Vector4f b = Vector4f::Random();
Vector4f x = A.inverse() * b;

Note that the function inverse() is defined in the LU module. See MatrixBase::inverse() for more details.

Cholesky (for positive definite matrices)

If the matrix $ A $ is positive definite, then the best method is to use a Cholesky decomposition. Such positive definite matrices often arise when solving overdetermined problems in a least square sense (see below). Eigen offers two different Cholesky decompositions: a $ LL^T $ decomposition where L is a lower triangular matrix, and a $ LDL^T $ decomposition where L is lower triangular with unit diagonal and D is a diagonal matrix. The latter avoids square roots and is therefore slightly more stable than the former one.

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

when the value of the right hand side $ b $ is not needed anymore, then it is faster to use the in place API, e.g.:

A.llt().solveInPlace(b);

In that case the value of $ b $ is replaced by the result $ x $.

If the linear problem has to solved for various right hand side $ b_i $, but same matrix $ A $, then it is highly recommended to perform the decomposition of $ A $ only once, e.g.:

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

LU decomposition (for most cases)

If the matrix $ A $ does not fit in any of the previous categories, or if you are unsure about the numerical stability of your problem, then you can use the LU solver based on a decomposition of the same name : see the section LU below. Actually, Eigen's LU module does not implement a standard LU decomposition, but rather a so-called LU decomposition with full pivoting and rank update which has much better numerical stability. The API of the LU solver is the same than the Cholesky one, except that there is no in place variant:

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

Again, the LU decomposition can be stored to be reused or to perform other kernel operations:

// ...
LU<MatrixXf> luOfA(A);
luOfA.solve(b, &x);
// ...

See the section LU below.

See also:
class LU, LU::solve(), LU module

SVD solver (for singular matrices and special cases)

Finally, Eigen also offer a solver based on a singular value decomposition (SVD). Again, the API is the same than with Cholesky or LU:

#include <Eigen/SVD>
MatrixXf A = MatrixXf::Random(20,20);
VectorXf b = VectorXf::Random(20);
VectorXf x;
A.svd().solve(b, &x);
SVD<MatrixXf> svdOfA(A);
svdOfA.solve(b, &x);
See also:
class SVD, SVD::solve(), SVD 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, &x);

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::LUDecomposition<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 LU

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