HOME

TheInfoList



OR:

In linear algebra, the Strassen algorithm, named after Volker Strassen, is an
algorithm for matrix multiplication Because matrix multiplication is such a central operation in many numerical algorithms, much work has been invested in making matrix multiplication algorithms efficient. Applications of matrix multiplication in computational problems are found in ...
. It is faster than the standard matrix multiplication algorithm for large matrices, with a better asymptotic complexity, although the naive algorithm is often better for smaller matrices. The Strassen algorithm is slower than the fastest known algorithms for extremely large matrices, but such galactic algorithms are not useful in practice, as they are much slower for matrices of practical size. For small matrices even faster algorithms exist. Strassen's algorithm works for any
ring Ring may refer to: * Ring (jewellery), a round band, usually made of metal, worn as ornamental jewelry * To make a sound with a bell, and the sound made by a bell :(hence) to initiate a telephone connection Arts, entertainment and media Film and ...
, such as plus/multiply, but not all semirings, such as min-plus or
boolean algebra In mathematics and mathematical logic, Boolean algebra is a branch of algebra. It differs from elementary algebra in two ways. First, the values of the variables are the truth values ''true'' and ''false'', usually denoted 1 and 0, whereas i ...
, where the naive algorithm still works, and so called combinatorial matrix multiplication.


History

Volker Strassen first published this algorithm in 1969 and thereby proved that the n^3 general matrix multiplication algorithm wasn't optimal. The Strassen algorithm's publication resulted in more research about matrix multiplication that led to both asymptotically lower bounds and improved computational upper bounds.


Algorithm

Let A, B be two square matrices over a
ring Ring may refer to: * Ring (jewellery), a round band, usually made of metal, worn as ornamental jewelry * To make a sound with a bell, and the sound made by a bell :(hence) to initiate a telephone connection Arts, entertainment and media Film and ...
\mathcal, for example matrices whose entries are integers or the real numbers. The goal of matrix multiplication is to calculate the matrix product C = AB. The following exposition of the algorithm assumes that all of these matrices have sizes that are powers of two (i.e., A, \, B, \, C \in \operatorname_ (\mathcal)), but this is only conceptually necessary -- if the matrices A, B are not of type 2^n \times 2^n, the "missing" rows and columns can be filled with zeros to obtain matrices with sizes of powers of two -- though real implementations of the algorithm do not do this in practice. The Strassen algorithm partitions A, B and C into equally sized
block matrices In mathematics, a block matrix or a partitioned matrix is a matrix that is '' interpreted'' as having been broken into sections called blocks or submatrices. Intuitively, a matrix interpreted as a block matrix can be visualized as the original ma ...
: A = \begin A_ & A_ \\ A_ & A_ \end, \quad B = \begin B_ & B_ \\ B_ & B_ \end, \quad C = \begin C_ & C_ \\ C_ & C_ \end, \quad with A_, B_, C_ \in \operatorname_ (\mathcal). The naive algorithm would be: : \begin C_ & C_ \\ C_ & C_ \end = \begin A_ B_ + A_ B_ & A_ B_ + A_ B_ \\ A_ B_ + A_ B_ & A_ B_ + A_ B_ \end. This construction does not reduce the number of multiplications: 8 multiplications of matrix blocks are still needed to calculate the C_ matrices, the same number of multiplications needed when using standard matrix multiplication. The Strassen algorithm defines instead new matrices: : \begin M_1 &= (A_ + A_) (B_ + B_); \\ M_2 &= (A_ + A_) B_; \\ M_3 &= A_ (B_ - B_); \\ M_4 &= A_ (B_ - B_); \\ M_5 &= (A_ + A_) B_; \\ M_6 &= (A_ - A_) (B_ + B_); \\ M_7 &= (A_ - A_) (B_ + B_), \\ \end using only 7 multiplications (one for each M_k) instead of 8. We may now express the C_ in terms of M_k: : \begin C_ & C_ \\ C_ & C_ \end = \begin M_1 + M_4 - M_5 + M_7 & M_3 + M_5 \\ M_2 + M_4 & M_1 - M_2 + M_3 + M_6 \end. We recursively iterate this division process until the
submatrices In mathematics, a matrix (plural matrices) is a rectangular array or table of numbers, symbols, or expressions, arranged in rows and columns, which is used to represent a mathematical object or a property of such an object. For example, \beg ...
degenerate into numbers (elements of the ring \mathcal). If, as mentioned above, the original matrix had a size that was not a power of 2, then the resulting product will have zero rows and columns just like A and B, and these will then be stripped at this point to obtain the (smaller) matrix C we really wanted. Practical implementations of Strassen's algorithm switch to standard methods of matrix multiplication for small enough submatrices, for which those algorithms are more efficient. The particular crossover point for which Strassen's algorithm is more efficient depends on the specific implementation and hardware. Earlier authors had estimated that Strassen's algorithm is faster for matrices with widths from 32 to 128 for optimized implementations. However, it has been observed that this crossover point has been increasing in recent years, and a 2010 study found that even a single step of Strassen's algorithm is often not beneficial on current architectures, compared to a highly optimized traditional multiplication, until matrix sizes exceed 1000 or more, and even for matrix sizes of several thousand the benefit is typically marginal at best (around 10% or less). A more recent study (2016) observed benefits for matrices as small as 512 and a benefit around 20%.


