Linear algebra and decompositions

This page explains how to solve linear systems, compute various decompositions such as LU, QR, SVD, eigendecompositions... After reading this page, don't miss our catalogue of dense matrix decompositions.

Basic linear solving

The problem: You have a system of equations, that you have written as a single matrix equation

\[ Ax \: = \: b \]

Where A and b are matrices (b could be a vector, as a special case). You want to find a solution x.

The solution: You can choose between various decompositions, depending on the properties of your matrix A, and depending on whether you favor speed or accuracy. However, let's start with an example that works in all cases, and is a good compromise:

Example:

Output:

// SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 #include <iostream> #include <Eigen/Dense> int main() { Eigen::Matrix3f A; Eigen::Vector3f b; A << 1, 2, 3, 4, 5, 6, 7, 8, 10; b << 3, 3, 4; std::cout << "Here is the matrix A:\n" << A << std::endl; std::cout << "Here is the vector b:\n" << b << std::endl; Eigen::Vector3f x = A.colPivHouseholderQr().solve(b); std::cout << "The solution is:\n" << x << std::endl; }

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

In this example, the colPivHouseholderQr() method returns an object of class ColPivHouseholderQR. Since here the matrix is of type Matrix3f, this line could have been replaced by:

ColPivHouseholderQR<Matrix3f> dec(A);
Vector3f x = dec.solve(b);

Here, ColPivHouseholderQR is a QR decomposition with column pivoting. It's a good compromise for this tutorial, as it works for all matrices while being quite fast. Here is a table of some other decompositions that you can choose from, depending on your matrix, the problem you are trying to solve, and the trade-off you want to make:

Decomposition	Method	Requirements on the matrix	Speed (small-to-medium)	Speed (large)	Robustness*
LLT	llt()	Positive definite	+++	+++	+
LDLT	ldlt()	Positive or negative semidefinite	+++	+	++
BunchKaufman	bunchKaufman()	Self-adjoint	+++	++	++
PartialPivLU	partialPivLu()	Invertible	++	++	+
HouseholderQR	householderQr()	None	++	++	+
ColPivHouseholderQR	colPivHouseholderQr()	None	+	-	+++
CompleteOrthogonalDecomposition	completeOrthogonalDecomposition()	None	+	-	+++
RandCompleteOrthogonalDecomposition	randCompleteOrthogonalDecomposition()	None	-	+	+++
BDCSVD	bdcSvd()	None	-	-	+++
JacobiSVD	jacobiSvd()	None	-	- - -	+++

* The Robustness column indicates how well the decomposition handles ill-conditioned or rank-deficient matrices. All decompositions give excellent accuracy when their requirements on the matrix are met and the problem is well-conditioned.

To get an overview of the true relative speed of the different decompositions, check this benchmark .

All of these decompositions offer a solve() method that works as in the above example.

Practical recommendations:

