namespace Eigen { /** \page TutorialAdvancedLinearAlgebra Tutorial 3/4 - Advanced linear algebra \ingroup Tutorial
\ref index "Overview" | \ref TutorialCore "Core features" | \ref TutorialGeometry "Geometry" | \b Advanced \b linear \b algebra | \ref TutorialSparse "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. \b Table \b of \b contents - \ref TutorialAdvSolvers - \ref TutorialAdvLU - \ref TutorialAdvCholesky - \ref TutorialAdvQR - \ref TutorialAdvEigenProblems \section TutorialAdvSolvers Solving linear problems This part of the tutorial focuses on solving systems of linear equations. Such systems can be written in the form \f$ A \mathbf{x} = \mathbf{b} \f$, where both \f$ A \f$ and \f$ \mathbf{b} \f$ are known, and \f$ \mathbf{x} \f$ is the unknown. Moreover, \f$ A \f$ is assumed to be a square matrix. The equation \f$ A \mathbf{x} = \mathbf{b} \f$ has a unique solution if \f$ A \f$ is invertible. If the matrix is not invertible, then the equation may have no or infinitely many solutions. All solvers assume that \f$ A \f$ is invertible, unless noted otherwise. Eigen offers various algorithms to this problem. The choice of algorithm mainly depends on the nature of the matrix \f$ A \f$, such as its shape, size and numerical properties. - The \ref TutorialAdvSolvers_LU "LU decomposition" (with partial pivoting) is a general-purpose algorithm which works for most problems. - Use the \ref TutorialAdvSolvers_Cholesky "Cholesky decomposition" if the matrix \f$ A \f$ is positive definite. - Use a special \ref TutorialAdvSolvers_Triangular "triangular solver" if the matrix \f$ A \f$ is upper or lower triangular. - Use of the \ref TutorialAdvSolvers_Inverse "matrix inverse" is not recommended in general, but may be appropriate in special cases, for instance if you want to solve several systems with the same matrix \f$ A \f$ and that matrix is small. - \ref TutorialAdvSolvers_Misc "Other solvers" (%LU decomposition with full pivoting, the singular value decomposition) are provided for special cases, such as when \f$ A \f$ is not invertible. The methods described here can be used whenever an expression involve the product of an inverse matrix with a vector or another matrix: \f$ A^{-1} \mathbf{v} \f$ or \f$ A^{-1} B \f$. \subsection TutorialAdvSolvers_LU LU decomposition (with partial pivoting) This is a general-purpose algorithm which performs well in most cases (provided the matrix \f$ A \f$ 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: \f[ \begin{aligned} x + 2y + 3z &= 3 \\ 4x + 5y + 6z &= 3 \\ 7x + 8y + 10z &= 4. \end{aligned} \f] The following program solves this system:
\include Tutorial_PartialLU_solve.cpp output: \include Tutorial_PartialLU_solve.out
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 \f$ B \f$ and then solve the equation \f$ A X = B \f$. 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:
\include Tutorial_solve_multiple_rhs.cpp output: \include Tutorial_solve_multiple_rhs.out
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.
\include Tutorial_solve_reuse_decomposition.cpp output: \include Tutorial_solve_reuse_decomposition.out
\b Warning: All this code presumes that the matrix \f$ A \f$ is invertible, so that the system \f$ A \mathbf{x} = \mathbf{b} \f$ has a unique solution. If the matrix \f$ A \f$ is not invertible, then the system \f$ A \mathbf{x} = \mathbf{b} \f$ 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 \f$ 8x + 10y + 12z = 7 \f$, which implies \f$ 4x + 5y + 6z = 3\frac12 \f$, in contradiction with the second equation. If we try to solve this inconsistent system with Eigen, we find:
\include Tutorial_solve_singular.cpp output: \include Tutorial_solve_singular.out
The %LU decomposition with \b full pivoting (class FullPivLU) and the singular value decomposition (class SVD) may be helpful in this case, as explained in the section \ref TutorialAdvSolvers_Misc below. \sa LU_Module, MatrixBase::partialPivLu(), PartialPivLU::solve(), class PartialPivLU. \subsection TutorialAdvSolvers_Cholesky Cholesky decomposition If the matrix \f$ A \f$ is \b symmetric \b positive \b 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 \f$ A \mathbf{x} = \mathbf{b} \f$ in which the matrix \f$ A \f$ has more rows than columns, so that there are more equations than unknowns. Typically, there is no vector \f$ \mathbf{x} \f$ which satisfies all the equation. Instead, we look for the least-square solution, that is, the vector \f$ \mathbf{x} \f$ for which \f$ \| A \mathbf{x} - \mathbf{b} \|_2 \f$ is minimal. You can find this vector by solving the equation \f$ A^T \! A \mathbf{x} = A^T \mathbf{b} \f$. If the matrix \f$ A \f$ has full rank, then \f$ A^T \! A \f$ 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 \f$ A \mathbf{x} = \mathbf{b} \f$. Eigen offers two different Cholesky decompositions: the LLT class provides a \f$ LL^T \f$ decomposition where L is a lower triangular matrix, and the LDLT class provides a \f$ LDL^T \f$ 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. \code #include 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 \endcode The LLT and LDLT classes also provide an \em in \em place API for the case where the value of the right hand-side \f$ b \f$ is not needed anymore. \code A.llt().solveInPlace(b); \endcode This code replaces the vector \f$ b \f$ by the result \f$ x \f$. As before, you can reuse the factorization if you have to solve the same linear problem with different right-hand sides, e.g.: \code // ... LLT lltOfA(A); lltOfA.solveInPlace(b0); lltOfA.solveInPlace(b1); // ... \endcode \sa Cholesky_Module, MatrixBase::llt(), MatrixBase::ldlt(), LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT. \subsection TutorialAdvSolvers_Triangular Triangular solver If the matrix \f$ A \f$ 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:
\include Tutorial_solve_triangular.cpp output: \include Tutorial_solve_triangular.out
The MatrixBase::triangularView() function constructs an object of the class TriangularView, and TriangularView::solve() then solves the system. There is also an \e in \e place variant:
\include Tutorial_solve_triangular_inplace.cpp output: \include Tutorial_solve_triangular_inplace.out
\sa MatrixBase::triangularView(), TriangularView::solve(), TriangularView::solveInPlace(), TriangularView class. \subsection TutorialAdvSolvers_Inverse Direct inversion (for small matrices) The solution of the system \f$ A \mathbf{x} = \mathbf{b} \f$ is given by \f$ \mathbf{x} = A^{-1} \mathbf{b} \f$. 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 \f$ A \f$ 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:
\include Tutorial_solve_matrix_inverse.cpp output: \include Tutorial_solve_matrix_inverse.out
Note that the function inverse() is defined in the \ref LU_Module. \sa MatrixBase::inverse(). \subsection TutorialAdvSolvers_Misc 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: \code #include // ... MatrixXf A = MatrixXf::Random(20,20); VectorXf b = VectorXf::Random(20); VectorXf x = A.svd().solve(b); SVD svdOfA(A); x = svdOfA.solve(b); \endcode %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. \code #include // ... MatrixXf A = MatrixXf::Random(20,20); VectorXf b = VectorXf::Random(20); VectorXf x = A.lu().solve(b); FullPivLU luOfA(A); x = luOfA.solve(b); \endcode See the section \ref TutorialAdvLU below. \sa class SVD, SVD::solve(), SVD_Module, class FullPivLU, LU::solve(), LU_Module. top\section TutorialAdvLU 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 \link MatrixBase::lu() lu() \endlink, which is the easiest way if you're going to use the LU decomposition only once, as in \code #include MatrixXf A = MatrixXf::Random(20,20); VectorXf b = VectorXf::Random(20); VectorXf x = A.lu().solve(b); \endcode Alternatively, you can construct a named LU decomposition, which allows you to reuse it for more than one operation: \code #include MatrixXf A = MatrixXf::Random(20,20); Eigen::FullPivLU 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; } \endcode \sa LU_Module, LU::solve(), class FullPivLU top\section TutorialAdvCholesky Cholesky todo \sa Cholesky_Module, LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT top\section TutorialAdvQR QR todo \sa QR_Module, class QR top\section TutorialAdvEigenProblems Eigen value problems todo \sa class SelfAdjointEigenSolver, class EigenSolver */ }