diff options
Diffstat (limited to 'doc/C05_TutorialLinearAlgebra.dox')
| -rw-r--r-- | doc/C05_TutorialLinearAlgebra.dox | 310 |
1 files changed, 0 insertions, 310 deletions
diff --git a/doc/C05_TutorialLinearAlgebra.dox b/doc/C05_TutorialLinearAlgebra.dox deleted file mode 100644 index af2096bf3..000000000 --- a/doc/C05_TutorialLinearAlgebra.dox +++ /dev/null @@ -1,310 +0,0 @@ -namespace Eigen { - -/** \page TutorialAdvancedLinearAlgebra Tutorial 3/4 - Advanced linear algebra - \ingroup Tutorial - -<div class="eimainmenu">\ref index "Overview" - | \ref TutorialCore "Core features" - | \ref TutorialGeometry "Geometry" - | \b Advanced \b linear \b algebra - | \ref TutorialSparse "Sparse matrix" -</div> - -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: - -<table class="tutorial_code"><tr><td> -\include Tutorial_PartialLU_solve.cpp -</td><td> -output: \include Tutorial_PartialLU_solve.out -</td></tr></table> - -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: - -<table class="tutorial_code"><tr><td> -\include Tutorial_solve_multiple_rhs.cpp -</td><td> -output: \include Tutorial_solve_multiple_rhs.out -</td></tr></table> - -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. - -<table class="tutorial_code"><tr><td> -\include Tutorial_solve_reuse_decomposition.cpp -</td><td> -output: \include Tutorial_solve_reuse_decomposition.out -</td></tr></table> - -\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: - -<table class="tutorial_code"><tr><td> -\include Tutorial_solve_singular.cpp -</td><td> -output: \include Tutorial_solve_singular.out -</td></tr></table> - -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 <Eigen/Cholesky> -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<MatrixXf> 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: - -<table class="tutorial_code"><tr><td> -\include Tutorial_solve_triangular.cpp -</td><td> -output: \include Tutorial_solve_triangular.out -</td></tr></table> - -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: - -<table class="tutorial_code"><tr><td> -\include Tutorial_solve_triangular_inplace.cpp -</td><td> -output: \include Tutorial_solve_triangular_inplace.out -</td></tr></table> - -\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: - -<table class="tutorial_code"><tr><td> -\include Tutorial_solve_matrix_inverse.cpp -</td><td> -output: \include Tutorial_solve_matrix_inverse.out -</td></tr></table> - -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 <Eigen/SVD> -// ... -MatrixXf A = MatrixXf::Random(20,20); -VectorXf b = VectorXf::Random(20); -VectorXf x = A.svd().solve(b); -SVD<MatrixXf> 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 <Eigen/LU> -// ... -MatrixXf A = MatrixXf::Random(20,20); -VectorXf b = VectorXf::Random(20); -VectorXf x = A.lu().solve(b); -FullPivLU<MatrixXf> 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. - - - -<a href="#" class="top">top</a>\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 <Eigen/LU> -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 <Eigen/LU> -MatrixXf A = MatrixXf::Random(20,20); -Eigen::FullPivLU<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; -} -\endcode - -\sa LU_Module, LU::solve(), class FullPivLU - -<a href="#" class="top">top</a>\section TutorialAdvCholesky Cholesky -todo - -\sa Cholesky_Module, LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT - -<a href="#" class="top">top</a>\section TutorialAdvQR QR -todo - -\sa QR_Module, class QR - -<a href="#" class="top">top</a>\section TutorialAdvEigenProblems Eigen value problems -todo - -\sa class SelfAdjointEigenSolver, class EigenSolver - -*/ - -} |