Matrix multiplication

Matrix multiplication combines an M by K matrix A with a K by N matrix B to produce an M by N matrix C, in which every entry of C is the dot product of one row of A with one column of B. In numerical computing the operation is normally called GEMM, for general matrix multiply, after the naming convention of the BLAS standard, and it is the single largest consumer of arithmetic in modern deep learning. Fully connected layers, recurrent layers, and convolutional layers all reduce to GEMM calls, and NVIDIA's own performance documentation treats GEMM as the fundamental building block of neural network computation ^[1].

Two very different research threads run through the subject. The first is hardware. Multiplying an M by K matrix by a K by N matrix costs 2 * M * N * K floating-point operations while touching only MK + KN + M*N values, so the ratio of arithmetic to memory traffic grows with the size of the matrices ^[1]. That property makes the operation an unusually good fit for wide parallel processors whose arithmetic units outrun their memory systems, which is one reason GPUs displaced CPUs for neural network training, why NVIDIA added dedicated tensor cores in 2017, and why Google built its tensor processing unit around a hardwired matrix unit.

The second thread is algorithmic, and it is much older. The schoolbook method needs n^3 scalar multiplications for two n by n matrices. Volker Strassen showed in 1969 that fewer suffice, and the search for the true minimum has continued ever since, most recently with help from machine learning systems such as AlphaTensor and AlphaEvolve. The two threads rarely meet: the algorithms that win on paper are almost never the ones running inside a training cluster.

The operation and its cost

For square matrices the textbook algorithm performs n^3 multiplications and roughly the same number of additions, giving O(n^3) arithmetic complexity. Because any correct algorithm must at least read both inputs, no method can do better than O(n^2), and the gap between those two exponents is the open problem discussed later in this article.

The Basic Linear Algebra Subprograms (BLAS) specification, maintained at Netlib, organizes dense linear algebra into three levels: Level 1 covers vector operations, Level 2 covers matrix-vector operations, and Level 3 covers matrix-matrix operations ^[2]. GEMM sits in Level 3, the level with enough arithmetic per byte of data moved to keep a modern processor busy, and it is the routine vendor libraries optimize hardest. Numerical stacks such as NumPy, PyTorch, and JAX route their dense linear algebra to a BLAS-compatible GEMM implementation.

The standard GEMM interface computes alpha * op(A) * op(B) + beta * C, where op() optionally transposes or conjugate-transposes an input. NVIDIA's cuBLAS implements exactly this form and supports FP16, BF16, FP32, FP64, FP8 in both E4M3 and E5M2 encodings, FP4 in the E2M1 encoding, integer and complex types, plus compute modes including TF32 and a BF16x9 emulation path ^[3].

Why neural networks are mostly matrix multiplication

A dense layer applies a weight matrix to a batch of activations, which is a matrix multiply by definition. Recurrent layers such as RNNs, LSTMs, and GRUs, and convolutional layers, are also built on GEMM ^[1]. The transformer architecture, which dominates current AI workloads, is built almost entirely from matrix products: the query, key, and value projections, the attention scores themselves, the output projection, and the two linear layers of the position-wise feed-forward block ^[4].

The share is easy to underestimate. Ivanov and colleagues profiled a BERT encoder layer at batch size 8 and sequence length 512 and classified every operator into tensor contractions (matrix multiplies), statistical normalizations, and element-wise operations. Tensor contractions accounted for 99.80 percent of the floating-point operations but only 61.0 percent of the measured runtime ^[5].

Scroll sideways for more →

That mismatch is the practical reason kernel fusion matters. Almost all the arithmetic lives in GEMM, but almost 40 percent of the time goes to softmax, layer normalization, and pointwise work that is limited by memory bandwidth rather than math. FlashAttention and similar fused kernels exist to close that gap.

The dependence is not absolute. Zhu and colleagues showed in 2024 that matrix multiplication can be removed from the attention and feed-forward blocks of a language model entirely, using ternary weights and element-wise operations instead, with performance competitive with standard transformers at scales up to 2.7 billion parameters, up to 61 percent less memory during training, and more than a tenfold reduction in memory during GPU inference ^[6]. That line of work remains a research direction rather than production practice.

Hardware built around the operation

