381 - A19 - Matrix Factorization

145 days ago by Professor381

John Travis

Mississippi College

When solving $Ax=b$ for a general system of equations, one often might want to use Gaussian Elimination with scaled partial pivoting.  However, one might want to solve such a system several times with a constant matrix $A$ but with a variety of "forcing term" vectors for b.  Using the standard approach would replicate a significant number of calculations in the reduction of A to triangular form.  There is however a way to save the intermediate steps of Gaussian Elimination so that only a final forward and backward solve would be necessary and not all of the elimination steps.

In general, one might be able to factor $A = LU$ where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix.  Then to solve $Ax=b$ simply replace and solve

$LUx = b$

using a two-step approach:

Solve $Ly=b$ for $y$ using forward solving.

Solve $Ux=y$ for $x$ using backward solving.

Since each of these require only $O(n^2)$ calculations, then repeated solves with various values for $b$ offer significant computational savings.
# Matrix Factorization using Gaussian Elimination clear A = [2 -1 3 0; -1 1 -2 -3; 0 0 4 1; 5 0 0 4] n = size(A)(1); L = zeros(n,n); for step = 1:n pivot = A(step,step); L(step:n,step) = A(step:n,step); A(step,step:n) = A(step,step:n)/pivot; for k=step+1:n multiplier = -A(k,step); A(k,step:n) = A(k,step:n) + multiplier*A(step,step:n); # kth row, all columns end end disp('The Factorization: ') L U=A disp('and multiplying these together gives the original A: ') A = L*U 
        
A =

 2 -1 3 0
 -1 1 -2 -3
 0 0 4 1
 5 0 0 4

The Factorization: 
L =

 2 0 0 0
 -1 0.5 0 0
 0 0 4 0
 5 2.5 -5 20.25

U =

 1 -0.5 1.5 0
 0 1 -1 -6
 0 0 1 0.25
 0 0 0 1

and multiplying these together gives the original A: 
A =

 2 -1 3 0
 -1 1 -2 -3
 0 0 4 1
 5 0 0 4
A =

 2 -1 3 0
 -1 1 -2 -3
 0 0 4 1
 5 0 0 4

The Factorization: 
L =

 2 0 0 0
 -1 0.5 0 0
 0 0 4 0
 5 2.5 -5 20.25

U =

 1 -0.5 1.5 0
 0 1 -1 -6
 0 0 1 0.25
 0 0 0 1

and multiplying these together gives the original A: 
A =

 2 -1 3 0
 -1 1 -2 -3
 0 0 4 1
 5 0 0 4
# Now, for a given vector b, solve Ax=b using LUx = b by forward and backward substitution clear x # clear y # if this cell is to be evaluated again then uncomment b = [8; -2; 7; 4] # First, forward substitution for Ly = b y(1) = b(1)/L(1,1); # the first equation for k = 2:n # each successive equation temp = 0; for j = 1:k-1 # accumulate the previously determined terms for row k temp += L(k,j)*y(j); end y(k) = (b(k)-temp)/L(k,k); end disp('Note the temporary values in the vector y:') y = y' # Finally, backward substitution for Ux = y x(n) = y(n)/U(n,n); # the last equation for k = n:-1:1 # each successive equation temp = 0; for j = n:-1:k+1 # accululate the previously determined terms for row k temp += U(k,j)*x(j); end x(k) = (y(k)-temp)/U(k,k); end disp('With the final answer for Ax=b as x equals:') x = x' disp('and to check, look at b-Ax:') b - A*x 
        
b =

 8
 -2
 7
 4

Note the temporary values in the vector y:
y =

 4
 4
 1.75
 -0.851852

With the final answer for Ax=b as x equals:
x =

 1.48148
 0.851852
 1.96296
 -0.851852

and to check, look at b-Ax:
ans =

 0
 0
 0
 8.88178e-16
b =

 8
 -2
 7
 4

Note the temporary values in the vector y:
y =

 4
 4
 1.75
 -0.851852

With the final answer for Ax=b as x equals:
x =

 1.48148
 0.851852
 1.96296
 -0.851852

and to check, look at b-Ax:
ans =

 0
 0
 0
 8.88178e-16
