Accuracy of Computations

Gauss’s Method lends itself to computerization. The code below illustrates. It operates on an $n \times n$ matrix named a, doing row combinations using the first row, then the second row, etc.

for(row=1; row<=n-1; row++){
  for(row_below=row+1; row_below<=n; row_below++){
    multiplier=a[row_below,row]/a[row,row];
    for(col=row; col<=n; col++){
      a[row_below,col]-=multiplier*a[row,col];
    }
  }
}

This is in the C language. The for(row=1; row<=n-1; row++){ .. } loop initializes row at $1$ and then iterates while row is less than or equal to $n - 1$ , each time through incrementing row by one with the ++ operation. The other non-obvious language construct is that the -= in the innermost loop has the effect of a[row_below,col]=-1*multiplier*a[row,col]+a[row_below,col].

While that code is a first take on mechanizing Gauss’s Method, it is naive. For one thing, it assumes that the entry in the row,row position is nonzero. So one way that it needs to be extended is to cover the case where finding a zero in that location leads to a row swap or to the conclusion that the matrix is singular.

We could add some if statements to cover those cases but we will instead consider another way in which this code is naive. It is prone to pitfalls arising from the computer’s reliance on floating point arithmetic.

For example, above we have seen that we must handle a singular system as a separate case. But systems that are nearly singular also require care. Consider this one (the extra digits are in the ninth significant place).

$\begin{matrix} \begin{array}{rcrcr} x & + & 2 y & = & 3 \\ 1.000 000 01 x & + & 2 y & = & 3.000 000 01 \end{array} & (*) \end{matrix}$

By eye we easily spot the solution $x = 1$ , $y = 1$ . A computer has more trouble. If it represents real numbers to eight significant places, called single precision, then it will represent the second equation internally as $1.000 000 0 x + 2 y = 3.000 000 0$ , losing the digits in the ninth place. Instead of reporting the correct solution, this computer will think that the two equations are equal and it will report that the system is singular.

For some intuition about how the computer could come up with something that far off, consider this graph of the system.

Two almost indistinguishable descending lines are drawn on coordinate axes, with the point (1,1) marked. The diagram illustrates how nearly coincident equations can make the intersection sensitive to small coefficient changes.

We cannot tell the two lines apart; this system is nearly singular in the sense that the two lines are nearly the same line. This gives the system ( $*$ ) the property that a small change in an equation can cause a large change in the solution. For instance, changing the $3.000 000 01$ to $3.000 000 03$ changes the intersection point from $(1, 1)$ to $(3, 0)$ . The solution changes radically depending on the ninth digit, which explains why an eight-place computer has trouble. A problem that is very sensitive to inaccuracy or uncertainties in the input values is ill-conditioned.

The above example gives one way in which a system can be difficult to solve on a computer. It has the advantage that the picture of nearly-equal lines gives a memorable insight into one way for numerical difficulties to happen. Unfortunately this insight isn’t useful when we wish to solve some large system. We typically will not understand the geometry of an arbitrary large system.

There are other ways that a computer’s results may be unreliable, besides that the angle between some of the linear surfaces is small. For example, consider this system (from [Hamming]).

$\begin{matrix} \begin{array}{rcrcr} 0.001 x & + & y & = & 1 \\ x & - & y & = & 0 \end{array} & (* *) \end{matrix}$

The second equation gives $x = y$ , so $x = y = 1 / 1.001$ and thus both variables have values that are just less than $1$ . A computer using two digits represents the system internally in this way (we will do this example in two-digit floating point arithmetic for clarity but inventing a similar one with eight or more digits is easy).

$\begin{array}{rcrcr} (1.0 \times 10^{- 3}) \cdot x & + & (1.0 \times 10^{0}) \cdot y & = & 1.0 \times 10^{0} \\ (1.0 \times 10^{0}) \cdot x & - & (1.0 \times 10^{0}) \cdot y & = & 0.0 \times 10^{0} \end{array}$

The row reduction step $- 1000 ρ_{1} + ρ_{2}$ produces a second equation $- 1001 y = - 1000$ , which this computer rounds to two places as $(- 1.0 \times 10^{3}) y = - 1.0 \times 10^{3}$ . The computer decides from the second equation that $y = 1$ and with that it concludes from the first equation that $x = 0$ . The $y$ value is close but the $x$ is bad— the ratio of the actual answer to the computer’s answer is infinite. In short, another cause of unreliable output is the computer’s reliance on floating point arithmetic when the system-solving code leads to using leading entries that are small.

An experienced programmer may respond by using double precision, which retains sixteen significant digits, or perhaps using some even larger size. This will indeed solve many problems. However, double precision has greater memory requirements and besides we can obviously tweak the above to give the same trouble in the seventeenth digit, so double precision isn’t a panacea. We need a strategy to minimize numerical trouble as well as some guidance about how far we can trust the reported solutions.

