Efficiently multiplying a list of vectors with a matrix in Blender

There are many situations when writing Blender add-ons where you might want to multiply a vector with a matrix. A common example is a coordinate transformation, for example to convert object coordinates to world coordinates.

If you only need to convert one vector, it would be no problem to do this directly in Blender:

object_space_vector = Vector((1,0,0,1)) 

world_space_vector = object_space_vector @ context.active_object.matrix_world

(Note that the world matrix is 4x4, so it will not just scale and rotate, but also translate a vector, but then this vector need to be 4d as well! The 4th component, often referred to as 'w', should be 1 if you want to perform that translation, like for positions, and 0 if not, so for normals.)

The challenge

But what if we needed to perform such a transformation on a million vectors? Matrix multiplication itself  isn't exactly cheap, requiring 16 multiplications and 12 additions in the 4d case, and loops in Python are exceptionally slow.

Luckily Blender comes bundled with the numpy module and there are efficient ways to get properties, including vertex coordinates into numpy arrays (as illustrated in this article for example), so we only need a way to use numpy functionality to efficiently multiply a list of vectors with a matrix.

Choices

It is always good to have choices and numpy in fact offers several ways to achieve the same result.

The @ operator

The first is matrix multiplication. If vector is a list of vectors with shape = n,4 and matrix is a 4x4 matrix, we can perform the multiplication of each of the vectors in the list with this expression:

(matrix @ vector.T).T

Because numpy's @ operator expects the first dimensions of each of the operands to match, we have to transpose our list of 4d vectors first, and transpose again at the end to get the results in the right shape. This is a bit unintuitive, but at least those transpose operations are very cheap: they just create a view with a different shape without copying anything.

The alternative approach would be to invert the matrix and write vector @ matrix instead, but inverting a 4x4 transformation matrix is not simply np.linalg.inv(matrix) because we would need to extract the translation and scaling parts first, the invert, then recombine. That's too complicated to my taste, so I go with this approach.

The dot() function

The dot function essentially does the exact same thing, and needs the same transpositions:

np.dot(matrix, vector.T).T

This function is supposed to be more general, in the sense it can deal with tensors too, but that might impact performance. We have a look at that in a minute.

The Einstein summation notation 

The most versatile tool in numpy's drawer is the einsum() function, which lets us define the order of every implicit loop over the list of vectors and the rows and columns of the matrix: 

np.einsum("ni,ji->nj", vector, matrix)

This may look complicated but saves us the double transposition shenanigans, and reads something like: 

use index n to index each vector

    use index i for each element of a vector

        multiply that vector element with an element of the matrix:

            use index j for each row of the matrix

                use index i for each column of the row

        sum the results for the row (because i is used to index the matrix and the vector)

    store each of those row sums as the columns of the nth result vector             

This may take some time to wrap your head around, but it is very versatile and should be fast too.

Performance

Speaking of fast, lets have a look at some benchmarks:

n	power of two	pure python	dot	matmul	einsum
16384	14	0.038	0.000	0.000	0.000
32768	15	0.073	0.000	0.000	0.001
65536	16	0.143	0.000	0.000	0.001
131072	17	0.268	0.000	0.000	0.002
262144	18	0.543	0.000	0.000	0.005
524288	19	1.103	0.001	0.000	0.010
1048576	20	2.202	0.004	0.006	0.022
2097152	21		0.010	0.016	0.037
4194304	22		0.023	0.021	0.070
8388608	23		0.042	0.034	0.137
16777216	24		0.088	0.067	0.264
33554432	25		0.159	0.131	0.522
67108864	26		0.322	0.262	1.045
134217728	27		0.625	0.511	2.068
268435456	28		1.251	1.027
536870912	29		2.481	2.060

For comparison I have also included a pure Python reference implementation to show how slow that is. It is included in the code on Github.

We measure all the way from 16K 4d vectors to over 500 million of those. Tests that are slow were taken out of the race to save time.

The first thing we notice is that all methods scale with the number of vectors: if we double the number of vectors, the time to multiply them all by the same matrix approximately doubles.

The second thing we notice is that Python is really slow: It already takes more than 2 seconds just to multiply a million vectors. This is what we expected of course, so not really remarkable.

An finally, that the @ operator ('matmul') is fastest of them all. So the Einstein summation might be versatile and the dot() function more general but that does come with some overhead.

We can visualize that in the following graph:

(Note the logarithmic scales, so the y-axis top line is 1 second, next one down 0.1 second, etc)

Apart from some irregularities for shorter lists of vectors, the lines are pretty straight and clearly the @ operator is the winner. It might not look all that much in the graph because of the logarithmic scales, but it consistently is 20% faster than the dot() function, and 20% is nothing to be sniffed at in my book.

Caveats & conclusion

Of course your mileage may vary and the exact numbers aren't all that relevant. These timings were measured on a AMD Ryzen 7 7700X with 64 GB memory, and if you have a different CPU and/or amount of memory these timings will differ, but the trends are expected to be the same.

Also, all code uses 64-bit floating point operations and for Blender vertex coordinates that would be overkill. Using 32-bit floating point might save you quite some time, even on a 64-bit CPU, because the vectors would take only half the space, increasing cache efficiency quite a bit, but I didn't tabulate those timings.

In the end, the fastest and simplest to code is probably using the @ operator.