7 Gaussian Elimination and LU Factorization - IIT

1 7 Gaussian Elimination and LU FactorizationIn this final section on matrix Factorization methods for solvingAx=bwe want totake a closer look at Gaussian Elimination (probably the best known method for solvingsystems of linear equations).The basic idea is to use left-multiplication ofA Cm mby (elementary) lowertriangular matrices ,L1, L2, .. , Lm 1to convertAto upper triangular form, ,Lm 1Lm 2.. L2L1 =eLA= that the product of lower triangular matrices is a lower triangular matrix, andthe inverse of a lower triangular matrix is also lower triangular. Therefore, LA=U A=LU,whereL= L 1. This approach can be viewed astriangular Why Would We Want to Do This?

2 Consider the systemAx=bwith LU factorizationA=LU. Then we haveL Ux =y= we can perform (a now familiar) 2-step solution the lower triangular systemLy=bforyby forward the upper triangular systemUx=yforxby back , consider the problemAX=B( , many different right-hand sides thatare associated with the same system matrix). In this case we need to compute thefactorizationA=LUonly once, and thenAX=B LU X=B,and we proceed as many forward substitutions (in parallel). X=Yby many back substitutions (in parallel).In order to appreciate the usefulness of this approach note that the operations countfor the matrix Factorization isO(23m3), while that for forward and back substitution isO(m2).

3 ExampleTake the matrixA= 1 1 12 3 54 6 8 59and compute its LU Factorization by applying elementary lower triangular transforma-tion chooseL1such that left-multiplication corresponds to subtracting multiples ofrow 1 from the rows below such that the entries in the first column ofAare zeroed out(cf. the first homework assignment). ThusL1A= 1 0 0 2 1 0 4 0 1 1 1 12 3 54 6 8 = 1 1 10 1 30 2 4 .Next, we repeat this operation analogously forL2(in order to zero what is left incolumn 2 of the matrix on the right-hand side above):L2(L1A) = 1 0 00 1 00 2 1 1 1 10 1 30 2 4 = 1 1 10 1 30 0 2 = (L2L1) 1=L 11L 12withL 11= 1 0 02 1 04 0 1 andL 12= 1 0 00 1 00 2 1 ,so thatL= 1 0 02 1 04 2 1.

