Skip to content

Conversation

@mdhaber
Copy link
Contributor

@mdhaber mdhaber commented Oct 1, 2022

Reference issue

Followup to gh-16002

What does this implement/fix?

gh-16002 added the ability to specify a covariance matrix via its inverse.
This PR adds the ability to specify a diagonal covariance matrix via its Cholesky decomposition.

@mdhaber mdhaber added scipy.stats enhancement A new feature or improvement labels Oct 1, 2022
Copy link
Member

@tupui tupui left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM, thanks Matt 👍

Do we plan to add more decomposition options?

@tupui tupui merged commit 9ae8fd0 into scipy:main Oct 2, 2022
@tupui tupui added this to the 1.10.0 milestone Oct 2, 2022
@mdhaber
Copy link
Contributor Author

mdhaber commented Oct 2, 2022

Do we plan to add more decomposition options?

I have two more on deck (that I've already written): eigendecomposition and LDL. LDL is the one that should be the default if we want to be able to instantiate a Covariance object directly because it is faster than eigendecomposition and more robust than Cholesky. Eigendecomposition is important because it's probably the most accurate. It would be easy to add others - basically, any matrix that has special structure for which we have a dedicated solver (e.g. solveh_banded, solve_toeplitz) could be useful. It may also be useful to add a sparse option, and we could add something for rotation matrices to satisfy gh-15675. But I can let others have fun with these; I'm just paving the way and adding the ones I think are most fundamental.

After that, @siddhantwahal has wanted (for a long time, see gh-11772) for rvs to re-use the decomposition. So if the user provides a Covariance object, we'll use it (instead of np.multivariate_normal) to generate random variates. You just take independent, standard normal random variate and color them (inverse of whiten).

Also, with a few modifications, these classes/multivariate_normal can broadcast over ND stacks of covariance matrices (and ND mean). This is something that @OriolAbril wanted to see for ArviZ/PyMC. I already made the required changes and wrote tests in my own branch; it's just a matter of putting them in a reasonably-sized PR. (All of these have essentially been copy-pasted from mdhaber#88, which is why they are coming one after the other.)

@tupui
Copy link
Member

tupui commented Oct 3, 2022

That's all great. I also just realizing that this work can also help qmc.MultivariateNormalQMC.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement A new feature or improvement scipy.stats

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants