User:Fawly/Computational complexity of 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 many fields including scientific computing and pattern recognition and in seemingly unrelated problems such as counting the paths through a graph. Many different algorithms have been designed for multiplying matrices on different types of hardware, including parallel and distributed systems, where the computational work is spread over multiple processors (perhaps over a network).

Directly applying the mathematical definition of matrix multiplication gives an algorithm that takes time on the order of n³ field operations to multiply two n × n matrices over that field (Θ(n³) in big O notation). Better asymptotic bounds on the time required to multiply matrices have been known since the Strassen's algorithm in the 1960s, but it is still unknown what the optimal time is (i.e., what the complexity of the problem is). As of December 2020, the matrix multiplication algorithm with best asymptotic complexity runs in O(n^2.3728596) time, given by Josh Alman and Virginia Vassilevska Williams, however this algorithm is a galactic algorithm because of the large constants and cannot be realized practically.

Algorithms exist that provide better running times than the straightforward ones. The first to be discovered was Strassen's algorithm, devised by Volker Strassen in 1969 and often referred to as "fast matrix multiplication". It is based on a way of multiplying two 2 × 2-matrices which requires only 7 multiplications (instead of the usual 8), at the expense of several additional addition and subtraction operations. Applying this recursively gives an algorithm with a multiplicative cost of $O(n^{\log _{2}7})\approx O(n^{2.807})$. Strassen's algorithm is more complex, and the numerical stability is reduced compared to the naïve algorithm, but it is faster in cases where n > 100 or so and appears in several libraries, such as BLAS. It is very useful for large matrices over exact domains such as finite fields, where numerical stability is not an issue.

The matrix multiplication algorithm that results of the definition requires, in the worst case, ⁠$n^{3}$⁠ multiplications of scalars and ⁠$(n-1)n^{2}$⁠ additions for computing the product of two square n×n matrices. Its computational complexity is therefore ⁠$O(n^{3})$⁠, in a model of computation for which the scalar operations require a constant time (in practice, this is the case for floating point numbers, but not for integers).

Rather surprisingly, this complexity is not optimal, as shown in 1969 by Volker Strassen, who provided an algorithm, now called Strassen's algorithm, with a complexity of $O(n^{\log _{2}7})\approx O(n^{2.8074}).$ The exponent appearing in the complexity of matrix multiplication has been improved several times, leading to the Coppersmith–Winograd algorithm with a complexity of O(n^2.3755) (1990). This algorithm has been slightly improved in 2010 by Stothers to a complexity of O(n^2.3737), in 2013 by Virginia Vassilevska Williams to O(n^2.3729), and in 2014 by François Le Gall to O(n^2.3728639). This was further refined in 2020 by Josh Alman and Virginia Vassilevska Williams to a final (up to date) complexity of O(n^2.3728596).

The greatest lower bound for the exponent of matrix multiplication algorithm is generally called ⁠$\omega$⁠. It is trivially true that ⁠$\omega \geq 2$⁠, because one must read the ⁠$n^{2}$⁠ elements of a matrix in order to multiply it with another matrix. Thus ⁠$2\leq \omega <2.373$⁠. It is unknown whether ⁠$\omega >2$⁠. The largest known lower bound for matrix-multiplication complexity is Ω(n² log(n)), for a restricted kind of arithmetic circuits, and is due to Ran Raz.

Improvement of estimates of exponent ω over time for the computational complexity of matrix multiplication $O(n^{\omega })$.

Freivalds' algorithm is a simple Monte Carlo algorithm that, given matrices A, B and C, verifies in Θ(n²) time if AB = C.

Matrix multiplication exponent

