Uncertainty estimation in high-dimensional inference problems without sampling?

Question

I'm working on a high-dimensional inference problem (around 2000 model parameters) for which we are able to robustly perform MAP estimation by finding the global maximum of the log-posterior using a combination of gradient-based optimisation and a genetic algorithm.

I'd very much like to be able to make some estimate of the uncertainties on the model parameters in addition to finding the MAP estimate.

We are able to efficiently calculate the gradient of the log-posterior with respect to the parameters, so long-term we're aiming to use Hamiltonian MCMC to do some sampling, but for now I'm interested in non-sampling based estimates.

The only approach I know of is to calculate the inverse of the Hessian at the mode to approximate the posterior as multivariate normal, but even this seems infeasible for such a large system, since even if we calculate the $\sim 4\times10^{6}$ elements of the Hessian I'm sure we couldn't find its inverse.

Can anyone suggest what kind of approaches are typically used in cases like this?

Thanks!

EDIT - additional information about the problem

Background
This is an inverse-problem related to a large physics experiment. We have a 2D triangular mesh which describes some physical fields, and our model parameters are the physical values of those fields at each vertex of the mesh. The mesh has about 650 vertices, and we model 3 fields, so that's where our 2000 model parameters come from.

Our experimental data is from instruments which do not measure these fields directly, but quantities that are complicated non-linear functions of the fields. For each of the different instruments we have a forward-model which maps the model parameters to predictions of the experimental data, and a comparison between the prediction and the measurement yields a log-likelihood.

We then sum up the log-likelihoods from all these different instruments, and also add some log-prior values which apply some physical constraints to the fields.

Consequently I doubt this 'model' falls neatly into a category - we don't have a choice of what the model is, it is dictated by how the actual instruments function that collect our experimental data.

Data set
The data set is composed of 500x500 images, and there is one image for each camera so total data points is 500x500x4 = $10^6$.

Error model
We take all errors in the problem to be Gaussian at the moment. At some point I might try to move over to a student-t error model just for some extra flexibility, but things still seem to function well with just Gaussians.

Likelihood example
This is a plasma physics experiment, and the vast majority of our data comes from cameras pointed at the plasma with particular filters in front of the lenses to look only at specific parts of the light spectrum.

To reproduce the data there are two steps; first we have to model the light that comes from the plasma on the mesh, then we have to model that light back to a camera image.

Modelling the light that comes from the plasma unfortunately depends on what are effectively rate coefficients, which say how much light is emitted by different processes given the fields. These rates are predicted by some expensive numerical models, so we have to store their output on grids, and then interpolate to look up values. The rate function data is only ever computed once - we store it then build a spline from it when the code starts up, and then that spline gets used for all the function evaluations.

Suppose $R_1$ and $R_2$ are the rate functions (which we evaluate by interpolation), then the emission at the $i$'th vertex of the mesh $\mathcal{E}_i$ is given by $$ \mathcal{E}_i = R_1(x_i, y_i) + z_i R_2(x_i, y_i) $$ where $(x,y,z)$ are the 3 fields we model on the mesh. Getting the vector of emissions to a camera image is easy, it's just multiplication with a matrix $\mathbf{G}$ which encodes what parts of the mesh each camera pixel looks through.

Since the errors are Gaussian the log-likelihood for this particular camera is then $$ \mathcal{L} = -\frac{1}{2} (\mathbf{G}\vec{\mathcal{E}} - \vec{d})^{\top}\mathbf{\Sigma}^{-1} (\mathbf{G}\vec{\mathcal{E}} - \vec{d}) $$

where $\vec{d}$ is the camera data. The total log-likelihood is a sum of 4 of the above expressions but for different cameras, which all have different versions of the rate functions $R_1, R_2$ because they're looking at different parts of the light spectrum.

