## Unit5.4.3Cholesky factorization algorithm (right-looking variant)

The most common algorithm for computing the Cholesky factorization of a given HPD matrix $A$ is derived as follows:

• Consider $A = L L^H \text{,}$ where $L$ is lower triangular.

Partition

\begin{equation} A = \FlaTwoByTwoSingleLine { \alpha_{11} }{ \star } { a_{21} }{ A_{22} } \quad \mbox{and} \quad L = \FlaTwoByTwoSingleLine { \lambda_{11} }{ 0 } { l_{21} }{ L_{22} }.\label{chapter05-cholesky-eqn-PartA}\tag{5.4.1} \end{equation}

Since $A$ is HPD, we know that

• By substituting these partitioned matrices into $A = L L^H$ we find that

\begin{equation*} \begin{array}{rcl} \FlaTwoByTwoSingleLine { \alpha_{11} }{ \star } { a_{21} }{ A_{22} } \amp=\amp \FlaTwoByTwoSingleLine { \lambda_{11} }{ 0 } { l_{21} }{ L_{22} } \FlaTwoByTwoSingleLine { \lambda_{11} }{ 0 } { l_{21} }{ L_{22} }^H = \FlaTwoByTwoSingleLine { \lambda_{11} }{ 0 } { l_{21} }{ L_{22} } \FlaTwoByTwoSingleLine { \bar \lambda_{11} }{ l_{21}^H } { 0 }{ L_{22} ^H} \\ \amp=\amp \FlaTwoByTwoSingleLine { \vert \lambda_{11} \vert^2 }{ \star } { \bar \lambda_{11} l_{21} } { l_{21} l_{21}^H + L_{22}^{\phantom{{}^{1}}} L_{22}^H }, \end{array} \end{equation*}

from which we conclude that

\begin{equation*} \FlaTwoByTwoSingleLineNoPar { \alpha_{11} = \vert \lambda_{11}\vert^2 }{ \star } { a_{21} = \lambda_{11} l_{21} }{ A_{22} = l_{21} l_{21}^H + L_{22} L_{22}^H } \end{equation*}

or, equivalently,

\begin{equation*} \FlaTwoByTwoSingleLineNoPar { \lambda_{11} = \pm \sqrt{ \alpha_{11} } }{ \star } { l_{21} = a_{21} / \bar \lambda_{11} } { L_{22}^{\phantom{{}^{1}}} = \Chol{ A_{22} - l_{21} l_{21}^H } } \ . \end{equation*}
• These equalities motivate the following algorithm for overwriting the lower triangular part of $A$ with the Cholesky factor of $A \text{:}$

• Partition $A \rightarrow \FlaTwoByTwoSingleLine { \alpha_{11} }{ \star } { a_{21} }{ A_{22} } \text{.}$

• Overwrite $\alpha_{11} \becomes \lambda_{11} = \sqrt{\alpha_{11}} \text{.}$ (Picking $\lambda_{11} = \sqrt{\alpha_{11}}$ makes it positive and real, and ensures uniqueness.)

• Overwrite $a_{21} \becomes l_{21} = a_{21} / \lambda_{11} \text{.}$

• Overwrite $A_{22} \becomes A_{22} - l_{21} l_{21}^H$ (updating only the lower triangular part of $A_{22}$). This operation is called a symmetric rank-1 update.

• Continue by computing the Cholesky factor of $A_{22} \text{.}$

The resulting algorithm is often called the "right-looking" variant and is summarized in Figure 5.4.3.1.

###### Homework5.4.3.1.

Give the approximate cost incurred by the algorithm in Figure 5.4.3.1 when applied to an $n \times n$ matrix.

$\frac{1}{3} n^3$ flops.

Solution

The cost of the Cholesky factorization of $A \in \C^{n \times n}$ can be analyzed as follows: In Figure 5.4.3.1 during the $k$th iteration (starting $k$ at zero) $A_{00}$ is $k \times k \text{.}$ Thus, the operations in that iteration cost

• $\alpha_{11} \becomes \sqrt{ \alpha_{11} } \text{:}$ this cost is negligible when $k$ is large.

• $a_{21} \becomes a_{21} / \alpha_{11} \text{:}$ approximately $(n-k-1)$ flops. This operation is typically implemented as $( 1 / \alpha_{11} ) a_{21} \text{.}$

• $A_{22} \becomes A_{22} - a_{21} a_{21}^H$ (updating only the lower triangular part of $A_{22}$): approximately $(n-k-1)^2$ flops.

Thus, the total cost in flops is given by

\begin{equation*} \begin{array}{l} C_{\rm Chol}( n ) \\ ~~~\approx~~~~ \lt \mbox{ sum over all iterations } \gt \\ \begin{array}[t]{c} \underbrace{ \sum_{k=0}^{n-1} (n-k-1)^2 } \\ (\mbox{Due to update of } A_{22} ) \end{array} + \begin{array}[t]{c} \underbrace{ \sum_{k=0}^{n-1} (n-k-1) } \\ (\mbox{Due to update of } a_{21} ) \end{array} \\ ~~~=~~~~ \lt \mbox{ change of variables } j = n-k-1 \gt \\ \sum_{j=0}^{n-1} j^2 + \sum_{j=0}^{n-1} j \\ ~~~\approx~~~~ \lt \sum_{j=0}^{n-1} j^2 \approx n^3/3; \sum_{j=0}^{n-1} j \approx n^2/2 \gt \\ \frac{1}{3}n^3 + \frac{1}{2} n^2 \\ ~~~ \approx ~~~~ \lt \mbox{ remove lower order term } \gt \\ \frac{1}{3}n^3. \end{array} \end{equation*}
###### Remark5.4.3.2.

Comparing the cost of the Cholesky factorization to that of the LU factorization in Homework 5.2.2.1, we see that taking advantage of symmetry cuts the cost approximately in half.

###### Homework5.4.3.2.

Implement the algorithm given in Figure 5.4.3.1 as

function [ A_out ] = Chol_right_looking( A )


by completing the code in Assignments/Week05/matlab/Chol_right_looking.m. Input is a HPD $m \times n$ matrix $A$ with only the lower triangular part stored. Output is the matrix $A$ that has its lower triangular part overwritten with the Cholesky factor. You may want to use Assignments/Week05/matlab/test_Chol_right_looking.m to check your implementation.

Solution

See Assignments/Week05/answers/Chol_right_looking.m. (Assignments/Week05/answers/Chol_right_looking.m)