template<typename _MatrixType, int _UpLo>
Eigen::LLT class

Standard Cholesky decomposition (LL^T) of a matrix and associated features.

Template parameters
_MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
_UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper. The other triangular part won't be read.

Contents

This class performs a LL^T Cholesky decomposition of a symmetric, positive definite matrix A such that A = LL^* = U^*U, where L is lower triangular.

While the Cholesky decomposition is particularly useful to solve selfadjoint problems like D^*D x = b, for that purpose, we recommend the Cholesky decomposition without square root which is more stable and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other situations like generalised eigen problems with hermitian matrices.

Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices, use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations has a solution.

Example:

MatrixXd A(3,3);
A << 4,-1,2, -1,6,0, 2,0,5;
cout << "The matrix A is" << endl << A << endl;

LLT<MatrixXd> lltOfA(A); // compute the Cholesky decomposition of A
MatrixXd L = lltOfA.matrixL(); // retrieve factor L  in the decomposition
// The previous two lines can also be written as "L = A.llt().matrixL()"

cout << "The Cholesky factor L is" << endl << L << endl;
cout << "To check this, let us compute L * L.transpose()" << endl;
cout << L * L.transpose() << endl;
cout << "This should equal the matrix A" << endl;

Output:

The matrix A is
 4 -1  2
-1  6  0
 2  0  5
The Cholesky factor L is
    2     0     0
 -0.5   2.4     0
    1 0.209  1.99
To check this, let us compute L * L.transpose()
 4 -1  2
-1  6  0
 2  0  5
This should equal the matrix A

Performance: for best performance, it is recommended to use a column-major storage format with the Lower triangular part (the default), or, equivalently, a row-major storage format with the Upper triangular part. Otherwise, you might get a 20% slowdown for the full factorization step, and rank-updates can be up to 3 times slower.

This class supports the inplace decomposition mechanism.

Note that during the decomposition, only the lower (or upper, as defined by _UpLo) triangular part of A is considered. Therefore, the strict lower part does not have to store correct values.

Base classes

template<typename Derived>
class SolverBase
A base class for matrix decomposition and solvers.

Constructors, destructors, conversion operators

LLT()
Default Constructor.
LLT(Index size) explicit
Default Constructor with memory preallocation.
template<typename InputType>
LLT(EigenBase<InputType>& matrix) explicit
Constructs a LLT factorization from a given matrix.

Public functions

auto adjoint() const -> const LLT&
template<typename InputType>
auto compute(const EigenBase<InputType>& a) -> LLT<MatrixType, _UpLo>&
auto info() const -> ComputationInfo
Reports whether previous computation was successful.
auto matrixL() const -> Traits::MatrixL
auto matrixLLT() const -> const MatrixType&
auto matrixU() const -> Traits::MatrixU
template<typename VectorType>
auto rankUpdate(const VectorType& v, const RealScalar& sigma) -> LLT<_MatrixType, _UpLo>&
auto rcond() const -> RealScalar
auto reconstructedMatrix() const -> MatrixType
template<typename Rhs>
auto solve(const MatrixBase<Rhs>& b) const -> const Solve<LLT, Rhs>

Function documentation

template<typename _MatrixType, int _UpLo>
Eigen::LLT<_MatrixType, _UpLo>::LLT()

Default Constructor.

The default constructor is useful in cases in which the user intends to perform decompositions via LLT::compute(const MatrixType&).

template<typename _MatrixType, int _UpLo>
Eigen::LLT<_MatrixType, _UpLo>::LLT(Index size) explicit

Default Constructor with memory preallocation.

Like the default constructor but with preallocation of the internal data according to the specified problem size.

template<typename _MatrixType, int _UpLo> template<typename InputType>
Eigen::LLT<_MatrixType, _UpLo>::LLT(EigenBase<InputType>& matrix) explicit

Constructs a LLT factorization from a given matrix.

This overloaded constructor is provided for inplace decomposition when MatrixType is a Eigen::Ref.

template<typename _MatrixType, int _UpLo>
const LLT& Eigen::LLT<_MatrixType, _UpLo>::adjoint() const

