Non-Hermitian Innovations Covariance

[charlesknipp]

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]

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]

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]

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.

[THargreaves]

This will largely be resolved by #30 and #36. All that is needed to close this issue is to replace the calls to (S + S') / 2 with this cleaner LinearAlgebra._hermitianpart!.

[THargreaves]

I plan to fix this issue in the next few days. My proposal is as follows that I would appreciate any thoughts on.

This would be combined with an implementation of the square-root Kalman filter for if high numerical stability is required.

1. We add a joseph_form::bool field to KalmanFilter.

If this is true (default false), we use the more numerically stable Joseph form formula for the posterior covariance instead of the optimised version we currently use.

The Joseph form is guaranteed to be symmetric by definition and I think it is more likely to be PD (but it could maybe still fail)

2. We add a pd_correct_strategy::PDCorrectStrategy field to KalmanFilter.

This determines what to do when Cholesky fails. I envision this being in the form of a try-catch so that we don't bother doing a correction unless the cholesky first fails, saving compute most of the time. The options would be

a) fail (default) — raise an error if Cholesky fails

b) diagonal loading — force symmetry, then add a small unit of the identity (fast but not well motivated)

c) eigenvalue thresholding — compute the eigen decomp and threshold all eigenvalues at some size (better motivated but quite slow)

For (c) ideally we would want to threshold at 0.0 since there's nothing wrong with a semi-definite covariance, but we'd then have to using PDMats to do the Cholesky

julia> cholesky(zeros(2, 2))
ERROR: PosDefException: matrix is not positive definite; Factorization failed.