Strassen Algorithm

Strassen Algorithm is a recursive divide-and-conquer algorithm in computation and algorithms and algebraic complexity theory that multiplies two n × n matrices in sub-cubic asymptotic time O(n^{log₂ 7}) ≈ O(n^2.807355). Published by German mathematician Volker Strassen in a landmark four-page 1969 paper titled "Gaussian Elimination is Not Optimal," the algorithm overturned the long-standing assumption in computational complexity theory that multiplying two n × n matrices inherently requires n³ scalar multiplications (or 2n³ − n² total arithmetic operations) as defined by the row-by-column inner-product formula.

Analogous to how the Karatsuba algorithm reduces two-limb integer multiplication from 4 scalar products to 3, Strassen discovered a bilinear rank-7 decomposition of the 2 × 2 matrix multiplication tensor that replaces 8 block matrix multiplications with 7 block multiplications and 18 block additions/subtractions (later reduced to 15 additions by Shmuel Winograd in 1971). Because matrix multiplication underpins block LU decomposition, Gaussian elimination, matrix inversion, determinant evaluation, and transitive closure on graphs, Strassen's algorithm reduced the asymptotic complexity of all fundamental dense linear algebra operations to O(n^ω) where ω ≤ log₂ 7 ≈ 2.8074, initiating the modern search for the asymptotic matrix multiplication exponent ω.

Naive 2 × 2 Block Multiplication Versus Rank-7 Bilinear Decomposition

Let A, B ∈ R^n×n be square matrices over an arbitrary ring R (not necessarily commutative, allowing recursive application to submatrices) with dimension n = 2^k. Partition A, B, and their product C = AB into four (n/2) × (n/2) quadrant blocks:

Standard block multiplication computes the four quadrants of C via C_ij = A_i1B_1j + A_i2B_2j, which invokes 8 recursive multiplications of (n/2) × (n/2) matrices and 4 block additions, yielding the recurrence T(n) = 8T(n/2) + Θ(n²) = Θ(n^{log₂ 8}) = Θ(n³).

Strassen observed that by linearly combining blocks of A and B before multiplying, all four quadrants of C can be reconstructed from just seven intermediate (n/2) × (n/2) matrix products M₁, M₂, ..., M₇:

1. M₁ = (A₁₁ + A₂₂)(B₁₁ + B₂₂)
2. M₂ = (A₂₁ + A₂₂)B₁₁
3. M₃ = A₁₁(B₁₂ − B₂₂)
4. M₄ = A₂₂(B₂₁ − B₁₁)
5. M₅ = (A₁₁ + A₁₂)B₂₂
6. M₆ = (A₂₁ − A₁₁)(B₁₁ + B₁₂)
7. M₇ = (A₁₂ − A₂₂)(B₂₁ + B₂₂)

Notice that A_ij always appears on the left and B_kl always appears on the right in every product M₁..M₇, so the identities hold without assuming commutativity of matrix multiplication (A_ijB_kl ≠ B_klA_ij).

Reconstruction of the Output Quadrants and Winograd's 15-Addition Variant

Expanding the seven products M₁ through M₇ and cancelling cross terms yields the four output quadrants of C = AB:

• C₁₁ = M₁ + M₄ − M₅ + M₇ = A₁₁B₁₁ + A₁₂B₂₁
• C₁₂ = M₃ + M₅ = A₁₁B₁₂ + A₁₂B₂₂
• C₂₁ = M₂ + M₄ = A₂₁B₁₁ + A₂₂B₂₁
• C₂₂ = M₁ − M₂ + M₃ + M₆ = A₂₁B₁₂ + A₂₂B₂₂

Counting block additions and subtractions: forming the operands of M₁..M₇ requires 10 block additions/subtractions, and assembling C₁₁..C₂₂ requires 8 block additions/subtractions, totaling 18 block additions of size (n/2) × (n/2). In 1971, Shmuel Winograd proved that 7 multiplications is the mathematical minimum for 2 × 2 matrices (the tensor rank of ⟨2, 2, 2⟩ is 7) and introduced the Strassen–Winograd variant, which shares intermediate sums to reduce the count from 18 to 15 block additions (the theoretical lower bound for rank-7 schedules established by Robert Probert in 1976).

Master Theorem Solution and Exact Operation Counts

Applying the Master Theorem to the Strassen–Winograd recurrence for total arithmetic operations T(n) on n = 2^k:

Solving this recurrence in closed form yields:

