Pivoting Strategies

Full (standard) Gaussian elimination

function A = GE(A)

INPUT: A is an n x m matrix

OUTPUT:  A an n x m upper-triangular matrix, or a message of failure

STEP 1: For i = 1,2,...,n-1 do STEPS 2-5
  STEP 2: Let p >= i be the smallest integer such that A(p,i) ~= 0. 
  STEP 3: If p cannot be found then
    DISPLAY('Method failed: matrix is rank deficient')
    OUTPUT(A);
    STOP.
  STEP 4: If p > i do Ri <-> Rp on A
  STEP 5: For j = i+1,i+2,...,n do STEP 6
    STEP 6: Do R_j - A(j,i)/A(i,i) R_i --> R_j on A
STEP 7: If A(n,n) = 0
  DISPLAY('Method failed: matrix is rank deficient')
  OUTPUT(A)
STEP 8: OUTPUT(A); STOP.

Backward substitution

function x = Backsub(U)

INPUT: U is an n x n + 1 upper-triangular matrix with non-zero diagonal entries

OUTPUT: the solution to U(1:n,1:n)x=U(1:n,n+1)

STEP 1: Set x = U(:,n+1)
STEP 2: 
    If U(n,n) = 0 then
        OUTPUT('Method failed: singular matrix')
        STOP.
    Else set x(n) = U(n,n+1)/U(n,n)
STEP 3: For i = n-1,...,1 do STEP 4
    STEP 4:
      If U(i,i) = 0 then
        OUTPUT('Method failed: singular matrix')
        STOP.
      Else set
        x(i) = [U(i,n+1) - SUM( U(i,j)x(j), j= i+1,...,n )]
        x(i) = x(i)/U(i,i)
STEP 5: OUTPUT(x); STOP.

A = rand(5); b = rand(5,1);
GE([A,b])

ans =

    0.8147    0.0975    0.1576    0.1419    0.6557    0.7577
         0    0.1701    0.7954    0.2640   -0.6933   -0.0993
         0         0   -1.5541    0.0682    2.9146    0.5846
         0    0.0000    0.0000   -0.8441   -3.2039   -1.0748
         0   -0.0000   -0.0000         0   -0.1376   -0.3364

A = rand(3); b = rand(3,1);
U = GE([A,b]);
x = Backsub(U)

x =

   -1.7588
   -1.1531
    1.9134

A*x-b

ans =

   1.0e-15 *

         0
   -0.1110
   -0.8882

One can encounter problems when the diagonal entries of $U$ become very small.

A = rand(100); b = rand(100,1);
A(1,:) = A(2,:) + 0.0001*A(1,:);
U = GE([A,b]);
min(abs(diag(U)))  % the smallest value of the diagonal of U
x = Backsub(U);
max(abs(A*x-b))

ans =

   1.4312e-05


ans =

   1.1014e-09

An example

Apply Gaussian elimination in 5-digit arithmetic to solve $Ax=b$ where

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

Assume $c > 10^6$ and $c = fl(c)$.

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

It follows that $fl(2-c) = fl(1-c) = -c$ because $c > 10^6$

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

Backward subsitution gives

$$x_2 = 1, \quad x_1 = c(1-x_2) = 0$$

$$Ax = b, \quad A = \begin{bmatrix} c^{-1} & 1 \\ 1 & 2 \end{bmatrix}, \quad b = \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \quad x = \begin{bmatrix} 0 \\ 1 \end{bmatrix}.$$ $$\begin{bmatrix} c^{-1} & 1 \\ 1 & 2 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \end{bmatrix} \neq b$$

If we swap the rows first $R_1 \leftrightarrow R_2$:

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

Because we are using 5-digit arithmetic and $c > 10^6$, $fl(1-2c^{-1}) = 1$ so we have:

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

Backward subsitution gives $x_2 = 1$ and $x_1 = 1 - 2x_2 = -1$

$$Ax = b, \quad A = \begin{bmatrix} c^{-1} & 1 \\ 1 & 2 \end{bmatrix}, \quad b = \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \quad x = \begin{bmatrix} 0 \\ 1 \end{bmatrix}.$$ $$\begin{bmatrix} c^{-1} & 1 \\ 1 & 2 \end{bmatrix} \begin{bmatrix} -1 \\ 1 \end{bmatrix} = \begin{bmatrix} fl(1-c^{-1}) \\ 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix} =b$$

Partial pivoting

To help remedy this issue, we need to choose a new way of deciding when to interchange rows. Rather than only interchanging if $A_{ii} = 0$, we can interchange rows to increase the size of the elements on the diagonal. This is called partial pivoting.

One step of 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}|$.

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

Perform $R_1 \leftrightarrow R_p$.

Do one step of standard Gaussian elimination.

To find the index of element that is the largest in absolute value:

n = 7; a = rand(1,n); MAX = 0; iMax = 0;
for i = 1:n
    if abs(a(i)) > MAX
        MAX = abs(a(i));
        iMax = i;
    end
end
a
iMax

a =

    0.1549    0.5555    0.7905    0.4439    0.9958    0.4366    0.3044


iMax =

     5

function iMAX = argmax(a) % returns zero if all entries in a are zero
    MAX = 0;
    iMAX = 0;
    for i = 1:length(a)
        if abs(a(i)) > MAX
            MAX = abs(a(i));
            iMAX = i;
        end
    end
end

Gaussian elimination with partial pivoting

function A = GEpp(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: For i = 1,2,...,n-1 do STEPS 2-6
  STEP 2: Set p = argmax(A(i:n,i))
  STEP 3: If p = 0 then
    DISPLAY('Method failed: matrix is rank deficient')
    OUTPUT(A);
    STOP.
  STEP 4: Set p = p + i - 1  %Adjust p to correspond to the right row of A
  STEP 5: Do Ri <-> Rp on A
  STEP 6: For j = i+1,i+2,...,n do STEP 5
    STEP 7: Do R_j - A(j,i)/A(i,i) R_i --> R_j on A
STEP 8: If A(n,n) = 0
  DISPLAY('Method failed: matrix is rank deficient')
  OUTPUT(A)
STEP 9: OUTPUT(A); STOP.

argmax([1,2])
argmax([2,3,5])
argmax([0,0])

ans =

     2


ans =

     3


ans =

     0

n = 200; A = rand(n); b = rand(n,1);
A(1,:) = A(2,:) + 0.0000000001*A(1,:);
U = GE([A,b]);
%min(abs(diag(U)))  % the smallest value of the diagonal of U
x = Backsub(U);
max(abs(A*x-b))

U = GEpp([A,b]);
%min(abs(diag(U)))  % the smallest value of the diagonal of U
x = Backsub(U);
max(abs(A*x-b))

ans =

    0.0048


ans =

   2.7072e-05

Example

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

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