Condition Number Estimate

And Topics from Chapter 2

To estimate the condition number of a matrix using the algorithm in part (a) of Computer Problem 2.4 (page 101, with a brief discussion on page 59 of the text), you need to solve a lower triangular system (Utv = c, where your algorithm must determine c as it goes) and then an upper triangular system (Lty = v). In general, if M is a lower triangular matrix (so M[i][j] = 0 for j > i), and the system is Mx = b, then row i of M represents the equation

M[i][0]*x[0] + M[i][1]*x[1] + ... + M[i][i-1]*x[i-1] + M[i][i]*x[i] = b[i]

Hence, if

sum = M[i][0]*x[0] + M[i][1]*x[1] + ... + M[i][i-1]*x[i-1]
then,

x[i] = (b[i] - sum) / M[i][i]

Most often in implementations of algorithms for solving systems, the vector b[i] is stored in the last column of the Matrix M, so b[i] would be M[i][n] where n = the number of rows of the coefficient matrix (which is also the number of columns of the coefficient matrix but the augmented matrix has one additional column - so in Matrix.cc the index of the last column is either numRows or numCols - 1). In problem 2.4, you don't know what b is -- you want to choose +1 or -1 for each entry so that the component of the solution x is as large (in magnitude) as possible. So, for each row, your program should

  1. Compute the sum (or as many of you did compute the negative of the sum and make appropriate adjustments to the sign in step 2).
  2. Choose the entry b[i] (or c[i] using the notation of the book -- but this does not need to be named anywhere in the program!!!) to make b[i] - sum have the largest magnitude -- so if sum is positive choose b[i] to be -1 otherwise choose b[i] to be 1.

An equation similar to the above can be written for solving an upper triangular system. In problem 2.4, the upper triangular system has ones on the diagonal so the last step of dividing by the entry on the diagonal is not needed (just as in solving the original system in my code there is no division when forward substituting using L).

Dealing with the Row Switches: The most efficient and easiest way to deal with the fact that there have been "virtual" row switches (which are kept track of in the row vector) is to use the LU matrix as it is!!!! You just need to think about what the indices are for the transpose of U and of L. As an example, suppose that in doing Gaussian elimination on a 3 by 3 matrix, the first and second row were switched in the first step then the second and third row were switched in the next step. Then (using C++/Java indices), row[0] = 1, row[1] = 2, and row[2] = 0. Suppose that physically LU (this is the way the matrix is physically stored in the computer) is as follows: 

          1  2  3
          4  5  6
          7  8  9
So, logically it is 
          4  5  6     (logical row 0 is physically in row 1 of LU)
          7  8  9
          1  2  3
So, the logical transpose is 
          4  7  1
          5  8  2
          6  9  3
 
and hence, Ut is

 
          4  0  0
          5  8  0
          6  9  3


and Lt is

          1  7  1
          0  1  2
          0  0  1


So, Ut[i][j] = LU[_______________][___________________]
(fill in the blank -- if you can't see it with i's and j's first use
some actual numbers).  Similarly, the entries in Lt can
be gotten directly from LU using correct indexing and can be used
in the backsubstitution algorithm to get a solution vector y.


To complete the approximation of the norm of A inverse you need to
solve Az = y (where y is your solution above).  To use the forback()
function already in the Matrix class (this is the one that SHOULD be
used when an LU decomposition has already been obtained), you need
to be sure your solution is in the last column of LU.  You can
put it there as you solve or copy it into that column -- so
LU[i] gets y[i]).


Copying LU or its transpose into a new matrix:  
There were many variations in your programs that involved making one
or more new matrices -- many of you got carried away and used
the vector row in too many places or the wrong places. 
If you found the logical transpose (as above --
the one that splits into the U and L) and worked with it (so, after
you found the transpose you DO NOT use the vector row in applying
the forward and backward substitution algorithms), the
components of your solution are not in the same order as the original.  For
example, if you get a solution vector y, then y[i] does not
go with the original row i -- it is a solution of original row[i].  
Hence when you copy your solution into
the last column of the original LU, LU[row[i]][numRows] gets y[i].  



Topics from Chapter 2 that will be on Test #2