Enhancing Sample Quality through
Minimum Energy Importance Weights

Bayesian inference is a popular tool for solving many statistical and engineering problems. However, in general, the closed-form derivation of the posterior $\pi$ is typically unattainable in practice, and thus approximation by Monte Carlo methods are often employed. Markov chain Monte Carlo (MCMC) and Importance Sampling (IS) are two popular techniques, but each with its own weakness.

The key idea of MCMC (Metropolis et al., 1953; Hastings, 1970) is to carefully construct a Markov chain that has the target posterior $\pi$ as the equilibrium distribution. Though easy to implement in practice, its asymptotic convergence is difficult to analyze and no finite-sample guarantee is known. To offset the random initialization and the auto-correlation within the MCMC samples, in practice it is recommended to discard the burn-in samples followed by thinning to yield the final sample set for the posterior approximation (Robert and Casella, 2013). Hence, at least tens of thousands Markov chain iterations are required for obtaining sufficient samples, but this becomes a major computational bottleneck when the posterior $\pi$ is expensive to evaluate, e.g. the big data setting (Bardenet et al., 2017), the calibration of complex computer code (Joseph et al., 2019), etc., in which it would take days or even weeks to finish running the MCMC.

On the other hand IS (Robert and Casella, 2013) instead draw samples from a simple proposal distribution $q$ and correct for the distributional mismatch by assigning the importance weight that is proportional to $\pi(\cdot)/q(\cdot)$. Improved convergence can be achieved if one can draw Quasi-Monte Carlo samples (Niederreiter, 1992) from the proposal $q$ (Aistleitner and Dick, 2015). However, the performance of IS is sensitive to the choice of the proposal because a poor choice could possibly lead to an estimator with infinite variance (Owen, 2013). It is challenging to design a good proposal $q$ for the high-dimensional complex posterior given the limited class of proposals that are analytically tractable.

It is natural to consider combining the advantages of MCMC and IS: applying IS to correct for the distributional bias of the finite MCMC samples, which is believed to match more closely to the target posterior $\pi$ than the samples from the simple proposal typically used by IS. This idea is explored in Markov Chain Importance Sampling (Rudolf et al., 2020; Schuster and Klebanov, 2021), which also allows for recycling the burn-in and thinned samples. However, it cannot be extended for the adaptive variant of MCMC that enjoys faster mixing in practice (Andrieu and Thoms, 2008; Vihola, 2012). On the other hand, Liu and Lee (2017) proposed the Black-Box Importance Sampling (BBIS) that can calculate the importance weights for samples simulated from arbitrary unknown black-box mechanisms. The importance weights are computed by minimizing over the kernelized Stein discrepancy (Liu et al., 2016; Chwialkowski et al., 2016; Oates et al., 2017), a goodness-of-fit measure for testing whether the samples are generated from the target distribution. However, computing the kernelized Stein discrepancy requires knowing the score function, i.e., the gradient of the log-density function $\nabla_{x}\log\pi(x)$.

Unfortunately, the score function is not always available in many engineering applications, e.g., calibration of black-box computer codes (Ngo and Scordelis, 1967; Mak et al., 2018), where only input/output pairs are accessible. A much more common situation is where the score function exists but computing it is not straightforward, e.g., calibration for system of differential equations (Section 6). See Sengupta et al. (2014) for the existing numerical approaches on the gradient computation within the dynamical models. The adjoint method is suggested to be the computational and memory efficient approach, which involves solving another system of differential equations known as the sensitivity equation. Consequently, the expense of computing the gradient is at least on par with the cost of solving the system. This would certainly be a computational bottleneck for large-scale system, not to mention the potential numerical instability from solving the sensitivity equation.

To address this shortcoming of kernelized Stein discrepancy-based BBIS algorithm, we propose BBIS-MED, a novel Black-Box Importance Sampling method using minimum energy design (MED) (Joseph et al., 2015) that only requires the unnormalized (log) density function, making it applicable for almost all Bayesian inference problems. BBIS-MED computes the importance weights by minimizing over an energy criterion, which has demonstrated its success as a good criterion for finding representative samples of any distribution (Joseph et al., 2019). Moreover, extensive simulations are provided to demonstrate that our proposed BBIS-MED is orders of magnitude faster than the state-of-the-art BBIS algorithm, e.g., for a 55 dimensional logistic regression example in Subsection 5.3, BBIS-MED only takes 6 seconds to compute the importance weights for 5,000 samples, in contrast to the 1.3 hours required by the BBIS algorithm of Liu and Lee (2017).

