Offensive Security Researcher & Quantum Cryptography Analyst

Strassen’s Algorithm for Matrix Multiplication: Divide-and-Conquer Approach

Abdul Wahab Junaid
Strassen’s Algorithm for Matrix Multiplication (Divide-and-Conquer Approach)

Matrix multiplication is one of those operations I use everywhere — graphics, machine learning, physics simulations, cryptography — so its efficiency matters a lot to me. The naive way of multiplying two $n \times n$ matrices takes $O(n^3)$ time, and for a long time that was assumed to be close to optimal. Strassen’s algorithm surprised the field by showing I can do better using a clever divide-and-conquer trick that reduces the number of recursive multiplications needed, at the cost of a few extra additions. It was one of the first results to show that “obvious” algorithms aren’t always the best ones.

History and Background

Strassen’s algorithm was published by the German mathematician Volker Strassen in 1969, in a paper titled “Gaussian Elimination is not Optimal.” Before Strassen’s work, it was widely believed that $O(n^3)$ was essentially unavoidable for matrix multiplication. Strassen’s breakthrough — cutting the required number of multiplications from 8 down to 7 in each recursive step — showed that this belief was false and kicked off an entire subfield of research into fast matrix multiplication. Since then, further improvements have been made, including the Coppersmith–Winograd algorithm (1990) and its successors, though these later algorithms tend to have such large constant factors that they’re mostly of theoretical interest, while Strassen’s algorithm remains practical for real-world use.

Problem Statement

Given two $n \times n$ matrices $A$ and $B$, I want to compute their product $C = A \times B$ faster than the standard $O(n^3)$ algorithm, ideally by reducing the number of scalar multiplications required, since multiplication is typically more expensive than addition.

Core Concepts

How It Works

Given two $n \times n$ matrices $A$ and $B$ (assuming $n$ is a power of 2 for simplicity), I split each into four $n/2 \times n/2$ blocks:

$$ A = \begin{pmatrix} A_{11} & A_{12} \ A_{21} & A_{22} \end{pmatrix}, \quad B = \begin{pmatrix} B_{11} & B_{12} \ B_{21} & B_{22} \end{pmatrix} $$

  1. I compute 7 products (instead of the naive 8) using specific linear combinations of the submatrices.
  2. I combine these 7 products using additions and subtractions to form the four quadrants of the result matrix $C$.
  3. I apply this process recursively to each of the 7 products, until the submatrices are small enough (e.g., $1\times1$ or some small threshold) to multiply directly.
  4. I assemble the final result matrix $C$ from its four quadrants.

Working Principle

The core insight is that scalar (or submatrix) multiplication is more computationally expensive than addition, especially as matrix sizes grow, since multiplication recurses while addition doesn’t. By trading one multiplication for several extra additions, Strassen reduces the total number of recursive multiplication calls from 8 to 7 per level. Because this saving compounds recursively across $\log_2 n$ levels, the overall complexity improves from $n^3$ to $n^{\log_2 7} \approx n^{2.807}$, which is asymptotically better for large $n$.

Mathematical Foundation

Strassen’s algorithm defines the following seven products:

$$ \begin{aligned} M_1 &= (A_{11} + A_{22})(B_{11} + B_{22}) \ M_2 &= (A_{21} + A_{22}) B_{11} \ M_3 &= A_{11} (B_{12} – B_{22}) \ M_4 &= A_{22} (B_{21} – B_{11}) \ M_5 &= (A_{11} + A_{12}) B_{22} \ M_6 &= (A_{21} – A_{11})(B_{11} + B_{12}) \ M_7 &= (A_{12} – A_{22})(B_{21} + B_{22}) \end{aligned} $$

The four quadrants of the result $C$ are then:

$$ \begin{aligned} C_{11} &= M_1 + M_4 – M_5 + M_7 \ C_{12} &= M_3 + M_5 \ C_{21} &= M_2 + M_4 \ C_{22} &= M_1 – M_2 + M_3 + M_6 \end{aligned} $$

I can verify $C_{11}$ algebraically: the standard formula requires $C_{11} = A_{11}B_{11} + A_{12}B_{21}$. Expanding $M_1 + M_4 – M_5 + M_7$ in terms of the $A$ and $B$ blocks and simplifying confirms it equals $A_{11}B_{11} + A_{12}B_{21}$, though I won’t expand the full algebra here since it’s a lengthy but mechanical verification.

The recurrence for the running time is:

$$ T(n) = 7T\left(\frac{n}{2}\right) + \Theta(n^2) $$

