Consider solving the following augmented system
[22c2c112].Assume we are working in 5 digit aritmetic and c≥106. We perform GE with partial pivoting. Initially, we do no row swaps. After one step, R2−1/2R1→R2 we have
[22c2c0fl(1−c)fl(2−c)].In 5-digit arithmetic fl(1−c)=−c and fl(2−c)=−c.
Performing backward subsitution, we get x2=−c/(−c)=1 x1=(2c−2cx2)/2=0
But
[22c11][01]=[2c1]≠[2c2]If we instead interchange rows because of the large elements in the first row, we'd perform
[11222c2c].R2−2R1→R2[1120fl(2c−2)fl(2c−4)].In 5-digit arithmetic fl(2c−2)=2c and fl(2c−4)=2c.
In 5-digit arithmetic
[22c11][11]=[fl(2+2c)2]=[2c2]Let A be an n×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=[a11a12⋯a1n0a22⋯a2n⋮⋱⋮00⋯ak+1,k+1ak+1,k+2⋯ak+1,n00⋯ak+2,k+1ak+2,k+1ak+2,n⋮⋮⋮⋮⋮00⋯an,kan,k+1an,n]With partial pivoting, we swap rows so as to replace ak+1,k+1 with the largest element (in absolute value) among the set {ak+1,k+1,ak+2,k+1,…,an,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
si=max1≤j≤n|ai,j|.The number si 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×(n+1) matrix.
Compute α=max1≤i≤n|Ai,1|/si.
Then choose the smallest value of p such that |Ap,1|/sp=α.
Perform R1↔Rp and s1↔sp
Then row reduce, below A11.
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'))'
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))
Complete pivoting involves both row and column interchanges. Let A be an n×n matrix. After k steps of Gaussian elimination, A has zeros below the first k diagonal entries:
A=[a11a12⋯a1n0a22⋯a2n⋮⋱⋮00⋯ak+1,k+1ak+1,k+2⋯ak+1,n00⋯ak+2,k+1ak+2,k+1ak+2,n⋮⋮⋮⋮⋮00⋯an,kan,k+1an,n]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.
Perform Gaussian elimination with partial pivoting to solve Ax=b where
A=[1−112−21030],b=[111].In 5-digit arithmetic with c>106
[22c2c112]unknown = [x1x2]C1↔C2In 5-digit arithmetic, this gives the solution of the system
[22c11][11]=[fl(2+2c)2]=[2c2]In 5-digit arithmetic with c>106
[c−111121]unknown = [x1x2]R1↔R2C1↔C2This is the solution (in 5-digit arithmetic)
[c−1112][−11]=[fl(1−c−1)1]=[11]