Matrix multiplication cost

Why this matters

Matrices appear in games (moving and rotating shapes), in image filters, in machine learning, and in scientific simulations. A program that multiplies two large matrices can run in a fraction of a second or take minutes, depending on the method it uses and how big the matrices are.

This lesson does three things. It counts the work of the usual method and shows why the count grows so fast. It shows Strassen's method, a divide and conquer method that multiplies large matrices with fewer multiplications. It ends with the cases where the extra bookkeeping of Strassen's method costs more than it saves.

The lesson uses the splitting idea from the earlier divide and conquer lessons in this module. The Master theorem lesson gives the formula used at the end of the Strassen section.

What a matrix is

A matrix is a grid of numbers with rows and columns. A matrix with 2 rows and 3 columns is called a 2 by 3 matrix. Each number in the grid is an entry.

In C++, a matrix is a std::vector of rows, and each row is a std::vector of numbers. The entry in row i and column j is m[i][j]. Counting starts at 0, so m[0][0] is the top left entry.

A = 1 2        B = 5 6
    3 4            7 8

Here A and B are 2 by 2 matrices. In C++ they are written as {{1, 2}, {3, 4}} and {{5, 6}, {7, 8}}.

The rule for multiplying

To multiply two matrices A and B, each entry of the result C is built from one row of A and one column of B. Take the first number of the row and multiply it by the first number of the column. Do the same for the second pair, the third pair, and so on. Add all the products together.

For the entry in row 0 and column 0 of C, the row is 1 2 and the column is 5 7:

C[0][0] = 1 x 5 + 2 x 7 = 5 + 14 = 19

The other three entries follow the same pattern:

C[0][1] = 1 x 6 + 2 x 8 = 22     (row 0 of A, column 1 of B)
C[1][0] = 3 x 5 + 4 x 7 = 43     (row 1 of A, column 0 of B)
C[1][1] = 3 x 6 + 4 x 8 = 50     (row 1 of A, column 1 of B)

So the product is C = {{19, 22}, {43, 50}}.

For two matrices of size n by n, each entry of C uses one row of A with n numbers and one column of B with n numbers. That means n multiplications for each entry. The result has n by n entries.

Counting the multiplications

The rule above gives a simple count. There are n by n entries in C, and each entry needs n multiplications. The total is:

n x n x n = n^3 multiplications

The additions follow the same pattern. Each entry needs n - 1 additions to add its n products, so the total is close to n^3 as well.

We write this cost as O(n^3). The notation O, called big O, describes how the work grows when n grows. It keeps the part that grows fastest and drops the constant factors. For matrix multiplication, the important fact is that doubling n multiplies the work by 8, because 2 cubed is 8.

The next program runs the usual method on small matrices and counts each multiplication. It uses a counter that increases on every multiplication.

C++ · editable

#include <iostream>
#include <vector>

using Matrix = std::vector<std::vector<long long>>;

long long multiplications = 0;

Matrix naive_multiply(const Matrix& a, const Matrix& b) {
    const std::size_t n = a.size();
    Matrix c(n, std::vector<long long>(n, 0));
    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = 0; j < n; ++j) {
            for (std::size_t k = 0; k < n; ++k) {
                c[i][j] += a[i][k] * b[k][j];
                ++multiplications;
            }
        }
    }
    return c;
}

int main() {
    Matrix a = {{1, 2}, {3, 4}};
    Matrix b = {{5, 6}, {7, 8}};
    Matrix c = naive_multiply(a, b);
    std::cout << "C = [[" << c[0][0] << ", " << c[0][1] << "], [" << c[1][0] << ", " << c[1][1] << "]]\n";
    std::cout << "multiplications: " << multiplications << "\n\n";

std::cout << "n   multiplications\n";
    for (std::size_t n = 2; n <= 5; ++n) {
        Matrix x(n, std::vector<long long>(n, 1));
        multiplications = 0;
        naive_multiply(x, x);
        std::cout << n << "   " << multiplications << "\n";
    }
}

The output is:

C = [[19, 22], [43, 50]]
multiplications: 8
 
