[SciPy-User] fast small matrix multiplication with cython?

[Skipper Seabold]

On Tue, Dec 7, 2010 at 9:54 AM, Charles R Harris
<charlesr.harris at gmail.com> wrote:
>
>
> On Mon, Dec 6, 2010 at 11:56 PM, Fernando Perez <fperez.net at gmail.com>
> wrote:
>>
>> Hi Skipper,
>>
>> On Mon, Dec 6, 2010 at 2:34 PM, Skipper Seabold <jsseabold at gmail.com>
>> wrote:
>> > I'm wondering if anyone might have a look at my cython code that does
>> > matrix multiplication and see where I can speed it up or offer some
>> > pointers/reading.  I'm new to Cython and my knowledge of C is pretty
>> > basic based on trial and (mostly) error, so I am sure the code is
>> > still very naive.
>>
>> a few years ago I had a similar problem, and I ended up getting a very
>> significant speedup by hand-coding a very unsafe, but very fast pure C
>> extension just to compute these inner products.  This was basically a
>> replacement for dot() that would only work with double precision
>> inputs of compatible dimensions and would happily segfault with
>> anything else, but it ran very fast.  The inner loop is implemented
>> completely naively, but it still beats calls to BLAS (even linked with
>> ATLAS) for small matrix dimensions (my case was also up to ~ 15x15).
>>
>> I'm attaching the code in case you find it useful, please keep in mind
>> I haven't compiled it in years, so it may have bit-rotted a little.
>>
>
> Blas adds quite a bit of overhead for multiplying small matrices, but so
> does calling from python. For implementing Kalman filters it might be better
> to write a whole Kalman class so that operations can be combined at the c
> level.
>
> Skipper, what kind of Kalman filter are you trying to implement?
>

It's just a linear Gaussian filter.  I use it to get the loglikelihood
of a univariate ARMA series with exact initial conditions.  As it
stands it is fairly inflexible, but if I can make it fast I would like
to generalize it.

There is a fair amount of scratch work in here, and some attempts at
generalized state space models, but all the action for my purposes is
in KalmanFilter.loglike

http://bazaar.launchpad.net/~jsseabold/statsmodels/statsmodels-skipper/annotate/head%3A/scikits/statsmodels/tsa/kalmanf/kalmanf.py#L505

It's not terribly slow, but I have to maximize the likelihood using
numerical derivatives, so it's getting called quite a few times.  A
1000 observation ARMA(2,2) series takes about 5-6 seconds on my
machine with fmin_l_bfgs_b.

Skipper