Fix to ensure consistency between insupport() and logpdf() for Inverse/Wishart

[brian-j-smith]

This pull request is meant to fix inconsistencies that can occur between insupport() results and logpdf() calculations for Wishart and InverseWishart distributions.

The issue:

• Let D denote a Wishart or InverseWishart distribution object, and X a real matrix to evaluate.
• logpdf(D, X) is intended to return -Inf if X is not in the support of D, and a numeric value otherwise.
• Support is tested with insupport(D, X).
• insupport() uses the Cholesky factorization of X, and logpdf() calculations involving X are non-Cholesky based.
• If X is close to being singular, Cholesky-based calculations of determinants, inverses, etc. can differ from non-Cholesky-based ones.
• This can lead to situations where insupport() returns true, but logpdf() calculations throw errors (e.g. the non-Cholesky-based inverse cannot be calculated, even though the Cholesky factorization can).
• In a nutshell, logpdf() can result in a terminal error instead of a numeric value or -Inf, as intended.

The solution proposed in this pull:

• Use Cholesky-based methods to compute the determinant and inverse of X in logpdf() - logdet(cholfact(X)) and inv(cholfact(X)) are the log-determinant and inverse of X, respectively.
• Note that trycholfact() is added and hasCholesky() modified a bit to allow the Cholesky factorization to be computed only once in a call to logpdf(). Determinants and inverses are then computed from that one factorization.
• Currently, a call to logpdf() can result in unrelated Cholesky factorization, determinant, and inverse calculations.
• Having them all be based on the factorization should make calculations consistent, faster, and more precise.

Feel free to let me know if there are any questions about or suggestions for the pull.

[coveralls]

Coverage decreased (-0.0%) when pulling 050ff40 on brian-j-smith:master into f4ed9f3 on JuliaStats:master.

[jiahao]

cc: @alanedelman

[brian-j-smith]

In re-reading my original post, I think it might sound more theoretical than intended…

The basic problem is that the insupport methods are based on Cholesky factorizations, whereas logpdf calculations are non-Cholesky-based. That inconsistency can lead to terminal errors (as detailed above). My pull gets rid of the inconsistencies.

In practice, I am running into the errors in my MCMC package (Mamba) when using the logpdf functions to sample covariance matrices that have Wishart or InverseWishart priors.

I have been testing the pull changes over the past week, and can confirm that my errors are gone and results are correct. The changes should be a net benefit to other users (fewer errors, faster and more precise calculations). I can’t think of any downsides, but would be happy to hear from anyone who has a different perspective. I also tried to be mindful of code readability and maintainability; e.g., only increased the total line count by one.

@jiahao or @lindahua: Do you have any specific questions about the pull or its merge-worthiness? Eliminating the errors is a priority for me, but I can talk about computation complexity advantages if you like.

Best.

[lindahua]

Looks good to me

[johnmyleswhite]

Yup, looks good to me as well