7
$\begingroup$

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.

$\endgroup$
8
$\begingroup$

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.

$\endgroup$
  • $\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$ @J.M. Thanks, I didn't know that $\endgroup$ – Nick Alger Nov 19 '19 at 18:35
6
$\begingroup$

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.

$\endgroup$
  • $\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
4
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.