Add sparse kalman filter

[guilhermebodin]

@raphaelsaavedra, There is a project going on where @AndreFCR and @marinadietze have a system with ~200 states. It currently takes about 50 minutes to estimate the system for 25k observations with some relaxation on the gradient norm. Let's see how low they can go.

[codecov]

Codecov Report

Merging #290 (cc6a454) into master (2ca5237) will decrease coverage by 1.22%.
The diff coverage is 61.98%.

@@            Coverage Diff             @@
##           master     #290      +/-   ##
==========================================
- Coverage   93.88%   92.66%   -1.23%     
==========================================
  Files          33       34       +1     
  Lines        3043     3164     +121     
==========================================
+ Hits         2857     2932      +75     
- Misses        186      232      +46     

Impacted Files	Coverage Δ
src/filters/sparse_univariate_kalman_filter.jl	61.98% <61.98%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 2ca5237...cc6a454. Read the comment docs.

[raphaelsaavedra]

That's really nice!

Looking at this PR, the main thing that comes to my mind is there seems to be a lot of duplicate code from the non-sparse types. Is it really necessary to have a new type specific for the sparse state? First, it looks that the types in SparseUnivariateKalmanState aren't sparse, just some calculations that happen in the filter, and second, if we want to allow SparseVector and SparseMatrixCSC, we could relax the types of the non-sparse state types to AbstractVector and AbstractMatrix.

But these are just my thoughts after taking a quick look, it might not be so simple.

[guilhermebodin]

@raphaelsaavedra agreed, I am going to try to relax the type on the StateSpaceSystems and KalmanStates. My main concern was that not all AbstractArrays work with LinearAlgebra.BLAS.gemm! but mul! is works with every type and it is just as fast for non sparse matrices.

using BenchmarkTools
A = rand(200, 200)
B = spzeros(200, 200)
C = similar(A)
B[1:100] .= 1
Fl = Float64


@btime mul!($C, $A, $(B'));
B_prime = zeros(200, 200)
B_prime[1:200] .= 1
@btime mul!($C, $A, $(B_prime'));
@btime LinearAlgebra.BLAS.gemm!('N', 'T', one(Fl), $A, $B_prime, zero(Fl), $C);

[guilhermebodin]

Well I gave it a try and the implementation to consider every vector ended up in much more allocations and bigger compile times

julia> @time begin
           y = randn(100)
           model = BasicStructural(y, 12)
           fit!(model)
       end
  3.616171 seconds (115.63 M allocations: 1.942 GiB, 4.11% gc time)
BasicStructural

for comparison, this is the one in master (and the one with the sparse filter on a separated file)

julia> @time begin
           y = randn(100)
               model = BasicStructural(y, 12)
           fit!(model)
       end
  0.157086 seconds (20.66 k allocations: 1.698 MiB)
BasicStructural

I think we are going with the sparse filter in a separated file with a lot of duplicated code for now

[marinadietze]

@guilhermebodin, @AndreFCR tried to run our model (the one with trend, slope and two stochastic seasonalities with s = 24 and s = 168) and the running time with the sparse kalman filter was considerably higher. We will try to investigate more why this is happening and return to you.

[guilhermebodin]

@marinadietze the speed ups were identified on julia 1.6