ALAFF LU factorization with partial pivoting

Subsection 5.3.3 LU factorization with partial pivoting

fit width

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

$\begin{equation*} \left[ A, p \right] \becomes \mbox{LUpiv}{ A }, \end{equation*}$

🔗

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

$\begin{equation*} \maxi( x ) \end{equation*}$

🔗

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

$\begin{equation*} \FlaThreeByThreeBR {A_{00}}{a_{01}}{A_{02}} {0}{\alpha_{11}}{a_{12}^T} {0}{a_{21}}{A_{22}} \end{equation*}$

🔗

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

$\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*}$

🔗

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.

🔗

$\begin{equation*} \newcommand{\routinename}{[ A, p ] = \mbox{LUpiv-right-looking}( A )} \newcommand{\guard}{ n( A_{TL} ) \lt n( A ) } \newcommand{\partitionings}{ A \rightarrow \FlaTwoByTwo{A_{TL}}{A_{TR}} {A_{BL}}{A_{BR}}, p \rightarrow \FlaTwoByOne{ p_T }{ p_B } } \newcommand{\partitionsizes}{ A_{TL} {\rm ~is~} 0 \times 0 , p_T \mbox{ has } 0 \mbox{ elements} } \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}}, \FlaTwoByOne{ p_T }{ p_B } \rightarrow \FlaThreeByOneB{ p_0 }{ \pi_1 }{ p_2 } } \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}}, \FlaTwoByOne{ p_T }{ p_B } \leftarrow \FlaThreeByOneT{ p_0 }{ \pi_1 }{ p_2 } } \newcommand{\update}{ \begin{array}{l} \pi_1 := \maxi\left( \begin{array}{c} \alpha_{11} \\ \hline a_{21} \end{array} \right) \\ \FlaThreeByThreeBR{A_{00}}{a_{01}}{A_{02}} {a_{10}^T}{ \alpha_{11} }{ a_{12}^T } {A_{20}}{ a_{21} }{ A_{22} } \becomes \FlaTwoByTwo {I}{0} {0}{P( \pi_1 )} \FlaThreeByThreeBR{A_{00}}{a_{01}}{A_{02}} {a_{10}^T}{ \alpha_{11} }{ a_{12}^T } {A_{20}}{ a_{21} }{ A_{22} } \\ a_{21} := a_{21} / \alpha_{11} \\ A_{22} := A_{22} - a_{21} a_{12}^T \end{array} } \FlaAlgorithm \end{equation*}$

Figure 5.3.3.1. Right-looking LU factorization algorithm with partial pivoting.

🔗

fit width

🔗

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

$\begin{equation*} L_{n-1} P_{n-1} \cdots L_0 P_0 A = U \end{equation*}$

🔗

or, equivalently,

$\begin{equation*} A = P_0^T L_{0}^{-1} \cdots P_{n-1}^T L_{n-1}^{-1} U . \end{equation*}$

🔗

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

$\begin{equation*} A = P_0 L_{0}^{-1} \cdots P_{n-1} L_{n-1}^{-1} U . \end{equation*}$

🔗

What we will finally show is that there are Gauss transforms

${L}_0^\star, \ldots {L}_{n-1}^\star$ such that

$\begin{equation*} A = P_0 \cdots P_{n-1} \begin{array}[t]{c} \underbrace{{L}_{0}^\star \cdots {L}_{n-1}^\star} \\ L \end{array} U \end{equation*}$

🔗

or, equivalently,

$\begin{equation*} \widetilde P( p ) A = P_{n-1} \cdots P_0 A = \begin{array}[t]{c} \underbrace{{L}_{0}^\star \cdots {L}_{n-1}^\star} \\ L \end{array} U, \end{equation*}$

🔗

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

$\begin{equation*} \left( \begin{array}{c c } \alpha_{11} \amp a_{12}^T \\ a_{21} \amp A_{22} \end{array} \right) \end{equation*}$

🔗

but also all of

$\begin{equation*} \left( \begin{array}{c |c c} a_{10}^T \amp \alpha_{11} \amp a_{12}^T \\ A_{20} \amp a_{21} \amp A_{22} \end{array} \right), \end{equation*}$

🔗

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{.}$

🔗

Homework 5.3.3.1.

Implement the algorithm given in Figure 5.3.3.1 as

function [ A_out ] = LUpiv_right_looking( A )

by completing the code in Assignments/Week05/matlab/LUpiv_right_looking.m. Input is an

$m \times n$ matrix

$A\text{.}$ Output is the matrix

$A$ that has been overwritten by the LU factorization and pivot vector

$p \text{.}$ You may want to use Assignments/Week05/matlab/test_LUpiv_right_looking.m to check your implementation.

The following utility functions may come in handy:

which we hope are self explanatory.

Solution

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