f2-10

2.10QR Decomposition 91Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).enddo 11 x(i)=sum/p(i) enddo 12 do14i=n,1,-1 SolveLT·x=y. sum=x(i) do13k=i+1,n sum=sum-a(k,i)*x(k) enddo 13 x(i)=sum/p(i) enddo 14 return END Atypicaluseof choldcandcholslisintheinversionofcovariancematricesdescribing thefitofdatatoamodel;see,e.g., §15.6. Inthis,andmanyotherapplications,oneoftenneeds L−1. The lower triangle of this matrix can be efficiently found from the output of choldc: do13i=1,n a(i,i)=1./p(i) do12j=i+1,n sum=0. do11k=i,j-1 sum=sum-a(j,k)*a(k,i) enddo 11 a(j,i)=sum/p(j) enddo 12 enddo 13 CITED REFERENCES AND FURTHER READING: Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag), Chapter I/1. Gill,P.E., Murray, W.,andWright, M.H. 1991, NumericalLinearAlgebraandOptimization , vol. 1 (Redwood City, CA: Addison-Wesley), §4.9.2. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §5.3.5. Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §4.2. 2.10 QR Decomposition There is another matrix factorization that is sometimes very useful, the so-called QR decomposition , A=Q·R (2.10.1 ) HereRis upper triangular, while Qis orthogonal, that is, QT·Q=1 (2.10.2 ) whereQTis the transpose matrix of Q. Although the decomposition exists for a general rectangularmatrix,weshallrestrictourtreatmenttothecasewhenallthematricesaresquare,with dimensions N×N. 92 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Like the other matrix factorizations we have met ( LU, SVD, Cholesky), QRdecompo- sition can be used to solve systems of linear equations. To solve A·x=b (2.10.3 ) first form QT·band then solve R·x=QT·b (2.10.4 ) by backsubstitution. Since QRdecomposition involves about twice as many operations as LUdecomposition, it is not used for typical systems of linear equations. However, we will meet special cases where QRis the method of choice. The standard algorithm for the QRdecomposition involves successive Householder transformations (to be discussed later in §11.2). We write a Householder matrix in the form 1−u⊗u/cwhere c=1 2u·u. Anappropriate Householder matrix applied to agiven matrix can zero all elements in a column of the matrix situated below a chosen element. Thus wearrangeforthefirstHouseholder matrix Q 1tozeroallelementsinthefirstcolumnof Abelow the first element. Similarly Q2zeroes all elements in the second column below the second element, and so on up to Qn−1. Thus R=Qn−1···Q1·A (2.10.5 ) Since the Householder matrices are orthogonal, Q=(Qn−1···Q1)−1=Q1···Qn−1 (2.10.6 ) In most applications we don’t need to form Qexplicitly; we instead store it in the factored form (2.10.6). Pivoting is not usually necessary unless the matrix Ais very close to singular. Ageneral QRalgorithmforrectangularmatricesincludingpivotingisgivenin [1]. Forsquare matrices, an implementation is the following: SUBROUTINE qrdcmp(a,n,np,c,d,sing) INTEGER n,npREAL a(np,np),c(n),d(n) LOGICAL sing Constructs the QR decomposition of a(1:n,1:n) , with physical dimension np. The upper triangular matrix Ris returned in the upper triangle of a, except for the diagonal elements ofRwhich are returned in d(1:n) . The orthogonal matrix Qis represented as a product of n−1Householder matrices Q1...Qn−1,w h e r eQj=1−uj⊗uj/cj.T h e ith component ofujis zero for i=1,...,j −1while the nonzero components are returned in a(i,j) for i=j ,...,n .sing returns as true if singularity is encountered during the decomposition, but the decomposition is still completed in this case. INTEGER i,j,k REAL scale,sigma,sum,tausing=.false.do 17k=1,n-1 scale=0. do11i=k,n scale=max(scale,abs(a(i,k))) enddo 11 if(scale.eq.0.)then Singular case. sing=.true.c(k)=0. d(k)=0. else FormQ kandQk·A. do12i=k,n a(i,k)=a(i,k)/scale enddo 12 sum=0.do 13i=k,n sum=sum+a(i,k)**2 enddo 13 sigma=sign(sqrt(sum),a(k,k)) a(k,k)=a(k,k)+sigma 2.10QR Decomposition 93Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).c(k)=sigma*a(k,k) d(k)=-scale*sigma do16j=k+1,n sum=0.do 14i=k,n sum=sum+a(i,k)*a(i,j) enddo 14 tau=sum/c(k)do 15i=k,n a(i,j)=a(i,j)-tau*a(i,k) enddo 15 enddo 16 endif enddo 17 d(n)=a(n,n) if(d(n).eq.0.)sing=.true. return END Thenextroutine, qrsolv,isused tosolvelinearsystems. Inmanyapplications onlythe part (2.10.4) of the algorithm is needed, so we separate it off into its own routine rsolv. SUBROUTINE qrsolv(a,n,np,c,d,b) INTEGER n,np REAL a(np,np),b(n),c(n),d(n) C USES rsolv Solves the set of nlinear equations A·x=b,w h e r e ais a matrix with physical dimension np. a,c,a n d dare input as the output of the routine qrdcmp and are not modified. b(1:n) is input as the right-hand side vector, and is overwritten with the solution vector on output. INTEGER i,jREAL sum,tau do 13j=1,n-1 FormQT·b. sum=0.do 11i=j,n sum=sum+a(i,j)*b(i) enddo 11 tau=sum/c(j)do 12i=j,n b(i)=b(i)-tau*a(i,j) enddo 12 enddo 13 call rsolv(a,n,np,d,b) SolveR·x=QT·b. return END SUBROUTINE rsolv(a,n,np,d,b) INTEGER n,np REAL a(np,np),b(n),d(n) Solves the set of nlinear equations R·x=b,w h e r eRis an upper triangular matrix stored inaandd.aanddare input as the output of the routine qrdcmp and are not modified. b(1:n) is input as the right-hand side vector, and is overwritten with the solution vector on output. INTEGER i,j REAL sum b(n)=b(n)/d(n)do 12i=n-1,1,-1 sum=0. do11j=i+1,n sum=sum+a(i,j)*b(j) enddo 11 b(i)=(b(i)-sum)/d(i) 94 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).enddo 12 return END See[2]for details on how to use QRdecomposition for constructing orthogonal bases, and for solving least-squares problems. (We prefer to use SVD, §2.6, for these purposes, because of its greater diagnostic capability in pathological cases.) Updatinga QR decomposition Somenumericalalgorithmsinvolvesolvingasuccessionoflinearsystemseachofwhich differs only slightly from its predecessor. Instead of doing O(N3)operations each time to solve the equations from scratch, one can often update a matrix factorization in O(N2) operations and use the new factorization to solve the next set of linear equations. The LU decomposition is complicated to update because of pivoting. However, QRturns out to be quite simple for a very common kind of update, A→A+s⊗t (2.10.7 ) (compare equation 2.7.1). In practice it is more convenient to work withthe equivalent form A=Q·R→A/prime=Q/prime·R/prime=Q·(R+u⊗v)( 2.10.8 ) One can go back and forth between equations (2.10.7) and (2.10.8) using the fact that Q is orthogonal, giving t=vand either s=Q·uoru=QT·s (2.10.9 ) Thealgorithm [2]hastwophases. Inthefirstweapply N−1Jacobirotations( §11.1)to reduceR+u⊗vto upper Hessenberg form. Another N−1Jacobi rotations transform this upper Hessenberg matrix to the new upper triangular matrix R/prime. The matrix Q/primeis simply the product of Qwith the 2(N−1)Jacobi rotations. In applications we usually want QT, and the algorithm can easily be rearranged to work with this matrix instead of with Q. SUBROUTINE qrupdt(r,qt,n,np,u,v) INTEGER n,npREAL r(np,np),qt(np,np),u(np),v(np) C USES rotate Given the QR decomposition of some n×nmatrix, calculates the QR decomposition of the matrix Q·(R+u⊗v). The matrices randqthave physical dimension np.N o t e t h a t QTis input and returned in qt. INTEGER i,j,kdo 11k=n,1,-1 Find largest ksuch that u(k)/negationslash=0. if(u(k).ne.0.)goto 1 enddo 11 k=1 1d o 12i=k-1,1,-1 Transform R+u⊗vto upper Hes- senberg. call rotate(r,qt,n,np,i,u(i),-u(i+1)) if(u(i).eq.0.)then u(i)=abs(u(i+1)) else if(abs(u(i)).gt.abs(u(i+1)))then u(i)=abs(u(i))*sqrt(1.+(u(i+1)/u(i))**2) else u(i)=abs(u(i+1))*sqrt(1.+(u(i)/u(i+1))**2) endif enddo 12 do13j=1,n r(1,j)=r(1,j)+u(1)*v(j) enddo 13 do14i=1,k-1 Transform upper Hessenberg matrix to upper triangular. call rotate(r,qt,n,np,i,r(i,i),-r(i+1,i)) enddo 14 2.11IsMatrix Inversionan N3Process? 95Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).return END SUBROUTINE rotate(r,qt,n,np,i,a,b) INTEGER n,np,iREAL a,b,r(np,np),qt(np,np) Given n×nmatrices randqtof physical dimension np, carry out a Jacobi rotation on rows i andi+1 of each matrix. aandbare the parameters of the rotation: cosθ=a/√ a2+b2, sinθ=b/√ a2+b2. INTEGER j REAL c,fact,s,w,y if(a.eq.0.)then Avoid unnecessary overflow or underflow. c=0.s=sign(1.,b) else if(abs(a).gt.abs(b))then fact=b/ac=sign(1./sqrt(1.+fact**2),a) s=fact*c else fact=a/bs=sign(1./sqrt(1.+fact**2),b) c=fact*s endifdo 11j=i,n Premultiply rby Jacobi rotation. y=r(i,j) w=r(i+1,j)r(i,j)=c*y-s*wr(i+1,j)=s*y+c*w enddo 11 do12j=1,n Premultiply qtby Jacobi rotation. y=qt(i,j) w=qt(i+1,j) qt(i,j)=c*y-s*wqt(i+1,j)=s*y+c*w enddo 12 returnEND We will make use of QRdecomposition, and its updating, in §9.7. CITED REFERENCES AND FURTHER READING: Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag), Chapter I/8. [1] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §§5.2, 5.3, 12.6. [2] 2.11 Is Matrix Inversion an N3Process? We close this chapter with a little entertainment, a bit of algorithmicprestidig- itation which probes more deeply into the subject of matrix inversion. We start with a seemingly simple question: