[Numpy-discussion] speeding up operations on small vectors

[Christoph Groth]

Dear numpy experts,

I could not find a satisfying solution to the following problem, so I
thought I would ask:

In one part of a large program I have to deal a lot with small (2d or
3d) vectors and matrices, performing simple linear algebra operations
with them (dot products and matrix multiplications).  For several
reasons I want the code to be in python.

The most natural choice seemed to be to use numpy.  However, the
constant time cost when dealing with numpy arrays seems to be immense,
as demonstrated by this toy program:

****************************************************************
import numpy as np
from time import time

def points_tuple(radius):
    rr = radius**2
    def inside(point):
        return point[0]**2 + point[1]**2 < rr

    M = ((1, 0), (0, 1))
    for x in xrange(-radius, radius + 1):
        for y in xrange(-radius, radius + 1):
            point = (M[0][0] * x + M[0][1] * y,
                     M[1][0] * x + M[1][1] * y)
            if inside(point):
                yield point

def points_numpy(radius):
    rr = radius**2
    def inside(point):
        return np.dot(point, point) < rr

    M = np.identity(2, dtype=int)
    for x in xrange(-radius, radius + 1):
        for y in xrange(-radius, radius + 1):
            point = np.dot(M, (x, y))
            if inside(point):
                yield point

def main():
    r = 200
    for func in [points_tuple, points_numpy]:
        t = time()
        for point in func(r):
            pass
        print func.__name__, time() - t, 'seconds'

if __name__ == '__main__':
    main()
****************************************************************

On my trusty old machine the output (python 2.6, numpy 1.5.1) is:

points_tuple 0.36815404892 seconds
points_numpy 6.20338392258 seconds

I do not need C performance here, but the latter is definitely too slow.

In the real program it's not so easy (but possible) to use tuples
because the code is dimension-independent.  Before considering writing
an own "small vector" module, I'd like about other possible solutions.
Other people must have stumbled across this before!

Thanks,
Christoph