Graphics processors were not designed for neural networks, and their early advantage came from raw parallel throughput rather than any matrix-specific circuitry. That changed with the Volta architecture in 2017. Each V100 streaming multiprocessor carries 8 tensor cores, 640 across the chip, and each one computes D = A * B + C on 4 by 4 matrices with FP16 inputs, a full-precision product, and FP32 accumulation. A single tensor core performs 64 mixed-precision fused multiply-adds per clock, giving 1,024 floating-point operations per clock per SM and up to 125 tensor TFLOPS for the whole GPU ^[7]. CUDA exposes them through the warp matrix multiply accumulate (WMMA) API, which issues 16 by 16 by 16 tiles across a warp ^[7].

Google took the more radical route. The first TPU, deployed in Google data centers from 2015 and described publicly in 2017, was built around a matrix multiply unit containing 65,536 8-bit multiply-accumulators in a systolic array, rated at 92 TeraOps per second with 28 MiB of software-managed on-chip memory. The paper reported roughly 15x to 30x the inference speed of contemporary CPUs and GPUs at 30x to 80x better performance per watt ^[8]. The design persists: current TPU documentation describes each TensorCore as one or more matrix multiply units (MXUs) plus a vector unit and a scalar unit, with each MXU built from a 256 by 256 systolic array of multiply-accumulators on TPU v6e and TPU7x, and 128 by 128 on earlier versions. All multiplies take bfloat16 inputs while accumulation runs in FP32 ^[9].

Other AI accelerators, including designs from Cerebras and Groq as well as FPGA based systems, face the same core problem: keeping large multiply-accumulate arrays supplied with data.

Arithmetic intensity and the memory wall

NVIDIA's performance guides frame kernel behavior in terms of arithmetic intensity, defined for GEMM as (M * N * K) / (MK + NK + M*N), compared against the processor's ops-to-byte ratio ^[1]. An operation whose intensity exceeds that ratio is math limited; below it, memory bandwidth decides ^[10]. Large square GEMMs sit comfortably in the math-limited regime, which is exactly why they can use tensor cores well. Skinny GEMMs, such as the single-token matrix-vector products of autoregressive decoding, do not, which is why KV cache sizing and batching dominate inference performance discussions.

The balance has moved steadily against memory. The V100 paired about 900 GB/s of memory bandwidth with 125 FP16 tensor TFLOPS; the A100 pairs up to 2,039 GB/s from 80 GB of HBM2 with 156 TF32 TFLOPS and 312 dense FP16 TFLOPS across 108 SMs ^[10]. Peak math grew faster than memory bandwidth, so tiling, data layout, and reuse became more important rather than less. GEMM kernels in cuBLAS use thread block tiles from 256x128 down to 64x64, with the larger tiles achieving better reuse and higher efficiency, and performance improves when M, N, and K are aligned to multiples of 16 bytes, or 128 bytes on A100 ^[1].

Writing these kernels by hand is impractical, so most work happens through libraries and compilers. CUTLASS provides C++ templates for GEMM across Volta through Blackwell along with the CuTe layout abstraction ^[11], cuDNN supplies the equivalent for neural network layers, and Triton offers a Python-level path to custom kernels.

Numerical precision

Because matrix multiplication is a sum of products, it tolerates low-precision inputs better than most numerical work, provided the accumulation is done in higher precision. That observation drives essentially all of the format engineering in AI hardware.

Mixed precision training was formalized by Micikevicius and colleagues in 2017: store weights, activations, and gradients in IEEE half precision, keep an FP32 master copy of the weights that accumulates optimizer updates, and scale the loss so that small gradient values survive the narrower FP16 range. The technique preserved accuracy on convolutional, recurrent, and generative models with more than 100 million parameters while roughly halving memory use ^[12].

Scroll sideways for more →

The FP8 interchange-format proposal defined the two encodings and reported that FP8 training matched the quality of 16-bit training on convolutional networks, recurrent networks, and transformers up to 175 billion parameters ^[14]. NVFP4 pushes further by shrinking the scaling block to 16 elements with an E4M3 scale factor, against the 32-element blocks and power-of-two E8M0 scales of the community MXFP4 format; NVIDIA reported pretraining a 12-billion-parameter hybrid Mamba-Transformer on 10 trillion tokens in NVFP4 with validation loss closely matching an FP8 baseline ^[16]. The general pattern is that inputs get narrower every generation while accumulators stay wide, and that quantization error is controlled by finer-grained scaling rather than by more mantissa bits.