The paper is organized as follows. Section 2 reviews the minimum energy designs (Joseph et al., 2015, 2019) that motivates the BBIS-MED algorithm. Section 3 proposes the computational efficient BBIS-MED algorithm that only requires the unnormalized density function. Section 4 outlines some existing BBIS algorithms from the literature, and extensive numerical comparisons to our proposed BBIS-MED is provided in Section 5. Section 6 presents a successful application of BBIS-MED on a real world Bayesian model calibration problem where the score function is not easy to compute. We conclude the article with some remarks in Section 7.

Let us start by defining the minimum energy designs (MED) in Joseph et al. (2015, 2019), which forms the foundation of our proposed Black-Box Importance Sampling algorithm.

MED is motivated from the optimal configuration of charged particles in physics. It considers the design points $\mathcal{D}_{n}$ as some positively charged particles and aims to place them in the proper position for minimizing the total potential energy. Under the charge function $\rho(x)=1/\pi^{1/(2p)}(x)$, Joseph et al. (2019) show that for $k\to\infty$, the MED converges asymptotically to the target distribution $\pi$ as $n\to\infty$.

MED is also known as the minimal Riesz energy points studied in the sphere-packing literature (Borodachov et al., 2008a, b, 2019) by recognizing that $d(x_{i},x_{j})^{-1}$ is the Riesz $1$-kernel. Again Borodachov et al. (2019) shows that using the charge function $\rho(x)=1/\pi^{1/(2p)}(x)$ yields the optimal configuration that has the target distribution $\pi$ as the limiting distribution, while only requiring $k>p$. Similar convergence result is also presented in Wang (2011). Given that the optimal configuration that minimizes (1) converges to $\pi$ as $n\to\infty$ when the charge function $\rho(x)=1/\pi^{1/(2p)}(x)$, in the limit we should have

where $\mu$ is any distribution supported on $\mathcal{X}$. The MED formulation (1) can be generalized to arbitrary symmetric proper kernel $K:\mathcal{X}\times\mathcal{X}\to\mathbb{R}\cup\{+\infty\}$, but less is known about the limiting distribution of the optimal configuration beyond the Riesz kernel (Borodachov et al., 2019).

Let $\{x_{i}\}_{i=1}^{n}$ be a set of samples in $\mathbb{R}^{p}$ that is simulated from an arbitrary (unknown) mechanism. The Black-Box Importance Sampling (BBIS) algorithm aims to find the normalized weights $\{w_{i}\}_{i=1}^{n}$, i.e., $w_{i}\geq 0$ and $\sum_{i=1}^{n}w_{i}=1$, such that the weighted samples $\{(x_{i},w_{i})\}_{i=1}^{n}$ give a good approximation to the target distribution $\pi$. Unlike Liu and Lee (2017) that use kernelized Stein discrepancy to measure the similarity between $\{(x_{i},w_{i})\}_{i=1}^{n}$ and $\pi$, we instead consider the energy criterion in (2) that require knowing $\pi$ only up to a constant of proportionality, making it applicable to almost all Bayesian problems.

More specifically, the Monte Carlo approximation for the expectation in (2), the continuous energy criterion, using weighted samples $\{(x_{i},w_{i})\}_{i=1}^{n}$ is

Compared to the MED formulation in (1), the approximation using weighted samples in (3) includes the self-interaction terms, i.e., terms with $i=j$, in the summation. Without these self-interaction terms, the minimum value of (3) is always obtained at 0 with weights $w_{i}=1$ for some $i$ and $w_{j}=0$ for all $j\neq i$. This is expected since if there is only one particle in the system, there will be no interaction (repulsion), and the total potential energy will always be minimized at 0. Hence, the self-interaction terms are required when weighted samples are used for approximating the energy criterion. Moreover, the self-interaction terms can also serve as a penalization to avoid the weights collapsing to only a few samples.

However, after including the self-interaction terms, (3) always yields infinity for any normalized weights $\{w_{i}\}_{i=1}^{n}$ since $d(x_{i},x_{j})=0$ for $i=j$. One fix is to add some small positive term $\delta$ to the distance computation to bound it away from 0, i.e, we consider the following formulation,

