f2-8

82 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).SUBROUTINE asolve(n,b,x,itrnsp) INTEGER n,itrnsp,ija,NMAX,i DOUBLE PRECISION x(n),b(n),sa PARAMETER (NMAX=1000)COMMON /mat/ sa(NMAX),ija(NMAX) The matrix is stored somewhere. do 11i=1,n x(i)=b(i)/sa(i) The matrix /tildewideAis the diagonal part of A,s t o r e di n the first nelements of sa.S i n c et h et r a n s p o s e matrix has the same diagonal, the flag itrnsp is not used.enddo 11 return END CITED REFERENCES AND FURTHER READING: Tewarson, R.P. 1973, Sparse Matrices (New York: Academic Press). [1] Jacobs,D.A.H.(ed.)1977, TheStateoftheArtinNumericalAnalysis (London:AcademicPress), Chapter I.3 (by J.K. Reid). [2] George,A.,andLiu,J.W.H.1981, ComputerSolutionofLargeSparsePositiveDefiniteSystems (Englewood Cliffs, NJ: Prentice-Hall). [3] NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.). [4] IMSL Math/Library UsersManual (IMSL Inc., 2500CityWest Boulevard, HoustonTX77042).[5] Eisenstat,S.C.,Gursky,M.C.,Schultz,M.H.,andSherman,A.H.1977, YaleSparseMatrixPack- age,TechnicalReports112and114(YaleUniversityDepartmentofComputerScience).[6] Knuth,D.E.1968, FundamentalAlgorithms ,vol.1ofTheArtofComputerProgramming (Reading, MA: Addison-Wesley), §2.2.6. [7] Kincaid,D.R., Respess, J.R., Young, D.M., andGrimes, R.G. 1982, ACMTransactions onMath- ematical Software , vol. 8, pp. 302–322. [8] PCGPAK User’s Guide (New Haven: Scientific Computing Associates, Inc.). [9] Bentley, J. 1986, Programming Pearls (Reading, MA: Addison-Wesley), §9. [10] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), Chapters 4 and 10, particularly §§10.2–10.3. [11] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 8. [12] Baker, L. 1991, More C Tools for Scientists and Engineers (New York: McGraw-Hill). [13] Fletcher,R.1976,in NumericalAnalysisDundee1975 ,Lecture NotesinMathematics, vol.506, A. Dold and B Eckmann, eds. (Berlin: Springer-Verlag), pp. 73–89. [14] Saad, Y., and Schulz, M. 1986, SIAM Journal on Scientific and Statistical Computing , vol. 7, pp. 856–869. [15] Bunch, J.R., and Rose, D.J. (eds.) 1976, Sparse Matrix Computations (New York: Academic Press). Duff, I.S., and Stewart, G.W. (eds.) 1979, Sparse Matrix Proceedings 1978 (Philadelphia: S.I.A.M.). 2.8 Vandermonde Matrices and Toeplitz Matrices In§2.4 the case of a tridiagonal matrix was treated specially, because that particular type of linear system admits a solution in only of order Noperations, rather than of order N3for the general linear problem. When such particular types 2.8VandermondeMatricesandToeplitzMatrices 83Sample 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).exist, it is important to know about them. Your computational savings, should you ever happen to be working on a problem that involves the right kind of particulartype, can be enormous. This section treats two special types of matrices that can be solved in of order N 2operations, not as good as tridiagonal, but a lot better than the general case. (Other than the operations count, these two types having nothing in common.) Matrices of the first type, termed Vandermonde matrices , occur in some problems havingtodowith thefittingofpolynomials,thereconstructionofdistributionsfrom their moments, and also other contexts. In this book, for example, a Vandermonde problem crops up in §3.5. Matrices of the second type, termed Toeplitz matrices , tend to occur in problems involving deconvolution and signal processing. In this book, a Toeplitz problem is encountered in §13.7. These are not the onlyspecial types of matrices worth knowing about. The Hilbert matrices , whose components are of the form aij=1 /(i+j−1),i , j = 1,...,Ncan be inverted by an exact integer algorithm, and are very difficultto invertinanyotherway,sincetheyarenotoriouslyill-conditioned(see [1]fordetails). The Sherman-Morrison and Woodbury formulas, discussed in §2.7, can sometimes be used to convert new special forms into old ones. Reference [2]gives some other special forms. We have not found these additional forms to arise as frequently as the two that we now discuss. Vandermonde Matrices A Vandermonde matrix of size N×Nis completely determined by Narbitrary numbers x1,x2,...,x N, in terms of which its N2components are the integer powers xj−1 i,i , j =1,...,N. Evidently there are two possible such forms, depending on whether we view the i’s as rows, j’s as columns, or vice versa. In the former case, we get a linear system of equations that looks like this,  1x 1x2 1···xN−1 1 1x2x2 2···xN−1 2 ............ 1x Nx2 N···xN−1 N · c 1 c2 ... cN = y 1 y2 ... yN (2.8.1 ) Performing the matrix multiplication, you will see that this equation solves for the unknown coefficients c iwhich fit a polynomial to the Npairs of abscissas and ordinates (xj,yj). Precisely this problem will arise in §3.5, and the routine given there will solve (2.8.1) by the method that we are about to describe. The alternative identification of rows and columns leads to the set of equations  11 ··· 1 x1 x2··· xN x2 1 x22··· x2 N ··· xN−1 1 xN−1 2 ···xN−1 N · w1 w2 w3 ··· wN = q1 q2 q3 ··· qN (2.8.2 ) Write this out and you will see that it relates to the problem of moments : Given the values ofNpoints xi, find the unknown weights wi, assigned so as to match the given values qjof the first Nmoments. (For more on this problem, consult [3].) The routine given in this section solves (2.8.2). 84 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).The method of solution of both (2.8.1) and (2.8.2) is closely related to Lagrange’s polynomialinterpolationformula,whichwewillnotformallymeetuntil §3.1below. Notwith- standing, the following derivation should be comprehensible: LetPj(x)be the polynomial of degree N−1defined by Pj(x)=N/productdisplay n=1 (n /negationslash=j)x−xn xj−xn=N/summationdisplay k=1Ajkxk−1(2.8.3 ) Here the meaning of the last equality is to define the components of the matrix Aijas the coefficients that arise when the product is multiplied out and like terms collected. The polynomial Pj(x)is a function of xgenerally. But you will notice that it is specifically designed so that it takes on a value of zero at all xiwithi/negationslash=j, and has a value of unity at x=xj. In other words, Pj(xi)=δij=N/summationdisplay k=1Ajkxk−1 i (2.8.4 ) But(2.8.4)saysthat Ajkisexactlytheinverseofthematrixofcomponents xk−1 i,which appears in (2.8.2), with the subscript as the column index. Therefore the solution of (2.8.2)is just that matrix inverse times the right-hand side, w j=N/summationdisplay k=1Ajkqk (2.8.5 ) Asforthetransposeproblem(2.8.1),wecanusethefactthattheinverseofthetranspose is the transpose of the inverse, so cj=N/summationdisplay k=1Akjyk (2.8.6 ) The routine in §3.5 implements this. It remains to find a good way of multiplying out the monomial terms in (2.8.3), in order toget thecomponents of Ajk. Thisisessentially abookkeeping problem, and wewillletyou read the routine itself to see how it can be solved. One trick is to define a master P(x)by P(x)≡N/productdisplay n=1(x−xn)( 2.8.7 ) workoutitscoefficients,andthenobtainthenumeratorsanddenominatorsofthespecific Pj’s viasyntheticdivisionbytheonesupernumeraryterm. (See §5.3formoreonsyntheticdivision.) Since each such division is only a process of order N, the total procedure is of order N2. You should be warned that Vandermonde systems are notoriously ill-conditioned, by their very nature. (As an aside anticipating §5.8, the reason is the same as that which makes Chebyshev fitting so impressively accurate: there exist high-order polynomials that are verygood uniform fits to zero. Hence roundoff error can introduce rather substantial coefficients ofthe leading termsof these polynomials.) Itisagood ideaalways tocompute Vandermonde problems in double precision. The routine for (2.8.2) which follows is due to G.B. Rybicki. SUBROUTINE vander(x,w,q,n) INTEGER n,NMAXDOUBLE PRECISION q(n),w(n),x(n) PARAMETER (NMAX=100) Solves the Vandermonde linear system/summationtext N i=1xk−1 iwi=qk(k=1 ,...,N ). Input consists of the vectors x(1:n) andq(1:n) ; the vector w(1:n) is output. 2.8VandermondeMatricesandToeplitzMatrices 85Sample 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).Parameters: NMAX is the maximum expected value of n. INTEGER i,j,k DOUBLE PRECISION b,s,t,xx,c(NMAX) if(n.eq.1)then w(1)=q(1) else do11i=1,n Initialize array. c(i)=0.d0 enddo 11 c(n)=-x(1) Coefficients of the master polynomial are found by recur- sion. do13i=2,n xx=-x(i)do 12j=n+1-i,n-1 c(j)=c(j)+xx*c(j+1) enddo 12 c(n)=c(n)+xx enddo 13 do15i=1,n Each subfactor in turn xx=x(i)t=1.d0 b=1.d0 s=q(n)do 14k=n,2,-1 is synthetically divided, b=c(k)+xx*b s=s+q(k-1)*b matrix-multiplied by the right-hand side, t=xx*t+b enddo 14 w(i)=s/t and supplied with a denominator. enddo 15 endif return END ToeplitzMatrices AnN×NToeplitz matrix is specified by giving 2N−1numbers Rk,k =−N+ 1,...,−1,0,1,...,N −1. Those numbers are then emplaced as matrix elements constant along the (upper-left to lower-right) diagonals of the matrix:  R 0 R−1R−2···R−(N−2)R−(N−1) R1 R0 R−1···R−(N−3)R−(N−2) R2 R1 R0···R−(N−4)R−(N−3) ··· ··· RN−2RN−3RN−4··· R0 R−1 RN−1RN−2RN−3··· R1 R0 (2.8.8 ) The linear Toeplitz problem can thus be written as N/summationdisplay j=1Ri−jxj=yi (i=1,...,N )( 2.8.9 ) where the xj’s,j=1,...,N, are the unknowns to be solved for. The Toeplitz matrix is symmetric if Rk=R−kfor all k. Levinson [4]developed an algorithm for fast solution of the symmetric Toeplitz problem, by a bordering method , that is, a recursive procedure that solves the M-dimensional Toeplitz problem M/summationdisplay j=1Ri−jx(M) j =yi (i=1,...,M )( 2.8.10 ) 86 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).inturnfor M=1,2,...untilM=N,thedesiredresult,isfinallyreached. Thevector x(M) j is the result at the Mth stage, and becomes the desired answer only when Nis reached. Levinson’s method is well documented in standard texts (e.g., [5]). The useful fact that the method generalizes to the nonsymmetric case seems to be less well known. At some risk of excessive detail, we therefore give a derivation here, due to G.B. Rybicki. Infollowingarecursionfromstep Mtostep M+1wefindthatourdeveloping solution x(M)changes in this way: M/summationdisplay j=1Ri−jx(M) j =yi i=1,...,M (2.8.11 ) becomes M/summationdisplay j=1Ri−jx(M+1) j +Ri−(M+1)x(M+1) M+1 =yi i=1,...,M +1 ( 2.8.12 ) By eliminating yiwe find M/summationdisplay j=1Ri−j/parenleftBigg x(M) j−x(M+1) j x(M+1) M+1/parenrightBigg =Ri−(M+1) i=1,...,M (2.8.13 ) or by letting i→M+1−iandj→M+1−j, M/summationdisplay j=1Rj−iG(M) j =R−i (2.8.14 ) where G(M) j≡x(M) M+1−j−x(M+1) M+1−j x(M+1) M+1(2.8.15 ) To put this another way, x(M+1) M+1−j=x(M) M+1−j−x(M+1) M+1G(M) j j=1,...,M (2.8.16 ) Thus, if we can use recursion to find the order Mquantities x(M)andG(M)andthe single order M+1quantity x(M+1) M+1, then all of the other x(M+1) jwill follow. Fortunately, the quantity x(M+1) M+1follows from equation (2.8.12) with i=M+1, M/summationdisplay j=1RM+1−jx(M+1) j +R0x(M+1) M+1 =yM+1 (2.8.17 ) For the unknown order M +1quantities x(M+1) jwe can substitute the previous order quantities in Gsince G(M) M+1−j=x(M) j−x(M+1) j x(M+1) M+1(2.8.18 ) The result of this operation is x(M+1) M+1 =/summationtextM j=1RM+1−jx(M) j−yM+1 /summationtextM j=1RM+1−jG(M) M+1−j−R0(2.8.19 ) The only remaining problem is to develop a recursion relation for G. Before we do that, however, we should point out that there are actually two distinct sets of solutions to theoriginal linear problem for a nonsymmetric matrix, namely right-hand solutions (which we 2.8VandermondeMatricesandToeplitzMatrices 87Sample 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).have been discussing) and left-hand solutions zi. The formalism for the left-hand solutions differs only in that we deal with the equations M/summationdisplay j=1Rj−iz(M) j =yi i=1,...,M (2.8.20 ) Then, the same sequence of operations on this set leads to M/summationdisplay j=1Ri−jH(M) j =Ri (2.8.21 ) where H(M) j≡z(M) M+1−j−z(M+1) M+1−j z(M+1) M+1(2.8.22 ) (compare with 2.8.14 – 2.8.15). The reason for mentioning the left-hand solutions now is that, by equation (2.8.21), the Hjsatisfy exactly the same equation as the xjexcept for the substitution yi→Rion the right-hand side. Therefore we can quickly deduce from equation (2.8.19) that H(M+1) M+1 =/summationtextM j=1RM+1−jH(M) j−RM+1 /summationtextM j=1RM+1−jG(M) M+1−j−R0(2.8.23 ) By the same token, Gsatisfies the same equation as z, except for the substitution yi→R−i. This gives G(M+1) M+1 =/summationtextM j=1Rj−M−1G(M) j−R−M−1 /summationtextM j=1Rj−M−1H(M) M+1−j−R0(2.8.24 ) Thesame“morphism”alsoturnsequation(2.8.16),anditspartnerfor z,intothefinalequations G(M+1) j =G(M) j−G(M+1) M+1H(M) M+1−j H(M+1) j =H(M) j−H(M+1) M+1G(M) M+1−j(2.8.25 ) Now, starting with the initial values x(1) 1=y1/R 0 G(1) 1=R−1/R 0 H(1) 1 =R1/R 0 (2.8.26 ) we can recurse away. At each stage Mwe use equations (2.8.23) and (2.8.24) to find H(M+1) M+1,G(M+1) M+1,andthenequation(2.8.25)tofindtheothercomponentsof H(M+1),G(M+1). From there the vectors x(M+1)and/or z(M+1)are easily calculated. Theprogram belowdoesthis. Itincorporatesthesecond equation in(2.8.25)intheform H(M+1) M+1−j=H(M) M+1−j−H(M+1) M+1G(M) j (2.8.27 ) so that the computation can be done “in place.” Notice that the above algorithm fails if R0=0. In fact, because the bordering method does not allow pivoting, the algorithm will fail if any of the diagonal principal minors of theoriginal Toeplitz matrix vanish. (Compare with discussion of the tridiagonal algorithm in§2.4.) If the algorithm fails, your matrix is not necessarily singular — you might just have to solve your problem by a slower and more general algorithm such as LUdecomposition with pivoting. The routine that implements equations (2.8.23)–(2.8.27) is also due to Rybicki. Note that the routine’s r(n+j)is equal to R jabove, so that subscripts on the rarray vary from 1to2N−1. 88 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).SUBROUTINE toeplz(r,x,y,n) INTEGER n,NMAX REAL r(2*n-1),x(n),y(n) PARAMETER (NMAX=100) Solves the Toeplitz system/summationtextN j=1R(N+i−j)xj=yi(i=1 ,...,N ). The Toeplitz matrix need not be symmetric. yandrare input arrays of length nand2*n-1 , respectively. x is the output array, of length n. Parameter: NMAX is the maximum anticipated value of n. INTEGER j,k,m,m1,m2REAL pp,pt1,pt2,qq,qt1,qt2,sd,sgd,sgn,shn,sxn, * g(NMAX),h(NMAX) if(r(n).eq.0.) goto 99x(1)=y(1)/r(n) Initialize for the recursion. if(n.eq.1)return g(1)=r(n-1)/r(n)h(1)=r(n+1)/r(n)do 15m=1,n Main loop over the recursion. m1=m+1 sxn=-y(m1) Compute numerator and denominator for x, sd=-r(n) do11j=1,m sxn=sxn+r(n+m1-j)*x(j)sd=sd+r(n+m1-j)*g(m-j+1) enddo 11 if(sd.eq.0.)goto 99x(m1)=sxn/sd whence x. do 12j=1,m x(j)=x(j)-x(m1)*g(m-j+1) enddo 12 if(m1.eq.n)returnsgn=-r(n-m1) Compute numerator and denominator for Gand H, shn=-r(n+m1) sgd=-r(n)do 13j=1,m sgn=sgn+r(n+j-m1)*g(j) shn=shn+r(n+m1-j)*h(j)sgd=sgd+r(n+j-m1)*h(m-j+1) enddo 13 if(sd.eq.0..or.sgd.eq.0.)goto 99g(m1)=sgn/sgd whence Gand H. h(m1)=shn/sdk=m m2=(m+1)/2 pp=g(m1)qq=h(m1) do 14j=1,m2 pt1=g(j)pt2=g(k)qt1=h(j) qt2=h(k) g(j)=pt1-pp*qt2g(k)=pt2-pp*qt1 h(j)=qt1-qq*pt2 h(k)=qt2-qq*pt1k=k-1 enddo 14 enddo 15 Back for another recurrence. pause ’never get here in toeplz’ 99 pause ’singular principal minor in toeplz’ END Ifyouareinthebusinessofsolving verylargeToeplitzsystems,youshouldfindoutabout so-called “new, fast” algorithms, which require only on the order of N(logN)2operations, 2.9CholeskyDecomposition 89Sample 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).compared to N2for Levinson’s method. These methods are too complicated to include here. Papers by Bunch [6]and de Hoog [7]will give entry to the literature. CITED REFERENCES AND FURTHER READING: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), Chapter 5 [also treats some other special forms]. Forsythe, G.E., and Moler, C.B. 1967, Computer Solution of Linear Algebraic Systems (Engle- wood Cliffs, NJ: Prentice-Hall), §19. [1] Westlake,J.R.1968, AHandbookofNumericalMatrixInversionandSolutionofLinearEquations (New York: Wiley). [2] von Mises, R. 1964, Mathematical Theory of Probability and Statistics (New York: Academic Press), pp. 394ff. [3] Levinson, N., Appendix B of N. Wiener, 1949, Extrapolation, Interpolation and Smoothing of Stationary Time Series (New York: Wiley). [4] Robinson,E.A.,andTreitel,S.1980, GeophysicalSignalAnalysis (EnglewoodCliffs,NJ:Prentice- Hall), pp. 163ff. [5] Bunch,J.R. 1985, SIAM JournalonScientific andStatistical Computing ,vol.6,pp.349–364.[6] de Hoog, F. 1987, Linear Algebra and Its Applications , vol. 88/89, pp. 123–138. [7] 2.9 Cholesky Decomposition If a square matrix Ahappens to be symmetric and positive definite, then it has a special, more efficient, triangular decomposition. Symmetric means that aij=ajifor i, j =1,...,N, whilepositive definite means that v·A·v>0for all vectors v (2.9.1 ) (In Chapter 11 we will see that positive definite has the equivalent interpretation that Ahas all positive eigenvalues.) While symmetric, positive definite matrices are rather special, theyoccur quite frequently in some applications, so their special factorization, called Cholesky decomposition ,isgoodtoknowabout. Whenyoucanuseit,Choleskydecompositionisabout a factor of two faster than alternative methods for solving linear equations. Instead of seeking arbitrary lower and upper triangular factors LandU, Cholesky decomposition constructs a lower triangular matrix Lwhose transpose L Tcan itself serve as the upper triangular part. In other words we replace equation (2.3.1) by L·LT=A (2.9.2 ) This factorization is sometimes referred to as “taking the square root” of the matrix A. The components of LTare of course related to those of Lby LT ij=Lji (2.9.3 ) Writing out equation (2.9.2) in components, one readily obtains the analogs of equations (2.3.12)–(2.3.13), Lii=/parenleftBigg aii−i−1/summationdisplay k=1L2 ik/parenrightBigg1/2 (2.9.4 ) and Lji=1 Lii/parenleftBigg aij−i−1/summationdisplay k=1LikLjk/parenrightBigg j=i+1,i+2,...,N (2.9.5 )