Matrix multiplication cost in Python

Why this matters

A matrix is a grid of numbers. Matrices show up wherever a whole table of numbers is transformed in one step. A 3D game multiplies a matrix for every object to move it, rotate it and place it on the screen. A machine learning model multiplies large weight matrices on every training step. A graph with n nodes can be stored as an n by n table, and multiplying that table by itself counts the routes of length two between nodes.

When the matrices are big, the cost of one multiplication decides whether a program finishes in a second or in an hour. This lesson measures that cost. It starts with the method you would write by hand, counts the work it does, and then shows Strassen's method, which is a divide and conquer method. It ends with the cases where the cleverer method is not worth using.

Plain Python runs every loop step through the interpreter, the program that reads and runs your Python code. That makes pure Python loops slow for big matrices. Libraries such as NumPy do the same multiplication in compiled code tuned for the machine. This lesson uses plain Python so that every step is visible. A real project with big matrices should use a library.

The lesson uses the result of the Master theorem, which is covered in its own lesson. Here, only the cost of matrix multiplication matters.

What a matrix is

A matrix has rows and columns. Its size is written as rows by columns. A 2 by 3 matrix has 2 rows and 3 columns. In this lesson every matrix is square with size n by n, so there are n rows and n columns.

Each number in the grid is called an entry. In Python, a matrix is a list of lists: each inner list is one row. The entry in row i and column j is m[i][j]. Python counts from 0, so the top left entry is m[0][0].

Multiplying two matrices A and B gives a new matrix C of the same size. The entry in row i and column j of C is a dot product: take row i of A and column j of B, multiply the pairs of numbers in the same position, and add the results. Here is a 2 by 2 example done by hand, with A equal to [[1, 2], [3, 4]] and B equal to [[5, 6], [7, 8]].

• Entry C[0][0] is 1 times 5, plus 2 times 7, which is 5 + 14 = 19.
• Entry C[0][1] is 1 times 6, plus 2 times 8, which is 6 + 16 = 22.
• Entry C[1][0] is 3 times 5, plus 4 times 7, which is 15 + 28 = 43.
• Entry C[1][1] is 3 times 6, plus 4 times 8, which is 18 + 32 = 50.

So C is [[19, 22], [43, 50]]. The program later in this lesson prints the same matrix.

Each entry needs two multiplications here. There are four entries, so the total is eight multiplications.

The naive method and its cost

The naive method is the rule above written as three nested loops. The outer loop picks a row of the result. The middle loop picks a value of k, which is the position inside the dot product. The inner loop goes across the columns of the result row. Each pass of the inner loop adds one product into one entry of C.

The output is:

[[19, 22], [43, 50]]

Changing the order of the three loops changes the order in which the work happens, not the amount of work.

The line c = [[0] * n for _ in range(n)] creates the result. The part [0] * n makes a row of n zeros. The list comprehension repeats that step n times, once for each row, so each row is a separate list. The common mistake section shows why this matters.

Counting the work

Count the multiplications. The innermost line runs once for each combination of i, k and j. There are n values for each of the three indexes, so the line runs n times n times n, which is n cubed. Computer scientists write this as O(n^3). O is Big-O notation: it describes how the number of steps grows as the input grows, and it drops constant factors and small terms.

The program below counts the multiplications for several sizes. It uses random entries with a fixed seed, so the numbers are the same on every run. The random.Random(1) object gives the same sequence each time it is created with the seed 1. The counts do not depend on the entry values, so they are the same on every machine.

The output is:

n=2 multiplications=8 n^3=8
n=4 multiplications=64 n^3=64
n=8 multiplications=512 n^3=512
n=16 multiplications=4096 n^3=4096
n=32 multiplications=32768 n^3=32768

Each time n doubles, the work grows by a factor of 8, because 2 cubed is 8. For a 1000 by 1000 matrix, the naive method does 1000 cubed, which is 1,000,000,000 multiplications. On the machine used for the timings later in this lesson, the naive method takes about 70 ms for a 128 by 128 product. Since the work grows with n cubed, a 1000 by 1000 product takes roughly 33 seconds with the same code. (1000 cubed divided by 128 cubed is about 476, and 476 times 70 ms is about 33 seconds.) Memory effects at that size would change the number.

