# Does a symmetric positive definite matrix also have $\mathbf{A} = \mathbf{L}^T\mathbf{L}$ (where $\mathbf{L}$ is a lower triangular matrix)?

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.

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.

• 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. – sunshine Nov 6 '19 at 0:37
• @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) – Nick Alger Nov 6 '19 at 5:12
• What you call the "anti-diagonal permutation matrix" is more commonly known as the "exchange matrix". – J. M. Nov 14 '19 at 1:14
• @J.M. Thanks, I didn't know that – Nick Alger Nov 19 '19 at 18:35

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:

1. Compute $$\lambda_{22} = \sqrt{\alpha_{22}}$$, a scalar square root.
2. Scale $$\mathbf l_{21} = \lambda_{22}^{-1} \mathbf a_{21}$$, a row scaling.
3. Update $$\mathbf { \tilde A_{11} } = \mathbf A_{11} - \mathbf l_{21}^T \mathbf l_{21}$$, a rank-1 outer product update.
4. 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.

• 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. – sunshine Nov 6 '19 at 0:41
• 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. – rchilton1980 Nov 6 '19 at 1:04
• Get it. Many thanks for Prof Richilton's reply. I thouroughly understand it. Best wishes.;-) – sunshine Nov 6 '19 at 1:27

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$$.

• Thanks for your concise and clear explanation.;-) – sunshine Nov 6 '19 at 0:22