文档库

最新最全的文档下载
当前位置:文档库 > Cholesky Decomposition

Cholesky Decomposition

2.9Cholesky Decomposition89

Sample 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 copying of machine-readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.wendangku.net/doc/7a42dc186bd97f192279e94a.html or call 1-800-872-7423 (North America only),or send email to directcustserv@http://www.wendangku.net/doc/7a42dc186bd97f192279e94a.html (outside North America).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.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),Chapter5[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,A Handbook of Numerical Matrix Inversion and Solution of Linear Equations(New Y ork:Wiley).[2]von Mises,R.1964,Mathematical Theory of Probability and Statistics(New Y ork:Academic Press),pp.394ff.[3]Levinson,N.,Appendix B of N.Wiener,1949,Extrapolation,Interpolation and Smoothing of Stationary Time Series(New Y ork:Wiley).[4]Robinson,E.A.,and T reitel,S.1980,Geophysical Signal Analysis(Englewood Cliffs,NJ:Prentice-Hall),pp.163ff.[5]Bunch,J.R.1985,SIAM Journal on Scienti?c and Statistical Computing,vol.6,pp.349–364.[6]de Hoog,F.1987,Linear Algebra and Its Applications,vol.88/89,pp.123–138.[7]2.9Cholesky Decomposition If a square matrix A happens to be symmetric and positive de?nite,then it has a special,more ef?cient,triangular decomposition.Symmetric means that a ij=a ji for i,j=1,...,N,while positive de?nite means that v·A·v>0for all vectors v(2.9.1)(In Chapter11we will see that positive de?nite has the equivalent interpretation that A has all positive eigenvalues.)While symmetric,positive de?nite matrices are rather special,they occur quite frequently in some applications,so their special factorization,called Cholesky decomposition,is good to know about.When you can use it,Cholesky decomposition is about a factor of two faster than alternative methods for solving linear equations.Instead of seeking arbitrary lower and upper triangular factors L and U,Cholesky decomposition constructs a lower triangular matrix L whose transpose L T can itself serve as the upper triangular part.In other words we replace equation(2.3.1)by L·L T=A(2.9.2)This factorization is sometimes referred to as“taking the square root”of the matrix A.The components of L T are of course related to those of L by L T ij=L ji(2.9.3)Writing out equation(2.9.2)in components,one readily obtains the analogs of equations(2.3.12)–(2.3.13), L ii=

a ii?

i?1

k=1

L2ik

1/2

(2.9.4)

and

L ji=1

L ii

a ij?

i?1

k=1

L ik L jk

j=i+1,i+2,...,N(2.9.5)

90Chapter2.Solution of Linear Algebraic Equations Sample 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 copying of machine-readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.wendangku.net/doc/7a42dc186bd97f192279e94a.html or call 1-800-872-7423 (North America only),or send email to directcustserv@http://www.wendangku.net/doc/7a42dc186bd97f192279e94a.html (outside North America).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 a ij with j≥i are referenced.(Since A is symmetric,these have complete information.)It is convenient,then,to have the factor L overwrite the subdiagonal(lower triangular but not including the diagonal)part of A,preserving the input upper triangular values of A.Only one extra vector of length N is needed to store the diagonal part of L.The operations count is N3/6executions of the inner loop(consisting of one multiply and one subtract),with also N square roots.As already mentioned,this is about a factor2better than LU decomposition of A(where its symmetry would be ignored).A straightforward implementation is SUBROUTINE choldc(a,n,np,p)INTEGER n,np REAL a(np,np),p(n)Given a positive-de?nite symmetric matrix a(1:n,1:n),with physical dimension np,this routine constructs its Cholesky decomposition,A=L·L T.On input,only the upper triangle of a need be given;it is not modi?ed.The Cholesky factor L is returned in the lower triangle of a,except for its diagonal elements which are returned in p(1:n).INTEGER i,j,k REAL sum do13i=1,n do12j=i,n sum=a(i,j)do11k=i-1,1,-1sum=sum-a(i,k)*a(j,k)enddo11if(i.eq.j)then if(sum.le.0.)pause’choldc failed’a,with rounding errors,is not positive de?nite.p(i)=sqrt(sum)else a(j,i)=sum/p(i)endif enddo12enddo13return END You might at this point wonder about pivoting.The pleasant answer is that Cholesky decomposition is extremely stable numerically,without any pivoting at all.Failure of choldc simply indicates that the matrix A(or,with roundoff error,another very nearby matrix)is not positive de?nite.In fact,choldc is an ef?cient way to test whether a symmetric matrix is positive de?nite.(In this application,you will want to replace the pause with 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,np REAL a(np,np),b(n),p(n),x(n)Solves the set of n linear equations A·x=b,where a is a positive-de?nite symmetric matrix with physical dimension np.a and p are input as the output of the routine choldc.Only the lower triangle of a is accessed.b(1:n)is input as the right-hand side vector.The solution vector is returned in x(1:n).a,n,np,and p are not modi?ed and can be left

in place for successive calls with di?erent right-hand sides b.b is not modi?ed unless you identify b and x in the calling sequence,which is allowed.

INTEGER i,k

REAL sum

do12i=1,n Solve L·y=b,storing y in x.

sum=b(i)

do11k=i-1,1,-1

sum=sum-a(i,k)*x(k)

2.10QR Decomposition91 Sample 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 copying of machine-readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.wendangku.net/doc/7a42dc186bd97f192279e94a.html or call 1-800-872-7423 (North America only),or send email to directcustserv@http://www.wendangku.net/doc/7a42dc186bd97f192279e94a.html (outside North America).enddo11x(i)=sum/p(i)enddo12do14i=n,1,-1Solve L T·x=y.sum=x(i)do13k=i+1,n sum=sum-a(k,i)*x(k)enddo13x(i)=sum/p(i)enddo14return END A typical use of choldc and cholsl is in the inversion of covariance matrices describing the?t of data to a model;see,e.g.,§15.6.In this,and many other applications,one often needs L?1.The lower triangle of this matrix can be ef?ciently found from the output of choldc:do13i=1,n a(i,i)=1./p(i)do12j=i+1,n sum=0.do11k=i,j-1sum=sum-a(j,k)*a(k,i)enddo11a(j,i)=sum/p(j)enddo12enddo13CITED REFERENCES AND FURTHER READING:Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.II of Handbook for Automatic Com-putation(New Y ork:Springer-Verlag),Chapter I/1.Gill,P.E.,Murray,W.,and Wright,M.H.1991,Numerical Linear Algebra and Optimization,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.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),§

4.2.2.10QR Decomposition There is another matrix factorization that is sometimes very useful,the so-called QR decomposition,A=Q·R(2.10.1)

Here R is upper triangular,while Q is orthogonal,that is,

Q T·Q=1(2.10.2)

where Q T is the transpose matrix of Q.Although the decomposition exists for a general rectangular matrix,we shall restrict our treatment to the case when all the matrices are square, with dimensions N×N.