Advanced Linear Algebra: Foundations to Frontiers

Subsection 5.5.1 Other LU factorization algorithms

There are actually five different (unblocked) algorithms for computing the LU factorization that were discovered over the course of the centuries. Here we show how to systematically derive all five. For details, we suggest Week 6 of our Massive Open Online Course titled "LAFF-On Programming for Correctness" [29].

Finding the algorithms starts with the following observations.

• Our algorithms will overwrite the matrix \(A \text{,}\) and hence we introduce \(\widehat A \) to denote the original contents of \(A \text{.}\) We will say that the precondition for the algorithm is that

  \begin{equation*} A = \widehat A \end{equation*}

  (\(A \) starts by containing the original contents of \(A \text{.}\))

• We wish to overwrite \(A \) with \(L \) and \(U \text{.}\) Thus, the postcondition for the algorithm (the state in which we wish to exit the algorithm) is that

  \begin{equation*} A = L \backslash U \wedge L U = \widehat A \end{equation*}

  (\(A \) is overwritten by \(L \) below the diagonal and \(U \) on and above the diagonal, where multiplying \(L \) and \(U\) yields the original matrix \(A\text{.}\))

• All the algorithms will march through the matrices from top-left to bottom-right, giving us the code skeleton in Figure 5.5.1.2. Since the computed \(L \) and \(U \) overwrite \(A \text{,}\) throughout they are partitioned comformal to (in the same way) as is \(A \text{.}\)

• Thus, before and after each iteration of the loop the matrices are viewed as quadrants:

  \begin{equation*} A \rightarrow \FlaTwoByTwo { A_{TL}}{ A_{TR}} { A_{BL}}{ A_{BR}} , L \rightarrow \FlaTwoByTwo { L_{TL}}{ 0 } { L_{BL}}{ L_{BR}} , \mbox{and} \quad U \rightarrow \FlaTwoByTwo { U_{TL}}{ U_{TR} } { 0}{ U_{BR}}. \end{equation*}

  where \(A_{TL} \text{,}\) \(L_{TL} \text{,}\) and \(U_{TL} \) are all square and equally sized.

• In terms of these exposed quadrants, in the end we wish for matrix \(A \) to contain

  \begin{equation*} \begin{array}{l} \FlaTwoByTwo { A_{TL} }{ A_{TR} } { A_{BL} }{ A_{BR} } = \FlaTwoByTwo { L \backslash U_{TL} }{ U_{TR} } { L_{BL} }{ L\backslash U_{BR} } \\ ~~~~~\wedge \FlaTwoByTwo { L_{TL} }{0} { L_{BL} }{ L_{BR} } \FlaTwoByTwo { U_{TL} }{U_{TR}} { 0 }{ U_{BR} } = \FlaTwoByTwo { \widehat A_{TL} }{ \widehat A_{TR} } { \widehat A_{BL} }{ \widehat A_{BR} } \end{array} \end{equation*}

• Manipulating this yields what we call the Partitioned Matrix Expression (PME), which can be viewed as a recursive definition of the LU factorization:

  \begin{equation*} \begin{array}{l} \FlaTwoByTwo { A_{TL} }{ A_{TR} } { A_{BL} }{ A_{BR} } = \FlaTwoByTwo { L \backslash U_{TL} }{ U_{TR} } { L_{BL} }{ L \backslash U_{BR} } \\ ~~~~~\wedge \begin{array}{c | c} { L_{TL} U_{TL} = \widehat A_{TL} } \amp { L_{TL} U_{TR} = \widehat A_{TR} } \\ \hline { L_{BL} U_{TL} = \widehat A_{BL} } \amp { L_{BR} U_{BR} = \widehat A_{BR} - L_{BL} U_{TR} }. \end{array} \end{array} \end{equation*}

