Skip to content

Commit

Permalink
Merge pull request #37 from fipelle/dev
Browse files Browse the repository at this point in the history 
Added support for UniformScaling R (useful for large dimensionality reduction problems)
  • Loading branch information
@fipelle
fipelle authored Mar 15, 2022
2 parents 8b6019a + 9035362 commit 1182033
Show file tree
Hide file tree
Showing 3 changed files with 59 additions and 22 deletions.
55 changes: 44 additions & 11 deletions src/kalman.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,19 +168,55 @@ function find_observed_data(settings::KalmanSettings, t::Int64)
end

"""
update_P_post!(P_post_old::SymMatrix, status::KalmanStatus, K_t::FloatMatrix, R_t::SubArray{Float64})
update_P_post!(P_post_old::Nothing, status::KalmanStatus, K_t::FloatMatrix, R_t::SubArray{Float64})
update_inv_F!(R::UniformScaling{Float64}, status::KalmanStatus, B_t::SubArray{Float64}, ind_not_missings::IntVector)
update_inv_F!(R::SymMatrix, status::KalmanStatus, B_t::SubArray{Float64}, ind_not_missings::IntVector)
Update status.inv_F.
"""
function update_inv_F!(R::UniformScaling{Float64}, status::KalmanStatus, B_t::SubArray{Float64}, ind_not_missings::IntVector)
F_t = status.buffer_m_n_obs'*B_t'; # I cannot use mul!(...) here since R is UniformScaling{Float64}
F_t += R;
status.inv_F = inv(Symmetric(F_t));
end

function update_inv_F!(R::SymMatrix, status::KalmanStatus, B_t::SubArray{Float64}, ind_not_missings::IntVector)
F_t = R[ind_not_missings, ind_not_missings];
mul!(F_t, status.buffer_m_n_obs', B_t', 1.0, 1.0);
status.inv_F = inv(Symmetric(F_t));
end