n   multiplications
2   8
3   27
4   64
5   125

The count for each size is n cubed: 2 cubed is 8, 3 cubed is 27, 4 cubed is 64, and 5 cubed is 125.

Why this is slow for large matrices

Take a matrix of size 1000 by 1000. The usual method needs 1000 cubed, which is one billion multiplications. Each multiplication also needs a memory read. A program that does one billion multiplications takes noticeable time even on a fast processor. A matrix of size 10000 by 10000 needs one trillion multiplications, which is 1000 times more work.

This is why people look for a method with fewer multiplications at large sizes.

Strassen's method

In 1969, Volker Strassen published a method that multiplies two 2 by 2 matrices with 7 multiplications instead of 8. The saving is small for 2 by 2 matrices. The idea matters for large matrices, because each big block can be handled by the same trick, and the saving repeats at each level of splitting.

Splitting into blocks

A block is a smaller square piece of a matrix. Split an n by n matrix into four blocks of size n/2 by n/2. Call them A11 (top left), A12 (top right), A21 (bottom left), and A22 (bottom right). Split B the same way.

The usual method computes each block of C from block products:

C11 = A11 x B11 + A12 x B21
C12 = A11 x B12 + A12 x B22
C21 = A21 x B11 + A22 x B21
C22 = A21 x B12 + A22 x B22

That is 8 block products. Each block product is a multiplication of two n/2 by n/2 matrices, so the total work is 8 times the work for half size. Since 8 times one eighth of n cubed is n cubed again, the split alone does not save anything.

Strassen's method uses 7 block products, built from sums and differences of blocks:

M1 = (A11 + A22) x (B11 + B22)
M2 = (A21 + A22) x B11
M3 = A11 x (B12 - B22)
M4 = A22 x (B21 - B11)
M5 = (A11 + A12) x B22
M6 = (A21 - A11) x (B11 + B12)
M7 = (A12 - A22) x (B21 + B22)
 
C11 = M1 + M4 - M5 + M7
C12 = M3 + M5
C21 = M2 + M4
C22 = M1 - M2 + M3 + M6

These formulas are fixed. You can check them by expanding the products and comparing terms. They are hard to find by hand, which is why the method is known by its inventor's name.

The formulas with numbers

Start with 1 by 1 blocks, which are plain numbers. Take A = {{1, 2}, {3, 4}} and B = {{5, 6}, {7, 8}}. The blocks are A11 = 1, A12 = 2, A21 = 3, A22 = 4, B11 = 5, B12 = 6, B21 = 7, and B22 = 8.

M1 = (1 + 4) x (5 + 8) = 5 x 13 = 65
M2 = (3 + 4) x 5       = 7 x 5  = 35
M3 = 1 x (6 - 8)       = 1 x -2 = -2
M4 = 4 x (7 - 5)       = 4 x 2  = 8
M5 = (1 + 2) x 8       = 3 x 8  = 24
M6 = (3 - 1) x (5 + 6) = 2 x 11 = 22
M7 = (2 - 4) x (7 + 8) = -2 x 15 = -30
 
C11 = 65 + 8 - 24 - 30 = 19
C12 = -2 + 24 = 22
C21 = 35 + 8 = 43
C22 = 65 - 35 - 2 + 22 = 50

The result is {{19, 22}, {43, 50}}, the same as the usual method. The count of multiplications dropped from 8 to 7. The count of additions and subtractions rose. The usual method needs 4 additions for this 2 by 2 case. Strassen's formulas need 18: 10 to form the sums of blocks before the products, and 8 to combine the products into C.

Why the saving repeats

When the blocks are large, each block product is itself a matrix multiplication of half size. Strassen's method can be applied to it again, and again, until the blocks are small. If the size is a power of two, written n = 2^k, each level of the recursion turns one product into 7 products of half the size. After k levels there are 7 to the power k products of size 1. Since n = 2^k, 7^k equals n raised to the power log2(7), which is about n^2.81.

So the work is about n^2.81 instead of n^3. For large n the gap grows: n^2.81 is much smaller than n^3 when n is large.

