# Sparse matrix implementation of the Kalman Filter?

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.

• I think this question would be better if you included the minimum amount of mathematics to describe the problem. Many people here are familiar with linear algebra, but not with the underlying Kalman filtering process. Describing the H matrix (whatever it may be), and the equations that involve it which you are trying to solve, should lead to a better answer. – Bill Barth Apr 16 '13 at 16:12
• You are perhaps right. However, Kalman filtering schemes are a large topic unto themselves. It would be too much to ask to someone to learn how Kalman Filters work from my question and from that devise an answer. This would be research paper level work (I assume so anyway). I think anyone who would be in a position to answer the question would not need additional details. – Bogdanovist Apr 17 '13 at 2:52

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.

• If I were you I'd start looking at EnKF for this problem. – Brian Borchers Apr 18 '13 at 3:49

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.)

• Hello Stephen. Thanks for joining scicomp. I understand that the question was quite general, so that a specific answer is impossible to give. However, also in view of the benefit of fellow visitors, you maybe can give more information on what your algorithm actually does and provide a stable link, eg. via the doi, to the reference. – Jan Jun 5 '13 at 16:10
• Just posting the abstract here would be a good way of "explaining" the link to Stack Exchange standards. – dmckee Jun 5 '13 at 19:10
• Whilst I agree with Jan, EE's who deal with KFs (like me!) generally don't read meteorological publications. The fact that the paper is recommended as an answer for a notoriously subtle problem is motivation enough for me to hunt it out. – Damien Jun 5 '13 at 23:02
• The ensemble Kalman ﬁlter (EnKF) can be interpreted in the context of linear regression theory. The ﬁlter equations are equivalent to the normal equations for a weighted least-squares estimate that minimizes a quadratic functional. Solving the normal equations is numerically unreliable and subject to large errors when the problem is ill-conditioned. A numerically reliable and efﬁcient algorithm is presented, based on the minimization of an alternative functional. The method relies on orthogonal rotations, is highly parallel and does not ‘square’ matrices to compute the analysis update. – Stephen Thomas Sep 5 '13 at 22:00

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.

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

• The abstract is: – Stephen Thomas Sep 5 '13 at 21:49
• Did you mean this to be a comment or edit to your other answer? – Jed Brown Sep 6 '13 at 7:23