Strassen's Algorithm — Catastrophic Cancellation in Seismic

For decades, everyone assumed matrix multiplication required exactly n^3 multiplications for n×n matrices. In 1969, Volker Strassen proved them wrong by showing you only need 7 recursive multiplications instead of 8 for the 2×2 case. Applied recursively, this changes the exponent from log_2(8)=3 to log_2(7)≈2.807.

Strassen's algorithm matters less for its practical use (numerical stability issues and large constants make it impractical below n≈512) than for what it proved: the O(n^3) barrier is not a fundamental law. It opened the door to algorithms approaching O(n^2) — Coppersmith-Winograd (1987) achieved O(n^2.376), and subsequent work has continued pushing toward O(n^2). As of 2026, the theoretical lower bound is unknown. The entire sub-cubic matrix multiplication literature traces back to Strassen's 1969 insight.

Standard Matrix Multiplication

For C = A × B where A, B are n×n matrices: C[i][j] = Σ A[i][k] × B[k][j] for k in 0..n-1

This requires n³ multiplications and n²(n-1) additions → O(n³).

Strassen's 7-Multiplication Formula

For 2×2 matrices A and B split into quadrants: A = [[a,b],[c,d]], B = [[e,f],[g,h]]

Result: C = [[M1+M4-M5+M7, M3+M5],[M2+M4, M1-M2+M3+M6]]

def strassen(A, B):
    n = len(A)
    if n == 1:
        return [[A[0][0] * B[0][0]]]
    # Split into quadrants
    mid = n // 2
    a = [r[:mid] for r in A[:mid]]
    b = [r[mid:] for r in A[:mid]]
    c = [r[:mid] for r in A[mid:]]
    d = [r[mid:] for r in A[mid:]]
    e = [r[:mid] for r in B[:mid]]
    f = [r[mid:] for r in B[:mid]]
    g = [r[:mid] for r in B[mid:]]
    h = [r[mid:] for r in B[mid:]]

    def add(X, Y): return [[X[i][j]+Y[i][j] for j in range(len(X[0]))] for i in range(len(X))]
    def sub(X, Y): return [[X[i][j]-Y[i][j] for j in range(len(X[0]))] for i in range(len(X))]

    M1 = strassen(add(a,d), add(e,h))
    M2 = strassen(add(c,d), e)
    M3 = strassen(a, sub(f,h))
    M4 = strassen(d, sub(g,e))
    M5 = strassen(add(a,b), h)
    M6 = strassen(sub(c,a), add(e,f))
    M7 = strassen(sub(b,d), add(g,h))

    C11 = add(sub(add(M1,M4),M5), M7)
    C12 = add(M3, M5)
    C21 = add(M2, M4)
    C22 = add(sub(add(M1,M3),M2), M6)

    n2 = len(C11)
    C = [[0]*(2*n2) for _ in range(2*n2)]
    for i in range(n2):
        for j in range(n2):
            C[i][j]       = C11[i][j]
            C[i][j+n2]    = C12[i][j]
            C[i+n2][j]    = C21[i][j]
            C[i+n2][j+n2] = C22[i][j]
    return C

A = [[1,2],[3,4]]
B = [[5,6],[7,8]]
print(strassen(A,B))  # [[19,22],[43,50]]

Master Theorem Analysis

T(n) = 7T(n/2) + O(n²)

a=7, b=2, f(n)=n² log_b(a) = log₂(7) ≈ 2.807

Since f(n) = O(n^(2.807-ε)), we are in Master Theorem Case 1: T(n) = Θ(n^log₂7) ≈ Θ(n^2.807)

Saving just one multiplication per level of recursion changed the exponent from 3 to 2.807.

Crossover Point

Strassen is only faster than naive multiplication for large n (typically n > 128-512 depending on implementation). The large number of additions and the overhead of recursion make it slower for small matrices. Practical implementations switch to naive multiplication below a threshold.

Numerical Stability and When to Avoid Strassen

Strassen's recursion amplifies floating-point errors. Each of the 7 products involves addition of potentially large intermediate values, leading to catastrophic cancellation. For example, M6 = (c-a)(e+f): if c and a are close, the subtraction loses significant digits. At deeper recursion, these errors compound. In double precision, errors can reach 1e-8 relative for n=1000, while naive stays near machine epsilon.