Becausematrix multiplication is such a central operation in manynumerical algorithms, much work has been invested in makingmatrix multiplication algorithms efficient. Applications of matrix multiplication in computational problems are found in many fields includingscientific computing andpattern recognition and in seemingly unrelated problems such as counting the paths through agraph.^[1] Many different algorithms have been designed for multiplying matrices on different types of hardware, includingparallel anddistributed 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 thattakes time on the order ofn³field operations to multiply twon ×n matrices over that field (Θ(n³) inbig O notation). Better asymptotic bounds on the time required to multiply matrices have been known since theStrassen's algorithm in the 1960s, but the optimal time (that is, thecomputational complexity of matrix multiplication) remains unknown. As of September 2025^[update], the best bound on theasymptotic complexity of a matrix multiplication algorithm isO(n^2.371339) time, given by Alman, Duan,Williams, Xu, Xu, and Zhou.^[2] However, this algorithm is agalactic algorithm because of the large constants and cannot be realized practically.

Thedefinition of matrix multiplication is that ifC =AB for ann ×m matrixA and anm ×p matrixB, thenC is ann ×p matrix with entries

${\displaystyle c_{ij}=\sum _{k=1}^m{a}_{ik}{b}_{kj}.}$

From this, a simple algorithm can be constructed which loops over the indicesi from 1 throughn andj from 1 throughp, computing the above using a nested loop:

This algorithm takestimeΘ(nmp) (inasymptotic notation).^[1] A common simplification for the purpose ofalgorithm analysis is to assume that the inputs are all square matrices of sizen ×n, in which case the running time isΘ(n³), i.e., cubic in the size of the dimension.^[3]

The three loops in iterative matrix multiplication can be arbitrarily swapped with each other without an effect on correctness or asymptotic running time. However, the order can have a considerable impact on practical performance due to thememory access patterns andcache use of the algorithm;^[1]which order is best also depends on whether the matrices are stored inrow-major order, column-major order, or a mix of both.

In particular, in the idealized case of afully associative cache consisting ofM bytes andb bytes per cache line (i.e.⁠M/b⁠ cache lines), the above algorithm is sub-optimal forA andB stored in row-major order. Whenn >⁠M/b⁠, every iteration of the inner loop (a simultaneous sweep through a row ofA and a column ofB) incurs a cache miss when accessing an element ofB. This means that the algorithm incursΘ(n³) cache misses in the worst case. As of 2010^[update], the speed of memories compared to that of processors is such that the cache misses, rather than the actual calculations, dominate the running time for sizable matrices.^[4]

The optimal variant of the iterative algorithm forA andB in row-major layout is atiled version, where the matrix is implicitly divided into square tiles of size√M by√M:^[4][5]

In the idealized cache model, this algorithm incurs onlyΘ(⁠n³/b√M⁠) cache misses; the divisorb√M amounts to several orders of magnitude on modern machines, so that the actual calculations dominate the running time, rather than the cache misses.^[4]

An alternative to the iterative algorithm is thedivide-and-conquer algorithm for matrix multiplication. This relies on theblock partitioning

which works for all square matrices whose dimensions are powers of two, i.e., the shapes are2ⁿ × 2ⁿ for somen. The matrix product is now

which consists of eight multiplications of pairs of submatrices, followed by an addition step. The divide-and-conquer algorithm computes the smaller multiplicationsrecursively, using thescalar multiplicationc₁₁ =a₁₁b₁₁ as its base case.

The complexity of this algorithm as a function ofn is given by the recurrence^[3]

${\displaystyle T(1)=\mathrm{\Theta }(1);}$

${\displaystyle T(n)=8T(n/2)+\mathrm{\Theta }(n^2),}$

accounting for the eight recursive calls on matrices of sizen/2 andΘ(n²) to sum the four pairs of resulting matrices element-wise. Application of themaster theorem for divide-and-conquer recurrences shows this recursion to have the solutionΘ(n³), the same as the iterative algorithm.^[3]

A variant of this algorithm that works for matrices of arbitrary shapes and is faster in practice^[4] splits matrices in two instead of four submatrices, as follows.^[6]Splitting a matrix now means dividing it into two parts of equal size, or as close to equal sizes as possible in the case of odd dimensions.

The cache miss rate of recursive matrix multiplication is the same as that of atiled iterative version, but unlike that algorithm, the recursive algorithm iscache-oblivious:^[6] there is no tuning parameter required to get optimal cache performance, and it behaves well in amultiprogramming environment where cache sizes are effectively dynamic due to other processes taking up cache space.^[4](The simple iterative algorithm is cache-oblivious as well, but much slower in practice if the matrix layout is not adapted to the algorithm.)

The number of cache misses incurred by this algorithm, on a machine withM lines of ideal cache, each of sizeb bytes, is bounded by^[6]: 13

${\displaystyle \mathrm{\Theta }\left(m+n+p+\frac{mn+np+mp}{b}+\frac{mnp}{b\sqrt{M}}\right)}$

