3.3: Partial Pivoting

When performing Gaussian elimination, the diagonal element that one uses during the elimination procedure is called the pivot. To obtain the correct multiple, one uses the pivot as the divisor to the elements below the pivot. Gaussian elimination in this form will fail if the pivot is zero. In this situation, a row interchange must be performed.

Even if the pivot is not identically zero, a small value can result in big roundoff errors. For very large matrices, one can easily lose all accuracy in the solution. To avoid these round-off errors arising from small pivots, row interchanges are made, and this technique is called partial pivoting (partial pivoting is in contrast to complete pivoting, where both rows and columns are interchanged). We will illustrate by example the LU decomposition using partial pivoting. Consider

\[A=\left(\begin{array}{rrr} -2 & 2 & -1 \\ 6 & -6 & 7 \\ 3 & -8 & 4 \end{array}\right) \nonumber \]

We interchange rows to place the largest element (in absolute value) in the pivot, or \(a_{11}\), position. That is,

\[\mathrm{A} \rightarrow\left(\begin{array}{rrr} 6 & -6 & 7 \\ -2 & 2 & -1 \\ 3 & -8 & 4 \end{array}\right)=\mathrm{P}_{12} \mathrm{~A} \nonumber \]

where

\[P_{12}=\left(\begin{array}{lll} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right) \nonumber \]

is a permutation matrix that when multiplied on the left interchanges the first and second rows of a matrix. Note that \(P_{12}^{-1}=P_{12}\). The elimination step is then

\[\mathrm{P}_{12} \mathrm{~A} \rightarrow\left(\begin{array}{rrc} 6 & -6 & 7 \\ 0 & 0 & 4 / 3 \\ 0 & -5 & 1 / 2 \end{array}\right)=\mathrm{M}_{1} \mathrm{P}_{12} \mathrm{~A} \nonumber \]

where

\[M_{1}=\left(\begin{array}{ccc} 1 & 0 & 0 \\ 1 / 3 & 1 & 0 \\ -1 / 2 & 0 & 1 \end{array}\right) \nonumber \]

The final step requires one more row interchange:

\[\mathrm{M}_{1} \mathrm{P}_{12} \mathrm{~A} \rightarrow\left(\begin{array}{rrc} 6 & -6 & 7 \\ 0 & -5 & 1 / 2 \\ 0 & 0 & 4 / 3 \end{array}\right)=\mathrm{P}_{23} \mathrm{M}_{1} \mathrm{P}_{12} \mathrm{~A}=\mathrm{U} \nonumber \]

Since the permutation matrices given by P are their own inverses, we can write our result as

\[\left(P_{23} M_{1} P_{23}\right) P_{23} P_{12} A=U \nonumber \]

Multiplication of \(M\) on the left by \(P\) interchanges rows while multiplication on the right by \(P\) interchanges columns. That is,

\[P_{23}\left(\begin{array}{cll} 1 & 0 & 0 \\ 1 / 3 & 1 & 0 \\ -1 / 2 & 0 & 1 \end{array}\right) P_{23}=\left(\begin{array}{ccc} 1 & 0 & 0 \\ -1 / 2 & 0 & 1 \\ 1 / 3 & 1 & 0 \end{array}\right) P_{23}=\left(\begin{array}{ccc} 1 & 0 & 0 \\ -1 / 2 & 1 & 0 \\ 1 / 3 & 0 & 1 \end{array}\right) \nonumber \]

The net result on \(M_{1}\) is an interchange of the nondiagonal elements \(1 / 3\) and \(-1 / 2\).

We can then multiply by the inverse of \(\left(\mathrm{P}_{23} \mathrm{M}_{1} \mathrm{P}_{23}\right)\) to obtain

\[P_{23} P_{12} A=\left(P_{23} M_{1} P_{23}\right)^{-1} U \nonumber \]

which we write as

\[\mathrm{PA}=\mathrm{LU} \nonumber \]

Instead of L, MATLAB will write this as

\[\mathrm{A}=\left(\mathrm{P}^{-1} \mathrm{~L}\right) \mathrm{U} \nonumber \]

For convenience, we will just denote \(\left(\mathrm{P}^{-1} \mathrm{~L}\right)\) by \(\mathrm{L}\), but by \(\mathrm{L}\) here we mean a permutated lower triangular matrix.

For example, in MATLAB, to solve \(A \mathbf{x}=\mathbf{b}\) for \(\mathbf{x}\) using Gaussian elimination, one types

\[\mathrm{x}=\mathrm{A} \backslash \mathrm{b} \nonumber \]

where \solves for \(\mathbf{x}\) using the most efficient algorithm available, depending on the form of \(\mathrm{A}\). If \(\mathrm{A}\) is a general \(n \times n\) matrix, then first the LU decomposition of \(\mathrm{A}\) is found using partial pivoting, and then \(x\) is determined from permuted forward and backward substitution. If \(A\) is upper or lower triangular, then forward or backward substitution (or their permuted version) is used directly.

If there are many different right-hand-sides, one would first directly find the LU decomposition of A using a function call, and then solve using . That is, one would iterate for different \(\mathbf{b}^{\prime}\) s the following expressions:

\[\begin{aligned} &{[\mathrm{LU}]=l u(\mathrm{~A})} \\ &\mathrm{y}=\mathrm{L} \backslash \mathrm{b} \\ &\mathrm{x}=\mathrm{U} \backslash y \end{aligned} \nonumber \]

where the second and third lines can be shortened to

\[\mathrm{x}=\mathrm{U} \backslash(\mathrm{L} \backslash \mathrm{b}) \nonumber \]

where the parenthesis are required. In lecture, I will demonstrate these solutions in MATLAB using the matrix \(A=[-2,2,-1 ; 6,-6,7 ; 3,-8,4] ;\) which is the example in the notes.