Cholesky, Doolittle and Crout Factorizations
Source inspiration: (Mathew 2000-2019).
LU-Factorization Background
A nonsingular matrix \(A \in \mathbb{R}^{n\times n}\) has an LU factorization if \[ A = LU, \] where \(L\) is lower triangular and \(U\) is upper triangular.
Common normalizations are:
- Doolittle form: diagonal entries of \(L\) are all \(1\).
- Crout form: diagonal entries of \(U\) are all \(1\).
- Cholesky form (for symmetric positive definite \(A\)): \(A=L_cL_c^\top=R^\top R\).
If \(A=LU\), then solving \(A\mathbf{x}=\mathbf{b}\) is equivalent to \[ L\mathbf{y}=\mathbf{b},\qquad U\mathbf{x}=\mathbf{y}. \] So one forward substitution and one back substitution solve the original system.
Formulas
For \(A=D-L-U\) splitting (diagonal part \(D\), strict lower part \(L\), strict upper part \(U\)), the three practical forms are:
\[ A=L_DU_D \quad \text{(Doolittle)}, ag{1} \]
\[ A=L_CU_C \quad \text{(Crout)}, ag{2} \]
\[ A=L_cL_c^\top=R^\top R \quad \text{(Cholesky, SPD matrices)}. ag{3} \]
Example 1 - One Matrix, Three Factorizations
Factorize
\[ A= \begin{pmatrix} 2&1&1&3&2\\ 1&2&2&1&1\\ 1&2&9&1&5\\ 3&1&1&7&1\\ 2&1&5&1&8 \end{pmatrix}. \]
Doolittle (\(\operatorname{diag}(L_D)=1\))
\[ L_D= \begin{pmatrix} 1&0&0&0&0\\ frac12&1&0&0&0\\ frac12&1&1&0&0\\ frac32&-\tfrac13&0&1&0\\ 1&0&\tfrac47&-\tfrac67&1 \end{pmatrix}, \qquad U_D= \begin{pmatrix} 2&1&1&3&2\\ 0&\tfrac32&\tfrac32&-\tfrac12&0\\ 0&0&7&0&4\\ 0&0&0&\tfrac73&-2\\ 0&0&0&0&2 \end{pmatrix}. \]
Crout (\(\operatorname{diag}(U_C)=1\))
\[ L_C= \begin{pmatrix} 2&0&0&0&0\\ 1&\tfrac32&0&0&0\\ 1&\tfrac32&7&0&0\\ 3&-\tfrac12&0&\tfrac73&0\\ 2&0&4&-2&2 \end{pmatrix}, \qquad U_C= \begin{pmatrix} 1&\tfrac12&\tfrac12&\tfrac32&1\\ 0&1&1&-\tfrac13&0\\ 0&0&1&0&\tfrac47\\ 0&0&0&1&-\tfrac67\\ 0&0&0&0&1 \end{pmatrix}. \]
Cholesky (\(A=R^\top R\))
One valid upper-triangular Cholesky factor is
\[ R= \begin{pmatrix} \sqrt2&\tfrac1{\sqrt2}&\tfrac1{\sqrt2}&\tfrac3{\sqrt2}&\sqrt2\\ 0&\sqrt{\tfrac32}&\sqrt{\tfrac32}&-\tfrac1{\sqrt6}&0\\ 0&0&\sqrt7&0&\tfrac4{\sqrt7}\\ 0&0&0&\sqrt{\tfrac73}&-2\sqrt{\tfrac37}\\ 0&0&0&0&\sqrt2 \end{pmatrix}. \]
Verification: \[ L_DU_D=A,\qquad L_CU_C=A,\qquad R^\top R=A. \]
Example 2 - Solve \(A\mathbf{x}=\mathbf{b}\) with Doolittle Factors
Use the same \(A\) and
\[ \mathbf{b}= \begin{pmatrix} -2\\4\\3\\-5\\1 \end{pmatrix}. \]
Step 1: Solve \(L_D\mathbf{y}=\mathbf{b}\).
\[ \mathbf{y}= \begin{pmatrix} -2\\5\\-1\\-\tfrac13\\\tfrac{23}{7} \end{pmatrix}. \]
Step 2: Solve \(U_D\mathbf{x}=\mathbf{y}\).
\[ \mathbf{x}= \begin{pmatrix} -\tfrac{629}{98}\\[2pt] frac{237}{49}\\[2pt] -\tfrac{53}{49}\\[2pt] frac{62}{49}\\[2pt] frac{23}{14} \end{pmatrix} \approx \begin{pmatrix} -6.41837\\4.83673\\-1.08163\\1.26531\\1.64286 \end{pmatrix}. \]
Verification:
| Row | Left-hand side \(A\mathbf{x}\) | Right-hand side |
|---|---|---|
| 1 | \(2x_1+x_2+x_3+3x_4+2x_5=-2\) | \(-2\;\checkmark\) |
| 2 | \(x_1+2x_2+2x_3+x_4+x_5=4\) | \(4\;\checkmark\) |
| 3 | \(x_1+2x_2+9x_3+x_4+5x_5=3\) | \(3\;\checkmark\) |
| 4 | \(3x_1+x_2+x_3+7x_4+x_5=-5\) | \(-5\;\checkmark\) |
| 5 | \(2x_1+x_2+5x_3+x_4+8x_5=1\) | \(1\;\checkmark\) |