Algorithms exist that provide better running times than the straightforward ones. The first to be discovered wasStrassen's algorithm, devised byVolker 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${\displaystyle O(n^{{\log }_{2}7})\approx O(n^{2.807})}$. Strassen's algorithm is more complex, and thenumerical stability is reduced compared to the naïve algorithm,^[7] but it is faster in cases wheren > 100 or so^[1] and appears in several libraries, such asBLAS.^[8] It is very useful for large matrices over exact domains such asfinite fields, where numerical stability is not an issue.

Since Strassen's algorithm is actually used in practical numerical software and computer algebra systems, improving on the constants hidden in thebig-O notation has its merits. A table that compares key aspects of the improved version based on recursive multiplication of 2×2-block matrices via 7 block matrix multiplications follows. As usual,${\displaystyle n}$ gives the dimensions of the matrix and${\displaystyle M}$ designates the memory size.

Progress for Strassen-like recursive 2x2-block matrix multiplication

Year	Reference	#matrix multiplications per step	#matrix additions per step	total arithmetic operations	total I/O-complexity
1969	Strassen[9]	7	18	${\displaystyle 7n^{{\log }_{2}7}-6n^2}$	${\displaystyle 6{\left(\frac{\sqrt{3}n}{\sqrt{M}}\right)}^{{\log }_{2}7}\cdot M-18n^2+3M}$
1971	Winograd[10]	7	15	${\displaystyle 6n^{{\log }_{2}7}-5n^2}$	${\displaystyle 5{\left(\frac{\sqrt{3}n}{\sqrt{M}}\right)}^{{\log }_{2}7}\cdot M-15n^2+3M}$
2017	Karstadt, Schwartz[11]	7	12	${\displaystyle 5n^{{\log }_{2}7}-4n^2+3n^2{\log }_{2}n}$	${\displaystyle 4{\left(\frac{\sqrt{3}n}{\sqrt{M}}\right)}^{{\log }_{2}7}\cdot M-12n^2+3n^2\cdot {\log }_2\left(\frac{\sqrt{2}n}{\sqrt{M}}\right)+5M}$
2023	Schwartz, Vaknin[12]	7	12	${\displaystyle 5n^{{\log }_{2}7}-4n^2+1.5n^2{\log }_{2}n}$	${\displaystyle 4{\left(\frac{\sqrt{3}n}{\sqrt{M}}\right)}^{{\log }_{2}7}\cdot M-12n^2+1.5n^2\cdot {\log }_2\left(\frac{\sqrt{2}n}{\sqrt{M}}\right)+5M}$

It is known that a Strassen-like algorithm with a 2×2-block matrix step requires at least 7 block matrix multiplications. In 1976 Probert^[13] showed that such an algorithm requires at least 15 additions (including subtractions); however, a hidden assumption was that the blocks and the 2×2-block matrix are represented in the same basis. Karstadt and Schwartz computed in different bases and traded 3 additions for less expensive basis transformations. They also proved that one cannot go below 12 additions per step using different bases. In subsequent work Beniamini et el.^[14] applied this base-change trick to more general decompositions than 2×2-block matrices and improved the leading constant for their run times.

It is an open question intheoretical computer science how well Strassen's algorithm can be improved in terms ofasymptotic complexity. Thematrix multiplication exponent, usually denoted${\displaystyle \omega }$, is the smallest real number for which any${\displaystyle n\times n}$ matrix over a field can be multiplied together using${\displaystyle n^{\omega +o(1)}}$ field operations. The current best bound on${\displaystyle \omega }$ is, by Alman, Duan,Williams, Xu, Xu, and Zhou.^[2] This algorithm, like all other recent algorithms in this line of research, is a generalization of the Coppersmith–Winograd algorithm, which was given byDon Coppersmith andShmuel Winograd in 1990.^[15] The conceptual idea of these algorithms is similar to Strassen's algorithm: a way is devised for multiplying twok ×k-matrices with fewer thank³ multiplications, and this technique is applied recursively. However, the constant coefficient hidden by thebig-O notation is so large that these algorithms are only worthwhile for matrices that are too large to handle on present-day computers.^[16][17] Victor Pan proposed so-called feasible sub-cubic matrix multiplication algorithms with an exponent slightly above 2.77, but in return with a much smaller hidden constant coefficient.^[18]

Freivalds' algorithm is a simpleMonte Carlo algorithm that, given matricesA,B andC, verifies inΘ(n²) time ifAB =C.

