Skip to main content

Section 4.3 Pivoting

The elementary implementation of Gaussian Elimination in Section 4.2 suffers two main problems:
  1. It fails if the element A[j,j] is zero when we evaluate L[i,j].
  2. There might be a loss of significant digits in the difference A[i,k]-L[i,j]*A[j,k] if L[i,j] is very big.

For instance, let us apply the algorithm to
\begin{equation*} A=\begin{pmatrix} \epsilon&1\\ 1&1\\ \end{pmatrix}. \end{equation*}
After the first step, the element A[2,2] becomes \(1-1/\epsilon\text{.}\) If \(\epsilon\) is small enough, then \(1/\epsilon\) is so much larger than 1 that \(1-1/\epsilon=-1/\epsilon\,\text{,}\) namely we get the same decomposition as if the original matrix were
\begin{equation*} \begin{pmatrix} \epsilon&1\\ 1&0\\ \end{pmatrix}. \end{equation*}
A concrete example. Let us make concrete calculations in an example in the floating point system \(D_3\) with base 10 and three significant digits that we introduced in Section 1.2.

Consider the system
\begin{equation*} Ax=b \end{equation*}
\begin{equation*} \begin{pmatrix} 0.00025&1\\ 1&1\\ \end{pmatrix} \begin{pmatrix} x\\ y\\ \end{pmatrix} = \begin{pmatrix} 1\\ 2\\ \end{pmatrix} \end{equation*}

Its exact solution is \(x=\frac{4000}{3999},y=\frac{3998}{3999}\text{.}\) Notice that these numbers are both equal to \(1.00\) in \(D_3.\)

Using the LU decomposition, though, we get
\begin{equation*} L = \begin{pmatrix} \phantom{400}1&0\\ 4000&1\\ \end{pmatrix} \end{equation*}
and
\begin{equation*} U = \begin{pmatrix} 0.00025&1\\ 0&1-4000\\ \end{pmatrix} = \begin{pmatrix} 0.00025&1\\ 0&-4000\\ \end{pmatrix} \end{equation*}
and the right-hand-side becomes
\begin{equation*} L^{-1}b = \begin{pmatrix}\phantom{-400}1&0\\ -4000&1\\ \end{pmatrix}\begin{pmatrix}1\\ 2\\ \end{pmatrix} =\begin{pmatrix}1\\ 2-4000\\ \end{pmatrix}=\begin{pmatrix}1\\ -4000\\ \end{pmatrix} \end{equation*}
(recall that we are in a 3-digits system!)

In other words, the system is equivalent to
\begin{equation*} \left\{ \begin{array}{c} 0.00025x+y=1\\ -4000y=-4000\\ \end{array} \right. \end{equation*}

so that \(y=1.00\) and therefore \(x=0.00\text{.}\)

Notice that we got a relative error of 100% in \(x\text{!!!}\)

Pivoting. This kind of problems can be fixed by rearranging at every step the rows of \(A\) so that, when we evaluate L[i,j], we always divide by the largest possible number.

We illustrate this method by applying it to the system of the previous example:
\begin{equation*} \begin{pmatrix} 0.00025&1\\ 1&1\\ \end{pmatrix}. \begin{pmatrix} x\\ y\\ \end{pmatrix} = \begin{pmatrix} 1\\ 2\\ \end{pmatrix} \end{equation*}
As discussed above, dividing by \(0.00025\) is not a good idea, so we re-write it as
\begin{equation*} \begin{pmatrix} 1&1\\ 0.00025&1\\ \end{pmatrix}. \begin{pmatrix} x\\ y\\ \end{pmatrix} = \begin{pmatrix} 2\\ 1\\ \end{pmatrix} \end{equation*}
Clearly nothing has changed: the solution of a system does not depend on the order we use to write its equations!

Now we use LU decomposition:
\begin{equation*} L= \begin{pmatrix} 1&0\\ 0.00025&1\\ \end{pmatrix},\; U= \begin{pmatrix} 1&1\\ 0&1\\ \end{pmatrix} \end{equation*}
In other words, the system becomes
\begin{equation*} \begin{pmatrix} 1&1\\ 0&1\\ \end{pmatrix}. \begin{pmatrix} x\\ y\\ \end{pmatrix} = L^{-1} \begin{pmatrix} 2\\ 1\\ \end{pmatrix} = \begin{pmatrix} 2\\ 1\\ \end{pmatrix} \end{equation*}
(remember that we are in a 3 digits floating-point system!)

whose solution is clearly \(x=y=1.00\text{.}\)

The error we experienced before pivoting disappeared completely!

Pivoting code. Below is an implementation of the LU method with pivoting. As in the last implementation of the LU method in , we use a single for loop. Notice that, since while pivoting we need to exchange rows in P, A and L, we need to define from the beginning all matrices. We define L to be the zero matrix and P to be the unit matrix.