chap31.htm

程序设计的艺术,此书是英文版的,经久不衰,非常的经典

<h1><a name="096f_1a8a">31.2 Strassen's algorithm for matrix multiplication<a name="096f_1a8a"></h1><P>
<a name="096f_1a88"><a name="096f_1a89">This section presents Strassen's remarkable recursive algorithm for multiplying<I> n<FONT FACE="Times New Roman" SIZE=3> </I><img src="../images/mult.gif"><I> n </I></FONT>matrices that runs in <img src="../images/bound.gif">(<I>n</I><SUP>lg 7</SUP>) = <I>O</I>(<I>n</I><SUP>2.81</SUP>) time. For sufficiently large<I> n,</I> therefore, it outperforms the naive <img src="../images/bound.gif">(<I>n</I><SUP>3</SUP>) matrix-multiplication algorithm <FONT FACE="Courier New" SIZE=2>MATRIX</FONT>-<FONT FACE="Courier New" SIZE=2>MULTIPLY</FONT> from Section 26.1.<P>





<h2>An overview of the algorithm</h2><P>
<a name="0970_1a8a">Strassen<FONT FACE="CG Times (W1)" SIZE=2>'</FONT>s algorithm can be viewed as an application of a familiar design technique: divide and conquer. Suppose we wish to compute the product <I>C</I> = <I>AB</I>, where each of <I>A</I>, <I>B</I>, and <I>C</I> are <I>n<FONT FACE="Times New Roman" SIZE=3> </I><img src="../images/mult.gif"><I> n</I></FONT> matrices. Assuming that <I>n </I>is an exact power of 2, we divide each of <I>A, B,</I> and <I>C</I> into four <I>n/2<FONT FACE="Times New Roman" SIZE=3> <img src="../images/mult.gif"> n/2 </I></FONT>matrices, rewriting the equation <I>C</I> = <I>AB</I> as follows:<P>
<img src="739_a.gif"><P>
<h4><a name="0970_1a8b">(31.9)<a name="0970_1a8b"></sub></sup></h4><P>
(Exercise 31.2-2 deals with the situation in which <I>n</I> is not an exact power of 2.) For convenience, the submatrices of <I>A</I> are labeled alphabetically from left to right, whereas those of <I>B</I> are labeled from top to bottom, in agreement with the way matrix multiplication is performed. Equation (31.9) corresponds to the four equations<P>
<pre><I>r </I>= <I>ae</I> + <I>bf</I> ,</sub></sup></pre><P>
<h4><a name="0970_1a8c">(31.10)<a name="0970_1a8c"></sub></sup></h4><P>
<pre><I>s = ag </I>+<I> bh ,</I></sub></sup></pre><P>
<h4><a name="0970_1a8d">(31.11)<a name="0970_1a8d"></sub></sup></h4><P>
<pre><I>t = ce </I>+<I> df ,</I></sub></sup></pre><P>
<h4><a name="0970_1a8e">(31.12)<a name="0970_1a8e"></sub></sup></h4><P>
<pre><I>u = cg </I>+<I> dh .</I></sub></sup></pre><P>
<h4><a name="0970_1a8f">(31.13)<a name="0970_1a8f"></sub></sup></h4><P>
Each of these four equations specifies two multiplications of <I>n/2 </I><img src="../images/mult.gif"><I><FONT FACE="Times New Roman" SIZE=3> n/2 </I></FONT>matrices and the addition of their <I>n/2<FONT FACE="Times New Roman" SIZE=3> </I><img src="../images/mult.gif"><I> n/2</I></FONT> products. Using these equations to define a straightforward divide-and-conquer strategy, we derive the following recurrence for the time <I>T</I>(<I>n</I>) to multiply two <I>n<FONT FACE="Times New Roman" SIZE=3> </I><img src="../images/mult.gif"><I> n</I></FONT> matrices:<P>
<pre><I>T</I>(<I>n</I>) = 8<I>T</I>(<I>n</I>/2) + <img src="../images/bound.gif">(<I>n</I><SUP>2</SUP>).</sub></sup></pre><P>
<h4><a name="0970_1a90">(31.14)<a name="0970_1a90"></sub></sup></h4><P>
Unfortunately, recurrence (31.14) has the solution<I> T</I>(<I>n</I>) = <img src="../images/bound.gif">(<I>n</I><SUP>3</SUP>), and thus this method is no faster than the ordinary one.<P>
Strassen discovered a different recursive approach that requires only 7 recursive multiplications of <I>n/2<FONT FACE="Times New Roman" SIZE=3> </I><img src="../images/mult.gif"><I> n</I></FONT>/2<I> </I>matrices and <img src="../images/bound.gif">(<I>n</I><SUP>2</SUP>) scalar additions and subtractions, yielding the recurrence<P>
<pre><I>T</I>(<I>n</I>) = 7<I>T</I>(<I>n</I>/2) + <img src="../images/bound.gif">(<I>n</I><SUP>2</SUP>)</sub></sup></pre><P>
<pre>= <img src="../images/bound.gif">(<I>n</I><SUP>lg 7</SUP>)</sub></sup></pre><P>
<pre>= <I>O</I>(<I>n</I><SUP>2.81</SUP>) .</sub></sup></pre><P>
<h4><a name="0970_1a91">(31.15)<a name="0970_1a91"></sub></sup></h4><P>
Strassen's method has four steps:<P>
1.     Divide the input matrices <I>A</I> and <I>B</I> into <I>n</I>/2<I><FONT FACE="Times New Roman" SIZE=3> </I><img src="../images/mult.gif"><I> n</I></FONT>/2 submatrices, as in equation (31.9).<P>
2.     Using <img src="../images/bound.gif">(<I>n</I><SUP>2</SUP>) scalar additions and subtractions, compute 14 <I>n</I>/2<I><FONT FACE="Times New Roman" SIZE=3> <img src="../images/mult.gif"> n</I></FONT>/2<I> </I>matrices <I>A</I><SUB>1</SUB>, <I>B</I><SUB>1</SUB>,<I> A</I><SUB>2</SUB>,<I> B</I><SUB>2</SUB>,<I> . . . </I>,<I> A</I><SUB>7</SUB>,<I> B</I><SUB>7</sub>.<P>
3.     Recursively compute the seven matrix products <I>P<SUB>i</SUB> = A<SUB>i</SUB>B<SUB>i</I></SUB> for <I>i </I>= 1, 2, . . . , 7.<P>
4.     Compute the desired submatrices <I>r, s, t, u </I>of the result matrix <I>C </I>by adding and/or subtracting various combinations of the <I>P<SUB>i</I> </SUB>matrices, using only <img src="../images/bound.gif">(<I>n<SUP>2</I></SUP>) scalar additions and subtractions.<P>
Such a procedure satisfies the recurrence (31.15). All that we have to do now is fill in the missing details.<P>
<P>







