f2-9

2.9CholeskyDecomposition 89Sample 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).compared to N2for Levinson’s method. These methods are too complicated to include here. Papers by Bunch [6]and de Hoog [7]will give entry to the literature. CITED REFERENCES AND FURTHER READING: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), Chapter 5 [also treats some other special forms]. Forsythe, G.E., and Moler, C.B. 1967, Computer Solution of Linear Algebraic Systems (Engle- wood Cliffs, NJ: Prentice-Hall), §19. [1] Westlake,J.R.1968, AHandbookofNumericalMatrixInversionandSolutionofLinearEquations (New York: Wiley). [2] von Mises, R. 1964, Mathematical Theory of Probability and Statistics (New York: Academic Press), pp. 394ff. [3] Levinson, N., Appendix B of N. Wiener, 1949, Extrapolation, Interpolation and Smoothing of Stationary Time Series (New York: Wiley). [4] Robinson,E.A.,andTreitel,S.1980, GeophysicalSignalAnalysis (EnglewoodCliffs,NJ:Prentice- Hall), pp. 163ff. [5] Bunch,J.R. 1985, SIAM JournalonScientific andStatistical Computing ,vol.6,pp.349–364.[6] de Hoog, F. 1987, Linear Algebra and Its Applications , vol. 88/89, pp. 123–138. [7] 2.9 Cholesky Decomposition If a square matrix Ahappens to be symmetric and positive definite, then it has a special, more efficient, triangular decomposition. Symmetric means that aij=ajifor i, j =1 ,...,N, whilepositive definite means that v·A·v>0for all vectors v (2.9.1 ) (In Chapter 11 we will see that positive definite has the equivalent interpretation that Ahas all positive eigenvalues.) While symmetric, positive definite matrices are rather special, theyoccur quite frequently in some applications, so their special factorization, called Cholesky decomposition ,isgoodtoknowabout. Whenyoucanuseit,Choleskydecompositionisabout a factor of two faster than alternative methods for solving linear equations. Instead of seeking arbitrary lower and upper triangular factors LandU, Cholesky decomposition constructs a lower triangular matrix Lwhose transpose L Tcan itself serve as the upper triangular part. In other words we replace equation (2.3.1) by L·LT=A (2.9.2 ) This factorization is sometimes referred to as “taking the square root” of the matrix A. The components of LTare of course related to those of Lby LT ij=Lji (2.9.3 ) Writing out equation (2.9.2) in components, one readily obtains the analogs of equations (2.3.12)–(2.3.13), Lii=/parenleftBigg aii−i−1/summationdisplay k=1L2 ik/parenrightBigg1/2 (2.9.4 ) and Lji=1 Lii/parenleftBigg aij−i−1/summationdisplay k=1LikLjk/parenrightBigg j=i+1 ,i+2 ,...,N (2.9.5 ) 90 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).If you apply equations (2.9.4) and (2.9.5) in the order i=1 ,2,...,N, you will see that the L’s that occur on the right-hand side are already determined by the time they are needed. Also, only components aijwith j≥iare referenced. (Since Ais symmetric, these have complete information.) It is convenient, then, to have the factor Loverwrite the subdiagonal (lower triangular but not including the diagonal) part of A, preserving the input uppertriangularvaluesof A. Onlyoneextravectoroflength Nisneededtostorethediagonal part ofL. The operations count is N3/6executions of the inner loop (consisting of one multiply and one subtract), with also Nsquare roots. As already mentioned, this is about a factor 2 better than LUdecomposition of A(where its symmetry would be ignored). A straightforward implementation is SUBROUTINE choldc(a,n,np,p) INTEGER n,npREAL a(np,np),p(n) Given a positive-definite symmetric matrix a(1:n,1:n) , with physical dimension np,t h i s routineconstructsitsCholeskydecomposition, A=L·LT. Oninput,onlytheuppertriangle ofaneedbegiven;itisnotmodified. TheCholeskyfactor Lisreturnedinthelowertriangle ofa, except for its diagonal elements which are returned in p(1:n). INTEGER i,j,kREAL sum do 13i=1,n do12j=i,n sum=a(i,j)do 11k=i-1,1,-1 sum=sum-a(i,k)*a(j,k) enddo 11 if(i.eq.j)then if(sum.le.0.)pause ’choldc failed’ a , with rounding errors, is not positivedefinite. p(i)=sqrt(sum) else a(j,i)=sum/p(i) endif enddo 12 enddo 13 returnEND You might at this point wonder about pivoting. The pleasant answer is that Cholesky decomposition isextremelystablenumerically, withoutanypivoting atall. Failureof choldc simply indicates that the matrix A(or, with roundoff error, another very nearby matrix) is not positive definite. In fact, choldcis an efficient way to test whethera symmetric matrix is positive definite. (In this application, you will want to replace the pausewith some less drastic signaling method.) Once your matrix is decomposed, the triangular factor can be used to solve a linear equation by backsubstitution. The straightforward implementation of this is SUBROUTINE cholsl(a,n,np,p,b,x) INTEGER n,npREAL a(np,np),b(n),p(n),x(n) Solves the set of nlinear equations A·x=b,w h e r e ais a positive-definite symmetric matrixwithphysicaldimension np.aandpareinputastheoutputoftheroutine choldc. Onlythelowertriangleof aisaccessed. b(1:n)isinputastheright-handsidevector. The solution vector is returned in x(1:n).a,n,np,a n d pare not modified and can be left in place for successive calls with different right-hand sides b.bis not modified unless you identify bandxin the calling sequence, which is allowed. INTEGER i,kREAL sum do 12i=1,n SolveL·y=b,s t o r i n g yinx. sum=b(i)do 11k=i-1,1,-1 sum=sum-a(i,k)*x(k) 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.