Positive definite corrections in get_u_cov

[odunbar]

EnsembleKalmanProcesses.jl/src/EnsembleKalmanProcess.jl

Line 612 in 0db3339

Σ = use_prior_cov ? get_u_cov_prior(ekp) : get_u_cov_final(ekp)

I often find that the additive inflation breaks due to non-positive definiteness when estimating covariances, also the get_u_cov_final() may not break the subspace property.

We can resolve both of these issues by adding either

1. the function posdef_correct() around the get_u_cov_XXX(ekp) when estimating the covariances,
2. the function posdef_correct() within the definition of get_u_cov() and get_u_cov_prior themselves. (so this change would apply more broadly)

I prefer option 2, as i can't see why you wouldn't want positive definiteness. Though it would slow down the call to this function. However, get_u_cov is only called within the source code in the additive inflation step, and once per iteration for EKS & UKI. The function is also occasionally used by users to obtain diagnostics.

NB posdef_correct currently lives in the LearningRateSchedulers file for some reason

[odunbar]

Note - one other thing is that Julia Distributions is built on needing PD covariances, even though PSD is "OK"...
https://discourse.julialang.org/t/multivariate-normal-with-positive-semi-definite-covariance-matrix/3029
So another solution might be to just take the square-root of "C" outside (also OK for positive semi-definite matrices) and then multiply this by MvNormal(0,I)
Though I do recognise that these square-roots might have some numerical issues too.

[eviatarbach]

Is the problem with the additive inflation just that MvNormal expects a PD matrix but EKI can return PSD matrices (when the ensemble is smaller than the number of parameters)? I suspect that (besides the extra computational costs) enforcing PD may have unexpected consequences, since the matrices are actually being modified. I would favor only enforcing PD in the additive inflation function.

[odunbar]

OK no problem. Where I am at with this:

The actual bug is that we get P Semi-D matrices C, from (e.g. as you suggest, the empirical "current" ensemble) and this causes the MvNormal(0,C) to fail the EKP step.
we could
(i) fix this by doing sqrt(C) * MvNormal(0,I) - this doesn't break the subspace property, and doesn't cause failure
(ii) fix this by updating with a PD matrix Cbar (possibly in conjunction with (i) this does break the subspace property
(iii) make all calls to get_cov has a PD matrix [I have since realized this doesn't actually work, and I think you suggest that too above]