Math 236  •  6 March 2001

Partial pivoting in Gaussian elimination

This notebook illustrates how the straightforward Gauss-Jordan algorithm can produce highly inaccurate results due to roundoff errors in the intermediate steps.  And it shows how to ameliorate such roundoff errors by the method of partial pivoting.

Functions for doing elementary row operations

swap::usage  =  "swap[mat, i, k] swaps the ith and kth rows of \
matrix mat.";
swap[mat_?MatrixQ,  i_Integer,  k_Integer]  :=  
      Module[{A  =  mat},  
        A[[{i, k}]]  =  
          A[[{k, i}]];  A]

scale::usage  =  "scale[mat, c, i] multiplies the ith row of \
matrix mat by scalar c.";
scale[mat_?MatrixQ,  c_?NumericQ,  i_Integer]  :=  
      Module[{A  =  mat},  
        A[[i]]  =  
          c A[[i]];  A]

addrow::usage  =  "addrow[mat, c, k, i] replaces the ith row of \
matrix mat by the sum of that row and scalar c times the kth row.";
addrow[mat_List,  c_?NumericQ,  k_Integer,  i_Integer]  :=  
      Module[{A  =  mat},  
        A[[i]]  =  A[[i]]  +  c A[[k]]; A]

The linear system

The linear system to be solved here is:
	0.001 (Subscript[x, 1])  +  (Subscript[x, 2])  =  
        1, 
               (Subscript[x, 1]) +  
          (Subscript[x, 2])  =  2. 
The augmented matrix of this system is:
	A = (({{0.001, 1, 1}, {1, 1, 2}})).

Gauss-Jordan algorithm with exact arithmetic

A = (({{(1/(10^3)), 1, 1}, {1, 1, 2}}));

Carry out Gauss-Jordan algorithm with exact, rational arithmetic

A  =  
      scale[A,  1/A[[1, 1]],  1];  
    A  //  MatrixForm

A  =  
      addrow[A,  -A[[2, 1]]/
          A[[1, 1]],  1,  2];  
    A  //  MatrixForm

A  =  
      scale[A,  1/A[[2, 2]],  2];  
    A   //  MatrixForm

A  =  
      addrow[A,  -A[[1, 2]]/
          A[[2, 2]],  2,  1];  
    A  //  MatrixForm

Get numerical approximation of answer

A  =  N[A];  A  //  MatrixForm

To display more digits, change the default print routine so as to display numbers with full machine precision of 16 digits:

