As we know, for a symmetric positive definite (SPD) matrix $\mathbf{A}$, there is a theorem about the Cholesky factorization $\mathbf{A}= \mathbf{L}\mathbf{L}^T$, where $\mathbf{L}$ is a lower triangular matrix. I am a little curious whether the factorization $\mathbf{A} = \mathbf{L}^T\mathbf{L}$ exists for a SPD matrix $\mathbf{A}$, where $\mathbf{L}$ is still a lower triangular matrix. Textbooks just give the first theorem not the second form. Though a teacher said that the second form also exists, I still have some doubts about that. Any hints and suggestions are welcome.
$\begingroup$
$\endgroup$
$\begingroup$
$\endgroup$
Let $P$ be the anti-diagonal permutation matrix, $$P = \begin{bmatrix} & & & 1 \\ & & 1 \\ & 1 \\ 1 \end{bmatrix}$$ so that $PAP$ is the version of $A$ with rows and columns reversed. The first $P$ swaps the rows, and the last $P$ swaps the columns. We have the Cholesky decomposition $$PAP=LL^T$$ which implies $$A = (PLP)(PLP)^T,$$ since $P^{-1}=P$. But now $PLP$ is upper triangular, so this is a factorization of $A$ into a product of an upper triangular matrix times a lower triangular matrix.
-
$\begingroup$ Thanks Prof Nick, I get it. I also want to ask that though there are many factorizations for SPD matrix, but why matlab just choose this form $A=LL^T$, where L is lower triangular. Is this form the most classic form？so matlab choose this one and many monographs also introduce this one? And the monographs do not introdce a little about other factorizations. $\endgroup$ – sunshine Nov 6 '19 at 0:37
-
$\begingroup$ @Zhen-WeiSun I'm not sure why books present it this way. But since the Cholesky factorization of a SPD matrix is unique, there are really only two possible symmetric triangular factorizations: $LL^T$ and $L^TL$. These factorizations are more useful than the LU factorization because they can be used to draw samples from a Gaussian with covariance $A^{-1}$ by computing $L^{-1}z$, $z\sim N(0,I)$. They are cheaper than the QR and eigenvalue decompositions because they often have less fill-in when $A$ is sparse. (also, I'm not a professor, but I'm flattered that you thought I was) $\endgroup$ – Nick Alger Nov 6 '19 at 5:12
-
1$\begingroup$ What you call the "anti-diagonal permutation matrix" is more commonly known as the "exchange matrix". $\endgroup$ – J. M. Nov 14 '19 at 1:14
-
$\begingroup$
$\endgroup$
Yes, for an SPD matrix $\mathbf A$ there are a variety of Cholesky-like decompositions, you can derive the $\mathbf A = \mathbf L^T \mathbf L$ variant by first writing down an educated/structured guess..
$\begin{bmatrix} \mathbf A_{11} & \mathbf a_{21}^T \\ \mathbf a_{21} & \alpha_{22} \end{bmatrix} = \begin{bmatrix} \mathbf L_{11}^T & \mathbf l_{21}^T \\ \mathbf 0 & \lambda_{22} \end{bmatrix} \begin{bmatrix} \mathbf L_{11} & \mathbf 0 \\ \mathbf l_{21} & \lambda_{22} \end{bmatrix} $
.. where the bold items are matrices/row vectors and the greek items are scalars. Next, multiply out the right side and equate it block-by-block to the left side, to deduce the following relationships between $\mathbf A$ and $\mathbf L$:
- Equate (1,1) blocks: $\mathbf A_{11} = \mathbf L_{11}^T \mathbf L_{11} + \mathbf l_{21}^T \mathbf l_{21}$
- Equate (2,1) blocks: $\mathbf a_{21} = \lambda_{22} \mathbf l_{21}$
- Equate (2,2) blocks: $\alpha_{22} = \lambda_{22} \lambda_{22}$
When properly sequenced (work from known quantities towards unknown ones), these relationships define the algorithm:
- Compute $\lambda_{22} = \sqrt{\alpha_{22}}$, a scalar square root.
- Scale $\mathbf l_{21} = \lambda_{22}^{-1} \mathbf a_{21}$, a row scaling.
- Update $ \mathbf { \tilde A_{11} } = \mathbf A_{11} - \mathbf l_{21}^T \mathbf l_{21}$, a rank-1 outer product update.
- Factor $\mathbf { \tilde A_{11} } = \mathbf L_{11}^T \mathbf L_{11}$, tail-recursion into the upper left submatrix.
LAPACK uses a similar algorithm whenever you apply Choleskly decomposition [potrf] to "upper" triangular storage (it forms $\mathbf A = \mathbf U \mathbf U^T$). All four of the decompositions ($\mathbf L \mathbf L^T$, $\mathbf L^T \mathbf L$, $\mathbf U \mathbf U^T$, $\mathbf U^T \mathbf U$) are possible and can be derived using similar ideas.
-
$\begingroup$ Thanks Prof Rchilton for your clear and detailed explanation, I get it. I also want to ask that though there are many factorizations for SPD matrix, but why matlab just choose this form $A=LL^T$, where L is lower triangular. Is this form the most classic form？so matlab choose this one and many monographs also introduce this one? And the monographs do not introdce a little about other factorizations. $\endgroup$ – sunshine Nov 6 '19 at 0:41
-
$\begingroup$ As you can see from the answers, all these variants are basically equivalent : triangular factorizations, that furnish the inverse operator, and properly exploit the symmetry/positivity. I guess the $\mathbf L \mathbf L^T$ variant is considered the canonical choice .. probably because it's the one that is most similar in structure/implementation to the $\mathbf L \mathbf U$ decomposition. That being said, there's nothing special about $\mathbf L \mathbf U$ either! It too can be rearranged into other/similar triangular factorizations, just like Cholesky. $\endgroup$ – rchilton1980 Nov 6 '19 at 1:04
-
$\begingroup$ Get it. Many thanks for Prof Richilton's reply. I thouroughly understand it. Best wishes.;-) $\endgroup$ – sunshine Nov 6 '19 at 1:27
$\begingroup$
$\endgroup$
Since a SPD matrix is invertible, we can make the Cholesky decomposition $A^{-1} = PP^T$.
Since $A$ is non-singular, so is $P$, and the inverse of a triangular matrix is triangular, so writing $L = P^{-1}$ we have $A^{-1}$ = $L^{-1}L^{-T}$.
Inverting both sides gives $A = L^TL$.
-
$\begingroup$ Thanks for your concise and clear explanation.;-) $\endgroup$ – sunshine Nov 6 '19 at 0:22