4 RemarkNote thatLalways is aunit lower triangularmatrix, , it has ones on thediagonal. Moreover,Lis always obtained as above, , the multipliers are accumulatedinto the lower triangular part with a change of claims made above can be verified as follows. First, we note that the multipliersinLkare of the form`jk=ajkakk,j=k+ 1, .. , m,so thatLk= `k+1, `m,k1 .Now, let`k= `k+1, `m,k .60 ThenLk=I `ke k, and therefore(I `ke k) =Lk(I+`ke k) L 1k=I `ke k`ke k=I,since the inner producte k`k= 0 because the only nonzero entry inek(the 1 in thek-thposition) does not hit any nonzero entries in`kwhich start in thek+ 1-st , for anykwe haveL 1k= `k+1, `m,k1 as addition,L 1kL 1k+1= (I+`ke k)(I+`k+1e k+1)=I+`ke k+`k+1e k+1+`ke k`k+1 =0ek+1= `k+1, `k+2,k+ `m,k`m,k+11 ,and in general we haveL=L 11.

5 L 1m 1= 1`2, `3, `m,1`m,2.. `m,m 11 .We can summarize the Factorization inAlgorithm(LU Factorization )InitializeU=A,L=Ifork= 1 :m 1forj=k+ 1 :m61L(j, k) =U(j, k)/U(k, k)U(j, k:m) =U(j, k:m) L(j, k)U(k, k:m) practice one can actually store bothLandUin the original matrixAsince it is known that the diagonal ofLconsists of all LU Factorization is the cheapest Factorization algorithm. Its operations countcan be verified to beO(23m3).However,LU Factorization cannot be guaranteed to be stable. The following exam-ples illustrate this fundamental problem is given if we encounter azero pivotas inA= 1 1 12 2 54 6 8 = L1A= 1 1 10 0 30 2 4.

6 Now the (2,2) position contains a zero and the algorithm will break down since it willattempt to divide by more subtle example is the followingbackward instability. TakeA= 1 1 12 2 + 54 6 8 with small . If = 1 then we have the initial example in this chapter, and for = 0we get the previous Factorization will result inL1A= 1 1 10 30 2 4 andL2L1A= 1 1 10 30 0 4 6 = multipliers wereL= 1 0 02 1 042 1 .Now we assume that a right-hand sidebis given asb= 100 and we attempt to solveAx= is on the order of machine accuracy, then the 4 in the entry 4 6 inUisinsignificant. Therefore, we have U= 1 1 10 30 0 6 and L=L,which leads to L U= 1 1 12 2 + 54 6 4 6= fact, the product is significantly different fromA.

7 Thus, using Land Uwe arenot able to solve a nearby problem , and thus the LU Factorization method isnotbackward we use the Factorization based on Land Uwith the above right-hand sideb, thenwe obtain x= 112 23 223 23 112 2 23 .Whereas if we were to use the exact factorizationA=LU, then we get the exactanswerx= 4 72 322 3 2 12 3 73 23 23 .RemarkEven though Land Uare close toLandU, the product L Uis not close toLU=Aand the computed solution xis PivotingExampleThe breakdown of the algorithm in our earlier example withL1A= 1 1 10 0 30 2 3 can be prevented by simplyswapping rows, , instead of trying to applyL2toL1 Awe first createP L1A= 1 0 00 0 10 1 0 L1A= 1 1 10 2 30 0 3 and are generally, stability problems can be avoided by swapping rows before applyingLk, , we performLm 1Pm 1.

8 L2P2L1P1A= strategy we use for swapping rows in stepkis to find the largest element in columnkbelow (and including) the diagonal the so-calledpivot element and swap its rowwith rowk. This process is referred to aspartial(row)pivoting. Partial column pivotingand complete (row and column) pivoting are also possible, but not very again the matrixA= 1 1 12 2 + 54 6 8 The largest element in the first column is the 4 in the (3,1) position. This is our firstpivot, and we swap rows 1 and 3. ThereforeP1A= 0 0 10 1 01 0 0 1 1 12 2 + 54 6 8 = 4 6 82 2 + 51 1 1 ,and thenL1P1A= 1 0 0 121 0 140 1 4 6 82 2 + 51 1 1 = 4 680 1 10 12 1.

9 Now we need to pick the second pivot element. For sufficiently small (in fact, unless12< <32), we pick 1 as the largest element in the second column below the firstrow. Therefore, the second permutation matrix is just the identity, and we haveP2L1P1A= 1 0 00 1 00 0 1 4 680 1 10 12 1 = 4 680 1 10 12 1 .To complete the Elimination phase, we need to perform the Elimination in the secondcolumn:L2P2L1P1A= 1 0 00 1 0012( 1)1 4 680 1 10 12 1 = 4 680 1 10 03 2 2( 1) = lower triangular matrixLis given byL=L 11L 12= 100121014 12( 1)1 ,and assuming that 1 1 we get L= 1 0 0121 014121 and U= 4 6 80 1 10 0 32.

10 64If we now check the computed Factorization L U, then we see L U= 4 6 82 2 51 1 1 =P A,which is just a permuted version of the original matrixAwith permutation matrixP=P2P1= 0 0 10 1 01 0 0 .Thus, this approach was backward , since we have the factorizationP A=LU, we can solve the linear systemAx=basP Ax=Pb LUx=Pb,and apply the usual two-step the lower triangular systemLy= the upper triangular systemUx= yieldsx= 7+4 3+2 22 3 2 12 3 73 23 23 .If we use the rounded factors Land Uinstead, then the computed solution is x= 73 23 23 ,which is the exact answer to the problem (see also the Maple ).