<h2>Determining the submatrix products</h2><P>
It is not clear exactly how Strassen discovered the submatrix products that are the key to making his algorithm work. Here, we reconstruct one plausible discovery method.<P>
Let us guess that each matrix product <I>P<SUB>i</I></SUB> can be written in the form<P>
<pre><I>P<SUB>i  </I></SUB>=<I>  A<SUB>i</SUB>B<SUB>i</I></sub></sup></pre><P>
<pre>=  (<img src="../images/alpha12.gif"><I><SUB>i</I>1</SUB><I>a</I> + <img src="../images/alpha12.gif"><I><SUB>i</I>2</SUB><I>b</I> + <img src="../images/alpha12.gif"><I><SUB>i</I>3</SUB><I>c</I> + <img src="../images/alpha12.gif"><I><SUB>i</I>4</SUB><I>d</I>) <img src="../images/dot10.gif"> (<img src="../images/beta14.gif"><I><SUB>i</I>1</SUB><I>e</I> + <img src="../images/beta14.gif"><I><SUB>i</I>2</SUB><I>f</I> + <img src="../images/beta14.gif"><I><SUB>i</I>3</SUB><I>g</I> + <img src="../images/beta14.gif"><I><SUB>i</I>4</SUB><I>h</I>) ,</sub></sup></pre><P>
<h4><a name="0971_0001">(31.16)<a name="0971_0001"></sub></sup></h4><P>
where the coefficients <img src="../images/alpha12.gif"><I><SUB>ij</I></SUB>, <img src="../images/beta14.gif"><I><SUB>ij</I> </SUB>are all drawn from the set {-1, 0, 1}. That is, we guess that each product is computed by adding or subtracting some of the submatrices of <I>A,</I> adding or subtracting some of the submatrices of <I>B</I>, and then multiplying the two results together. While more general strategies are possible, this simple one turns out to work.<P>
If we form all of our products in this manner, then we can use this method recursively without assuming commutativity of multiplication, since each product has all of the <I>A</I> submatrices on the left and all of the <I>B</I> submatrices on the right. This property is essential for the recursive application of this method, since matrix multiplication is not commutative.<P>
For convenience, we shall use 4 <img src="../images/mult.gif"> 4 matrices to represent linear combinations of products of submatrices, where each product combines one submatrix of <I>A</I> with one submatrix of <I>B</I> as in equation (31.16). For example, we can rewrite equation (31.10) as<P>
<img src="740_a.gif"><P>
<img src="741_a.gif"><P>
The last expression uses an abbreviated notation in which "+" represents +1, "<img src="../images/dot10.gif"><FONT FACE="CG Times (W1)" SIZE=2></FONT> represents 0, and <FONT FACE="CG Times (W1)" SIZE=2>"</FONT>-" represents -1. (From here on, we omit the row and column labels.) Using this notation, we have the following equations for the other submatrices of the result matrix <I>C</I>:<P>
<img src="741_b.gif"><P>
We begin our search for a faster matrix-multiplication algorithm by observing that the submatrix <I>s</I> can be computed as <I>s</I> = <I>P</I><SUB>1</SUB> + <I>P</I><SUB>2</SUB>, where <I>P</I><SUB>1 </SUB>and <I>P</I><SUB>2</SUB> are computed using one matrix multiplication each:<P>
<img src="741_c.gif"><P>
The matrix <I>t</I> can be computed in a similar manner as <I>t</I> = <I>P</I><SUB>3</SUB> + <I>P</I><SUB>4</SUB>, where<P>
<img src="742_a.gif"><P>
Let us define an <I><B>essential term</I></B> to be one of the eight terms appearing on the right-hand side of one of the equations (31.10)--(31.13). We have now used 4 products to compute the two submatrices <I>s</I> and <I>t</I> whose essential terms are <I>ag</I>, <I>bh</I>, <I>ce</I>, and <I>df</I> . Note that <I>P</I><SUB>1</SUB> computes the essential term <I>ag</I>, <I>P</I><SUB>2</SUB> computes the essential term <I>bh</I>, <I>P</I><SUB>3</SUB> computes the essential term <I>ce</I>, and <I>P</I><SUB>4</SUB> computes the essential term <I>df</I>. Thus, it remains for us to compute the remaining two submatrices <I>r</I> and <I>u</I>, whose essential terms are the diagonal terms <I>ae</I>, <I>bf</I>, <I>cg</I>, and <I>dh</I>, without using more than 3 additional products. We now try the innovation <I>P</I><SUB>5</SUB> in order to compute two essential terms at once:<P>
<img src="742_b.gif"><P>
In addition to computing both of the essential terms <I>ae</I> and <I>dh</I>, <I>P</I><SUB>5 </SUB>computes the inessential terms <I>ah</I> and <I>de</I>, which need to be cancelled somehow. We can use <I>P</I><SUB>4</SUB> and <I>P</I><SUB>2</SUB> to cancel them, but two other inessential terms then appear:<P>
<img src="742_c.gif"><P>
By adding an additional product<P>
<img src="743_a.gif"><P>
however, we obtain<P>
<img src="743_b.gif"><P>
We can obtain <I>u</I> in a similar manner from <I>P</I><SUB>5</SUB> by using <I>P</I><SUB>1</SUB> and <I>P</I><SUB>3</SUB> to move the inessential terms of <I>P</I><SUB>5</SUB> in a different direction:<P>
<img src="743_c.gif"><P>
By subtracting an additional product<P>
<img src="743_d.gif"><P>
we now obtain<P>
<img src="743_e.gif"><P>
The 7 submatrix products <I>P</I><SUB>1</SUB>, <I>P</I><SUB>2</SUB>, . . . , <I>P</I><SUB>7</SUB> can thus be used to compute the product <I>C</I> = <I>AB</I>, which completes the description of Strassen's method.<P>
<P>







