c11-0 (Numerical Recipes in C), страница 2

. . λN )(11.0.10)This is a particular case of a similarity transform of the matrix A,A→Z−1 · A · Z(11.0.11)for some transformation matrix Z. Similarity transformations play a crucial rolein the computation of eigenvalues, because they leave the eigenvalues of a matrixunchanged. This is easily seen fromdet Z−1 · A · Z − λ1 = det Z−1 · (A − λ1) · Z= det |Z| det |A − λ1| det Z−1 (11.0.12)= det |A − λ1|Equation (11.0.10) shows that any matrix with complete eigenvectors (which includesall normal matrices and “most” random nonnormal ones) can be diagonalized by asimilarity transformation, that the columns of the transformation matrix that effectsthe diagonalization are the right eigenvectors, and that the rows of its inverse arethe left eigenvectors.For real, symmetric matrices, the eigenvectors are real and orthonormal, so thetransformation matrix is orthogonal. The similarity transformation is then also anorthogonal transformation of the formA→ZT · A · Z(11.0.13)While real nonsymmetric matrices can be diagonalized in their usual case of completeeigenvectors, the transformation matrix is not necessarily real.

It turns out, however,that a real similarity transformation can “almost” do the job. It can reduce thematrix down to a form with little two-by-two blocks along the diagonal, all otherelements zero. Each two-by-two block corresponds to a complex-conjugate pairof complex eigenvalues. We will see this idea exploited in some routines givenlater in the chapter.The “grand strategy” of virtually all modern eigensystem routines is to nudgethe matrix A towards diagonal form by a sequence of similarity transformations,A→P−11 · A · P1→−1−1P−13 · P2 · P1 · A · P1 · P2 · P3→−1P−12 · P1 · A · P1 · P2→etc.(11.0.14)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).Diagonalization of a Matrix460Chapter 11.EigensystemsIf we get all the way to diagonal form, then the eigenvectors are the columns ofthe accumulated transformationXR = P1 · P2 · P3 · . . .(11.0.15)A = FL · FRor equivalentlyF−1L · A = FR(11.0.16)If we now multiply back together the factors in the reverse order, and use the secondequation in (11.0.16) we getFR · FL = F−1L · A · FL(11.0.17)which we recognize as having effected a similarity transformation on A with thetransformation matrix being FL ! In §11.3 and §11.6 we will discuss the QR methodwhich exploits this idea.Factorization methods also do not converge exactly in a finite number oftransformations.

But the better ones do converge rapidly and reliably, and, whenfollowing an appropriate initial reduction by simple similarity transformations, theyare the methods of choice.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).Sometimes we do not want to go all the way to diagonal form.

For example, ifwe are interested only in eigenvalues, not eigenvectors, it is enough to transformthe matrix A to be triangular, with all elements below (or above) the diagonal zero.In this case the diagonal elements are already the eigenvalues, as you can see bymentally evaluating (11.0.2) using expansion by minors.There are two rather different sets of techniques for implementing the grandstrategy (11.0.14). It turns out that they work rather well in combination, so mostmodern eigensystem routines use both. The first set of techniques constructs individual Pi ’s as explicit “atomic” transformations designed to perform specific tasks, forexample zeroing a particular off-diagonal element (Jacobi transformation, §11.1), ora whole particular row or column (Householder transformation, §11.2; eliminationmethod, §11.5).

In general, a finite sequence of these simple transformations cannotcompletely diagonalize a matrix. There are then two choices: either use the finitesequence of transformations to go most of the way (e.g., to some special form liketridiagonal or Hessenberg, see §11.2 and §11.5 below) and follow up with the secondset of techniques about to be mentioned; or else iterate the finite sequence of simpletransformations over and over until the deviation of the matrix from diagonal isnegligibly small. This latter approach is conceptually simplest, so we will discussit in the next section; however, for N greater than ∼ 10, it is computationallyinefficient by a roughly constant factor ∼ 5.The second set of techniques, called factorization methods, is more subtle.Suppose that the matrix A can be factored into a left factor FL and a right factorFR .

Then11.0 Introduction461“Eigenpackages of Canned Eigenroutines”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).You have probably gathered by now that the solution of eigensystems is a fairlycomplicated business.

It is. It is one of the few subjects covered in this book forwhich we do not recommend that you avoid canned routines. On the contrary, thepurpose of this chapter is precisely to give you some appreciation of what is goingon inside such canned routines, so that you can make intelligent choices about usingthem, and intelligent diagnoses when something goes wrong.You will find that almost all canned routines in use nowadays trace their ancestryback to routines published in Wilkinson and Reinsch’s Handbook for AutomaticComputation, Vol.

II, Linear Algebra [2]. This excellent reference, containing papersby a number of authors, is the Bible of the field. A public-domain implementationof the Handbook routines in FORTRAN is the EISPACK set of programs [3]. Theroutines in this chapter are translations of either the Handbook or EISPACK routines,so understanding these will take you a lot of the way towards understanding thosecanonical packages.IMSL [4] and NAG [5] each provide proprietary implementations, in FORTRAN,of what are essentially the Handbook routines.A good “eigenpackage” will provide separate routines, or separate paths throughsequences of routines, for the following desired calculations:• all eigenvalues and no eigenvectors• all eigenvalues and some corresponding eigenvectors• all eigenvalues and all corresponding eigenvectorsThe purpose of these distinctions is to save compute time and storage; it is wastefulto calculate eigenvectors that you don’t need.

Often one is interested only inthe eigenvectors corresponding to the largest few eigenvalues, or largest few inmagnitude, or few that are negative. The method usually used to calculate “some”eigenvectors is typically more efficient than calculating all eigenvectors if you desirefewer than about a quarter of the eigenvectors.A good eigenpackage also provides separate paths for each of the abovecalculations for each of the following special forms of the matrix:• real, symmetric, tridiagonal• real, symmetric, banded (only a small number of sub- and superdiagonalsare nonzero)• real, symmetric• real, nonsymmetric• complex, Hermitian• complex, non-HermitianAgain, the purpose of these distinctions is to save time and storage by using the leastgeneral routine that will serve in any particular application.In this chapter, as a bare introduction, we give good routines for the followingpaths:• all eigenvalues and eigenvectors of a real, symmetric, tridiagonal matrix(§11.3)• all eigenvalues and eigenvectors of a real, symmetric, matrix (§11.1–§11.3)• all eigenvalues and eigenvectors of a complex, Hermitian matrix(§11.4)• all eigenvalues and no eigenvectors of a real, nonsymmetric matrix462Chapter 11.Eigensystems(§11.5–§11.6)We also discuss, in §11.7, how to obtain some eigenvectors of nonsymmetricmatrices by the method of inverse iteration.Generalized and Nonlinear Eigenvalue ProblemsA · x = λB · x(11.0.18)where A and B are both matrices.

Most such problems, where B is nonsingular,can be handled by the equivalent(B−1 · A) · x = λx(11.0.19)Often A and B are symmetric and B is positive definite. The matrix B−1 · A in(11.0.19) is not symmetric, but we can recover a symmetric eigenvalue problemby using the Cholesky decomposition B = L · LT of §2.9. Multiplying equation(11.0.18) by L−1 , we getwhereC · (LT · x) = λ(LT · x)(11.0.20)C = L−1 · A · (L−1 )T(11.0.21)The matrix C is symmetric and its eigenvalues are the same as those of the originalproblem (11.0.18); its eigenfunctions are LT · x.