Computational complexity of matrix multiplication

If A, B are {{math|n × n}} matrices over a field, then their product AB is also an {{math|n × n}} matrix over that field, defined entrywise as

:$$

(AB)_{ij} = \sum_{k = 1}^n A_{ik} B_{kj}.

{{For|implementation techniques (in particular parallel and distributed algorithms)|Matrix multiplication algorithm}}

The simplest approach to computing the product of two {{math|n × n}} matrices A and B is to compute the arithmetic expressions coming from the definition of matrix multiplication. In pseudocode:

input A and B, both n by n matrices

initialize C to be an n by n matrix of all zeros

for i from 1 to n:

for j from 1 to n:

for k from 1 to n:

C[i][j] = C[i][j] + A[i][k]*B[k][j]

output C (as A*B)

This algorithm requires, in the worst case, {{tmath|n^3}} multiplications of scalars and {{tmath|n^3 - n^2}} additions for computing the product of two square {{math|n×n}} matrices. Its computational complexity is therefore {{tmath|O(n^3)}}, in a model of computation where field operations (addition and multiplication) take constant time (in practice, this is the case for floating point numbers, but not necessarily for integers).

{{Main|Strassen algorithm}}

Strassen's algorithm improves on naive matrix multiplication through a divide-and-conquer approach. The key observation is that multiplying two {{math|2 × 2}} matrices can be done with only 7 multiplications, instead of the usual 8 (at the expense of 11 additional addition and subtraction operations). This means that, treating the input {{math|n×n}} matrices as block {{math|2 × 2}} matrices, the task of multiplying {{math|n×n}} matrices can be reduced to 7 subproblems of multiplying {{math|n/2×n/2}} matrices. Applying this recursively gives an algorithm needing $O(\; n^\{\backslash log\_\{2\}7\})\; \backslash approx\; O(n^\{2.807\})$ field operations.

Unlike algorithms with faster asymptotic complexity, Strassen's algorithm is used in practice. The numerical stability is reduced compared to the naive algorithm,{{cite journal | last1=Miller | first1=Webb | title=Computational complexity and numerical stability | citeseerx = 10.1.1.148.9947 | year=1975 | journal=SIAM News | volume=4 | issue=2 | pages=97–107 | doi=10.1137/0204009}} but it is faster in cases where {{math|n > 100}} or so{{cite book |first=Steven |last=Skiena |date=2012 |author-link=Steven Skiena |title=The Algorithm Design Manual |url=https://archive.org/details/algorithmdesignm00skie_772 |url-access=limited |publisher=Springer |pages=[https://archive.org/details/algorithmdesignm00skie_772/page/n56 45]–46, 401–403 |doi=10.1007/978-1-84800-070-4_4|chapter=Sorting and Searching |isbn=978-1-84800-069-8 }} and appears in several libraries, such as BLAS.{{cite book |last1=Press |first1=William H. |last2=Flannery |first2=Brian P. |last3=Teukolsky |first3=Saul A. |author3-link=Saul Teukolsky |last4=Vetterling |first4=William T. |title=Numerical Recipes: The Art of Scientific Computing |publisher=Cambridge University Press |edition=3rd |isbn=978-0-521-88068-8 |year=2007 |page=[https://archive.org/details/numericalrecipes00pres_033/page/n131 108]|title-link=Numerical Recipes }} Fast matrix multiplication algorithms cannot achieve component-wise stability, but some can be shown to exhibit norm-wise stability.{{cite journal | last1=Ballard | first1=Grey | last2=Benson | first2=Austin R. | last3=Druinsky | first3=Alex | last4=Lipshitz | first4=Benjamin | last5=Schwartz | first5=Oded | title=Improving the numerical stability of fast matrix multiplication | year=2016 | journal=SIAM Journal on Matrix Analysis and Applications | volume=37 | issue=4 | pages=1382–1418 | doi=10.1137/15M1032168 | arxiv=1507.00687| s2cid=2853388 }} It is very useful for large matrices over exact domains such as finite fields, where numerical stability is not an issue.

File:MatrixMultComplexity svg.svg

File:MatrixMultComplexity1990 svg.svg

