f2-11

2.11IsMatrix Inversionan N3Process? 95Sample 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).return END SUBROUTINE rotate(r,qt,n,np,i,a,b) INTEGER n,np,iREAL a,b,r(np,np),qt(np,np) Given n×nmatrices randqtof physical dimension np, carry out a Jacobi rotation on rows i andi+1 of each matrix. aandbare the parameters of the rotation: cosθ=a/√ a2+b2, sinθ=b/√ a2+b2. INTEGER j REAL c,fact,s,w,y if(a.eq.0.)then Avoid unnecessary overflow or underflow. c=0.s=sign(1.,b) else if(abs(a).gt.abs(b))then fact=b/ac=sign(1./sqrt(1.+fact**2),a) s=fact*c else fact=a/bs=sign(1./sqrt(1.+fact**2),b) c=fact*s endifdo 11j=i,n Premultiply rby Jacobi rotation. y=r(i,j) w=r(i+1,j)r(i,j)=c*y-s*wr(i+1,j)=s*y+c*w enddo 11 do12j=1,n Premultiply qtby Jacobi rotation. y=qt(i,j) w=qt(i+1,j) qt(i,j)=c*y-s*wqt(i+1,j)=s*y+c*w enddo 12 returnEND We will make use of QRdecomposition, and its updating, in §9.7. CITED REFERENCES AND FURTHER READING: Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag), Chapter I/8. [1] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §§5.2, 5.3, 12.6. [2] 2.11 Is Matrix Inversion an N3Process? We close this chapter with a little entertainment, a bit of algorithmicprestidig- itation which probes more deeply into the subject of matrix inversion. We start with a seemingly simple question: 96 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).How many individual multiplications does it take to perform the matrix mul- tiplication of two 2×2matrices, /parenleftbigg a11 a12 a21 a22/parenrightbigg ·/parenleftbigg b11 b12 b21 b22/parenrightbigg =/parenleftbigg c11 c12 c21 c22/parenrightbigg (2.11.1 ) Eight, right? Here they are written explicitly: c11=a11×b11+a12×b21 c12=a11×b12+a12×b22 c21=a21×b11+a22×b21 c22=a21×b12+a22×b22(2.11.2 ) Do you think that one can write formulas for the c’s that involve only seven multiplications? (Try it yourself, before reading on.) Suchasetofformulaswas,infact,discoveredbyStrassen [1]. Theformulasare: Q1≡(a11+a22)×(b11+b22) Q2≡(a21+a22)×b11 Q3≡a11×(b12−b22) Q4≡a22×(−b11+b21) Q5≡(a11+a12)×b22 Q6≡(−a11+a21)×(b11+b12) Q7≡(a12−a22)×(b21+b22)(2.11.3 ) in terms of which c11=Q1+Q4−Q5+Q7 c21=Q2+Q4 c12=Q3+Q5 c22=Q1+Q3−Q2+Q6(2.11.4 ) What’s the use of this? There is one fewer multiplication than in equation (2.11.2), but many more additions and subtractions. It is not clear that anything has been gained. But notice that in (2.11.3) the a’s and b’s are never commuted. Therefore(2.11.3)and(2.11.4)arevalidwhenthe a’sand b’sarethemselvesmatrices. The problem of multiplying two very large matrices (of order N=2mfor some integerm) can now be broken down recursively by partitioning the matrices into quarters, sixteenths, etc. And note the key point: The savings is not just a factor “7/8”; it is that factor at eachhierarchical level of the recursion. In total it reduces the process of matrix multiplication to order Nlog27instead of N3. 2.11Is MatrixInversionan N3Process? 97Sample 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).What about all the extra additions in (2.11.3)–(2.11.4)? Don’t they outweigh the advantage of the fewer multiplications? For large N, it turns out that there are six times as many additions as multiplications implied by (2.11.3)–(2.11.4). But, ifNis very large, this constant factor is no match for the change in the exponent from N3toNlog27. Withthis“fast”matrixmultiplication,Strassenalsoobtainedasurprisingresult for matrix inversion [1]. Suppose that the matrices /parenleftbigg a11 a12 a21 a22/parenrightbigg and/parenleftbigg c11 c12 c21 c22/parenrightbigg (2.11.5 ) areinversesofeachother. Thenthe c’scanbeobtainedfromthe a’sbythefollowing operations (compare equations 2.7.22 and 2.7.25): R1=Inverse (a11) R2=a21×R1 R3=R1×a12 R4=a21×R3 R5=R4−a22 R6=Inverse (R5) c12=R3×R6 c21=R6×R2 R7=R3×c21 c11=R1−R7 c22=−R6(2.11.6 ) In(2.11.6)the“inverse”operatoroccursjusttwice. Itistobeinterpretedasthe reciprocal if the a’s and c’s are scalars, but as matrix inversion if the a’s and c’s are themselvessubmatrices. Imaginedoingtheinversionofaverylargematrix,oforder N=2m, recursively by partitions in half. At each step, halving the order doubles the number of inverse operations. But this means that there are only Ndivisions in all! So divisions don’t dominate in the recursive use of (2.11.6). Equation (2.11.6) is dominated,infact,byits 6multiplications. Since thesecanbedonebyan Nlog27 algorithm, so can the matrix inversion! Thisisfun,butlet’slookatpracticalities: Ifyouestimatehowlarge Nhastobe beforethedifferencebetweenexponent3andexponent log27=2 .807issubstantial enoughto outweighthe bookkeepingoverhead,arisingfromthe complicatednature of the recursive Strassen algorithm, you will find that LUdecomposition is in no immediate danger of becoming obsolete. If, on the other hand, you like this kind of fun, then try these: (1) Can you multiplythecomplexnumbers (a+ib)and (c+id)inonlythreerealmultiplications? [Answer: see §5.4.] (2) Can you evaluate a general fourth-degree polynomial in 98 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).xfor many different values of xwith only threemultiplications per evaluation? [Answer: see §5.3.] CITED REFERENCES AND FURTHER READING: Strassen, V. 1969, Numerische Mathematik , vol. 13, pp. 354–356. [1] Kronsj¨o, L. 1987, Algorithms: Their Complexity and Efficiency , 2nd ed. (New York: Wiley). Winograd, S. 1971, Linear Algebra and Its Applications , vol. 4, pp. 381–388. Pan, V. Ya. 1980, SIAM Journal on Computing , vol. 9, pp. 321–342. Pan, V. 1984, How to Multiply Matrices Faster , Lecture Notes in Computer Science, vol. 179 (New York: Springer-Verlag) Pan,V.1984, SIAMReview ,vol.26,pp.393–415.[Morerecentresultsthatshowthatanexponent of 2.496 can be achieved — theoretically!]