MATH2071: LAB #7: FactorizationsIntro duction Exercise 1Orthogonal Matrices Exercise 2A Set of Questions Exercise 3The Gram Schmidt Method Exercise 4Gram-Schmidt Factorization Exercise 5Householder Matrices Exercise 6Householder Factorization Exercise 7The QR Method for Linear Systems Exercise 8The singular value decomposition Exercise 9Exercise 101 IntroductionWe have seen one factorization–PLU–that can be used to solve a linear system provided that the system issquare, and that it is nonsingular, and that it is not too badly conditioned. However, if we want to handleproblems with a bad condition number, or that are singular, or even rectangular, we are going to need tocome up with a different approach. In this lab, we will look at two versions of the QR factorization:A = QRwhere Q is an “orthogonal” matrix, and R is upper triangular, and also the “singular value decomposition”(SVD) that factors a matrix into the formA = USVT.We will see that the QR factorization can be used to solve the same problems that the PLU factorizationhandles, but can be extended to do several other tasks for which the PLU factorization is useless. In situationswhere we are concerned about controlling roundoff error, the QR factorization can be more accurate thanthe PLU factorization, and is even a little easier to use when solving linear systems, since the inverse of Qis simply QT.We will also see that the SVD can be used to compress data, to identify the range and nullspace of alinear operator, and even to solve non-square systems.You may find it convenient to print the pdf version of this lab rather than the web page itself. This labwill take three sessions.2 Orthogonal MatricesDefinition: An “orthogonal matrix” is a real matrix whose inverse is equal to its transpose.By convention, an orthogonal matrix is usually denoted by the symbol Q. The definition of an orthogonalmatrix immediately implies thatQQT= QTQ = I.From the definition of an orthogonal matrix, and from the fact that the L2vector norm of x can be definedby:||x||2=q(xTx)and the fact that(Ax)T= xTATyou should be able to deduce that, for orthogonal matrices,||Qx||2= ||x||21If I multiply x by Q and its L2norm doesn’t change, then Qx must lie on the circle whose radius is x. Inother words, the only thing that has happened is that Qx has been rotated around the origin by some angle,or reflected through the origin. That means that an orthogonal matrix represents a rotation or reflection.When matrices are complex, the term “unitary” is an analog to “orthogonal.” A matrix is unitary ifUUH= I, where the H superscript refers to the “Hermitian” or “conjugate-transpose” of the matrix. InMatlab, the prime op erator constructs the Hermitian and the dot-prime op erator constructs the transpose.A real matrix that is unitary is actually orthogonal.3 A Set of QuestionsSuppose we have N vectors xi, each of dimension M. Here are some common and related questions:• Can we determine if the vectors are linearly independent?• Can we determine the dimension of the space the vectors span? (This is the rank.)• Can we construct an orthonormal basis for the space spanned by the vectors?• Can we construct an orthonormal basis for the space of vectors perpendicular to all the xi?• Given a new vector xN+1, can we determine whether this vector lies in the space spanned by the othervectors?• If a new vector lies in the space, what specific linear combination of the old vectors is equal to the newone?• If we create an M by N matrix using the vectors as columns, what is the rank of the matrix?• How can we solve a rectangular system of linear equations?Since we actually plan to do these calculations numerically we also need to b e concerned about how wedeal with numerical error. For instance, in exact arithmetic, it’s enough to say a matrix is not singular. Innumerical calculations, we would like to be able to say that a particular matrix, while not singular, is ”nearlysingular”.4 The Gram Schmidt MethodThe “Gram Schmidt method” can be thought of as a process that analyzes a set of vectors X, producingan “equivalent” (and possibly smaller) set of vectors Q that span the same space, have unit L2norm, andare pairwise orthogonal. Note that if we can produce such a set of vectors, then we can easily answer manyof the questions in our set of problems. In particular, the size of the set Q tells us whether the set X waslinearly independent, and the dimension of the space spanned and the rank of a matrix constructed from thevectors. And of course, the vectors in Q are the orthonormal basis we were seeking. So you should believethat being able to compute the vectors Q is a valuable ability. In fact, the Gram-Schmidt method can helpus with all the problems on our list.We will be using the Gram-Schmidt process to factor a matrix, but the process itself pops up repeatedlywhenever sets of vectors must be made orthogonal. You may see, for example, it in Krylov space methodsfor iteratively solving systems of equations and for eigenvalue problems.In words, the Gram-Schmidt process goes like this:1. Start with no vectors in Q;2. Consider x, the next (possibly the first) vector in X;23. For every vector qialready in Q, compute the projection of qion x (i.e., qi· x) and subtract it fromx to get a vector orthogonal to qi;4. Compute the norm of (what’s left of) x. If the norm is zero (or too small), discard x from Q; otherwise,divide x by its norm and move it from X to Q.5. If there are more vectors in X, return to step 2.6. When you get here, X is empty and you have an orthonormal set of vectors Q spanning the same spaceas the original set X.Here is a sketch of the Gram Schmidt process as an algorithm. Assume that nxis the number of x vectors:nq= 0 % nqwill become the number of q vectorsfor k = 1 to nxfor ` = 1 to nq% if nq= 0 the loop is skippedr`k= q`· xkxk= xk- r`kq`endrkk=√xk· xkif rkk> 0nq= nq+ 1qnq= xk/rkkendendYou should be able to match this algorithm to the previous verbal description. Can you see how the L2norm of a vector is being computed?Exercise 1:(a) Implement the Gram-Schmidt process in an m-file called gram schmidt.m. Your function shouldhave the signature:function Q = gram_schmidt ( X )% Q = gram_schmidt ( X )% more comments% your name and the dateAssume the x vectors are stored as columns of the matrix X. Be sure you understand this datastructure because it is a potential source of confusion. In the