A basic improvement on the naive code above is to not determine the factor to use for row combinations by simply taking the entry in the row,row position, but rather to look at all of the entries in the row column below the row,row entry and take one that is likely to give reliable results because it is not too small. This is partial pivoting.

For example, to solve the troublesome system ( $* *$ ) above we start by looking at both equations for a best entry to use, and take the $1$ in the second equation as more likely to give good results. The combination step of $- .001 ρ_{2} + ρ_{1}$ gives a first equation of $1.001 y = 1$ , which the computer will represent as $(1.0 \times 10^{0}) y = 1.0 \times 10^{0}$ , leading to the conclusion that $y = 1$ and, after back-substitution, that $x = 1$ , both of which are close to right. We can adapt the code from above to do this.

for(row=1; row<=n-1; row++){
/* find the largest entry in this column (in row max) */
  max=row;
  for(row_below=row+1; row_below<=n; row_below++){
    if (abs(a[row_below,row]) > abs(a[max,row]));
      max = row_below;
  }
/* swap rows to move that best entry up */
  for(col=row; col<=n; col++){
    temp=a[row,col];
    a[row,col]=a[max,col];
    a[max,col]=temp;
  }
/* proceed as before */
  for(row_below=row+1; row_below<=n; row_below++){
    multiplier=a[row_below,row]/a[row,row];
     for(col=row; col<=n; col++){
       a[row_below,col]-=multiplier*a[row,col];
    }
  }
}

A full analysis of the best way to implement Gauss’s Method is beyond the scope of this book (see [Wilkinson 1965]), but the method recommended by most experts first finds the best entry among the candidates and then scales it to a number that is less likely to give trouble. This is scaled partial pivoting.

In addition to returning a result that is likely to be reliable, most well-done code will return a conditioning number that describes the factor by which uncertainties in the input numbers could be magnified to become inaccuracies in the results returned (see [Rice]).

The lesson is that just because Gauss’s Method always works in theory, and just because computer code correctly implements that method, doesn’t mean that the answer is reliable. In practice, always use a package where experts have worked hard to counter what can go wrong.

Exercises

Exercise 1 Worked answer

Using two decimal places, add $253$ and $2 / 3$ .

Back to Exercise 1

Answer. Scientific notation is convenient for expressing the two-place restriction: $.25 \times 10^{2} + .67 \times 10^{0} = .25 \times 10^{2}$ . Note that adding the $2 / 3$ has no effect on the total.
Exercise 2 Worked answer

This intersect-the-lines problem contrasts with the example discussed above.

$\begin{array}{rcrcr} x & + & 2 y & = & 3 \\ 3 x & - & 2 y & = & 1 \end{array}$

Illustrate that in this system some small change in the numbers will produce only a small change in the solution by changing the constant in the bottom equation to $1.008$ and solving. Compare it to the solution of the unchanged system.

Back to Exercise 2

Answer. The reduction

$⟶_{\begin{array}{l} \end{array}}^{- 3 ρ_{1} + ρ_{2}} \begin{array}{rcrcr} x & + & 2 y & = & 3 \\ - 8 y & = & - 7.992 \end{array}$

gives $(x, y) = (1.002, 0.999)$ . So for this system a small change in the constant produces only a small change in the solution.
Exercise 3 Worked answer

Consider this system ([Rice]).

$\begin{array}{rcrcr} 0.000 3 x & + & 1.556 y & = & 1.559 \\ 0.345 4 x & - & 2.346 y & = & 1.108 \end{array}$
1. Solve it.
2. Solve it by rounding at each step to four digits.
Back to Exercise 3

Answer.
1. The fully accurate solution is that $x = 10$ and $y = 1$ .
2. The four-digit reduction
  
  $⟶_{\begin{array}{l} \end{array}}^{- (.3454 / .0003) ρ_{1} + ρ_{2}} (\begin{array}{ccc} .0003 & 1.556 & 1.569 \\ 0 & - 1794 & - 1805 \end{array})$
  
  gives the conclusion that $y = 1.006$ , which is not bad, and that $x = 12.21$ . Of course, this is twenty percent different than the correct answer.
Exercise 4 Worked answer

Rounding inside the computer often has an effect on the result. Assume that your machine has eight significant digits.
1. Show that the machine will compute $(2 / 3) + ((2 / 3) - (1 / 3))$ as unequal to $((2 / 3) + (2 / 3)) - (1 / 3)$ . Thus, computer arithmetic is not associative.
2. Compare the computer’s version of $(1 / 3) x + y = 0$ and $(2 / 3) x + 2 y = 0$ . Is twice the first equation the same as the second?
Back to Exercise 4