In 2022,DeepMind introduced AlphaTensor, aneural network that used a single-player game analogy to invent thousands of matrix multiplication algorithms, including some previously discovered by humans and some that were not.^[19] Operations were restricted to thenon-commutative ground field^{[clarification needed]} (normal arithmetic) andfinite field${\displaystyle \mathbb{Z}/2\mathbb{Z}}$ (mod 2 arithmetic). The best "practical" (explicit low-rank decomposition of a matrix multiplication tensor) algorithm found ran in O(n^2.778).^[20] Finding low-rank decompositions of such tensors (and beyond) is NP-hard; optimal multiplication even for 3×3 matricesremains unknown, even in commutative field.^[20] On 4×4 matrices, AlphaTensor unexpectedly discovered a solution with 47 multiplication steps, an improvement over the 49 required with Strassen’s algorithm of 1969, albeit restricted to mod 2 arithmetic. Similarly, AlphaTensor solved 5×5 matrices with 96 rather than Strassen's 98 steps. Based on the surprising discovery that such improvements exist, other researchers were quickly able to find a similar independent 4×4 algorithm, and separately tweaked Deepmind's 96-step 5×5 algorithm down to 95 steps in mod 2 arithmetic and to 97^[21] in normal arithmetic.^[22] Some algorithms were completely new: for example, (4, 5, 5) was improved to 76 steps from a baseline of 80 in both normal and mod 2 arithmetic.

Thedivide-and-conquer algorithm sketched earlier can beparallelized in two ways forshared-memory multiprocessors. These are based on the fact that the eight recursive matrix multiplications in

can be performed independently of each other, as can the four summations (although the algorithm needs to "join" the multiplications before doing the summations). Exploiting the full parallelism of the problem, one obtains an algorithm that can be expressed infork–join stylepseudocode:^[23]

Here,fork is a keyword that signal a computation may be run in parallel with the rest of the function call, whilejoin waits for all previously "forked" computations to complete.partition achieves its goal by pointer manipulation only.

This algorithm has acritical path length ofΘ(log²n) steps, meaning it takes that much time on an ideal machine with an infinite number of processors; therefore, it has a maximum possiblespeedup ofΘ(n³/log²n) on any real computer. The algorithm isn't practical due to the communication cost inherent in moving data to and from the temporary matrixT, but a more practical variant achievesΘ(n²) speedup, without using a temporary matrix.^[23]

On modern architectures with hierarchical memory, the cost of loading and storing input matrix elements tends to dominate the cost of arithmetic. On a single machine this is the amount of data transferred between RAM and cache, while on a distributed memory multi-node machine it is the amount transferred between nodes; in either case it is called thecommunication bandwidth. The naïve algorithm using three nested loops usesΩ(n³) communication bandwidth.

Cannon's algorithm, also known as the2D algorithm, is acommunication-avoiding algorithm that partitions each input matrix into a block matrix whose elements are submatrices of size√M/3 by√M/3, whereM is the size of fast memory.^[24] The naïve algorithm is then used over the block matrices, computing products of submatrices entirely in fast memory. This reduces communication bandwidth toO(n³/√M), which is asymptotically optimal (for algorithms performingΩ(n³) computation).^[25][26]

In a distributed setting withp processors arranged in a√p by√p 2D mesh, one submatrix of the result can be assigned to each processor, and the product can be computed with each processor transmittingO(n²/√p) words, which is asymptotically optimal assuming that each node stores the minimumO(n²/p) elements.^[26] This can be improved by the3D algorithm, which arranges the processors in a 3D cube mesh, assigning every product of two input submatrices to a single processor. The result submatrices are then generated by performing a reduction over each row.^[27] This algorithm transmitsO(n²/p^2/3) words per processor, which is asymptotically optimal.^[26] However, this requires replicating each input matrix elementp^1/3 times, and so requires a factor ofp^1/3 more memory than is needed to store the inputs. This algorithm can be combined with Strassen to further reduce runtime.^[27] "2.5D" algorithms provide a continuous tradeoff between memory usage and communication bandwidth.^[28] On modern distributed computing environments such asMapReduce, specialized multiplication algorithms have been developed.^[29]

There are a variety of algorithms for multiplication onmeshes. For multiplication of twon×n on a standard two-dimensional mesh using the 2DCannon's algorithm, one can complete the multiplication in 3n-2 steps although this is reduced to half this number for repeated computations.^[30] The standard array is inefficient because the data from the two matrices does not arrive simultaneously and it must be padded with zeroes.

The result is even faster on a two-layered cross-wired mesh, where only 2n-1 steps are needed.^[31] The performance improves further for repeated computations leading to 100% efficiency.^[32] The cross-wired mesh array may be seen as a special case of a non-planar (i.e. multilayered) processing structure.^[33]

In a 3D mesh withn³ processing elements, two matrices can be multiplied in${\displaystyle \mathcal{O}(\log n)}$ using the DNS algorithm.^[34]

• Computational complexity of mathematical operations
• Computational complexity of matrix multiplication
• CYK algorithm § Valiant's algorithm
• Matrix chain multiplication
• Method of Four Russians
• Multiplication algorithm
• Sparse matrix–vector multiplication

• Numerical linear algebra
• Matrix theory
• Matrix multiplication algorithms

• CS1: long volume value
• Articles with short description
• Short description matches Wikidata
• Articles containing potentially dated statements from September 2025
• All articles containing potentially dated statements
• Articles containing potentially dated statements from 2010
• Wikipedia articles needing clarification from May 2025
• Articles with example pseudocode