Use cov(u,dims=2) not cov(u,u,dims=2)

[odunbar]

EnsembleKalmanProcesses.jl/src/EnsembleKalmanProcess.jl

Line 238 in 3ba08bb

return cov(u, u, dims = 2)

In some examples (e.g. with additive inflation) I have encountered non-positive empirical covariance matrices from calls such as the above.

Even though this should be equal to cov(u,dims=2) in some instances of u the call to cov(u,u,dims=2) produces order(machine precision) asymmetry, and thus isposdef(cov(u,u,dims=2))=false. meanwhile, cov(u,dims=2) seems to always return symmetric matrices and is a safer option as isposdef(cov(u,dims=2))=true.

e.g.

julia> y = rand( MvNormal(I(5)),10);

julia> cyy = cov(y,y,dims=2);

julia> cy = cov(y,dims=2);

julia> cyy - permutedims(cyy,(2,1))
5×5 Matrix{Float64}:
 0.0           0.0          0.0           0.0          -1.38778e-17
 0.0           0.0          0.0           0.0           2.77556e-17
 0.0           0.0          0.0           0.0          -5.55112e-17
 0.0           0.0          0.0           0.0           2.77556e-17
 1.38778e-17  -2.77556e-17  5.55112e-17  -2.77556e-17   0.0

julia> cy - permutedims(cy,(2,1))
5×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

And so the importantly, cyy fails the checks:

julia> cy ≈ cyy
true

julia> isposdef(cy)
true

julia> isposdef(cyy)
false