Vendor throughput figures track the format ladder directly. NVIDIA rates the H100 SXM at 67 teraFLOPS of ordinary FP32, against 989 teraFLOPS of TF32, 1,979 teraFLOPS of BF16 or FP16, and 3,958 teraFLOPS of FP8 tensor-core throughput, with the tensor-core numbers quoted with sparsity and dense rates at half those values ^[17].

Fast algorithms and the exponent omega

Strassen's 1969 paper "Gaussian elimination is not optimal" in Numerische Mathematik showed that two 2 by 2 matrices can be multiplied with 7 scalar multiplications instead of 8, and that applying the scheme recursively gives an algorithm with exponent log2 7, about 2.8074 ^[18][19]. The result launched a search for the matrix multiplication exponent omega, defined as the smallest real number such that two n by n matrices can be multiplied in O(n^(omega+e)) field operations for every e greater than zero.

Scroll sideways for more →

Years for the entries after 1990 are those of the first preprint, and the pre-2020 bounds follow the standard tabulation of the problem's history ^[19]. Every bound since 1986 has come from refinements of the laser method, applied in recent work to higher powers of the Coppersmith-Winograd tensor ^[20][21]. Duan, Wu, and Zhou introduced asymmetric hashing to address combination loss in that analysis and reached 2.371866 ^[21]; Vassilevska Williams, Xu, Xu, and Zhou improved the variant to 2.371552 ^[22]; and the 2024 paper by Alman, Duan, Vassilevska Williams, Xu, Xu, and Zhou broke the symmetry requirement between two of the three dimensions to reach 2.371339, which remains the record ^[23].

None of these record-setting algorithms is used in practice. The constant factors hidden inside the asymptotic notation are so large that the methods only pay off for matrices too big for present-day computers to handle, which is why they are called galactic algorithms ^[19]. Strassen's method is the genuine exception, becoming faster than the schoolbook algorithm somewhere around n greater than 100 and appearing in several numerical libraries ^[19]. Deep learning kernels, however, are built the conventional way: cuBLAS and CUTLASS organize GEMM as tiled multiply-accumulate work mapped onto tensor-core instructions rather than as Strassen-style recursion ^[3][11].

Machine-discovered algorithms

In October 2022, Google DeepMind published AlphaTensor in Nature, framing the discovery of a matrix multiplication algorithm as a single-player game in which an agent based on AlphaZero decomposes the matrix multiplication tensor into rank-one terms ^[24]. For 4 by 4 matrices over the field with two elements it found a decomposition of rank 47, beating the 49 multiplications of Strassen's algorithm applied twice and marking the first improvement on Strassen in a finite field in 50 years ^[24][25]. For multiplying a 4 by 5 matrix by a 5 by 5 matrix it found a rank-76 decomposition against a previous best of 80 and a schoolbook count of 100 ^[24][25]. DeepMind also reported that specializing AlphaTensor's search for particular hardware produced algorithms that multiplied large matrices 10 to 20 percent faster than commonly used implementations on an NVIDIA V100 GPU and a Google TPU v2 ^[25].

The results did not stand unchallenged for long. Three days after the Nature paper appeared, Manuel Kauers and Jakob Moosbauer posted an algorithm multiplying 5 by 5 matrices over the field with two elements in 95 multiplications, one fewer than the 96 AlphaTensor had announced ^[26]. They followed in December 2022 with a flip graph method, a random walk over a graph of equivalent multiplication schemes, which improved counts for the (4,4,5) and (5,5,5) formats over both characteristic two and general ground fields ^[27]. The episode is a reasonable illustration of how the field works: the reinforcement learning result was real, and conventional search caught up almost immediately.

DeepMind returned to the problem in May 2025 with AlphaEvolve, an evolutionary coding agent that uses Gemini models to propose and refine program mutations under automated evaluation. It improved 14 different matrix multiplication targets, including a procedure for multiplying two 4 by 4 complex-valued matrices with 48 scalar multiplications, which the authors described as the first improvement over Strassen's algorithm in that setting after 56 years ^[28][29]. The paper notes that the discovered algorithms use complex-valued multiplications, which can be applied to exact multiplication of complex or real-valued matrices.