• If your matrix is symmetric positive definite, use LLT. It is the fastest and is perfectly accurate for this class of problems. If your matrix is only positive or negative semidefinite, use LDLT, whose symmetric pivoting keeps the factorization stable in the presence of tiny or zero pivots. LDLT produces a purely diagonal D and so cannot handle indefinite matrices; for a self-adjoint matrix that is indefinite, use BunchKaufman, which produces a block-diagonal D at LDLT-like speed and also reveals the inertia.
• For a general invertible matrix, PartialPivLU is the best choice. It is fast (uses cache-friendly blocking) and reliable for the vast majority of problems.
• For least squares problems (over- or under-determined systems), CompleteOrthogonalDecomposition is the recommended default. Like the SVD, it robustly computes the minimum-norm solution for rank-deficient and under-determined problems, but at the cost of a QR decomposition rather than an SVD. Use ColPivHouseholderQR if you only need least squares for full-rank overdetermined systems and don't need the minimum-norm property.
• For large least squares problems, consider RandCompleteOrthogonalDecomposition. It computes the same kind of factorization as CompleteOrthogonalDecomposition, but its rank-revealing QR step selects pivots from a small random sketch, which lets nearly all of the work be done by level-3 BLAS kernels instead of the memory-bound column scan used by classical column pivoting. On small matrices the sketch overhead dominates and plain CompleteOrthogonalDecomposition is preferable.
• SVD decompositions (BDCSVD, JacobiSVD) are the most robust but also the slowest. Use these when you need singular values/vectors, not just the solution.
• HouseholderQR is the fastest option for full-rank least squares problems, but it does not reveal rank and cannot compute minimum-norm solutions for rank-deficient problems.
• FullPivLU and FullPivHouseholderQR use complete pivoting, which rules out blocking and so makes them substantially slower, increasingly so as the matrix grows. What it buys – a smaller growth factor for LU, a more reliable rank decision for QR – rarely repays that cost on its own, so PartialPivLU and ColPivHouseholderQR are the better default; LAPACK made the same call, offering no complete-pivoting LU or QR driver. Where complete pivoting does earn its cost is rank determination, since column pivoting is reliable in practice but not guaranteed. If your application depends on getting the rank right, the full-pivoting classes are a reasonable choice, especially for small matrices, where the lost blocking costs comparatively little.

Here's an example showing the use of LLT for a symmetric positive definite system, also demonstrating that using a general matrix (not a vector) as right hand side is possible:

Example:

Output:

// SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 #include <iostream> #include <Eigen/Dense> int main() { Eigen::Matrix2f A, b; A << 2, -1, -1, 3; b << 1, 2, 3, 1; std::cout << "Here is the matrix A:\n" << A << std::endl; std::cout << "Here is the right hand side b:\n" << b << std::endl; Eigen::Matrix2f x = A.ldlt().solve(b); std::cout << "The solution is:\n" << x << std::endl; }

Here is the matrix A: 2 -1 -1 3 Here is the right hand side b: 1 2 3 1 The solution is: 1.2 1.4 1.4 0.8

For a much more complete table comparing all decompositions supported by Eigen (notice that Eigen supports many other decompositions), see our special page on this topic.

Least squares solving

The recommended method to solve under- or over-determined linear systems in the least squares sense is CompleteOrthogonalDecomposition. Like the SVD, it robustly computes the minimum-norm least squares solution, correctly handling rank-deficient and under-determined problems, but it is significantly faster since it is based on a rank-revealing QR decomposition rather than a full SVD. For large problems, RandCompleteOrthogonalDecomposition computes the same kind of factorization with a randomized blocked rank-revealing QR that spends nearly all of its time in level-3 BLAS kernels.

If you also need the singular values or vectors themselves (not just the least squares solution), use BDCSVD, which scales well for large problems and automatically falls back to JacobiSVD for smaller ones.

Here is an example using the SVD:

Example:

Output:

// SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 #include <iostream> #include <Eigen/Dense> int main() { Eigen::MatrixXf A = Eigen::MatrixXf::Random(3, 2); std::cout << "Here is the matrix A:\n" << A << std::endl; Eigen::VectorXf b = Eigen::VectorXf::Random(3); std::cout << "Here is the right hand side b:\n" << b << std::endl; std::cout << "The least-squares solution is:\n" << A.bdcSvd<Eigen::ComputeThinU | Eigen::ComputeThinV>().solve(b) << std::endl; }

Here is the matrix A: 0.68 0.597 -0.211 0.823 0.566 -0.605 Here is the right hand side b: -0.33 0.536 -0.444 The least-squares solution is: -0.67 0.314

If you know more about the problem, faster methods are available. If your matrix is full rank, HouseholderQR is the fastest method. If your matrix is full rank and well conditioned, using the Cholesky decomposition (LLT) on the normal equations can be faster still. Our page on least squares solving has more details.

Checking if a matrix is singular

Only you know what error margin you want to allow for a solution to be considered valid. So Eigen lets you do this computation for yourself, if you want to, as in this example:

Example:

Output:

// SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 #include <iostream> #include <Eigen/Dense> using Eigen::MatrixXd; int main() { MatrixXd A = MatrixXd::Random(100, 100); MatrixXd b = MatrixXd::Random(100, 50); MatrixXd x = A.fullPivLu().solve(b); double relative_error = (A * x - b).norm() / b.norm(); // norm() is L2 norm std::cout << "The relative error is:\n" << relative_error << std::endl; }

The relative error is: 2.24987e-14

Computing eigenvalues and eigenvectors

You need an eigendecomposition here, see available such decompositions on this page. Make sure to check if your matrix is self-adjoint, as is often the case in these problems. Here's an example using SelfAdjointEigenSolver, it could easily be adapted to general matrices using EigenSolver or ComplexEigenSolver.

The computation of eigenvalues and eigenvectors does not necessarily converge, but such failure to converge is very rare. The call to info() is to check for this possibility.

Example:

Output:

// SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 #include <iostream> #include <Eigen/Dense> int main() { Eigen::Matrix2f A; A << 1, 2, 2, 3; std::cout << "Here is the matrix A:\n" << A << std::endl; Eigen::SelfAdjointEigenSolver<Eigen::Matrix2f> eigensolver(A); if (eigensolver.info() != Eigen::Success) abort(); std::cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << std::endl; std::cout << "Here's a matrix whose columns are eigenvectors of A \n" << "corresponding to these eigenvalues:\n" << eigensolver.eigenvectors() << std::endl; }

Here is the matrix A: 1 2 2 3 The eigenvalues of A are: -0.236 4.24 Here's a matrix whose columns are eigenvectors of A corresponding to these eigenvalues: -0.851 -0.526 0.526 -0.851

Computing inverse and determinant

First of all, make sure that you really want this. While inverse and determinant are fundamental mathematical concepts, in numerical linear algebra they are not as useful as in pure mathematics. Inverse computations are often advantageously replaced by solve() operations, and the determinant is often not a good way of checking if a matrix is invertible.

However, for very small matrices, the above may not be true, and inverse and determinant can be very useful.

While certain decompositions, such as PartialPivLU and FullPivLU, offer inverse() and determinant() methods, you can also call inverse() and determinant() directly on a matrix. If your matrix is of a very small fixed size (at most 4x4) this allows Eigen to avoid performing a LU decomposition, and instead use formulas that are more efficient on such small matrices.

Here is an example:

Example:	Output:
// SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 #include <iostream> #include <Eigen/Dense> int main() { Eigen::Matrix3f A; A << 1, 2, 1, 2, 1, 0, -1, 1, 2; std::cout << "Here is the matrix A:\n" << A << std::endl; std::cout << "The determinant of A is " << A.determinant() << std::endl; std::cout << "The inverse of A is:\n" << A.inverse() << std::endl; }	Here is the matrix A: 1 2 1 2 1 0 -1 1 2 The determinant of A is -3 The inverse of A is: -0.667 1 0.333 1.33 -1 -0.667 -1 1 1

Separating the computation from the construction

In the above examples, the decomposition was computed at the same time that the decomposition object was constructed. There are however situations where you might want to separate these two things, for example if you don't know, at the time of the construction, the matrix that you will want to decompose; or if you want to reuse an existing decomposition object.

What makes this possible is that:

• all decompositions have a default constructor,
• all decompositions have a compute(matrix) method that does the computation, and that may be called again on an already-computed decomposition, reinitializing it.

For example:

Example:

Output:

// SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 #include <iostream> #include <Eigen/Dense> int main() { Eigen::Matrix2f A, b; Eigen::LLT<Eigen::Matrix2f> llt; A << 2, -1, -1, 3; b << 1, 2, 3, 1; std::cout << "Here is the matrix A:\n" << A << std::endl; std::cout << "Here is the right hand side b:\n" << b << std::endl; std::cout << "Computing LLT decomposition..." << std::endl; llt.compute(A); std::cout << "The solution is:\n" << llt.solve(b) << std::endl; A(1, 1)++; std::cout << "The matrix A is now:\n" << A << std::endl; std::cout << "Computing LLT decomposition..." << std::endl; llt.compute(A); std::cout << "The solution is now:\n" << llt.solve(b) << std::endl; }

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