Splitting a matrix into blocks

Divide and conquer works on matrices by splitting each one into four smaller squares, called blocks. If the matrix has size n by n, each block has size n/2 by n/2. Name the blocks by their position: A11 is the top left block, A12 is the top right, A21 is the bottom left and A22 is the bottom right. Do the same for B.

In Python, a block is taken with a slice. A slice like row[0:2] gives the values at indexes 0 and 1. The end index 2 is not included.

The product C can be written in terms of the blocks, using the same dot product rule as before but with blocks instead of single numbers.

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

Every product on the right is the multiplication of two matrices of size n/2 by n/2. That is the same problem as the original, only half as large. So a divide and conquer method solves eight smaller products and then adds the results.

Let T(n) be the work for an n by n matrix. The work is eight products of half size, plus the cost of the additions, which is about n squared. Written as a recurrence, T(n) = 8 T(n/2) + n squared. The Master theorem says that this recurrence gives n cubed again. Splitting into blocks changes the shape of the work, but not its size. To get a real gain, the method must do fewer block products. Strassen found a way to do that.

Strassen's trick: seven products instead of eight

Volker Strassen published his method in 1969. The idea is to compute seven products of blocks, not eight. Each of the seven products is built from sums and differences of blocks. The four blocks of C are then built from those seven products, using only additions and subtractions.

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

This looks like a trick with no meaning, so check it with numbers. Take the same 2 by 2 matrices as before, and treat each number as a block of size 1 by 1. Then A11 is 1, A12 is 2, A21 is 3, A22 is 4, B11 is 5, B12 is 6, B21 is 7 and B22 is 8. The seven products are:

• M1 = (1 + 4) times (5 + 8) = 5 times 13 = 65
• M2 = (3 + 4) times 5 = 7 times 5 = 35
• M3 = 1 times (6 - 8) = 1 times -2 = -2
• M4 = 4 times (7 - 5) = 4 times 2 = 8
• M5 = (1 + 2) times 8 = 3 times 8 = 24
• M6 = (3 - 1) times (5 + 6) = 2 times 11 = 22
• M7 = (2 - 4) times (7 + 8) = -2 times 15 = -30

Now the four entries of C:

• C11 = 65 + 8 - 24 - 30 = 19
• C12 = -2 + 24 = 22
• C21 = 35 + 8 = 43
• C22 = 65 - 35 - 2 + 22 = 50

These are the same four numbers as the naive method gave. The seven products used 7 multiplications, where the naive method used 8. For 1 by 1 blocks the saving is one multiplication. The saving is worth much more when the blocks are large, because each block product is itself split again with the same trick.

The cost of the new method

Each level of the recursion makes 7 smaller products instead of 8. Written as a recurrence, T(n) = 7 T(n/2) + n squared. The Master theorem gives the growth rate n raised to log base 2 of 7, which is about n raised to 2.807. This is written O(n^2.807). The exponent is smaller than 3, so for large n the method does less work.

In a test run of the same recursion with the cutoff (the block size where the recursion stops) set to 1, the count of single multiplications was 7, 49, 343 and 2401 for sizes 2, 4, 8 and 16. Each of those is a power of 7, as the recurrence predicts.

For a 1024 by 1024 matrix, the difference is large. Using only the multiplications and ignoring the additions, n cubed is about 1.07 billion, and 1024 raised to 2.807 is about 0.28 billion. The ratio is about 3.8. The exact numbers come from the recurrence, so treat them as a guide.

The recursion also has a depth limit to know about. Each level halves the size, so a 1024 by 1024 matrix goes 1024, 512, 256, 128, 64 before the cutoff, which is 4 levels. Python stops a recursion that goes more than 1000 calls deep. The limit is 1000 by default, which sys.getrecursionlimit() reports. This method never comes close to it.

