JuliaGNSS / KalmanFilters.jl

Various Kalman Filters: KF, UKF, AUKF and their Square root variant
Other
41 stars 6 forks source link

Forcing positive diagonals when performing QR decomposition for SRKF #16

Open THargreaves opened 2 weeks ago

THargreaves commented 2 weeks ago

First off, thank you very much to the developers of this package! I really love the design and there is clear attention to detail.

I've been using the SRKF implementation and noticed a small issue. During calc_cross_cov_innovation_posterior the new matrices are computed using a QR decomposition. The Julia/LAPLACK implementations of QR decomposition do not enforce that R has positive entries on the diagonal. This means that when S, P_post are extracted from the block diagonals, they may not be valid L matrices for Cholesky decompositions (these must have positive diagonals).

This causes trouble when computing the incremental marginal likelihood via logpdf(MvNormal(ỹ, PDMat(S)), y). This is because logdet of Cholesky is computed as the some of logs of the diagonal elements (assuming they are positive, and no checks are put in place to ensure they actually are).

I've included an example below.

using Distributions
using KalmanFilters
using LinearAlgebra
using PDMats

x = [0.5623846194937526, 0.6019375204345232]
P = Cholesky{Float64,Matrix{Float64}}([1.2197182702938882 1.2659926501405483; 1.5441543654342051 0.15004572573480882], 'U', 0)
H = [0.3850845754724803 0.19625883331364058; 0.5445983221231325 0.2997283029235214]
y = [0.6250101801493309, 0.07868380451798473]
R = LinearAlgebra.Cholesky{Float64, Matrix{Float64}}([0.5689422644561306 0.05362293003382934; 0.030508351240219524 0.08648515775194665], 'U', 0)

ỹ = KalmanFilters.calc_innovation(H, x, y)
PHᵀ, S, P_post = KalmanFilters.calc_cross_cov_innovation_posterior(P, H, R, nothing)
logpdf(MvNormal(ỹ, PDMat(S)), y)
#> ERROR: DomainError with -0.9166852531274825:
#> log was called with a negative real argument but will only return a complex result if called with a complex argument. Try log(Complex(x)).

I believe this issue also exists with the time update and perhaps elsewhere in the SR implementations.

I'm happy to make changes to the code to fix these issues but I wanted to check in with the devs first. I'm also happy to just make adjustments in my code to "correct" the Cholesky decomposition to have positive diagonals if you don't want the (small) overhead of enforcing positivity.

zsoerenm commented 1 week ago

~~I'm not totally sure this should be part of KalmanFilters.jl. Kalman filter doesn't require the diagonal elements to be positive right, or should it? If you have a requirement for that I suggest that you call the functions on your own (like calc_upper_triangular_of_qr etc and enforce the diagonal elements to be positive). If you think they should be positive, I'm happy to change my mind.~~ Actually, I think you are right. They are cholesky decomposition, so they should be positive. So yes, we should ensure that they are positive. Would you mind opening a PR?