数值计算代写 Numerical Computing代写

Numerical Computing, Spring 2022

数值计算代写 Is this the same mathematically?Do you get nearly the same output as previously on the example from the previous question?

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

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

instead of AG’s

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 Q is orthogonal¹ if Q^T Q = I (and hence Q^T = Q⁻¹ and so QQ^T = I, where I is the identity matrix).

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

H = I − 2uu^T

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 ⁽ⁿ⁾P ⁽ⁿ⁻¹⁾ · · P ⁽³⁾P ⁽²⁾P ⁽¹⁾, 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 min_x “Ax − b“, all we have to do is write

where c ∈ Rⁿ and d ∈ R^m⁻ⁿ, 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: 数值计算代写

±||a_k||e₁, 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 ±||a_k|| 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)Howdoes the number of flops (floating point operations) required by house, which was originally O(mn²), 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 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 Q^T b, but here we want Q = QI or Qˆ = QIˆ, not Q^T = Q^T 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 ||Q^T 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.²

(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 ||A−QR|| or ||A−QˆRˆ|| and the “orthogonality” residual norm ||Q^T 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 a_ijare 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.

合作平台：essay代写论文代写写手招聘英国留学生代写