doc/C05_TutorialLinearAlgebra.dox - mirror - Git at Google

 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
 (&le;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

 */

 }
	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

	*/

	}