Prior example
We have various priors which effectively just set certain upper and lower bounds on various quantities, but these tend not to act too strongly on the problem. We do have one prior that acts strongly, which effectively applies Laplacian-type smoothing to the fields. It also takes a Gaussian form: $$ \text{log-prior} = -\frac{1}{2}\vec{x}^{\top}\mathbf{S}\vec{x} -\frac{1}{2}\vec{y}^{\top}\mathbf{S}\vec{y} -\frac{1}{2}\vec{z}^{\top}\mathbf{S}\vec{z} $$

Answer 1

First of all, I think your statistical model is wrong. I change your notation to one more familiar to statisticians, thus let

$$\mathbf{d}=\mathbf{y}=(y_1,\dots,y_N),\ N=10^6$$

be your vector of observations (data), and

$$\begin{align} \mathbf{x}&=\boldsymbol{\theta}=(\theta_1,\dots,\theta_p) \\ \mathbf{y}&=\boldsymbol{\phi}=(\phi_1,\dots,\phi_p) \\ \mathbf{z}&=\boldsymbol{\rho}=(\rho_1,\dots,\rho_p), \ p \approx 650 \\ \end{align}$$

your vectors of parameters, of total dimension $d=3p \approx 2000$. Then, if I understood correctly, you assume a model

$$ \mathbf{y} = \mathbf{G}\mathbf{r_1}(\boldsymbol{\theta}, \boldsymbol{\phi})+\boldsymbol{\rho}\mathbf{G}\mathbf{r_2}(\boldsymbol{\theta}, \boldsymbol{\phi}))+\boldsymbol{\epsilon},\ \boldsymbol{\epsilon}\sim\mathcal{N}(0,I_N) $$

where $\mathbf{G}$ is the $N\times d$ spline interpolation matrix.

This is clearly wrong. There's no way the errors at different points in the image from the same camera, and at the same point in images from different cameras, are independent. You should look into spatial statistics and models such as generalized least squares, semivariogram estimation, kriging, Gaussian Processes, etc.

---

Having said that, since your question is not whether the model is a good approximation of the actual data generating process, but how to estimate such a model, I'll show you a few options to do that.

HMC

2000 parameters is not a very large model, unless you're training this thing on a laptop. The dataset is bigger ($10^6$ data points), but still, if you have access to cloud instances or machines with GPUs, frameworks such as Pyro or Tensorflow Probability will make short work of such a problem. Thus, you could simply use GPU-powered Hamiltonian Monte Carlo.

Pros: "exact" inference, in the limit of a infinite number of samples from the chain.

Cons: no tight bound on the estimation error, multiple convergence diagnostic metrics exist, but none is ideal.

Large sample approximation

With an abuse of notation, let's denote by $\theta$ the vector obtained by concatenating your three vectors of parameters. Then, using the Bayesian central limit theorem (Bernstein-von Mises), you could approximate $p(\theta\vert \mathbf{y})$ with $\mathcal{N}(\hat{\theta_0}_n,I_n^{-1}(\theta_0))$, where $\theta_0$ is the "true" parameter value, $\hat{\theta_0}_n$ is the MLE estimate of $\theta_0$ and $I_n^{-1}(\theta_0)$ is the Fisher information matrix evaluated at $ \theta_0$. Of course, $\theta_0$ being unknown, we'll use $I_n^{-1}(\hat{\theta_0}_n)$ instead. The validity of the Bernstein-von Mises theorem depends on a few hypotheses which you can find, e.e g., here: in your case, assuming that $R_1,R_2$ are smooth and differentiable, the theorem is valid, because the support of a Gaussian prior is the whole parameter space. Or, better, it would be valid, if your data were actually i.i.d. as you assume, but I don't believe they are, as I explained in the beginning.

Pros: especially useful in the $p<<N$ case. Guaranteed to converge to the right answer, in the iid setting, when the likelihood is smooth and differentiable and the prior is nonzero in a neighborhood of $\theta_0$.

