6. Midterm #1 ProgrammingΒΆ

This is the take-home programming part of Midterm #1. You are not allowed to collaborate with any other student, failure to comply will result in a score of zero on the entire exam.

Please turn in your code along with your solution. Submissions should be made on Sakai.

Problem

Write a routine for estimating the condition number \(\kappa(A)\) of a matrix \(A\) using the \(L_1\)-norm. You will need to compute \(\lVert A\rVert_1\), which should be easy, and estimate \(\lVert A^{-1}\rVert_1\), which is more challenging. Note that we only want to estimate \(\lVert A^{-1}\rVert_1\), not compute it exactly (up to finite precision). One way to estimate \(\lVert A^{-1}\rVert_1\) is to choose a vector \(y\) such that the ratio \(\lVert z\rVert/\lVert y\rVert\) is large, where \(z\) is the solution to \(Az=y\). Try two different approaches to choosing \(y\):

  1. Choose \(y\) as the solution to the system \(A^Ty = c\), where \(c\) is a vector each of whose components is \(\pm 1\), with the sign for each component chosen by the following heuristic. Using the factorization \(A = LU\), the system \(A^Ty=c\) is solved in two stages, successively solving the triangular systems \(U^Tv = c\) and \(L^Ty = v\). At each step of the first triangular solution, choose the corresponding component of \(c\) to be either \(1\) or \(-1\), depending on which will make the resulting component of \(v\) larger in magnitude. (You will need to write a custom triangular solution routine to implement this.) Then solve the second triangular system in the usual way for \(y\).

    The idea here is that any ill-conditioning in \(A\) will be reflected in \(U\), resulting in a relatively large \(v\). The relatively well-conditioned unit triangular matrix \(L\) will then preserve this relationship, resulting in a relatively large \(y\).

  2. Choose five different vectors \(y\) randomly and use the one producing the largest ratio \(\lVert z\rVert/\lVert y\rVert\).

Test both the approaches described in parts (a) and (b) on each of the following two matrices:

\[\begin{split}A_1 &=& \enspace \left[\begin{array}{ccc} -10 & 7 & 0 \\ 5 & -1 & 5 \\ -3 & 2 & 6 \end{array}\right] \\ A_2 &=& \enspace \left[\begin{array}{ccc} 92 & 66 & 25 \\ -73 & 78 & 24 \\ -80 & 37 & 10 \end{array}\right]\end{split}\]

Report your results using these two methods. Check the quality of your estimates using the function np.linalg.cond(A,1) and report your results.

For solving the linear system \(Ax = b\), you can use the following code:

import numpy as np
x = np.linalg.solve(A,b)

Note: The only place where you should not use this routine is when solving the system \(U^Tv = c\), where you need to adjust each component of \(c\) depending on the magnitude of the corresponding component of \(v\).

For computing the \(LU\) decomposition of a matrix \(A\), you can use the following code:

from scipy.linalg import lu
A = np.array([[10.,-7.,0.],[5.,-1.,5.],[-3.,2.,6.]])
P,L,U = lu(A,permute_l=False)

Both the matrices \(A_1\) and \(A_2\) have been chosen such that the permutation matrix \(P = I_{3\times 3}\), the identity matrix. So you do not need to worry about pivoting in this problem. Finally, to generate a random vector \(y\) using NumPy, you can use the following code:

>>> np.random.rand(3,1)
array([[ 0.16823198],
       [ 0.91167363],
       [ 0.34837631]])

Needless to say, you cannot use the function np.linalg.norm in this problem to compute the norm of a matrix. However, you can (and should) use it to compute the norm of a vector. Please use double precision, which is the default precision in Python.