64 Chapter 2.Solution of Linear Algebraic Equations small wi's and then using equation (2.6.7)is very often better(in the sense of the residual A.x-b being smaller)than both the direct-method solution and the SVD solution where the small wi's are left nonzero. It may seem paradoxical that this can be so,since zeroing a singular value corresponds to throwing away one linear combination of the set of equations that we are trying to solve.The resolution of the paradox is that we are throwing away precisely a combination ofequations that is so corrupted by roundofferror as to be at best useless;usually it is worse than useless since it"pulls"the solution vector way off towards infinity along some direction that is almost a nullspace vector.In doing this,it compounds the roundoff problem and makes the residual A.x-b larger. 81 SVD cannot be applied blindly,then.You have to exercise some discretion in deciding at what threshold to zero the small wi's,and/or you have to have some idea nted for what size of computed residual A.x-b is acceptable. As an example,here is a "backsubstitution"routine svbksb for evaluating equation(2.6.7)and obtaining a solution vector x from a right-hand side b,given that the SVD of a matrix A has already been calculated by a call to svdcmp.Note ERICA that this routine presumes that you have already zeroed the small w,'s.It does not do this for you.If you haven't zeroed the small wi's,then this routine is just as 令 ill-conditioned as any direct method,and you are misusing SVD. #include "nrutil.h" void svbksb(float **u,float w[,float **v,int m,int n,float b[],float x) Solves A.X=B for a vector X,where A is specified by the arrays u[1..m][1..n],w[1..n], 9 Program v[1..n][1..n]as returned by svdcmp.m and n are the dimensions of a,and will be equal for square matrices.b[1..m]is the input right-hand side.x [1..n]is the output solution vector. No input quantities are destroyed,so the routine may be called sequentially with different b's. int jj,j,i; float s,*tmpi tmp=vector(1,n); for(j=1;j<=n;j++){ Calculate UTB. s=0.0; if (w[j]){ Nonzero result only if wj is nonzero. for(i=1;i<=m;1++)s+=u[1][j门*b[i]； v@cambri s/=w[j]; This is the divide by wj tmp[j]=s; 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING (ISBN 10621 43108 d for(j=1;j<=n;j++){ Matrix multiply by V to get answer. s=0.0; for (jj=1;jj<=n;jj++)s 膜 +v[j][jj]*tmp[jj]; x[j]=s; (outside North Amer Software. free_vector(tmp,1,n); Note that a typical use of svdcmp and svbksb superficially resembles the typical use of ludcmp and lubksb:In both cases,you decompose the left-hand matrix A just once,and then can use the decomposition either once or many times with different right-hand sides.The crucial difference is the "editing"of the singular values before svbksb is called:64 Chapter 2. Solution of Linear Algebraic Equations 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). small wj ’s and then using equation (2.6.7) is very often better (in the sense of the residual |A · x − b| being smaller) than both the direct-method solution and the SVD solution where the small wj ’s are left nonzero. It may seem paradoxical that this can be so, since zeroing a singular value corresponds to throwing away one linear combination of the set of equations that we are trying to solve. The resolution of the paradox is that we are throwing away precisely a combination of equations that is so corrupted by roundoff error as to be at best useless; usually it is worse than useless since it “pulls” the solution vector way off towards infinity along some direction that is almost a nullspace vector. In doing this, it compounds the roundoff problem and makes the residual |A · x − b| larger. SVD cannot be applied blindly, then. You have to exercise some discretion in deciding at what threshold to zero the small wj ’s, and/or you have to have some idea what size of computed residual |A · x − b| is acceptable. As an example, here is a “backsubstitution” routine svbksb for evaluating equation (2.6.7) and obtaining a solution vector x from a right-hand side b, given that the SVD of a matrix A has already been calculated by a call to svdcmp. Note that this routine presumes that you have already zeroed the small wj ’s. It does not do this for you. If you haven’t zeroed the small wj ’s, then this routine is just as ill-conditioned as any direct method, and you are misusing SVD. #include "nrutil.h" void svbksb(float **u, float w[], float **v, int m, int n, float b[], float x[]) Solves A·X = B for a vector X, where A is specified by the arrays u[1..m][1..n], w[1..n], v[1..n][1..n] as returned by svdcmp. m and n are the dimensions of a, and will be equal for square matrices. b[1..m] is the input right-hand side. x[1..n] is the output solution vector. No input quantities are destroyed, so the routine may be called sequentially with different b’s. { int jj,j,i; float s,*tmp; tmp=vector(1,n); for (j=1;j<=n;j++) { Calculate U T B. s=0.0; if (w[j]) { Nonzero result only if wj is nonzero. for (i=1;i<=m;i++) s += u[i][j]*b[i]; s /= w[j]; This is the divide by wj . } tmp[j]=s; } for (j=1;j<=n;j++) { Matrix multiply by V to get answer. s=0.0; for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj]; x[j]=s; } free_vector(tmp,1,n); } Note that a typical use of svdcmp and svbksb superficially resembles the typical use of ludcmp and lubksb: In both cases, you decompose the left-hand matrix A just once, and then can use the decomposition either once or many times with different right-hand sides. The crucial difference is the “editing” of the singular values before svbksb is called: