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.
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])
A = rand(3); b = rand(3,1);
U = GE([A,b]);
x = Backsub(U)
A*x-b
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))
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)$.
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$$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$
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
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
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])
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))
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}.$$