Fast, small, and repetitive matrix multiplication in Python
Question: I'm looking for a way to very quickly multiply together many 4x4 matrices
using Python/Cython/Numpy, can anyone give any suggestions?

To show my current attempt, I have an algorithm which needs to compute

    
    
    A_1 * A_2 * A_3 * ... * A_N 
    

where every

    
    
    A_i != A_j
    

An example of it in Python:

    
    
    means = array([0.0, 0.0, 34.28, 0.0, 0.0, 3.4])
    stds = array([ 4.839339, 4.839339, 4.092728, 0.141421, 0.141421, 0.141421])
    
    def fn():
        steps = means+stds*numpy.random.normal(size=(60,6))
        A = identity(4)
        for step in steps:
            A = dot(A, transform_step_to_4by4(step))
    %timeit fn()
    
    1000 loops, best of 3: 570 us per loop
    

Implementing this algorithm in Cython/Numpy is approximately 100x slower than
the equivalent code using Eigen/C++ with all optimizations. I really don't
want to use C++, though.


Answer: If you are having to make a Python function call to produce each of the
matrices you want to multiply, then you are basically screwed performance-
wise. But if you can vectorize the `transform_step_to_4by4` function, and have
it return an array with shape `(n, 4, 4)` then you can save some time using
`matrix_multiply`:

    
    
    import numpy as np
    from numpy.core.umath_tests import matrix_multiply
    
    matrices = np.random.rand(64, 4, 4) - 0.5
    
    def mat_loop_reduce(m):
        ret = m[0]
        for x in m[1:]:
            ret = np.dot(ret, x)
        return ret
    
    def mat_reduce(m):
        while len(m) % 2 == 0:
            m = matrix_multiply(m[::2], m[1::2])
        return mat_loop_reduce(m)
    
    In [2]: %timeit mat_reduce(matrices)
    1000 loops, best of 3: 287 us per loop
    
    In [3]: %timeit mat_loop_reduce(matrices)
    1000 loops, best of 3: 721 us per loop
    
    In [4]: np.allclose(mat_loop_reduce(matrices), mat_reduce(matrices))
    Out[4]: True
    

You now have log(n) Python calls instead of n, good for a 2.5x speed-up, which
will get close to 10x for n = 1024. Apparently `matrix_multiply` is a ufunc,
and as such has a `.reduce` method, which would allow your code to run no
loops in Python. I have not been able to make it run though, keep getting an
arcane error:

    
    
    In [7]: matrix_multiply.reduce(matrices)
    ------------------------------------------------------------
    Traceback (most recent call last):
      File "<ipython console>", line 1, in <module>
    RuntimeError: Reduction not defined on ufunc with signature
    

