From 94393ddfff9ae66e71c5bc397e96794d56d0f062 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 18 Apr 2018 16:30:43 -0600 Subject: [PATCH] Update documentation of class Householder. --- include/deal.II/lac/householder.h | 60 +++++++++++++++++++++++-------- 1 file changed, 46 insertions(+), 14 deletions(-) diff --git a/include/deal.II/lac/householder.h b/include/deal.II/lac/householder.h index cb550657cb..3ae06ad9f8 100644 --- a/include/deal.II/lac/householder.h +++ b/include/deal.II/lac/householder.h @@ -39,10 +39,32 @@ template class Vector; /** * QR-decomposition of a full matrix. * - * Whenever an object of this class is created, it copies the matrix given and - * computes its QR-decomposition by Householder algorithm. Then, the function - * least_squares() can be used to compute the vector x minimizing - * ||Ax-b|| for a given vector b. + * This class computes the QR-decomposition of given matrix by the + * Householder algorithm. Then, the function least_squares() can be + * used to compute the vector $x$ minimizing $\|Ax-b\|$ for a given + * vector $b$. The QR decomposition of $A$ is useful for this purpose + * because the minimizer is given by the equation + * $x=(A^TA)^{-1}A^Tb=(R^TQ^TQR)^{-1}R^TQ^Tb$ which is easy to compute + * because $Q$ is an orthogonal matrix, and consequently + * $Q^TQ=I$. Thus, + * $x=(R^TR)^{-1}R^TQ^Tb=R^{-1}R^{-T}R^TQ^Tb=R^{-1}Q^Tb$. Furthermore, + * $R$ is triangular, so applying $R^{-1}$ to a vector only involves a + * backward or forward solve. + * + * + *

Implementation details

+ * + * The class does not in fact store the $Q$ and $R$ factors explicitly + * as matrices. It does store $R$, but the $Q$ factor is stored as the + * product of Householder reflections of the form $Q_i = I-v_i v_i^T$ + * where the vectors $v_i$ are so that they can be stored in the + * lower-triangular part of an underlying matrix object, whereas $R$ + * is stored in the upper triangular part. + * + * The $v_i$ vectors and the $R$ matrix now are in conflict because they + * both want to use the diagonal entry of the matrix, but we can only + * store one in these positions, of course. Consequently, the entries + * $(v_i)_i$ are stored separately in the `diagonal` member variable. * * @note Instantiations for this template are provided for @ and * @; others can be generated in application programs (see the @@ -57,7 +79,7 @@ class Householder : private FullMatrix { public: /** - * Declare type of container size. + * Declare type of container size type. */ typedef types::global_dof_index size_type; @@ -67,21 +89,23 @@ public: Householder () = default; /** - * Create an object holding the QR-decomposition of a matrix. + * Create an object holding the QR-decomposition of the matrix $A$. */ template - Householder (const FullMatrix &); + Householder (const FullMatrix &A); /** - * Compute the QR-decomposition of another matrix. + * Compute the QR-decomposition of the given matrix $A$. + * + * This overwrites any previously computed QR decomposition. */ template void - initialize (const FullMatrix &); + initialize (const FullMatrix &A); /** * Solve the least-squares problem for the right hand side src. The - * return value is the Euclidean norm of the approximation error. + * returned scalar value is the Euclidean norm of the approximation error. * * @arg @c dst contains the solution of the least squares problem on return. * @@ -94,7 +118,7 @@ public: const Vector &src) const; /** - * This function does the same as the one for BlockVectors. + * This function does the same as the previous one, but for BlockVectors. */ template double least_squares (BlockVector &dst, @@ -105,15 +129,21 @@ public: * interface. */ template - void vmult (VectorType &dst, const VectorType &src) const; + void vmult (VectorType &dst, + const VectorType &src) const; + /** + * A wrapper to least_squares() that implements multiplication with + * the transpose matrix. + */ template void Tvmult (VectorType &dst, const VectorType &src) const; private: /** - * Storage for the diagonal elements of the orthogonal transformation. + * Storage for the diagonal elements of the orthogonal + * transformation. See the class documentation for more information. */ std::vector diagonal; }; @@ -181,6 +211,7 @@ Householder::initialize(const FullMatrix &M) } + template template Householder::Householder(const FullMatrix &M) @@ -189,6 +220,7 @@ Householder::Householder(const FullMatrix &M) } + template template double @@ -208,7 +240,7 @@ Householder::least_squares (Vector &dst, // m > n, m = src.n, n = dst.n // Multiply Q_n ... Q_2 Q_1 src - // Where Q_i = I-v_i v_i^T + // Where Q_i = I - v_i v_i^T for (size_type j=0; j