102 Chapter 2.Solution of Linear Algebraic Equations We will make use of QR decomposition,and its updating,in 89.7. CITED REFERENCES AND FURTHER READING: Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.Il of Handbook for Automatic Com- putation (New York:Springer-Verlag),Chapter 1/8.[1] Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),885.2,5.3,12.6.[2] 2.11 Is Matrix Inversion an N3 Process? We close this chapter with a little entertainment,a bit of algorithmic prestidig- itation which probes more deeply into the subject of matrix inversion.We start with a seemingly simple question: How many individual multiplications does it take to perform the matrix mul- tiplication of two 2 x 2 matrices. University Press. THE a12 (2.11.1) 121 a22 Programs Eight,right?Here they are written explicitly: 爱 OF SCIENTIFIC( 61 c11=a11×b11+a12×b21 c12=a11×b12+a12×b22 (2.11.2) c21=a21×b11+a22×b21 C22=a21×b12+a22×b22 Do you think that one can write formulas for the c's that involve only seven 是留会 10-621 4310-5 multiplications?(Try it yourself,before reading on.) Numerical Recipes Such a set of formulas was,in fact,discovered by Strassen [11.The formulas are: (outside Q1≡(a11+a22）×(b11+b22) North Software. Q2≡(a21+a22)×b11 Q3≡a11×(b12-b22) Q4三a22×(-b11+b21) (2.11.3) Q5三(a11+a12）×b22 Q6≡(-a11+a21)×(b11+b2) Qr≡(a12-a22)×(b21+b22)
102 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 machinereadable 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). We will make use of QR decomposition, 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 Computation (New York: Springer-Verlag), Chapter I/8. [1] Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins University Press), §§5.2, 5.3, 12.6. [2] 2.11 Is Matrix Inversion an N3 Process? We close this chapter with a little entertainment, a bit of algorithmic prestidigitation which probes more deeply into the subject of matrix inversion. We start with a seemingly simple question: How many individual multiplications does it take to perform the matrix multiplication of two 2 × 2 matrices, a11 a12 a21 a22 · b11 b12 b21 b22 = c11 c12 c21 c22 (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.) Such a set of formulas was, in fact, discovered by Strassen [1]. The formulas are: 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)
2.11 Is Matrix Inversion an N3 Process? 103 in terms of which c11=Q1+Q4-Q5+Q7 c21=Q2+Q4 (2.11.4) c12=Q3+Q5 c22=Q1+Q3-Q2+Q6 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)are valid when the a's and b's are themselves matrices. The problem of multiplying two very large matrices(of order N =2m for some 不 integer m)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 each hierarchical level of the recursion.In total it reduces RECIPES I the process of matrix multiplication to order Nlog7 instead of N3. 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 University Press. 斋 six times as many additions as multiplications implied by (2.11.3)(2.11.4).But, if N is very large,this constant factor is no match for the change in the exponent from N3 to Nlog2 7 With this "fast"matrix multiplication,Strassen also obtained a surprising result Programs for matrix inversion[1].Suppose that the matrices 兴 OF SCIENTIFIC( a11 a12 and C11 C12 (2.11.5) 6 a21 a22 C21 C22 are inverses of each other.Then the c's can be obtained from the a's by the following operations (compare equations 2.7.22 and 2.7.25): COMPUTING (ISBN R1 Inverse(a11) @cambri 1988-1992 by Numerical Recipes 10-621 R2=a21×R Rg=R1×a12 R4=a21×R3 (outside 膜 R5=R4-a22 North Software. R6=Inverse(Rs) (2.11.6) C12 R3 X R6 C21=R6 X R2 R7=R3×c21 c11=R1-R7 c22=-R6
2.11 Is Matrix Inversion an N 3 Process? 103 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 machinereadable 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). 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) are valid when the a’s and b’s are themselves matrices. The problem of multiplying two very large matrices (of order N = 2 m for some integer m) 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 each hierarchical level of the recursion. In total it reduces the process of matrix multiplication to order N log2 7 instead of N3. 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, if N is very large, this constant factor is no match for the change in the exponent from N3 to Nlog2 7. With this “fast” matrix multiplication, Strassen also obtained a surprising result for matrix inversion [1]. Suppose that the matrices a11 a12 a21 a22 and c11 c12 c21 c22 (2.11.5) are inverses of each other. Then the c’s can be obtained from the a’s by the following 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)
104 Chapter 2.Solution of Linear Algebraic Equations In(2.11.6)the"inverse"operator occurs just twice.It is to be interpreted as the reciprocal if the a's and c's are scalars,but as matrix inversion if the a's and c's are themselves submatrices.Imagine doing the inversion of a very large matrix,of order 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 N divisions in all!So divisions don't dominate in the recursive use of(2.11.6).Equation (2.11.6) is dominated,in fact,by its 6 multiplications.Since these can be done by an Nlog27 algorithm,so can the matrix inversion! This is fun,but let's look at practicalities:If you estimate how large N has to be before the difference between exponent 3 and exponent log2 7=2.807 is substantial enough to outweigh the bookkeeping overhead.arising from the complicated nature of the recursive Strassen algorithm,you will find that LU decomposition 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 83g 18881992 1-800 multiply the complex numbers (a+ib)and(c+id)in only three real multiplications? [Answer:see $5.4.](2)Can you evaluate a general fourth-degree polynomial in to any from NUMERICAL RECIPES IN x for many different values of x with only three multiplications per evaluation? Answer: see 85.3.1 CITED REFERENCES AND FURTHER READING: (North America 州bMe se to make one paper THE Strassen.V.1969,Numerische Mathematik,vol.13,pp.354-356.[1] ART Kronsjo,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 Programs 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,SIAM Review,vol.26,pp.393-415.[More recent results that show that an exponent of 2.496 can be achieved-theoretically!] Ito directcustserv@cambridge.org (outside North America). 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) Software
104 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 machinereadable 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). In (2.11.6) the “inverse” operator occurs just twice. It is to be interpreted as the reciprocal if the a’s and c’s are scalars, but as matrix inversion if the a’s and c’s are themselves submatrices. Imagine doing the inversion of a very large matrix, of order 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 N divisions in all! So divisions don’t dominate in the recursive use of (2.11.6). Equation (2.11.6) is dominated, in fact, by its 6 multiplications. Since these can be done by an N log2 7 algorithm, so can the matrix inversion! This is fun, but let’s look at practicalities: If you estimate how large N has to be before the difference between exponent 3 and exponent log 2 7=2.807 is substantial enough to outweigh the bookkeeping overhead, arising from the complicated nature of the recursive Strassen algorithm, you will find that LU decomposition 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 multiply the complex numbers(a+ib) and (c+id)in only three real multiplications? [Answer: see §5.4.] (2) Can you evaluate a general fourth-degree polynomial in x for many different values of x with only three multiplications 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, SIAM Review, vol. 26, pp. 393–415. [More recent results that show that an exponent of 2.496 can be achieved — theoretically!]