Winograd form

It is possible to reduce the number of matrix additions by instead using the following form discovered by Winograd: \begin a & b \\ c & d \end \begin A & C \\ B & D \end = \begin aA + bB & w + v + (a + b - c - d)D \\ w + u + d(B + C - A - D) & w + u + v \end where u = (c - a)(C - D), v = (c + d)(C - A), w = aA + (c + d - a)(A + D - C). This reduces the number of matrix additions and subtractions from 18 to 15. The number of matrix multiplications is still 7, and the asymptotic complexity is the same.


Asymptotic complexity

The outline of the algorithm above showed that one can get away with just 7, instead of the traditional 8, matrix-matrix multiplications for the sub-blocks of the matrix. On the other hand, one has to do additions and subtractions of blocks, though this is of no concern for the overall complexity: Adding matrices of size N/2 requires only (N/2)^2 operations whereas multiplication is substantially more expensive (traditionally 2 (N/2)^3 addition or multiplication operations). The question then is how many operations exactly one needs for Strassen's algorithms, and how this compares with the standard matrix multiplication that takes approximately 2 N^3 (where N = 2^n) arithmetic operations, i.e. an asymptotic complexity \Theta (N^3). The number of additions and multiplications required in the Strassen algorithm can be calculated as follows: let f(n) be the number of operations for a 2^n \times 2^n matrix. Then by recursive application of the Strassen algorithm, we see that f(n) = 7 f(n-1) + l 4^n, for some constant l that depends on the number of additions performed at each application of the algorithm. Hence f(n) = (7 + o(1))^n, i.e., the asymptotic complexity for multiplying matrices of size N = 2^n using the Strassen algorithm is O( +o(1)n) = O(N^) \approx O(N^). The reduction in the number of arithmetic operations however comes at the price of a somewhat reduced
numerical stability In the mathematical subfield of numerical analysis, numerical stability is a generally desirable property of numerical algorithms. The precise definition of stability depends on the context. One is numerical linear algebra and the other is algorit ...
, and the algorithm also requires significantly more memory compared to the naive algorithm. Both initial matrices must have their dimensions expanded to the next power of 2, which results in storing up to four times as many elements, and the seven auxiliary matrices each contain a quarter of the elements in the expanded ones. Strassen's algorithm needs to be compared to the "naive" way of doing the matrix multiplication that would require 8 instead of 7 multiplications of sub-blocks. This would then give rise to the complexity one expects from the standard approach: O(8^) = O(N^) = O(N^3). The comparison of these two algorithms shows that ''asymptotically'', Strassen's algorithm is faster: There exists a size N_\text so that matrices that are larger are more efficiently multiplied with Strassen's algorithm than the "traditional" way. However, the asymptotic statement does not imply that Strassen's algorithm is ''always'' faster even for small matrices, and in practice this is in fact not the case: For small matrices, the cost of the additional additions of matrix blocks outweighs the savings in the number of multiplications. There are also other factors not captured by the analysis above, such as the difference in cost on today's hardware between loading data from memory onto processors vs. the cost of actually doing operations on this data. As a consequence of these sorts of considerations, Strassen's algorithm is typically only used on "large" matrices. This kind of effect is even more pronounced with alternative algorithms such as the one by Coppersmith and Winograd: While ''asymptotically'' even faster, the cross-over point N_\text is so large that the algorithm is not generally used on matrices one encounters in practice.