• Now, consider the code skeleton for the LU factorization in Figure 5.5.1.2. At the top of the loop (right after the \(\color{blue} {\bf while}\)), we want to maintain certain contents in matrix \(A \text{.}\) Since we are in a loop, we haven't yet overwritten \(A \) with the final result. Instead, some progress toward this final result has been made. The way we can find what the state of \(A \) is that we would like to maintain is to take the PME and delete subexpressions. For example, consider the following condition on the contents of \(A \text{:}\)

  \begin{equation*} \begin{array}{l} \FlaTwoByTwo { A_{TL} }{ A_{TR} } { A_{BL} }{ A_{BR} } = \FlaTwoByTwo { L \backslash U_{TL} }{ U_{TR} } { L_{BL} }{ \widehat A_{BR} - L_{BL} U_{TR} } \\ ~~~~~\wedge \begin{array}{c | c} { L_{TL} U_{TL} = \widehat A_{TL} } \amp { L_{TL} U_{TR} = \widehat A_{TR} } \\ \hline { L_{BL} U_{TL} = \widehat A_{BL} } \amp { \phantom{L_{BR} U_{BR} = \widehat A_{BR} - L_{BL} U_{TR}} }. \end{array} \end{array} \end{equation*}

  What we are saying is that \(A_{TL} \text{,}\) \(A_{TR} \text{,}\) and \(A_{BL} \) have been completely updated with the corresponding parts of \(L \) and \(U\text{,}\) and \(A_{BR}\) has been partially updated. This is exactly the state that the right-looking algorithm that we discussed in Subsection 5.2.2 maintains! What is left is to factor \(A_{BR} \text{,}\) since it contains \(\widehat A_{BR} - L_{BL} U_{TR} \text{,}\) and \(\widehat A_{BR} - L_{BL} U_{TR} = L_{BR} U_{BR} \text{.}\)

• By carefully analyzing the order in which computation must occur (in compiler lingo: by performing a dependence analysis), we can identify five states that can be maintained at the top of the loop, by deleting subexpressions from the PME. These are called loop invariants. There are five for LU factorization:

  \begin{equation*} \small \begin{array}{c} \begin{array}{c | c} \FlaTwoByTwo { A_{TL} }{ A_{TR} } { A_{BL} }{ A_{BR} } = \FlaTwoByTwo { L \backslash U_{TL} }{ \widehat A_{TR} } { \widehat A_{BL} }{ \widehat A_{BR} } \amp \FlaTwoByTwo { A_{TL} }{ A_{TR} } { A_{BL} }{ A_{BR} } = \FlaTwoByTwo { L \backslash U_{TL} }{ U_{TR} } { \widehat A_{BL} }{ \widehat A_{BR} } \\ \mbox{Invariant 1} \amp \mbox{Invariant 2} \\ \hline \FlaTwoByTwo { A_{TL} }{ A_{TR} } { A_{BL} }{ A_{BR} } = \FlaTwoByTwo { L \backslash U_{TL} }{ \widehat A_{TR} } { L_{BL} }{ \widehat A_{BR} } \amp \FlaTwoByTwo { A_{TL} }{ A_{TR} } { A_{BL} }{ A_{BR} } = \FlaTwoByTwo { L \backslash U_{TL} }{ U_{TR} } { L_{BL} }{ \widehat A_{BR} } \\ \mbox{Invariant 3} \amp \mbox{Invariant 4} \\ \hline \end{array} \\ \FlaTwoByTwo { A_{TL} }{ A_{TR} } { A_{BL} }{ A_{BR} } = \FlaTwoByTwo { L \backslash U_{TL} }{ U_{TR} } { L_{BL} }{ \widehat A_{BR}-L_{BL} U_{TR} } \\ \mbox{Invariant 5} \end{array} \end{equation*}

• Key to figuring out what updates must occur in the loop for each of the variants is to look at how the matrices are repartitioned at the top and bottom of the loop body.

For each of the five algorithms for LU factorization, we will derive the loop invariant, and then derive the algorithm from the loop invariant.