Linear Filtering in High Dimensional State Space

Question

I am working with Gaussian Linear State Space models of the form:

$$y_t=F_t\Theta_t+v_t$$ $$\Theta_t=G_t\Theta_{t-1}+w_t$$ $$v_t \sim N(0, V_t)$$ $$w_t \sim N(0, W_t)$$

Where $y_t$ is my observed time-series (usually multivariate) and $\Theta_t$ is the state of the system at time $t$. $F_t$ and $G_t$ are the system and observation equations respectively.

I have in general been inferring these models using Kalman Filtering/Smoothing but I am now working with data/models that have very large state-spaces (~500 dimensions) and I find that the implementations of the Kalman Filter I have tried take too long to run.

Of note - the matrix $G_t$ and $F_t$ are often quite sparse in my models.

I am looking for advice regarding either extensions or alternatives to the Kalman filter that allow me to infer this type of Dynamic Linear Model. My initial thoughts are to use a reduced dimensional representation of the state space but I am neither sure of how to implement this or what exactly it would look like.

I would very much appreciate references or links to available code (R or C++ preferred).

Thank you!

Answer 1

A few thoughts:

1. Use the Woodbury-Matrix formula, if it helps, when you're computing the Kalman gain matrix. Usually that inverse matrix in there takes the longest to calculate. With this you might be able to write that inverse in terms of the inverse of a smaller matrix. This is why I asked about your observation dimension.

2. If $F_t$, $G_t$, $V_t$ and $W_t$ are block-diagonal, or sparse as you say, you might be able "decouple". I remember vaguely if that these are all block-diagonal, and if your first state distribution has uncorrelated components, they stay uncorrelated after every iteration of the Kalman filter. So you can just fit individual state space models to these components. Proving this is straightforward. You just write the Kalman filter equations in terms of block matrices, and check to make sure the off-diagonal blocks are $\mathbf{0}$ matrices.

3. Use the "information filter." I've never used it myself, but apparently you just keep track of the inverse covariance matrix, or precision matrix, at every time step. You never have to invert anything. But it has other drawbacks. I forget what they were.

Answer 2

You don't necessarily need to do anything beyond exploit the structure of your problem.

In particular, don't compute, or store or in any other way chew up time and space, dealing with the parts of the calculation that will just be zero. Figure out where the things you need to change will change and only store and update those. Also, some people over-compress the steps of the calculation ("look, I wrote it in one less line of algebra"), but in the process, can mask some of the opportunities to exploit redundant calculation. Some algorithms don't explicitly compute the Kalman gain as a step, for example, but there might be a speedup there (I've seen some algorithms implicitly computing it three times).

In addition there are some other obvious things to avoid

• don't multiply by 1.

• don't add zero.

I don't mean "check for this calculation in code and skip it" (that's almost useless) - I mean work out where you would have done this calculation using your knowledge of the structure of the problem and never go there in code at all.

Further:

• don't compute the same thing twice. (At worst, you could consider memoization, but much better to avoid even considering calculating what you know you'll already know; don't do in software what your brain can easily anticipate and you can simply avoid altogether)

  There's often a variety of symmetries to exploit as well. e.g. If you're already computing one thing, don't then turn around and compute its transpose; you already have it (this is a subset of computing the same thing twice).

  Often there can be some kind of block structure (or similar structures), sometimes with some efficiencies between (such as one block being closely related to another) or within blocks (e.g. symmetries or skew-symmetries)

Some time thinking about the avoidable calculations you're doing can make a huge difference. Do the biggest things first (I mean the things that will avoid a lot of calculation, not the things that take the most time to implement); you may find you don't need to get every ounce of gain when you can get 85% of it for less than half the work.

But by the same token, don't overoptimize to the extent that you're getting in the way of the things the hardware and numerical libraries can do well (this depends on the precise nature of your hardware and libraries). Some of this can be subtle, but don't knock yourself out over tiny gains (or even negative gains if you're not careful)