AbstractPDMat restriction causes problems

[oxinabox]

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]

Related: #940

[nickrobinson251]

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

[cscherrer]

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]

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]

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]

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.