The full program

The program below is a complete Strassen implementation with three extras. It stops the recursion at blocks of size 64 and uses the naive method below that size, because very small blocks cost more in overhead than they save. It pads any size up to a power of two with zeros, so a 100 by 100 matrix is computed as a 128 by 128 matrix and then cut back down. It also checks that the two matrices have the same size.

The output is:

2x2 naive:    [[19, 22], [43, 50]]
2x2 strassen: [[19, 22], [43, 50]]
n=3 multiply matches naive: True
n=5 multiply matches naive: True
n=100 multiply matches naive: True
n=128 multiply matches naive: True
error: expected two square matrices of the same size, got 2 and 3

Read the code from the top. naive is the base case. block takes one quarter of a matrix with slicing. Each slice makes a new list, so the recursive calls never change the caller's matrices. strassen makes the four input blocks, builds the sums and differences for the seven products, makes the seven recursive calls, and then builds the four result blocks. join puts the four result blocks back into one matrix, by joining each pair of rows with + on lists. multiply is the entry point that a caller should use. It checks the sizes, pads to a power of two and cuts the answer back to the original size.

Measuring it

The table below comes from one run of each size on WSL2 Linux with Python 3.14.4. The size 128 and 256 rows show the fastest of the runs listed. The size 512 row comes from a single run. Times depend on the machine and on what else is running, so treat the numbers as a sense of scale, not as fixed values.

At 128 and 256 the two methods take about the same time, and Strassen is a few percent slower. At 512 Strassen is about 1.1 times faster. That result comes from one run, so it may not hold on another run.

Each Strassen level builds many new lists for the sums and differences. At these sizes, that work takes about as long as the multiplications that Strassen saves.

When the extra work is not worth it

Strassen saves multiplications, but it adds work of another kind. Here are the costs, one by one.

Additions. Each level of the recursion does 18 block additions or subtractions. Ten of them build the operands of the seven products, and eight build the four result blocks. Each addition of two blocks of size n/2 costs about n squared divided by 4 entry operations. For small sizes, the additions cost about as much as the multiplications saved.

Memory. Each level creates new lists for the sums, the differences and the seven products. Every list holds references to int objects, so each block takes space for its own list structure as well as its numbers. The naive method creates one result matrix and no other lists.

Size. As the timing table shows, below some size the naive method is faster. The cutoff of 64 in the program is a choice. In a test run at size 256, a cutoff of 8 took 808.6 ms, while a cutoff of 64 took 573.6 ms. The best of two runs was used for both. A cutoff that is too small makes the program slower.

Accuracy with decimal numbers. The naive method adds products of numbers in a fixed order. Strassen subtracts sums of numbers that may be close to each other, which can lose digits when the numbers are decimals stored as floats. The result can then be slightly less accurate than the naive answer. With whole numbers stored as ints, there is no rounding, so the results match exactly.

Code that is hard to check. The seven products and the four combining formulas must all be right. A wrong sign in one formula gives a matrix that looks plausible but is wrong. The naive method is three loops that are easy to check by reading.

Here is what the choice means in a real program. If a program multiplies 4 by 4 matrices in a loop, Strassen would make it slower and use more memory, so the naive method is the right choice. If a program multiplies a few 2048 by 2048 matrices, plain Python will be slow whichever method it uses. A numerical library is the usual choice there, because it runs the work in compiled code tuned for the machine. Writing a fast matrix library is a separate skill, and this lesson only explains why the choice is not automatic.

Common mistakes

Assuming Strassen is always faster

The timing table shows that at size 256 the naive method is slightly faster on this machine. A method with fewer multiplications is not always faster in practice. Measure with the sizes you need before choosing.

Passing matrices that are not a power of two

The core strassen function assumes the size is a power of two, because it splits the size by 2 at each level. If the cutoff is set to 1 and a 3 by 3 list is passed in, the split gives 1 by 1 blocks and a 2 by 2 result. In one test run with the cutoff set to 1 and a 3 by 3 input, the function returned a 2 by 2 list and printed no error. The multiply function in the program avoids this by padding to a power of two first.

