ALAFF Gaussian elimination with row exchanges

Subsection 5.3.1 Gaussian elimination with row exchanges

Homework 5.3.1.1.

Perform Gaussian elimination as explained in Subsection 5.2.1 to solve

\begin{equation*} \left( \begin{array}{c c} 0 \amp 1 \\ 1 \amp 0 \end{array} \right) \left( \begin{array}{c} \chi_0 \\ \chi_1 \end{array} \right) = \left( \begin{array}{c} 2 \\ 1 \end{array} \right) \end{equation*}

Solution

The appended system is given by

\begin{equation*} \left( \begin{array}{c c | c} 0 \amp 1 \amp 2\\ 1 \amp 0 \amp 1 \end{array} \right) . \end{equation*}

In the first step, the multiplier is computed as \(\lambda_{1,0} = 1/0 \) and the algorithm fails. Yet, it is clear that the (unique) solution is

\begin{equation*} \left( \begin{array}{c} \chi_0 \\ \chi_1 \end{array} \right) = \left( \begin{array}{c} 1 \\ 2 \end{array} \right). \end{equation*}

The point of the exercise: Gaussian elemination and, equivalently, LU factorization as we have discussed so far can fail if a "divide by zero" is encountered. The element on the diagonal used to compute the multipliers in a current iteration of the outer-most loop is called the pivot (element). Thus, if a zero pivot is encountered, the algorithms fail. Even if the pivot is merely small (in magnitude), as we will discuss in a future week, roundoff error encountered when performing floating point operations will likely make the computation "numerically unstable," which is the topic of next week's material.

The simple observation is that the rows of the matrix (and corresponding right-hand side element) correspond to linear equations that must be simultaneously solved. Reordering these does not change the solution. Reordering in advance so that no zero pivot is encountered is problematic, since pivots are generally updated by prior computation. However, when a zero pivot is encountered, the row in which it appears can simply be swapped with another row so that the pivot is replaced with a nonzero element (which then becomes the pivot). In exact arithmetic, it suffices to ensure that the pivot is nonzero after swapping. As mentioned, in the presence of roundoff error, any element that is small in magnitude can create problems. For this reason, we will swap rows so that the element with the largest magnitude (among the elements in the "current" column below the diagonal) becomes the pivot. This is known as partial pivoting or row pivoting.

Homework 5.3.1.2.

When performing Gaussian elimination as explained in Subsection 5.2.1 to solve

\begin{equation*} \left( \begin{array}{c c} 10^{-k} \amp 1 \\ 1 \amp 0 \end{array} \right) \left( \begin{array}{c} \chi_0 \\ \chi_1 \end{array} \right) = \left( \begin{array}{c} 1 \\ 1 \end{array} \right), \end{equation*}

set

\begin{equation*} 1 - 10^k \end{equation*}

\begin{equation*} - 10^k \end{equation*}

(since we will assume \(k \) to be large and hence \(1 \) is very small to relative to \(10^k \)). With this modification (which simulates roundoff error that may be encountered when performing floating point computation), what is the answer?

Next, solve

\begin{equation*} \left( \begin{array}{c c} 1 \amp 0 \\ 10^{-k} \amp 1 \end{array} \right) \left( \begin{array}{c} \chi_0 \\ \chi_1 \end{array} \right) = \left( \begin{array}{c} 1 \\ 1 \end{array} \right). \end{equation*}

What do you observe?

Solution

The appended system is given by

\begin{equation*} \left( \begin{array}{c c | c} 10^{-k} \amp 1 \amp 1\\ 1 \amp 0 \amp 1 \end{array} \right) . \end{equation*}

In the first step, the multiplier is computed as \(\lambda_{1,0} = 10^k \) and the updated appended system becomes

\begin{equation*} \left( \begin{array}{c c | c} 10^{-k} \amp 1 \amp 1\\ 0 \amp - 10^k \amp 1 - 10^k \end{array} \right) \end{equation*}

which is rounded to

\begin{equation*} \left( \begin{array}{c c | c} 10^{-k} \amp 1 \amp 1\\ 0 \amp - 10^k \amp -10^k \end{array} \right) . \end{equation*}

We then compute

\begin{equation*} \chi_1 = (-10^k) / (-10^k) = 1 \end{equation*}

and

\begin{equation*} \chi_0 = (1 - \chi_1) / 10^{-k} = (1-1)/ 10^{-k} = 0. \end{equation*}

If we instead start with the equivalent system

\begin{equation*} \left( \begin{array}{c c | c} 1 \amp 0 \amp 1 \\ 10^{-k} \amp 1 \amp 1 \end{array} \right) . \end{equation*}

the appended system after one step becomes

\begin{equation*} \left( \begin{array}{c c | c} 1 \amp 0 \amp 1 \\ 0 \amp 1 \amp 1 - 10^{-k} \end{array} \right) \end{equation*}

which yields the solution

\begin{equation*} \left( \begin{array}{c} \chi_0 \\ \chi_1 \end{array} \right) = \left( \begin{array}{c} 1 \\ 1 - 10^{-k} \end{array} \right). \end{equation*}

which becomes

\begin{equation*} \left( \begin{array}{c} \chi_0 \\ \chi_1 \end{array} \right) = \left( \begin{array}{c} 1 \\ 1 \end{array} \right). \end{equation*}

as \(k \) gets large.

What this illustrates is how a large multiple of a row being added to another row can wipe out information in that second row. After one step of Gaussian elimination, the system becomes equivalent to one that started with

\begin{equation*} \left( \begin{array}{c c | c} 10^{-k} \amp 1 \amp 1\\ 1 \amp 0 \amp 0 \end{array} \right) . \end{equation*}