c11-1 (Numerical Recipes in C)

46311.1 Jacobi Transformations of a Symmetric MatrixCITED REFERENCES AND FURTHER READING:Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),Chapter 6. [1]Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). [2]IMSL Math/Library Users Manual (IMSL Inc., 2500 CityWest Boulevard, Houston TX 77042).

[4]NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.),Chapter F02. [5]Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), §7.7. [6]Wilkinson, J.H. 1965, The Algebraic Eigenvalue Problem (New York: Oxford University Press). [7]Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 13.Horn, R.A., and Johnson, C.R.

1985, Matrix Analysis (Cambridge: Cambridge University Press).11.1 Jacobi Transformations of a SymmetricMatrixThe Jacobi method consists of a sequence of orthogonal similarity transformations of the form of equation (11.0.14). Each transformation (a Jacobi rotation) isjust a plane rotation designed to annihilate one of the off-diagonal matrix elements.Successive transformations undo previously set zeros, but the off-diagonal elementsnevertheless get smaller and smaller, until the matrix is diagonal to machine precision.

Accumulating the product of the transformations as you go gives the matrixof eigenvectors, equation (11.0.15), while the elements of the final diagonal matrixare the eigenvalues.The Jacobi method is absolutely foolproof for all real symmetric matrices.

Formatrices of order greater than about 10, say, the algorithm is slower, by a significantconstant factor, than the QR method we shall give in §11.3. However, the Jacobialgorithm is much simpler than the more efficient methods. We thus recommend itfor matrices of moderate order, where expense is not a major consideration.The basic Jacobi rotation Ppq is a matrix of the formPpq=1···c...··· s.1 ..−s··· c···(11.1.1)1Here all the diagonal elements are unity except for the two elements c in rows (andcolumns) p and q.

All off-diagonal elements are zero except the two elements s and−s. The numbers c and s are the cosine and sine of a rotation angle φ, so c2 + s2 = 1.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 machinereadable files (including this one) to any servercomputer, is strictly prohibited.

To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide, 2nd ed., vol. 6 ofLecture Notes in Computer Science (New York: Springer-Verlag). [3]464Chapter 11.EigensystemsA plane rotation such as (11.1.1) is used to transform the matrix A according toA0 = PTpq · A · Ppq(11.1.2)Multiplying out equation (11.1.2) and using the symmetry of A, we get the explicitformulasa0rp = carp − sarqa0rq = carq + sarpr 6= p, r 6= q(11.1.4)a0pp = c2 app + s2 aqq − 2scapq(11.1.5)a0qqa0pq= s app + c aqq + 2scapq(11.1.6)= (c − s )apq + sc(app − aqq )(11.1.7)2222The idea of the Jacobi method is to try to zero the off-diagonal elements by aseries of plane rotations.

Accordingly, to set a0pq = 0, equation (11.1.7) gives thefollowing expression for the rotation angle φθ ≡ cot 2φ ≡aqq − appc2 − s2=2sc2apq(11.1.8)If we let t ≡ s/c, the definition of θ can be rewrittent2 + 2tθ − 1 = 0(11.1.9)The smaller root of this equation corresponds to a rotation angle less than π/4in magnitude; this choice at each stage gives the most stable reduction. Using theform of the quadratic formula with the discriminant in the denominator, we canwrite this smaller root ast=sgn(θ)√|θ| + θ2 + 1(11.1.10)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Now, PTpq · A changes only rows p and q of A, while A · Ppq changes only columnsp and q.

Notice that the subscripts p and q do not denote components of Ppq , butrather label which kind of rotation the matrix is, i.e., which rows and columns itaffects. Thus the changed elements of A in (11.1.2) are only in the p and q rowsand columns indicated below:· · · a01p · · · a01q · · ·......  .. ....  0 ap1 · · · a0pp · · · a0pq · · · a0pn ...... (11.1.3)A0 =  ...... 000  a0 q1 · · · aqp · · · aqq · · · aqn  .......  ..... 00· · · anp · · · anq · · ·11.1 Jacobi Transformations of a Symmetric Matrix465If θ is so large that θ2 would overflow on the computer, we set t = 1/(2θ).

Itnow follows thatc= √1t2 + 1s = tc(11.1.11)(11.1.12)a0pq = 0(11.1.13)The idea in the remaining equations is to set the new quantity equal to the oldquantity plus a small correction. Thus we can use (11.1.7) and (11.1.13) to eliminateaqq from (11.1.5), givinga0pp = app − tapq(11.1.14)a0qq = aqq + tapqa0rp = arp − s(arq + τ arp )(11.1.15)(11.1.16)a0rq = arq + s(arp − τ arq )(11.1.17)Similarly,where τ (= tan φ/2) is defined byτ≡s1+c(11.1.18)One can see the convergence of the Jacobi method by considering the sum ofthe squares of the off-diagonal elementsS=X|ars |2(11.1.19)r6=sEquations (11.1.4)–(11.1.7) imply thatS 0 = S − 2|apq |2(11.1.20)(Since the transformation is orthogonal, the sum of the squares of the diagonalelements increases correspondingly by 2|apq |2 .) The sequence of S’s thus decreasesmonotonically.

Since the sequence is bounded below by zero, and since we canchoose apq to be whatever element we want, the sequence can be made to convergeto zero.Eventually one obtains a matrix D that is diagonal to machine precision. Thediagonal elements give the eigenvalues of the original matrix A, sinceD = VT · A · V(11.1.21)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).When we actually use equations (11.1.4)–(11.1.7) numerically, we rewrite themto minimize roundoff error.

Equation (11.1.7) is replaced by466Chapter 11.EigensystemswhereV = P1 · P2 · P3 · · ·(11.1.22)the Pi ’s being the successive Jacobi rotation matrices. The columns of V are theeigenvectors (since A · V = V · D). They can be computed by applying(11.1.23)at each stage of calculation, where initially V is the identity matrix. In detail,equation (11.1.23) is0vrs= vrs(s 6= p, s 6= q)0vrp= cvrp − svrq0vrq(11.1.24)= svrp + cvrqWe 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 inwhich the elements are to be annihilated.