c11-2 (Numerical Recipes in C)

11.2 Reduction of a Symmetric Matrix to Tridiagonal Form469}}CITED REFERENCES AND FURTHER READING:Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), §8.4.Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide, 2nd ed., vol. 6 ofLecture Notes in Computer Science (New York: Springer-Verlag). [1]Wilkinson, J.H., and Reinsch, C.

1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). [2]11.2 Reduction of a Symmetric Matrixto Tridiagonal Form: Givens andHouseholder ReductionsAs already mentioned, the optimum strategy for finding eigenvalues andeigenvectors is, first, to reduce the matrix to a simple form, only then beginning aniterative procedure. For symmetric matrices, the preferred simple form is tridiagonal.The Givens reduction is a modification of the Jacobi method. Instead of trying toreduce the matrix all the way to diagonal form, we are content to stop when thematrix is tridiagonal.

This allows the procedure to be carried out in a finite numberof steps, unlike the Jacobi method, which requires iteration to convergence.Givens MethodFor the Givens method, we choose the rotation angle in equation (11.1.1) soas to zero an element that is not at one of the four “corners,” i.e., not app , apq ,or aqq in equation (11.1.3). Specifically, we first choose P23 to annihilate a31(and, by symmetry, a13 ). Then we choose P24 to annihilate a41 . In general, wechoose the sequenceP23 , P24 , . . . , P2n ; P34 , .

. . , P3n ; . . . ; Pn−1,nwhere Pjk annihilates ak,j−1. The method works because elements such as a0rp anda0rq , with r 6= p r 6= q, are linear combinations of the old quantities arp and arq , bySample 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).for (j=i+1;j<=n;j++)if (d[j] >= p) p=d[k=j];if (k != i) {d[k]=d[i];d[i]=p;for (j=1;j<=n;j++) {p=v[j][i];v[j][i]=v[j][k];v[j][k]=p;}}470Chapter 11.EigensystemsHouseholder MethodThe Householder algorithm reduces an n × n symmetric matrix A to tridiagonalform by n − 2 orthogonal transformations.

Each transformation annihilates therequired part of a whole column and whole corresponding row. The basic ingredientis a Householder matrix P, which has the formP = 1 − 2w · wT(11.2.1)where w is a real vector with |w|2 = 1. (In the present notation, the outer or matrixproduct of two vectors, a and b is written a · bT , while the inner or scalar product ofthe vectors is written as aT · b.) The matrix P is orthogonal, becauseP2 = (1 − 2w · wT ) · (1 − 2w · wT )= 1 − 4w · wT + 4w · (wT · w) · wT(11.2.2)=1Therefore P = P−1 .

But PT = P, and so PT = P−1 , proving orthogonality.Rewrite P asP =1 −u · uTH(11.2.3)where the scalar H isH≡and u can now be any vector.column of A. Choose1 2|u|2(11.2.4)Suppose x is the vector composed of the firstu = x ∓ |x|e1(11.2.5)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).equation (11.1.4). Thus, if arp and arq have already been set to zero, they remainzero as the reduction proceeds.

Evidently, of order n2 /2 rotations are required,and the number of multiplications in a straightforward implementation is of order4n3 /3, not counting those for keeping track of the product of the transformationmatrices, required for the eigenvectors.The Householder method, to be discussed next, is just as stable as the Givensreduction and it is a factor of 2 more efficient, so the Givens method is not generallyused. Recent work (see [1]) has shown that the Givens reduction can be reformulatedto reduce the number of operations by a factor of 2, and also avoid the necessityof taking square roots. This appears to make the algorithm competitive with theHouseholder reduction. However, this “fast Givens” reduction has to be monitoredto avoid overflows, and the variables have to be periodically rescaled.

There doesnot seem to be any compelling reason to prefer the Givens reduction over theHouseholder method.47111.2 Reduction of a Symmetric Matrix to Tridiagonal Formwhere e1 is the unit vector [1, 0, . . ., 0]T , and the choice of signs will be madelater. Thenu· (x ∓ |x|e1 )T · xH2u · (|x|2 ∓ |x|x1)=x−2|x|2 ∓ 2|x|x1= x−uP·x=x−= ±|x|e1This shows that the Householder matrix P acts on a given vector x to zero all itselements except the first one.To reduce a symmetric matrix A to tridiagonal form, we choose the vector xfor the first Householder matrix to be the lower n − 1 elements of the first column.Then the lower n − 2 elements will be zeroed:100P1 · A =  ...0···0(n−1)0P1 a11 a  21  ·  a31   ..  .a13···irrelevanta1nan10=a12a11 k 0 .. .a12a13···a1nirrelevant(11.2.7)0Here we have written the matrices in partitioned form, with (n−1)P denoting aHouseholder matrix with dimensions (n − 1) × (n − 1).

The quantity k is simplyplus or minus the magnitude of the vector [a21 , . . . , an1]T .The complete orthogonal transformation is nowa11 kA0 = P · A · P =  0 . ..0We have used the fact that PT = P.k0···irrelevant0(11.2.8)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).(11.2.6)472Chapter 11.EigensystemsThe identity block in the upper left corner insures that the tridiagonalization achievedin the first step will not be spoiled by this one, while the (n − 2)-dimensionalHouseholder matrix (n−2)P2 creates one additional row and column of the tridiagonaloutput. Clearly, a sequence of n − 2 such transformations will reduce the matrixA to tridiagonal form.Instead of actually carrying out the matrix multiplications in P · A · P, wecompute a vectorp≡ThenA·uH(11.2.10)u · uT) = A − p · uTHA0 = P · A · P = A − p · uT − u · pT + 2Ku · uTA · P = A · (1 −where the scalar K is defined byK=If we writeuT · p2H(11.2.11)q ≡ p − Ku(11.2.12)then we haveA0 = A − q · uT − u · qT(11.2.13)This is the computationally useful formula.Following [2], the routine for Householder reduction given below actually startsin the nth column of A, not the first as in the explanation above.

In detail, theequations are as follows: At stage m (m = 1, 2, . . . , n − 2) the vector u has the formuT = [ai1 , ai2, . . . , ai,i−2, ai,i−1 ±√σ, 0, . . . , 0](11.2.14)Herei ≡ n − m + 1 = n, n − 1, . . . , 3(11.2.15)and the quantity σ (|x|2 in our earlier notation) isσ = (ai1 )2 + · · · + (ai,i−1 )2(11.2.16)We choose the sign of σ in (11.2.14) to be the same as the sign of ai,i−1 to lessenroundoff error.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.