<h2>Discussion</h2><P>
The large constant hidden in the running time of Strassen<FONT FACE="CG Times (W1)" SIZE=2>'</FONT>s algorithm makes it impractical unless the matrices are large (<I>n</I> at least 45 or so) and dense (few zero entries). For small matrices, the straightforward algorithm is preferable, and for large, sparse matrices, there are special sparsematrix algorithms that beat Strassen's in practice. Thus, Strassen<FONT FACE="CG Times (W1)" SIZE=2>'</FONT>s method is largely of theoretical interest.<P>
By using advanced techniques beyond the scope of this text, one can in fact multiply <I>n</I> <img src="../images/mult.gif"> <I>n</I> matrices in better than <img src="../images/bound.gif">(<I>n</I><SUP>1g 7</SUP>) time. The current best upper bound is approximately <I>O</I>(<I>n</I><SUP>2.376</SUP>). The best lower bound known is just the obvious <img src="../images/omega12.gif">(<I>n</I><SUP>2</SUP>) bound (obvious because we have to fill in <I>n</I><SUP>2 </SUP>elements of the product matrix). Thus, we currently do not know how hard matrix multiplication really is.<P>
Strassen<FONT FACE="CG Times (W1)" SIZE=2>'</FONT>s algorithm does not require that the matrix entries be real numbers. All that matters is that the number system form an algebraic ring. If the matrix entries do not form a ring, however, sometimes other techniques can be brought to bear to allow his method to apply. These issues are discussed more fully in the next section.<P>
<P>







<h2><a name="0973_1a8f">Exercises<a name="0973_1a8f"></h2><P>
<a name="0973_1a90">31.2-1<a name="0973_1a90"><P>
Use Strassen<FONT FACE="CG Times (W1)" SIZE=2>'</FONT>s algorithm to compute the matrix product<P>
<img src="744_a.gif"><P>
Show your work.<P>