nearly complete page 6 / linear algebra + examples

fix the previous/next links

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
 
 **/
 }

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
 
 */
 

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
+
 */
 
 }

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
 
 */
 

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;
+}

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;
+}

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;
 }

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;
+}

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;
+}

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
+}

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;
+}

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;
+}