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

[Skipper Seabold]

On Tue, Dec 7, 2010 at 12:05 PM,  <josef.pktd at gmail.com> wrote:
> On Tue, Dec 7, 2010 at 11:53 AM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
>>
>>
>> On Tue, Dec 7, 2010 at 9:47 AM, Skipper Seabold <jsseabold at gmail.com> wrote:
>>>
>>> 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.
>>>
>>
>> Just a guess here, but the numerical derivative bit makes it sounds like you
>> are implementing a generalized Kalman filter. Have you looked at unscented
>> Kalman filters?
>
> It's still a linear filter, non-linear optimization comes in because
> the exact loglikelihood function for ARMA is non-linear in the
> coefficients.
> (There might be a way to calculate the derivative in the same loop,
> but that's a different issue.)
>

Right.  The derivative is of the whole likelihood function with
respect to the parameters that make up the system matrices in the
state equation.  You can calculate the derivatives as part of the
recursions but the literature I've seen suggests that numerical
derivatives of the state equation matrices are the way to go.

Skipper