c2-6 (779464), страница 3
Текст из файла (страница 3)
Occasionally, however, there might be column degeneracies in A. Inthis case you will need to zero some small wj values after all. The correspondingcolumn in V gives the linear combination of x’s that is then ill-determined even bythe supposedly overdetermined set.Sometimes, although you do not need to zero any wj ’s for computationalreasons, you may nevertheless want to take note of any that are unusually small:Their corresponding columns in Vare linear combinations of x’s which are insensitiveto your data.
In fact, you may then wish to zero these wj ’s, to reduce the number offree parameters in the fit. These matters are discussed more fully in Chapter 15.Constructing an Orthonormal BasisSuppose that you have N vectors in an M -dimensional vector space, withN ≤ M . Then the N vectors span some subspace of the full vector space.Often you want to construct an orthonormal set of N vectors that span the samesubspace. The textbook way to do this is by Gram-Schmidt orthogonalization,starting with one vector and then expanding the subspace one dimension at atime.
Numerically, however, because of the build-up of roundoff errors, naiveGram-Schmidt orthogonalization is terrible.The right way to construct an orthonormal basis for a subspace is by SVD:Form an M × N matrix A whose N columns are your vectors. Run the matrixthrough svdcmp.
The columns of the matrix U (which in fact replaces A on outputfrom svdcmp) are your desired orthonormal basis vectors.You might also want to check the output wj ’s for zero values. If any occur,then the spanned subspace was not, in fact, N dimensional; the columns of Ucorresponding to zero wj ’s should be discarded from the orthonormal basis set.(QR factorization, discussed in §2.10, also constructs an orthonormal basis,see [5].)Approximation of MatricesNote that equation (2.6.1) can be rewritten to express any matrix Aij as a sumof outer products of columns of U and rows of VT , with the “weighting factors”being the singular values wj ,Aij =NXk=1wk Uik Vjk(2.6.13)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). x = 2.6 Singular Value Decomposition67SVD AlgorithmHere is the algorithm for constructing the singular value decomposition of anymatrix. See §11.2–§11.3, and also [4-5] , for discussion relating to the underlyingmethod.#include <math.h>#include "nrutil.h"void svdcmp(float **a, int m, int n, float w[], float **v)Given a matrix a[1..m][1..n], this routine computes its singular value decomposition, A =U ·W ·V T . The matrix U replaces a on output.
The diagonal matrix of singular values W is output as a vector w[1..n]. The matrix V (not the transpose V T ) is output as v[1..n][1..n].{float pythag(float a, float b);int flag,i,its,j,jj,k,l,nm;float anorm,c,f,g,h,s,scale,x,y,z,*rv1;rv1=vector(1,n);g=scale=anorm=0.0;Householder reduction to bidiagonal form.for (i=1;i<=n;i++) {l=i+1;rv1[i]=scale*g;g=s=scale=0.0;if (i <= m) {for (k=i;k<=m;k++) scale += fabs(a[k][i]);if (scale) {for (k=i;k<=m;k++) {a[k][i] /= scale;s += a[k][i]*a[k][i];}f=a[i][i];g = -SIGN(sqrt(s),f);h=f*g-s;a[i][i]=f-g;for (j=l;j<=n;j++) {for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];f=s/h;for (k=i;k<=m;k++) a[k][j] += f*a[k][i];}for (k=i;k<=m;k++) a[k][i] *= scale;}}w[i]=scale *g;g=s=scale=0.0;if (i <= m && i != n) {for (k=l;k<=n;k++) scale += fabs(a[i][k]);if (scale) {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 ever encounter a situation where most of the singular values wj of amatrix A are very small, then A will be well-approximated by only a few terms in thesum (2.6.13).
This means that you have to store only a few columns of U and V (thesame k ones) and you will be able to recover, with good accuracy, the whole matrix.Note also that it is very efficient to multiply such an approximated matrix by avector x: You just dot x with each of the stored columns of V, multiply the resultingscalar by the corresponding wk , and accumulate that multiple of the correspondingcolumn of U.
If your matrix is approximated by a small number K of singularvalues, then this computation of A · x takes only about K(M + N ) multiplications,instead of M N for the full matrix.68Chapter 2.Solution of Linear Algebraic Equations}}anorm=FMAX(anorm,(fabs(w[i])+fabs(rv1[i])));}for (i=n;i>=1;i--) {Accumulation of right-hand transformations.if (i < n) {if (g) {for (j=l;j<=n;j++)Double division to avoid possible underflow.v[j][i]=(a[i][j]/a[i][l])/g;for (j=l;j<=n;j++) {for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];for (k=l;k<=n;k++) v[k][j] += s*v[k][i];}}for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;}v[i][i]=1.0;g=rv1[i];l=i;}for (i=IMIN(m,n);i>=1;i--) {Accumulation of left-hand transformations.l=i+1;g=w[i];for (j=l;j<=n;j++) a[i][j]=0.0;if (g) {g=1.0/g;for (j=l;j<=n;j++) {for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];f=(s/a[i][i])*g;for (k=i;k<=m;k++) a[k][j] += f*a[k][i];}for (j=i;j<=m;j++) a[j][i] *= g;} else for (j=i;j<=m;j++) a[j][i]=0.0;++a[i][i];}for (k=n;k>=1;k--) {Diagonalization of the bidiagonal form: Loop overfor (its=1;its<=30;its++) {singular values, and over allowed iterations.flag=1;for (l=k;l>=1;l--) {Test for splitting.nm=l-1;Note that rv1[1] is always zero.if ((float)(fabs(rv1[l])+anorm) == anorm) {flag=0;break;}if ((float)(fabs(w[nm])+anorm) == anorm) break;}if (flag) {c=0.0;Cancellation of rv1[l], if l > 1.s=1.0;for (i=l;i<=k;i++) {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).for (k=l;k<=n;k++) {a[i][k] /= scale;s += a[i][k]*a[i][k];}f=a[i][l];g = -SIGN(sqrt(s),f);h=f*g-s;a[i][l]=f-g;for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;for (j=l;j<=m;j++) {for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];for (k=l;k<=n;k++) a[j][k] += s*rv1[k];}for (k=l;k<=n;k++) a[i][k] *= scale;2.6 Singular Value Decomposition69}}z=w[k];if (l == k) {Convergence.if (z < 0.0) {Singular value is made nonnegative.w[k] = -z;for (j=1;j<=n;j++) v[j][k] = -v[j][k];}break;}if (its == 30) nrerror("no convergence in 30 svdcmp iterations");x=w[l];Shift from bottom 2-by-2 minor.nm=k-1;y=w[nm];g=rv1[nm];h=rv1[k];f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);g=pythag(f,1.0);f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;c=s=1.0;Next QR transformation:for (j=l;j<=nm;j++) {i=j+1;g=rv1[i];y=w[i];h=s*g;g=c*g;z=pythag(f,h);rv1[j]=z;c=f/z;s=h/z;f=x*c+g*s;g = g*c-x*s;h=y*s;y *= c;for (jj=1;jj<=n;jj++) {x=v[jj][j];z=v[jj][i];v[jj][j]=x*c+z*s;v[jj][i]=z*c-x*s;}z=pythag(f,h);w[j]=z;Rotation can be arbitrary if z = 0.if (z) {z=1.0/z;c=f*z;s=h*z;}f=c*g+s*y;x=c*y-s*g;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).f=s*rv1[i];rv1[i]=c*rv1[i];if ((float)(fabs(f)+anorm) == anorm) break;g=w[i];h=pythag(f,g);w[i]=h;h=1.0/h;c=g*h;s = -f*h;for (j=1;j<=m;j++) {y=a[j][nm];z=a[j][i];a[j][nm]=y*c+z*s;a[j][i]=z*c-y*s;}70Chapter 2.Solution of Linear Algebraic Equationsfor (jj=1;jj<=m;jj++) {y=a[jj][j];z=a[jj][i];a[jj][j]=y*c+z*s;a[jj][i]=z*c-y*s;}}}free_vector(rv1,1,n);}#include <math.h>#include "nrutil.h"float pythag(float a, float b)Computes (a2 + b2 )1/2 without destructive underflow or overflow.{float absa,absb;absa=fabs(a);absb=fabs(b);if (absa > absb) return absa*sqrt(1.0+SQR(absb/absa));else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));}(Double precision versions of svdcmp, svbksb, and pythag, named dsvdcmp,dsvbksb, and dpythag, are used by the routine ratlsq in §5.13.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
