# 381 - A19 - Matrix Factorization

## 283 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