LU Decomposition 

We are going to learn how to write A as the product of a lower triangular matrix times a lower triangular matrix.  In the process we will learn more about how elementary row operations work.  This first look at LU decomposition will hold when there are no row switches necessary in the Gausselimination process. 

We will now illustrate how to perform elementary row operations by apply the operations to the identity matrix and the pre-multiplying by the new matrix.  The act of pre-multiplying is multiplying on the left by the new matrix.   

>	with(linalg):

 

>	with(LinearAlgebra):

 

>	A:=<<2,1,-1>|<4,1,0>|<2,2,2>>;

 

>	I3:=IdentityMatrix(3);

 

>	E[1]:=addrow(I3,1,2,-1/2);

 

>	A1:=E[1].A;

 

Multiplication on the left by the elementary matrix E[1] is the same as the first step in the Gauss Elimination process.  The E[1] matrix was the same as the Gauss Elimination process  but carried out on the Identity matrix.  

>	E[2]:=addrow(I3,1,3,1/2);

 

>	A2:=E[2].A1;

 

We can do a whole row of operations by combining them just as in the way we do Gauss elimination . 

>	Ea[1]:=addrow(addrow(I3,1,2,-1/2),1,3,1/2);

 

>	A2a:=Ea[1].A;

 

We have obtained the same results as doing the operations one step at a time.   Pay close attention to what is being done in each step.  You are not going to be doing the problems by using the elementary matrices but they will help in understanting the process.  We not need to get rid of the 2 in the third row and second column.  Again we will use the elimentary matrices to accomplish this  

>	E[3]:=addrow(I3,2,3,2);

 

>	U:=E[3].A2a;

 

We can now look at what we really have that is  

= U  or = U.  It would appear that we were now half way to our LU decomposition.  The next thing we would like to look at is can we find the inverses of , , , and in an easy manner.  The answer must be yes or why do this.  Finding the inverse of Elementary matrices is a very important concept. Elementary matrices have ones on the diagonal and one non-zero entry in an entry below the diagonal.  The matrix Ea[1] can also be work like an elementary matrix but we will not go any into this in detail.   

>	IE1:=inverse(E[1]);

 

>	IE2:=inverse(E[2]);

 

>	IEa1:=inverse(Ea[1]);

 

>	IE3:=inverse(E[3]);

 

The inverses are found by changing the sign of the multiplier in the matrices.  This can be done in general.  We will look at this after we find the lower trangular matrix.  

Now that we know how to find the inverses we will look at what this really means.  Remember that order really does make a difference when we talk about multiplication of matrices.   

If = U  then =   or if we use the alternate we would have 

= which then leads to = .  The last thing we need to look at is the product of the lower triangular matrices to make sure they are lower triangular and can we find an easy way to find this product.  Remember we are looking for no more work than when we did Gauss elimination.  You may not think that is true at this time but we will look at it closer in a few minutes.   

Because we used the addrow command and then had to use the inverse command from the same package we will have to use the old way of multiplying matrices.  The good news is that we will not be using this to find the L in the LU decomposition.  Only to show you why the technique works.     

>	L:=evalm(IE1&*IE2&*IE3);

 

>	La:=evalm(inverse(Ea[1])&*inverse(E[3]));

 

What is the pattern?   

You need to notice that the multipliers with the signs changed from the Gauss Elimination steps will give you the numbers in the lower triangular matrix L as we will now discuss in more detail.  

Let us now do Gauss elimination steps and look at the multipliers to find the L matrix. 

>	A1:=addrow(addrow(A,1,2,-1/2),1,3,1/2);

 

>	U:=addrow(A1,2,3,2);

 

>	L:=diag(1,1,1);

 

The first addrow command in A1 tells us that 1/2 goes in the 2,1 position in the L matrix.   The second command  in A1 tells us that -1/2 goes in the 3,1 position in the L matrix.   Finally the last 2 in the U command tell us that -2 will go in the 3,2 position in the L matrix.  This can be accomplished by the following commands.  

>	L[2,1]:=1/2;

 

>	L[3,1]:=-1/2;

 

>	L[3,2]:=-2;

 

>

