Example:Output:
\include TutorialLinAlgExSolveColPivHouseholderQR.cpp \verbinclude TutorialLinAlgExSolveColPivHouseholderQR.out
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: \code ColPivHouseholderQR dec(A); Vector3f x = dec.solve(b); \endcode 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)
Accuracy
PartialPivLU partialPivLu() Invertible ++ ++ +
FullPivLU fullPivLu() None - - - +++
HouseholderQR householderQr() None ++ ++ +
ColPivHouseholderQR colPivHouseholderQr() None + - +++
FullPivHouseholderQR fullPivHouseholderQr() None - - - +++
CompleteOrthogonalDecomposition completeOrthogonalDecomposition() None + - +++
LLT llt() Positive definite +++ +++ +
LDLT ldlt() Positive or negative
semidefinite
+++ + ++
BDCSVD bdcSvd() None - - +++
JacobiSVD jacobiSvd() None - - - - +++
To get an overview of the true relative speed of the different decompositions, check this \link DenseDecompositionBenchmark benchmark \endlink. All of these decompositions offer a solve() method that works as in the above example. If you know more about the properties of your matrix, you can use the above table to select the best method. For example, a good choice for solving linear systems with a non-symmetric matrix of full rank is PartialPivLU. If you know that your matrix is also symmetric and positive definite, the above table says that a very good choice is the LLT or LDLT decomposition. Here's an example, also demonstrating that using a general matrix (not a vector) as right hand side is possible:
Example:Output:
\include TutorialLinAlgExSolveLDLT.cpp \verbinclude TutorialLinAlgExSolveLDLT.out
For a \ref TopicLinearAlgebraDecompositions "much more complete table" comparing all decompositions supported by Eigen (notice that Eigen supports many other decompositions), see our special page on \ref TopicLinearAlgebraDecompositions "this topic". \section TutorialLinAlgLeastsquares Least squares solving The most general and accurate method to solve under- or over-determined linear systems in the least squares sense, is the SVD decomposition. Eigen provides two implementations. The recommended one is the BDCSVD class, which scales well for large problems and automatically falls back to the JacobiSVD class for smaller problems. For both classes, their solve() method solved the linear system in the least-squares sense. Here is an example:
Example:Output:
\include TutorialLinAlgSVDSolve.cpp \verbinclude TutorialLinAlgSVDSolve.out
An alternative to the SVD, which is usually faster and about as accurate, is CompleteOrthogonalDecomposition. Again, if you know more about the problem, the table above contains methods that are potentially faster. If your matrix is full rank, HouseHolderQR is the method of choice. If your matrix is full rank and well conditioned, using the Cholesky decomposition (LLT) on the matrix of the normal equations can be faster still. Our page on \link LeastSquares least squares solving \endlink has more details. \section TutorialLinAlgSolutionExists 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:
\include TutorialLinAlgExComputeSolveError.cpp \verbinclude TutorialLinAlgExComputeSolveError.out
\section TutorialLinAlgEigensolving Computing eigenvalues and eigenvectors You need an eigendecomposition here, see available such decompositions on \ref TopicLinearAlgebraDecompositions "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: