Consider solving the following augmented system

$$\left[\begin{array}{cc|c} 2 & 2c & 2c \\ 1 & 1 & 2 \end{array} \right].$$

Assume we are working in $5$ digit aritmetic and $c \geq 10^6$. We perform GE with partial pivoting. Initially, we do no row swaps. After one step, $R_2- 1/2 R_1 \to R_2$ we have

$$\left[\begin{array}{cc|c} 2 & 2c & 2c \\ 0 & fl(1-c) & fl(2-c) \end{array} \right].$$

In 5-digit arithmetic $fl(1- c) = -c$ and $fl(2-c) = -c$.

$$\left[\begin{array}{cc|c} 2 & 2c & 2c \\ 0 & -c & -c \end{array} \right].$$

Performing backward subsitution, we get $$x_2 = -c/(-c) = 1$$ $$x_1 = (2c - 2c x_2)/2 = 0$$

But

$$\left[\begin{array}{cc|c} 2 & 2c \\ 1 & 1 \end{array} \right] \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 2c \\ 1 \end{bmatrix} \neq \begin{bmatrix} 2c \\ 2 \end{bmatrix}$$

If we instead interchange rows because of the large elements in the first row, we'd perform

$$\left[\begin{array}{cc|c} 1 & 1 & 2 \\ 2 & 2c & 2c \end{array} \right].$$ $$R_2 - 2 R_1 \to R_2$$ $$\left[\begin{array}{cc|c} 1 & 1 & 2 \\ 0 & fl(2c-2) & fl(2c -4) \end{array} \right].$$

In 5-digit arithmetic $fl(2c- 2) = 2c$ and $fl(2c-4) = 2c$.

$$\left[\begin{array}{cc|c} 1 & 1 & 2 \\ 0 & 2c & 2c \end{array} \right]$$ $$x_2 = 2c/2c = 1$$ $$x_1 = 2 - x_2 = 1$$

In 5-digit arithmetic

$$\left[\begin{array}{cc|c} 2 & 2c \\ 1 & 1 \end{array} \right] \begin{bmatrix} 1 \\ 1 \end{bmatrix} = \begin{bmatrix} fl(2+2c) \\ 2 \end{bmatrix} = \begin{bmatrix} 2c \\ 2 \end{bmatrix}$$

Scaled partial pivoting

Let $A$ be an $n \times n$ matrix. After $k$ steps of Gaussian elimination, $A$ has zeros below the first $k$ diagonal entries. With partial pivoting, we look at the entries that are in the $(k+1)$th column and rows $k+1$ through $n$:

$$A = \begin{bmatrix} a_{11} & a_{12} &&& \cdots &&& a_{1n}\\ 0 & a_{22} &&&\cdots&&& a_{2n}\\ \vdots &\ddots &&&&&& \vdots \\ 0 &0 & \cdots & a_{k+1,k+1} & a_{k+1,k+2} &\cdots && a_{k+1,n}\\ 0 &0 & \cdots & a_{k+2,k+1} & a_{k+2,k+1} &&& a_{k+2,n}\\ \vdots & \vdots && \vdots & \vdots &&& \vdots\\ 0 &0 & \cdots & a_{n,k} & a_{n,k+1} &&& a_{n,n}\\ \end{bmatrix}$$

With partial pivoting, we swap rows so as to replace $a_{k+1,k+1}$ with the largest element (in absolute value) among the set $\{a_{k+1,k+1}, a_{k+2,k+1},\ldots,a_{n,k+1}\}$.

Often, what is more important is the relative size of the first non-zero entry in the row to the largest entry (in absolute value). Before any row operations define

$$s_i = \max_{1 \leq j \leq n} |a_{i,j}|.$$

The number $s_i$ is computed once, is attached to row $i$, and will need to be swapped when row $i$ is swapped.

One step of scaled partial pivoting is as follows:

Let $A$ be an $n \times (n+1)$ matrix.

Compute $\alpha = \max_{1 \leq i \leq n} |A_{i,1}|/s_i$.

Then choose the smallest value of $p$ such that $|A_{p,1}|/s_p = \alpha$.

Perform $R_1 \leftrightarrow R_p$ and $s_1 \leftrightarrow s_p$

Then row reduce, below $A_{11}$.

Full Gaussian elimination with scaled partial pivoting

function A = GEspp(A)

INPUT: A is an n x m matrix

OUTPUT:  A an n x m upper-triangular matrix, or Inf if the method failed

STEP 1: Define s = max(abs(A(1:n,1:n)')' % the scale factors
STEP 2: For i = 1,2,...,n-1 do STEPS 3-8
  STEP 3: Set p = argmax(A(i:n,i)./s(i:n)) %This is a concise way to find the entry with the largest relative size
  STEP 4: If p = 0 then
    DISPLAY('Method failed: matrix is rank deficient')
    OUTPUT(A);
    STOP.
  STEP 5: Set p = p + i - 1  %Adjust p to correspond to the right row of A
  STEP 6: Do Ri <-> Rp on A  %if you swap rows
  STEP 7: Do si <-> sp on s  %you need to swap scale factors
  STEP 8: For j = i+1,i+2,...,n do STEP 9
    STEP 9: Do R_j - A(j,i)/A(i,i) R_i --> R_j on A
STEP 10: If A(n,n) = 0
  DISPLAY('Method failed: matrix is rank deficient')
  OUTPUT(A)
STEP 11: OUTPUT(A); STOP.

A = [5,2,3; 1,1,1; 4,3,1];
A
max(abs(A'))'

A =

     5     2     3
     1     1     1
     4     3     1


ans =

     5
     1
     4

Scaled partial pivoting is slightly better:

n = 100; A = rand(n); b = rand(n,1);
A(1,:) = (1:n).^4.*A(2,:) + A(1,:); % create a "difficult" matrix
U = GEpp([A,b]);
x = Backsub(U);
max(abs(x-A\b))

U = GEspp([A,b]);
x = Backsub(U);
max(abs(x-A\b))

ans =

   2.5175e-08


ans =

   1.0271e-08

Full (complete) pivoting

Complete pivoting involves both row and column interchanges. Let $A$ be an $n \times n$ matrix. After $k$ steps of Gaussian elimination, $A$ has zeros below the first $k$ diagonal entries:

$$A = \begin{bmatrix} a_{11} & a_{12} &&& \cdots &&& a_{1n}\\ 0 & a_{22} &&&\cdots&&& a_{2n}\\ \vdots &\ddots &&&&&& \vdots \\ 0 &0 & \cdots & a_{k+1,k+1} & a_{k+1,k+2} &\cdots && a_{k+1,n}\\ 0 &0 & \cdots & a_{k+2,k+1} & a_{k+2,k+1} &&& a_{k+2,n}\\ \vdots & \vdots && \vdots & \vdots &&& \vdots\\ 0 &0 & \cdots & a_{n,k} & a_{n,k+1} &&& a_{n,n}\\ \end{bmatrix}$$

We swap rows and columns to move the largest element (in absolute value) in the entire bottom right block to $(k+1,k+1)$ entry.

Example

Perform Gaussian elimination with partial pivoting to solve $A x = b$ where

$$A = \begin{bmatrix} 1 & -1 & 1 \\ 2 & -2 & 1 \\ 0 & 3 & 0 \end{bmatrix}, \quad b = \begin{bmatrix} 1 \\ 1\\ 1 \end{bmatrix}.$$

$$\left[ \begin{array}{ccc|c} 1 & -1 & 1 & 1 \\ 2 & -2 & 1 & 1 \\ 0 & 3 & 0 & 1 \end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_1\\x_2\\x_3 \end{bmatrix}$$

$$R_1 \leftrightarrow R_3$$

$$\left[ \begin{array}{ccc|c} 0 & 3 & 0 & 1 \\ 2 & -2 & 1 & 1\\ 1 & -1 & 1 & 1 \end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_1\\x_2\\x_3\end{bmatrix}$$

$$C_1 \leftrightarrow C_2$$

$$\left[ \begin{array}{ccc|c} 3 & 0 & 0 & 1 \\ -2 & 2 & 1 & 1\\ -1 & 1 & 1 & 1 \end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_2\\x_1\\x_3\end{bmatrix}$$

$$R_2 + 2/3 R_1 \to R_2$$ $$R_3 + 1/3 R_1 \to R_3$$

$$\left[ \begin{array}{ccc|c} 3 & 0 & 0 & 1 \\ 0 & 2 & 1 & 5/3\\ 0 & 1 & 1 & 4/3 \end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_2\\x_1\\x_3\end{bmatrix}$$

$$R_3 - 1/2 R_2 \to R_3$$

$$\left[ \begin{array}{ccc|c} 3 & 0 & 0 & 1 \\ 0 & 2 & 1 & 5/3\\ 0 & 0 & 1/2 & 1/2 \end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_2\\x_1\\x_3\end{bmatrix}$$

$$x_3 = 1, \quad x_1 = (5/3 - x_3)/2 = 1/3, \quad x_2 = 1/3$$

Example (where partial pivoting failed)

In 5-digit arithmetic with $c > 10^6$

$$\left[\begin{array}{cc|c} 2 & 2c & 2c \\ 1 & 1 & 2 \end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_1\\x_2\end{bmatrix}$$ $$C_1 \leftrightarrow C_2$$

$$\left[\begin{array}{cc|c} 2c & 2 & 2c \\ 1 & 1 & 2 \end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_2\\x_1\end{bmatrix}$$ $$R_2 - R_1/(2c) \to R_2$$

$$\left[\begin{array}{cc|c} 2c & 2 & 2c \\ 0 & fl(1-c^{-1}) & 1 \end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_2\\x_1\end{bmatrix}$$ $$fl(1-c^{-1}) = 1$$ $$\left[\begin{array}{cc|c} 2c & 2 & 2c \\ 0 & 1 & 1 \end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_2\\x_1\end{bmatrix}$$ $$x_ 1 =1, \quad x_2 = fl((2c - 2)/2c) = 1$$

In 5-digit arithmetic, this gives the solution of the system

$$\left[\begin{array}{cc|c} 2 & 2c\\ 1 & 1 \end{array} \right] \begin{bmatrix} 1 \\ 1 \end{bmatrix} = \begin{bmatrix} fl(2 + 2c) \\ 2 \end{bmatrix} = \begin{bmatrix} 2c \\ 2 \end{bmatrix}$$

Example (where standard Gaussian elimination failed)

In 5-digit arithmetic with $c > 10^6$

$$\left[\begin{array}{cc|c} c^{-1} & 1 & 1 \\ 1 & 2 & 1\end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_1\\x_2\end{bmatrix}$$ $$R_1 \leftrightarrow R_2 \quad C_1 \leftrightarrow C_2$$

$$\left[\begin{array}{cc|c} 2 & 1 & 1 \\ 1 & c^{-1} & 1\end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_2\\x_1\end{bmatrix}$$ $$R_2 - 1/2 R_1 \to R_2$$

$$\left[\begin{array}{cc|c} 2 & 1 & 1 \\ 0 & fl(c^{-1}-1/2) & 1/2\end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_2\\x_1\end{bmatrix}$$ $$fl(c^{-1} -1/2) = -1/2$$

$$\left[\begin{array}{cc|c} 2 & 1 & 1 \\ 0 & -1/2 & 1/2\end{array} \right] \quad \text{unknown = } \begin{bmatrix} x_2\\x_1\end{bmatrix}$$ $$x_1 = -1, \quad x_2 = fl( (1-x_1)/2) = 1$$

This is the solution (in 5-digit arithmetic)

$$\left[\begin{array}{cc|c} c^{-1} & 1 \\ 1 & 2 \end{array} \right] \begin{bmatrix} -1\\1\end{bmatrix} = \begin{bmatrix} fl(1-c^{-1}) \\ 1\end{bmatrix} = \begin{bmatrix} 1\\ 1\end{bmatrix}$$