where the Riesz kernel in (3) is replaced by the inverse multiquadric kernel $K(x_{i},x_{j})=(d^{2}(x_{i},x_{j})+\delta)^{-1/2}$ for some small $\delta>0$. Unfortunately, there is no known result on the limiting distribution for the optimal configuration of MED for the inverse multiquadric kernel. However, given that the inverse multiquadric kernel with small $\delta$ closely mimics the Riesz kernel except when $d(x_{i},x_{j})\approx 0$, one can assume that the difference between the limiting distribution of the optimal configuration with respect to the inverse multiquadric kernel and the limiting distribution with respect to the Riesz kernel should be very small. Hence, utilizing the limiting behavior of MED with respect to the Riesz kernel discussed in Section 2, one can believe that for $\rho(x)=1/\pi^{1/(2p)}(x)$, the minimal value of (4) with respect to the weights $\{w_{i}\}_{i=1}^{n}$ is obtained when the weighted samples $\{(x_{i},w_{i})\}_{i=1}^{n}$ are approximately following $\pi$. Recall that $\{w_{i}\}_{i=1}^{n}$ serve as the self-normalized importance weights that correct for the mismatch between the target distribution $\pi$ and the unknown underlying distribution of $\{x_{i}\}_{i=1}^{n}$.

This leads to the following optimization problem in matrix form for computing the black-box importance weights $w^{*}=(w^{*}_{1},\ldots,w^{*}_{n})$ that matches any samples $\{x_{i}\}_{i=1}^{n}$ to the target distribution $\pi$ via minimizing the energy criterion,

where $\Delta^{n}$ is the $n$-dimensional probability simplex and

where $\gamma(x)=\pi(x)/Z$ is the unnormalized density function and $C=Z^{k/p}$ is some constant depending on the unknown normalizing constant $Z$. Since the constant $C$ can be factored out without changing the optimum solution for (5), our proposed BBIS algorithm only requires the unnormalized log-density function, which is the most elementary prerequisite for any sampling method. Moreover, the optimization problem in (5) is a convex quadratic programming problem (see Appendix B.1) that can be solved efficiently with theoretical convergence guarantee by the simplex constrained convex optimization techniques, including mirror descent (Nemirovskij and Yudin, 1983; Beck and Teboulle, 2003), projected gradient descent (Duchi et al., 2008), and Frank-Wolfe (Jaggi, 2013). We employ projected gradient descent for the simulations presented in this paper.

Last, let us discuss the choices for the parameters in (5).

• •

  For the distance function $d(x_{i},x_{j})$, we use the Mahalanobis distance,

  $\displaystyle d^{2}(x_{i},x_{j})=(x_{i}-x_{j})^{T}\Sigma^{-1}(x_{i}-x_{j}),$

  where $\Sigma$ is the covariance matrix estimated from the samples $\{x_{i}\}_{i=1}^{n}$. It has been shown in (Joseph et al., 2019) that for computing MED of correlated distributions, better empirical performance is observed from Mahalanobis distance compared to Euclidean distance.

• •

  For $\delta$, we find that after standardization of the samples $\{x_{i}\}_{i=1}^{n}$ to be component-wise zero mean and unit variance, using $\delta=0.01$ yields robust empirical performance on problems across different dimensions.

• •

  For $k$, the theoretical result of MED from Section 2 suggests that one should use $k>p$, but empirically the performance is not good (see Appendix B.2). Our guess is that the theoretical results no longer hold when the self-interaction terms are included in the energy criterion. Using $k=1$ gives the most robust performance from extensive numerical simulations.

Given that our proposed BBIS algorithm is motivated from MED, for the rest of paper we will refer it as BBIS-MED and the resulting weights $w^{*}$ as the minimum energy importance weights. Though the theoretical guarantee of BBIS-MED requires future study, which is known to be a very difficult mathematical problem, extensive simulations are provided in the later Sections to demonstrate that any arbitrary simulated samples $\{x_{i}\}_{i=1}^{n}$ with the minimum energy importance weights $w^{*}$ provide better approximation for the target distribution $\pi$ than the unweighted samples $\{x_{i}\}_{i=1}^{n}$ alone.

Let us now review some existing BBIS algorithms in the literature that we will compare BBIS-MED to empirically in Section 5.

In this section we provide extensive simulations comparing our proposed BBIS-MED to the two other BBIS algorithms discussed in Section 4: BBIS-KSD and BBIS-KDE.

Now we demonstrate the application of our proposed BBIS-MED algorithm on a real world model calibration problem where the score function of the posterior is difficult to compute. Vapor phase infiltration (VPI) process is one important class of chemical processes that can be used for transforming organic polymer membranes into hybrid organic-inorganic membranes (Leng and Losego, 2017). These hybrid membranes enjoys 10 times energy efficiency than distillation for the chemical separation process, e.g. purification of clean water. Understanding the chemical behavior of the VPI process is critical for designing those membranes. Ren et al. (2021) proposed a reaction-diffusion partial differential equations (PDEs) model $g(t;\theta)$ for studying VPI, where $t$ denotes the time in seconds and $\theta$ denotes the model parameters. See Appendix C for the model details. Figure 7 shows the evaluation of the PDEs model $g$ on some randomly chosen $\theta$’s, and none is aligned well with the experimental observation. Given that only temporal observations are available from experiment, i.e., spatial information is aggregated, this limits the choice of methods we can use to calibrate $\theta$.

