Numerical Computing, Spring 2022
Homework Assignment 5
QR Factorization via Householder Reduction
Be sure to read A&G, section 6.2 and 6.3 (you can skip the subsections on Gram-Schmidt and the SVD) as well as my notes on QR factorization before starting this homework. Please note that this homework uses the notation in my notes, with
Copyright By PowCoder代写 加微信 powcoder
A=QR, whereR=Rˆ ism×nandRˆisn×n, (1) 0
instead of AG’s
A=QR0, whereRisn×n. InbothcasesQism×m,andifwewriteQ=[Qˆ Q ̃],withQˆm×n,we
have, in my notation,
A = QR = [Qˆ Q ̃]Rˆ = QˆRˆ 0
which is the reduced (economy-size) QR factorization, A = QˆRˆ, in my notation. In AG notation, the reduced factorization is A = QˆR.
1. (10 pts) Recall that a square matrix Q is orthogonal1 if QT Q = I (and hence QT = Q−1 and so QQT = I, where I is the identity matrix).
(a) Show that, given u with ∥u∥ = 1, the Householder reflector H = I − 2uuT
is an orthogonal matrix. (The norm is always the 2-norm when we are talking about least squares, QR and orthogonal matrices.)
(b) Show that the matrix
where P is an orthogonal matrix, is also orthogonal.
1As I mentioned, “orthonormal” would be a better name, but we are stuck with the standard name “orthogonal”.
(c) Show that the product P(n)P(n−1)···P(3)P(2)P(1), where each P(k) is orthogonal, is also orthogonal.
2. (10 pts) The equation in the middle of p. 188 of A&G which I’ve marked with a sticky note shows how we can use the QR factorization to solve least squares problems. Rewriting this in my notation, this becomes
∥ b − A x ∥ = ∥ Q T ( b − A x ) ∥ = ∥ Q T b − Q T Q R x ∥ = ∥ Q T b − Rˆ x ∥ . 0
Now to solve the least squares problem minx ∥Ax − b∥, all we have to
do is write
Q T b = Qˆ T b = Qˆ T b = c Q ̃ T Q ̃ T b d
where c ∈ Rn and d ∈ Rm−n, and then solve the square triangular system Rˆx = c. Explain exactly what is wrong with the following statement using the reduced QR factorization A = QˆR instead:
∥b − Ax∥ = ∥QˆT (b − Ax)∥ = ∥QˆT b − QˆT QˆRˆx∥ = ∥QˆT b − Rˆx∥.
If this were true it would mean every least squares problems has a solution with a zero residual, since we could just solve the triangular system Rˆx = QˆT b to make the right-hand side zero.
3. (10pts)Modifythematlabfunctionhouseinthetext(p.198)sothat after computing the kth Householder reflector, instead of multiplying it onto subcolumns k,…,n, you multiply it only onto subcolumns k+1,…,n, since we know what the answer will be for subcolumn k: ±∥ak∥e1, which can be stored explicitly. To figure out which sign to use, you need to read the algorithm description in the text. (Actually, you only need to store ±∥ak∥ in the diagonal entry A(k,k), because the zeros below it will be overwritten by the reflector information anyway.) Submit the modified code and, using this input, submit output showing that you still get approximately the same answer A (it might not be exactly the same to 16 digits, but it should be close: the easiest way to check is to compute, say, norm(Aout1 – Aout2), which should be close to machine precision.
4. (20 pts) Going back to the original function house, let’s modify it differently, changing
A(k:m,k:n)-2*u*(u’*A(k:m,k:n));
A(k:m,k:n)-2*(u*u’)*A(k:m,k:n));
(a) Is this the same mathematically? Do you get nearly the same output as previously on the example from the previous question?
(b) How does the number of flops (floating point operations) required by house, which was originally O(mn2), change, in terms of m and n? Work this out before running the program, explaining your reasoning carefully.
(c) Confirm your answer to the number of flops computationally, using tic…toc or etime, on random problems A=randn(m,n) with m, n large enough that you can see the difference. To an- swer this systematically, plot the original and the new running time as a function of n, fixing m = 10n, for both the original and the modified program, and explain how this supports your answer in part (b). As always, be sure your plot has appropriate legends, axis labels and title, and use a plot command such as plot(timevec,’r-*’) which is easy to read – or semilogy or loglog if one of those seems better.
5. (20 pts) The function house does not actually compute the Q or Qˆ of the full or reduced (economy-sized) QR factorization; instead it applies n Householder reflectors to compute the R part of the QR factorization (which overwrites the upper triangle of A), and it returns the u vectors defining those Householder reflectors by overwriting the part of A that lies below the diagonal, and, as this is not quite enough space to store all those vectors, putting the first component of each u vector in the output argument q.
(a) Write another matlab function myqr that has the same basic
syntax as matlab’s qr, so that [Q,R]=myqr(A) returns the full
Q, R factors while [Qhat,Rhat]=myqr(A,0) returns the reduced
factors Qˆ, R. It should call house (not matlab’s qr!) to do
the Householder reduction. Then to obtain Q or Qˆ you need to ˆ
efficiently apply the Householder transformations to I or I, in
the correct order: see p. 5 of the notes. You can conveniently
compute Iˆ by eye(m,n). It may help to look at the code lsol
on p. 198, which computes QˆT b, the first n rows of QT b, but here ˆˆTT
wewantQ=QIorQ=QI,notQ =Q I.Whenyoureturn R or Rˆ, you will need to explicitly set the desired zero entries to
zero, since house overwrites the lower triangular part of A with the Householder reflector information.
(b) Check that you got the right answer by computing the factoriza- tion residual norm ∥A − QR∥ or ∥A − QˆRˆ∥, and the “orthogonal- ity” residual norm ∥QT Q − I∥ or ∥QˆT Qˆ − I∥. As long as these are both around machine precision, and R or Rˆ is triangular, you have a “good” answer. However, you cannot expect it to be close to the output of matlab’s qr, first because the signs of the diag- onal entries of R or Rˆ might be different, and secondly because of conditioning issues.2
(c) How many flops does your code require in each case (full QR and reduced QR), in terms of m and n? As in the previous question, work this out before running the program, explaining your reasoning.
(d) Then confirm your answer computationally on large enough prob- lems that you can confirm your claim. Compare the running time to matlab’s qr, which will certainly be faster: how much faster? Is the ratio of the running times of myqr and qr constant as you increase n (again keeping m = 10n)? As previously, plot the running times on a nicely labelled graph.
6. (10 pts) Modify the original function house to compute u using the wrong choice of sign, and compare the resulting output of myqr on a problem that you choose to create cancellation in the computation of the first Householder reflector. What effect does this have on the residual norm ∥A−QR∥ or ∥A−QˆRˆ∥ and the “orthogonality” residual norm ∥QT Q−I∥ or ∥QˆT Qˆ−I∥? Try varying your input choice to create more or less cancellation and see what effect this has on the residuals. As with all the questions, submit both code and output, and in this case, also the input that you chose for A.
7. (20 pts) Suppose that A is an m × n matrix with the property that its entries aij are zero for i ≥ j +p (and i ≤ m), for some integer p (the “lower bandwidth” of the matrix). Thus, the number of nonzero subdiagonals of the matrix, including the main diagonal, is p.
2These are subtle and we won’t go into them, but if you are interested a good reference is Lecture 16 of the book by Trefethen and Bau.
(a) ModifytheHouseholderalgorithmandtheoriginalfunctionhouse so that it works efficiently on such a matrix A: specifically, it should not access the zero entries below the p nonzero subdiag- onals. You should provide p as an additional input parameter. Test your code on this input (for which p = 5).
(b) How many flops does your solution require in terms of m, n and p? As in the other questions, work this out before running the program, and explain your reasoning.
(c) Confirmyouranswercomputationally,usingetimeonlargeenough problems that you can confirm your claim. To construct test problems A use triu(randn(m,n),-(p-1)), with m = 10n and p = 10, and plot the running time as a function of n.
As well as the answers to the questions, remember to submit supporting code (with comments) and output.
Remember that if you work with another student, in addition to ac- knowledging that fact you must each write your own codes and your own solutions which cannot be the same.
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com