I am wondering if it is possible to specify a guide that is a high-dimensional MVG with a sparse precision matrix. I am dealing with a high-dimensional regression problem in which the number of parameters is ~ 1 million and so there is no shot at modeling a full covariance matrix. However, I have a strong prior on which parameters are correlated with each other.
Due to the enormous size of the precision / covariance matrices, Pyro would need to have some way of representing the sparse matrix efficiently.
specifying the precision matrix would probably be the least of your problems, since you would also need to be able to sample from the distribution (which involves the inverse square root of the precision) and compute the log_prob (which would involve e.g. computing the log determinant of the precision). all of this is possible but it wouldn’t be a small lift. your best bet is probably to use AutoLowRankMultivariateNormal
Hi Martin, I came across your earlier response about using sparse precision matrices for high-dimensional Gaussian guides in Pyro.
I recently released a PyTorch package implementing selected inversion for sparse block matrices, with autograd, sparse solves, Gaussian sampling, log-determinants, and selected marginal covariances.
I’m aware that Pyro has its own structured Gaussian tools, including AutoGaussian, and that my implementation is not a scalable replacement for those.
I’m curious whether a standalone differentiable selected-inverse primitive would ever be useful in this setting, or whether Pyro’s existing machinery already covers the relevant use cases.
it could be but presumably you’d either need a concise way to specify the precision matrix or you’d need some algorithm that does so based on the structure of the model. and that’s where most of the hard work would be?
My library operates one layer below that: it assumes the block-sparse precision matrix and its sparsity pattern are already specified, then provides differentiable selected covariances, solves, log-determinants, and sampling. It doesn’t automatically infer those structures from a Pyro model.
I can see how extracting and parameterizing the precision structure from an arbitrary model, as AutoGaussian does, would be a big part of a general-purpose integration.
My intended use case is narrower: models where the sparse precision structure is already known, such as GMRFs or graph-structured Gaussian variational posteriors.
Do you think there’s room for a reusable numerical backend in that setting, or do Pyro’s existing Gaussian backends already handle those cases adequately?
haven’t thought about this much but it seems like the simplest start place would be to consider models with, say, random variables a,b,c,d,e,f, … and then allow specifying options like “fully diagonal but arbitrary c-f and b-d factors” in precision matrix or “fully diagonal but arbitrary c-f and b-d factors as well as full-rank a-a factor” etc etc. one question is how often will such ansaetze be sufficiently sparse for any of this to matter?
numpyro allows things along these lines for mass matrices:
*+ dense_mass=\[("x", "y")\]: use a dense mass matrix for the joint*
*(x, y) and a diagonal mass matrix for z*
*+ dense_mass=\[\] (equivalent to dense_mass=False): use a diagonal mass*
*matrix for the joint (x, y, z)*
*+ dense_mass=\[("x", "y", "z")\] (equivalent to full_mass=True):*
*use a dense mass matrix for the joint (x, y, z)*
*+ dense_mass=\[("x",), ("y",), ("z")\]: use dense mass matrices for*
*each of x, y, and z (i.e. block-diagonal with 3 blocks)*