Returns the adjoint of *this, that is, a const reference to the decomposition itself as the underlying matrix is self-adjoint.

This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as: x = decomposition.adjoint().solve(b)

template<typename _MatrixType, int _UpLo> template<typename InputType>
LLT<MatrixType, _UpLo>& Eigen::LLT<_MatrixType, _UpLo>::compute(const EigenBase<InputType>& a)

Returns a reference to *this

Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of matrix

Example:

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

using namespace std;
using namespace Eigen;

int main()
{
   Matrix2f A, b;
   LLT<Matrix2f> llt;
   A << 2, -1, -1, 3;
   b << 1, 2, 3, 1;
   cout << "Here is the matrix A:\n" << A << endl;
   cout << "Here is the right hand side b:\n" << b << endl;
   cout << "Computing LLT decomposition..." << endl;
   llt.compute(A);
   cout << "The solution is:\n" << llt.solve(b) << endl;
   A(1,1)++;
   cout << "The matrix A is now:\n" << A << endl;
   cout << "Computing LLT decomposition..." << endl;
   llt.compute(A);
   cout << "The solution is now:\n" << llt.solve(b) << endl;
}

Output:

Here is the matrix A:
 2 -1
-1  3
Here is the right hand side b:
1 2
3 1
Computing LLT decomposition...
The solution is:
1.2 1.4
1.4 0.8
The matrix A is now:
 2 -1
-1  4
Computing LLT decomposition...
The solution is now:
    1  1.29
    1 0.571

template<typename _MatrixType, int _UpLo>
ComputationInfo Eigen::LLT<_MatrixType, _UpLo>::info() const

Reports whether previous computation was successful.

Returns Success if computation was successful, NumericalIssue if the matrix.appears not to be positive definite.

template<typename _MatrixType, int _UpLo>
Traits::MatrixL Eigen::LLT<_MatrixType, _UpLo>::matrixL() const

Returns a view of the lower triangular matrix L

template<typename _MatrixType, int _UpLo>
const MatrixType& Eigen::LLT<_MatrixType, _UpLo>::matrixLLT() const

Returns the LLT decomposition matrix

TODO: document the storage layout

template<typename _MatrixType, int _UpLo>
Traits::MatrixU Eigen::LLT<_MatrixType, _UpLo>::matrixU() const

Returns a view of the upper triangular matrix U

template<typename _MatrixType, int _UpLo> template<typename VectorType>
LLT<_MatrixType, _UpLo>& Eigen::LLT<_MatrixType, _UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma)

Performs a rank one update (or dowdate) of the current decomposition. If A = LL^* before the rank one update, then after it we have LL^* = A + sigma * v v^* where v must be a vector of same dimension.

template<typename _MatrixType, int _UpLo>
RealScalar Eigen::LLT<_MatrixType, _UpLo>::rcond() const

Returns an estimate of the reciprocal condition number of the matrix of which *this is the Cholesky decomposition.

template<typename _MatrixType, int _UpLo>
MatrixType Eigen::LLT<_MatrixType, _UpLo>::reconstructedMatrix() const

Returns the matrix represented by the decomposition, i.e., it returns the product: L L^*. This function is provided for debug purpose.

template<typename _MatrixType, int _UpLo> template<typename Rhs>
const Solve<LLT, Rhs> Eigen::LLT<_MatrixType, _UpLo>::solve(const MatrixBase<Rhs>& b) const

Returns the solution x of $ A x = b $ using the current decomposition of A.

Since this LLT class assumes anyway that the matrix A is invertible, the solution theoretically exists and is unique regardless of b.

Example:

typedef Matrix<float,Dynamic,2> DataMatrix;
// let's generate some samples on the 3D plane of equation z = 2x+3y (with some noise)
DataMatrix samples = DataMatrix::Random(12,2);
VectorXf elevations = 2*samples.col(0) + 3*samples.col(1) + VectorXf::Random(12)*0.1;
// and let's solve samples * [x y]^T = elevations in least square sense:
Matrix<float,2,1> xy
 = (samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations));
cout << xy << endl;

Output:

2.02
2.97