Rank or bilinear complexity

The bilinear complexity or rank of a
bilinear map In mathematics, a bilinear map is a function combining elements of two vector spaces to yield an element of a third vector space, and is linear in each of its arguments. Matrix multiplication is an example. Definition Vector spaces Let V, ...
is an important concept in the asymptotic complexity of matrix multiplication. The rank of a bilinear map \phi:\mathbf A \times \mathbf B \rightarrow \mathbf C over a field F is defined as (somewhat of an
abuse of notation In mathematics, abuse of notation occurs when an author uses a mathematical notation in a way that is not entirely formally correct, but which might help simplify the exposition or suggest the correct intuition (while possibly minimizing errors a ...
) :R(\phi/\mathbf F) = \min \left\ In other words, the rank of a bilinear map is the length of its shortest bilinear computation. The existence of Strassen's algorithm shows that the rank of 2 \times 2 matrix multiplication is no more than seven. To see this, let us express this algorithm (alongside the standard algorithm) as such a bilinear computation. In the case of matrices, the dual spaces A* and B* consist of maps into the field F induced by a scalar double-dot product, (i.e. in this case the sum of all the entries of a Hadamard product.) It can be shown that the total number of elementary multiplications L required for matrix multiplication is tightly asymptotically bound to the rank R, i.e. L = \Theta(R), or more specifically, since the constants are known, R / 2 \le L \le R. One useful property of the rank is that it is submultiplicative for tensor products, and this enables one to show that 2^n \times 2^n \times 2^n matrix multiplication can be accomplished with no more than 7n elementary multiplications for any n. (This n-fold tensor product of the 2 \times 2 \times 2 matrix multiplication map with itself — an n-th tensor power—is realized by the recursive step in the algorithm shown.)


Cache behavior

Strassen's algorithm is cache oblivious. Analysis of its
cache Cache, caching, or caché may refer to: Places United States * Cache, Idaho, an unincorporated community * Cache, Illinois, an unincorporated community * Cache, Oklahoma, a city in Comanche County * Cache, Utah, Cache County, Utah * Cache Count ...
behavior algorithm has shown it to incur :\Theta \left(1 + \frac + \frac \right) cache misses during its execution, assuming an idealized cache of size M (i.e. with M / b lines of length b).


Implementation considerations

The description above states that the matrices are square, and the size is a power of two, and that padding should be used if needed. This restriction allows the matrices to be split in half, recursively, until limit of scalar multiplication is reached. The restriction simplifies the explanation, and analysis of complexity, but is not actually necessary; and in fact, padding the matrix as described will increase the computation time and can easily eliminate the fairly narrow time savings obtained by using the method in the first place. A good implementation will observe the following: * It is not necessary or desirable to use the Strassen algorithm down to the limit of scalars. Compared to conventional matrix multiplication, the algorithm adds a considerable O(n^) workload in addition/subtractions; so below a certain size, it will be better to use conventional multiplication. Thus, for instance, a 1600 \times 1600 does not need to be padded to 2048 \times 2048, since it could be subdivided down to 25 \times 25 matrices and conventional multiplication can then be used at that level. * The method can indeed be applied to square matrices of any dimension. If the dimension is even, they are split in half as described. If the dimension is odd, zero padding by one row and one column is applied first. Such padding can be applied on-the-fly and lazily, and the extra rows and columns discarded as the result is formed. For instance, suppose the matrices are 199 \times 199. They can be split so that the upper-left portion is 100 \times 100 and the lower-right is 99 \times 99. Wherever the operations require it, dimensions of 99 are zero padded to 100 first. Note, for instance, that the product M_2 is only used in the lower row of the output, so is only required to be 99 rows high; and thus the left factor A_ + A_ used to generate it need only be 99 rows high; accordingly, there is no need to pad that sum to 100 rows; it is only necessary to pad A_ to 100 columns to match A_. Furthermore, there is no need for the matrices to be square. Non-square matrices can be split in half using the same methods, yielding smaller non-square matrices. If the matrices are sufficiently non-square it will be worthwhile reducing the initial operation to more square products, using simple methods which are essentially O(n^), for instance: * A product of size N \times N\ast \times 10N/math> can be done as 20 separate \times N\ast \times N/math> operations, arranged to form the result; * A product of size \times 10N\ast 0N \times N/math> can be done as 10 separate \times N\ast \times N/math> operations, summed to form the result. These techniques will make the implementation more complicated, compared to simply padding to a power-of-two square; however, it is a reasonable assumption that anyone undertaking an implementation of Strassen, rather than conventional multiplication, will place a higher priority on computational efficiency than on simplicity of the implementation. In practice, Strassen's algorithm can be implemented to attain better performance than conventional multiplication even for small matrices, for matrices that are not at all square, and without requiring workspace beyond buffers that are already needed for a high-performance conventional multiplication.


See also

*
Computational complexity of mathematical operations The following tables list the computational complexity of various algorithms for common mathematical operations. Here, complexity refers to the time complexity of performing computations on a multitape Turing machine. See big O notation for an ...
*
Gauss–Jordan elimination In mathematics, Gaussian elimination, also known as row reduction, is an algorithm for solving systems of linear equations. It consists of a sequence of operations performed on the corresponding matrix of coefficients. This method can also be used ...
*
Coppersmith–Winograd algorithm In theoretical computer science, the computational complexity of matrix multiplication dictates how quickly the operation of matrix multiplication can be performed. Matrix multiplication algorithms are a central subroutine in theoretical and num ...
* Z-order matrix representation * Karatsuba algorithm, for multiplying ''n''-digit integers in O(n^) instead of in O(n^2) time * Toom-Cook algorithm, a faster generalization of the Karatsuba algorithm that permits recursive divide-and-conquer decomposition into more than 2 blocks at a time * Gauss's complex multiplication algorithm multiplies two complex numbers using 3 real multiplications instead of 4


References

*
Thomas H. Cormen Thomas H. Cormen is the co-author of ''Introduction to Algorithms'', along with Charles Leiserson, Ron Rivest, and Cliff Stein. In 2013, he published a new book titled '' Algorithms Unlocked''. He is a professor of computer science at Dartmou ...
,
Charles E. Leiserson Charles Eric Leiserson is a computer scientist, specializing in the theory of parallel computing and distributed computing, and particularly practical applications thereof. As part of this effort, he developed the Cilk multithreaded language. H ...
,
Ronald L. Rivest Ronald Linn Rivest (; born May 6, 1947) is a cryptographer and an Institute Professor at MIT. He is a member of MIT's Department of Electrical Engineering and Computer Science (EECS) and a member of MIT's Computer Science and Artificial Int ...
, and
Clifford Stein Clifford Seth Stein (born December 14, 1965), a computer scientist, is a professor of industrial engineering and operations research at Columbia University in New York, NY, where he also holds an appointment in the Department of Computer Scie ...
. '' Introduction to Algorithms'', Second Edition. MIT Press and McGraw-Hill, 2001. . Chapter 28: Section 28.2: Strassen's algorithm for matrix multiplication, pp. 735–741. *


External links

* (also includes formulas for fast matrix inversion) *Tyler J. Earnest,
Strassen's Algorithm on the Cell Broadband Engine
' {{Numerical linear algebra Matrix multiplication algorithms