"""
update_P_post!(P_post_old::SymMatrix, R_t::UniformScaling{Float64}, status::KalmanStatus, K_t::FloatMatrix, ind_not_missings::IntVector)
update_P_post!(P_post_old::Nothing, R_t::UniformScaling{Float64}, status::KalmanStatus, K_t::FloatMatrix, ind_not_missings::IntVector)
update_P_post!(P_post_old::SymMatrix, R::SymMatrix, status::KalmanStatus, K_t::FloatMatrix, ind_not_missings::IntVector)
update_P_post!(P_post_old::Nothing, R::SymMatrix, status::KalmanStatus, K_t::FloatMatrix, ind_not_missings::IntVector)
Update status.P_post.
"""
function update_P_post!(P_post_old::SymMatrix, status::KalmanStatus, K_t::FloatMatrix, R_t::SubArray{Float64})
function update_P_post!(P_post_old::SymMatrix, R_t::UniformScaling{Float64}, status::KalmanStatus, K_t::FloatMatrix, ind_not_missings::IntVector)
mul!(status.buffer_m_m, status.L, status.P_prior);
mul!(status.P_post.data, status.buffer_m_m, status.L');
mul!(status.buffer_m_n_obs, K_t, R_t);
mul!(status.P_post.data, status.buffer_m_n_obs, K_t', 1.0, 1.0);
end

function update_P_post!(P_post_old::Nothing, R_t::UniformScaling{Float64}, status::KalmanStatus, K_t::FloatMatrix, ind_not_missings::IntVector)
mul!(status.buffer_m_m, status.L, status.P_prior);
status.P_post = Symmetric(status.buffer_m_m*status.L');
mul!(status.buffer_m_n_obs, K_t, R_t);
mul!(status.P_post.data, status.buffer_m_n_obs, K_t', 1.0, 1.0);
end

function update_P_post!(P_post_old::SymMatrix, R::SymMatrix, status::KalmanStatus, K_t::FloatMatrix, ind_not_missings::IntVector)
R_t = @view R[ind_not_missings, ind_not_missings];
mul!(status.buffer_m_m, status.L, status.P_prior);
mul!(status.P_post.data, status.buffer_m_m, status.L');
mul!(status.buffer_m_n_obs, K_t, R_t);
mul!(status.P_post.data, status.buffer_m_n_obs, K_t', 1.0, 1.0);
end

function update_P_post!(P_post_old::Nothing, status::KalmanStatus, K_t::FloatMatrix, R_t::SubArray{Float64})
function update_P_post!(P_post_old::Nothing, R::SymMatrix, status::KalmanStatus, K_t::FloatMatrix, ind_not_missings::IntVector)
R_t = @view R[ind_not_missings, ind_not_missings];
mul!(status.buffer_m_m, status.L, status.P_prior);
status.P_post = Symmetric(status.buffer_m_m*status.L');
mul!(status.buffer_m_n_obs, K_t, R_t);
Expand Down Expand Up @@ -224,9 +260,8 @@ Kalman filter a-posteriori update. All measurements are not observed at time t.
"""
function aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::IntVector)

# Convenient views
# Convenient view
B_t = @view settings.B[ind_not_missings, :];
R_t = @view settings.R[ind_not_missings, ind_not_missings];

# Forecast error
status.e = settings.Y.data[ind_not_missings, status.t];
Expand All @@ -237,9 +272,7 @@ function aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_mi
status.buffer_m_n_obs = status.P_prior*B_t';

# Inverse of the forecast error covariance matrix
F_t = settings.R[ind_not_missings, ind_not_missings];
mul!(F_t, status.buffer_m_n_obs', B_t', 1.0, 1.0);
status.inv_F = inv(Symmetric(F_t));
update_inv_F!(settings.R, status, B_t, ind_not_missings);

# Kalman gain
K_t = status.buffer_m_n_obs*status.inv_F;
Expand All @@ -253,8 +286,8 @@ function aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_mi
mul!(status.X_post, K_t, status.e, 1.0, 1.0);

# A posteriori estimates: P_post (P_post is updated using the Joseph form)
update_P_post!(status.P_post, status, K_t, R_t);

update_P_post!(status.P_post, settings.R, status, K_t, ind_not_missings);
# Update log likelihood
if settings.compute_loglik == true
update_loglik!(status);
Expand Down
14 changes: 7 additions & 7 deletions src/types.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ where ``e_{t} \\sim N(0_{nx1}, R)`` and ``U_{t} \\sim N(0_{mx1}, Q)``.
struct KalmanSettings
Y::MSeries
B::FloatMatrix
R::SymMatrix
R::Union{UniformScaling{Float64}, SymMatrix}
C::FloatMatrix
D::FloatMatrix
Q::SymMatrix
Expand All @@ -81,7 +81,7 @@ KalmanSettings constructors
=#

"""
KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, Q::SymMatrix; kwargs...)
KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::Union{UniformScaling{Float64}, SymMatrix}, C::FloatMatrix, Q::SymMatrix; kwargs...)
`KalmanSettings` constructor.
Expand All @@ -99,7 +99,7 @@ KalmanSettings constructors
# Notes
This particular constructor sets `D` to be an identity matrix, `X0` to be a vector of zeros and computes `P0` via `solve_discrete_lyapunov(C, Q)`.
"""
function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, Q::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)
function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::Union{UniformScaling{Float64}, SymMatrix}, C::FloatMatrix, Q::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)

# Compute default value for missing parameters
n, T = size(Y);
Expand All @@ -113,7 +113,7 @@ function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix,
end

"""
KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix; kwargs...)
KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::Union{UniformScaling{Float64}, SymMatrix}, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix; kwargs...)
`KalmanSettings` constructor.
Expand All @@ -132,7 +132,7 @@ end
# Notes
This particular constructor sets `X0` to be a vector of zeros and computes `P0` via `solve_discrete_lyapunov(C, Q)`.
"""
function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)
function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::Union{UniformScaling{Float64}, SymMatrix}, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)

# Compute default value for missing parameters
n, T = size(Y);
Expand All @@ -146,7 +146,7 @@ function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix,
end

"""
KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, X0::FloatVector, P0::SymMatrix; kwargs...)
KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::Union{UniformScaling{Float64}, SymMatrix}, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, X0::FloatVector, P0::SymMatrix; kwargs...)
`KalmanSettings` constructor.
Expand All @@ -164,7 +164,7 @@ end
- `compute_loglik`: Boolean (`true` for computing the loglikelihood in the Kalman filter - default: `true`)
- `store_history`: Boolean (`true` to store the history of the filter and smoother - default: `true`)
"""
function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, X0::FloatVector, P0::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)
function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::Union{UniformScaling{Float64}, SymMatrix}, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, X0::FloatVector, P0::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)

# Compute default value for missing parameters
n, T = size(Y);
Expand Down
12 changes: 8 additions & 4 deletions test/kalman.jl
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,10 @@
"""
ksettings_input_test(ksettings::KalmanSettings, Y::JArray, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, X0::FloatArray, P0::SymMatrix, DQD::SymMatrix, n::Int64, T::Int64, m::Int64; compute_loglik::Bool=true, store_history::Bool=true)
ksettings_input_test(ksettings::KalmanSettings, Y::JArray, B::FloatMatrix, R::Union{UniformScaling{Float64}, SymMatrix}, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, X0::FloatArray, P0::SymMatrix, DQD::SymMatrix, n::Int64, T::Int64, m::Int64; compute_loglik::Bool=true, store_history::Bool=true)
Return true if the entries of ksettings are correct (false otherwise).
"""
function ksettings_input_test(ksettings::KalmanSettings, Y::JArray, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, X0::FloatArray, P0::SymMatrix, DQD::SymMatrix, n::Int64, T::Int64, m::Int64; compute_loglik::Bool=true, store_history::Bool=true)
function ksettings_input_test(ksettings::KalmanSettings, Y::JArray, B::FloatMatrix, R::Union{UniformScaling{Float64}, SymMatrix}, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, X0::FloatArray, P0::SymMatrix, DQD::SymMatrix, n::Int64, T::Int64, m::Int64; compute_loglik::Bool=true, store_history::Bool=true)

return ~(false in [ksettings.Y.data === Y;
ksettings.B == B;
ksettings.R == R;
Expand Down Expand Up @@ -85,11 +86,11 @@ function test_kalman_output(ksettings::KalmanSettings, kstatus::SizedKalmanStatu
end

"""
kalman_test(Y::JArray, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, benchmark_data::Tuple)
kalman_test(Y::JArray, B::FloatMatrix, R::Union{UniformScaling{Float64}, SymMatrix}, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, benchmark_data::Tuple)
Run a series of tests to check whether the kalman.jl functions work.
"""
function kalman_test(Y::JArray, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, benchmark_data::Tuple)
function kalman_test(Y::JArray, B::FloatMatrix, R::Union{UniformScaling{Float64}, SymMatrix}, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, benchmark_data::Tuple)

# Benchmark data
benchmark_n, benchmark_T, benchmark_m, benchmark_X0, benchmark_P0, benchmark_DQD, benchmark_X_prior, benchmark_P_prior,
Expand Down Expand Up @@ -199,6 +200,7 @@ end

# Run tests
kalman_test(Y, B, R, C, D, Q, benchmark_data);
kalman_test(Y, B, 1e-4I, C, D, Q, benchmark_data); # same test with UniformScaling R
end

@testset "multivariate model" begin
Expand Down Expand Up @@ -249,6 +251,7 @@ end

# Run tests
kalman_test(Y, B, R, C, D, Q, benchmark_data);
kalman_test(Y, B, 1e-4I, C, D, Q, benchmark_data); # same test with UniformScaling R
end

@testset "multivariate model (non-diagonal)" begin
Expand Down Expand Up @@ -300,4 +303,5 @@ end

# Run tests
kalman_test(Y, B, R, C, D, Q, benchmark_data);
kalman_test(Y, B, 1e-4I, C, D, Q, benchmark_data); # same test with UniformScaling R
end

0 comments on commit 1182033

Please sign in to comment.