Answer.
1. For the first one, first, $(2 / 3) - (1 / 3)$ is $.666 666 67 - .333 333 33 = .333 333 34$ and so $(2 / 3) + ((2 / 3) - (1 / 3)) = .666 666 67 + .333 333 34 = 1.000 000 0$ .
  
  For the other one, first $((2 / 3) + (2 / 3)) = .666 666 67 + .666 666 67 = 1.333 333 3$ and so $((2 / 3) + (2 / 3)) - (1 / 3) = 1.333 333 3 - .333 333 33 = .999 999 97$ .
2. The first equation is $.333 333 33 \cdot x + 1.000 000 0 \cdot y = 0$ while the second is $.666 666 67 \cdot x + 2.000 000 0 \cdot y = 0$ .
Exercise 5 Worked answer

Ill-conditioning is not only dependent on the matrix of coefficients. This example [Hamming] shows that it can arise from an interaction between the left and right sides of the system. Let $ε$ be a small real.

$\begin{array}{rcrcrcr} 3 x & + & 2 y & + & z & = & 6 \\ 2 x & + & 2 ε y & + & 2 ε z & = & 2 + 4 ε \\ x & + & 2 ε y & - & ε z & = & 1 + ε \end{array}$
1. Solve the system by hand. Notice that the $ε$ ’s divide out only because there is an exact cancellation of the integer parts on the right side as well as on the left.
2. Solve the system by hand, rounding to two decimal places, and with $ε = 0.001$ .
Back to Exercise 5

Answer.
1. This calculation
  
  $\begin{aligned} ⟶_{\begin{matrix} - (1 / 3) ρ_{1} + ρ_{3} \end{matrix}}^{- (2 / 3) ρ_{1} + ρ_{2}} (\begin{array}{cccc} 3 & 2 & 1 & 6 \\ 0 & - (4 / 3) + 2 ε & - (2 / 3) + 2 ε & - 2 + 4 ε \\ 0 & - (2 / 3) + 2 ε & - (1 / 3) - ε & - 1 + ε \end{array}) \\ ⟶_{\begin{matrix} \end{matrix}}^{- (1 / 2) ρ_{2} + ρ_{3}} (\begin{array}{cccc} 3 & 2 & 1 & 6 \\ 0 & - (4 / 3) + 2 ε & - (2 / 3) + 2 ε & - 2 + 4 ε \\ 0 & ε & - 2 ε & - ε \end{array}) \end{aligned}$
  
  gives a third equation of $y - 2 z = - 1$ . Substituting into the second equation gives $((- 10 / 3) + 6 ε) \cdot z = (- 10 / 3) + 6 ε$ so $z = 1$ and thus $y = 1$ . With those, the first equation says that $x = 1$ .
2. As above, scientific notation is convenient to express the restriction on the numbe of digits.
  
  The solution with two digits retained is $z = 2.1$ , $y = 2.6$ , and $x = - .43$ .
  
  $\begin{matrix} (\begin{array}{cccc} .30 \times 10^{1} & .20 \times 10^{1} & .10 \times 10^{1} & .60 \times 10^{1} \\ .10 \times 10^{1} & .20 \times 10^{- 3} & .20 \times 10^{- 3} & .20 \times 10^{1} \\ .30 \times 10^{1} & .20 \times 10^{- 3} & - .10 \times 10^{- 3} & .10 \times 10^{1} \end{array}) \\ \begin{aligned} ⟶_{\begin{matrix} - (1 / 3) ρ_{1} + ρ_{3} \end{matrix}}^{- (2 / 3) ρ_{1} + ρ_{2}} (\begin{array}{cccc} .30 \times 10^{1} & .20 \times 10^{1} & .10 \times 10^{1} & .60 \times 10^{1} \\ 0 & - .13 \times 10^{1} & - .67 \times 10^{0} & - .20 \times 10^{1} \\ 0 & - .67 \times 10^{0} & - .33 \times 10^{0} & - .10 \times 10^{1} \end{array}) \\ ⟶_{\begin{matrix} \end{matrix}}^{- (.67 / 1.3) ρ_{2} + ρ_{3}} (\begin{array}{cccc} .30 \times 10^{1} & .20 \times 10^{1} & .10 \times 10^{1} & .60 \times 10^{1} \\ 0 & - .13 \times 10^{1} & - .67 \times 10^{0} & - .20 \times 10^{1} \\ 0 & 0 & .15 \times 10^{- 2} & .31 \times 10^{- 2} \end{array}) \end{aligned} \end{matrix}$

References cited in this section

Hamming

Richard W. Hamming, Introduction to Applied Numerical Analysis, Hemisphere Publishing, 1971.

Wilkinson 1965

The Algebraic Eigenvalue Problem, J. H. Wilkinson, Oxford University Press, 1965

Rice

John R. Rice, Numerical Mathods, Software, and Analysis, second edition, Academic Press, 1993.