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]

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 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