From 76873e9969607c44c937bdbfb69acb82b5113fe0 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 24 Jun 2020 11:20:25 -0600 Subject: [PATCH] Augment documentation of project_onto_orthogonal_tensors(). --- include/deal.II/base/tensor.h | 68 ++++++++++++++++++++++++++++++----- 1 file changed, 59 insertions(+), 9 deletions(-) diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index 956bc7e674..713e08a5d0 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -2640,16 +2640,66 @@ cofactor(const Tensor<2, dim, Number> &t) /** - * Return the nearest orthogonal matrix by - * combining the products of the SVD decomposition: $\mathbf U \mathbf{V}^T$, - * where $\mathbf U$ and $\mathbf V$ are computed from the SVD decomposition: - * $\mathbf U \mathbf S \mathbf V^T$, - * effectively replacing $\mathbf S$ with the identity matrix. - * @param A The tensor which to find the closest orthogonal tensor to. - * @pre @p Number must be either `float` or `double`. + * Return the nearest orthogonal matrix + * $\hat {\mathbf A}=\mathbf U \mathbf{V}^T$ by + * combining the products of the singular value decomposition (SVD) + * ${\mathbf A}=\mathbf U \mathbf S \mathbf V^T$ for a given input + * ${\mathbf A}$, effectively replacing $\mathbf S$ with the identity matrix. + * + * This is a (nonlinear) + * [projection + * operation](https://en.wikipedia.org/wiki/Projection_(mathematics)) since when + * applied twice, we have $\hat{\hat{\mathbf A}}=\hat{\mathbf A}$ as is easy to + * see. (That is because the SVD of $\hat {\mathbf A}$ is simply + * $\mathbf U \mathbf I \mathbf{V}^T$.) Furthermore, $\hat {\mathbf A}$ is + * really an orthogonal matrix because orthogonal matrices have to satisfy + * ${\hat {\mathbf A}}^T \hat {\mathbf A}={\mathbf I}$, which here implies + * that + * @f{align*}{ + * {\hat {\mathbf A}}^T \hat {\mathbf A} + * &= + * \left(\mathbf U \mathbf{V}^T\right)^T\left(\mathbf U \mathbf{V}^T\right) + * \\ + * &= + * \mathbf V \mathbf{U}^T + * \mathbf U \mathbf{V}^T + * \\ + * &= + * \mathbf V \left(\mathbf{U}^T + * \mathbf U\right) \mathbf{V}^T + * \\ + * &= + * \mathbf V \mathbf I \mathbf{V}^T + * \\ + * &= + * \mathbf V \mathbf{V}^T + * \\ + * &= + * \mathbf I + * @f} + * due to the fact that the $\mathbf U$ and $\mathbf V$ factors that come out + * of the SVD are themselves orthogonal matrices. + * + * @param A The tensor for which to find the closest orthogonal tensor. + * @tparam Number The type used to store the entries of the tensor. + * Must be either `float` or `double`. * @pre In order to use this function, this program must be linked with the - * LAPACK library. - * @pre @p A must not be singular. + * LAPACK library. + * @pre @p A must not be singular. This is because, conceptually, the problem + * to be solved here is trying to find a matrix $\hat{\mathbf A}$ that + * minimizes some kind of distance from $\mathbf A$ while satisfying the + * quadratic constraint + * ${\hat {\mathbf A}}^T \hat {\mathbf A}={\mathbf I}$. This is not so + * dissimilar to the kind of problem where one wants to find a vector + * $\hat{\mathbf x}\in{\mathbb R}^n$ that minimizes the quadratic objective + * function $\|\hat {\mathbf x} - \mathbf x\|^2$ for a given $\mathbf x$ + * subject to the constraint $\|\mathbf x\|^2=1$ -- in other + * words, we are seeking the point $\hat{\mathbf x}$ on the unit sphere + * that is closest to $\mathbf x$. This problem has a solution for all + * $\mathbf x$ except if $\mathbf x=0$. The corresponding condition + * for the problem we are considering here is that $\mathbf A$ must not + * have a zero eigenvalue. + * * @relatesalso Tensor */ template -- 2.39.5