﻿ 数值计算代写 Numerical Computing代写 – 天才代写

# 数值计算代写 Numerical Computing代写

2022-05-15 11:48 星期日 所属： CS代写 浏览：410 ## Numerical Computing, Spring 2022

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

A QR,  where R = is m × n and Rˆ is n × n, (1)

A = Q , where R is n × n.  数值计算代写

In both cases Q is m × m, and if we write Q = [Qˆ Qˆ] have, in my notation,with Qˆ  m × n, we

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.(10pts) Recall that a square matrix Qis orthogonal1 if QTQ= I(and hence QT= Q−1 and so QQT= I, where Iis the identity matrix).

(a)Show that, given u with “u” = 1, the Householderreflector

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 thematrix where P is an orthogonal matrix, is also orthogonal.

1 As I mentioned, “orthonormal” would be a better name, but we are stuck with the standard name “orthogonal”.

(c)Showthat the product P (n)P (n1)  · · 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  数值计算代写 Now to solve the least squares problem minx Ax − b“, all we have to do is write

where c ∈ Rn and d ∈ Rmn, 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: 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) Modify the matlab function house in the text (p. 198) so that 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));

to

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)Howdoes 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 ..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.(20pts) The function house does not actually compute the Qor 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 Qor 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 we want Q QI  or Qˆ = QIˆ, not QT QT I.  When you return 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- tionresidual 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)Thenconfirm your answer computationally on large enough prob- lems that you can confirm your  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 onthe residual norm ||AQR|| or ||AQˆRˆ|| and the “orthogonality” residual norm ||QT QI|| 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 aijare 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.  数值计算代写

2 These 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)Modify the Householder algorithm and the original function house so that it works efficiently on such a matrix A: specifically, it shouldnot 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)Confirmyour answer computationally, using etime on large enough 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. 