The extra additions cost about n^2 at each level. The recursion has about log2(n) levels, and the work is dominated by the multiplications at the bottom.

The program below uses Strassen's method on powers of two. It counts the multiplications in two ways. The first run uses pure recursion down to 1 by 1 blocks, so the counts match the 7 to the power k formula. The second run stops at blocks of size 64 and uses the usual method below that size, which is what a practical program does.

C++ · editable

#include <chrono>
#include <cstddef>
#include <iostream>
#include <random>
#include <vector>

using Matrix = std::vector<std::vector<long long>>;

long long scalar_multiplications = 0;

Matrix naive_multiply(const Matrix& a, const Matrix& b) {
    const std::size_t n = a.size();
    Matrix c(n, std::vector<long long>(n, 0));
    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = 0; j < n; ++j) {
            for (std::size_t k = 0; k < n; ++k) {
                c[i][j] += a[i][k] * b[k][j];
                ++scalar_multiplications;
            }
        }
    }
    return c;
}

Matrix add(const Matrix& a, const Matrix& b, long long sign = 1) {
    const std::size_t n = a.size();
    Matrix c(n, std::vector<long long>(n));
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j)
            c[i][j] = a[i][j] + sign * b[i][j];
    return c;
}

Matrix block(const Matrix& m, std::size_t row, std::size_t col, std::size_t size) {
    Matrix out(size, std::vector<long long>(size));
    for (std::size_t i = 0; i < size; ++i)
        for (std::size_t j = 0; j < size; ++j)
            out[i][j] = m[row + i][col + j];
    return out;
}

void place(Matrix& dest, const Matrix& part, std::size_t row, std::size_t col) {
    for (std::size_t i = 0; i < part.size(); ++i)
        for (std::size_t j = 0; j < part.size(); ++j)
            dest[row + i][col + j] = part[i][j];
}

// n must be a power of two. Blocks of size <= cutoff use the naive method.
Matrix strassen(const Matrix& a, const Matrix& b, std::size_t cutoff) {
    const std::size_t n = a.size();
    if (n <= cutoff) {
        return naive_multiply(a, b);
    }
    const std::size_t h = n / 2;
    const Matrix a11 = block(a, 0, 0, h), a12 = block(a, 0, h, h);
    const Matrix a21 = block(a, h, 0, h), a22 = block(a, h, h, h);
    const Matrix b11 = block(b, 0, 0, h), b12 = block(b, 0, h, h);
    const Matrix b21 = block(b, h, 0, h), b22 = block(b, h, h, h);

const Matrix m1 = strassen(add(a11, a22), add(b11, b22), cutoff);
    const Matrix m2 = strassen(add(a21, a22), b11, cutoff);
    const Matrix m3 = strassen(a11, add(b12, b22, -1), cutoff);
    const Matrix m4 = strassen(a22, add(b21, b11, -1), cutoff);
    const Matrix m5 = strassen(add(a11, a12), b22, cutoff);
    const Matrix m6 = strassen(add(a21, a11, -1), add(b11, b12), cutoff);
    const Matrix m7 = strassen(add(a12, a22, -1), add(b21, b22), cutoff);

Matrix c(n, std::vector<long long>(n));
    place(c, add(add(add(m1, m4), m5, -1), m7), 0, 0);
    place(c, add(m3, m5), 0, h);
    place(c, add(m2, m4), h, 0);
    place(c, add(add(add(m1, m3), m2, -1), m6), h, h);
    return c;
}

