f2-5

2.5IterativeImprovementofaSolutiontoLinearEquations 47Sample 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).Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), Example 5.4.3, p. 166. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §9.11. Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag), Chapter I/6. [1] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §4.3. 2.5 Iterative Improvement of a Solution to Linear Equations Obviously it is not easy to obtain greater precision for the solution of a linear set than the precision of your computer’s floating-point word. Unfortunately, forlarge sets of linear equations, it is not always easy to obtain precision equal to, or even comparable to, the computer’s limit. In direct methods of solution, roundoff errors accumulate, and they are magnified to the extent that your matrix is closeto singular. You can easily lose two or three significant figures for matrices which (you thought) were farfrom singular. Ifthishappenstoyou,thereisa neattricktorestorethefullmachineprecision, callediterative improvement of thesolution. Thetheoryis verystraightforward(see Figure 2.5.1): Suppose that a vector xis the exact solution of the linear set A·x=b (2.5.1 ) You don’t, however, know x. You only know some slightly wrong solution x+δx, where δxistheunknownerror. Whenmultipliedbythematrix A,yourslightlywrong solutiongivesaproductslightlydiscrepantfromthedesiredright-handside b,namely A·(x+δx)=b+δb (2.5.2 ) Subtracting (2.5.1) from (2.5.2) gives A·δx=δb (2.5.3 ) But (2.5.2)can also be solved,trivially,for δb. Substitutingthis into (2.5.3)gives A·δx=A·(x+δx)−b (2.5.4 ) In this equation, the whole right-hand side is known, since x+δxis the wrong solution that you want to improve. It is essential to calculate the right-hand sidein double precision, since there will be a lot of cancellation in the subtraction of b. Then, we need only solve (2.5.4)for the error δx, then subtract this from the wrong solution to get an improved solution. An important extra benefit occurs if we obtained the original solution by LU decomposition. Inthiscase we alreadyhavethe LUdecomposedformof A, andall we need do to solve (2.5.4)is compute the right-handside and backsubstitute! The code to do all this is concise and straightforward: 48 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).A A−1δxx + δx x bb + δb δb Figure 2.5.1. Iterative improvement of the solution to A·x=b. Thefirst guessx+δxis multiplied by Ato produce b+δb. The known vector bis subtracted, giving δb. The linear set with this right-hand side is inverted, giving δx. This is subtracted from the first guess giving an improved solution x. SUBROUTINE mprove(a,alud,n,np,indx,b,x) INTEGER n,np,indx(n),NMAXREAL a(np,np),alud(np,np),b(n),x(n) PARAMETER (NMAX=500) Maximum anticipated value of n. C USES lubksb Improves a solution vector x(1:n) of the linear set of equations A·X=B.T h e m a t r i x a(1:n,1:n) , and the vectors b(1:n) andx(1:n) are input, as is the dimension n.A l s o input is alud ,t h eLU decomposition of aas returned by ludcmp , and the vector indx also returned by that routine. On output, only x(1:n) is modified, to an improved set of values. INTEGER i,j REAL r(NMAX) DOUBLE PRECISION sdpdo 12i=1,n Calculate the right-hand side, accumulating the resid- ual in double precision. sdp=-b(i) do11j=1,n sdp=sdp+dble(a(i,j))*dble(x(j)) enddo 11 r(i)=sdp enddo 12 call lubksb(alud,n,np,indx,r) Solve for the error term, do13i=1,n and subtract it from the old solution. x(i)=x(i)-r(i) enddo 13 return END You should note that the routine ludcmpin§2.3 destroys the input matrix as itLUdecomposesit. Since iterative improvementrequires boththe original matrix andits LUdecomposition,youwillneedtocopy Abeforecalling ludcmp. Likewise lubksbdestroysbin obtaining x, so make a copy of balso. If you don ’t mind this extra storage, iterative improvement is highlyrecommended: It is a process of order only N2operations (multiply vector by matrix, and backsubstitute —see discussion following equation 2.3.7); it never hurts; and it can really give you your money’s worth if it saves an otherwise ruined solution on which you have already spent of order N3operations. 2.5IterativeImprovementofaSolutiontoLinearEquations 49Sample 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).You can call mproveseveral times in succession if you want. Unless you are starting quite far from the true solution, one call is generally enough; but a secondcall to verify convergence can be reassuring. More onIterative Improvement It is illuminating (and will be useful later in the book) to give a somewhat more solid analyticalfoundation forequation(2.5.4),andalsotogivesomeadditionalresults. Implicitinthe previous discussion was the notion that the solution vector x+δxhas an error term; but we neglected the fact that the LUdecomposition of Ais itself not exact. A different analytical approach starts with some matrix B 0that is assumed to be an approximate inverse of the matrix A, so thatB0·Ais approximately the identity matrix 1. Define theresidual matrix RofB0as R≡1−B0·A (2.5.5 ) which is supposed to be “small”(we will be more precise below). Note that therefore B0·A=1−R (2.5.6 ) Next consider the following formal manipulation: A−1=A−1·(B−1 0·B0)=(A−1·B−1 0)·B0=(B0·A)−1·B0 =(1−R)−1·B0=(1+R+R2+R3+···)·B0(2.5.7 ) We can de fine the nth partial sum of the last expression by Bn≡(1+R+···+Rn)·B0 (2.5.8 ) so thatB∞→A−1, if the limit exists. It now is straightforward to verify that equation (2.5.8) satis fies some interesting recurrence relations. As regards solving A·x=b, wherexandbare vectors, de fine xn≡Bn·b (2.5.9 ) Then it is easy to show that xn+1=xn+B0·(b−A·xn)( 2.5.10 ) This is immediately recognizable as equation (2.5.4), with −δx=xn+1−xn, and with B0 taking the role of A−1. We see, therefore, that equation (2.5.4) does not require that the LU decomposition of Abe exact, butonly that theimplied residual Rbe small. Inrough terms,if the residual is smaller than the square root of your computer ’s roundoff error, then after one applicationofequation(2.5.10)(thatis,goingfrom x0≡B0·btox1)thefirstneglectedterm, of orderR2, will be smaller than the roundoff error. Equation (2.5.10), like equation (2.5.4), moreover, canbe appliedmore thanonce, sinceitusesonly B0,andnot anyofthehigher B’s. Amuchmoresurprisingrecurrencewhichfollowsfromequation(2.5.8)isonethatmore thandoublesthe order nat each stage: B2n+1=2Bn−Bn·A·Bn n=0,1,3,7,... (2.5.11 ) Repeated application of equation (2.5.11), from a suitable starting matrix B0, converges quadratically to the unknown inverse matrix A−1(see§9.4 for the de finition of “quadrati- cally”). Equation(2.5.11) goesby various names,including Schultz’s Method andHotelling’s Method; see Pan and Reif [1]for references. In fact, equation (2.5.11) is simply the iterative Newton-Raphson method of root- finding ( §9.4) applied to matrix inversion. Before you get too excited about equation (2.5.11), however, you should notice that it involves two full matrixmultiplications at each iteration. Each matrix multiplication involvesN 3adds and multiplies. But we already saw in §§2.1–2.3 that direct inversion of Arequires only N3adds and N3multiplies in toto. Equation (2.5.11) is therefore practical only when special circumstances allow it to be evaluated much more rapidly than is the case for generalmatrices. We will meet such circumstances later, in §13.10. 50 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).In the spirit of delayed grati fication, let us nevertheless pursue the two related issues: When does the series in equation (2.5.7) converge; and what is a suitable initial guess B0(if, for example, an initial LUdecomposition is not feasible)? We can de fine the norm of a matrix as the largest ampli fication of length that it is able to induce on a vector, /bardblR/bardbl≡ maxv/negationslash=0|R·v| |v|(2.5.12 ) Ifweletequation(2.5.7)actonsomearbitraryright-handside b,asonewantsamatrixinverse to do, it is obvious that a suf ficient condition for convergence is /bardblR/bardbl<1( 2.5.13 ) Pan and Reif [1]point out that a suitable initial guess for B0is any suf ficiently smallconstant /epsilon1times the matrix transpose of A, that is, B0=/epsilon1ATorR=1−/epsilon1AT·A (2.5.14 ) ToseewhythisissoinvolvesconceptsfromChapter11;wegivehereonlythebriefestsketch: AT·Ais a symmetric, positive de finite matrix, so it has real, positive eigenvalues. In its diagonal representation, Rtakes the form R=diag (1−/epsilon1λ1,1−/epsilon1λ2,..., 1−/epsilon1λ N)( 2.5.15 ) where all the λi’s are positive. Evidently any /epsilon1satisfying 0</epsilon1< 2/(max iλi)will give /bardblR/bardbl<1. It is not dif ficult to show that the optimal choice for /epsilon1, giving the most rapid convergence for equation (2.5.11), is /epsilon1=2/(max iλi+m i n iλi)( 2.5.16 ) Rarely does one know the eigenvalues of AT·Ain equation (2.5.16). Pan and Reif derive several interesting bounds, which are computable directly from A. The following choices guarantee the convergence of Bnasn→∞, /epsilon1≤1/slashbigg/summationdisplay j,ka2 jkor /epsilon1≤1/slashbigg/parenleftbigg max i/summationdisplay j|aij|× max j/summationdisplay i|aij|/parenrightbigg (2.5.17 ) The latter expression is truly a remarkable formula, which Pan and Reif derive by noting that the vector norm in equation (2.5.12) need not be the usual L2norm, but can instead be either theL∞(max) norm, or the L1(absolute value) norm. See their work for details. Another approach, with which we have had some success, is to estimate the largest eigenvaluestatistically,bycalculating si≡|A·vi|2forseveralunitvector vi’swithrandomly chosendirectionsin N-space. Thelargesteigenvalue λcanthenbebounded bythemaximum of2m a x siand 2NVar(si)/µ(si), where Var and µdenote the sample variance and mean, respectively. CITED REFERENCES AND FURTHER READING: Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis , 2nd ed. (Reading, MA: Addison- Wesley), §2.3.4, p. 55. Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), p. 74. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §5.5.6, p. 183. Forsythe, G.E., and Moler, C.B. 1967, Computer Solution of Linear Algebraic Systems (Engle- wood Cliffs, NJ: Prentice-Hall), Chapter 13. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §9.5, p. 437. Pan,V.,andReif,J.1985,inProceedingsoftheSeventeenthAnnualACMSymposiumonTheory of Computing (New York: Association for Computing Machinery). [1] 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”answer that 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 )