# Now, let's put this all together in functions function ans = forwardsolve(mat,rhs) # presumes mat is lower n = size(rhs)(1); y(1) = rhs(1)/mat(1,1); # the first equation for k = 2:n # each successive equation temp = 0; for j = 1:k-1 # accumulate the previously determined terms for row k temp += mat(k,j)*y(j); end y(k) = (rhs(k)-temp)/mat(k,k); end ans = y'; endfunction function ans = backwardsolve(mat,rhs) # presumes mat is upper n = size(rhs)(1); x(n) = rhs(n)/mat(n,n); # the last equation for k = n:-1:1 # each successive equation temp = 0; for j = n:-1:k+1 # accululate the previously determined terms for row k temp += mat(k,j)*x(j); end x(k) = (rhs(k)-temp)/mat(k,k); end ans = x'; endfunction 
       
disp('Solving using') b = [8; -2; 7; 4] y = forwardsolve(L,b) disp('yields') x = backwardsolve(U,y) disp('with b-Ax = ') b-A*x disp('Solving using') b = [-9; -2; -7; 5] y = forwardsolve(L,b) disp('yields') x = backwardsolve(U,y) disp('with b-Ax = ') b-A*x 
        
Solving using
b =

 8
 -2
 7
 4

y =

 4
 4
 1.75
 -0.851852

yields
x =

 1.48148
 0.851852
 1.96296
 -0.851852

with b-Ax = 
ans =

 0
 0
 0
 8.88178e-16

Solving using
b =

 -9
 -2
 -7
 5

y =

 -4.5
 -13
 -1.75
 2.53086

yields
x =

 -1.02469
 -0.197531
 -2.38272
 2.53086

with b-Ax = 
ans =

 0
 0
 0
 -7.10543e-15
Solving using
b =

 8
 -2
 7
 4

y =

 4
 4
 1.75
 -0.851852

yields
x =

 1.48148
 0.851852
 1.96296
 -0.851852

with b-Ax = 
ans =

 0
 0
 0
 8.88178e-16

Solving using
b =

 -9
 -2
 -7
 5

y =

 -4.5
 -13
 -1.75
 2.53086

yields
x =

 -1.02469
 -0.197531
 -2.38272
 2.53086

with b-Ax = 
ans =

 0
 0
 0
 -7.10543e-15
# We can even find the inverse matrix column by column by solving Ax=e_k, where e_k is the kth unit vector. e1 = [1;0;0;0]; e2 = [0;1;0;0]; e3 = [0;0;1;0]; e4 = [0;0;0;1]; inv1 = backwardsolve(U,forwardsolve(L,e1)); inv2 = backwardsolve(U,forwardsolve(L,e2)); inv3 = backwardsolve(U,forwardsolve(L,e3)); inv4 = backwardsolve(U,forwardsolve(L,e4)); invA = [inv1 inv2 inv3 inv4] disp('Check our result by comparing to the built-in function for inverses:') format short invA-inv(A) 
        
invA =

 0.197531 0.197531 -0.0493827 0.160494
 -0.419753 0.580247 0.604938 0.283951
 0.0617284 0.0617284 0.234568 -0.0123457
 -0.246914 -0.246914 0.0617284 0.0493827

Check our result by comparing to the built-in function for inverses:
ans =

   1.1102e-16   8.3267e-17   0.0000e+00   0.0000e+00
   1.1102e-16   1.1102e-16   0.0000e+00  -5.5511e-17
   6.9389e-18   6.9389e-18  -2.7756e-17   0.0000e+00
  -2.7756e-17  -2.7756e-17   6.9389e-18   0.0000e+00
invA =

 0.197531 0.197531 -0.0493827 0.160494
 -0.419753 0.580247 0.604938 0.283951
 0.0617284 0.0617284 0.234568 -0.0123457
 -0.246914 -0.246914 0.0617284 0.0493827

Check our result by comparing to the built-in function for inverses:
ans =

   1.1102e-16   8.3267e-17   0.0000e+00   0.0000e+00
   1.1102e-16   1.1102e-16   0.0000e+00  -5.5511e-17
   6.9389e-18   6.9389e-18  -2.7756e-17   0.0000e+00
  -2.7756e-17  -2.7756e-17   6.9389e-18   0.0000e+00