Cons: The biggest con, as you noted, is the need to invert the Fisher information matrix. Also, I wouldn't know how to judge the accuracy of the approximation empirically, short of using a MCMC sampler to draw samples from $p(\theta\vert \mathbf{y})$. Of course, this would defeat the utility of using B-vM in the first place.

Variational inference

In this case, rather than finding the exact $p(\theta\vert \mathbf{y})$ (which would require the computation of a $d-$dimensional integral), we choose to approximate $p$ with $q_{\phi}(\theta)$, where $q$ belongs to the parametric family $\mathcal{Q}_{\phi}$ indexed by the parameter vector $\phi$. We look for $\phi^*$ s.t. some measure of discrepancy between $q$ and $p$ is minimzed. Choosing this measure to be the KL divergence, we obtain the Variational Inference method:

$$\DeclareMathOperator*{\argmin}{arg\,min} \phi^*=\argmin_{\phi\in\Phi}D_{KL}(q_{\phi}(\theta)||p(\theta\vert\mathbf{y}))$$

Requirements on $q_{\phi}(\theta)$:

• it should be differentiable with respect to $\phi$, so that we can apply methods for large scale optimization, such as Stochastic Gradient Descent, to solve the minimization problem.
• it should be flexible enough that it can approximate accurately $p(\theta\vert\mathbf{y})$ for some value of $\phi$, but also simple enough that it's easy to sample from. This is because estimating the KL divergence (our optimization objective) requires estimating an expectation w.r.t $q$.

You might choose $q_{\phi}(\theta)$ to be fully factorized, i.e., the product of $d$ univariate probability distributions:

$$ q_{\phi}(\theta)=\prod_{i=1}^d q_{\phi_i}(\theta_i)$$

this is the so-called mean-field Variational Bayes method. One can prove (see, e.g., Chapter 10 of this book) that the optimal solution for each of the factors $q_{\phi_j}(\theta_j)$ is

$$ \log{q_j^*(\theta_j)} = \mathbb{E}_{i\neq j}[\log{p(\mathbf{y},\theta)}] + \text{const.}$$

where $p(\mathbf{y},\theta)$ is the joint distribution of parameters and data (in your case, it's the product of your Gaussian likelihood and the Gaussian priors over the parameters) and the expectation is with respect to the other variational univariate distributions $q_1^*(\theta_1),\dots,q_{j-1}^*(\theta_{j-1}),q_{j+1}^*(\theta_{j+1}),\dots,q_{d}^*(\theta_{d})$. Of course, since the solution for one of the factors depends on all the other factors, we must apply an iterative procedure, initializing all the distributions $q_{i}(\theta_{i})$ to some initial guess and then iteratively updating them one at a time with the equation above. Note that instead of computing the expectation above as a $(d-1)-$dimensional integral, which would be prohibitive in your case where the priors and the likelihood aren't conjugate, you could use Monte Carlo estimation to approximate the expectation.

The mean-field Variational Bayes algorithm is not the only possible VI algorithm you could use: the Variational Autoencoder presented in Kingma & Welling, 2014, "Auto-encoding Variational Bayes" is an interesting alternative, where, rather than assuming a fully factorized form for $q$, and then deriving a closed-form expression for the $q_i$, $q$ is assumed to be multivariate Gaussian, but with possibly different parameters at each of the $N$ data points. To amortize the cost of inference, a neural network is used to map the input space to the variational parameters space. See the paper for a detailed description of the algorithm: VAE implementations are again available in all the major Deep Learning frameworks.

Answer 2

you may want to check out some of the "bayesX" software and possibly also the "inla" software. both of these are likely to have some ideas that you can try. Google it

both rely very heavily on exploiting sparsity in the parameterisation of the precision matrix (I.e conditional independence, markov type model) - and have inversion algorithms designed for this. most of the examples are based on either multi level or auto regressive guassian models. should be fairly similar to the example you posted