For original 18-addition Strassen, the leading constant is 7n^{log₂ 7} − 6n². In practical implementations, recursion is not continued down to 1 × 1 scalars; instead, recursion terminates at a crossover threshold n₀ ∈ [512, 2048], below which submatrices are multiplied using cache-optimized Level-3 BLAS (dgemm). Stopping after r = 1, 2, or 3 recursive levels reduces total floating-point operations by (7/8)^r (12.5%, 23.4%, or 33.0% fewer flops) while preserving hardware vectorization. For non-power-of-two dimensions n, dynamic peeling or even-dimension padding avoids wasting work on zero rows.

Finite-Precision Stability and Subsequent Fast Matrix Multiplication Bounds

Because Strassen's formulas add and subtract submatrices from different rows and columns (such as A₁₁ + A₂₂), cancellation can amplify rounding errors when matrix entries vary widely in magnitude. Nicholas Higham (1990) proved that while classical matrix multiplication satisfies a componentwise backward error bound |Ĉ − AB| ≤ nε_mach|A||B| in IEEE 754 arithmetic, Strassen's algorithm satisfies only a normwise bound:

Thus, Strassen's algorithm is weakly stable: it preserves accuracy for well-scaled dense matrices or 1–3 levels of recursion, and is exact over finite fields GF(p) or arbitrary-precision integers (accelerated at the scalar level by the fast Fourier transform). Following Strassen (ω < 2.8074), Victor Pan (1978, ω < 2.796), Don Coppersmith and Shmuel Winograd (1990, ω < 2.375477), Virginia Vassilevska Williams (2012, ω < 2.37287), and Josh Alman and Vassilevska Williams (2021–2024, ω < 2.371339) pushed the theoretical exponent closer to the conjectured limit ω = 2, though Strassen–Winograd remains the primary sub-cubic algorithm with small enough constants for practical deployment.

How does Strassen's algorithm reduce the complexity of matrix inversion and Gaussian elimination?

In the same 1969 paper, Strassen showed that an invertible n × n matrix A partitioned into 2 × 2 blocks can be inverted via its Schur complement S = A₂₂ − A₂₁A₁₁⁻¹A₁₂ using 2 recursive (n/2) × (n/2) inversions and 6 (n/2) × (n/2) matrix multiplications. This recurrence I(n) = 2I(n/2) + O(n^ω) solves to I(n) = O(n^ω), reducing LU decomposition, Gaussian elimination, and determinant computation to O(n^2.8074).

Why don't standard BLAS libraries always use Strassen's algorithm by default?

Modern CPUs and GPUs have high arithmetic throughput relative to memory bandwidth. Strassen's algorithm requires temporary workspace arrays to hold intermediate block sums (increasing memory traffic) and replaces componentwise stability with normwise stability. Consequently, BLAS libraries use classical cache-blocked dgemm by default and invoke 1–3 levels of Strassen–Winograd only when n ≥ 1024–2048 and workspace memory is available.

What is the minimum known number of multiplications for 3 × 3 matrices?

Whereas the minimal bilinear rank for 2 × 2 matrices is known to be 7 (achieved by Strassen), the minimal rank for 3 × 3 matrices (naive cost 27 multiplications) remains an open problem in algebraic complexity: Julian Laderman constructed a 23-multiplication scheme in 1976 (log₃ 23 ≈ 2.854, slightly slower than Strassen's log₂ 7 ≈ 2.807), and the true minimum rank of 3 × 3 matrix multiplication is known to lie between 19 and 23.

1. Strassen, V. (1969). "Gaussian elimination is not optimal." Numerische Mathematik, 13(4), 354–356. https://doi.org/10.1007/BF02165411
2. Winograd, S. (1971). "On multiplication of 2 × 2 matrices." Linear Algebra and its Applications, 4(4), 381–388. https://doi.org/10.1016/0024-3795(71)90009-790009-7)
3. Higham, N. J. (1990). "Exploiting fast matrix multiplication within the level 3 BLAS." ACM Transactions on Mathematical Software, 16(4), 352–368. https://doi.org/10.1145/98267.98290
4. Coppersmith, D., & Winograd, S. (1990). "Matrix multiplication via arithmetic progressions." Journal of Symbolic Computation, 9(3), 251–280. https://doi.org/10.1016/S0747-7171(08)80013-280013-2)
5. Ballard, G., Demmel, J., Holtz, O., & Schwartz, O. (2012). "Graph expansion and communication costs of fast matrix multiplication." Journal of the ACM, 59(6), 32:1–32:23. https://doi.org/10.1145/2395116.2395121