c15-4 (779588), страница 2
Текст из файла (страница 2)
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).computation of the matrix inverse. In theory, since AT · A is positive definite,Cholesky decomposition is the most efficient way to solve the normal equations.However, in practice most of the computing time is spent in looping over the datato form the equations, and Gauss-Jordan is quite adequate.We need to warn you that the solution of a least-squares problem directly fromthe normal equations is rather susceptible to roundoff error.
An alternative, andpreferred, technique involves QR decomposition (§2.10, §11.3, and §11.6) of thedesign matrix A. This is essentially what we did at the end of §15.2 for fitting data toa straight line, but without invoking all the machinery of QR to derive the necessaryformulas. Later in this section, we will discuss other difficulties in the least-squaresproblem, for which the cure is singular value decomposition (SVD), of which we givean implementation.
It turns out that SVD also fixes the roundoff problem, so it is ourrecommended technique for all but “easy” least-squares problems. It is for these easyproblems that the following routine, which solves the normal equations, is intended.The routine below introduces one bookkeeping trick that is quite useful inpractical work. Frequently it is a matter of “art” to decide which parameters akin a model should be fit from the data set, and which should be held constant atfixed values, for example values predicted by a theory or measured in a previousexperiment.
One wants, therefore, to have a convenient means for “freezing”and “unfreezing” the parameters ak . In the following routine the total number ofparameters ak is denoted ma (called M above). As input to the routine, you supplyan array ia[1..ma], whose components are either zero or nonzero (e.g., 1).
Zerosindicate that you want the corresponding elements of the parameter vector a[1..ma]to be held fixed at their input values. Nonzeros indicate parameters that should befitted for. On output, any frozen parameters will have their variances, and all theircovariances, set to zero in the covariance matrix.15.4 General Linear Least Squares675}for (j=2;j<=mfit;j++)Fill in above the diagonal from symmetry.for (k=1;k<j;k++)covar[k][j]=covar[j][k];gaussj(covar,mfit,beta,1);Matrix solution.for (j=0,l=1;l<=ma;l++)if (ia[l]) a[l]=beta[++j][1];Partition solution to appropriate coefficients*chisq=0.0;a.for (i=1;i<=ndat;i++) {Evaluate χ2 of the fit.(*funcs)(x[i],afunc,ma);for (sum=0.0,j=1;j<=ma;j++) sum += a[j]*afunc[j];*chisq += SQR((y[i]-sum)/sig[i]);}covsrt(covar,ma,ia,mfit);Sort covariance matrix to true order of fittingfree_vector(afunc,1,ma);coefficients.free_matrix(beta,1,ma,1,1);}That last call to a function covsrt is only for the purpose of spreading thecovariances back into the full ma × ma covariance matrix, in the proper rows andcolumns and with zero variances and covariances set for variables which wereheld frozen.The function covsrt is as follows.#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}void covsrt(float **covar, int ma, int ia[], int mfit)Expand in storage the covariance matrix covar, so as to take into account parameters that arebeing held fixed.
(For the latter, return zero covariances.){int i,j,k;float swap;for (i=mfit+1;i<=ma;i++)for (j=1;j<=i;j++) covar[i][j]=covar[j][i]=0.0;k=mfit;for (j=ma;j>=1;j--) {if (ia[j]) {for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j])for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i])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).(*funcs)(x[i],afunc,ma);ym=y[i];if (mfit < ma) {Subtract off dependences on known piecesfor (j=1;j<=ma;j++)of the fitting function.if (!ia[j]) ym -= a[j]*afunc[j];}sig2i=1.0/SQR(sig[i]);for (j=0,l=1;l<=ma;l++) {if (ia[l]) {wt=afunc[l]*sig2i;for (j++,k=0,m=1;m<=l;m++)if (ia[m]) covar[j][++k] += wt*afunc[m];beta[j][1] += ym*wt;}}676Chapter 15.Modeling of DataSolution by Use of Singular Value Decompositionfinda2χ2 = |A · a − b|that minimizes(15.4.16)Comparing to equation (2.6.9), we see that this is precisely the problem that routinessvdcmp and svbksb are designed to solve.
The solution, which is given by equation(2.6.12), can be rewritten as follows: If U and V enter the SVD decompositionof A according to equation (2.6.1), as computed by svdcmp, then let the vectorsU(i) i = 1, . . . , M denote the columns of U (each one a vector of length N ); andlet the vectors V(i) ; i = 1, . . . , M denote the columns of V (each one a vectorof length M ). Then the solution (2.6.12) of the least-squares problem (15.4.16)can be written asa=M XU(i) · bi=1wiV(i)(15.4.17)where the wi are, as in §2.6, the singular values calculated by svdcmp.Equation (15.4.17) says that the fitted parameters a are linear combinations ofthe columns of V, with coefficients obtained by forming dot products of the columnsSample 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).In some applications, the normal equations are perfectly adequate for linearleast-squares problems. However, in many cases the normal equations are very closeto singular. A zero pivot element may be encountered during the solution of thelinear equations (e.g., in gaussj), in which case you get no solution at all.
Or avery small pivot may occur, in which case you typically get fitted parameters akwith very large magnitudes that are delicately (and unstably) balanced to cancel outalmost precisely when the fitted function is evaluated.Why does this commonly occur? The reason is that, more often than experimenters would like to admit, data do not clearly distinguish between two or more ofthe basis functions provided. If two such functions, or two different combinationsof functions, happen to fit the data about equally well — or equally badly — thenthe matrix [α], unable to distinguish between them, neatly folds up its tent andbecomes singular.
There is a certain mathematical irony in the fact that least-squaresproblems are both overdetermined (number of data points greater than number ofparameters) and underdetermined (ambiguous combinations of parameters exist);but that is how it frequently is. The ambiguities can be extremely hard to noticea priori in complicated problems.Enter singular value decomposition (SVD).
This would be a good time for youto review the material in §2.6, which we will not repeat here. In the case of anoverdetermined system, SVD produces a solution that is the best approximation inthe least-squares sense, cf. equation (2.6.10). That is exactly what we want. Inthe case of an underdetermined system, SVD produces a solution whose values (forus, the ak ’s) are smallest in the least-squares sense, cf. equation (2.6.8).
That isalso what we want: When some combination of basis functions is irrelevant to thefit, that combination will be driven down to a small, innocuous, value, rather thanpushed up to delicately canceling infinities.In terms of the design matrix A (equation 15.4.4) and the vector b (equation15.4.5), minimization of χ2 in (15.4.3) can be written as15.4 General Linear Least Squares677i=1Here each ± is followed by a standard deviation.
The amazing fact is that,decomposed in this fashion, the standard deviations are all mutually independent(uncorrelated). Therefore they can be added together in root-mean-square fashion.What is going on is that the vectors V(i) are the principal axes of the error ellipsoidof the fitted parameters a (see §15.6).It follows that the variance in the estimate of a parameter aj is given by2MM XX1Vji2[V(i) ]j =σ (aj ) =wi2wi2i=1(15.4.19)i=1whose result should be identical with (15.4.14). As before, you should not besurprised at the formula for the covariances, here given without proof,Cov(aj , ak ) =M XVji Vkii=1wi2(15.4.20)We introduced this subsection by noting that the normal equations can failby encountering a zero pivot. We have not yet, however, mentioned how SVDovercomes this problem.
The answer is: If any singular value wi is zero, itsreciprocal in equation (15.4.18) should be set to zero, not infinity. (Compare thediscussion preceding equation 2.6.7.) This corresponds to adding to the fittedparameters a a zero multiple, rather than some random large multiple, of any linearcombination of basis functions that are degenerate in the fit. It is a good thing to do!Moreover, if a singular value wi is nonzero but very small, you should alsodefine its reciprocal to be zero, since its apparent value is probably an artifact ofroundoff error, not a meaningful number.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