int main() {
    std::cout << "Strassen, pure recursion (cutoff 1)\n";
    std::cout << "n    naive mults   strassen mults\n";
    for (std::size_t n : {2, 4, 8, 16, 32, 64}) {
        Matrix x(n, std::vector<long long>(n, 1)), y = x;
        scalar_multiplications = 0;
        naive_multiply(x, y);
        long long naive_count = scalar_multiplications;
        scalar_multiplications = 0;
        strassen(x, y, 1);
        std::cout << n << "    " << naive_count << "    " << scalar_multiplications << "\n";
    }

std::cout << "\nTiming, cutoff 64 (one run each, -O2)\n";
    std::mt19937 rng(7);
    std::uniform_int_distribution<long long> dist(-9, 9);
    for (std::size_t n : {128, 256, 512}) {
        Matrix x(n, std::vector<long long>(n)), y = x;
        for (auto& row : x) for (auto& v : row) v = dist(rng);
        for (auto& row : y) for (auto& v : row) v = dist(rng);
        auto t0 = std::chrono::steady_clock::now();
        Matrix c1 = naive_multiply(x, y);
        auto t1 = std::chrono::steady_clock::now();
        Matrix c2 = strassen(x, y, 64);
        auto t2 = std::chrono::steady_clock::now();
        auto ms = [](auto d) { return std::chrono::duration<double, std::milli>(d).count(); };
        std::cout << "n=" << n << "  naive " << ms(t1 - t0) << " ms, strassen " << ms(t2 - t1)
                  << " ms, same result: " << (c1 == c2 ? "yes" : "no") << "\n";
    }
}

The output on the machine used for this lesson, built with g++ -std=c++20 -Wall -O2, was:

Strassen, pure recursion (cutoff 1)
n    naive mults   strassen mults
2    8    7
4    64    49
8    512    343
16    4096    2401
32    32768    16807
64    262144    117649
 
Timing, cutoff 64 (one run each, -O2)
n=128  naive 1.88624 ms, strassen 1.99572 ms, same result: yes
n=256  naive 12.1493 ms, strassen 13.2111 ms, same result: yes
n=512  naive 138.101 ms, strassen 68.9411 ms, same result: yes

The count column for Strassen is 7 to the power k. Each row doubles n and multiplies the Strassen count by 7, while the naive count is multiplied by 8.

The timing numbers vary from run to run. Another run, with the naive loops in a different order, measured the 512 case at 36 ms for the naive method and 42 ms for Strassen. The timings only show that the answer depends on the machine, the loop order and the size. They do not give a fixed crossover size, the size where one method starts to beat the other.

When the extra work is not worth it

Strassen's method saves multiplications and pays for it in other work. Four costs decide whether the trade is worth it.

Extra additions and copies. Each level forms 18 sums and differences for every 2 by 2 split. The code also copies blocks into new matrices. These steps take about n^2 work per level. At small sizes this work is large compared with the multiplications saved, so the usual method is faster.

Memory. Every level creates several temporary matrices. A program that runs out of memory fails before the time saving can matter.

Cache use. A processor keeps recently used numbers in a small fast memory called the cache. The usual method can be written so that it reuses numbers while they are still in the cache. Strassen's block copies spread the work over many temporary matrices, which can reduce that reuse. The result depends on the machine.

Rounding error. Strassen's method adds and subtracts intermediate sums before multiplying. For decimal numbers, these steps can lose more precision than the usual method. Programs that need exact results with whole numbers do not have this problem. Programs that use decimal numbers should compare the error before they choose the method.

Two more points matter in practice.

If the size is not a power of two, the usual approach is to pad the matrices with zeros up to the next power of two. Padding has a cost. A 100 by 100 matrix padded to 128 by 128 makes the usual method do 128 cubed divided by 100 cubed, which is about 2.1 times more multiplications. The padding must be paid for before the saving can help.

Production numerical libraries spend most of their effort on a tuned version of the usual method. That version splits the work into small blocks that fit in the cache. Strassen-type methods appear in some libraries for very large sizes. For a real program, the right step is to measure the code on the target machine with the target sizes, as the timing above does.

Common mistakes

Using Strassen on a size that is not a power of two. The split n / 2 assumes that n is even at each level. When n is 3, n / 2 is 1, so the four blocks are 1 by 1. The code then uses only the top left 2 by 2 part of the 3 by 3 matrix.

We called strassen from the program above on a 3 by 3 identity matrix, with cutoff 1. It printed:

1 2 0 
4 5 0 
0 0 0

The product of a matrix and the identity should be the same matrix, 1 2 3 / 4 5 6 / 7 8 9. The last row and column are lost, and no warning appears. A correct program pads to 4 by 4 first, or uses the usual method for sizes that are not powers of two.