Finally, you can tell the decomposition constructor to preallocate storage for decomposing matrices of a given size, so that when you subsequently decompose such matrices, no dynamic memory allocation is performed (of course, if you are using fixed-size matrices, no dynamic memory allocation happens at all). This is done by just passing the size to the decomposition constructor, as in this example:

HouseholderQR<MatrixXf> qr(50,50);
MatrixXf A = MatrixXf::Random(50,50);
qr.compute(A); // no dynamic memory allocation

Rank-revealing decompositions

Certain decompositions are rank-revealing, i.e. are able to compute the rank of a matrix. These are typically also the decompositions that behave best in the face of a non-full-rank matrix (which in the square case means a singular matrix). On this table you can see for all our decompositions whether they are rank-revealing or not.

Rank-revealing decompositions offer at least a rank() method. They can also offer convenience methods such as isInvertible(), and FullPivLU is the only one that additionally computes the kernel (null-space) and the image (column-space) of the matrix, through kernel() and image(). ColPivHouseholderQR and CompleteOrthogonalDecomposition are faster rank-revealing factorizations providing rank() and isInvertible() but no null-space basis. CompleteOrthogonalDecomposition and RandCompleteOrthogonalDecomposition additionally offer pseudoInverse(); they are the only decompositions in Eigen that do, so they are also the cheapest route to a Moore-Penrose pseudo-inverse, which would otherwise have to be assembled by hand from an SVD. Note that solve() should still be preferred over multiplying by an explicitly formed pseudo-inverse. Here is an example using FullPivLU:

Example:

Output:

// SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 #include <iostream> #include <Eigen/Dense> int main() { Eigen::Matrix3f A; A << 1, 2, 5, 2, 1, 4, 3, 0, 3; std::cout << "Here is the matrix A:\n" << A << std::endl; Eigen::FullPivLU<Eigen::Matrix3f> lu_decomp(A); std::cout << "The rank of A is " << lu_decomp.rank() << std::endl; std::cout << "Here is a matrix whose columns form a basis of the null-space of A:\n" << lu_decomp.kernel() << std::endl; std::cout << "Here is a matrix whose columns form a basis of the column-space of A:\n" << lu_decomp.image(A) << std::endl; // yes, have to pass the original A }

Here is the matrix A: 1 2 5 2 1 4 3 0 3 The rank of A is 2 Here is a matrix whose columns form a basis of the null-space of A: 0.5 1 -0.5 Here is a matrix whose columns form a basis of the column-space of A: 5 1 4 2 3 3

Of course, any rank computation depends on the choice of an arbitrary threshold, since practically no floating-point matrix is exactly rank-deficient. Eigen picks a sensible default threshold, which depends on the decomposition but is typically the diagonal size times machine epsilon. While this is the best default we could pick, only you know what is the right threshold for your application. You can set this by calling setThreshold() on your decomposition object before calling rank() or any other method that needs to use such a threshold. The decomposition itself, i.e. the compute() method, is independent of the threshold. You don't need to recompute the decomposition after you've changed the threshold.

Example:

Output:

// SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 #include <iostream> #include <Eigen/Dense> int main() { Eigen::Matrix2d A; A << 2, 1, 2, 0.9999999999; Eigen::FullPivLU<Eigen::Matrix2d> lu(A); std::cout << "By default, the rank of A is found to be " << lu.rank() << std::endl; lu.setThreshold(1e-5); std::cout << "With threshold 1e-5, the rank of A is found to be " << lu.rank() << std::endl; }

By default, the rank of A is found to be 2 With threshold 1e-5, the rank of A is found to be 1