Advanced Linear Algebra: Foundations to Frontiers

Subsection 5.3.3 LU factorization with partial pivoting

Having introduced our notation for permutation matrices, we can now define the LU factorization with partial pivoting: Given an \(m \times n \) matrix \(A \text{,}\) we wish to compute

• vector \(p \) of \(n \) integers that indicates how rows are pivoting as the algorithm proceeds,

• a unit lower trapezoidal matrix \(L \text{,}\) and

• an upper triangular matrix \(U \)

so that \(\widetilde P( p ) A = L U \text{.}\) We represent this operation by

where upon completion \(A \) has been overwritten by \(\{ L \backslash U \} \text{,}\) which indicates that \(U \) overwrites the upper triangular part of \(A \) and \(L \) is stored in the strictly lower triangular part of \(A \text{.}\)

Let us start with revisiting the derivation of the right-looking LU factorization in Subsection 5.2.2. The first step is to find a first permutation matrix \(\widetilde P( \pi_1 ) \) such that the element on the diagonal in the first column is maximal in value. (Mathematically, any nonzero value works. We will see that ensuring that the multiplier is less than one in magnitude reduces the potential for accumulation of error.) For this, we will introduce the function

which, given a vector \(x \text{,}\) returns the index of the element in \(x \) with maximal magnitude (absolute value). The algorithm then proceeds as follows:

• Partition \(A \text{,}\) \(L \) as follows:

  \begin{equation*} A \rightarrow \FlaTwoByTwoSingleLine { \alpha_{11} }{ a_{12}^T } { a_{21} }{ A_{22} }, \quad \mbox{and} \quad L \rightarrow \FlaTwoByTwoSingleLine { 1 }{ 0 } {l_{21}}{L_{22}} \end{equation*}

• Compute \(\pi_1 = \maxi\left( \begin{array}{c} \alpha_{11} \\ \hline a_{21} \end{array} \right) \text{.}\)

• Permute the rows: \(\FlaTwoByTwoSingleLine { \alpha_{11} }{ a_{12}^T } { a_{21} }{ A_{22} } \becomes \widetilde P( \pi_1 ) \FlaTwoByTwoSingleLine { \alpha_{11} }{ a_{12}^T } { a_{21} }{ A_{22} } \text{.}\)

• Compute \(l_{21} \becomes a_{21} / \alpha_{11} \text{.}\)

• Update \(A_{22} \becomes A_{22} - l_{21} a_{12}^T \text{.}\)

This completes the introduction of zeroes below the diagonal of the first column.

Now, more generally, assume that the computation has proceeded to the point where matrix \(A \) has been overwritten by

where \(A_{00} \) is upper triangular. If no pivoting was added one would compute \(l_{21} \becomes a_{21} / \alpha_{11} \) followed by the update

Now, instead one performs the steps

• Compute

  \begin{equation*} \pi_1 := \maxi\left( \begin{array}{c} \alpha_{11} \\ \hline a_{21} \end{array} \right). \end{equation*}

• Permute the rows:

  \begin{equation*} \FlaThreeByThreeBR{A_{00}}{a_{01}}{A_{02}} {0}{ \alpha_{11} }{ a_{12}^T } {0}{ a_{21} }{ A_{22} } \becomes \FlaTwoByTwo {I}{0} {0}{\widetilde P( \pi_1 )} \FlaThreeByThreeBR{A_{00}}{a_{01}}{A_{02}} {0}{ \alpha_{11} }{ a_{12}^T } {0}{ a_{21} }{ A_{22} } \end{equation*}

• Update

  \begin{equation*} l_{21} \becomes a_{21} / \alpha_{11}\text{.} \end{equation*}

• Update

  \begin{equation*} \begin{array}{l} \FlaThreeByThreeBR {A_{00}}{a_{01}}{A_{02}} {0}{\alpha_{11}}{a_{12}^T} {0}{a_{21}}{A_{22}} \becomes \FlaThreeByThreeBR {I}{0}{0} {0}{1}{0} {0}{-l_{21}}{I} \FlaThreeByThreeBR {A_{00}}{a_{01}}{A_{02}} {0}{\alpha_{11}}{a_{12}^T} {0}{a_{21}}{A_{22}}\\ ~~~~~= \FlaThreeByThreeBR {A_{00}}{a_{01}}{A_{02}} {0}{\alpha_{11}}{a_{12}^T} {0}{0}{A_{22} - l_{21} a_{12}^T}. \end{array} \end{equation*}

This algorithm is summarized in Figure 5.3.3.1. In that algorithm, the lower triangular matrix \(L \) is accumulated below the diagonal.

What this algorithm computes is a sequence of Gauss transforms \(L_{0}, \ldots , L_{n-1} \) and permutations \(P_0, \ldots, P_{n-1} \) such that

or, equivalently,

Actually, since \(P_k = \left( \begin{array}{c | c} I_{k \times k} \amp 0 \\ \hline 0 \amp \widetilde P( \pi ) \end{array} \right)\) for some \(\pi \) , we know that \(P_k^T = P_k \) and hence

What we will finally show is that there are Gauss transforms \({L}_0^\star, \ldots {L}_{n-1}^\star \) such that

or, equivalently,

which is what we set out to compute.

Here is the insight. If only we know how to order the rows of \(A \) and right-hand side \(b \) correctly, then we would not have to pivot. But we only know how to pivot as the computation unfolds. Recall that the multipliers can overwrite the elements they zero in Gaussian elimination and do so when we formulate it as an LU factorization. By not only pivoting the elements of

but also all of

we are moving the computed multipliers with the rows that are being swapped. It is for this reason that we end up computing the LU factorization of the permuted matrix: \(\widetilde P( p ) A \text{.}\)

function [ A_out ] = LUpiv_right_looking( A )