AlphaEvolve's applied results are arguably more consequential than its algebraic ones. DeepMind reported that a heuristic it found for the Borg cluster scheduler recovers on average 0.7 percent of the company's worldwide compute, that it sped up a matrix multiplication kernel used in Gemini training by 23 percent for a 1 percent reduction in overall training time, and that it improved a FlashAttention kernel implementation by up to 32.5 percent ^[28].

Recent developments

Hardware in 2024 through 2026 has continued narrowing the multiplier while widening the accumulator. The Blackwell architecture packs 208 billion transistors on a custom TSMC 4NP process and adds a second-generation Transformer Engine using fine-grain micro-tensor scaling to enable 4-bit floating point AI ^[30]. NVIDIA rates the 72-GPU GB200 NVL72 system at 1,440 petaFLOPS of NVFP4 tensor-core throughput with sparsity and 720 dense, against 360 petaFLOPS at FP16 or BF16 with sparsity, or 180 dense ^[31]. The successor Rubin generation carries what NVIDIA calls enhanced fifth-generation tensor cores tuned for NVFP4 and FP8 arithmetic, quoted at up to 50 petaFLOPS of NVFP4 inference, with the Vera Rubin NVL72 system quoted at 3,600 petaFLOPS ^[32].

Google announced its seventh-generation TPU, Ironwood, in April 2025, quoting peak FP8 compute of 4,614 teraFLOPS per chip, 192 GB of high bandwidth memory per chip at 7.37 TB/s (six times the capacity and 4.5 times the bandwidth of Trillium), and configurations up to 9,216 chips totaling 42.5 FP8 exaFLOPS ^[33].

On the software side, DeepSeek released DeepGEMM as an open tensor core kernel library covering FP8, FP4, and BF16 GEMMs for SM90 and SM100 GPUs, reporting up to 1,550 TFLOPS on an H800 ^[34]. NVIDIA's cuBLAS has meanwhile added emulation compute modes, listed alongside TF32 and including a BF16x9 path and fixed-point emulation, which run a GEMM at one precision using tensor-core instructions built for another ^[3].

See also

• Tensor core
• Systolic array
• AlphaTensor
• AlphaEvolve
• TPU
• Mixed precision training

References

