search

CliMA / EnsembleKalmanProcesses.jl

Implements Optimization and approximate uncertainty quantification algorithms, Ensemble Kalman Inversion, and Ensemble Kalman Processes.
https://clima.github.io/EnsembleKalmanProcesses.jl/dev
Apache License 2.0
80 stars 18 forks source link

SparseInversion not as robust as Inversion #134

Open ilopezgp opened 2 years ago

ilopezgp commented 2 years ago

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
jinlong83 commented 2 years ago

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.

ilopezgp commented 2 years ago

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?

jinlong83 commented 2 years ago

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).

odunbar commented 2 years ago

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 commented 8 months ago

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