Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Multivariate Gaussian Problem with Assumptions on Covariance #791

Open
janrpeters opened this issue Nov 2, 2018 · 16 comments
Open

Multivariate Gaussian Problem with Assumptions on Covariance #791

janrpeters opened this issue Nov 2, 2018 · 16 comments

Comments

@janrpeters
Copy link
Contributor

I am working on a few tutorials in machine learning as well as papers. One recurrent theme in these algorithms is that there are linear transformations applied to Gaussian variables, e.g.,

δu = rand(MvNormal([0,0],𝐖));
δx = 𝐉*δu;

However, in most cases, we aim to transform the covariance matrix and use

δx = rand(MvNormal([0,0],𝐉'𝐖𝐉));

instead (that happens in really well-known things such as Kalman filters). It also causes a problem with Julia's implementation of Gaussians as ,𝐉'𝐖𝐉 is not always positive definite and with

𝐉 = [-0.707805 0.706408; 0.706408 0.707805];
𝐖 = diagm(0=>[1,0.1])/20;

our Method 1 would work but Method 2 would fail as 𝐉'𝐖𝐉 is only semi-positive definite. This case of high practical relevance! May I appeal to changing the requirement from positive definite to semi-positive definite?

I would like to use Distributions.jl directly instead of having to create a wrapper around it!

@simonbyrne
Copy link
Member

𝐉'_𝐖_𝐉 is only semi-positive definite

Actually, it's because it's not quite symmetric (if you wrap it in a Symmetric block, it is perfectly fine).

If we proceed with #307, we could allow something like

dx = MvNormal([0,0],𝐖)
du = 𝐉' * dx

Though what that would do internally isn't quite clear (i.e. how should the covariance component of du be stored?)

@simonbyrne
Copy link
Member

This is actually a good argument for getting rid of PDMats, so as to make it easier for custom factorizations (see #688 (comment)).

If all you're doing is sampling, then the most efficient way to store it is as a 𝐉 and the Cholesky factor of 𝐖 (which would have already been precomputed for dx).

If you want to compute pdfs, then you would need to do some more factorization of 𝐉 (probably an RQ decomposition, since then you could multiply the upper triangular bits together to another upper triangular part).

@janrpeters
Copy link
Contributor Author

How is 𝐉'𝐖𝐉, that is,

2×2 Array{Float64,2}:
0.0275445 -0.0225
-0.0225 0.0274556

not symmetric?

@simonbyrne
Copy link
Member

Floating point rounding error:

julia> A = 𝐉'* 𝐖 * 𝐉;

julia> A-A'
2×2 Array{Float64,2}:
 0.0          -3.46945e-18
 3.46945e-18   0.0

@janrpeters
Copy link
Contributor Author

Oh wow! Thanks!

Could we maybe always wrap this around? At least for errors of rounding magnitude? I am pretty sure that I am not the only one who ended up one afternoon cursing and wondering whether something was wrong in my math as it happened only for 1/230 cases...

@janrpeters
Copy link
Contributor Author

Allover, I am so in love with Julia & Distributions.jl. But getting such kinks out would be really good, I want to make my Julia code so nice that I can avoid doing pseudo-code and doing without a Symmetric wrapper would be really nice in this respect.

@simonbyrne
Copy link
Member

Could we maybe always wrap this around? At least for errors of rounding magnitude? I am pretty sure that I am not the only one who ended up one afternoon cursing and wondering whether something was wrong in my math as it happened only for 1/230 cases...

Unfortunately it's far from trivial to figure out what an appropriate tolerance should be. We should at least give a more informative error though (i.e. suggest the Symmetric wrapper).

@matbesancon
Copy link
Member

@simonbyrne is there something to do on this one? Add a @warn somewhere?

@simonbyrne
Copy link
Member

The warning would probably go in the PDMats constructor.

@SebastianCallh
Copy link

I am using Optim together with Distributions for Laplace approximations and run into the same issues where the estimated covariance matrix very often have rounding errors and is not exactly symmetric. This causes the call to MvNormal to crash the program with PosDefException: matrix is not Hermitian; Cholesky factorization failed.
What is the best way of dealing with this situation at the moment? I cannot find any mention of the Symmetric wrapper in the documentation so I guess things have changed since this issues was created.

@devmotion
Copy link
Member

I found https://github.com/timholy/PositiveFactorizations.jl to be helpful for dealing with numerical issues arising from estimated covariance matrices. You should be able to compute the (approximate) Cholesky decomposition with PositiveFactorizations manually, and then construct a PDMat that you can use a covariance matrix for MvNormal.

@btmit
Copy link

btmit commented Apr 9, 2021

Symmetry issues aside, the fact remains that MvNormal fails with positive semi-definite matrices. Should this be considered a bug?

@devmotion
Copy link
Member

PDMatsExtras adds support for psd matrices to PDMats. You should be able to use it to wrap your psd covariance matrix in a format that MvNormal accepts.

@btmit
Copy link

btmit commented Apr 21, 2021

@devmotion: quick question on usage of PDMatsExtras. Here's a short example of what I'm running into:

using PDMats, PDMatsExtras

m = rand(2,2);  # some 2x2 matrix
pd = PDMat(a*a');  # some PD matrix

e = [0, 1];
psd = PSDMat(Symmetric(a*Diagonal(e)*a'));  # some PSD matrix
eigmin(psd)  # check PSD

X_A_Xt(pd, a)  # all good
X_A_Xt(psd, a)  #! not defined

Should I be importing both PDMats and PDMatsExtras like I am or is there a different intended usage to get the extra methods on these sorts of functions?

@devmotion
Copy link
Member

devmotion commented Apr 21, 2021

I am not familiar with PDMatsExtras (have neither used it nor contributed to it) but to me it seems this is a bug in PDMatsExtras and the definitions in https://github.com/invenia/PDMatsExtras.jl/blob/2c689308b846d69cfbced00d4108f77105a60de3/src/psd_mat.jl#L109-L132 should be qualified (i.e., they should define PDMats.X_A_Xt etc.). Currently they do not implement the PDMats interface but rather define different internal methods such as PDMatsExtras.X_A_Xt. It would be good to report this issue in PDMatsExtras.

@btmit
Copy link

btmit commented Apr 21, 2021

Thanks and sorry for the making the assumption you were a contributor.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

6 participants