It is an open question in theoretical computer science how well Strassen's algorithm can be improved in terms of asymptotic complexity. The matrix multiplication exponent, usually denoted $\omega$, is the smallest real number for which any $n\times n$ matrix over a field can be multiplied together using $n^{\omega +o(1)}$ field operations. The current best bound on $\omega$ is $\omega <2.3728596$, by Josh Alman and Virginia Vassilevska Williams. This algorithm, like all other recent algorithms in this line of research, is a generalization of the Coppersmith–Winograd algorithm, which was given by Don Coppersmith and Shmuel Winograd in 1990, was the best matrix multiplication algorithm until 2010, and has an asymptotic complexity of O(n^2.375477). The conceptual idea of these algorithms are similar to Strassen's algorithm: a way is devised for multiplying two k × k-matrices with fewer than k³ multiplications, and this technique is applied recursively. However, the constant coefficient hidden by the Big O notation is so large that these algorithms are only worthwhile for matrices that are too large to handle on present-day computers.

Since any algorithm for multiplying two n × n-matrices has to process all 2n² entries, there is an asymptotic lower bound of Ω(n²) operations. Raz proved a lower bound of Ω(n² log(n)) for bounded coefficient arithmetic circuits over the real or complex numbers.

Group theory reformulation of matrix multiplication algorithms

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. Several of their conjectures have since been disproven by Blasiak, Cohn, Church, Grochow, Naslund, Sawin, and Umans using the Slice Rank method. 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.

Related complexities

The importance of the computational complexity of matrix multiplication relies on the facts that many algorithmic problems may be solved by means of matrix computation, and most problems on matrices have a complexity which is either the same as that of matrix multiplication (up to a multiplicative constant), or may be expressed in term of the complexity of matrix multiplication or its exponent $\omega .$

There are several advantages of expressing complexities in terms of the exponent $\omega$ of matrix multiplication. Firstly, if $\omega$ is improved, this will automatically improve the known upper bound of complexity of many algorithms. Secondly, in practical implementations, one never uses the matrix multiplication algorithm that has the best asymptotical complexity, because the constant hidden behind the big O notation is too large for making the algorithm competitive for sizes of matrices that can be manipulated in a computer. Thus expressing complexities in terms of $\omega$ provide a more realistic complexity, since it remains valid whichever algorithm is chosen for matrix computation.

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 $\omega$ include characteristic polynomial, eigenvalues (but not eigenvectors), Hermite normal form, and Smith normal form.

Matrix inversion, determinant and Gaussian elimination

In his 1969 paper, where he proved the complexity $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^{\omega })$ for some $\omega \geq 2$

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

${\begin{bmatrix}{A}&{B}\\{C}&{D}\end{bmatrix}}.$

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 A and ${D}-{CA}^{-1}{B}$ are invertible.

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

$I(2n)\leq 2I(n)+6M(n)+4A(n).$

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

${\begin{aligned}I(2^{k})&\leq 2I(2^{k-1})+6M(2^{k-1})+4A(2^{k-1})\\&\leq 2^{2}I(2^{k-2})+6(M(2^{k-1})+2M(2^{k-2}))+4(A(2^{k-1})+2A(2^{k-2}))\\&\ldots \end{aligned}}$

If $M(n)\leq cn^{\omega },$ and $\alpha =2^{\omega }\geq 4,$ one gets eventually

${\begin{aligned}I(2^{k})&\leq 2^{k}I(1)+6c(\alpha ^{k-1}+2\alpha ^{k-2}+\cdots +2^{k-1}\alpha ^{0})+k2^{k+1}\\&\leq 2^{k}+6c{\frac {\alpha ^{k}-2^{k}}{\alpha -2}}+k2^{k+1}\\&\leq d(2^{k})^{\omega }.\end{aligned}}$

for some constant 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 A is invertible, the equality

${\begin{bmatrix}{A}&{B}\\{C}&{D}\end{bmatrix}}={\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 {\begin{bmatrix}{A}&{B}\\{C}&{D}\end{bmatrix}}=\det(A)\det(D-CA^{-1}B).$

See also

• Computational complexity of mathematical operations
• CYK algorithm, §Valiant's algorithm
• Matrix chain multiplication
• Sparse matrix-vector multiplication