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.