EnsembleKalmanProcesses.jl icon indicating copy to clipboard operation
EnsembleKalmanProcesses.jl copied to clipboard

SparseInversion not as robust as Inversion

Open ilopezgp opened this issue 2 years ago • 5 comments

I tried using the SparseInversion algorithm instead of Inversion for an integration test in CEDMF.jl, and I run into a SingularException in L111.

The configuration used for the integration test is

γ = 1.0
threshold_eki = true
threshold_value = 5e-2
reg = 1e-2

In general, we are missing documentation for how to choose reg and how to make the convex optimization part of the update robust.

Stack trace:

Stacktrace:
--
  | [1] checknonsingular
  | @ /central/software/julia/1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/factorization.jl:19 [inlined]
  | [2] checknonsingular
  | @ /central/software/julia/1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/factorization.jl:21 [inlined]
  | [3] #lu!#146
  | @ /central/software/julia/1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/lu.jl:82 [inlined]
  | [4] #lu#153
  | @ /central/software/julia/1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/lu.jl:279 [inlined]
  | [5] lu (repeats 2 times)
  | @ /central/software/julia/1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/lu.jl:278 [inlined]
  | [6] \(A::Matrix{Float64}, B::Diagonal{Float64, Vector{Float64}})
  | @ LinearAlgebra /central/software/julia/1.7.0/share/julia/stdlib/v1.7/LinearAlgebra/src/generic.jl:1142
  | [7] sparse_eki_update(ekp::EnsembleKalmanProcess{Float64, Int64, SparseInversion{Float64, Int64}}, u::Matrix{Float64}, g::Matrix{Float64}, y::Matrix{Float64}, obs_noise_cov::Matrix{Float64})
  | @ EnsembleKalmanProcesses /central/scratch/esm/slurm-buildkite/calibrateedmf-ci/depot/cpu/packages/EnsembleKalmanProcesses/OjLuD/src/SparseEnsembleKalmanInversion.jl:111
  | [8] (::EnsembleKalmanProcesses.var"#failsafe_update#26")(ekp::EnsembleKalmanProcess{Float64, Int64, SparseInversion{Float64, Int64}}, u::Matrix{Float64}, g::Matrix{Float64}, y::Matrix{Float64}, obs_noise_cov::Matrix{Float64}, failed_ens::Vector{Int64})
  | @ EnsembleKalmanProcesses /central/scratch/esm/slurm-buildkite/calibrateedmf-ci/depot/cpu/packages/EnsembleKalmanProcesses/OjLuD/src/SparseEnsembleKalmanInversion.jl:24
  | [9] update_ensemble!(ekp::EnsembleKalmanProcess{Float64, Int64, SparseInversion{Float64, Int64}}, g::Matrix{Float64}; cov_threshold::Float64, Δt_new::Float64, deterministic_forward_map::Bool, failed_ens::Nothing)
  | @ EnsembleKalmanProcesses /central/scratch/esm/slurm-buildkite/calibrateedmf-ci/depot/cpu/packages/EnsembleKalmanProcesses/OjLuD/src/SparseEnsembleKalmanInversion.jl:191

ilopezgp avatar Mar 29 '22 00:03 ilopezgp

The value of γ may be too small in this case? When the sparsity constraint is too strong, we would quickly have collapsed ensemble that could trigger singularity issue in convex optimization. I'll investigate this issue a little bit more to see if we can better guarantee the robustness.

jinlong83 avatar Mar 29 '22 07:03 jinlong83

The value of γ may be too small in this case? When the sparsity constraint is too strong, we would quickly have collapsed ensemble that could trigger singularity issue in convex optimization. I'll investigate this issue a little bit more to see if we can better guarantee the robustness.

That would be great, thanks! Is there any guideline we could document for the value of reg, or any heuristic to choose its value?

ilopezgp avatar Mar 29 '22 15:03 ilopezgp

The value of γ may be too small in this case? When the sparsity constraint is too strong, we would quickly have collapsed ensemble that could trigger singularity issue in convex optimization. I'll investigate this issue a little bit more to see if we can better guarantee the robustness.

That would be great, thanks! Is there any guideline we could document for the value of reg, or any heuristic to choose its value?

The value of reg corresponds to ensemble inflation and it would be ideal to keep its value small. The value of γ needs a little bit more trial and error, as we want to have a γ that leads to a sparse solution but would not significantly increase the original loss (which can be visualized in more details by checking the agreement with observation data).

jinlong83 avatar Mar 30 '22 03:03 jinlong83

Added a small note in the documentation for this https://clima.github.io/EnsembleKalmanProcesses.jl/dev/ensemble_kalman_inversion/#Sparsity-Inducing-Ensemble-Kalman-Inversion

odunbar avatar Sep 21 '22 21:09 odunbar

As a continuation here, we have removed the DataMisfitController() tests as these have started leading to test failures.

odunbar avatar Feb 05 '24 23:02 odunbar