Mismatched sizes

If the two matrices have different sizes, the loops read past the end of a list. In one run that called strassen directly with a 4 by 4 matrix and a 2 by 2 matrix, the program stopped with this message:

IndexError: list index out of range

The message points to a position, not to the real cause. The multiply function checks the sizes first and raises a ValueError with a clear message, as shown in the program output above.

Building a matrix with * on a list of lists

[[0] * n] * n looks like it makes an n by n matrix, but it makes n references to one row. Changing one row changes all of them:

The output is:

bad: [[7, 0, 0], [7, 0, 0], [7, 0, 0]]
good: [[7, 0, 0], [0, 0, 0], [0, 0, 0]]

Use a list comprehension, as the code in this lesson does, so that each row is a separate list.

Python integers do not overflow

Java's int type holds numbers up to about 2 billion, so a product can overflow silently. A Python int has no fixed size. It grows as needed:

The output is:

product: 10000000000
big: 1267650600228229401496703205376

The limit to watch for is the float type. A float keeps about 53 bits of precision, so whole numbers above 2 to the power of 53 can lose their last digit when they are stored as floats:

The output is:

float: 9007199254740992.0
int sum: 18014398509481986
float sum: 1.8014398509481984e+16

Keep the entries of a matrix as ints when they are whole numbers. Do not convert them to floats.

A cutoff that is too small

Setting the cutoff to 1 makes the recursion go all the way down to single numbers. That creates a very large number of tiny lists. Keep the cutoff at a size where the naive method is still fast, such as 64, and measure before changing it.

Practice

1. How many multiplications does the naive method use to multiply two 6 by 6 matrices? Use the formula n cubed.

Show the answer

216. The size is n equal to 6, so n cubed is 6 times 6 times 6, which is 216.

217. Using the Strassen formulas with the cutoff set to 1, how many single-number multiplications does an 8 by 8 product use?

Show the answer

343. Each level of the recursion makes seven products, so the count is 7 cubed, which is 343. The test run in the cost section gave 7, 49, 343 and 2401 for sizes 2, 4, 8 and 16.

344. Multiply [[2, 1], [0, 3]] by [[1, 4], [2, 5]] with the naive method and write down the result.

Show the answer

[[4, 13], [6, 15]]. Entry C[0][0] is 2 times 1 plus 1 times 2, which is 4. Entry C[0][1] is 2 times 4 plus 1 times 5, which is 13. Entry C[1][0] is 0 times 1 plus 3 times 2, which is 6. Entry C[1][1] is 0 times 4 plus 3 times 5, which is 15.

4. The naive method takes about 70 ms for a 128 by 128 product on the machine in the timing table. Roughly how long would a 256 by 256 product take with the naive method? Ignore memory effects.

Show the answer

About 560 ms. Doubling n multiplies the work of the naive method by 2 cubed, which is 8, and 8 times 70 ms is 560 ms. The table shows 568.9 ms for size 256, which is close. Times on your machine will differ.

5. Explain why the multiply function in the program raises an error for a 2 by 2 matrix and a 3 by 3 matrix, instead of computing a result.

Show the answer

The two matrices must have the same size, because the block method splits both of them at the same place. The program checks the size of the second matrix against the first and stops with a message before any work starts. If it did not check, the naive loops would run without an error and give a wrong answer, because they would read only the top left 2 by 2 part of the 3 by 3 matrix.

Summary

• The naive method uses n cubed multiplications for n by n matrices, which is O(n^3).
• Splitting into four blocks and doing eight block products still gives n cubed, so the split alone does not help.
• Strassen's method uses seven block products and 18 additions or subtractions per level, which gives about n raised to 2.807.
• Strassen is not faster for small matrices, uses more memory and can lose accuracy with decimal numbers, so the choice depends on the size and the goal.
• In the timing run on this machine, Strassen was faster only at 512 and slower at 128 and 256.
• Check the sizes first, keep whole numbers as ints, and measure before choosing a method.