The matrix multiplication exponent, usually denoted {{math|ω}}, is the smallest real number for which any two $n\backslash times\; n$ matrices over a field can be multiplied together using $n^\{\backslash omega\; +\; o(1)\}$ field operations. This notation is commonly used in algorithms research, so that algorithms using matrix multiplication as a subroutine have bounds on running time that can update as bounds on {{math|ω}} improve.

Using a naive lower bound and schoolbook matrix multiplication for the upper bound, one can straightforwardly conclude that {{math|2 ≤ ω ≤ 3}}. Whether {{math|1=ω = 2}} is a major open question in theoretical computer science, and there is a line of research developing matrix multiplication algorithms to get improved bounds on {{math|ω}}.

All recent algorithms in this line of research use the laser method, a generalization of the Coppersmith–Winograd algorithm, which was given by Don Coppersmith and Shmuel Winograd in 1990 and was the best matrix multiplication algorithm until 2010.{{cite journal|doi=10.1016/S0747-7171(08)80013-2 |title=Matrix multiplication via arithmetic progressions |url=http://www.cs.umd.edu/~gasarch/TOPICS/ramsey/matrixmult.pdf |year=1990 |last1=Coppersmith |first1=Don |last2=Winograd |first2=Shmuel |journal=Journal of Symbolic Computation |volume=9|issue=3|pages=251|doi-access=free }} The conceptual idea of these algorithms is similar to Strassen's algorithm: a method is devised for multiplying two {{math|k × k}}-matrices with fewer than {{math|k³}} multiplications, and this technique is applied recursively. The laser method has limitations to its power: Ambainis, Filmus and Le Gall prove that it cannot be used to show that {{math|ω < 2.3725}} by analyzing higher and higher tensor powers of a certain identity of Coppersmith and Winograd and neither {{math|ω < 2.3078}} for a wide class of variants of this approach.{{Cite book |last1=Ambainis |first1=Andris |title=Proceedings of the forty-seventh annual ACM symposium on Theory of Computing |last2=Filmus |first2=Yuval |last3=Le Gall |first3=François |date=2015-06-14 |publisher=Association for Computing Machinery |isbn=978-1-4503-3536-2 |series=STOC '15 |location=Portland, Oregon, USA |pages=585–593 |chapter=Fast Matrix Multiplication |doi=10.1145/2746539.2746554 |chapter-url=https://doi.org/10.1145/2746539.2746554 |arxiv=1411.5414 |s2cid=8332797}} In 2022 Duan, Wu and Zhou devised a variant breaking the first of the two barriers with {{math|ω < 2.37188}}, they do so by identifying a source of potential optimization in the laser method termed combination loss for which they compensate using an asymmetric version of the hashing method in the Coppersmith–Winograd algorithm.

Nonetheless, the above are classical examples of galactic algorithms. On the opposite, the above Strassen's algorithm of 1969 and Pan's algorithm of 1978, whose respective exponents are slightly above and below 2.78, have constant coefficients that make them feasible.{{cite journal | last1=Laderman | first1=Julian | last2=Pan | first2=Victor |last3=Sha | first3=Xuan-He | title=On practical algorithms for accelerated matrix multiplication | year=1992 | journal=Linear Algebra and Its Applications | volume=162-164 | pages=557–588 | doi=10.1016/0024-3795(92)90393-O}}

Henry Cohn, Robert Kleinberg, Balázs Szegedy and Chris Umans put methods such as the Strassen and Coppersmith–Winograd algorithms in an entirely different group-theoretic context, by utilising triples of subsets of finite groups which satisfy a disjointness property called the triple product property (TPP). They also give conjectures that, if true, would imply that there are matrix multiplication algorithms with essentially quadratic complexity. This implies that the optimal exponent of matrix multiplication is 2, which most researchers believe is indeed the case. One such conjecture is that families of wreath products of Abelian groups with symmetric groups realise families of subset triples with a simultaneous version of the TPP.{{Cite book | last1 = Cohn | first1 = H. | last2 = Kleinberg | first2 = R. | last3 = Szegedy | first3 = B. | last4 = Umans | first4 = C. | chapter = Group-theoretic Algorithms for Matrix Multiplication | doi = 10.1109/SFCS.2005.39 | title = 46th Annual IEEE Symposium on Foundations of Computer Science (FOCS'05) | pages = 379 | year = 2005 | isbn = 0-7695-2468-0 | s2cid = 41278294 | url = https://authors.library.caltech.edu/23966/ }}{{cite book |first1=Henry |last1=Cohn |first2=Chris |last2=Umans |chapter=A Group-theoretic Approach to Fast Matrix Multiplication |arxiv=math.GR/0307321 |title=Proceedings of the 44th Annual IEEE Symposium on Foundations of Computer Science, 11–14 October 2003 |year=2003 |publisher=IEEE Computer Society |pages=438–449 |doi=10.1109/SFCS.2003.1238217 |isbn=0-7695-2040-5 |s2cid=5890100 }} Several of their conjectures have since been disproven by Blasiak, Cohn, Church, Grochow, Naslund, Sawin, and Umans using the Slice Rank method.{{Cite book | last1 = Blasiak | first1 = J. | last2 = Cohn | first2 = H. | last3 = Church | first3 = T. | last4 = Grochow | first4 = J. | last5 = Naslund | first5= E. | last6 = Sawin | first6 = W. | last7=Umans | first7= C.| chapter= On cap sets and the group-theoretic approach to matrix multiplication | doi = 10.19086/da.1245 | title = Discrete Analysis | year = 2017 | page = 1245 | s2cid = 9687868 | url = http://discreteanalysisjournal.com/article/1245-on-cap-sets-and-the-group-theoretic-approach-to-matrix-multiplication}} Further, Alon, Shpilka and Chris Umans have recently shown that some of these conjectures implying fast matrix multiplication are incompatible with another plausible conjecture, the sunflower conjecture,{{cite journal |journal=Electronic Colloquium on Computational Complexity |date=April 2011 |author1-link=Noga Alon |last1=Alon |first1=N. |last2=Shpilka |first2=A. |last3=Umans |first3=C. |url=http://eccc.hpi-web.de/report/2011/067/ |title=On Sunflowers and Matrix Multiplication |id=TR11-067 }} which in turn is related to the cap set problem.

There is a trivial lower bound of {{tmath|\omega \ge 2}}. Since any algorithm for multiplying two {{math|n × n}}-matrices has to process all {{math|2n²}} entries, there is a trivial asymptotic lower bound of {{math|Ω(n²)}} operations for any matrix multiplication algorithm. Thus {{tmath|2\le \omega < 2.37188}}. It is unknown whether {{tmath|\omega > 2}}. The best known lower bound for matrix-multiplication complexity is {{math|Ω(n² log(n))}}, for bounded coefficient arithmetic circuits over the real or complex numbers, and is due to Ran Raz.{{cite book | last1 = Raz | first1 = Ran | title = Proceedings of the thiry-fourth annual ACM symposium on Theory of computing | chapter = On the complexity of matrix product | author-link = Ran Raz | year = 2002| pages = 144–151 | doi = 10.1145/509907.509932 | isbn = 1581134959 | s2cid = 9582328 }}

The exponent ω is defined to be a limit point, in that it is the infimum of the exponent over all matrix multiplication algorithms. It is known that this limit point is not achieved. In other words, under the model of computation typically studied, there is no matrix multiplication algorithm that uses precisely {{math|O(n^ω)}} operations; there must be an additional factor of {{math|n^o(1)}}.

Similar techniques also apply to rectangular matrix multiplication. The central object of study is $\backslash omega(k)$, which is the smallest $c$ such that one can multiply a matrix of size $n\backslash times\; \backslash lceil\; n^k\backslash rceil$ with a matrix of size $\backslash lceil\; n^k\backslash rceil\; \backslash times\; n$ with $O(n^\{c\; +\; o(1)\})$ arithmetic operations. A result in algebraic complexity states that multiplying matrices of size $n\backslash times\; \backslash lceil\; n^k\backslash rceil$ and $\backslash lceil\; n^k\backslash rceil\; \backslash times\; n$ requires the same number of arithmetic operations as multiplying matrices of size $n\backslash times\; \backslash lceil\; n^k\backslash rceil$ and $n\; \backslash times\; n$ and of size $n\; \backslash times\; n$ and $n\backslash times\; \backslash lceil\; n^k\backslash rceil$, so this encompasses the complexity of rectangular matrix multiplication.{{cite conference

| last1 = Gall | first1 = Francois Le

| last2 = Urrutia | first2 = Florent

| editor-last = Czumaj | editor-first = Artur

| arxiv = 1708.05622

| contribution = Improved rectangular matrix multiplication using powers of the Coppersmith-Winograd tensor

| doi = 10.1137/1.9781611975031.67

| pages = 1029–1046

| publisher = Society for Industrial and Applied Mathematics

| title = Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2018, New Orleans, LA, USA, January 7–10, 2018

| year = 2018| isbn = 978-1-61197-503-1

}} This generalizes the square matrix multiplication exponent, since $\backslash omega(1)\; =\; \backslash omega$.

Since the output of the matrix multiplication problem is size $n^2$, we have $\backslash omega(k)\; \backslash geq\; 2$ for all values of $k$. If one can prove for some values of $k$ between 0 and 1 that $\backslash omega(k)\; \backslash leq\; 2$, then such a result shows that $\backslash omega(k)\; =\; 2$ for those $k$. The largest k such that $\backslash omega(k)\; =\; 2$ is known as the dual matrix multiplication exponent, usually denoted α. α is referred to as the "dual" because showing that $\backslash alpha\; =\; 1$ is equivalent to showing that $\backslash omega\; =\; 2$. Like the matrix multiplication exponent, the dual matrix multiplication exponent sometimes appears in the complexity of algorithms in numerical linear algebra and optimization.{{Cite journal|last1=Cohen|first1=Michael B.|last2=Lee|first2=Yin Tat|last3=Song|first3=Zhao|date=2021-01-05|title=Solving Linear Programs in the Current Matrix Multiplication Time|url=https://doi.org/10.1145/3424305|journal=Journal of the ACM|volume=68|issue=1|pages=3:1–3:39|doi=10.1145/3424305|issn=0004-5411|arxiv=1810.07896|s2cid=231955576 }}

The first bound on α is by Coppersmith in 1982, who showed that $\backslash alpha\; >\; 0.17227$.{{Cite journal|last=Coppersmith|first=D.|date=1982-08-01|title=Rapid Multiplication of Rectangular Matrices|url=https://epubs.siam.org/doi/10.1137/0211037|journal=SIAM Journal on Computing|volume=11|issue=3|pages=467–471|doi=10.1137/0211037|issn=0097-5397|url-access=subscription}} The current best peer-reviewed bound on α is $\backslash alpha\; \backslash geq\; 0.321334$, given by Williams, Xu, Xu, and Zhou.

{{further|Computational complexity of mathematical operations#Matrix algebra}}

Problems that have the same asymptotic complexity as matrix multiplication include determinant, matrix inversion, Gaussian elimination (see next section). Problems with complexity that is expressible in terms of $\backslash omega$ include characteristic polynomial, eigenvalues (but not eigenvectors), Hermite normal form, and Smith normal form.{{citation needed|date=March 2018}}

In his 1969 paper, where he proved the complexity $O(n^\{\backslash log\_2\; 7\})\; \backslash approx\; O(n^\{2.807\})$ for matrix computation, Strassen proved also that matrix inversion, determinant and Gaussian elimination have, up to a multiplicative constant, the same computational complexity as matrix multiplication. The proof does not make any assumptions on matrix multiplication that is used, except that its complexity is $O(n^\backslash omega)$ for some $\backslash omega\; \backslash ge\; 2$

The starting point of Strassen's proof is using block matrix multiplication. Specifically, a matrix of even dimension {{math|2n×2n}} may be partitioned in four {{math|n×n}} blocks

:

Under this form, its inverse is

:$$

\begin{bmatrix} {A} & {B} \\ {C} & {D} \end{bmatrix}^{-1} =

\begin{bmatrix}

{A}^{-1}+{A}^{-1}{B}({D}-{CA}^{-1}{B})^{-1}{CA}^{-1} & -{A}^{-1}{B}({D}-{CA}^{-1}{B})^{-1}

\\ -({D}-{CA}^{-1}{B})^{-1}{CA}^{-1} & ({D}-{CA}^{-1}{B})^{-1}

\end{bmatrix},

provided that {{mvar|A}} and $\{D\}-\{CA\}^\{-1\}\{B\}$ are invertible.

Thus, the inverse of a {{math|2n×2n}} matrix may be computed with two inversions, six multiplications and four additions or additive inverses of {{math|n×n}} matrices. It follows that, denoting respectively by {{math|I(n)}}, {{math|M(n)}} and {{math|1= A(n) = n²}} the number of operations needed for inverting, multiplying and adding {{math|n×n}} matrices, one has

:$I(2n)\; \backslash le\; 2I(n)\; +\; 6M(n)+\; 4\; A(n).$

If $n=2^k,$ one may apply this formula recursively:

:$\backslash begin\{align\}$

I(2^k) &\le 2I(2^{k-1}) + 6M(2^{k-1})+ 4 A(2^{k-1})\\

&\le 2^2I(2^{k-2}) + 6(M(2^{k-1})+2M(2^{k-2})) + 4(A(2^{k-1}) + 2A(2^{k-2}))\\

&\,\,\,\vdots

\end{align}

If $M(n)\backslash le\; cn^\backslash omega,$ and $\backslash alpha=2^\backslash omega\backslash ge\; 4,$ one gets eventually

:$\backslash begin\{align\}$

I(2^k) &\le 2^k I(1) + 6c(\alpha^{k-1}+2\alpha^{k-2} + \cdots +2^{k-1}\alpha^0) + k 2^{k+1}\\

&\le 2^k + 6c\frac{\alpha^k-2^k}{\alpha-2} + k 2^{k+1}\\

&\le d(2^k)^\omega

\end{align}

for some constant {{mvar|d}}.

For matrices whose dimension is not a power of two, the same complexity is reached by increasing the dimension of the matrix to a power of two, by padding the matrix with rows and columns whose entries are 1 on the diagonal and 0 elsewhere.

This proves the asserted complexity for matrices such that all submatrices that have to be inverted are indeed invertible. This complexity is thus proved for almost all matrices, as a matrix with randomly chosen entries is invertible with probability one.

The same argument applies to LU decomposition, as, if the matrix {{mvar|A}} is invertible, the equality

:

= \begin{bmatrix}I & 0\\CA^{-1}&I\end{bmatrix}\,\begin{bmatrix}A&B\\0&D-CA^{-1}B\end{bmatrix}

defines a block LU decomposition that may be applied recursively to $A$ and $D-CA^\{-1\}B,$ for getting eventually a true LU decomposition of the original matrix.

The argument applies also for the determinant, since it results from the block LU decomposition that

:

\det(A)\det(D-CA^{-1}B).

Related to the problem of minimizing the number of arithmetic operations is minimizing the number of multiplications, which is typically a more costly operation than addition. A $O(n^\backslash omega)$ algorithm for matrix multiplication must necessarily only use $O(n^\backslash omega)$ multiplication operations, but these algorithms are impractical. Improving from the naive $n^3$ multiplications for schoolbook multiplication, $4\backslash times\; 4$ matrices in $\backslash mathbb\{Z\}/2\backslash mathbb\{Z\}$ can be done with 47 multiplications,See [https://www.nature.com/articles/s41586-022-05172-4/figures/6 Extended Data Fig. 1: Algorithm for multiplying 4 × 4 matrices in modular arithmetic ($\backslash mathbb\{Z\}\_\{2\}$)) with 47 multiplications] in {{Cite journal |title=Discovering faster matrix multiplication algorithms with reinforcement learning | year=2022 |language=en |doi=10.1038/s41586-022-05172-4| pmid=36198780 | last1=Fawzi | first1=A. | last2=Balog | first2=M. | last3=Huang | first3=A. | last4=Hubert | first4=T. | last5=Romera-Paredes | first5=B. | last6=Barekatain | first6=M. | last7=Novikov | first7=A. | last8=r Ruiz | first8=F. J. | last9=Schrittwieser | first9=J. | last10=Swirszcz | first10=G. | last11=Silver | first11=D. | last12=Hassabis | first12=D. | last13=Kohli | first13=P. | journal=Nature | volume=610 | issue=7930 | pages=47–53 | pmc=9534758 | bibcode=2022Natur.610...47F }} $3\backslash times\; 3$ matrix multiplication over a commutative ring can be done in 21 multiplications{{cite journal

| last = Rosowski | first = Andreas

| arxiv = 1904.07683

| doi = 10.1016/j.jsc.2022.05.002

| journal = Journal of Symbolic Computation

| mr = 4433063

| pages = 302–321

| title = Fast commutative matrix algorithms

| volume = 114

| year = 2023}}{{cite journal |last1=Makarov |first1=O. M. |title=An algorithm for multiplying 3×3 matrices |journal=Zhurnal Vychislitel'noi Matematiki I Matematicheskoi Fiziki |volume=26 |issue=2 |year=1986 |pages=293–294 |url=https://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=zvmmf&paperid=4056&option_lang=eng |access-date=5 October 2022}}

:Also in {{cite journal |doi=10.1016/0041-5553(86)90203-X |title=An algorithm for multiplying 3×3 matrices |year=1986 |last1=Makarov |first1=O. M. |journal=USSR Computational Mathematics and Mathematical Physics |volume=26 |pages=179–180 }} (23 if non-commutative{{Cite journal |last=Laderman |first=Julian D. |date=1976 |title=A noncommutative algorithm for multiplying 3×3 matrices using 23 multiplications |url=https://www.ams.org/bull/1976-82-01/S0002-9904-1976-13988-2/ |journal=Bulletin of the American Mathematical Society |language=en |volume=82 |issue=1 |pages=126–128 |doi=10.1090/S0002-9904-1976-13988-2 |issn=0002-9904|doi-access=free }}). The lower bound of multiplications needed is 2mn+2n−m−2 (multiplication of n×m-matrices with m×n-matrices using the substitution method, m⩾n⩾3), which means n=3 case requires at least 19 multiplications and n=4 at least 34.{{Cite journal |last=Bläser |first=Markus |date=February 2003 |title=On the complexity of the multiplication of matrices of small formats |journal=Journal of Complexity |language=en |volume=19 |issue=1 |pages=43–60 |doi=10.1016/S0885-064X(02)00007-9|doi-access=free }} For n=2 optimal 7 multiplications 15 additions are minimal, compared to only 4 additions for 8 multiplications.{{Cite journal |last=Winograd |first=S. |date=1971-10-01 |title=On multiplication of 2 × 2 matrices |journal=Linear Algebra and Its Applications |language=en |volume=4 |issue=4 |pages=381–388 |doi=10.1016/0024-3795(71)90009-7 |issn=0024-3795|doi-access=free }}{{Cite book |last=L. |first=Probert, R. |url=http://worldcat.org/oclc/1124200063 |title=On the complexity of matrix multiplication |date=1973 |publisher=University of Waterloo |oclc=1124200063}}

• Computational complexity of mathematical operations
• CYK algorithm, §Valiant's algorithm
• Freivalds' algorithm, a simple Monte Carlo algorithm that, given matrices {{mvar|A}}, {{mvar|B}} and {{mvar|C}}, verifies in {{math|Θ(n²)}} time if {{math|AB {{=}} C}}.
• Matrix chain multiplication
• Matrix multiplication, for abstract definitions
• Matrix multiplication algorithm, for practical implementation details
• Sparse matrix–vector multiplication

{{reflist|30em}}

• [https://fmm.univ-lille.fr/ Yet another catalogue of fast matrix multiplication algorithms]
• {{cite journal |last1=Fawzi |first1=A. |last2=Balog |first2=M. |last3=Huang |first3=A. |first4=T. |last4=Hubert |first5=B. |last5=Romera-Paredes |first6=M. |last6=Barekatain |first7=A. |last7=Novikov |first8=F.J.R. |last8=Ruiz |first9=J. |last9=Schrittwieser |first10=G. |last10=Swirszcz |first11=D. |last11=Silver |first12=D. |last12=Hassabis |first13=P. |last13=Kohli |title=Discovering faster matrix multiplication algorithms with reinforcement learning |journal=Nature |volume=610 |issue= 7930|pages=47–53 |date=2022 |doi=10.1038/s41586-022-05172-4 |pmid=36198780 |pmc=9534758 |bibcode=2022Natur.610...47F |url=}}

Category:Computer arithmetic algorithms

Category:Computational complexity theory

Category:Matrix theory

Category:Unsolved problems in computer science