2.6 Singular Value Decomposition 67 If you ever encounter a situation where most of the singular values wi of a matrix A are very small,then A will be well-approximated by only a few terms in the sum(2.6.13).This means that you have to store only a few columns of U and V(the same 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 a vector x:You just dot x with each of the stored columns of V.multiply the resulting scalar by the corresponding w,and accumulate that multiple of the corresponding column of U.If your matrix is approximated by a small number K of singular values,then this computation of A.x takes only about K(M+N)multiplications. instead of MN for the full matrix. SVD Algorithm Here is the algorithm for constructing the singular value decomposition of any matrix.See $11.2-811.3,and also [4-51,for discussion relating to the underlying , 鱼 100 method from NUMERICAL RECIPES I 18881892 #include <math.h> #include "nrutil.h" (Nort void svdcmp(float a,int m,int n,float w[],float **v) America server computer, to make one paper /Cambridge University Press. THE Given a matrix a[1..m][1..n],this routine computes its singular value decomposition,A U.W.VT.The matrix U replaces a on output.The diagonal matrix of singular values W is out- ART put as a vector w[1..n].The matrix V(not the transpose VT)is output as v[1..n][1..n] Programs float pythag(float a,float b); int flag,i,its,j,jj,k,1,nm; float anorm,c,f,g,h,s,scale,x,y,Z,*rvi; rvi-vector(1,n); g=scale=anorm=0.0; Householder reduction to bidiagonal form to dir for(i=1:1<=n:1++){ 1=1+1; rv1[i]=scale*gi g=s=scale=0.0: rectcustsen OF SCIENTIFIC COMPUTING (ISBN 0-521 1f(1<=m)[ for (k=i;k<=m;k++)scale +fabs(a[k][i]); if (scale){ for(k=1;k<=m;k++)[ a[k][i]/scale; v@cambridge.org 1988-1992 by Numerical Recipes s+=a[k][1]*a[k][i]; Numerical Recipes books or -43108-5 f=a[i][i]; g =-SIGN(sqrt(s),f); h=f*g-s; a[i][i]=f-g; Software. for(j=1;j<=n;j++){ for(s=0.0,k=i;k<=m;k+)s+a[k][1]*a[k][j]; (outside North America) f=s/hi for (k=i;k<=m;k++)a[k][i]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 &k i !n){ for (k=l;k<=n;k++)scale +fabs(a[i][k]); if (scale){2.6 Singular Value Decomposition 67 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). If you ever encounter a situation where most of the singular values wj of a matrix A are very small, then A will be well-approximated by only a few terms in the sum (2.6.13). This means that you have to store only a few columns of U and V (the same 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 a vector x: You just dot x with each of the stored columns of V, multiply the resulting scalar by the corresponding wk, and accumulate that multiple of the corresponding column of U. If your matrix is approximated by a small number K of singular values, then this computation of A · x takes only about K(M + N) multiplications, instead of MN for the full matrix. SVD Algorithm Here is the algorithm for constructing the singular value decomposition of any matrix. See §11.2–§11.3, and also [4-5], for discussion relating to the underlying method. #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 out￾put 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) {