1. ^NVIDIA, "Matrix Multiplication Background User's Guide", NVIDIA Deep Learning Performance Documentation. docs.nvidia.com/...rformance-matrix-multiplication
2. ^Netlib, "BLAS (Basic Linear Algebra Subprograms)". netlib.org/blas
3. ^NVIDIA, "cuBLAS Documentation", CUDA Toolkit Documentation. docs.nvidia.com/...cublas
4. ^Ashish Vaswani et al., "Attention Is All You Need", arXiv:1706.03762, 12 June 2017. arxiv.org/...1706.03762
5. ^Andrei Ivanov, Nikoli Dryden, Tal Ben-Nun, Shigang Li, Torsten Hoefler, "Data Movement Is All You Need: A Case Study on Optimizing Transformers", arXiv:2007.00072, 30 June 2020. arxiv.org/...2007.00072
6. ^Rui-Jie Zhu et al., "Scalable MatMul-free Language Modeling", arXiv:2406.02528, 4 June 2024. arxiv.org/...2406.02528
7. ^NVIDIA, "Programming Tensor Cores in CUDA 9", NVIDIA Technical Blog. developer.nvidia.com/...amming-tensor-cores-cuda-9
8. ^Norman P. Jouppi et al., "In-Datacenter Performance Analysis of a Tensor Processing Unit", arXiv:1704.04760, 16 April 2017. arxiv.org/...1704.04760
9. ^Google Cloud, "TPU architecture", Cloud TPU documentation. docs.cloud.google.com/...system-architecture-tpu-vm
10. ^NVIDIA, "GPU Performance Background User's Guide", NVIDIA Deep Learning Performance Documentation. docs.nvidia.com/...dl-performance-gpu-background
11. ^NVIDIA, "CUTLASS" repository README, GitHub. github.com/...cutlass
12. ^Paulius Micikevicius et al., "Mixed Precision Training", arXiv:1710.03740, 10 October 2017 (ICLR 2018). arxiv.org/...1710.03740
13. ^NVIDIA, "Accelerating AI Training with NVIDIA TF32 Tensor Cores", NVIDIA Technical Blog. developer.nvidia.com/...ing-with-tf32-tensor-cores
14. ^Paulius Micikevicius et al., "FP8 Formats for Deep Learning", arXiv:2209.05433, 12 September 2022. arxiv.org/...2209.05433
15. ^NVIDIA, "Floating-Point 8: An Introduction to Efficient, Lower-Precision AI Training", NVIDIA Technical Blog. developer.nvidia.com/...ower-precision-ai-training
16. ^NVIDIA, "NVFP4 Trains with Precision of 16-Bit and Speed and Efficiency of 4-Bit", NVIDIA Technical Blog. developer.nvidia.com/...ed-and-efficiency-of-4-bit
17. ^NVIDIA, "NVIDIA H100 Tensor Core GPU" product page and specification table. nvidia.com/...h100
18. ^Volker Strassen, "Gaussian elimination is not optimal", Numerische Mathematik, volume 13, pages 354-356, August 1969. link.springer.com/...BF02165411
19. ^Wikipedia, "Computational complexity of matrix multiplication". en.wikipedia.org/...exity_of_matrix_multiplication
20. ^Josh Alman, Virginia Vassilevska Williams, "A Refined Laser Method and Faster Matrix Multiplication", arXiv:2010.05846, 12 October 2020. arxiv.org/...2010.05846
21. ^Ran Duan, Hongxun Wu, Renfei Zhou, "Faster Matrix Multiplication via Asymmetric Hashing", arXiv:2210.10173, 18 October 2022. arxiv.org/...2210.10173
22. ^Virginia Vassilevska Williams, Yinzhan Xu, Zixuan Xu, Renfei Zhou, "New Bounds for Matrix Multiplication: from Alpha to Omega", arXiv:2307.07970, 16 July 2023. arxiv.org/...2307.07970
23. ^Josh Alman, Ran Duan, Virginia Vassilevska Williams, Yinzhan Xu, Zixuan Xu, Renfei Zhou, "More Asymmetry Yields Faster Matrix Multiplication", arXiv:2404.16349, 25 April 2024. arxiv.org/...2404.16349
24. ^Alhussein Fawzi et al., "Discovering faster matrix multiplication algorithms with reinforcement learning", Nature, volume 610, issue 7930, pages 47-53, 5 October 2022. pmc.ncbi.nlm.nih.gov/...PMC9534758
25. ^Google DeepMind, "Discovering novel algorithms with AlphaTensor", 5 October 2022. deepmind.google/...vel-algorithms-with-alphatensor
26. ^Manuel Kauers, Jakob Moosbauer, "The FBHHRBNRSSSHK-Algorithm for Multiplication in Z_2^{5x5} is still not the end of the story", arXiv:2210.04045, 8 October 2022. arxiv.org/...2210.04045
27. ^Manuel Kauers, Jakob Moosbauer, "Flip Graphs for Matrix Multiplication", arXiv:2212.01175, 2 December 2022. arxiv.org/...2212.01175
28. ^Google DeepMind, "AlphaEvolve: A Gemini-powered coding agent for designing advanced algorithms", 14 May 2025. deepmind.google/...r-designing-advanced-algorithms
29. ^Alexander Novikov et al., "AlphaEvolve: A coding agent for scientific and algorithmic discovery", arXiv:2506.13131, 16 June 2025. arxiv.org/...2506.13131
30. ^NVIDIA, "NVIDIA Blackwell Architecture" technology page. nvidia.com/...blackwell-architecture
31. ^NVIDIA, "NVIDIA GB200 NVL72" product page and specification table. nvidia.com/...gb200-nvl72
32. ^NVIDIA, "NVIDIA Tensor Cores" product page. nvidia.com/...tensor-cores
33. ^Google, "Ironwood: The first Google TPU for the age of inference", The Keyword, 9 April 2025. blog.google/...ironwood-tpu-age-of-inference
34. ^DeepSeek, "DeepGEMM" repository README, GitHub. github.com/...DeepGEMM