Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Forcing positive diagonals when performing QR decomposition for SRKF #16

Closed
THargreaves opened this issue Oct 3, 2024 · 2 comments · Fixed by #17
Closed

Forcing positive diagonals when performing QR decomposition for SRKF #16

THargreaves opened this issue Oct 3, 2024 · 2 comments · Fixed by #17

Comments

@THargreaves
Copy link
Contributor

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
Copy link
Member

zsoerenm commented Oct 7, 2024

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?

@zsoerenm
Copy link
Member

Fixed in #17

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants