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.

SuperHive summer sale last day

 

The summer sale is closing. So this is your last chance at getting a serious discount on road signs in my shop, but there is loads of other good stuff on sale as well.

US Highway Shield Generator

 

After publishing the US Road Sign Generator I realized people might want a way to generate those highway shields that show the interstate numbers, so I created just that: All lightweight (both geometry and texture) and of course utilizing the common 'Highway Gothic' typeface.

Its available on Superhive (formerly BlenderMarket). Check it out if you want to go on a road trip!

SuperHive summer sale

Yup, the traditional summer sale is there again.  So if you are in the market for some road signs you may want to check my shop, but there is loads of other good stuff on sale as well.

Blender 5.2: Speeding up ImBuf access

In Blender 5.2, the ImBuf object now features a method that can be used as a context manager, yielding a Python memoryview object. This opens up a fast and convenient way to use an ImBuf to draw things using the blf module, and then copy the result to a Blender Image object's pixels attribute.

However, we cannot simply write:

with myimbuf.with_buffer() as mv:
    myimage.pixels[:] = mv[:]

This fails because Python's memoryview slice assignments have limitations when working with non-flat array configurations, which is what an Imbuf always is (width x height x color channels).  Instead, we have to cast the memoryview to a flat array. 

Furthermore, because either the source or the destination needs to be in a compatible byte format, we cannot cast a multidimensional float view directly to a flat float view in one go. We must resort to a rather inelegant double cast coupled with the foreach_set method:

myimage.pixels.foreach_set(mv.cast("b").cast("f"))

This first casts everything to a flat array of bytes, then to a flat array of floats, what is what foreach_set expects.

It works, and both cast() and foreach_set() are highly performant (cast() does not even touch the data, just the metadata), but it does feel a bit clunky. This isn't Blender's fault, though—it has been a long-standing behavior regarding multidimensional array manipulation discussed on the Python Core Development Issue Tracker.

Note that we now have a fast way to move data from an ImBuf to an Image, but ImBuf objects aren't all that common: we can create one and write text into it with the functions in the blf module, but as far as I can tell the is currently no way to copy a framebuffer to an Image or an ImBuf, nor is there a way to use the gpu module directly to draw into an ImBuf, so quickly creating an Image still requires a laborious route that involves copying the data from the framebuffer pixel-by-pixel to the Image.

Hopefully framebuffers will get buffer protocol support in the future, but that's probably not an easy task as there the actual data lives on the GPU and needs to be transferred. So let's keep an eye on the dev logs. 😀

US Road Sign Generator

After publishing a few sets of US road signs, I went back a created a geometry nodes based road sign generator that can be used to generate the iconic green signs you see along highways in the US.

Configurable to up to three lines, will all sort of options and three different materials. All lightweight (both geometry and texture) and of course utilizing the common 'Highway Gothic' typeface.

Its available on Superhive (formerly BlenderMarket). Check it out if you want to give your drivers some directions!

US Road Signs - Warning Signs Roadwork

After publishing a set of general US road signs, I thought why not publish one with signs commonly seen around road works?

The set is conveniently organized as an asset library and available on Superhive (formerly BlenderMarket). Check it out if you want to keep your rendered road crews and passing driver safe!