c2-9 (779467), страница 3
Текст из файла (страница 3)
. . , N ). The Toeplitz matrix neednot be symmetric. y[1..n] and r[1..2*n-1] are input arrays; x[1..n] is the output array.{int j,k,m,m1,m2;float pp,pt1,pt2,qq,qt1,qt2,sd,sgd,sgn,shn,sxn;float *g,*h;96Chapter 2.Solution of Linear Algebraic EquationsPapers by Bunch [6] and de Hoog [7] will give entry to the literature.CITED REFERENCES AND FURTHER READING:Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), Chapter 5 [also treats some other special forms].Westlake, J.R. 1968, A Handbook of Numerical Matrix Inversion and Solution of Linear Equations(New York: Wiley).
[2]von Mises, R. 1964, Mathematical Theory of Probability and Statistics (New York: AcademicPress), pp. 394ff. [3]Levinson, N., Appendix B of N. Wiener, 1949, Extrapolation, Interpolation and Smoothing ofStationary Time Series (New York: Wiley). [4]Robinson, E.A., and Treitel, S. 1980, Geophysical Signal Analysis (Englewood Cliffs, NJ: PrenticeHall), pp. 163ff. [5]Bunch, J.R. 1985, SIAM Journal on Scientific and Statistical 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 DecompositionIf a square matrix A happens to be symmetric and positive definite, then it has aspecial, more efficient, triangular decomposition. Symmetric means that aij = aji fori, j = 1, . . . , N , while positive definite means thatv · A · v > 0 for all vectors v(2.9.1)(In Chapter 11 we will see that positive definite has the equivalent interpretation that A hasall positive eigenvalues.) While symmetric, positive definite matrices are rather special, theyoccur quite frequently in some applications, so their special factorization, called Choleskydecomposition, is good to know about. When you can use it, Cholesky decomposition is abouta factor of two faster than alternative methods for solving linear equations.Instead of seeking arbitrary lower and upper triangular factors L and U, Choleskydecomposition constructs a lower triangular matrix L whose transpose LT can itself serve asthe upper triangular part.
In other words we replace equation (2.3.1) byL · LT = A(2.9.2)This factorization is sometimes referred to as “taking the square root” of the matrix A. Thecomponents of LT are of course related to those of L byLTij = 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),!1/2i−1XLii = aii −L2ik(2.9.4)k=1andLji1=Liiaij −i−1Xk=1!Lik Ljkj = i + 1, i + 2, . . .
, N(2.9.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).Forsythe, G.E., and Moler, C.B.
1967, Computer Solution of Linear Algebraic Systems (Englewood Cliffs, NJ: Prentice-Hall), §19. [1]2.9 Cholesky Decomposition97#include <math.h>void choldc(float **a, int n, float p[])Given a positive-definite symmetric matrix a[1..n][1..n], this routine constructs its Choleskydecomposition, A = L · LT . On input, only the upper triangle of a need be given; it is notmodified. The Cholesky factor L is returned in the lower triangle of a, except for its diagonalelements which are returned in p[1..n].{void nrerror(char error_text[]);int i,j,k;float sum;for (i=1;i<=n;i++) {for (j=i;j<=n;j++) {for (sum=a[i][j],k=i-1;k>=1;k--) sum -= a[i][k]*a[j][k];if (i == j) {if (sum <= 0.0)a, with rounding errors, is not positive definite.nrerror("choldc failed");p[i]=sqrt(sum);} else a[j][i]=sum/p[i];}}}You might at this point wonder about pivoting.
The pleasant answer is that Choleskydecomposition is extremely stable numerically, without any pivoting at all. Failure of choldcsimply indicates that the matrix A (or, with roundoff error, another very nearby matrix) isnot positive definite. In fact, choldc is an efficient way to test whether a symmetric matrixis positive definite.
(In this application, you will want to replace the call to nrerror withsome less drastic signaling method.)Once your matrix is decomposed, the triangular factor can be used to solve a linearequation by backsubstitution. The straightforward implementation of this isvoid cholsl(float **a, int n, float p[], float b[], float x[])Solves the set of n linear equations A · x = b, where a is a positive-definite symmetric matrix.a[1..n][1..n] and p[1..n] are input as the output of the routine choldc. Only the lowertriangle of a is accessed. b[1..n] is input as the right-hand side vector.
The solution vector isreturned in x[1..n]. a, n, and p are not modified and can be left in place for successive callswith different right-hand sides b. b is not modified unless you identify b and x in the callingsequence, which is allowed.{int i,k;float sum;for (i=1;i<=n;i++) {Solve L · y = b, storing y in x.for (sum=b[i],k=i-1;k>=1;k--) sum -= a[i][k]*x[k];x[i]=sum/p[i];}for (i=n;i>=1;i--) {Solve LT · x = y.for (sum=x[i],k=i+1;k<=n;k++) sum -= a[k][i]*x[k];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).If you apply equations (2.9.4) and (2.9.5) in the order i = 1, 2, . . . , N , you will seethat the L’s that occur on the right-hand side are already determined by the time they areneeded. Also, only components aij with j ≥ i are referenced. (Since A is symmetric,these have complete information.) It is convenient, then, to have the factor L overwrite thesubdiagonal (lower triangular but not including the diagonal) part of A, preserving the inputupper triangular values of A. Only one extra vector of length N is needed to store the diagonalpart of L. The operations count is N 3 /6 executions of the inner loop (consisting of onemultiply and one subtract), with also N square roots.
As already mentioned, this is about afactor 2 better than LU decomposition of A (where its symmetry would be ignored).A straightforward implementation is98Chapter 2.Solution of Linear Algebraic Equationsx[i]=sum/p[i];}}for (i=1;i<=n;i++) {a[i][i]=1.0/p[i];for (j=i+1;j<=n;j++) {sum=0.0;for (k=i;k<j;k++) sum -= a[j][k]*a[k][i];a[j][i]=sum/p[j];}}CITED REFERENCES AND FURTHER READING:Wilkinson, J.H., and Reinsch, C.
1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag), Chapter I/1.Gill, P.E., Murray, W., and Wright, M.H. 1991, Numerical Linear Algebra and Optimization, vol. 1(Redwood City, CA: Addison-Wesley), §4.9.2.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§5.3.5.Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed.
(Baltimore: Johns HopkinsUniversity Press), §4.2.2.10 QR DecompositionThere is another matrix factorization that is sometimes very useful, the so-called QRdecomposition,A = Q·R(2.10.1)Here R is upper triangular, while Q is orthogonal, that is,QT · Q = 1(2.10.2)where QT is the transpose matrix of Q. Although the decomposition exists for a generalrectangular matrix, we shall restrict our treatment to the case when all the matrices are square,with dimensions N × N .Like the other matrix factorizations we have met (LU , SVD, Cholesky), QR decomposition can be used to solve systems of linear equations.
To solveA·x =b(2.10.3)R · x = QT · b(2.10.4)first form Q · b and then solveTby backsubstitution. Since QR decomposition involves about twice as many operations asLU decomposition, it is not used for typical systems of linear equations. However, we willmeet special cases where QR is the method 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).A typical use of choldc and cholsl is in the inversion of covariance matrices describingthe fit of data to a model; see, e.g., §15.6. In this, and many other applications, one often needsL−1 . The lower triangle of this matrix can be efficiently found from the output of choldc:.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
