Linear Algebra Powertool

The LU Decomposition

OnLine 3.5

Worksheet by Russell Blyth.

> restart: with(linalg): with(plots): with(plottools):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

Outline

This worksheet investigates the LU decomposition of a matrix.

The basic objectives are:

• 1) Compute the LU decomposition, if possible.
• 2) See how elementary matrices are involved in the process.
• 3) Understand how the LU decomposition is useful in solving systems of equations.

Finding an LU decomposition

First define a random 3 x 4 matrix.

> A := randmatrix(3,4,entries=rand(-5..5));

Let's keep track of the elementary matrices corresponding to the row operations used, as well as the cumulative product of the elementary matrices, which we can do by starting with the matrix A augmented by the appropriate identity matrix:

> Iden:=evalm(array(identity, 1..3,1..3));
AAug := augment(A,Iden);

Now we apply elementary row operations to the augmented matrix. We hope to add multiples of rows to lower rows to get a row echelon form. We exhibit the elementary matrices at each step.

(Note that the matrix A is not truly "random" - every time the worksheet is started, it produces the same matrix A. Thus we know in advance which row operations are required in this worksheet. This is handy for teaching purposes. If you want a truly "random" matrix, you need to precede the randmatrix command with a "randomize()", which generates a seed for the random functions in Maple, based on the computer clock.)

> AAugR := addrow(AAug,1,2,1/4);
E := addrow(Iden,1,2,1/4);

> AAugR := addrow(AAugR,1,3,-1);
E := addrow(Iden,1,3,-1);

> AAugR := addrow(AAugR,2,3,5*4/17);
E := addrow(Iden,2,3,5*4/17);

AAugR is in row echelon form. Note we don't bother to force leading ones. The matrix U can now be extracted:

> U := delcols(AAugR,5..7);

The matrix B referred to in class can also be extracted:

> B := delcols(AAugR,1..4);

The matrix L we seek is just the inverse of this B:

> L := inverse(B);

Now note the remarkable relationship between the entries of L, the elementary matrices, and the multipliers used in the row reductions! The matrix L can in fact be written down without actually computing the inverse of B.

Check that A = LU

> evalm(L &* U); evalm(A);

Exercise:

1) Repeat the above for a random 4 x 4 matrix AA (be careful not to reuse any of the variable names used previously, since some get reused below!)

Now let's see how to use the LU decomposition to solve a system AX = Y

Let's generate a RHS vector Y

> Y := randmatrix(3,1,entries=rand(-5..5));

Solve LZ = Y for Z; note this amounts to a simple forward-substitution (no further row reductions). We show the augmented matrix and have Maple compute.

> augment(L,Y);
Z:=linsolve(L,Y);

Then solve UX = Z for X; note this amounts to a simple back-substitution (no further row reductions). We show the augmented matrix and again have Maple compute the solution for us.

> augment(U,Z);
X := linsolve(U,Z);

(Unfortunately Maple does not always select the parameter the way we do.)

Exercise:

2) Generate a random 4 x 1 vector YY, and then repeat the steps above to solve (AA)(X) = YY.

>

3) What advantage would the LU decomposition have for computer computation of solutions to linear systems?

4) Use Maple's help command to find out about the LUdecomp command in the linalg package, Use the command to check your LU decomposition of AA in exercise 1.

>

>

>

>