Applying the Master Method with $a = 7$, $b = 2$, and $f(n) = \Theta(n^2)$: I compute $n^{\log_2 7} \approx n^{2.807}$. Since $f(n) = n^2$ is polynomially smaller than $n^{2.807}$, this falls into Case 1:

$$ T(n) = \Theta(n^{\log_2 7}) \approx \Theta(n^{2.807}) $$

Diagrams

flowchart TD
    A["A (n x n), B (n x n)"] --> B["Split into 2x2 blocks: A11, A12, A21, A22, B11, B12, B21, B22"]
    B --> C["Compute 7 products M1..M7 using combined sums/differences"]
    C --> D["Combine M1..M7 into result quadrants C11, C12, C21, C22"]
    D --> E["Assemble final result matrix C"]

Pseudocode

STRASSEN(A, B, n):
    if n == 1:
        return A[0][0] * B[0][0]

    // split A and B into quadrants
    A11, A12, A21, A22 = SPLIT(A)
    B11, B12, B21, B22 = SPLIT(B)

    M1 = STRASSEN(A11 + A22, B11 + B22, n/2)
    M2 = STRASSEN(A21 + A22, B11, n/2)
    M3 = STRASSEN(A11, B12 - B22, n/2)
    M4 = STRASSEN(A22, B21 - B11, n/2)
    M5 = STRASSEN(A11 + A12, B22, n/2)
    M6 = STRASSEN(A21 - A11, B11 + B12, n/2)
    M7 = STRASSEN(A12 - A22, B21 + B22, n/2)

    C11 = M1 + M4 - M5 + M7
    C12 = M3 + M5
    C21 = M2 + M4
    C22 = M1 - M2 + M3 + M6

    return COMBINE(C11, C12, C21, C22)

Step-by-Step Example

Let me multiply two $2\times 2$ matrices:

$$ A = \begin{pmatrix} 1 & 2 \ 3 & 4 \end{pmatrix}, \quad B = \begin{pmatrix} 5 & 6 \ 7 & 8 \end{pmatrix} $$

Here $A_{11}=1, A_{12}=2, A_{21}=3, A_{22}=4$, and $B_{11}=5, B_{12}=6, B_{21}=7, B_{22}=8$ (treating them as $1\times1$ blocks since the matrix itself is $2\times2$).

$$ \begin{aligned} M_1 &= (1+4)(5+8) = 5 \times 13 = 65 \ M_2 &= (3+4)\times 5 = 7 \times 5 = 35 \ M_3 &= 1 \times (6-8) = 1 \times (-2) = -2 \ M_4 &= 4 \times (7-5) = 4 \times 2 = 8 \ M_5 &= (1+2)\times 8 = 3 \times 8 = 24 \ M_6 &= (3-1)(5+6) = 2\times 11 = 22 \ M_7 &= (2-4)(7+8) = -2\times 15 = -30 \end{aligned} $$

Now the quadrants:

$$ \begin{aligned} C_{11} &= 65 + 8 – 24 + (-30) = 19 \ C_{12} &= -2 + 24 = 22 \ C_{21} &= 35 + 8 = 43 \ C_{22} &= 65 – 35 + (-2) + 22 = 50 \end{aligned} $$

So:

$$ C = \begin{pmatrix} 19 & 22 \ 43 & 50 \end{pmatrix} $$

I can double check this against standard matrix multiplication: $1\times5 + 2\times7 = 19$, $1\times6+2\times8=22$, $3\times5+4\times7=43$, $3\times6+4\times8=50$ — it matches exactly.

Time Complexity

Space Complexity

Strassen’s algorithm requires additional space for storing the 7 intermediate products and the various sums/differences of submatrices at each recursion level, giving a space complexity of $O(n^2)$ (same asymptotic order as the input/output matrices themselves), though with a larger constant factor than the naive algorithm due to the temporary matrices needed at each level of recursion.

Correctness Analysis

The correctness of Strassen’s algorithm follows directly from the algebraic identities defining $M_1$ through $M_7$ and the quadrant formulas — these can be verified by fully expanding both sides in terms of the original submatrix entries and confirming they match the standard matrix multiplication formulas for each quadrant. Since this holds at every level of recursion, and the base case (multiplying $1\times1$ “matrices,” i.e., scalars) is trivially correct, the algorithm is correct by induction on the recursion depth.

Advantages

Disadvantages

Applications

Implementation in C

#include <stdio.h>
#include <stdlib.h>

