diff options
| -rw-r--r-- | doc/C03_TutorialArrayClass.dox | 3 | ||||
| -rw-r--r-- | doc/C04_TutorialBlockOperations.dox | 5 | ||||
| -rw-r--r-- | doc/C05_TutorialAdvancedInitialization.dox | 7 | ||||
| -rw-r--r-- | doc/C06_TutorialLinearAlgebra.dox | 141 | ||||
| -rw-r--r-- | doc/examples/TutorialLinAlgComputeTwice.cpp | 23 | ||||
| -rw-r--r-- | doc/examples/TutorialLinAlgExComputeSolveError.cpp | 14 | ||||
| -rw-r--r-- | doc/examples/TutorialLinAlgExSolveColPivHouseholderQR.cpp | 6 | ||||
| -rw-r--r-- | doc/examples/TutorialLinAlgExSolveLDLT.cpp | 16 | ||||
| -rw-r--r-- | doc/examples/TutorialLinAlgInverseDeterminant.cpp | 16 | ||||
| -rw-r--r-- | doc/examples/TutorialLinAlgRankRevealing.cpp | 20 | ||||
| -rw-r--r-- | doc/examples/TutorialLinAlgSelfAdjointEigenSolver.cpp | 17 | ||||
| -rw-r--r-- | doc/examples/TutorialLinAlgSetThreshold.cpp | 19 |
12 files changed, 273 insertions, 14 deletions
diff --git a/doc/C03_TutorialArrayClass.dox b/doc/C03_TutorialArrayClass.dox index 2444e6ed8..4e3b4b08e 100644 --- a/doc/C03_TutorialArrayClass.dox +++ b/doc/C03_TutorialArrayClass.dox @@ -1,6 +1,6 @@ namespace Eigen { -/** \page TutorialArrayClass Tutorial page 3 - The Array Class +/** \page TutorialArrayClass Tutorial page 3 - The %Array class \ingroup Tutorial \li \b Previous: \ref TutorialMatrixArithmetic @@ -238,6 +238,7 @@ array3 = array1.abs2(); </table> </td></tr></table> +\li \b Next: \ref TutorialBlockOperations **/ } diff --git a/doc/C04_TutorialBlockOperations.dox b/doc/C04_TutorialBlockOperations.dox index 689828481..3f2916945 100644 --- a/doc/C04_TutorialBlockOperations.dox +++ b/doc/C04_TutorialBlockOperations.dox @@ -1,10 +1,10 @@ namespace Eigen { -/** \page TutorialBlockOperations Tutorial page 4 - Block operations +/** \page TutorialBlockOperations Tutorial page 4 - %Block operations \ingroup Tutorial \li \b Previous: \ref TutorialArrayClass -\li \b Next: (not yet written) +\li \b Next: \ref TutorialAdvancedInitialization This tutorial explains the essentials of Block operations together with many examples. @@ -288,6 +288,7 @@ Output: \include Tutorial_BlockOperations_vector.out </td></tr></table> +\li \b Next: \ref TutorialAdvancedInitialization */ diff --git a/doc/C05_TutorialAdvancedInitialization.dox b/doc/C05_TutorialAdvancedInitialization.dox index a658986c9..1e048d4b2 100644 --- a/doc/C05_TutorialAdvancedInitialization.dox +++ b/doc/C05_TutorialAdvancedInitialization.dox @@ -1,8 +1,11 @@ namespace Eigen { -/** \page TutorialAdvancedInitialization Tutorial - Advanced initialization +/** \page TutorialAdvancedInitialization Tutorial page 5 - Advanced initialization \ingroup Tutorial +\li \b Previous: \ref TutorialBlockOperations +\li \b Next: \ref TutorialLinearAlgebra + \section TutorialMatrixArithmCommaInitializer Comma initializer Eigen offers a comma initializer syntax which allows to set all the coefficients @@ -24,6 +27,8 @@ TODO mention using the comma initializer to fill a block xpr like m.row(i) << 1, TODO add more sections about Identity(), Zero(), Constant(), Random(), LinSpaced(). +\li \b Next: \ref TutorialLinearAlgebra + */ } diff --git a/doc/C06_TutorialLinearAlgebra.dox b/doc/C06_TutorialLinearAlgebra.dox index f31fcfc0f..7c851ec34 100644 --- a/doc/C06_TutorialLinearAlgebra.dox +++ b/doc/C06_TutorialLinearAlgebra.dox @@ -3,14 +3,14 @@ namespace Eigen { /** \page TutorialLinearAlgebra Tutorial page 6 - Linear algebra and decompositions \ingroup Tutorial -\li \b Previous: TODO +\li \b Previous: \ref TutorialAdvancedInitialization \li \b Next: TODO This tutorial explains how to solve linear systems, compute various decompositions such as LU, -QR, SVD, eigendecompositions... for more advanced topics, don't miss our special page on +QR, %SVD, eigendecompositions... for more advanced topics, don't miss our special page on \ref TopicLinearAlgebraDecompositions "this topic". -\section TutorialLinAlgBasicSolve How do I solve a system of linear equations? +\section TutorialLinAlgBasicSolve Basic linear solving \b The \b problem: You have a system of equations, that you have written as a single matrix equation \f[ Ax \: = \: b \f] @@ -26,10 +26,10 @@ and is a good compromise: </tr> </table> -In this example, the colPivHouseholderQr() method returns an object of class ColPivHouseholderQR. This line could -have been replaced by: +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); +ColPivHouseholderQR<Matrix3f> dec(A); Vector3f x = dec.solve(b); \endcode @@ -107,11 +107,138 @@ depending on your matrix and the trade-off you want to make: All of these decompositions offer a solve() method that works as in the above example. -For a much more complete table comparing all decompositions supported by Eigen (notice that Eigen +For example, if your matrix is positive definite, the above table says that a very good +choice is then the LDLT decomposition. Here's an example, also demonstrating that using a general +matrix (not a vector) as right hand side is possible. + +<table class="tutorial_code"> +<tr> + <td>\include TutorialLinAlgExSolveLDLT.cpp </td> + <td>output: \verbinclude TutorialLinAlgExSolveLDLT.out </td> +</tr> +</table> + +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 TutorialLinAlgSolutionExists Checking if a solution really exists + +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: + +<table class="tutorial_code"> +<tr> + <td>\include TutorialLinAlgExComputeSolveError.cpp </td> + <td>output: \verbinclude TutorialLinAlgExComputeSolveError.out </td> +</tr> +</table> + +\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. + +<table class="tutorial_code"> +<tr> + <td>\include TutorialLinAlgSelfAdjointEigenSolver.cpp </td> + <td>output: \verbinclude TutorialLinAlgSelfAdjointEigenSolver.out </td> +</tr> +</table> + +\section TutorialLinAlgEigensolving Computing inverse and determinant + +First of all, make sure that you really want this. While inverse and determinant are fundamental mathematical concepts, +in \em numerical linear algebra they are not as popular as in pure mathematics. Inverse computations are often +advantageously replaced by solve() operations, and the determinant is often \em not a good way of checking if a matrix +is invertible. + +However, for \em very \em small matrices, the above is not 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: +<table class="tutorial_code"> +<tr> + <td>\include TutorialLinAlgInverseDeterminant.cpp </td> + <td>output: \verbinclude TutorialLinAlgInverseDeterminant.out </td> +</tr> +</table> + +\section TutorialLinAlgLeastsquares Least squares solving +Eigen doesn't currently provide built-in linear least squares solving functions, but you can easily compute that yourself +from Eigen's decompositions. The most reliable way is to use a SVD (or better yet, JacobiSVD), and in the future +these classes will offer methods for least squares solving. Another, potentially faster way, is to use a LLT decomposition +of the normal matrix. In any case, just read any reference text on least squares, and it will be very easy for you +to implement any linear least squares computation on top of Eigen. + +\section TutorialLinAlgSeparateComputation 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: +\li all decompositions have a default constructor, +\li 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: + +<table class="tutorial_code"> +<tr> + <td>\include TutorialLinAlgComputeTwice.cpp </td> + <td>output: \verbinclude TutorialLinAlgComputeTwice.out </td> +</tr> +</table> + +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: +\code +HouseholderQR<MatrixXf> qr(50,50); +MatrixXf A = MatrixXf::Random(50,50); +qr.compute(A); // no dynamic memory allocation +\endcode + +\section TutorialLinAlgRankRevealing 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 \ref TopicLinearAlgebraDecompositions "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 some are also providing methods to compute the kernel (null-space) and image (column-space) of the matrix, as is the +case with FullPivLU: + +<table class="tutorial_code"> +<tr> + <td>\include TutorialLinAlgRankRevealing.cpp </td> + <td>output: \verbinclude TutorialLinAlgRankRevealing.out </td> +</tr> +</table> + +Of course, any rank computation depends on the choice of an arbitrary threshold, since practically no +floating-point matrix is \em 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 compute(), as in this example: + +<table class="tutorial_code"> +<tr> + <td>\include TutorialLinAlgSetThreshold.cpp </td> + <td>output: \verbinclude TutorialLinAlgSetThreshold.out </td> +</tr> +</table> + +\li \b Next: TODO */ diff --git a/doc/examples/TutorialLinAlgComputeTwice.cpp b/doc/examples/TutorialLinAlgComputeTwice.cpp new file mode 100644 index 000000000..06ba6461a --- /dev/null +++ b/doc/examples/TutorialLinAlgComputeTwice.cpp @@ -0,0 +1,23 @@ +#include <iostream> +#include <Eigen/Dense> + +using namespace std; +using namespace Eigen; + +int main() +{ + Matrix2f A, b; + LLT<Matrix2f> llt; + A << 2, -1, -1, 3; + b << 1, 2, 3, 1; + cout << "Here is the matrix A:\n" << A << endl; + cout << "Here is the right hand side b:\n" << b << endl; + cout << "Computing LLT decomposition..." << endl; + llt.compute(A); + cout << "The solution is:\n" << llt.solve(b) << endl; + A(1,1)++; + cout << "The matrix A is now:\n" << A << endl; + cout << "Computing LLT decomposition..." << endl; + llt.compute(A); + cout << "The solution is now:\n" << llt.solve(b) << endl; +} diff --git a/doc/examples/TutorialLinAlgExComputeSolveError.cpp b/doc/examples/TutorialLinAlgExComputeSolveError.cpp new file mode 100644 index 000000000..f362fb71a --- /dev/null +++ b/doc/examples/TutorialLinAlgExComputeSolveError.cpp @@ -0,0 +1,14 @@ +#include <iostream> +#include <Eigen/Dense> + +using namespace std; +using namespace Eigen; + +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 + cout << "The relative error is:\n" << relative_error << endl; +} diff --git a/doc/examples/TutorialLinAlgExSolveColPivHouseholderQR.cpp b/doc/examples/TutorialLinAlgExSolveColPivHouseholderQR.cpp index 29c22be41..3a99a94d7 100644 --- a/doc/examples/TutorialLinAlgExSolveColPivHouseholderQR.cpp +++ b/doc/examples/TutorialLinAlgExSolveColPivHouseholderQR.cpp @@ -10,8 +10,8 @@ int main() Vector3f b; A << 1,2,3, 4,5,6, 7,8,10; b << 3, 3, 4; - cout << "Here is the matrix A:" << endl << A << endl; - cout << "Here is the vector b:" << endl << b << endl; + cout << "Here is the matrix A:\n" << A << endl; + cout << "Here is the vector b:\n" << b << endl; Vector3f x = A.colPivHouseholderQr().solve(b); - cout << "The solution is:" << endl << x << endl; + cout << "The solution is:\n" << x << endl; } diff --git a/doc/examples/TutorialLinAlgExSolveLDLT.cpp b/doc/examples/TutorialLinAlgExSolveLDLT.cpp new file mode 100644 index 000000000..f8beacd27 --- /dev/null +++ b/doc/examples/TutorialLinAlgExSolveLDLT.cpp @@ -0,0 +1,16 @@ +#include <iostream> +#include <Eigen/Dense> + +using namespace std; +using namespace Eigen; + +int main() +{ + Matrix2f A, b; + A << 2, -1, -1, 3; + b << 1, 2, 3, 1; + cout << "Here is the matrix A:\n" << A << endl; + cout << "Here is the right hand side b:\n" << b << endl; + Matrix2f x = A.ldlt().solve(b); + cout << "The solution is:\n" << x << endl; +} diff --git a/doc/examples/TutorialLinAlgInverseDeterminant.cpp b/doc/examples/TutorialLinAlgInverseDeterminant.cpp new file mode 100644 index 000000000..43970ff05 --- /dev/null +++ b/doc/examples/TutorialLinAlgInverseDeterminant.cpp @@ -0,0 +1,16 @@ +#include <iostream> +#include <Eigen/Dense> + +using namespace std; +using namespace Eigen; + +int main() +{ + Matrix3f A; + A << 1, 2, 1, + 2, 1, 0, + -1, 1, 2; + cout << "Here is the matrix A:\n" << A << endl; + cout << "The determinant of A is " << A.determinant() << endl; + cout << "The inverse of A is:\n" << A.inverse() << endl; +}
\ No newline at end of file diff --git a/doc/examples/TutorialLinAlgRankRevealing.cpp b/doc/examples/TutorialLinAlgRankRevealing.cpp new file mode 100644 index 000000000..c5165077f --- /dev/null +++ b/doc/examples/TutorialLinAlgRankRevealing.cpp @@ -0,0 +1,20 @@ +#include <iostream> +#include <Eigen/Dense> + +using namespace std; +using namespace Eigen; + +int main() +{ + Matrix3f A; + A << 1, 2, 5, + 2, 1, 4, + 3, 0, 3; + cout << "Here is the matrix A:\n" << A << endl; + FullPivLU<Matrix3f> lu_decomp(A); + cout << "The rank of A is " << lu_decomp.rank() << endl; + cout << "Here is a matrix whose columns form a basis of the null-space of A:\n" + << lu_decomp.kernel() << endl; + cout << "Here is a matrix whose columns form a basis of the column-space of A:\n" + << lu_decomp.image(A) << endl; // yes, have to pass the original A +} diff --git a/doc/examples/TutorialLinAlgSelfAdjointEigenSolver.cpp b/doc/examples/TutorialLinAlgSelfAdjointEigenSolver.cpp new file mode 100644 index 000000000..e98444347 --- /dev/null +++ b/doc/examples/TutorialLinAlgSelfAdjointEigenSolver.cpp @@ -0,0 +1,17 @@ +#include <iostream> +#include <Eigen/Dense> + +using namespace std; +using namespace Eigen; + +int main() +{ + Matrix2f A; + A << 1, 2, 2, 3; + cout << "Here is the matrix A:\n" << A << endl; + SelfAdjointEigenSolver<Matrix2f> eigensolver(A); + cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl; + cout << "Here's a matrix whose columns are eigenvectors of A " + << "corresponding to these eigenvalues:\n" + << eigensolver.eigenvectors() << endl; +} diff --git a/doc/examples/TutorialLinAlgSetThreshold.cpp b/doc/examples/TutorialLinAlgSetThreshold.cpp new file mode 100644 index 000000000..e0927cf27 --- /dev/null +++ b/doc/examples/TutorialLinAlgSetThreshold.cpp @@ -0,0 +1,19 @@ +#include <iostream> +#include <Eigen/Dense> + +using namespace std; +using namespace Eigen; + +int main() +{ + Matrix2d A; + FullPivLU<Matrix2d> lu; + A << 2, 1, + 2, 0.9999999999; + lu.compute(A); + cout << "By default, the rank of A is found to be " << lu.rank() << endl; + cout << "Now recomputing the LU decomposition with threshold 1e-5" << endl; + lu.setThreshold(1e-5); + lu.compute(A); + cout << "The rank of A is found to be " << lu.rank() << endl; +} |