search

JuliaStats / Distributions.jl

A Julia package for probability distributions and associated functions.
Other
1.1k stars 413 forks source link

AbstractPDMat restriction causes problems #1219

Open oxinabox opened 3 years ago

oxinabox commented 3 years ago

So often the AbstractPDMat restriction on MvNormal and other Multivariate distributions causes pain. Discussed in slack with @timholy and @andreasnoack that i had a problem creating a Positive Definate Woodbury Matrix. Woodbury.jl wants it ot subtype Factorization, but for it to work with Distributions.jl it needs to subtype AbstractPDMat. A positive definate woodbury matrix is actually just a SymWoodbury with a few extra contraints about nonnegatives on diagonals. So ideally would actually just use SymWoodbury.

Other examples include if I wanted to work with a NamedDimsArray or KeyedArray, or some other AxisArray, and/or their respective factorizations. With very little effort one can get an indexed distribution as output -- if we could use indexed types within the internals.

I propose that the covariences be relaxed so some mininal set of operations that must be defined; and maybe to passing LinearAlgebra.isposdef The minimal set of operations would be something like inv. Then we can special case AbstractPDMat to use faster invquad. And we can special case Matrix inputs to be converted to PDMats. But other inputs can stay as they are. We can quite likely support any AbstractMatrix and even some Factorization subtypes, via ducktyping.

oxinabox commented 3 years ago

Related: https://github.com/JuliaStats/Distributions.jl/issues/940

nickrobinson251 commented 3 years ago

workaround for some cases: https://github.com/invenia/PDMatsExtras.jl

cscherrer commented 2 years ago

In my experience, there's essentially no benefit to carrying around the full covariance matrix; it works just fine to work in terms of the positive semidefinite Cholesky factor, or even better (for logpdf) the Cholesky factor of the precision.

Beyond the PDS requirement, PDMats has the unusual approach "allocate no matter what". Give it a Matrix, it will allocate for the Cholesky. Give it a Cholesky and it will allocate to store the full matrix. The latter is usually a waste, at least in a PPL context.

devmotion commented 2 years ago

Beyond the PDS requirement, PDMats has the unusual approach "allocate no matter what".

This pops up quite regularly, I think it would be good to add support for an AbstractPDMat type in PDMats that does not reconstruct the full matrix.

sethaxen commented 2 years ago

Agreed, most of the overloads that use .mat could be written to be more efficient or just a little less efficient by never computing the full matrix.

oxinabox commented 2 years ago

Another one came up, we have a MvGaussian that has a block diagonal structure. That is a pretty significant speedup for many operations; and a huge saving in memory.

BlockDiagonals.jl doesn't know anything about PDMats. I would rather not have to replicate it all in PDMatsExtra.jl just so i can fit it in to an AbstractPDMat.