f11-1

456 Chapter11. EigensystemsSample 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).CITED REFERENCES AND FURTHER READING: Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 6. [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag). [2] Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide , 2nd ed., vol. 6 of Lecture Notes in Computer Science (New York: Springer-Verlag). [3] IMSL Math/Library UsersManual (IMSL Inc., 2500CityWest Boulevard, HoustonTX77042).[4] NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.), Chapter F02. [5] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §7.7. [6] Wilkinson,J.H.1965, TheAlgebraicEigenvalueProblem (NewYork:OxfordUniversityPress).[7] Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 13. Horn,R.A.,andJohnson,C.R.1985, MatrixAnalysis (Cambridge:CambridgeUniversityPress). 11.1 Jacobi Transformations of a Symmetric Matrix The Jacobi method consists of a sequenceof orthogonalsimilarity transforma- tions of the form of equation (11.0.14). Each transformation (a Jacobi rotation )i s just a plane rotation designed to annihilate one of the off-diagonalmatrix elements. Successive transformationsundopreviouslyset zeros, but the off-diagonalelements nevertheless get smaller and smaller, until the matrix is diagonal to machine preci-sion. Accumulating the product of the transformations as you go gives the matrix of eigenvectors, equation (11.0.15), while the elements of the final diagonal matrix are the eigenvalues. The Jacobi methodis absolutelyfoolproofforall real symmetricmatrices. For matricesofordergreaterthanabout10,say,the algorithmis slower,bya significant constant factor, than the QRmethod we shall give in §11.3. However, the Jacobi algorithm is much simpler than the more efficient methods. We thus recommend it for matrices of moderate order, where expense is not a major consideration. The basic Jacobi rotation P pqis a matrix of the form Ppq= 1 ··· c··· s... 1... −s··· c ··· 1 (11.1.1 ) Here all the diagonal elements are unity except for the two elements cin rows (and columns) pandq. All off-diagonalelementsare zeroexceptthe two elements sand −s. Thenumbers candsarethecosineandsineofarotationangle φ,soc 2+s2=1. 11.1JacobiTransformationsofa SymmetricMatrix 457Sample 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).Aplanerotationsuchas(11.1.1)isusedtotransformthematrix Aaccordingto A/prime=PT pq·A·Ppq (11.1.2 ) Now,PT pq·Achanges only rows pandqofA, whileA·Ppqchanges only columns pandq. Notice that the subscripts pandqdo not denote components of Ppq,b u t rather label which kind of rotation the matrix is, i.e., which rows and columns itaffects. Thus the changed elements of Ain (11.1.2) are only in the pandqrows and columns indicated below: A /prime= ··· a /prime 1p··· a/prime 1q··· ............ a /prime p1··· a/prime pp··· a/prime pq··· a/prime pn ............ a /prime q1··· a/prime qp··· a/prime qq··· a/prime qn ............ ··· a/prime np··· a/prime nq··· (11.1.3 ) Multiplying out equation (11.1.2)and using the symmetry of A, we get the explicit formulas a /prime rp=carp−sa rq a/prime rq=carq+sa rpr/negationslash=p, r/negationslash=q (11.1.4) a/prime pp=c2app+s2aqq−2sca pq (11.1.5 ) a/prime qq=s2app+c2aqq+2sca pq (11.1.6 ) a/prime pq=(c2−s2)apq+sc(app−aqq)( 11.1.7 ) The idea of the Jacobi method is to try to zero the off-diagonal elements by a series of plane rotations. Accordingly, to set a/prime pq=0, equation (11.1.7) gives the following expression for the rotation angle φ θ≡cot 2φ≡c2−s2 2sc=aqq−app 2apq(11.1.8 ) If we let t≡s/c, the definition of θcan be rewritten t2+2tθ−1=0 ( 11.1.9 ) The smaller root of this equation correspondsto a rotation angle less than π/4 in magnitude; this choice at each stage gives the most stable reduction. Using the form of the quadratic formula with the discriminant in the denominator, we can write this smaller root as t=sgn(θ) |θ|+√ θ2+1(11.1.10 ) 458 Chapter11. EigensystemsSample 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θis so large that θ2would overflow on the computer, we set t=1/(2θ).I t now follows that c=1√ t2+1(11.1.11 ) s=tc (11.1.12 ) Whenweactuallyuseequations(11.1.4)–(11.1.7)numerically,werewritethem to minimize roundoff error. Equation (11.1.7) is replaced by a/prime pq=0 ( 11.1.13 ) The idea in the remaining equations is to set the new quantity equal to the old quantityplusasmallcorrection. Thuswecanuse(11.1.7)and(11.1.13)toeliminatea qqfrom (11.1.5), giving a/prime pp=app−tapq (11.1.14 ) Similarly, a/prime qq=aqq+tapq (11.1.15 ) a/prime rp=arp−s(arq+τa rp)( 11.1.16 ) a/prime rq=arq+s(arp−τa rq)( 11.1.17 ) where τ(= tan φ/2)is defined by τ≡s 1+c(11.1.18 ) One can see the convergence of the Jacobi method by considering the sum of the squares of the off-diagonal elements S=/summationdisplay r/negationslash=s|ars|2(11.1.19 ) Equations (11.1.4)–(11.1.7) imply that S/prime=S−2|apq|2(11.1.20 ) (Since the transformation is orthogonal, the sum of the squares of the diagonal elementsincreasescorrespondinglyby 2|apq|2.) Thesequenceof S’sthusdecreases monotonically. Since the sequence is bounded below by zero, and since we can choose apqto be whateverelement we want, the sequence can be made to converge to zero. Eventually one obtains a matrix Dthat is diagonal to machine precision. The diagonal elements give the eigenvalues of the original matrix A, since D=VT·A·V (11.1.21 ) 11.1JacobiTransformationsofa SymmetricMatrix 459Sample 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).where V=P1·P2·P3··· (11.1.22 ) thePi’s being the successive Jacobi rotation matrices. The columns of Vare the eigenvectors (since A·V=V·D). They can be computed by applying V/prime=V·Pi (11.1.23 ) at each stage of calculation, where initially Vis the identity matrix. In detail, equation (11.1.23) is v/prime rs=vrs (s/negationslash=p, s/negationslash=q) v/prime rp=cvrp−svrq v/prime rq=svrp+cvrq(11.1.24 ) We rewrite these equations in terms of τas in equations (11.1.16) and (11.1.17) to minimize roundoff. The only remaining question is the strategy one should adopt for the order in whichtheelementsaretobeannihilated. Jacobi’soriginalalgorithmof1846searched thewholeuppertriangleateachstageandsetthelargestoff-diagonalelementtozero.Thisis a reasonablestrategyforhandcalculation,butit is prohibitiveona computer sincethesearchalonemakeseachJacobirotationaprocessoforder N 2insteadof N. A better strategy for our purposes is the cyclic Jacobi method , where one annihilates elements in strict order. For example, one can simply proceed down the rows: P12,P13, ...,P1n; thenP23,P24, etc. One can show that convergence is generally quadratic for both the original or the cyclic Jacobi methods, for nondegenerate eigenvalues. One such set of n(n−1)/2Jacobi rotations is called asweep. The program below, based on the implementations in [1,2], uses two further refinements: •In the first three sweeps, we carry out the pqrotation only if |apq|>/epsilon1 for some threshold value /epsilon1=1 5S0 n2(11.1.25 ) where S0is the sum of the off-diagonal moduli, S0=/summationdisplay r<s|ars| (11.1.26 ) •After four sweeps, if |apq|/lessmuch| app|and|apq|/lessmuch| aqq|, we set |apq|=0 and skip the rotation. The criterion used in the comparison is |apq|< 10−(D+2)|app|,where Disthenumberofsignificantdecimaldigitsonthe machine, and similarly for |aqq|. 460 Chapter11. EigensystemsSample 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).In the following routine the n×nsymmetric matrix a(1:n,1:n) is stored in annp×nparray. On output, the superdiagonal elements of aare destroyed, but the diagonal and subdiagonal are unchanged and give full information on the original symmetric matrix a. The parameter dis a vectorof length np. On output, it returns the eigenvalues of ain its first nelements. During the computation, it contains the current diagonal of a. The matrix voutputs the normalized eigenvector belonging tod(k)in its kth column. The parameter nrotis the number of Jacobi rotations that were needed to achieve convergence. Typical matrices require 6 to 10 sweeps to achieve convergence,or 3n2to5n2 Jacobi rotations. Each rotation requires of order 4noperations, each consisting of a multiply and an add, so the total labor is of order 12n3to20n3operations. Calculation of the eigenvectors as well as the eigenvalues changes the operation count from 4nto6nper rotation, which is only a 50 percent overhead. SUBROUTINE jacobi(a,n,np,d,v,nrot) INTEGER n,np,nrot,NMAXREAL a(np,np),d(np),v(np,np) PARAMETER (NMAX=500) Computes all eigenvalues and eigenvectors of a real symmetric matrix a,w h i c hi so fs i ze n byn,storedinaphysical npbynparray. Onoutput, elements of aabovethediagonalare destroyed. dreturnstheeigenvaluesof ainitsfirst nelements. visamatrixwiththesame logical and physical dimensions as a, whose columns contain, on output, the normalized eigenvectors of a.nrotreturns the number of Jacobi rotations that were required. INTEGER i,ip,iq,j REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX) do12ip=1,n Initialize to the identity matrix. do11iq=1,n v(ip,iq)=0. enddo 11 v(ip,ip)=1. enddo 12 do13ip=1,n b(ip)=a(ip,ip) Initialize banddto the diagonal of a. d(ip)=b(ip) z(ip)=0. This vector will accumulate terms of the form tapq as in equation (11.1.14). enddo 13 nrot=0 do24i=1,50 sm=0. do15ip=1,n-1 Sum off-diagonal elements. do14iq=ip+1,n sm=sm+abs(a(ip,iq)) enddo 14 enddo 15 if(sm.eq.0.)return Thenormalreturn,whichreliesonquadraticconver- gence to machine underflow. if(i.lt.4)then tresh=0.2*sm/n**2 ...on the first three sweeps. else tresh=0. ...thereafter. endif do22ip=1,n-1 do21iq=ip+1,n g=100.*abs(a(ip,iq)) After four sweeps, skip the rotation ifthe off-diagonal element is small. if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip))) * .and.(abs(d(iq))+g.eq.abs(d(iq))))then a(ip,iq)=0. else if(abs(a(ip,iq)).gt.tresh)then h=d(iq)-d(ip) if(abs(h)+g.eq.abs(h))then 11.1JacobiTransformationsofa SymmetricMatrix 461Sample 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).t=a(ip,iq)/h t=1/(2θ) else theta=0.5*h/a(ip,iq) Equation (11.1.10). t=1./(abs(theta)+sqrt(1.+theta**2))if(theta.lt.0.)t=-t endif c=1./sqrt(1+t**2)s=t*ctau=s/(1.+c) h=t*a(ip,iq) z(ip)=z(ip)-hz(iq)=z(iq)+hd(ip)=d(ip)-h d(iq)=d(iq)+h a(ip,iq)=0.do 16j=1,ip-1 Case of rotations 1≤j<p. g=a(j,ip) h=a(j,iq)a(j,ip)=g-s*(h+g*tau)a(j,iq)=h+s*(g-h*tau) enddo 16 do17j=ip+1,iq-1 Case of rotations p<j<q. g=a(ip,j) h=a(j,iq) a(ip,j)=g-s*(h+g*tau)a(j,iq)=h+s*(g-h*tau) enddo 17 do18j=iq+1,n Case of rotations q<j ≤n. g=a(ip,j)h=a(iq,j) a(ip,j)=g-s*(h+g*tau) a(iq,j)=h+s*(g-h*tau) enddo 18 do19j=1,n g=v(j,ip) h=v(j,iq)v(j,ip)=g-s*(h+g*tau)v(j,iq)=h+s*(g-h*tau) enddo 19 nrot=nrot+1 endif enddo 21 enddo 22 do23ip=1,n b(ip)=b(ip)+z(ip) d(ip)=b(ip) Update dwith the sum of tapq, z(ip)=0. and reinitialize z. enddo 23 enddo 24 pause ’too many iterations in jacobi’returnEND Note that the above routine assumes that underflows are set to zero. On machines where this is not true, the program must be modified. The eigenvalues are not ordered on output. If sorting is desired, the following routine can be invoked to reorder the output of jacobior of later routines in this chapter. (The method, straight insertion, is N2rather than NlogN; but since you have just done an N3procedure to get the eigenvalues, you can afford yourself this little indulgence.) 462 Chapter11. EigensystemsSample 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).SUBROUTINE eigsrt(d,v,n,np) INTEGER n,np REAL d(np),v(np,np) Giventheeigenvalues dandeigenvectors vasoutputfrom jacobi(§11.1)or tqli(§11.3), this routine sorts the eigenvalues into descending order, and rearranges the columns of v correspondingly. The method is straight insertion. INTEGER i,j,kREAL pdo 13i=1,n-1 k=i p=d(i)do 11j=i+1,n if(d(j).ge.p)then k=j p=d(j) endif enddo 11 if(k.ne.i)then d(k)=d(i)d(i)=p do 12j=1,n p=v(j,i)v(j,i)=v(j,k) v(j,k)=p enddo 12 endif enddo 13 returnEND CITED REFERENCES AND FURTHER READING: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press),§8.4. Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide , 2nd ed., vol. 6 of Lecture Notes in Computer Science (New York: Springer-Verlag). [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag). [2] 11.2 Reduction of a Symmetric Matrix to Tridiagonal Form: Givens and Householder Reductions As already mentioned, the optimum strategy for finding eigenvalues and eigenvectors is, first, to reduce the matrix to a simple form, only then beginning an iterativeprocedure. Forsymmetricmatrices,thepreferredsimpleformistridiagonal. TheGivens reduction is a modification of the Jacobi method. Instead of trying to reduce the matrix all the way to diagonal form, we are content to stop when thematrix is tridiagonal. This allows the procedureto be carried out in a finite number of steps, unlike the Jacobi method, which requires iteration to convergence.