Wait, What?
The school formula for matrix multiplication is mathematically correct—and still leaves most of the real algorithmic problem untouched.
If A is m×k and B is k×n, each output entry is a dot product. That definition gives a correct O(mkn) algorithm, but professional performance depends on loop order, memory layout, blocking, vectorisation, parallelism, numerical precision, library kernels and sometimes asymptotically faster algorithms such as Strassen’s method.
This article owns the algorithmic job of forming matrix products efficiently. The existing Numerical Linear Algebra Algorithms article owns factorisations, conditioning and linear-system methods. The Cache-Efficient Algorithms article owns the general locality model. Here those ideas are applied to one of computing’s most important kernels: matrix multiplication.
Quick Answer
Learn matrix multiplication through the route shape compatibility → dot products → triple-loop algorithm → correctness invariant → loop permutations → memory layout → cache locality → blocking/tiling → vectorisation → Level-3 BLAS GEMM → arithmetic intensity → multicore parallelism → GPU execution → Strassen’s seven products → recurrence and exponent → crossover size → numerical error → rectangular and batched GEMM → benchmarking and validation. A beginner should be able to compute and code the classical method. A professional should be able to explain why mathematically identical loop nests perform differently, know when to call an optimised GEMM rather than hand-write loops, and distinguish asymptotic arithmetic complexity from real data-movement cost.
1. Start With Shape Before Arithmetic
A matrix product AB is defined when the number of columns of A equals the number of rows of B. If A is m×k and B is k×n, the result C is m×n. Shape errors should be caught before any loop begins.
2. Each Output Cell Is a Dot Product
The entry C[i,j] equals the sum over p of A[i,p]×B[p,j]. This definition immediately exposes the three logical indices: output row i, output column j and reduction index p.
3. The Classical Algorithm Has Three Nested Loops
A direct implementation iterates over output rows, output columns and the shared dimension. For square n×n matrices, this uses Θ(n³) scalar multiply-add work. Before optimisation, learners should be able to state exactly what each loop index means.
4. The Partial-Sum Invariant Proves Correctness
After processing the first t positions of the shared dimension, the accumulator for C[i,j] should equal the dot product restricted to p<t. This invariant makes the algorithm’s correctness almost mechanical and connects directly to the existing Algorithm Correctness Proofs article.
5. Six Loop Orders Can Compute the Same Product
The indices i, j and p can be nested in different orders. Mathematically, each correct ordering can produce the same result. Operationally, they touch memory differently. That makes loop order an early lesson in the gap between an abstract algorithm and its machine execution.
6. Memory Layout Decides Which Direction Is Cheap
Row-major storage places consecutive elements of a row next to each other; column-major storage does the analogous thing for columns. Accessing data in storage order usually improves spatial locality and cache-line use.
7. Naïve Triple Loops Can Waste Memory Bandwidth
If a loop repeatedly walks down a column in row-major memory, each useful value may arrive with many nearby values that are not immediately used. The arithmetic count is unchanged, but cache misses and data movement rise.
8. Blocking Reuses Small Tiles Before Eviction
Blocked multiplication partitions matrices into submatrices sized so useful working data can remain in faster memory. Instead of streaming entire rows and columns repeatedly, the algorithm performs many operations on a tile while its data is still close to the processor.
9. Blocking Preserves the Algebra
Partitioned matrices obey the same multiplication rule at block level. The algorithm is not approximating; it is regrouping operations to improve locality. This is a powerful general optimisation pattern: change execution order while preserving the mathematical sum.
10. Vectorisation Exploits Several Arithmetic Lanes at Once
Modern CPUs can perform operations on multiple numbers with one vector instruction. High-performance kernels organise data and loop structure so these lanes stay busy, often using small register-resident micro-kernels inside larger cache tiles.
11. GEMM Is the Professional Primitive
The Basic Linear Algebra Subprograms define standard matrix operations. Level-3 BLAS contains matrix–matrix routines, including GEMM for general matrix multiplication. Netlib describes BLAS as portable building blocks used throughout high-quality linear algebra software. See Netlib BLAS and the GEMM reference.
12. GEMM Computes More Than Plain AB
The standard form is C ← αAB + βC, with options to transpose A or B. That general form lets many applications fuse scaling and accumulation into the same optimised kernel instead of creating extra temporary matrices.
13. Arithmetic Intensity Explains Why GEMM Can Run Fast
Large dense matrix multiplication performs many arithmetic operations per value loaded from memory when blocking is effective. This high arithmetic intensity makes GEMM well suited to modern CPUs and GPUs, which have enormous compute throughput relative to main-memory bandwidth.
14. GPUs Change the Scale, Not the Mathematical Contract
GPU implementations distribute tiles across many parallel threads and rely on specialised memory hierarchies and matrix hardware. NVIDIA’s current performance guide frames GEMM as a fundamental building block of deep-learning operations and explains the interaction between matrix dimensions and hardware utilisation. See NVIDIA Matrix Multiplication Background User’s Guide.
15. Parallelism Needs Enough Work Per Task
Splitting tiny matrices across many threads can lose to scheduling and synchronisation overhead. The existing Parallel Algorithms article owns work/span reasoning; matrix multiplication provides a concrete case where task granularity and data placement determine whether more processors help.
16. Strassen’s Insight Is Seven Products Instead of Eight
For 2×2 block matrices, classical block multiplication needs eight recursive block products. Strassen discovered formulas that use seven block products plus additional additions and subtractions. Recursing gives O(n^log₂7), approximately O(n^2.81), arithmetic operations for square matrices.
17. Fewer Multiplications Do Not Automatically Mean Faster
Strassen introduces extra additions, temporaries, irregular memory traffic and recursion overhead. Practical implementations therefore switch back to highly tuned classical kernels below a crossover size. The best crossover depends on hardware, dimensions, precision and library engineering.
18. Asymptotic Complexity and Practical Performance Are Different Questions
Fast matrix-multiplication research continues to improve theoretical exponents, but algorithms with spectacular asymptotic bounds can have enormous constants or impractical structure. A 2025 Berkeley Simons Institute lecture notes that the current best theoretical exponent is below 2.37134 while also emphasising communication cost and practical reality. The professional lesson is to keep theory and workload evidence in the same conversation.
19. Numerical Error Also Depends on the Algorithm
Floating-point addition and multiplication are not exact real arithmetic. Different parenthesisations change rounding. Fast algorithms such as Strassen may have different error properties from conventional multiplication. Research in SIAM’s Journal on Matrix Analysis and Applications examines how scaling and algorithm design can improve the numerical stability of fast multiplication. See Improving the Numerical Stability of Fast Matrix Multiplication.
20. Integer and Floating-Point Products Need Different Checks
For small integers within a safe range, exact comparison against a reference implementation is straightforward. For floating-point matrices, validation should use tolerances grounded in scale and precision rather than demanding bit-for-bit equality between different operation orders.
21. Rectangular Matrices Change the Workload
Real GEMMs are often tall-skinny, short-wide or highly rectangular. A kernel tuned for large square matrices may underperform. Performance analysis should use the actual m, n and k distributions rather than assuming n×n.
22. Small Matrices Create a Different Regime
For very small dimensions, function-call overhead, packing and kernel setup can dominate. Intel’s oneMKL has documented JIT-generated kernels aimed specifically at improving small GEMM workloads, a reminder that “use the fastest asymptotic method” is not a professional performance rule.
23. Batched GEMM Amortises Repeated Small Problems
Many applications need thousands of small independent matrix products. Batched APIs group them so launch and scheduling overhead can be reduced and hardware can process more work concurrently.
24. Benchmark With Warm-Up, Repetition and Correctness Checks
- Warm up JITs, caches and device kernels where relevant.
- Separate allocation and data transfer from the kernel if the question is kernel speed.
- Measure several matrix shapes, not one.
- Record data type and threading configuration.
- Check the numerical result every time performance code changes.
- Report throughput together with matrix dimensions and hardware.
25. Common Learning Failure States
- Multiplying entries position-by-position instead of row-by-column.
- Ignoring shape compatibility.
- Memorising O(n³) without knowing the three indices.
- Assuming all loop orders perform the same on hardware.
- Calling blocking an approximation.
- Writing hand loops when a tuned BLAS should be used.
- Assuming Strassen is faster for every size.
- Ignoring temporary memory and data movement.
- Comparing floating-point results with exact equality.
- Benchmarking only square matrices or only one machine.
26. A Beginner-to-Professional Learning Ladder
- Level 1: verify matrix shapes and compute one output entry.
- Level 2: implement the classical triple loop.
- Level 3: state and use the partial-sum invariant.
- Level 4: compare loop orders under row-major and column-major storage.
- Level 5: implement simple tiling and measure cache effects.
- Level 6: call GEMM correctly with transposition and scaling options.
- Level 7: explain arithmetic intensity, vectorisation and parallel granularity.
- Level 8: derive Strassen’s recurrence and crossover trade-off.
- Level 9: validate floating-point products and benchmark multiple shapes.
- Level 10: choose kernels, precision, batching and hardware strategy for a production workload.
27. Teach the Indices Visually Before Optimising
Give each learner a small A, B and blank C. Highlight one row of A, one column of B and the accumulating C cell. Ask them to predict which memory locations will be touched next. Only after the index roles are stable should loop order and blocking be introduced.
28. Use Faded Worked Examples for Loop Transformations
Start with a fully annotated i-j-p loop showing memory accesses. Next remove cache annotations. Then provide a partially blocked loop with tile bounds missing. Finally ask the learner to reconstruct and benchmark the blocked version. Programming-education research on faded worked examples with metacognitive scaffolding found strong benefits for novice problem solving. See The Effects of Worked-Out Example and Metacognitive Scaffolding on Problem-Solving Programming.
29. Retrieve Both Formula and Performance Model
Delayed practice should include two different questions: reconstruct the C[i,j] formula and explain why a blocked version can outperform a naïve version despite the same Θ(n³) arithmetic count. That separation prevents complexity notation from becoming a substitute for machine reasoning.
30. Immediate, Delayed and Transfer Checks
- Immediate: compute a 2×3 times 3×2 product.
- Code: implement and test a classical loop nest.
- Invariant: explain the accumulator after t reduction steps.
- Performance: predict which loop order improves locality for a stated storage layout.
- Delayed: derive the purpose of tiling from memory hierarchy constraints.
- Transfer: decide whether to use a hand kernel, BLAS GEMM, batching or GPU execution for four workloads.
- Professional: benchmark dimensions, precision, threading, transfers and error together.
AI Assistance Boundary
AI can generate loop variants, benchmark harnesses and diagrams of cache tiles. The learner should still be able to verify the product, explain each index, check numerical error, interpret performance counters and decide whether an optimisation preserved the mathematical contract.
Professional Direction
Advanced study includes cache-oblivious multiplication, register micro-kernels, packing, roofline analysis, communication lower bounds, distributed GEMM, mixed precision, tensor cores, Strassen–Winograd variants, fast rectangular matrix multiplication, sparse matrix products and automatic kernel generation.
Algorithm-learning rule: when two matrix-multiplication programs have the same arithmetic formula, do not assume they have the same cost. Count what moves through memory, what fits in cache, what vector units can consume, what synchronisation is required, and what numerical error the chosen execution order introduces.
eduKate Sengkang 