`project` in `projectSymPosDef` projects in the tangent space

[Affie]

ApproxManifoldProducts.jl/src/services/KernelEval.jl

Line 7 in 0f4b8db

issymmetric(_c) ? _c : project(SymmetricPositiveDefinite(s[1]),_c,_c)

project a matrix from the embedding onto the tangent space T_p\mathcal P(n) of the [SymmetricPositiveDefinite] matrices, i.e. the set of symmetric matrices.
ref

[mateuszbaran]

Note that you can't really project onto the SPD manifold. For example, if _c is diagm([-1, 1, 1]), what would you expect to get? diagm([0, 1, 1]) is not a valid point on SPD manifold. If you want to admit positive semi-definite matrices, then it's a manifold with corners we don't have (yet).

[dehann]

Hi @mateuszbaran , we thought project SPD would be a good candidate for "fixing" covariance matrices that are nearly symmetric. Turns out Distributions.jl is very strict about symmetry, I mean too strict -- as in up maximum machine precision strict. I've had cases where attempted fix 0.5*(C + C') did not pass Distributions.jl requirments, although the matrix was clearly full rank (large eigen values etc.).

We are not looking for drastic projections such as diagm([-1,1]). It's a good example though. Perhaps if we limit the smallest allowable eigen value at say 1e-9, then accept a SPD projection as good enough? We could also consider limiting the allowable asymmetry.

[mateuszbaran]

Hi! I see the problem now. I've never encountered such strictness in Distributions.jl but it could happen due to use of Cholesky decomposition. Maybe wrapping in Symmetric would fix the issue?

By the way, I've recently been working with tangent space normal distributions (for Riemannian CMA-ES) and the algorithm actually samples from eigendecomposition instead of Cholesky decomposition (even the original state-of-the-art Euclidean one), so maybe that could be worth considering. If Symmetric is not enough to fix the matrix, it's likely eigendecomposition would have to be performed anyway.