Sparse matrix implementation of the Kalman Filter?

Question

I have a Kalman Filter based modelling code that I have developed for a near-real time regional ionospheric mapping application. The code assimilates data from different sensors into a map (described by a set of basis functions) using a Kalman Filter.

I am trying to scale this up to a larger region and more sensors, however the matrix algebra part of the Kalman Filter is becoming very slow, due to the large matrices (thousands of rows/columns) involved. I suspect the best way to attack the runtime issue is to use the fact that these matrices are typically very sparse with 80% or more of the total elements zero. The reason for this is that each sensor has a bias parameter that is jointly estimated with the map coefficients. This shows up as a 1 in the column for that sensor in the Kalman H matrix, with zero in the columns for every other sensor and map co-efficient. There are hundreds of sensors each contributing 8-10 observations at each epoch, hence a lot of zeros.

I could look at implementing the components of the Kalman filter using sparse algorithms, specifically multiplication and inversion*, but I wonder if there is an even better approach that re-casts the Kalman filter in a different form more suitable for cases when the matrices are sparse? I know I could use an ensemble Kalman filter or something similar, but if possible I'd like to retain the optimality of the pure linear Kalman filter; the total data volume is not prohibitive, just the large sparse matrices that result from the linear model.

In terms of implementation, this is done in IDL, however the core matrix algebra is done via calls to external optimised LA libraries (specifically ATLAS).

*I know that an optimum Kalman filter implementation avoids inversion and instead uses a UD decomposition. I am considering trying to implement something like this, so that may be the answer, but I'm fishing for whether there is a better solution given the sparseness of the matrices.

Answer 1

With sparse matrices, it's frequently the case that although a matrix $A$ is sparse, $A^{-1}$ is dense. In such cases the Cholesky or $LU$ factorization of $A$ is more likely to be sparse (particularly if the rows/columns of $A$ are reordered to improve the sparsity pattern.) In most cases, if you want to exploit sparsity and aren't interested in using an iterative algorithm for solving systems involving $A$, then you're better off using some factorization of the matrix rather than explicitly computing $A^{-1}$.

For Kalman filtering in particular, rather than computing

$S_{k}=(H_{k}P_{k-1,k}H_{k}^{T}+R_{k})^{-1}$

you'd typically be better off working working with a factorization of $S_{k}^{-1}$. Since $S_{k}$ is symmetric and should be positive definite, you can use a Cholesky factorization or $LDL^{T}$ factorization to do this. You've told us that your $H_{k}$ matrix is sparse, but you haven't told us anything about whether $R_{k}$ is sparse or otherwise structured, and of course $P$ could be quite dense.

One reason that the Ensemble Kalman Filter (EnKF) and various particle filtering techniques are so popular is that for systems with an extremely large state vector, conventional Kalman filtering becomes very difficult. EnKF can be efficiently implemented for very large state vectors if $R_{k}$ is diagonal or nearly diagonal. These questions have been dealt with in great depth by people working in the field of data assimilation, so I'd suggest starting your research by reading about how they've dealt with these issues.

Answer 2

We have a robust algorithm for the Ensemble Kalman (and regular Kalman) filter. It is well suited to sparse matrices and parallel computation because it is based on orthogonal matrices and it related to the square root or UD algorithms.

Would be glad to send the paper

Thomas, S. J., J. Hacker and J. Anderson, (2009): A robust formulation of the ensemble Kalman filter Quart J. Royal Met. Soc, vol 135, 507-521,

(PDF from the publisher is free.)

Answer 3

Long time back I had a chance to work on dimensionality reduction which deals with processing data that falls in large sets. The basic idea behind it is that it processes the data through a few steps to orient it in the way that most information can be calculated from it.

It works pretty fine for matrices as well and is used largely. The best part is that you need not even program it as there are standard libraries available for it already. Major mathematical tools like Matlab and Mathematica also support this functionality straightforwardly.

There are two main algorithms that achieve this - Principal Component Analysis and Singular Value Decomposition.

What these algorithms actually achieve is finding that data that actually affects your reading by a significant margin. The internet is full of information on these algorithms. This will show you how Apache is doing it.

Answer 4

Apologies for not adding to the discussion for a while - but I would be glad to post the abstract of the paper - but also gladly go into the reasons why the computations are organized differently than standard KF of EnKF approaches