evalm(L);

 

>

L.U;

 

So you will remember what was just discussed we will go over it again.  Notice that we will never have to find an Elementary matrix to find either U or L.   

Look at the A1 command above in inner most addrow command said to multiply -1/2 time row one and add the results to row two.  The entry in the 2,1 position of the L matrix will then be 1/2 (change the sign of the multiplier 1/2 and place it in this position.  The outer command tells us to multiply 1/2 time row one and add it to row 3.  So we change the sign to -1/2 and place it in the 3,1 position of the L matrix. 

Finally the last operation tell us to multiply 2 times the second row of the A1 matrix and add it to the third row.  Again, we change the sign of the multiplier 2 to -2 and place it in the 3,2 position.  You can thus build the lower matrix L without doing any more work than Gauss elimination.  The upper matrix is the one that you obtained from Gauss elimination.  The building of the L matrix only works if you follow the Gauss Elimination process going from top to bottom and left to right.  Do not try to jump around.  Always check to see if your L*U = A.   

We now want to use the decomposition to solve Ax = b.    What we have is LUx = b or if we set 

Ux = y we can now write the system as Ly = b.     

>	b:=<<1,0,1>>;

 

>	y:=ForwardSubstitute(<L|b>);

 

>	x:=BackwardSubstitute(<U|y>);

 

Let us check our answer.  

>

A.x;

 

LU Decomposition for a Singular Matrix and Solution for a System with a Singular Coefficient matrix. (No row switches) 

>	B:=<<5,-3,-2,1>|<-5,3,2,-1>|<10,2,0,10>|<0,2,-1,1>>;

 

>	B1:=addrow(addrow(addrow(B,1,2,3/5),1,3,2/5),1,4,-1/5);

 

Even though there is no non-zero entry in the second through fourth column we can go the entry in the 2,3 position and zero out the entries.  Notice that we now are finished because we have an upper triangular matrix.   

>	U1:=addrow(addrow(B1,2,3,-1/2),2,4,-1);

 

>	U:=addrow(U1,3,4,-1/2);

 

(2.1)

 

>	L:=diag(1,1,1,1);

 

>

 

>	L[2,1]:=-3/5;

 

The entry in the second row first column is found by look at the inner most command of addrow(B,12,3/5).  The Gauss elimination step said to multiply row one by 3/5 and add it to row 2.  As in the first example to find the entry in the 2,1 position of the L matrix we change the sign of the multiplier 3/5 to -3/5 and place it in the 2,1 position.   

>	L[3,1]:=-2/5;

 

The second command the the B1 command line tells us to multiply 2/5 time the first row and all it to the third row.  Hence, the entry in the 3 row first column of the L matrix is -2/5.   

>	L[4,1]:=1/5;

 

The final command in the B1 command line tells us that we need to multiply row one by -1/5 and add it to row four.  This then tells us that the entry in the 4,1 position is 1/5.  Always change the sign of the multiply to obtain the entry in the L matrix.  

>	L[3,2]:=1/2;

 

>	L[4,2]:=1;

 

>	L[4,3]:=1/2;

 

(2.2)

 

>

evalm(L);

 

>	evalm(L&*U);

 

>

B;

 

(2.3)

 

We will now look at a system of equations Bx = c 

>	c:=<<5,1,0,5>>;

 

>	y:=ForwardSubstitute(<L|c>);

 

>	x:=BackwardSubstitute(<U|y>);

 

>

B.x;

 

>

 

Exercises 

The following exercises are to worked using Maple.  You are again to do the addrow commands to create the U matrix and then use the multipliers to obtain the L matrix.  There are no row switches in these problems.  1.8.1, 1.8.2, 1.8,4(This problem can not be done in Maple but must be done using text)  The inverse of A if it exists should be written in terms of L and U.  For part b a hint is to recall that L is lower triangular and has ones on the diagonal.  Also, recall how you take the determinate of the product of two matrices.  1.8.6  Use the E matrices to find the inverse now just the inverse command from Maple.  You can verify your answer by using the inverse command in Maple.  Problems are due September 26, 2007.  Again you can hand me a hard copy or e-mail me a copy.