[SciPy-Dev] Allow specifying inverse-variance matrix to multivariate_normal (scipy.stats)

[Adrian Price-Whelan]

Hi all --

scipy.stats.multivariate_normal currently accepts a mean and
covariance matrix, and always computes the inverse-variance matrix
internally and explicitly. But I often have covariance matrices with
inverses that can be computed through other linear algebra tricks or
identities. It would therefore be great if
scipy.stats.multivariate_normal would allow optionally passing in an
inverse-variance matrix, which is sometimes useful for (a)
performance, or (b) numerical precision issues with covariance
matrices with large values. To be clear: I'm just suggesting that the
MVN distribution could accept an inverse-variance matrix to bypass the
inverse calculation, I'm not suggesting implementing or supporting
linear algebra identities.

To get around this in the past, e.g., to compute the log-PDF values at
some positions, I've typically implemented by own custom function, or
have used the implementation in astroML
(https://github.com/astroML/astroML/blob/master/astroML/utils/utils.py#L22).
It would be nice to have this supported in scipy directly!

Thoughts?

(Note: I made an issue for this idea, but was told to start discussion
here instead https://github.com/scipy/scipy/issues/11053)

best,
- adrian

-- 
Adrian M. Price-Whelan
Flatiron Institute, NYC
http://adrn.github.io