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

[Skipper Seabold]

On Tue, Dec 7, 2010 at 7:37 PM, Charles R Harris
<charlesr.harris at gmail.com> wrote:
>
>
> On Tue, Dec 7, 2010 at 5:09 PM, <josef.pktd at gmail.com> wrote:
>>
>> On Tue, Dec 7, 2010 at 6:45 PM, Charles R Harris
>> <charlesr.harris at gmail.com> wrote:
>> >
>> >
>> > On Tue, Dec 7, 2010 at 10:39 AM, Skipper Seabold <jsseabold at gmail.com>
>> > wrote:
>> >>
>> >> On Tue, Dec 7, 2010 at 12:17 PM, Charles R Harris
>> >> <charlesr.harris at gmail.com> wrote:
>> >> >
>> >> >
>> >> > On Tue, Dec 7, 2010 at 10:05 AM, <josef.pktd at gmail.com> wrote:
>> >>
>> >> <snip>
>> >>
>> >> >> 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.)
>> >> >>
>> >> >
>> >> > The unscented Kalman filter is a better way to estimate the
>> >> > covariance
>> >> > of a
>> >> > non-linear process, think of it as a better integrator. If the
>> >> > propagation
>> >> > is easy to compute, which seems to be the case here, it will probably
>> >> > save
>> >> > you some time. You might even be able to use the basic idea and skip
>> >> > the
>> >> > Kalman part altogether.
>> >> >
>> >> > My general aim here is to optimize the algorithm first before getting
>> >> > caught
>> >> > up in the details of matrix multiplication in c. Premature
>> >> > optimization
>> >> > and
>> >> > all that.
>> >> >
>> >>
>> >> Hmm I haven't seen this mentioned much in what I've been reading or
>> >> the documentation on existing software for ARMA processes, so I never
>> >> thought much about it.  I will have a closer look.  Well, google turns
>> >> up this thread...
>> >>
>> >
>> > I've started reading up a bit on what you are doing and the application
>> > doesn't use extended Kalman filters, so the suggestion to use unscented
>> > Kalman filters is irrelevant. Sorry about that ;) I'm still wading
>> > through
>> > the various statistical notation thickets to see if there might be a
>> > better
>> > form to use for the problem but I don't see one at the moment.
>>
>> There are faster ways to get the likelihood for a simple ARMA process
>> than using a Kalman Filter. I think the main advantage and reason for
>> the popularity of Kalman Filter for this is that it is easier to
>> extend. So using too many tricks that are specific to the simple ARMA
>> might take away much of the advantage of getting a fast Kalman Filter.
>>
>> I didn't read much of the details for the Kalman Filter for this, but
>> that was my conclusion from the non-Kalman Filter literature.
>>

That's the idea.  The advantages of the KF are that it's inherently
structural and it's *very* general.  The ARMA case was just a jumping
off point, but has also proved to be a sticking point.  I'd like to
have a fast and general linear Gaussian KF available for larger state
space models, as it's the baseline workhorse for estimating linearized
large macroeconomic models at the moment.

>
> Well, there are five forms of the standard Kalman filter that I am somewhat
> familiar with and some are better suited to some applications than others.
> But at this point I don't see that there is any reason not to use the common
> form for the ARMA case. It would be interesting to see some profiling since
> the matrix inversions are likely to dominate as the number of variables go
> up.
>

If interested, I am using Durbin and Koopman's "Time Series Analysis
by State Space Methods" and Andrew Harvey's "Forecasting, Structural
Time Series Models, and the Kalman Filter" as my main references for
this.  The former is nice and concise but has a lot of details,
suggestions, and use cases.

I have looked some more and it does seem that the filter converges to
its steady state after maybe 2 % of the iterations depending on the
properties of the series, so for the ARMA case I can switch to the
fast recursions only updating the state (not quite sure on the time
savings yet), but I am moving away from my goal of a fast and general
KF implementation...

About the matrix inversions, the ARMA model right now is only
univariate, so there is no real inverting of matrices.  The suggestion
of Durbin and Koopman for larger, multivariate cases is to split it
into a series of univariate problems in order to avoid inversion.
They provide in the book some benchmarks on computational efficiency
in terms of multiplications needed based on their experience writing
http://www.ssfpack.com/index.html.

Skipper