LKJ distribution: check if matrix is positive definite.

[agoldin]

Without this check MCMC strays into neverland with negative eigenvalues for
covariance matrix. This is not good at all ( see
#873 ).

[springcoil]

Hey @agoldin this looks good to me - I think this is a good addition to PyMC3 especially since this kind of error has come up before - but is there any way you could add a quick test. Maybe deliberately trying a distribution with negative eigenvalues for this function and then proving that an exception is raised.

No disagreements about multivariate.py maybe @twiecki or @jsalvatier have some comments.

[agoldin]

I'll look where I can add tests over weekend. I am not yet very comfortable with github way of doing things. Thx.

[twiecki]

It's a bit surprising that T.all(eigh(X)[0] > 0) doesn't work. It raises a gradient error for you?

[agoldin]

The derivative of Eigh is defined via class EighGrad, which does not have gradient defined. Even if eigh appears only as an argument of switch() , guess_scaling wants to take second derivative of all the functions that logp() is composed of, even if those derivatives are later multiplied by 0. I'd imagine to fix it one should play with then (at this point it is too much for me). Instead I just defined a new operation and explicitly expressed that derivative of PositiveDefiniteMatrix() is zero. Not being able to use NUTS sampler is a problem, because it is so much better.

Anyway cholesky() should be faster than eigh(). It would be great to be able to use it on GPU, but even matrix inversion in Theano seems to be CPU only.

May be in far future, when pymc4 is based on tensor flow or something similar.

[twiecki]

@nouiz Do you think this makes sense from a theano perspective?

[agoldin]

@springcoil , pls check if this is a test you have in mind. Thanks.

[twiecki]

@agoldin This looks good to me.

[twiecki]

@agoldin can you replace the check in Wishart with this one?

[agoldin]

@twiecki , pls check.

[twiecki]

LGTM, thanks @agoldin!

[nouiz]

It seem good overall. But personnaly, I don't like a str that don't return the class name. I think it make things more confusing when you search. So personally, I would remove that method.

[twiecki]

Great, thanks for reviewing @nouiz!

[agoldin]

@nouiz , did not think of str. Will fix tonight.

[springcoil]

Sorry for missing this @agoldin - it looks good to me too. I'd fix the str though like @nouiz said.

[agoldin]

Thanks, @springcoil . I will not bother then to reduce confusion.