search

TuringLang / GeneralisedFilters.jl

Filtering and smoothing algorithms for state space models with analytic (or approximately/partially analytic) solutions
MIT License
0 stars 0 forks source link

Non-Hermitian Innovations Covariance #2

Open charlesknipp opened 1 month ago

charlesknipp commented 1 month ago

Dependent on RNG, the Kalman filter may fail for the RBPF test case. You can replicate this issue by using MersenneTwister instead of StableRNG.

The filtering algorithm raises a PosDefException when evaluating the log likelihood. This is caused by non-symmetry in the innovations covariance S. A quick fix would be to deploy the following:

S = LinearAlgebra._hermitianpart!(H * Σ * H') + R
K = Σ * H' / cholesky(S)

but this problem extends to other covariance matrices. So it may be worth investigating other instances which potentially fail a Cholesky decomposition.

On a semi-related note, it is common for some models (particularly in macroeconomics) to have rank deficient covariance matrices. These will also raise errors when taking a Cholesky decomposition. While this is not necessary for the Kalman filter to run, this will fail to generate an MvNormal for the state transition density.

charlesknipp commented 1 month ago

I usually catch the main issue by liberal use of PDMats, which is admittedly inefficient for particle methods. You can see how I deal with both problems here, but that makes heavy use of the PDMat.

devmotion commented 1 month ago

In my experience it's best to use and update a factorization in the Kalman filter, see eg https://en.wikipedia.org/wiki/Kalman_filter#Factored_form.

THargreaves commented 1 month ago

Ah yes. We've ran into this issue in the past and had just used the hack S = (S + S') / 2. Using LinearAlgebra._hermitianpart! looks much cleaner.

The rank-deficient covariance matrices is not something I had considered, but is definitely something worth thinking about in the long-run.

I soon plan to implement a square-root variant of the Kalman Filter as devmotion suggests so that should at least provide a stable (if a bit slower) option.