Passing the wrong number of arguments. The function takes a cutoff as its third argument. Calling it with only two arguments fails to compile. g++ 15.2.0 printed this:

missing_arg.cpp: In function ‘int main()’:
missing_arg.cpp:6:24: error: too few arguments to function ‘Matrix strassen(const Matrix&, const Matrix&, std::size_t)’
    6 |     Matrix c = strassen(a, a);
      |                ~~~~~~~~^~~~~~
missing_arg.cpp:3:8: note: declared here
    3 | Matrix strassen(const Matrix& a, const Matrix& b, std::size_t cutoff);
      |        ^~~~~~~~

The message names the missing argument and the signature. Add the cutoff, or give it a default value.

Copying a sign wrong. The formula for M3 uses B12 - B22. If you write B12 + B22, the program compiles and runs, and the product is wrong. Compare each formula with the list above when you write the code, and test the result against the usual method on a small random matrix.

Counting only multiplications. Strassen's method saves one multiplication in eight at the 2 by 2 level, but it adds 18 additions where the usual method adds 4. A count that ignores the additions gives a false picture of the work.

Assuming Strassen is always faster. In the timing run above, the usual method was faster at 128, the two were close at 256, and Strassen was faster at 512. In the other run, Strassen was slower at 512. The benefit appears only when the size is large enough to pay for the extra work. Measure before you choose.

Practice

Exercise 1. How many multiplications does the usual method use to multiply two 5 by 5 matrices?

Show the answer

It uses 5 cubed, which is 125 multiplications. The program above printed 125 for n = 5.

Exercise 2. Use pure Strassen recursion (cutoff 1) on two 8 by 8 matrices. How many scalar multiplications does it use?

Show the answer

Each level of the recursion turns one product into 7 products of half the size. Eight is 2 cubed, so there are 7 cubed products of size 1, which is 343. The program printed 343.

Exercise 3. Use Strassen with a cutoff of 4 on two 16 by 16 matrices. The blocks of size 4 or less use the usual method. How many scalar multiplications are used?

Show the answer

The first split turns one 16 by 16 product into 7 products of size 8. The second split turns each of those into 7 products of size 4. That gives 49 blocks of size 4. Each block of size 4 uses 4 cubed, which is 64, multiplications. The total is 49 times 64, which is 3136. Calling strassen with cutoff 4 on 16 by 16 matrices counts 3136 multiplications.

Exercise 4. A program must multiply two 100 by 100 matrices with the usual method. Another program pads them to 128 by 128 first. How many times more multiplications does the padded version use?

Show the answer

The padded version uses 128 cubed divided by 100 cubed. That is 2097152 divided by 1000000, which is 2.097152. The padded version uses about 2.1 times more multiplications, so padding alone makes this case slower.

Exercise 5. Why did the 3 by 3 identity test give a wrong answer without an error?

Show the answer

The split uses n / 2, which is 1 for n = 3. The four blocks are then 1 by 1, and the code only reads the top left 2 by 2 part of the matrix. The last row and column are never used, and the code writes zeros there. The code has no check that n is a power of two, so it runs without a warning.

Summary

The usual method computes each of the n by n entries with n multiplications, so it does n cubed multiplications, written O(n^3).
Strassen's method computes a 2 by 2 product with 7 block products instead of 8, at the cost of 18 additions and subtractions.
Applied again and again to power-of-two sizes, Strassen's method does about n to the power 2.81 multiplications.
Extra additions, copies, memory and rounding error can make Strassen slower at small and medium sizes.
Non-power-of-two sizes need padding, and the padding adds work.
Measure on the target machine before choosing a method, because the crossover size depends on the machine.

Matrix multiplication cost

Why this matters

What a matrix is

The rule for multiplying

Counting the multiplications

Why this is slow for large matrices

Strassen's method

Splitting into blocks

The formulas with numbers

Why the saving repeats

When the extra work is not worth it

Common mistakes

Practice

Summary

Stuck? Ask the tutor