One approach is to apply the Bayesian calibration (Kennedy and O’Hagan, 2001). For simplicity assume there is no model discrepancy, thus the measurement errors are independent, and we can model them using i.i.d. Gaussian with mean zero and variance $\sigma^{2}$, i.e., the likelihood for $y(t)$, the noisy observation at time $t$, is

For the prior, we consider $p(\sigma^{2})\propto 1/\sigma^{2}$ and $p(\theta)$ for the PDEs parameter (see Appendix C). It follows that the posterior is

where $\mathcal{T}=\{t_{1},\ldots,t_{N}\}$ are the 251 equally spaced out time points on $[0,125000]$ that we have both experimental observation and the PDEs model output. After integrating out $\sigma^{2}$ we have

Note that the PDEs model $g$ is in the log posterior, and hence the score function of $p(\theta|y)$ cannot be easily computed¹¹1Though for this particular PDEs example the score function can be computed via the adjoint method, it would still be computationally expensive. Moreover the purpose here is to demonstrate BBIS-MED can be easily applied without knowing the score function, which is common in many real world model calibration problems that the computer code is completely black-box.. Moreover, PDEs is expensive to solve, and each evaluation of the log posterior takes about 7 seconds on a 2 cores 2.3 GHz laptop using the deSolve package in R (Soetaert et al., 2010). Given the limited computational resource, we can only afford to run the adaptive MCMC (Vihola, 2012) for 2,500 iterations, but this is far from converging by comparing the densities of the post burn-in MCMC samples (first half of the chain is discarded for burn-in) to the histograms of the “groundtruth” posterior samples in Figure 8. The “groundtruth” posterior samples are obtained by running adaptive MCMC for 25,000 iterations with the chain initialized at the mean of the post burn-in MCMC samples and the initial kernel covariance is set to be the diagonal of the covariance of the post burn-in MCMC samples. Again the first 5,000 samples are discarded for burn-in and the rest 20,000 samples are used for the “groundtruth” posterior samples.

Fortunately, our proposed BBIS-MED algorithm could be the remedy. Applying BBIS-MED on the 2,500 MCMC proposal samples we obtain the weighted samples that provide a much better approximation for the posterior (red dashed lines versus the histogram’s in Figure 8). The improvement over the post burn-in MCMC samples can also be quantified by the energy distance (0.044 versus 0.067). Lastly, Figure 7 shows the PDEs output evaluated at the mean of the BBIS-MED weighted samples in red solid line, and it better matches to the experimental observations than the randomly chosen $\theta$’s.

In this paper we propose a novel Black-Box Importance Sampling algorithm, BBIS-MED, that is built upon the energy criterion (Joseph et al., 2015, 2019) motivated from the optimal configuration of charged particles in physics. Different from the kernelized Stein discrepancy based BBIS algorithm (Liu and Lee, 2017), our proposed BBIS-MED only require the unnormalized (log) density function, making it applicable to almost all Bayesian problems. Moreover, BBIS-MED also outperforms other BBIS algorithms that do not need the score function, e.g. the KDE-based BBIS algorithm (Delyon and Portier, 2016), and the unweighted baseline on extensive simulations. The order of magnitude reduction in computational time is yet another great benefit of our proposed BBIS-MED algorithm.

Though the theoretical guarantee of BBIS-MED requires further study, we observe its outstanding empirical performances on a wide range of problems, including strongly correlated distributions (Figure 1(b)), high dimensional multi-modal distributions (Figure 5), and the complicated posterior from real world applications (Figure 6 and Figure 8). Moreover, the comparison of BBIS-MED to the Monte Carlo baseline (Figure 3) also further strengthen the validity of our proposed method.

Although this paper focus on applying BBIS-MED to correct the bias from the finite MCMC samples, i.e., a simple yet efficient post-processing step to improve the approximation quality of any MCMC samples, BBIS-MED has much broader potential applications, including optimal MCMC thinning (Chopin and Ducrocq, 2021; Riabiz et al., 2022), result refinement of complex variational proposals (Wainwright et al., 2008; Rezende and Mohamed, 2015), covariance shift in model training (Sugiyama et al., 2007), extension to black-box importance weights for discrete variables (Appendix D), and many other statistical and machine learning problems.