Implement Square Root Filter

[cstjean]

After running the Kalman filter on some data, I'd like to sample from the posterior, but I get

> rand(estimates[end])
ArgumentError: matrix is not symmetric/Hermitian. This error can be avoided by calling chol(Hermitian(A)) which will ignore either the upper or lower triangle of the matrix.

Stacktrace:
 [1] chol(::Array{Float64,2}) at ./linalg/cholesky.jl:185
 [2] rand(::GaussianDistributions.Gaussian{Array{Float64,1},Array{Float64,2}}) at /home/cst-jean/.julia/v0.6/GaussianDistributions/src/GaussianDistributions.jl:68
 [3] include_string(::String, ::String) at ./loading.jl:522

I admit that I don't understand these issues as well as I should. Should I just follow their advice, and use P.μ + chol(Hermitian(P.Σ))'*randn(T, length(P.μ))? It seems to work.

[mschauer]

That, or use convert_hermitian(G::Gaussian) = Gaussian(G.μ, Hermitian(G.Σ)).

A different solution is make the Arguments to the Kalman-Filter-Type which are covariances actually Hermitians, but I don't know if this would work in your code.

[cstjean]

Cool, thank you.

A different solution is make the Arguments to the Kalman-Filter-Type which are covariances actually Hermitians, but I don't know if this would work in your code.

I tried to make the prior Hermitian, but it didn't propagate to the posteriors. Is that what you meant?

[mschauer]

Ah, I see that this is not yet possible. I'll have to a square root filter anyway which alleviates the problem.

[touste]

Hi, I have a similar issue with some of my data. My Q and R matrices are Hermitian, but the problem still arises.

It happens as soon as R has off-diagonal terms.

Here is a MWE:

Φ = Matrix{Float64}(I, 6, 6)
for i in 1:3
	Φ[i,i+3] = 1.0
end

b = zeros(6)
Q = Matrix{Float64}(I, 6, 6) *0.1

E = LinearEvolution(Φ, Gaussian(b, Q))

H = zeros(3, 6)
for i in 1:3
	H[i,i] = 1.0
end

R = zeros(3, 3)
R[1,1] = R[2,2] = R[3,3] = 0.1
R[1,2] = R[2,1] = 0.02
R[1,3] = R[3,1] = 0.01

O = LinearObservation(E, H, R)

x0 = zeros(Float64, 6)
P0 = copy(Q)

Y = [i => [i*1.0, i*2.0, i*4.0] for i in 1:10]

Xf, ll = kalmanfilter(O, 0 => Gaussian(x0, P0), Y)

Error stack:

ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.
Stacktrace:
 [1] checkpositivedefinite(::Int64) at /usr/share/julia/stdlib/v1.3/LinearAlgebra/src/factorization.jl:18
 [2] #cholesky!#125(::Bool, ::typeof(cholesky!), ::Array{Float64,2}, ::Val{false}) at /usr/share/julia/stdlib/v1.3/LinearAlgebra/src/cholesky.jl:255
 [3] #cholesky! at ./none:0 [inlined] (repeats 2 times)
 [4] #cholesky#129 at /usr/share/julia/stdlib/v1.3/LinearAlgebra/src/cholesky.jl:348 [inlined]
 [5] cholesky at /usr/share/julia/stdlib/v1.3/LinearAlgebra/src/cholesky.jl:348 [inlined] (repeats 2 times)
 [6] whiten(::Array{Float64,2}, ::Array{Float64,1}) at .julia/packages/GaussianDistributions/bzffd/src/GaussianDistributions.jl:77
 [7] sqmahal(::Gaussian{Array{Float64,1},Array{Float64,2}}, ::Array{Float64,1}) at .julia/packages/GaussianDistributions/bzffd/src/GaussianDistributions.jl:86
 [8] logpdf(::Gaussian{Array{Float64,1},Array{Float64,2}}, ::Array{Float64,1}) at .julia/packages/GaussianDistributions/bzffd/src/GaussianDistributions.jl:94
 [9] dyniterate(::Bind{LinearObservation{LinearEvolution{Array{Float64,2},Gaussian{Array{Float64,1},Array{Float64,2}}},Array{Float64,2},Array{Float64,2}},Array{Pair{Int64,Array{Float64,1}},1}}, ::Tuple{Int64,Kalman.Filter{Pair{Int64,Gaussian{Array{Float64,1},Array{Float64,2}}}}}) at .julia/packages/Kalman/yfnWQ/src/filter.jl:90
 [10] kalmanfilter(::LinearObservation{LinearEvolution{Array{Float64,2},Gaussian{Array{Float64,1},Array{Float64,2}}},Array{Float64,2},Array{Float64,2}}, ::Pair{Int64,Gaussian{Array{Float64,1},Array{Float64,2}}}, ::Array{Pair{Int64,Array{Float64,1}},1}) at .julia/packages/Kalman/yfnWQ/src/filter.jl:149

[mschauer]

Hi @touste
You could define

 Kalman.llikelihood(yres, S) = Kalman.logpdf(Kalman.Gaussian(zero(yres), (S+S')/2), yres)

to force symmetry.

But do you need actually the likelihood? If not, I could make a version of the filter which does not compute the likelihood, then this problem does not show up. You could also define Kalman.llikelihood(yres, S) = NaN to avoid the likelihood computation.

[touste]

Hi, thanks for the tip, it works perfectly!

Yes I need the likelihood, but I wonder why is this definition not the default one in your package? It seems that there can be plenty of cases where S is not symmetric.

[mschauer]

master now follows your suggestion