/* Allocates an n x n matrix of ints, dynamically. */
int** allocMatrix(int n) {
    int** mat = malloc(n * sizeof(int*));
    for (int i = 0; i < n; i++)
        mat[i] = calloc(n, sizeof(int));
    return mat;
}

void freeMatrix(int** mat, int n) {
    for (int i = 0; i < n; i++) free(mat[i]);
    free(mat);
}

/* Adds two n x n matrices: result = A + B */
int** addMatrix(int** A, int** B, int n) {
    int** result = allocMatrix(n);
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            result[i][j] = A[i][j] + B[i][j];
    return result;
}

/* Subtracts two n x n matrices: result = A - B */
int** subMatrix(int** A, int** B, int n) {
    int** result = allocMatrix(n);
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            result[i][j] = A[i][j] - B[i][j];
    return result;
}

/* Recursively multiplies two n x n matrices using Strassen's algorithm.
   Assumes n is a power of 2. Falls back to direct multiplication at n == 1. */
int** strassen(int** A, int** B, int n) {
    int** C = allocMatrix(n);

    if (n == 1) {
        C[0][0] = A[0][0] * B[0][0];
        return C;
    }

    int half = n / 2;
    int** A11 = allocMatrix(half); int** A12 = allocMatrix(half);
    int** A21 = allocMatrix(half); int** A22 = allocMatrix(half);
    int** B11 = allocMatrix(half); int** B12 = allocMatrix(half);
    int** B21 = allocMatrix(half); int** B22 = allocMatrix(half);

    /* split A and B into quadrants */
    for (int i = 0; i < half; i++) {
        for (int j = 0; j < half; j++) {
            A11[i][j] = A[i][j];
            A12[i][j] = A[i][j + half];
            A21[i][j] = A[i + half][j];
            A22[i][j] = A[i + half][j + half];

            B11[i][j] = B[i][j];
            B12[i][j] = B[i][j + half];
            B21[i][j] = B[i + half][j];
            B22[i][j] = B[i + half][j + half];
        }
    }

    /* seven recursive multiplications */
    int** t1 = addMatrix(A11, A22, half);
    int** t2 = addMatrix(B11, B22, half);
    int** M1 = strassen(t1, t2, half);

    int** t3 = addMatrix(A21, A22, half);
    int** M2 = strassen(t3, B11, half);

    int** t4 = subMatrix(B12, B22, half);
    int** M3 = strassen(A11, t4, half);

    int** t5 = subMatrix(B21, B11, half);
    int** M4 = strassen(A22, t5, half);

    int** t6 = addMatrix(A11, A12, half);
    int** M5 = strassen(t6, B22, half);

    int** t7 = subMatrix(A21, A11, half);
    int** t8 = addMatrix(B11, B12, half);
    int** M6 = strassen(t7, t8, half);

    int** t9 = subMatrix(A12, A22, half);
    int** t10 = addMatrix(B21, B22, half);
    int** M7 = strassen(t9, t10, half);

    /* combine into result quadrants */
    for (int i = 0; i < half; i++) {
        for (int j = 0; j < half; j++) {
            C[i][j] = M1[i][j] + M4[i][j] - M5[i][j] + M7[i][j];
            C[i][j + half] = M3[i][j] + M5[i][j];
            C[i + half][j] = M2[i][j] + M4[i][j];
            C[i + half][j + half] = M1[i][j] - M2[i][j] + M3[i][j] + M6[i][j];
        }
    }

    /* free all temporary matrices */
    int** temps[] = {A11,A12,A21,A22,B11,B12,B21,B22,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,M1,M2,M3,M4,M5,M6,M7};
    for (int k = 0; k < 25; k++) freeMatrix(temps[k], half);

    return C;
}

int main() {
    int n = 2;
    int** A = allocMatrix(n);
    int** B = allocMatrix(n);

    A[0][0]=1; A[0][1]=2; A[1][0]=3; A[1][1]=4;
    B[0][0]=5; B[0][1]=6; B[1][0]=7; B[1][1]=8;

    int** C = strassen(A, B, n);

    printf("Result matrix C:\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++)
            printf("%d ", C[i][j]);
        printf("\n");
    }

    freeMatrix(A, n);
    freeMatrix(B, n);
    freeMatrix(C, n);
    return 0;
}

Sample Input and Output

For the input matrices $A = \begin{pmatrix}1&2\3&4\end{pmatrix}$ and $B = \begin{pmatrix}5&6\7&8\end{pmatrix}$, the program outputs:

Result matrix C:
19 22 
43 50

This matches the hand-computed result from the step-by-step example above.

Optimization Techniques

Common Mistakes

Further Reading

Exit mobile version