Lab assignment L1

LU factorization

For due dates see the calendar web page. For hand in details the lab main page.

Objective

To practically implement and analyse the routines involved in LU factorization and thereby get a hands-on knowledge of the numerics of factorization, how to write numerical routines in Matlab and/or Numpy, particularly how to access sub-parts of vectors and matrices, and how to write clean, correct vectorized code.

Procedure

The goal is to compose a linear equation solver from the following parts. Here the solver should work with the explicit matrices as shown in Heath Section 2.4.3, not the for loops in alg 2.4.

We provide a fwdSubst function (below). Experiment with it alone for some triangular matrices.
Modify it to do a backSubst function.
Using Octave/Matlab built in [L, U] = lu(A) function (or NumPy's) solve some equation systems of your own design with your fwd and bkw subst functions per point 6 below.
We provide a function [M,L] = elimMat(A,k) which given matrix A, index k computes the elimination matrix M and L as specified on Heath p67 .
- For each row the elimMat function draw the vectors and matrices involved indicating dimensions. You can print the code and draw by hand next to the rows.
Write an myLU function which combines L1, L2 ... M2 M1 A into L and U
Show how to use your programs to solve Ax=b by method of LU and then solving the trianglular L and U systems( p. 68 )
1. Factorization: L U = A
2. fwdSubst to solve for y in: Ly=b
3. backSubst to solve for x in: Ux = y
Test your program with e.g. example 2.13 and one test of your own design. Write a m-file script to execute the the example and comment on what's happening in each step - you can have a program flag that when enabled shows salient vectors and matrices, or (cleaner code) simply two versions of your programs where one prints and the other not.
For your own test aim to solve a real world problem (ie like the pipeline temp problem in class, bridge problem, car problem or tomography problem in this lab). Mathematicallhy model the problem, derive the linear equations system and show how to solve it. Does it generalize to higher dimensions (ie like the pipeline can have more interior points, or the bridge more elements)?

It is not required to implement pivoting in this lab, but think about: What happens in your program if diagonal element 0? If the system is ill conditioned?

(NOT required) In particular, for those who want some challenge, think about the following problem: Given a matrix A and vector b, decide if the linear system Ax=b is consistent or not; if yes, then deliver one solution (an arbitrary one if there are multiple ones).

Test your program on the matrix, say,

A = [1, 2, 3, 4, 0; 0, 0, 0, 2, 3; 0, 0, 0, 2, 1; 0, 0, 0, 0, 0; 0, 0, 0, 0, 0];

b = [1; 2; 3; 0; 0];

Style

Show how to use matlab matrix operations to write the above in a clear compact form with a minimum number of for loops. (e.g. copying the for-loop based algorithms 2.2 and 2.3 will not yield good marks). Write your code in a readable, commented and suitably indented fashion. (There is a matlab mode for emacs you can download)

Hints

First understand and experiment with the supplied algorithms code below. It's helpful to do this in an interactive way - print/show the intermediate vectors and matrices. Modify it to do back substitution. Then write the back substitution and myLU routines. myLU can be written recursively much the same way as the fwdSubst, or written with one "for" loop. You need no loop in computing the elimination matrix and at most one loop (or recursion) in computing the LU and forward, backward substitutions. Achieving this depends mostly on skills in linear algebra to do operations on vectors and matrices as a whole instead of element by element and on understanding Matlab's matrix elements adressing as described in the first part of the tutorial. This way each function (elim matrix, lu and substitutions) can each be written in just 5-10 lines of code.

Code:

% elimMat   Compute elimination matrix for k'th row
% cmput 340

function [M,L] = elimMat(A,k)
  [p,q] = size(A);
  M = eye(p,q);
  e = M(k,:);
  m = zeros(q,1);
  m(k+1:n) = -A(k+1:q,k)/A(k,k);
  L = M - m*e;
  M = M + m*e; 





% fwdSubst  Forward substitution to solve Ly=b by tail recursion
% cmput 340,

function y = fwdSubst(L,b,k)

  [m,n]=size(L);
  if ~exist('k')  % If first call no k param given, but k=1
    k=1;
  end
  
  y=b(k)/L(k,k)
  if k < n % Recursion step
    l = [zeros(k,1);L(k+1:m,k)]
    y = [y;fwdSubst(L,b-y*l,k+1)]
end

Python


import numpy as np

# find the M and L which can be used to eliminate entries of a matrix
def elimMat(A, k):
    elimination_col = -A[k+1:, k] / A[k, k] # extract the column used for elimination

    # construct the M by adding in the elimination_col
    M = np.eye(*A.shape)
    M[k+1:, k] = elimination_col

    # inverse of M is the negative of the corners in this case
    L = np.eye(*A.shape)
    L[k+1:, k] = -elimination_col

    return M, L

def fwdSubst(L, b, k=0):
    _, cols = L.shape

    y = b[k]/L[k, k]
    if k < cols - 1: # recursion step
        l = np.concatenate([np.zeros(k+1), L[k+1:, k]])
        recursive_ys = fwdSubst(L, b - y * l, k + 1)

        return np.concatenate([[y], recursive_ys])
    else:
        return np.array([y])

Misc. Q and A:

Subject: Re: L1 myLu(...)
References: 
--text follows this line--
Paul Cartledge  writes:

> If I could get some clarification on the myLU(...) function, that would be
> great.
> 
> I believe it is wrong to assume that we will only be getting 3x3 matrices
> for this lab, so I am wondering how we are to send in parameters to
> myLU().

Obvioulsy your function should not just handle 3x3 matrices.
Look back to the reading list for lab 1, and make sure you
have read and understood the parts of the Matlab tutorial that
was covered. In particular, ch 3 has the answer for how to extract
properties (such as size) from a passed in matrix. Especially pay
attention to table 4. Then take a look at the sample fwdSubst function
supplied and see how parameters are used.

> 
> If there's a varying number of matrices to be sent in, is there a

Why would you like to pass in a varying number of matrices?

> way to store the intermediate matrices we need for computation obtained
> from elimMat() so that we can send them into myLU().

If you use recursion you have the benefit of a new set of intermediate
variables at each call. Again study the tutorial, then the sample routine.

/Martin

----------------------------------------

Subject: Re: Leading Diag Entries of Unreduced Portions being Zero
References: 
--text follows this line--
Don Rivers  writes:

> The textbook talks about the possibility of there being a leading zero in
> a part we've yet to reduce.  Is there an elegant way for us to interchange
> rows, 

1. Yes, this is essentially the same as piviting

2. You are not required to implement this in the lab, but in preparation
for the lab test you should be aware of the issue from your readings 
(as you obvioulsy are) so you can answer oral questions.

> or will all of the matrices that we will be tested on not have
> leading zeros in parts that we are working on?

On a test matrix with a leading 0 you will have to point out why 
your code doesn't handle it, identify the problem, and the particular
place in the code.

/Martin

------------------------------------

Subject: Re: NxN, or MxN?
References: 
--text follows this line--
Don Rivers  writes:

> I don't know if this has been answered, but are only going to be tested on
> Square Matrices, or on nonSquares as well?

Not required to handle non-square matrices.

> 
> I notice that matlab's lu() complains if you pass it a nonsquare matrix,
> so do we have a similar requirement?

Actually, it is possible to compute an L U factorization through
elimination even for an mxn matrix. It is just not implemented
in matlab 5 (but was fixed in matlab 6). 

/Martin