f2-6

2.6SingularValueDecomposition 51Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).2.6 Singular Value Decomposition Thereexistsaverypowerfulsetoftechniquesfordealingwithsetsofequations ormatricesthatareeithersingularorelsenumericallyveryclosetosingular. Inmanycases where Gaussian elimination and LUdecomposition fail to give satisfactory results, this set of techniques, known as singular value decomposition ,o rSVD, will diagnose for you precisely what the problem is. In some cases, SVD will not only diagnose the problem, it will also solve it, in the sense of giving you a useful numerical answer, although, as we shall see, not necessarily “the” answerthat you thought you should get. SVDisalsothemethodofchoiceforsolvingmost linearleast-squares problems. We will outline the relevant theory in this section, but defer detailed discussion ofthe use of SVD in this application to Chapter 15, whose subject is the parametric modeling of data. SVDmethodsarebasedonthefollowingtheoremoflinearalgebra,whoseproof isbeyondourscope: Any M×NmatrixAwhosenumberofrows Misgreaterthan or equal to its number of columns N, can be written as the product of an M×N column-orthogonal matrix U,a nN×Ndiagonal matrix Wwith positive or zero elements(the singularvalues ),andthetransposeofan N×Northogonalmatrix V. Thevariousshapes ofthesematrices will bemadeclearerbythefollowingtableau:  A = U ·  w1 w2 ··· ··· wN · VT  (2.6.1 ) The matrices UandVare each orthogonal in the sense that their columns are orthonormal, M/summationdisplay i=1UikUin=δkn1≤k≤N 1≤n≤N(2.6.2 ) N/summationdisplay j=1VjkVjn=δkn1≤k≤N 1≤n≤N(2.6.3 ) 52 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).or as a tableau,  UT ·  U =  VT ·  V  =  1  (2.6.4 ) SinceVis square, it is also row-orthonormal, V·VT=1. The SVD decomposition can also be carried out when M<N. In this case the singular values wjforj=M+1,...,Nare all zero, and the corresponding columns of Uare also zero. Equation (2.6.2) then holds only for k,n≤M. The decomposition (2.6.1) can always be done, no matter how singular the matrix is, and it is “almost” unique. That is to say, it is unique up to (i) makingthe same permutation of the columns of U, elements of W, and columns of V(or rows ofV T), or (ii) forminglinear combinationsof any columns of UandVwhose correspondingelementsof Whappentobeexactlyequal. Animportantconsequence of the permutation freedom is that for the case M<N, a numerical algorithm for the decomposition need not return zero wj’s for j=M+1,...,N; theN−M zero singular values can be scattered among all positions j=1,2,...,N. At the end of this section, we give a routine, svdcmp, that performs SVD on an arbitrary matrix A, replacing it by U(they are the same shape) and returning WandVseparately. The routine svdcmpis based on a routine by Forsythe et al.[1],whichis inturnbasedontheoriginalroutineofGolubandReinsch,found,in variousforms,in [2-4]and elsewhere. These referencesincludeextensivediscussion of the algorithm used. As much as we dislike the use of black-boxroutines, we are going to ask you to accept this one, since it would take us too far afield to coverits necessary background material here. Suffice it to say that the algorithm is very stable, andthatit is veryunusualforit everto misbehave. Most oftheconceptsthat enter the algorithm (Householder reduction to bidiagonal form, diagonalization byQRprocedure with shifts) will be discussed further in Chapter 11. Ifyouareassuspiciousofblackboxesasweare,youwillwanttoverifyyourself thatsvdcmpdoeswhatwesayitdoes. Thatisveryeasytodo: Generateanarbitrary matrixA, call the routine, and then verify by matrix multiplication that (2.6.1) and (2.6.4) are satisfied. Since these two equations are the only defining requirementsfor SVD, this procedureis (for the chosen A) a complete end-to-endcheck. Now let us find out what SVD is good for. 2.6SingularValueDecomposition 53Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SVD ofa Square Matrix Ifthematrix Ais square, N×Nsay,thenU,V,andWareall squarematrices ofthesamesize. Theirinversesarealsotrivialtocompute: UandVareorthogonal, so their inverses are equal to their transposes; Wis diagonal, so its inverse is the diagonalmatrixwhoseelementsarethereciprocalsoftheelements wj. From(2.6.1) it now follows immediately that the inverse of Ais A−1=V·[diag (1/w j)]·UT(2.6.5 ) The only thing that can go wrong with this construction is for one of the wj’s to be zero, or (numerically) for it to be so small that its value is dominated by roundoff error and therefore unknowable. If more than one of the wj’s have this problem, then the matrix is even more singular. So, first of all, SVD gives you a clear diagnosis of the situation. Formally,the conditionnumber of a matrixis definedas the ratioof the largest (in magnitude) of the wj’s to the smallest of the wj’s. A matrix is singular if its condition number is infinite, and it is ill-conditioned if its condition number is too large, that is, if its reciprocal approachesthe machine’s floating-pointprecision (for example, less than 10−6for single precision or 10−12for double). For singular matrices, the concepts of nullspace andrangeare important. Consider the familiar set of simultaneous equations A·x=b (2.6.6 ) whereAis a square matrix, bandxare vectors. Equation (2.6.6) defines Aas a linear mapping from the vector space xto the vector space b.I fAis singular, then there is some subspace of x, called the nullspace, that is mapped to zero, A·x=0. The dimension of the nullspace (the number of linearly independent vectors xthat can be found in it) is called the nullityofA. Now,there is also some subspaceof bthat can be “reached”by A, in the sense thatthereexistssome xwhichismappedthere. Thissubspaceof biscalledtherange ofA. Thedimensionoftherangeiscalledthe rankofA.I fAisnonsingular,thenits rangewillbeallofthevectorspace b,soitsrankis N.I fAissingular,thentherank will be less than N. Infact, the relevanttheoremis “rankplus nullity equals N.” What has this to do with SVD? SVD explicitly constructs orthonormal bases for the nullspace and range of a matrix. Specifically, the columns of Uwhose same-numberedelements wjarenonzeroareanorthonormalsetofbasisvectorsthat span the range; the columns of Vwhose same-numbered elements wjarezeroare an orthonormal basis for the nullspace. Now let’s haveanotherlook at solving the set of simultaneouslinear equations (2.6.6)in thecase that Ais singular. First, the set of homogeneous equations,where b=0, is solved immediately by SVD: Any column of Vwhose corresponding wj is zero yields a solution. When the vector bon the right-handside is not zero, the important question is whether it lies in the range of Aor not. If it does, then the singular set of equations doeshave a solution x; in fact it has more than one solution, since any vector in the nullspace (any column of Vwith a corresponding zero wj) can be added to x in any linear combination. 54 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Ifwewanttosingleoutoneparticularmemberofthissolution-setofvectorsas arepresentative,wemightwanttopicktheonewiththesmallestlength |x|2. Hereis howtofindthatvectorusingSVD:Simply replace 1/w jbyzeroif wj=0.(Itisnot veryoftenthat onegets toset ∞=0!) Thencompute(workingfromrighttoleft) x=V·[diag (1/w j)]·(UT·b)( 2.6.7 ) This will be the solution vector of smallest length; the columns of Vthat are in the nullspace complete the specification of the solution set. Proof: Consider |x+x/prime|, wherex/primelies in the nullspace. Then, if W−1denotes the modified inverse of Wwith some elements zeroed, |x+x/prime|=/vextendsingle/vextendsingleV·W−1·UT·b+x/prime/vextendsingle/vextendsingle =/vextendsingle/vextendsingleV·(W−1·UT·b+VT·x/prime)/vextendsingle/vextendsingle =/vextendsingle/vextendsingleW−1·UT·b+VT·x/prime/vextendsingle/vextendsingle(2.6.8 ) Here the first equality follows from (2.6.7), the second and third from the orthonor- mality of V. If you now examine the two terms that make up the sum on the right-handside,youwill seethatthefirst onehasnonzero jcomponentsonlywhere wj/negationslash=0,whilethesecondone,since x/primeisinthenullspace,hasnonzero jcomponents only where wj=0. Thereforethe minimum length obtains for x/prime=0, q.e.d. Ifbisnotintherangeofthesingularmatrix A,thenthesetofequations(2.6.6) has no solution. But here is some good news: If bis not in the range of A, then equation (2.6.7) can still be used to construct a “solution” vector x. This vector x will not exactly solve A·x=b. But, among all possible vectors x, it will do the closest possible job in the least squares sense. In other words (2.6.7) finds xwhichminimizes r≡|A·x−b| (2.6.9 ) The number ris called the residualof the solution. The proofis similar to (2.6.8): Supposewe modify xbyaddingsome arbitrary x/prime. ThenA·x−bis modified by adding some b/prime≡A·x/prime. Obviously b/primeis in the range of A. We then have /vextendsingle/vextendsingleA·x−b+b/prime/vextendsingle/vextendsingle=/vextendsingle/vextendsingle(U·W·VT)·(V·W−1·UT·b)−b+b/prime/vextendsingle/vextendsingle =/vextendsingle/vextendsingle(U·W·W−1·UT−1)·b+b/prime/vextendsingle/vextendsingle =/vextendsingle/vextendsingleU·/bracketleftbig (W·W−1−1)·UT·b+UT·b/prime/bracketrightbig/vextendsingle/vextendsingle =/vextendsingle/vextendsingle(W·W−1−1)·UT·b+UT·b/prime/vextendsingle/vextendsingle(2.6.10 ) Now, (W·W−1−1)isadiagonalmatrixwhichhasnonzero jcomponentsonlyfor wj=0, whileUTb/primehas nonzero jcomponentsonlyfor wj/negationslash=0, sinceb/primelies in the range ofA. Therefore the minimum obtains for b/prime=0, q.e.d. Figure 2.6.1 summarizes our discussion of SVD thus far. 2.6SingularValueDecomposition 55Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).A ⋅ x = b SVD “solution” of A ⋅ x = csolutions of A ⋅ x = c′ solutions of A ⋅ x = dnull spaceof A SVD solution of A ⋅ x = drange of A dc (b)(a)A xb c′ Figure 2.6.1. (a) A nonsingular matrix Amaps a vector space into one of the same dimension. The vectorxis mapped into b, so thatxsatisfies the equation A·x=b. (b) A singular matrix Amaps a vector space into one of lower dimensionality, here a plane into a line, called the “range”ofA. The “nullspace ”ofAismappedto zero. Thesolutions of A·x=dconsist ofany oneparticular solution plus any vector in the nullspace, here forming a line parallel to the nullspace. Singular value decomposition(SVD) selects the particular solution closest to zero, as shown. The point clies outside of the range ofA,s oA·x=chas no solution. SVD finds the least-squares best compromise solution, namely a solution of A·x=c /prime, as shown. In the discussion since equation (2.6.6),we have been pretendingthat a matrix either is singular or else isn ’t. That is of course true analytically. Numerically, however, the far more common situation is that some of the wj’s are very small but nonzero, so that the matrix is ill-conditioned. In that case, the direct solution methods of LUdecomposition or Gaussian elimination may actually give a formal solution to the set of equations (that is, a zero pivot may not be encountered); but thesolutionvectormayhavewildlylargecomponentswhosealgebraiccancellation, when multiplying by the matrix A, may give a very poor approximation to the right-hand vector b. In such cases, the solution vector xobtained by zeroingthe 56 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).small wj’s and then using equation (2.6.7) is very often better (in the sense of the residual |A·x−b|beingsmaller)than boththedirect-methodsolution andtheSVD 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 thatwe are trying to solve. The resolution of the paradox is that we are throwing away preciselyacombinationofequationsthatissocorruptedbyroundofferrorastobeat best useless; usually it is worse than useless since it “pulls”the solution vector way off towards in finity along some directionthat is almost a nullspace vector. In doing this, it compoundstheroundoffproblemandmakes the residual |A·x−b|larger. SVD cannot be applied blindly, then. You have to exercise some discretion in decidingatwhatthresholdtozerothesmall w j’s,and/oryouhavetohavesomeidea what size of computed residual |A·x−b|is acceptable. As an example, here is a “backsubstitution ”routine svbksbfor evaluating equation (2.6.7) and obtaining a solution vector xfrom a right-hand side b,g i v e n that the SVD of a matrix Ahas already been calculated by a call to svdcmp. Note that this routine presumes that youhave already zeroed the small wj’s. It does not do this for you. If you haven’tzeroed the small wj’s, then this routine is just as ill-conditioned as any direct method, and you are misusing SVD. SUBROUTINE svbksb(u,w,v,m,n,mp,np,b,x) INTEGER m,mp,n,np,NMAXREAL b(mp),u(mp,np),v(np,np),w(np),x(np)PARAMETER (NMAX=500) Maximum anticipated value of n. Solves A·X=Bfor a vector X,w h e r e Ais specified by the arrays u,w,vas returned by svdcmp .mandnare the logical dimensions of a, and will be equal for square matrices. mp andnpare the physical dimensions of a.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. INTEGER i,j,jjREAL s,tmp(NMAX) do 12j=1,n Calculate UTB. s=0.if(w(j).ne.0.)then Nonzero result only if w jis nonzero. do11i=1,m s=s+u(i,j)*b(i) enddo 11 s=s/w(j) This is the divide by wj. endif tmp(j)=s enddo 12 do14j=1,n Matrix multiply by Vto get answer. s=0. do13jj=1,n s=s+v(j,jj)*tmp(jj) enddo 13 x(j)=s enddo 14 return END Note that a typical use of svdcmpandsvbksbsuperficially resembles the typical use of ludcmpandlubksb: In both cases, you decompose the left-hand matrixAjust once, and then can use the decomposition either once or many times withdifferentright-handsides. Thecrucialdifferenceisthe “editing”ofthesingular 2.6SingularValueDecomposition 57Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).values before svbksbis called: REAL a(np,np),u(np,np),w(np),v(np,np),b(np),x(np) ... do12i=1,n Copy ainto ui fy o ud o n ’tw a n ti tt ob ed e s t r o y e d . do11j=1,n u(i,j)=a(i,j) enddo 11 enddo 12 call svdcmp(u,n,n,np,np,w,v) SVD the square matrix a. wmax=0. Will be the maximum singular value obtained. do13j=1,n if(w(j).gt.wmax)wmax=w(j) enddo 13 wmin=wmax*1.0e-6 This is where we set the threshold for singular values allowed to be nonzero. The constant is typical,but not universal. You have to experiment withyour own application.do 14j=1,n if(w(j).lt.wmin)w(j)=0. enddo 14 call svbksb(u,w,v,n,n,np,np,b,x) Now we can backsubstitute. SVD forFewerEquations than Unknowns If you have fewer linear equations Mthan unknowns N, then you are not expecting a unique solution. Usually there will be an N−Mdimensional family of solutions. If you want to find this whole solution space, then SVD can readily do the job. The SVD decomposition will yield N−Mzero or negligible wj’s, since M<N. There may be additional zero wj’s from any degeneracies in your M equations. Be sure that you find this manysmall wj’s, and zero them before calling svbksb,whichwillgiveyoutheparticularsolutionvector x. Asbefore,thecolumns ofVcorrespondingto zeroed wj’s are the basis vectors whose linear combinations, added to the particular solution, span the solution space. SVD forMore Equations than Unknowns This situation will occur in Chapter 15, when we wish to find the least-squares solution to an overdetermined set of linear equations. In tableau, the equations to be solved are  A · x = b (2.6.11 ) The proofs that we gave above for the square case apply without modi fication to the case ofmoreequationsthan unknowns. The least-squaressolutionvector xis 58 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).given by (2.6.7), which, with nonsquare matrices, looks like this,  x = V · diag( 1/w j) ·  UT · b  (2.6.12 ) In general, the matrix Wwill not be singular, and no w j’s will need to be set to zero. Occasionally, however, there might be column degeneracies in A.I n this case you will need to zero some small wjvalues after all. The corresponding column in Vgives the linear combination of x’s that is then ill-determined even by the supposedly overdetermined set. Sometimes, although you do not need to zero any wj’s forcomputational reasons, you may nevertheless want to take note of any that are unusually small: Theircorrespondingcolumnsin Varelinearcombinationsof x’swhichareinsensitive to yourdata. In fact, youmaythenwish to zerothese wj’s, to reducethenumberof free parametersin the fit. These matters are discussed morefully in Chapter 15. Constructing an OrthonormalBasis Suppose that you have Nvectors in an M-dimensional vector space, with N≤M. Then the Nvectors span some subspace of the full vector space. Often you want to construct an orthonormal set of Nvectors that span the same subspace. The textbook way to do this is by Gram-Schmidt orthogonalization, starting with one vector and then expanding the subspace one dimension at a time. Numerically, however, because of the build-up of roundoff errors, naive Gram-Schmidt orthogonalization is terrible. The right way to construct an orthonormal basis for a subspace is by SVD: Form an M×NmatrixAwhose Ncolumns are your vectors. Run the matrix through svdcmp. The columns of the matrix U(which in fact replaces Aon output from 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, Ndimensional; the columns of U correspondingto zero wj’s should be discardedfrom the orthonormalbasis set. (QR factorization, discussed in §2.10, also constructs an orthonormal basis, see[5].) ApproximationofMatrices Note that equation (2.6.1)can be rewritten to express any matrix Aijas a sum of outer products of columns of Uand rows of VT, with the “weighting factors ” being the singular values wj, Aij=N/summationdisplay k=1wkUikVjk (2.6.13 ) 2.6SingularValueDecomposition 59Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).If you ever encounter a situation where mostof the singular values wjof a matrixAareverysmall,then Awillbewell-approximatedbyonlyafewtermsinthe sum(2.6.13). Thismeansthatyouhavetostoreonlyafewcolumnsof UandV(the samekones)andyouwill beabletorecover,withgoodaccuracy,thewholematrix. Note also that it is very ef ficient to multiply such an approximatedmatrix by a vectorx: Youjustdot xwith eachofthestoredcolumnsof V, multiplytheresulting scalar by the corresponding wk, and accumulate that multiple of the corresponding column of U. If your matrix is approximated by a small number Kof singular values, then this computationof A·xtakes onlyabout K(M+N)multiplications, instead of MNfor the full matrix. SVD Algorithm Here is the algorithm for constructing the singular value decompositionof any matrix. See §11.2–§11.3, and also [4-5], for discussion relating to the underlying method. SUBROUTINE svdcmp(a,m,n,mp,np,w,v) INTEGER m,mp,n,np,NMAX REAL a(mp,np),v(np,np),w(np) PARAMETER (NMAX=500) Maximum anticipated value of n. C USES pythag Given a matrix a(1:m,1:n) , with physical dimensions mpbynp, this routine computes its singular value decomposition, A=U·W ·VT.T h em a t r i x Ureplaces aon output. The diagonal matrix of singular values Wis output as a vector w(1:n) .T h e m a t r i x V(not the transpose VT) is output as v(1:n,1:n) . INTEGER i,its,j,jj,k,l,nm REAL anorm,c,f,g,h,s,scale,x,y,z,rv1(NMAX),pythagg=0.0 Householder reduction to bidiagonal form. scale=0.0 anorm=0.0 do 25i=1,n l=i+1 rv1(i)=scale*g g=0.0s=0.0scale=0.0 if(i.le.m)then do 11k=i,m scale=scale+abs(a(k,i)) enddo 11 if(scale.ne.0.0)then do12k=i,m a(k,i)=a(k,i)/scale s=s+a(k,i)*a(k,i) enddo 12 f=a(i,i)g=-sign(sqrt(s),f) h=f*g-s a(i,i)=f-gdo 15j=l,n s=0.0 do13k=i,m s=s+a(k,i)*a(k,j) enddo 13 f=s/hdo 14k=i,m a(k,j)=a(k,j)+f*a(k,i) enddo 14 60 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).enddo 15 do16k=i,m a(k,i)=scale*a(k,i) enddo 16 endif endif w(i)=scale *gg=0.0s=0.0 scale=0.0 if((i.le.m).and.(i.ne.n))then do 17k=l,n scale=scale+abs(a(i,k)) enddo 17 if(scale.ne.0.0)then do18k=l,n a(i,k)=a(i,k)/scale s=s+a(i,k)*a(i,k) enddo 18 f=a(i,l) g=-sign(sqrt(s),f) h=f*g-sa(i,l)=f-g do 19k=l,n rv1(k)=a(i,k)/h enddo 19 do23j=l,m s=0.0 do21k=l,n s=s+a(j,k)*a(i,k) enddo 21 do22k=l,n a(j,k)=a(j,k)+s*rv1(k) enddo 22 enddo 23 do24k=l,n a(i,k)=scale*a(i,k) enddo 24 endif endifanorm=max(anorm,(abs(w(i))+abs(rv1(i)))) enddo 25 do32i=n,1,-1 Accumulation of right-hand transformations. if(i.lt.n)then if(g.ne.0.0)then do26j=l,n Double division to avoid possible underflow. v(j,i)=(a(i,j)/a(i,l))/g enddo 26 do29j=l,n s=0.0do 27k=l,n s=s+a(i,k)*v(k,j) enddo 27 do28k=l,n v(k,j)=v(k,j)+s*v(k,i) enddo 28 enddo 29 endif do31j=l,n v(i,j)=0.0 v(j,i)=0.0 enddo 31 endif v(i,i)=1.0 2.6SingularValueDecomposition 61Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).g=rv1(i) l=i enddo 32 do39i=min(m,n),1,-1 Accumulation of left-hand transformations. l=i+1 g=w(i) do33j=l,n a(i,j)=0.0 enddo 33 if(g.ne.0.0)then g=1.0/gdo 36j=l,n s=0.0 do34k=l,m s=s+a(k,i)*a(k,j) enddo 34 f=(s/a(i,i))*gdo 35k=i,m a(k,j)=a(k,j)+f*a(k,i) enddo 35 enddo 36 do37j=i,m a(j,i)=a(j,i)*g enddo 37 else do38j= i,m a(j,i)=0.0 enddo 38 endif a(i,i)=a(i,i)+1.0 enddo 39 do49k=n,1,-1 Diagonalization of the bidiagonal form: Loop over singular values, and over allowed iterations. do48its=1,30 do41l=k,1,-1 Test for splitting. nm=l-1 Note that rv1(1) is always zero. if((abs(rv1(l))+anorm).eq.anorm) goto 2if((abs(w(nm))+anorm).eq.anorm) goto 1 enddo 41 1 c=0.0 Cancellation of rv1(l) ,i fl>1. s=1.0do 43i=l,k f=s*rv1(i) rv1(i)=c*rv1(i)if((abs(f)+anorm).eq.anorm) goto 2g=w(i) h=pythag(f,g) w(i)=hh=1.0/h c= (g*h) s=-(f*h)do 42j=1,m y=a(j,nm) z=a(j,i) a(j,nm)=(y*c)+(z*s)a(j,i)=-(y*s)+(z*c) enddo 42 enddo 43 2 z=w(k) if(l.eq.k)then Convergence. if(z.lt.0.0)then Singular value is made nonnegative. w(k)=-zdo 44j=1,n v(j,k)=-v(j,k) enddo 44 62 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).endif goto 3 endif if(its.eq.30) pause ’no convergence in svdcmp’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))/xc=1.0 Next QR transformation: s=1.0 do 47j=l,nm i=j+1 g=rv1(i) y=w(i)h=s*gg=c*g z=pythag(f,h) rv1(j)=zc=f/z s=h/z f= (x*c)+(g*s)g=-(x*s)+(g*c)h=y*s y=y*c do 45jj=1,n x=v(jj,j) z=v(jj,i) v(jj,j)= (x*c)+(z*s)v(jj,i)=-(x*s)+(z*c) enddo 45 z=pythag(f,h)w(j)=z Rotation can be arbitrary if z=0. if(z.ne.0.0)then z=1.0/z c=f*z s=h*z endif f= (c*g)+(s*y) x=-(s*g)+(c*y)do 46jj=1,m y=a(jj,j) z=a(jj,i) a(jj,j)= (y*c)+(z*s)a(jj,i)=-(y*s)+(z*c) enddo 46 enddo 47 rv1(l)=0.0 rv1(k)=f w(k)=x enddo 48 3 continue enddo 49 return END FUNCTION pythag(a,b) REAL a,b,pythag Computes (a2+b2)1/2without destructive underflow or overflow. 2.7SparseLinearSystems 63Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).REAL absa,absb absa=abs(a) absb=abs(b) if(absa.gt.absb)then pythag=absa*sqrt(1.+(absb/absa)**2) else if(absb.eq.0.)then pythag=0. else pythag=absb*sqrt(1.+(absa/absb)**2) endif endifreturn END (Double precision versions of svdcmp,svbksb, and pythag, named dsvdcmp, dsvbksb, and dpythag, are used by the routine ratlsqin§5.13. You can easily maketheconversions,orelsegettheconvertedroutinesfromthe NumericalRecipes diskette.) CITED REFERENCES AND FURTHER READING: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §8.3 and Chapter 12. Lawson, C.L., and Hanson, R. 1974, Solving Least Squares Problems (Englewood Cliffs, NJ: Prentice-Hall), Chapter 18. Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 9. [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(NewYork: Springer-Verlag), Chapter I.10 by G.H. Golub andC. Reinsch. [2] Dongarra, J.J., et al. 1979, LINPACK User’s Guide (Philadelphia: S.I.A.M.), Chapter 11. [3] Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide , 2nd ed., vol. 6 of Lecture Notes in Computer Science (New York: Springer-Verlag). Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §6.7. [4] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §5.2.6. [5] 2.7 Sparse Linear Systems A system of linear equations is called sparseif only a relatively small number of its matrix elements aijare nonzero. It is wasteful to use general methods of linear algebra on such problems, because most of the O(N3)arithmetic operations devotedtosolvingthesetofequationsorinvertingthematrixinvolvezerooperands. Furthermore, you might wish to work problems so large as to tax your available memory space, and it is wasteful to reserve storage for unfruitful zero elements.Note that there are two distinct (and not always compatible) goals for any sparse matrix method: saving time and/or saving space. We have already considered one archetypal sparse form in §2.4, the band diagonal matrix. In the tridiagonal case, e.g., we saw that it was possible to save