Homework assigned 10/29, due 11/6

Use maple where appropriate for the operation counts

Problem 1

Let L_1 be a lower triangular matrix which is ones on the diagonal, zeros everywhere else, but possibly nonzero in the 1st column.

Generally, let L_k be a lower triangular matrix which is ones on the diagonal, zeros everywhere else, but possibly nonzero in the k_th column.

Prove that the inverse of L_k is a matrix with identical structure except the entries below the diagonal have their signs switched. (Do this proof with matrix notation, keeping track of indices.)

Problem 2

Let L_k be as described above.

Let L be a lower triangular matrix which is ones, zeros everywhere else, but the k+1st, k+2nd, ... nth columns can have non-zero entries.

Prove that L_k * L is a lower triangular matrix which is ones, zeros everywhere else, but the kth, k+1st, k+2nd, ... nth columns have non-zero entries which are exactly those non-zero entries directly inherited from L and L_k.

Problem 3

Do an operation count for the LU decomposition algorithm I gave in class to find out how many operations are involved for an n x n matrix A. (operations are +,-,*,/, and <,>.)

Problem 4

Write an algorithm to solve Ly = b, where L is lower triangular. Do an operation count to find out how many operations are involved for an n x n matrix L.

Problem 5

What was the operation count to solve Ux = y for an upper triangular n x n matrix U?

Problem 6

What was the operation count to solve Ax = b for a n x n matrix A using Gaussian elimination?

Problem 7

If you're going to do it only once, which way is faster?

Problem 8

Given 10 vectors b_1, ... b_10 and A is 5 x 5, how many operations would it take to find the 10 solutions using Gaussian elimination each time? Using the LU decomposition?

Problem 9

Given an n x n matrix A, how many problems would you have to solve before the LU decomposition approach would be faster than doing Gaussian eliminations each time?

Computer problem 1

In matlab, if you type "flops(0)", it resets the flop counter. If you type "flops" sometime after that, it tells you how many flops were used since you last reset it.

Write a *.m file which given a number n returns a vector "icount". The program will first find 100 random n x n matrices, 100 random vectors b, and will find the resulting 100 solutions. You then make icount(1) = flops. Reset the flop counter. Now find a single n x n matrix. Perform its LU decomposition. Find 100 random vectors b, and find the solution x by first solving L y = b using low_solve.m and then solving U x = y using up_solve.m. Then make icount(2) = flops.

I've provided the routine up_solve.m, you have to write the routine low_solve.m (Let me know if you fail in copying and cannibalizing the routine up_solve.m to do this.)

You now have a *.m file which given n will return icount. Using this routine, create a long vector as follows:

>> [icount] = expt(2);
>> y(1,1) = icount(1); y(2,1) = icount(2);
>> [icount] = expt(3);
>> y(1,2) = icount(1); y(2,2) = icount(2);
>> [icount] = expt(4);
>> y(1,3) = icount(1); y(2,3) = icount(2);
>> [icount] = expt(5);
>> y(1,4) = icount(1); y(2,4) = icount(2);

Repeat this until you've performed the routine for the 10 x 10 matrices. (Hint: you might want to do this with another *.m file, if you're quicker with an editor than with typing.)

Let x = 2:10. Plot y(1,:) against x. >> hold on. Plot y(2,:) against x using a different line-type.

Explain your results.

Computer problem 2

I mentioned that LU decompositions can also be done when Gaussian elimination preferred pivoting. For a 3 x 3 matrix A, Gaussian elimination with pivoting would produce two permutation matrices P_1 and P_2 and two gaussian matrices L_1 and L_2 such that
L_2 P_2 L_1 P_1 A = U
This can be rewritten as follows
L P A = L'_2 L'_1 P_2 P_1 A = U
if L'_2 = L_2, and L'_1 = P_2 L_1 inv(P_2). More generally, for a 4 x 4 matrix A, Gaussian elimination with pivoting produces
L_3 P_3 L_2 P_2 L_1 P_1 A = U
which can be rewritten as follows
L P A = L'_3 L'_2 L'_1 P_3 P_2 P_1 A = U
if L'_3 = L_3, L'_2 = P_3 L_2 inv(P_3), and L'_1 = P_3 P_2 L_1 inv(P_2) inv(P_3)

1) Verify the above statements. (This is just a formal manipulation.)

2) How do we know that the new matrices L'_1 and L'_2 are actually lower triangular? If you can't prove it, at least go through a 3 x 3 example by hand and verify that they really are.

3) Matlab will give us these matrices: [L,U,P] = LU(A) where PA = LU . To solve A x = b, this means it's good enough if we solve PA x = Pb = b'. So first multiply b by P to find a new right hand side. We are now on familiar territory. PA x = b' ==> LU x = b'. First solve L y = b', then solve U x = y.

4) Now, take the matrix A from exercise 6 on page 242. Take b = [0.3097; 0.2838; 0.4742; 0.4716]. First use my LU solver:

>> [L,U] = LU(A);
>> [y] = low_solve(L,b);
>> [x] = up_solve(U,y);
>> max(A*x-b)

Then use the matlab solver with pivoting:

>> [L,U,P] = lu(A);
>> bp = P*b;
>> [y] = low_solve(L,bp);
>> [x] = up_solve(U,y);
>> max(A*x-b)

What did you find? What did you learn? Could you have used the matlab version of the lu docomposition routine? (What happens if you'd typed "[L,U] = lu(A)"?)