MATH 191, SPRING 2023, WILKENING Programming Assignment 2: Due Wed, March 15
Combine the write-ups for parts 1 and 2 below into a single pdf file and submit it to grade- scope. Upload your code to gradescope as well. If your write-up is a jupyter notebook, upload it as a notebook in the “code” assignment on gradescope, and upload it as a pdf file in the “write-up” assignment.
Part 1: Here we take a closer look at LAPACK’s low level QR factorization routines. In python, you can access this code via
[QR,tau],R = scipy.linalg.qr(A,mode=”raw”).
The matrix QR contains both R and the Householder transformation data that we called v1, . . . , vn in lecture 16. It also returns the τj components from lecture 16. We can apply Q or QT to a vector or matrix using the LAPACK routine dormqr. For example, to apply QT to a vector b in python, you can do this:
y,work,info = scipy.linalg.lapack.dormqr(“L”, “T”, QR, tau, b, 1)
Here work is a “work array” from the dormqr fortran interface, and info is used to return an error message if one occurs. In matlab, it was once possible to access this raw form, but that feature seems to have been removed. So I wrote a simple implementation called qrRaw.m with the same functionality as scipy.linalg.qr and posted it to bCourses. I also wrote applyQT.m and applyQ.m to replace dormqr. E.g., here is the matlab command to toapplyQT tob
y = applyQT(QR,tau,b)
(a) Repeat HW5#10 using these routines. Here A = 2 −1 and b = −1. Your
2 −2 −1 taskistocomputeτ1 ∈R,τ2 ∈R,v1 ∈R3,v2 ∈R3 andR∈R2×2 suchthattheHouseholder
transformations H1 = I − τ1v1v1T and H2 = I − τ2v2v2T reduce A to upper triangular form, T R
Q A=H2H1A= 0 =R0.
Also solve Ax = b in the least-squares sense and compute the minimum value of ∥r∥ = ∥b − Ax∥. Compare the computed answer to your (or my) solution of HW5#10 to make sure you are calling the subroutines and interpreting the output correctly. Remember that LAPACK’s QR routine doesn’t actually store the leading 1 in the v vectors. You’ll probably want to set “format long” in matlab or “np.set printoptions(precision=14)” in python.
−1.4 0.16 −2.44 (b)Repeatpart(a)forA= 3.2 2.92 andb=0.72.
3.2 3.92 −0.28 1.6 −1.04 6.36
(c) One reason to use the v,τ format is that you can set some of the Hk’s to the identity matrix (by setting τk = 0) if the kth column of A(k) already points along ek. Check that
Date: March 3, 2023.
Code Help
2 MATH 191, SPRING 2023, WILKENING
matlab or python are doing this by seeing what v1 and τ1 are for these two 3 × 1 matrices A: 3 3
case 1: A = 0 , case 2: A = 0 . 0 1.0e−18
So this tiny perturbation leads to a reflection being used in case 2 instead of the identity in
case 1. Explain the values of v1 and τ1 that you get in case 2. (I.e., also do the calculation
with pencil/paper for that case.)
(d) Consider the 12 × 8 matrix Aij = 1 . Here i, j start at 1. If you use zero-based i+j −1
indexing, then Aij = 1 . This is a rectangular version of a Hilbert matrix. (Hilbert i+j +1
matrices are famously ill-conditioned.) Write a simple code to generate this rectangular Hilbert matrix and create a vector x0 of length 8 with entries that are random numbers uniformly distributed between 0 and 1. Then compute b = Ax0 on the computer and solve the least squares problem Ax = b using the QR factorization code above (call the result x1), and by solving the normal equations AT Ax = AT b (call the result x2). Make a table with columns given by [x0, x1, x2, x1 − x0, x2 − x0]. (The table is 8 × 5. Switch back to “format short” or “np.set printoptions(precision=3)” here.) Also compute the normwise relative errors ∥x1 − x0∥/∥x0∥ and ∥x2 − x0∥/∥x0∥ in the 2-norm.
Part 2: (Variant of Question 2.9 from page 95 of Demmel’s book, which is page 65 of demmel chap02.pdf). Let B be an (n + 1) × (n + 1) lower bi-diagonal matrix of the form
b1 a1 .
b .. B=2
.. .an−1
Derive an algorithm for computing κ1(B) = ∥B∥1∥B−1∥1 exactly (ignoring roundoff.) Your
algorithm should be as cheap as possible; it should be possible to do using no more than 2n
additions, n + 1 multiplications, n + 1 divisions, 2n + 1 absolute values, and 2n comparisons.
(Anything close to this is acceptable. Note that ∥B−1∥1 is the maximum absolute column
sumofB−1,andthejthcolumnofB−1,whichwedenotebyyj ∈Rn+1 for0≤j≤n,can
be obtained by forward solving Byj = ej with ej the jth elementary unit vector in Rn+1.
The jth absolute column sum is just ∥yj∥1 = n (yj)i, where we note that (yj)i = 0 for i=j
i < j. The challenge of this problem is finding a way to re-use intermediate calculations so that you can compute ∥yn∥1, ∥yn−1∥1, ∥yn−2∥1, . . . , ∥y0∥1 cheaply, in that order, keeping track of the largest one. You also have to compute ∥B∥1, but that is easy.)
Implement your algorithm with the calling interface
kappa1(a,b)
where a is a vector of length n+1, b has length n, and the return value is κ1(B). (Depending on whether you use matlab or python, either b or a is indexed awkwardly in the matrix above; make sure you correctly match up the entries of the input vectors with the variables of your derivation.) Use your method to compute κ1(Bn) for (n+1)×(n+1) matrices Bn with 2’s on the main diagonal and 1’s on the subdiagonal. Plot 3−κ1(Bn) versus n for 1 ≤ n ≤ 50 with the y-axis scaled logarithmically. Also report κ1(Bn) with n = 20 to 14 digits of precision. Repeat this calculation with 1’s on the main diagonal and 2’s on the subdiagonal. This time plot κ1(Bn) versus n for 1 ≤ n ≤ 50 with the y-axis scaled logarithmically, and again report κ1(Bn) for n = 20 to 14 digits. Include the derivation of your algorithm in the write-up.
Programming Help