Skip to main content

Subsection 5.4.3 Cholesky 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.

\begin{equation*} \newcommand{\routinename}{A = \mbox{Chol-right-looking}( A )} \newcommand{\guard}{ n( A_{TL} ) \lt n( A ) } \newcommand{\partitionings}{ A \rightarrow \FlaTwoByTwo{A_{TL}}{A_{TR}} {A_{BL}}{A_{BR}} } \newcommand{\partitionsizes}{ A_{TL} {\rm ~is~} 0 \times 0 } \newcommand{\repartitionings}{ \FlaTwoByTwo{A_{TL}}{A_{TR}} {A_{BL}}{A_{BR}} \rightarrow \FlaThreeByThreeBR{A_{00}}{a_{01}}{A_{02}} {a_{10}^T}{\alpha_{11}}{a_{12}^T} {A_{20}}{a_{21}}{A_{22}} } \newcommand{\moveboundaries}{ \FlaTwoByTwo{A_{TL}}{A_{TR}} {A_{BL}}{A_{BR}} \leftarrow \FlaThreeByThreeTL{A_{00}}{a_{01}}{A_{02}} {a_{10}^T}{\alpha_{11}}{a_{12}^T} {A_{20}}{a_{21}}{A_{22}} } \newcommand{\update}{ \begin{array}{l} \alpha_{11} := \lambda_{11} = \sqrt{ \alpha_{11} } \\ a_{21} := l_{21} = a_{21} / \alpha_{11} \\ A_{22} := A_{22} - a_{21} a_{21}^H ~~~~~~~~ \mbox{(syr: update only lower triangular part)} \end{array} } \FlaAlgorithm \end{equation*}
Figure 5.4.3.1. Cholesky factorization algorithm (right-looking variant). The operation "syr" refers to "symmetric rank-1 update", which performs a rank-1 update, updating only the lower triangular part of the matrix in this algorithm.
Homework 5.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.

Answer

\(\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*}
Remark 5.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.

Homework 5.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)

login