$PrePrint  =  NumberForm[#,  16] &;

A  //  MatrixForm

Restore the default print routine:

$PrePrint  =   .

{x1Actual,  x2Actual}  =  
      A[[All,  3]]

Thus the actual solution of the linear system is:
	(Subscript[x, 1])  =  
        1.001001001  …, 
(Subscript[x, 2])  =  
        0.998998998  … .

Gauss-Jordan algorithm with 3-digit precision

To see the difficulty with the straightforward, unadorned Gauss-Jordan algorithm, simulate a computer that can only do arithmetic on numbers with precision 3 decimal digits.

To force a floating-point number to have a precision of exactly 3 decimal digits, use the form x`3.  Thus 2.`3 means the number with value 2 but having precision only 3 decimal digits.  And 1.`3*^-3 means the floating-point number having mantissa 1 with precision 3 and base-10 exponent -3, that is, the number 1.00 × (10^(-3)) where the mantissa 1.00 has precision of 3 decimal digits (the exponent -3 is an exact integer).

The elementary row operations are the same ones as before, and in the same order:

A = (({{1.`3*^-3, 1.`3, 1.`3}, {1.`3, 1.`3, 2.`3}})); A  //  MatrixForm

For later reference, keep the value of the matrix after the first scale operation.

A  =  A1  =  
        scale[A,  1/A[[1, 1]],  
          1];  A  //  MatrixForm

Here comes the crucial step:

A  =  
      addrow[A,  -A[[2, 1]]/
          A[[1, 1]],  1,  2];  
    A  //  MatrixForm

The result just obtained is very different from what was obtained with exact rational arithmetic.  Why?

?addrow

The entry A[[2, 1]] after the addrow operation is the result of the following calculation:

c  =  -A1[[2, 1]]/
          A1[[1, 
            1]];
x  =  
        A1[[1,  
          2]];
c x
y  =  A1[[2,  2]]
y  +  c x

Now continue with the Gauss-Jordan algorithm.  The roundoff error just introduced is propagated in the remaining steps:

A  =  
      scale[A,  1/A[[2, 2]],  2];  
    A   //  MatrixForm

A  =  
      addrow[A,  -A[[1, 2]]/
          A[[2, 2]],  2,  1];  
    A  //  MatrixForm

{x1GJ,  x2GJ}  =  
      A[[All,  3]]

Thus the solution of the linear system obtained by straightforward Gauss-Jordan algorithm (in precision 3) is:
	(Subscript[x, 1])  =  0. , 
(Subscript[x, 2])  =  
        1.

Whereas the value of (Subscript[x, 2]) obtained here is as close to the actual value as can possibly be obtained with only 3-digit precision, the value of (Subscript[x, 1]) has 100% error!

Gauss-Jordan algorithm with partial pivoting in 3-digit precision

To try to eliminate the roundoff error in the straightforward Gauss-Jordan algorithm, use the strategy of partial pivoting:

Before scaling the leading nonzero entry in the current row, swap this row with the row at or below its position that has the entry of greatest absolute value in the column of the leading nonzero entry.

Apply this strategy, again with the numbers having a precision of 3 digits.

A = (({{1.`3*^-3, 1.`3, 1.`3}, {1.`3, 1.`3, 2.`3}})); A  //  MatrixForm

The leading nonzero entry in row 1 is 0.00100, but the entry 1.00 in column 1 below it has larger absolute value.  So according to the partial pivoting strategy, the first step is to swap the two rows:

A  =  swap[A,  1,  2];  A  //  MatrixForm

Now proceed as before.  (The scaling operation is perhaps unnecessary, but proceed blindly.)

A  =  
      scale[A,  1/A[[1, 1]],  1];  
    A  //  MatrixForm

Here comes the step that, without partial pivoting, introduced roundoff error:

A  =  
      addrow[A,  -A[[2, 1]]/
          A[[1, 1]],  1,  2];  
    A  //  MatrixForm

A  =  
      scale[A,  1/A[[2, 2]],  2];  
    A   //  MatrixForm

A  =  
      addrow[A,  -A[[1, 2]]/
          A[[2, 2]],  2,  1];  
    A  //  MatrixForm

{x1pp,  x2pp}  =  
      A[[All,  3]]

Thus the solution of the linear system obtained by the Gauss-Jordan algorithm with partial pivoting (in precision 3) is:
	(Subscript[x, 1])  =  1. , 
(Subscript[x, 2])  =  
        1.

The values of both (Subscript[x, 1]) and (Subscript[x, 2]) are as close to the actual values as can possibly be obtained with only 3-digit precision.

Partial pivoting may not be enough!

Change the linear system to the equivalent one obtained by multiplying the first equation by 1000:
	(Subscript[x, 1])  +  1000  (Subscript[x, 2])  =  
        1000, 
               (Subscript[x, 1]) +  
          (Subscript[x, 2])  =  2. 
The augmented matrix of this system is:
	A = (({{1, 1000, 1000}, {1, 1, 2}})).

Use again the Gauss-Jordan algorithm with partial pivoting, again with the numbers having a precision of 3 digits.

A = (({{1.`3, 1.`3*^3, 1.`3*^3}, {1.`3, 1.`3, 2.`3}})); A  //  MatrixForm

The entries in column 1 are the same, so according to the partial pivoting strategy, a swap is not done now.

A  =  
      scale[A,  1/A[[1, 1]],  1];  
    A  //  MatrixForm

Here comes the crucial step:

A  =  
      addrow[A,  -A[[2, 1]]/
          A[[1, 1]],  1,  2];  
    A  //  MatrixForm

And now we are back with the same roundoff error as was introduced for the original system without partial pivoting.  And once again this roundoff error propagates through the remaining steps to the final answer.

A  =  
      scale[A,  1/A[[2, 2]],  2];  
    A   //  MatrixForm

A  =  
      addrow[A,  -A[[1, 2]]/
          A[[2, 2]],  2,  1];  
    A  //  MatrixForm

{x1Newpp,  x2Newpp}  =  
      A[[All,  3]]

Thus the solution of the linear system obtained (in precision 3) by the Gauss-Jordan algorithm with partial pivoting—where, in fact, no actual partial pivoting was required—is:
	(Subscript[x, 1])  =  0. , 
(Subscript[x, 2])  =  
        1. 
Once again there is 100% error in (Subscript[x, 1]).

The moral is that it is not so much the sizes of the entries of the matrix themselves, but rather the size of their ratios to one another.  The key is to avoid (relatively speaking) "large" intermediate results.
A Gaussian elimination implementation that is to minimize the propagation of roundoff error needs to use a more sophisticated strategy than simple partial pivoting.  To learn about such sophisticated strategies, take a course in numerical analysis, or consult a book about this subject.

An example using machine precision

Perhaps you think the example above is too artificial in that it uses a simulated machine precision of only 3 decimal digits.  After all, real computers use 15 or 16 digits precision for ordinary floating-point calculations.

It's easy to modify the example to demonstrate the very same need for a strategy to minimize roundoff error.  Change it to:

(10^(-17)) (Subscript[x, 1])  +  (Subscript[x, 2])  =  
        1, 
                 (Subscript[x, 1]) +  
          (Subscript[x, 2])  =  2. 
The augmented matrix of this system is:
	A = (({{(10^(-17)), 1, 1}, {1, 1, 2}}))..

Gauss-Jordan algorithm with exact arithmetic

A = (({{(1/(10^7)), 1, 1}, {1, 1, 2}}));

A  =  
      scale[A,  1/A[[1, 1]],  1];  
    A  //  MatrixForm

A  =  
      addrow[A,  -A[[2, 1]]/
          A[[1, 1]],  1,  2];  
    A  //  MatrixForm

A  =  
      scale[A,  1/A[[2, 2]],  2];  
    A   //  MatrixForm

A  =  
      addrow[A,  -A[[1, 2]]/
          A[[2, 2]],  2,  1];  
    A  //  MatrixForm

A  =  N[A];  A  //  MatrixForm

$PrePrint  =  NumberForm[#,  16] &;
A  //  MatrixForm
$PrePrint  =   .

Gauss-Jordan algorithm with machine precision example—no partial pivoting

A = (({{1.*^-17, 1., 1.}, {1., 1., 2.}}))

A  =  
      scale[A,  1/A[[1, 1]],  1];  
    A  //  MatrixForm

A  =  
      addrow[A,  -A[[2, 1]]/
          A[[1, 1]],  1,  2];  
    A  //  MatrixForm

A  =  
      scale[A,  1/A[[2, 2]],  2];  
    A   //  MatrixForm

A  =  
      addrow[A,  -A[[1, 2]]/
          A[[2, 2]],  2,  1];  
    A  //  MatrixForm

Once again, despite using full machine precision, the error in (Subscript[x, 1]) is 100%.

Gauss-Jordan algorithm with machine precision—with partial pivoting

A = (({{1.*^-17, 1., 1.}, {1., 1., 2.}}));

A  =  swap[A,  1,  2];  A  //  MatrixForm

A  =  
      scale[A,  1/A[[1, 1]],  1];  
    A  //  MatrixForm

A  =  
      addrow[A,  -A[[2, 1]]/
          A[[1, 1]],  1,  2];  
    A  //  MatrixForm

A  =  
      scale[A,  1/A[[2, 2]],  2];  
    A   //  MatrixForm

A  =  
      addrow[A,  -A[[1, 2]]/
          A[[2, 2]],  2,  1];  
    A  //  MatrixForm

Through the use of partial pivoting, the roundoff error has been treated.