Simplest "Square Root" Kalman Smoother

Question

Hi everyone,

I've been working on a Kalman Filter/Smoother library for python, and I was looking for a "simple" Kalman Smoothing algorithm that allows me to propgate the Cholesky factorization of the covariance matrix rather than the whole thing. The benefit of doing so is to avoid numerical instabilities causing small eigenvalues to turn negative. I've seen a few papers from the 70s and 80s describing so-called "Square Root" Kalman Smoothers, but each of them involves some pretty complicated factorization that would make my code unreadable. My ultimate goal is to have something both easy to use and easy to read, so I ask the community: do you know of any concise Kalman Smoothing update equations?

EDIT: Barring a better solution, I've decided to use the the Bierman-Thornton U-D Filter and Smoother. It unfortunately is not a very clear algorithm, involving an uncommon matrix decomposition (UDU for upper triangular U and diagonal D) and requires I process each coordinate of the observed vector one at a time, but it seems reasonably well known and documented.

Answer 1

Hi, I have implemented the Algorithm described described in "Fixed-Interval Smoothing Algorithm Based on Singular Value Decomposition", Youmin Zhang and X. Rong Li, 1996, IEEE International Converence on Control Applications If you want I can provide the cleaned up Algorithm as a PDF Document. However, the Algorithm performs poor when the Measurement Noise is low. Numerical errors are the reason. I could improved by substituting the formula for the forward Kalman-gain by an equivalent one. However, does the achieved precision not exceed that of the simple RTS-Smoother. It seems to me that the Algorithm does not hold its promise.

Recently an alternative has appeared that I would like to try out one day, it is given in "Square Root Modiﬁed Bryson–Frazier Smoother", Richard G. Gibbs, IEEE Transactions on Automation and Control, 2011. As a great advantage it allows singular measurement noise covariance R. However it seems to end up in a complex algorithm. Does anyone have written down this Algorithm?

Answer 2

Hello,

The square-root Bryson-Frazier smoother seems simple:

x <-- x_initial
P <-- P_initial
# Note that each covariance matrix is lowerly triangulated square-root one.
for y in y_series:
    # Calculate the lowerly-triangulated matrix by using (Householder) orthogonalization as follows:
        |     R~    0     0  |    | R  H*P    0 |
        | F*K*R~    P~    0  | =  | 0  F*P  G*Q | * (orthogonal matrix)
        |   K*R~  C*P~   (*) |    | 0    P    0 |
    Calculate K and C by using back-substitution.
    e = y - H*x
    x <-- x + K*e
    Save x, K*e, and C in filter_series.
    x <-- F*x
    P <-- P~
lam = 0
for x, K*e, C in reversed(filter_series):
    x_smooth = x - C*lam
    lam = C*lam - K*e

If you want to avoid the back-substitutions, see @article{ParkPetKailathT1995a, year=1995, volume=40, number=4, pages="761--766", author="Park, P. and Kailath, T.", title="Square-Root {Bryson}-{Frazier} Smoothing Algorithms", journal=IEEE_J_AC}

Answer 3

Probably helpful: http://infoscience.epfl.ch/record/161468/files/cholupdate.pdf?version=1

Simplest "Square Root" Kalman Smoother

Subscription:

Related questions