f2-0

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 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).Chapter 2. Solution of Linear Algebraic Equations 2.0 Introduction A set of linear algebraic equations looks like this: a11x1+a12x2+a13x3+··· +a1NxN=b1 a21x1+a22x2+a23x3+··· +a2NxN=b2 a31x1+a32x2+a33x3+··· +a3NxN=b3 ··· ··· aM1x1+aM2x2+aM3x3+··· +aMN xN=bM(2.0.1 ) Here the Nunknowns xj,j=1 ,2,...,Nare related by Mequations. The coefficients aijwith i=1 ,2,...,Mand j=1 ,2,...,Nare known numbers, as are theright-hand side quantities bi,i=1 ,2,...,M. Nonsingularversus Singular Sets ofEquations IfN=Mthen there are as many equations as unknowns, and there is a good chance of solving for a unique solution set of xj’s. Analytically, there can fail to be a unique solution if one or more of the Mequations is a linear combination of the others, a condition called row degeneracy , or if all equations contain certain variables only in exactly the same linear combination, called column degeneracy . (For square matrices, a row degeneracy implies a column degeneracy, and vice versa.) A set of equations that is degenerate is called singular. We will consider singular matrices in some detail in §2.6. Numerically, at least two additional things can go wrong: While not exact linear combinations of each other, some of the equations may be so close to linearly dependentthat roundofferrors in the machine render them linearly dependent at some stage in the solution process. Inthis case your numerical procedure will fail, and it can tell you that it has failed. 22 2.0Introduction 23Sample 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).Accumulated roundoff errors in the solution process can swamp the true solution. This problem particularly emerges if Nis too large. The numerical procedure does not fail algorithmically. However, it returns a set of x’s that are wrong, as can be discoveredby direct substitution back intotheoriginalequations. Thecloserasetofequationsistobeingsingular,the more likely this is to happen, since increasingly close cancellations will occur during the solution. In fact, the preceding item can be viewed as the special case wherethe loss of significanceis unfortunatelytotal. Much of the sophistication of complicated “linear equation-solvingpackages” is devoted to the detection and/or correction of these two pathologies. As you work with large linear sets of equations, you will develop a feeling for when such sophistication is needed. It is difficult to give any firm guidelines, since there is no such thing as a “typical” linear problem. But here is a rough idea: Linear sets withNas large as 20 or 50 can be routinely solved in single precision (32 bit floating representations) without resorting to sophisticated methods, ifthe equations are not close to singular. With double precision (60 or 64 bits), this number can readilybe extended to Nas large as several hundred, after which point the limiting factor is generally machine time, not accuracy. Even larger linear sets, Nin the thousands or greater, can be solved when the coefficients are sparse (that is, mostly zero), by methods that take advantage of the sparseness. We discuss this further in §2.7. At the other end of the spectrum, one seems just as often to encounter linear problems which, by their underlying nature, are close to singular. In this case, you mightneed to resort to sophisticated methods even for the case of N=1 0(though rarely for N=5). Singular value decomposition ( §2.6) is a technique that can sometimes turn singular problems into nonsingular ones, in which case additional sophistication becomes unnecessary. Matrices Equation (2.0.1) can be written in matrix form as A·x=b (2.0.2 ) Heretheraiseddotdenotesmatrixmultiplication, Aisthematrixofcoefficients,and bis the right-hand side written as a column vector, A= a11 a12 ... a 1N a21 a22 ... a 2N ··· aM1 aM2 ... a MN b= b1 b2 ··· bM  (2.0.3 ) By convention, the first index on an element aijdenotes its row, the second index its column. A computer will store the matrix Aas a two-dimensional array. However, computer memory is numbered sequentially by its address, and so is intrinsically one-dimensional. Therefore the two-dimensional array Awill, at the hardware level, either be stored by columns in the order a11,a21,...,a M1,a 12,a22,...,a M2, ..., a 1N,a2N,...a MN 24 Chapter2. SolutionofLinearAlgebraicEquationsSample 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).17 22 34 48 5×65 73 89 91 10×111 121 130 142 15×16× × 18× 19× 20×21× 22× 23× 24× 25×26× 27× 28× 29× 30×31× 32× 33× 34× 35×36× 37× 38× 39× 40×17 m mpnp n Figure 2.0.1. A matrix of logical dimension mbynis stored in an array of physical dimension mpbynp. Locationsmarkedby “x”containextraneousinformationwhichmaybeleftoverfromsomeprevioususeof thephysicalarray. Circlednumbersshowtheactual orderingofthearrayincomputermemory,notusually relevant to the programmer. Note, however, that the logical array does not occupy consecutive memory locations. Tolocate an (i,j)element correctly, a subroutine must be told mpandnp, not just iandj. or elsestored by rows in the order a11,a12,...,a 1N,a 21,a22,...,a 2N, ..., a M1,a M2,...a MN FORTRAN always stores by columns, and user programs are generally allowed to exploit this fact to their advantage. By contrast, C,Pascal, and other languages generally store by rows. Note one confusing point in the terminology,that a matrix which is stored by columns (as in FORTRAN) has itsrow(i.e.,first) index changing mostrapidlyasonegoeslinearlythroughmemory,theoppositeofacar ’sodometer! For most purposes you don ’tneedto know what the order of storage is, since you reference an element by its two-dimensional address: a34=a(3,4). It is, however, essential that you understand the difference between an array ’sphysical dimensions and itslogical dimensions . When you pass an array to a subroutine, you must, in general, tell the subroutine bothof these dimensions. The distinction betweenthem is this: It may happenthat youhavea 4×4matrixstored in an array dimensioned as 10×10. This occurs most frequently in practice when you have dimensioned to the largest expected value of N, but are at the moment considering a value of Nsmaller than that largest possible one. In the example posed, the 16 elementsofthematrixdonotoccupy16consecutivememorylocations. Rathertheyare spread out among the 100 dimensioned locations of the array as if the whole 10×10matrix were filled. Figure 2.0.1 shows an additional example. Ifyouhaveasubroutinetoinvertamatrix,itscallmighttypicallylooklikethis: call matinv(a,ai,n,np) 2.0Introduction 25Sample 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).Here the subroutine has to be told both the logical size of the matrix that you want to invert (here n=4), and the physical size of the array in which it is stored (here np =1 0). Thisseemslikeatrivialpoint,andwearesorrytobelaborit. Butitturnsoutthat mostreportedfailuresofstandardlinearequationandmatrixmanipulationpackages are due to user errors in passing inappropriatelogical or physical dimensions! Tasks of ComputationalLinear Algebra We will consider the following tasks as falling in the general purview of this chapter: Solutionofthematrixequation A·x=bforanunknownvector x,whereA isasquarematrixofcoef ficients,raiseddotdenotesmatrixmultiplication, andbis a known right-hand side vector ( §2.1–§2.10). Solutionofmorethanonematrixequation A·xj=bj,forasetofvectors xj,j=1 ,2,...,eachcorrespondingtoadifferent,knownright-handside vectorbj. In this task the key simpli fication is that the matrix Ais held constant, while the right-handsides, the b’s, are changed( §2.1–§2.10). Calculationofthematrix A−1whichisthematrixinverseofasquarematrix A, i.e.,A·A−1=A−1·A=1, where1is the identity matrix (all zeros except for ones on the diagonal). This task is equivalent, for an N×N matrixA, to the previous task with Ndifferentbj’s(j=1 ,2,...,N ), namely the unit vectors ( bj=all zero elements except for 1 in the jth component). The corresponding x’s are then the columns of the matrix inverse of A(§2.1 and §2.3). Calculation of the determinant of a square matrix A(§2.3). IfM< N ,o ri f M =Nbut the equations are degenerate, then there are effectively fewer equations than unknowns. In this case there can be either no solution,orelsemorethanonesolutionvector x. Inthelatterevent,thesolutionspace consists of a particular solution xpadded to any linear combination of (typically) N−Mvectors (which are said to be in the nullspace of the matrix A). The task offinding the solution space of Ainvolves Singular value decomposition of a matrix A. This subject is treated in §2.6. In the opposite case there are more equations than unknowns, M> N. When this occurs there is, in general, no solution vector xto equation (2.0.1), and the set of equations is said to be overdetermined . It happens frequently, however, that the best“compromise ”solution is sought, the one that comes closest to satisfying all equations simultaneously. If closeness is de fined in the least-squares sense, i.e., that the sum of the squares of the differences between the left- and right-handsides ofequation(2.0.1)beminimized,thentheoverdeterminedlinearproblemreducestoa 26 Chapter2. Solutionof LinearAlgebraicEquationsSample 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).(usually) solvable linear problem, called the Linear least-squares problem. Thereducedsetofequationstobesolvedcanbewrittenasthe N×Nsetofequations (AT·A)·x=(AT·b)( 2.0.4 ) whereATdenotes the transpose of the matrix A. Equations (2.0.4) are called the normal equations of the linear least-squares problem. There is a close connection between singularvalue decompositionand the linear least-squares problem,and the latter is also discussed in §2.6. You should be warned that direct solution of the normalequations(2.0.4)isnotgenerallythebestwayto findleast-squaressolutions. Some other topics in this chapter include Iterative improvement of a solution ( §2.5) Various special forms: symmetric positive-de finite ( §2.9), tridiagonal (§2.4), band diagonal ( §2.4), Toeplitz ( §2.8), Vandermonde ( §2.8), sparse (§2.7) Strassen’s“fast matrix inversion ”(§2.11). Standard Subroutine Packages We cannot hope, in this chapter or in this book, to tell you everything there is to know about the tasks that have been de fined above. In many cases you will have no alternative but to use sophisticated black-box program packages. Several good onesareavailable. LINPACKwasdevelopedatArgonneNationalLaboratoriesand deserves particular mention because it is published, documented, and available forfreeuse. AsuccessortoLINPACK,LAPACK,isnowbecomingavailable. Packages available commercially include those in the IMSL and NAG libraries. Youshouldkeepinmindthatthesophisticatedpackagesaredesignedwithvery largelinear systems in mind. They thereforego to great effortto minimize not only the number of operations, but also the required storage. Routines for the varioustasks are usually provided in several versions, corresponding to several possible simplifications in the form of the input coef ficient matrix: symmetric, triangular, banded, positive de finite, etc. If you have a large matrix in one of these forms, you should certainly take advantage of the increased ef ficiency provided by these different routines, and not just use the form providedfor general matrices. There is also a great watershed dividing routines that are direct(i.e., execute in a predictable number of operations) from routines that are iterative(i.e., attempt to converge to the desired answer in however many steps are necessary). Iterativemethods become preferablewhen the battle against loss of signi ficance is in danger of being lost, either due to large Nor because the problem is close to singular. We will treat iterative methods only incompletely in this book, in §2.7 and in Chapters 18 and 19. These methods are important, but mostly beyond our scope. We will, however, discuss in detail a technique which is on the borderline between direct and iterative methods, namely the iterative improvementof a solution that has been obtained by direct methods ( §2.5). 2.1Gauss-JordanElimination 27Sample 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: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press). Gill,P.E., Murray, W.,andWright, M.H. 1991, NumericalLinearAlgebraandOptimization , vol. 1 (Redwood City, CA: Addison-Wesley). Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 4. Dongarra, J.J., et al. 1979, LINPACK User’s Guide (Philadelphia: S.I.A.M.). Coleman,T.F.,andVanLoan,C.1988, HandbookforMatrixComputations (Philadelphia:S.I.A.M.). Forsythe, G.E., and Moler, C.B. 1967, Computer Solution of Linear Algebraic Systems (Engle- wood Cliffs, NJ: Prentice-Hall). Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag). Westlake,J.R.1968, AHandbookofNumericalMatrixInversionandSolutionofLinearEquations (New York: Wiley). Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis , 2nd ed. (Reading, MA: Addison- Wesley), Chapter 2. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), Chapter 9. 2.1 Gauss-Jordan Elimination For inverting a matrix, Gauss-Jordan elimination is about as ef ficient as any other method. For solving sets of linear equations, Gauss-Jordan elimination produces boththe solution of the equations for one or more right-handside vectors b, and also the matrixinverse A−1. However,its principalweaknesses are (i) that it requires all the right-hand sides to be stored and manipulated at the same time, and (ii) that when the inverse matrix is notdesired, Gauss-Jordan is three times slower thanthebestalternativetechniqueforsolvingasinglelinearset( §2.3). Themethod ’s principal strength is that it is as stable as any other direct method, perhaps even a bit more stable when full pivoting is used (see below). If you come along later with an additional right-hand side vector, you can multiplyitbytheinversematrix,ofcourse. Thisdoesgiveananswer,butonethatis quite susceptible to roundofferror,not nearlyas goodas if the new vector had beenincluded with the set of right-hand side vectors in the first instance. Forthesereasons,Gauss-Jordaneliminationshouldusuallynotbeyourmethod offirst choice, either for solving linear equations or for matrix inversion. The decompositionmethodsin §2.3arebetter. WhydowegiveyouGauss-Jordanatall? Because it is straightforward, understandable, solid as a rock, and an exceptionallygood“psychological ”backupforthosetimesthatsomethingisgoingwrongandyou think itmightbe your linear-equation solver. Some people believe that the backup is more than psychological, that Gauss- Jordan elimination is an “independent ”numerical method. This turns out to be mostly myth. Except for the relatively minor differences in pivoting, described below, the actual sequence of operations performed in Gauss-Jordan elimination is very closely related to that performedby the routines in the next two sections.