Keywords: Gibbs measures, Large deviations, Monte Carlo, maximum mean discrepancy.

1 Introduction↩︎

Numerical integration, or quadrature, is of routine use in applied mathematics, such as in Bayesian inference [1]. For a target measure \(\pi\) on \(\mathbb{R}^d\), a quadrature is a measure \(\mu\) supported on a finite number \(n\) of atoms such that one can bound the worst-case integration error \[\label{e:worst95case95error} \mathsf{E}_\mathcal{F}(\mu) = \sup_{f\in\mathcal{F}} \left\vert \int f\mathrm{d}(\pi-\mu)\right\vert\tag{1}\] for a class of functions \(\mathcal{F} \subset L^1(\pi)\). A wealth of constructions have been investigated in moderate and large dimension. For instance, Monte Carlo methods correspond to taking \(\mu\) to be a random measure; they come with strong probabilistic controls on 1 when \(\mathcal{F}\) is a singleton [2]. On the other hand, quasi-Monte Carlo methods correspond to deterministic quadratures that come with uniform controls on a larger set \(\mathcal{F}\) of smooth functions [3].

More recently, an intermediate string of works has been approximating \(\pi\) by a minimizer of the so-called maximum mean discrepancy (MMD) \[\begin{align} \mu \mapsto I_{K}(\mu - \pi):= \iint K(x, y)\,\mathrm{d}(\mu -\pi)(x)\,\mathrm{d}(\mu - \pi)(y), \label{e:mmd} \end{align}\tag{2}\] where \(K\) is a suitable bivariate function called an interaction kernel, and the minimization is carried over measures supported on \(n\) points. The motivation is twofold. First, when \(K\) is bounded, one can show that 2 is the worst-case error 1 where the set \(\mathcal{F}\) of integrands is the unit ball of the reproducing kernel Hilbert space (RKHS) \(\mathcal{H}_K\) of kernel \(K\) [4], see 1. Throughout the paper, we will abusively denote by \(\mathrm{E}_{\mathcal{H}_K}\) the worst–case error in 1 when \(\mathcal{F}\) is the unit ball of \(\mathcal{H}_K\). Second, as long as the kernel mean embedding (also known as potential) \[\label{e:potential} U_K^\pi:z \mapsto \int K(z, x)\mathrm{d}\pi(x)\tag{3}\] can be numerically evaluated at any point of \(\mathbb{R}^d\), 2 can be evaluated as well, which makes its minimization amenable to optimization algorithms. At the price of the evaluation of 3 , MMD minimization seems to offer a compromise between fast scalable algorithms and small integration errors over an RKHS [5]–[14]; see [14] for an overview of existing methods. Very recently, [15] have also shown that, while minimizing the MMD is a nonconvex problem, stationary points (instead of minimizers) of the MMD actually already lead to efficient quadrature, while being more computationally tractable.

Randomized quadratures based on determinantal point processes [16], [17] have also been found to control the MMD [18], [19]. In the same spirit, we proposed in [20] a direct probabilistic relaxation of MMD minimization, where we draw the atoms \(y_1, \dots, y_n\) of our quadrature from the (Gibbs) measure on \(\left(\mathbb{R}^{d}\right)^{n}\) defined as \[\begin{align} \label{e:gibbs95measure} \mathrm{d}\mathbb{P}_{n, \beta_n}^{V}(y_1, \dots, y_n) = \frac{1}{Z_{n, \beta_n}^{V}}\exp\left(-\beta_n H_n(y_1, \dots, y_n)\right)\,\mathrm{d}y_1\,\dots\mathrm{d}y_n, \end{align}\tag{4}\] where \(\beta_n>0\) is a parameter called the inverse temperature, \[\begin{align} \label{e:H95n} H_n(y_1, \dots, y_n) = \frac{1}{2 n^2}\sum_{i \neq j} K(y_i, y_j) + \frac{1}{n}\sum_{i = 1}^{n}V(y_i), \end{align}\tag{5}\] includes an interaction term that is an empirical proxy for 2 , \(V:\mathbb{R}^d\mapsto \mathbb{R}\) is a potential that goes to infinity as \(\Vert y\Vert\rightarrow \infty\) sufficiently fast to confine the nodes \(y_1, \dots, y_n\), i.e. avoid the density in 4 to be maximal when the \(y_i\)s diverge to \(\infty\), and \(Z_{n, \beta_n}^{V}\) is a normalizing constant.

The family of distributions 4 encompasses well-known point processes. For instance, when \[K(x,y) = g(x,y) := \left\{ \begin{array}{ll} \frac{1}{\lvert x -y\rvert^{d-2}}, & \textrm{ if } d\ge 3,\\ -\log|x-y|, & \textrm{ if } d=2, \end{array} \right.\] the Gibbs measure 4 is called a Coulomb gas, or one-component plasma [21], [22]. It models the electrostatic interaction of unit-charge elementary particles in \(\mathbb{R}^{d}\) that are further confined by an external field \(V\). Alternatively, when \(d=1\) and \(K(x,y) = -\log|x-y|\), the measure 4 is called one-dimensional log-gas. It appears naturally as the distribution of the eigenvalues of some classical ensembles of random matrices [23].

In general, the potential \(V\) in 4 has to be confining enough to make the probability measure well-defined, i.e. to guarantee that \(0 < Z_{n,\beta_n}^V<\infty.\) A particularly interesting choice is \[\begin{align} \label{e:jellium} V(z) = V^{\pi}(z) := - U_{K}^{\pi}(z) + \Phi(z), \quad z\in\mathbb{R}^d, \end{align}\tag{6}\] where \(\Phi\) is a simple confining term that vanishes on the (known) support of \(\pi\), and \(U_{K}^{\pi}\) is the potential of \(\pi\) with respect to \(K\) given by 3 . In that case, it is easy to check through the Euler-Lagrange equations that \(\pi\) is the equilibrium measure of the system, in the sense that the empirical measure \(\mu_n\) of the particles \(y_1, \dots, y_n\) converges to \(\pi\) as \(n\) goes to infinity, making \(\mu_n\) a candidate quadrature rule in numerical integration w.r.t. \(\pi\). Additionally, for a smooth function \(f\), the Monte Carlo estimator \(\mu_n(f)\) is expected to have low mean-square error because the particles spread very regularly through \(\mathbb{R}^d\), avoiding the quadrature to miss any relevant oscillation by leaving an uncovered hole in \(\mathbb{R}^d\). Relatedly, in physics, when \(K=g\) and \(V=V^\pi\), a random configuration drawn from 4 is called jellium, and one can imagine positive electric charges interacting with a negatively charged background distribution \(\pi\).

One way to quantify the fast convergence to \(\pi\) of the empirical measure of \(y_1, \dots, y_n\) drawn from \(\eqref{e:gibbs95measure}\) with \(V=V^\pi\) in 6 is through large deviations principles (LDPs). In particular, [24] and [25] established that for a very general class of bounded interaction kernels \(K\), the worst-case integration error on the unit ball of the RKHS \(\mathcal{H}_K\) is small with overwhelming probability, i.e. for any \(r > 0\) and \(\epsilon > 0\), for large enough \(n\), \[\begin{align} \mathbb{P}_{n, \beta_n}^{V^{\pi}}\left(\mathsf{E}_{\mathcal{H}_K}\left(\frac{1}{n}\sum_{i = 1}^{n}\delta_{y_i}\right) > r\right) \leq \exp\left(- \left(\frac{1}{2}- \epsilon\right) \beta_n r^2\right), \label{e:ldp95gibbs95bounded} \end{align}\tag{7}\] Similarly, taking for \(K\) the (unbounded) Coulomb interaction 12 , they prove \[\begin{align} \mathbb{P}_{n, \beta_n}^{V^{\pi}}\left(\mathsf{E}_{\mathcal{BL}}\left(\frac{1}{n}\sum_{i = 1}^{n}\delta_{y_i}\right) > r\right) \leq \exp\left(- c \beta_n r^2\right), \label{e:ldp95gibbs} \end{align}\tag{8}\] for some constant \(c>0\), where \(\mathcal{BL}\) is the space of functions that are \(1\)-Lipschitz and bounded by \(1\). Throughout the paper, when the pushforward of the Gibbs measure \(\mathbb{P}_{n, \beta_n}^{V^{\pi}}\) by the application \[y_1, \dots, y_n \mapsto \frac{1}{n}\sum_{i = 1}^{n}\delta_{y_i}\] satisfies a LDP, we will abusively say that the sequence of empirical measures \(\xi_n = \frac{1}{n}\sum_{i = 1}^{n}\delta_{y_i}\) satisfies the same LDP. The results 7 and 8 are to compare to classical results for i.i.d. quadrature nodes, such as Sanov’s theorem, where \(\beta_n\) in the exponential is replaced by \(n\). The same rate \(n\) appears for Markov chain quadratures [26]. As soon as \(\beta_n\) grows faster than \(n\), the right-hand sides of 7 and 8 thus decrease faster than i.i.d. or Markov chain quadratures, justifying our use of the terms overwhelming probability. This fast convergence motivated our study of numerical quadratures based on the Gibbs measure \(\mathbb{P}_{n, \beta_n}^{V^{\pi}}\) in [20].

There are two main issues with this strategy, though. First, we do not know of any exact sampling algorithm for the Gibbs measure 4 , so that one has to rely on approximate samplers like Markov chain Monte Carlo (MCMC) methods. Second, and maybe more fundamentally, computing the potential 6 requires evaluating \(U_{K}^{\pi}\), which requires integrating w.r.t. \(\pi\), seemingly defeating our initial purpose. In this paper, we address this second issue, completely removing the need to evaluate \(U_K^\pi\) while preserving the better-than-Monte-Carlo integration guarantees of Gibbs measures.

Because the problem of evaluating \(U_K^\pi\) is common to all algorithms that minimize the MMD 2 , workarounds have naturally already been investigated. They come in two clusters. First, one can rely on Stein’s trick [27], which defines a kernel such that \(U_{K}^{\pi} \equiv 0\), but modifies the set of functions \(f\) that are efficiently integrated in an unclear way, see related work below. Alternatively, the majority of MMD minimization algorithms [6]–[8], [28] builds a Monte Carlo approximation to 3 , using a single draw of a random measure that approximates \(\pi\), e.g. using MCMC. If that random measure has a large number of atoms compared to the cardinality of the final quadrature \(n\), this two-level Monte Carlo technique reduces the range of potential applications by forcing the user to spend a significant part of their computational budget on preparing the distribution from which they will draw their final quadrature. Apart from [13] and [14], who bypass the issue by controlling the difference in MMD of their \(n\)-point quadrature to an \(n^2\)-point MCMC reference, to our knowledge there is no guarantee on the worst-case integration error 1 of any MMD minimization algorithm that includes this second level of approximation of 3 . Our contribution is to provide such a guarantee for Gibbs measures.

Formally, let \(\nu_n\) be a random \(n\)-atom measure on \(\mathbb{R}^d\), called the background. Conditionally on \(\nu_n\), let \(y_1, \dots, y_n\) be drawn from the Gibbs measure 4 with confining potential \[\label{e:random95potential} V^{\nu_n} = - U_K^{\nu_n} + \Phi,\tag{9}\] that is 3 where we replace the target measure \(\pi\) by the atomic measure \(\nu_n\). Following the vocabulary of large deviations theory [25], we call the conditional distribution \(\mathbb{P}_{n, \beta_n}^{\mathrm{Q}, \pi}\) of \(y_1, \dots, y_n\) given \(\nu_n\) a quenched Gibbs measure. In this paper, we show large deviations principles for this quenched Gibbs measure.

Our first result deals with a bounded interaction kernel \(K\) in the Gibbs measure. We show in 3 that, on any event such that \(\nu_n\) converges to \(\pi\), the distribution \(\mathbb{P}_{n, \beta_n}^{\mathrm{Q}, \pi}\) of the quadrature \(y_1, \dots, y_n\) almost surely satisfies the fast LDP 7 , just like the Gibbs measure with intractable potential 3 . Consequently, we avoid the evaluation of 3 –the bottleneck of most MMD minimization methods– at no statistical cost on the LDP. Note that it is computationally easy to build a \(\nu_n\) that converges almost surely to \(\pi\) using standard Monte Carlo techniques, such as self-normalized importance sampling.

Our second result covers the more difficult case of the (singular) Coulomb kernel \(K=g\) in 4 . This time, we show that the quenched Gibbs measure almost surely preserves the LDP 8 if the atoms of \(\nu_n\) are themselves drawn from a Gibbs measure with Coulomb interaction, under weak requirements on the equilibrium measure, and with importance weights compensating for the mismatch between that equilibrium measure and the target \(\pi\). Here again, we obtain the same statistical guarantee on the worst-case integration error without ever needing to evaluate the intractable potential 3 . Seeing Gibbs measures as an avatar of MMD minimization, we believe that ours is the first work that formally includes the approximation of the kernel mean embedding \(U_{K}^{\pi}\) in the analysis of its quadrature. We hope that this can inspire the analysis of non-Gibbs-based MMD minimization algorithms.

As a side remark for the reader versed into the literature on the Coulomb interaction, we provide en route a control on the uniform convergence of the approximation of the potential. For a well-chosen regularization \(K_\zeta\) of the Coulomb kernel, we show in 6 and 7 that \[\underset{z \in \mathbb{R}^{d}}{\sup\,} \left\vert\nu_n(K_{\zeta}(z, \cdot)) - U_{g}^{\pi}(z)\right\vert \rightarrow 0\] almost surely as \(n \rightarrow + \infty\). We believe this intermediate result is of independent interest since the extrema of the potential field generated by a Coulomb gas have recently drawn attention, with exact scalings obtained at first order in dimension \(2\) [29], [30], more on this below.

The paper is organized as follows. In the end of this section, we quickly survey related work. In 2, we rigorously introduce the objects of interest and our assumptions, and we state our main results. Sections 4 and 5 are devoted to the proofs of the main results. Numerical experiments, discussion on remaining practical issues and proofs of auxiliary lemmas are given in appendix.

Related work↩︎

Large deviations↩︎

Large deviations for Gibbs measures with pairwise interactions as in 4 have been thoroughly studied by [24], [25]. In our work, we need to further include a random environment, i.e. an empirical approximation to the potential.

Actually, [25] give a list of sufficient conditions to derive large deviations principles in varying environments, with applications to Gibbs measures. Yet, the latter applications deal with Gibbs measures for which we cannot pre-specify the equilibrium measure, an important requirement for Monte Carlo integration where the equilibrium measure should be (at least an approximation of) the target integration measure. Our proofs thus consist in checking the sufficient conditions of [25] for the potential 9 . Checking these conditions in our setting is however no sinecure, especially in the case of the Coulomb interaction kernel.

Extrema of the potential of a Coulomb gas↩︎

In dimension \(1\) and \(2\), the two-dimensional Coulomb interaction kernel \(g(x, y) = -\log \lvert x - y\rvert\), along with \(V(x) = \lvert x \rvert^2 /2\) and \(\beta_n=2n^2\), yields a Gibbs measure that also describes the distribution of the eigenvalues of random matrices, respectively the Gaussian unitary and Ginibre ensembles [23]. In this case, the interaction potential \(\sum_{i = 1}^{n} g(z, x_i)\) generated at point \(z\) by points \(x_1, \dots, x_n\) drawn from 4 coincides with the log of the modulus of the characteristic polynomial of the random matrix, evaluated in \(z\). The behavior of the maximum of this quantity has been a topic of interest, following the famous conjectures of [31] in the one-dimensional case. In dimension \(2\), for general \(V\) and \(\beta_n = \beta n^2\), it has been for instance shown by [29], [30], [32] that \[\begin{align} \underset{z \in D(x, r)}{\sup\,} \frac{n}{\log n}\left\{ \frac{1}{n}\sum_{i = 1}^{n}g(z, x_i)- \int g(z, x)\,\mathrm{d}\mu_V(x) \right\} \underset{n \rightarrow + \infty}{\rightarrow} \frac{1}{\sqrt{\beta}} \text{\, in probability,}\label{e:lln95max95pot} \end{align}\tag{10}\] where \(\mu_V\) is the compactly supported measure that satisfies 6 for a fixed \(V\) and \(D(x, r)\) is any disk strictly contained in the support of \(\mu_V\). Obtaining a similar behavior in dimension \(d\geq3\) is still an open question up to our knowledge, but we show as a byproduct of our proof in 6 that for general \(V\) and \(\beta_n\) growing faster than \(n\), \[\underset{z \in \mathbb{R}^{d}}{\sup\,}\left\lvert \frac{1}{n}\sum_{i = 1}^{n} K_{\zeta}(z, x_i) - \int g(z, x)\,\mathrm{d}\mu_V(x)\right\rvert \underset{n \rightarrow + \infty}{\rightarrow} 0, \quad\text{\mathbb{P}-almost surely,}\] where \(K_{\zeta}\) is a regularized version of \(g\). Regularization is indeed needed to keep the supremum of the absolute value finite. Obtaining the exact scaling for the leading order as in 10 requires fine estimates on the fluctuations of the system at the second order, a hot topic in the field for \(d \geq 3\).

Monte Carlo integration and KSDs↩︎

Computing \(U_{K}^{\pi}(z)\) has been up to our knowledge one of main limitations of MMD minimization algorithms. One standard way to bypass this computation is to rely on a large MCMC chain targeting \(\pi\) [8], which is fair if the computational bottleneck is the number \(n\) of evaluations of the test function \(f\). This can be the case for instance in computational biology when each evaluation of the likelihood of the model can range from a few minutes in cardiac electrophysiology models [33] to a few hours for population models in ecology [34]. Yet, this comes at the price of losing theoretical integration guarantees for Monte Carlo integration with respect to \(\pi\).

Kernel Stein Discrepancies (KSDs, [27]) have been introduced as a way to bypass this computation. Indeed, given a kernel \(K\) and a target measure \(\pi\), one can define a new kernel \(K_{\pi}\) with closed form such that \(U_{K_{\pi}}^{\pi}(z) = 0\) for all \(z\) under a mild regularity assumption on \(K\) when \(\pi\) vanishes. The latter can always be ensured even when \(\pi\) does not vanish up to multiplying the kernel \(K\) by a vanishing weight function at the boundary of the integration domain [35]. The kernel Stein discrepancy is then defined as the square root of the MMD in 2 with kernel \(K_{\pi}\). Minimizing this quantity has led to fast and practical sampling algorithms [10]–[12] along with decay rates in \(n^{-1/2}\) for the KSD with target \(\pi\) in some cases. Yet, replacing \(K\) by \(K_{\pi}\) replaces the set of test functions \(\mathcal{H}_K\) by \(\mathcal{H}_{K_\pi}\), and there is little knowledge about the kind of functions contained in \(\mathcal{H}_{K_\pi}\). There is evidence that the quadrature is influenced in that case by pathological test functions that oscillate wildly in regions where \(\pi\) puts little mass [36]. Even though one is able to show [37], [38] that under some assumptions on \(K\), \(\mathcal{H}_{K_\pi}\) is still large enough to imply convergence in Wassertein or bounded Lipschitz distance, we only have access so far to estimates of the form \[\mathsf{E}_{\mathcal{BL}}(\mu) \in \mathcal{O}\left(\left(\log \frac{1}{\mathrm{KSD}(\mu, \pi)}\right)^{-1/2}\right)\] as the \(\mathrm{KSD}\) decreases. Thus, the decay rates in \(\mathcal{H}_{K_{\pi}}\) only transfer as slow logarithmic decay rates for meaningful distances for Monte Carlo integration.

2 Definitions and main results↩︎

Let \(d \geq 1\) and \(\pi\) be a probability measure on \(\mathbb{R}^{d}\), which we call the target distribution. Recall that our goal is to get a computationally tractable empirical measure supported on \(n\) points that approximates \(\pi\).

As explained in 1, we will build our quadrature by balancing two effects: a repulsive pairwise interaction preventing points to be too close from one another, and a confinement term favoring regions of large mass under \(\pi\). These two effects are modeled by quantities known as energies and interaction kernels, and are combined in a Gibbs measure.

Throughout the paper, some particular kernels will play an important role.

For a kernel \(K\) satisfying 2, the interaction energy between two distributions coincides with the squared worst-case integration error for test functions in the reproducing kernel Hilbert space \(\mathcal{H}_K\) of kernel \(K\).

Non-singular kernels such as the Gaussian or Matern kernels are standard in kernel regression [39]. Other natural examples come from regularizing singular kernels. For instance, for \(\epsilon > 0\) and \(s > 0\), consider the regularized Riesz interaction \[\begin{align} \label{e:riesz95kernel95reg} K_{s, \epsilon}(x, y) = \frac{1}{\left(\lvert x -y \rvert^2 + \epsilon^2\right)^{s/2}}. \end{align}\tag{13}\] In particular, for \(n \geq 1\) and \(\zeta > 0\), we will also use the following non-singular regularization of the Coulomb interaction \(g\) in 12 , \[\begin{align} \label{e:coulomb95kernel95reg} K_{\zeta}(x, y) = \frac{1}{\left(\lvert x -y \rvert^2 + n^{-2\zeta}\right)^{\frac{d-2}{2}}}. \end{align}\tag{14}\]

We will work with the following assumptions on the confinement \(V\).

We quickly give basic facts concerning \(V\) and its link with the minimizers of the energies in 1, see for instance [8], [24] for references. Under 3, the energies of 1 are well defined when \(K\) either satisfies 2 or is the Coulomb kernel 12 , though the energy is possibly infinite in the Coulomb case. In the same way, \(I_{K}^{V}\) in 11  is known to have a unique minimizer over the set of probability distributions, called the equilibrium measure \(\mu_V\). Furthermore, \(\mu_V\) then has finite energy and compact support.

Although the so-called Euler-Lagrange equations characterize the equilibrium measure for general interaction kernels, the computation of the equilibrium measure is not explicit in general. Yet there are exceptions. First, when \(K = g\), a few choices of potential \(V\) do lead to explicit expressions. For instance, for \(V(x) = \frac{d-2}{2 R^d}\lvert x \rvert^{2}\), \(\mu_V\) is the uniform measure on \(B(0, R)\). This is the main reason why the Coulomb kernel \(g\) will play an important role in the upcoming results, although it involves some additional theoretical difficulties. Second, and of paramount important for our Monte Carlo ambitions, we can force the equilibrium measure to be a given distribution \(\pi\). Indeed, under 1 and when \(K\) either satisfies 2 or is the Coulomb kernel 12 , setting \[V^{\pi}(z) := \begin{cases} - U_{K}^{\pi}(z) \text{ for } z \in B(0, R)\\ - U_{K}^{\pi}(z) + \Phi(z) \text{ otherwise} \end{cases}, \label{e:euler95lagrange}\tag{15}\] with for instance \(\Phi(z) = \left[\lvert z\rvert^{2} - R^2\right]_{+},\) yields that \(V^{\pi}\) satisfies 3 and that \(\mu_{V^{\pi}} = \pi.\) We could make a more general choice of \(\Phi\) to confine the points in \(B(0, R)\), but the quadratic choice is quite canonical. Additionally, note that with \(V=V^\pi\), the energy \(I_{K}^{V^{\pi}}(\mu)\) in 1 writes as \[I_{K}^{V^{\pi}}(\mu) = \frac{1}{2}I_{K}(\mu - \pi) - \frac{1}{2}I_{K}(\pi) + \int_{z \notin B(0, R)} \Phi(z)\,\mathrm{d}\mu(z)\] when it is finite, so that minimizing this quantity is indeed natural if one wishes the distribution \(\mu\) to be close to \(\pi\). The reader should keep in mind that evaluating \(V^\pi\) in 15 requires evaluating \(U_K^\pi\), which is itself an integral w.r.t. \(\pi\); this a priori forbids the use of \(V^\pi\) in a Monte Carlo integration algorithm with target distribution \(\pi\).

We now define the probability measure from which we want the points to be drawn from, promoting configurations of points with small energy \(I_{K}^{V}\) up to diagonal terms.

We recall a few results from [24]. Under the assumptions of 3, \(Z_{n, \beta_n}^{V}\) is positive and finite, making the Gibbs measure \(\mathbb{P}_{n, \beta_n}^{V}\) a well-defined probability distribution.³ Moreover, if \(x_1, \dots, x_n \sim \mathbb{P}_{n, \beta_n}^{V}\) and \(\xi_n = \frac{1}{n} \sum_{i=1}^n \delta_{x_i}\), it is known that in the low-temperature regime \(\beta_n / n \longrightarrow + \infty\) which we denote by \(\beta_n \gg n\), \(\xi_n\) converges in law to \(\mu_V\), \(\mathbb{P}\)-almost surely, as \(n\) grows. Actually, in the same low-temperature regime, \(\xi_n\) satisfies a large deviations principle (LDP) with speed \(\beta_n\) and rate function \(\mu \mapsto I_{K}^{V}(\mu) - I_{K}^{V}(\mu_V)\), meaning in particular that for any \(r > 0\), there exists \(n_0\) such that for any \(n \geq n_0\), \[\begin{align} \mathbb{P}_{n, \beta_n}^{V}\left(\xi_n \notin B_{\mathbf{d}}(\mu_V, r)\right) \leq \exp\left(-\beta_n \underset{\mu \not\in B_{\mathbf{d}}(\mu_V, r)}{\inf\,} \{I_{K}^{V}(\mu) - I_{K}^{V}(\mu_V)\}\right)\label{e:bound95gibbs} \end{align}\tag{16}\] where \(B_{\mathbf{d}}(\mu_V, r) := \{ \nu \in \mathcal{P}(\mathbb{R}^d), \mathbf{d}(\mu_V, \nu) < r\}\), and \(\mathbf{d}\) is any distance on probability distributions compatible with the convergence in distribution. In particular this implies 8 and 7 .

Recall that our motivation of numerical integration w.r.t. \(\pi\) hints at the choice \(V = V^{\pi}\) as in 15 , which however involves an intractable integral with respect to \(\pi\). We will replace this integral itself by an easier-to-draw Monte Carlo estimator \(V_n\) of \(V^\pi\), using a quadrature \(\nu_n\). Conditionally on \(V_n\), we call the Gibbs measure with potential \(V_n\) the quenched Gibbs measure.

The atoms \(x_1, \dots, x_n\) of \(\nu_n\) in 4 act as a random environment (or background) that approximates \(\pi\). Then, the new points \(y_1, \dots, y_n \sim \mathbb{P}_{n, \beta_n}^{\mathrm{Q}, \pi}\) repel each other, while the random environment acts as a confinement, attracting the nodes \(y_1, \dots y_n\) towards regions of high probability under the target distribution \(\pi\). Crucially, all of this happens without ever needing to evaluate the intractable potential \(U_{K}^{\pi}\) in 15 . Our goal is now to show that for \(\mathbb{P}\)-almost-all background realization, the quenched empirical measure still converges to the target distribution \(\pi\) at the same speed as 16 . Without further assumption, we can already state our main result for bounded kernels.

The proof of the LDP is in 4. The first consequence is immediate, and the second one follows from ?? , along with the fact that \(I_{K}(\mu_{n}^{\mathrm{Q}} - \pi) \leq 2(I_{K}^{V_{\pi}}(\mu_{n}^{\mathrm{Q}}) - I_{K}^{V^{\pi}}(\pi))\); see [20]. Taking a step back, 3 yields tighter simultaneous confidence intervals for Monte Carlo integration of functions in the unit ball of \(\mathcal{H}_{K}\) than with points drawn i.i.d from \(\pi\) or using MCMC, yet only asymptotically. We stress the fact that we put no condition on the background \(\nu_n\), other than the almost sure convergence to \(\pi\). Our LDP is thus valid even in the case where \(\nu_n\) is the empirical measure of the history of a suitable MCMC chain up to time \(n\). This justifies a common practical heuristics of replacing the intractable potential by a Monte Carlo estimate obtained using a long MCMC run.

From 3, one might consider taking an inverse temperature \(\beta_n\) as large as possible in order to maximize the gain in the tail bound. Yet, one still has to produce samples from the quenched Gibbs measure with inverse temperature \(\beta_n\) in order to provide Monte Carlo estimators. Heuristically, this becomes more and more difficult as \(\beta_n\) becomes large as particles drawn from the Gibbs measure then try to get close to the deterministic minimizers of the energy \(H_{n}^{\mathrm{Q}}\) in 4. We thus expect a tradeoff between sampling accuracy or complexity and tail bound guarantees governed by \(\beta_n\). However, there is to our knowledge no algorithm to sample approximately from 4 that comes with theoretical considerations on this tradeoff, so we give only practical considerations in 6.1.

Motivated by simultaneous confidence intervals for the wider class of bounded Lipschitz functions, we now turn to our second main result, covering the case when \(K=g\) is the Coulomb kernel and \(d\geq 3\). We shall need another set of assumptions, always in conjunction with 1 that guarantees that \(\pi\) is compactly supported.

When \(V\) is the quadratic potential \(V(z) = \frac{d-2}{2R^d}\lvert z\rvert^2\), for instance, 4 and 5 are satisfied as soon as \(S_{\pi} \subset B(0, R)\) since \(\mu_V\) is the uniform measure on \(B(0, R)\). We shall also need weights in our Monte Carlo estimator of the potential in the case \(K =g\).

We have the following convergence for the weighted empirical measures.

Note that since computing \(w\) is invariant by multiplication of \(\pi'\) by a constant, \(\nu_n\) does not depend on the normalization constant of the target \(\pi\), and so neither does any linear statistic w.r.t \(\nu_n\). In the following theorem, we chose \(x_1, \dots, x_n\) as a realization of a Coulomb gas with external potential \(V\), and we show that for this choice of \(\nu_n\) in ?? , we have the desired LDP for the quenched Gibbs measure.

The proof is again in 4. The first consequence is immediate, and the second one follows from the Coulomb transport inequality [40]. Again 5 gives tighter asymptotic simultaneous asymptotic confidence intervals for Monte Carlo integration of bounded Lipschitz functions than with points drawn i.i.d from \(\pi\) or using MCMC. In that case however, we need the first approximation \(\nu_n\) to \(\pi\) to be built as well on a Coulomb gas, i.e. to exhibit fast converging properties, this is due to the fact that the quality of \(\nu_n\) interplays with the regularization that we need to make to deal with the singularity of \(g\). Moreover, the fact that the regularized Coulomb interaction 14 satisfies 2 will play an important role in the proof. In particular, in dimensions \(1\) and \(2\), the natural counterpart to 14 is given by \[K_{\zeta}'(x, y) = -\frac{1}{2}\log \left(\left\lvert x - y \right\rvert^2 + n^{-2\zeta}\right),\] which fails to satisfy 2 since \(I_{K}(\mu) > 0\) only for finite signed Borel measures with total mass \(0\) [8]. Still, it would be interesting to try to circumvent this issue in the proofs.

We will use a general result on LDPs for Gibbs measures from [25] to get 3 and 5. The case of non singular interactions is relatively easy to handle, whereas more control is needed in the Coulomb case because of the regularization. In particular, the key control will be given by 7. In the latter, we prove uniform convergence of the regularized potential generated by the Coulomb gas \(x_1, \dots, x_n\). We first begin with an unweighted version of the desired result, which we will use in the proof of the weighted version and which we believe to be interesting as discussed in 1.

We give a sketch of proof, to be fleshed out in 5.1. We first note that \(K_{\zeta}\) satisfies 2, so that \(K_{\zeta}(z, .)\) is in the RKHS \(\mathcal{H}_{K_\zeta}\) for all \(z\) and has squared norm \(n^{\zeta (d-2)}\) in \(\mathcal{H}_{K_\zeta}\). In particular, the duality formula ?? yields \[\underset{z \in \mathbb{R}^{d}}{\sup}\,\left\lvert U_{K_{\zeta}}^{\mu_{n}}(z) - U_{K_{\zeta}}^{\mu_V}(z)\right\rvert \leq \sqrt{n^{\zeta(d-2)} I_{K_{\zeta}}(\mu_{n}-\mu_V)}.\] Then, one can show (3) that \(I_{K_{\zeta}}\) is always smaller than \(I_g\). Since the points \(x_1, \dots, x_n\) are generated from a Coulomb gas with interaction \(g\), one should thus be able to use concentration properties of Gibbs measures (4) to show that \(n^{\zeta(d-2)} I_g\) goes to zero. However, \(g\) is not integrable at infinity, so that \(I_g(\mu_n - \mu_V)\) is infinite. We bypass this issue by convolving each Dirac \(\delta_{x_i}\) in \(\mu_n\) with a uniform distribution on a small ball of radius \(n^{-\epsilon}\) to make sense of the energy \(I_g\). Denoting by \(\mu_{n}^{(\epsilon)}\) this convolved empirical measure, we thus show that \[n^{\zeta (d-2)} I_{H_{\zeta}}(\mu_{n}^{(\epsilon)} - \mu_V) \leq n^{\zeta(d-2)}I_g(\mu_{n}^{(\epsilon)} - \mu_V) \rightarrow 0\] for a suitable choice of \(\epsilon\) and \(\zeta\). It then remains to show that the error terms vanishes when one replaces \(U_{K_{\zeta}}^{\mu_n} - U_{g}^{\mu_V}\) by \(U_{K_{\zeta}}^{\mu_{n}^{(\epsilon)}} - U_{K_{\zeta}}^{\mu_V}\), which we do in 1 and 2.

We now state the same result for the weighted measures, the proof of which relies on 6. In the proof of 5, we will use the following:

The constraint on \(\zeta\) is upper bounded by \(\mathcal{O}(1/d^2)\) both in 6 and 7, but it is not expected to be optimal at all.

We give again a sketch of proof, with details to be found in 7. The beginning of the proof is in the same vein as 6: we show using 2 and 5 that it is enough (i.e. the error terms vanish) to study \(U_{K_{\zeta}}^{\nu_{n}^{(\epsilon)}} - U_{K_\zeta}^{\pi}\) where \(\nu_{n}^{(\epsilon)}\) is the weighted empirical measure with each Dirac convolved by a uniform distribution on \(B(0, n^{-\epsilon})\). For that purpose, we need instead of \(\nu_{n}^{(\epsilon)}\) to consider \(\nu_n^{\epsilon, \eta}\) in which the weights are smoothed versions of W obtained by convolving with a Gaussian mollifier. Then, one wishes to decouple the weights from the unweighted measures \(\mu_{n}^{(\epsilon)}\) and \(\mu_V\) to apply once again the concentration properties for the Gibbs measure from which the points are generated. To do so, a key idea is to rephrase the energy \(I_{K}(\mu - \nu)\) in Fourier space as an \(L^2\)-norm between the characteristic functions of \(\mu\) and \(\nu\) with respect to a spectral measure induced by \(K\). For \(K = K_{\zeta}\) and \(K = g\), we give explicit expressions and useful estimates for this spectral measure in 7.1. We then show in 6 that we can split this \(L^2\)-norm in two terms: one vanishing term and the other one controlled by quantities of the form \(n^{2\zeta d} I_g(\mu_{n}^{(\epsilon)}-\mu_V)\) which converges to \(0\) again for suitable choices of \(\epsilon\) and \(\zeta\).

3 Discussion and perspectives↩︎

We have shown that a suitable random approximation of the potential can be integrated into large deviations results that support the use of Gibbs measures for Monte Carlo integration. Several open questions deserve further investigation.

Extensions↩︎

First, for exhaustiveness, generalizing 5 to dimensions \(1\) and \(2\) could be of interest. The difficulty is that the natural replacement for the regularized Coulomb kernel 14 in dimensions 1 and 2 is not positive semidefinite.

A second interesting question lies within the choice of the size of the background \(\nu_n\). Our 3 is only asymptotic, and does not bring to light the impact of the quality of \(\nu_n\); yet the latter matters in practice, as we experimentally show in 6. Quantifying this impact at fixed \(n\) would require a concentration inequality in the spirit of 3. For the record, one can easily show using the computations of [20] that in the setup of 3, there exists a constant \(c > 0\) such that for \(\mathbb{P}\)-almost all \(\omega\) and for any \(n \geq 2\) and \(r > n^{-1/2}\), \[\label{e:concentration95quenched} \mathbb{P}_{n, \beta_n}^{\mathrm{Q}}\left(\mathrm{E}_{\mathcal{H}_K}\left(\frac{1}{n}\sum_{i = 1}^{n}\delta_{y_i}\right) > r\right) \leq \exp\left(-c \beta_n r^2 + 2 \beta_n \,\mathrm{E}_{\mathcal{H}_K}\left(\nu_n(\omega)\right)\right).\tag{22}\] If \(\nu_n\) is the empirical measure of the history of a Markov chain of size \(M_n\), it is known [13] that \(\mathrm{E}_{\mathcal{H}_K}\left(\nu_n\right) = \mathcal{O}(M_{n}^{-1/2})\) with high probability. Therefore, there exists constants \(u_0, u_1\) such that for any \(r > u_0 (n^{-1/2} \vee M_{n}^{-1/4})\), with high probability, \[\mathbb{P}_{n, \beta_n}^{\mathrm{Q}}\left(\mathrm{E}_{\mathcal{H}_K}\left(\frac{1}{n}\sum_{i = 1}^{n}\delta_{y_i}\right) > r\right) \leq \exp\left(-u_1 \beta_n r^2 \right).\] When \(M_n = n\) for instance, this suggests a downgrade in the constraint on \(r\) compared to the concentration result of [20] with exact knowledge of the potential \(U_{K}^{\pi}\), while there is no downgrade when using a larger initial Markov chain \(M_n \geq n^2\). In line with [13], [14], this seems to point out that the background should be a significantly longer Markov chain than the cardinality of the final quadrature. Alternatively, one could choose a repulsive background with variance reduction, such as the Coulomb gas we used in 5. Still, these recommendations are to be taken with a pinch of salt as our bound 22 might be crude.

A fast CLT for linear statistics↩︎

A stronger argument in favor of using Gibbs measures for Monte Carlo integration would be a central limit theorem with a fast rate. While obtaining such a result for the singular Coulomb interaction is notoriously difficult [41], one can reasonably hope that this would be easier for a bounded interaction kernel.

Sampling from the Gibbs measure↩︎

Now that we have incorporated the approximation of the potential into a fast Monte Carlo result, the elephant in the room is our need for an MCMC algorithm to approximately sample from our Gibbs measures in reasonable time. The experiments given in 6 show that this can be an issue if one wants to see smaller integration errors or variances than MCMC in practice. In the absence of an exact sampler, understanding the sampling complexity of an MCMC algorithm like MALA for this specific target in terms of the number of iterations \(T\), the step size, the dimension \(d\), the number of particles \(n\), and the inverse temperature \(\beta_n\) would be important, and it would be interesting to compare the different methods we considered in 6 with fixed budget. Additionally, the recent work of [15] suggests that stationary points of the MMD are both computationally tractable and lead to efficient numerical integration of the functions in a large subset of the underlying RKHS. For integrating these particular functions, deterministic optimization of the MMD might thus be preferable. Alternatively, when MALA fails to reach the targeted Gibbs distribution because the gradient of the log target is too small, we may thus be in the vicinity of a stationary point of the MMD, which in turn would not be bad news in terms of numerical integration. Finally, sampling algorithms other than MALA, tailored to our Gibbs measures, could also achieve better performances. One could as well think of drawing points from a discrete Gibbs measure supported on a pre-drawn background \(\nu_n\), following the two-step rationale of coresets or kernel thinning [13].

4 Proofs: quenched large deviations↩︎

In [25], the author established a LDP for Gibbs measures in a quite general framework. In 4.1, we recall his results in the case when \(\frac{\beta_n}{n}\) converges to infinity (that we denote by \(\beta_n \gg n\)) and detail in particular the conditions that we need to check to apply the results in our context. In 4.2, we check these conditions in the case of non-singular kernels and conclude the proof of Theorem 3. To further prove 5, where the kernel is the Coulomb kernel \(g\), we will need 7, whose proof is the main technical input of the present paper and is postponed to 5.1.

4.1 Reminder on LDPs for Gibbs measures, after [25]↩︎

Equip the set \(\mathcal{P} (\mathbb{R}^d)\) of probability measures on \(\mathbb{R}^{d}\) with the topology of weak convergence. The following theorem is a general large deviations result for Gibbs measures, due to [25].

When Conditions [assum:ldp:stable] to [assum:ldp:regularity] are met, it is said that \(W\) is the macroscopic zero-temperature limit of \(W_n\). Our objective is to apply Theorem 8 to the discrete energy \(H_{n}^{\mathrm{Q}}\), where we recall from 4 that \[\label{e:hnq} H_{n}^{\mathrm{Q}}(y_1, \dots, y_n) = \frac{1}{2n^2}\sum_{i \neq j} K(y_i, y_j) + \frac{1}{n}\sum_{i = 1}^{n} v(y_i) - \frac{1}{n}\sum_{i = 1}^{n} U_{K}^{\nu_n}(y_i),\tag{23}\] with \(v = 0\) inside \(B(0, R)\) and \(v = \Phi\) outside and where \(\nu_n\) is a random measure. The role of the macroscopic zero-temperature limit \(W\) in Theorem 8 would naturally be played by the energy functional in 1, that is \[I_{K}^{V^{\pi}}(\mu) = \frac{1}{2}\iint K(x, y)\,\mathrm{d}\mu(x)\,\mathrm{d}\mu(y) + \int v\,\mathrm{d}\mu - \int U_{K}^{\pi}\,\mathrm{d}\mu = I_{K}^{v}(\mu) - \int U_{K}^{\pi}\,\mathrm{d}\mu\] when \(I_{K}^{V^{\pi}}(\mu) < \infty\).

The main difficulty in applying Theorem 8 is the dependence of \(H_n^{\mathrm{Q}}\) to the random measure \(\nu_n\). [25] does use 8 to prove large deviations for conditioned Gibbs measures (their Theorem 4.1) and particles in varying environment (their Theorem 4.2), but their framework cannot be directly applied to our case since they either restrict to particles supported on a compact subset, or use a different sign convention for the background distributions. Neither of these settings covers our case, since we wish for a Gibbs measure supported on \(\left(\mathbb{R}^{d}\right)^{n}\), and because the conditioned Gibbs measure considered in [25] does not converge to an easily predefinable target due to their sign convention. We will however rely on some standard arguments from [25] in the upcoming proofs. Relating their proof to our framework turned out to be relatively easy for non-singular interactions, but tricky in the Coulomb case, where our 7 will give the missing key result.

We start with a result that only considers the deterministic part of \(H_n^{\mathrm{Q}}\), and can be considered a corollary of 8.

For the specific choice of \(v=0\) in \(B(0,R)\) and \(v= \Phi\) outside, 3 is clearly satisfied, so that we can already apply 9 with \(K\) being either Coulomb or satisfying 2. Defining \(W_n\) and \(\nu\) as in 9, \(I_{K}^{v}\) is the zero-temperature macroscopic limit of \(W_n\), and we can rewrite the unnormalized quenched Gibbs measure we are interested in as \[\begin{gather} \exp\left(-\beta_n H_n^{\mathrm{Q}}(y_1, \dots, y_n)\right) \,\mathrm{d}y_1\,\dots\,\mathrm{d}y_n = \\ \exp\left(- \beta_n\left[ W_n(y_1, \dots, y_n) - \frac{1}{n}\sum_{i = 1}^{n} U_{K}^{\nu_n}(y_i)\right]\right)\,\mathrm{d}\nu(y_1)\,\dots\,\mathrm{d}\nu(y_n). \end{gather}\]

If we denote \(W_n^{\mathrm{Q}}(y_1, \dots, y_n) = W_n(y_1, \dots, y_n) - \frac{1}{n}\sum_{i = 1}^{n} U_{K}^{\nu_n}(y_i)\), which is symmetric and measurable, we want to show that \(I_{K}^{V^{\pi}}\) is the zero-temperature macroscopic limit of \(W_n^{\mathrm{Q}}\) \(\mathbb{P}\)-almost surely. Indeed, since the partition function of \(\mathbb{P}_{n, \beta_n}^{\mathrm{Q}, \pi}\) is well-defined \(\mathbb{P}\)-almost surely for \(n\) large enough (i.e. finite and positive, see 4), we then will be able to apply 8, concluding the proof of 3 and 5 respectively. It remains to check that \((W_n^{\mathrm{Q}}, I_{K}^{V^{\pi}})\) satisfies Condition [assum:ldp:stable] to [assum:ldp:regularity] \(\mathbb{P}\)-almost surely and for \(n\) large enough. This part of the proof will differ in the non-singular case and in the Coulomb case.

4.2 Proof of 3 (non-singular kernels)↩︎

Let \(K\) satisfy 2. We examine the conditions of 8 one by one.

First, \(0 \leq K \leq C\) by assumption, so that \(-U_{K}^{\nu_n} \geq -C\) and \(W_n^{\mathrm{Q}}\) satisfies Condition [assum:ldp:stable] since \(W_n\) satisfies it also.

Second, assume that \(\tilde{W}_n^{\mathrm{Q}}(\mu_j) \leq A\) for some constant \(A\). Then the \((\mu_j)_j\) have to be empirical measures \(\mu_j = \frac{1}{n_j} \sum_{i = 1}^{n} \delta_{y_{i, j}}\). Using again \(U_{K}^{\nu_j} \leq C\), this yields \(W_{n_j}(\mu_j) \leq A + C\). Since \(W_n\) satisfies Condition [assum:ldp:confining], this implies that \((\mu_j)_j\) is relatively compact, and thus \(W_n^{\mathrm{Q}}\) satisfies the same condition.

Third, assume that \(\mu_n \rightarrow \mu\) with \(\mu_n\) empirical measures. Since \(K\) is bounded and continuous and \(\nu_n \rightarrow \pi\) \(\mathbb{P}\)-almost surely, \(\iint K\,\mathrm{d}\nu_n\otimes \mu_n \rightarrow \int U_{K}^{\pi}\,\mathrm{d}\mu\) \(\mathbb{P}\)-almost surely by Fubini. Since \(W_{n}\) satisfies Condition [assum:ldp:lower95limit], \(W_n^{\mathrm{Q}}\) also satisfies it \(\mathbb{P}\)-almost surely.

Fourth, let \(\mu\) be such that \(I_{K}^{V^{\pi}}(\mu) < +\infty\). Then, since \(U_{K}^{\pi}\) is bounded, \(I_{K}^{v}(\mu) < +\infty\). By 9, \(W_n\) satisfies Condition [assum:ldp:regularity], so that there exists a sequence \((\mu_n) \in \mathcal{N}\) such that \(\mu_n \rightarrow \mu\) and \(\overline{\lim}\, I_{K}^{v}(\mu_n) \leq W(\mu)\) with \[\mathcal{N} = \{ \xi \in \mathcal{P}(\mathbb{R}^{d})\,:\, \mathrm{KL}(\xi, \nu) < +\infty,\, \overline{\lim}\,\mathbb{E}_{\xi^{\otimes n}}\left[W_n\right] \leq I_{K}^{v}(\xi)\}.\] Since \(-U_{K}^{\pi}\) is bounded and continuous by dominated convergence, \(\mu \mapsto -\int U_{K}^{\pi}\,\mathrm{d}\mu\) is continuous on \(\mathcal{P}(\mathbb{R}^{d})\), so that \(\overline{\lim} I_{K}^{V^{\pi}}(\mu_n) \leq I_{K}^{V^{\pi}}(\mu)\). Finally, for any \(\xi \in \mathcal{P}(\mathbb{R}^{d})\), \[\mathbb{E}_{\xi^{\otimes n}}\left[-\frac{1}{n}\sum_{i = 1}^{n} U_{K}^{\nu_n}(y_i)\right] = - \int U_{K}^{\nu_n}\,\mathrm{d}\xi = -\int U_{K}^{\xi}\,\mathrm{d}\nu_n\] by Fubini. Again, this yields \(\mathbb{E}_{\xi^{\otimes n}}\left[-\frac{1}{n}\sum_{i = 1}^{n} U_{K}^{\nu_n}(y_i)\right] \rightarrow - \int U_{K}^{\xi}\,\mathrm{d}\pi = - \int U_{K}^{\pi}\,\mathrm{d}\xi\) \(\mathbb{P}\)-almost surely, so that \(\overline{\lim}\,\mathbb{E}_{\xi^{\otimes n}}\left[W_n^{\mathrm{Q}}\right] \leq I_{K}^{V^{\pi}}(\xi)\) for any \(\xi \in \mathcal{N}\), and in particular for the \((\mu_n)_n\). Thus Condition [assum:ldp:regularity] is satisfied for \(W_n^{\mathrm{Q}}\), \(\mathbb{P}\)-almost surely. This concludes the proof of 3.

4.3 Proof of 5 (Coulomb kernel)↩︎

Similarly, when \(K=g\) is the Coulomb kernel and for our specific choice of \(v,\) we check that all the conditions to apply 8 are satisfied.

Since the interaction is neither bounded nor continuous, we loose the continuity and boundedness properties for the energies compared to 4.2, so our approach must differ a bit. We define \[B_n := \underset{z \in \mathbb{R}^{d}}{\sup}\left\lvert U_{K_{\zeta}}^{\nu_n}(z) - U_{g}^{\pi}(z)\right\rvert.\] We show 5 given 7, whose proof is the technical crux and is postponed to 5.1.

We first consider Condition [assum:ldp:lower95limit]. Assume that the sequence of empirical measures \(\mu_n = \frac{1}{n}\sum_{i = 1}^{n} \delta_{z_i}\) converges to \(\mu\). Due to the maximum principle, if \(U_{g}^{\pi} \leq M\) on the support of \(\pi\) then this holds everywhere [42]. Under the assumption that \(\pi\) has compact support and continuous density, \(z \mapsto U_{g}^{\pi}(z)\) is continuous [24], so that \(U_{g}^{\pi}\) is bounded on the support of \(\pi\) and thus everywhere. Thus, \(\mu \longmapsto -\int U_{g}^{\pi}\,\mathrm{d}\mu\) is lower semi-continuous and bounded from below. In particular we have \[\begin{align} \label{e:uniform95bound95pot95ldp} C \leq - \int U_{g}^{\pi}\,\mathrm{d}\mu \leq \underline{\lim}\, \left\{-\frac{1}{n} \sum_{i = 1}^{n} U_{g}^{\pi}(z_i)\right\} \leq \underline{\lim}\,\left\{-\frac{1}{n} \sum_{i = 1}^{n} U_{K_{\zeta}}^{\nu_n}(z_i) + B_n\right\}. \end{align}\tag{24}\] Furthermore, \[\begin{align} \label{e:sum95lim95infs} \underline{\lim}\, W_{n}(z_1, \dots, z_n) + \underline{\lim}\, \left\{-\frac{1}{n}\sum_{i = 1}^{n} U_{K_{\zeta}}^{\nu_n}(z_i)\right\} \leq \underline{\lim}\, W_{n}^{Q}(z_1, \dots, z_n). \end{align}\tag{25}\] Because \(W_n\) satisfies Condition [assum:ldp:lower95limit] by 9, and \(B_n \rightarrow 0\) \(\mathbb{P}\)-almost surely by 7, 24 and 25 imply that \(W_n^{\mathrm{Q}}\) also satisfies it, \(\mathbb{P}\)-almost surely.

We now consider Condition [assum:ldp:stable]. Using again 7 along with \[\begin{align} -\frac{1}{n}\sum_{i = 1}^{n} U_{K_{\zeta}}^{\nu_n}(z_i) \geq - \frac{1}{n} \sum_{i = 1}^{n} U_{g}^{\pi}(z_i) - B_n, \end{align}\] and the fact that \(U_{g}^{\pi}\) is bounded, \(-\frac{1}{n}\sum_{i = 1}^{n} U_{K_{\zeta}}^{\nu_n}(z_i)\) is bounded from below by a constant \(C\) uniformly in \(n\) for \(n\) large enough and \(\mathbb{P}\)-almost surely. Since \(W_n\) is also uniformly lower bounded, this stays true for \(W_n^{\mathrm{Q}}\).

We now consider Condition [assum:ldp:confining]. We proceed as for non-singular interactions: if \(\tilde{W}_{n_j}^{\text{Q}}(\mu_j) \leq A\), then \(W_{n_j}(\mu_j) \leq A + C\) since we just showed that \(-\frac{1}{n}\sum_{i = 1}^{n} U_{K_{\zeta}}^{\nu_n}(z_i)\) is bounded below for \(n\) large enough. So \((\mu_j)\) is relatively compact and \(W_n^{\mathrm{Q}}\) satisfies the assumption \(\mathbb{P}\)-almost surely for \(n\) large enough.

We finish by Condition [assum:ldp:regularity]. We already know that \(U_{g}^{\pi}\) is bounded and continuous so \(\overline{\lim}\,I_{g}^{V^{\pi}}(\mu_n) \leq I_{g}^{V^{\pi}}(\mu)\) by the same argument as for non-singular interactions. Then, for \(\xi \in \mathcal{P}(\mathbb{R}^{d})\), \[\begin{align} \mathbb{E}_{\xi^{\otimes n}}\left[-\frac{1}{n}\sum_{i = 1}^{n} U_{K_{\zeta}}^{\nu_n}(y_i)\right] = - \int U_{K_{\zeta}}^{\nu_n}\,\mathrm{d}\xi. \end{align}\] Using 7, we get \[\overline{\lim}\,\mathbb{E}_{\xi^{\otimes n}}\left[-\frac{1}{n}\sum_{i = 1}^{n} U_{K_{\zeta}}^{\nu_n}(y_i)\right] \leq - \int U_{g}^{\pi}\,\mathrm{d}\xi,\] \(\mathbb{P}\)-almost surely. Again, the inequality for the sum of superior limits is in the good direction and we obtain \(\overline{\lim}\,\mathbb{E}_{\xi^{\otimes n}}\left[W_n^{\mathrm{Q}}\right] \leq I_{g}^{V^{\pi}}(\xi)\) for any \(\xi \in \mathcal{N}\), the set in the Condition [assum:ldp:regularity] satisfied by \(W_n\). This concludes the proof of 5.

5 Proofs: uniform convergence for the Coulomb potential field↩︎

We use here the notations of 2, in particular the definitions of the Coulomb kernel \(g\) from 12 and its regularized version \(K_\zeta\) in 14 . We give the proofs of 6 and 7 along with the statements of the auxiliary lemmas we use. For conciseness, we delay the proofs of all lemmas to 7.

5.1 Proof of 6↩︎

Recall the proof strategy that we sketched after stating 6. We first bound the differences \(U_{g}^{\pi} - U_{K_{\zeta}}^{\pi}\) and \(U_{K_{\zeta}}^{\mu_n} - U_{K_{\zeta}}^{\mu_{n}^{(\epsilon)}}\) uniformly in \(z\), in 1 and 2. Then, we use 3 and 4 to conclude by the concentration properties of the Coulomb gas and with a suitable choice of regularization parameters. These intermediate results are then assembled into the proof of 6.

5.2 Proof of 7↩︎

In this section, we will always work under 1, 3, 4, and 5, even though most of the following lemmas do not require all of them to be true at once.

We begin in 5 by a Taylor expansion similar to 2, but rather dealing with the weights \(W\) from 20 , the motivation of this will become clear afterwards in the proof. Since \(W\) is only assumed to be continuous, we smooth \(W\) by convolving it with a Gaussian.

In the following lemma, we bound the difference of potentials for weighted measures. To do so, we rely on Cauchy–Schwarz to get bounds in term of the \(L_2\) norm of the difference of the characteristic functions of the measures. We then isolate the weights from the \(L_2\) norm of the difference of characteristic functions for unweighted measures, which is the interaction energy between the unweighted measures (see 7). Recall that we know how to bound the latter through concentration inequalities (4).

We will now conclude using the concentration as before and trying to optimize the different parameters to make all the error terms go to \(0\).

Acknowledgements↩︎

We thank Raphaël Butez for stimulating discussions around LDPs and Gibbs measures. We acknowledge support from ERC grant Blackjack ERC-2019-STG-851866, ANR grant Baccarat ANR-20-CHIA-0002 and Labex CEMPI (ANR-11-LABX-0007-01).

6 Experiments↩︎

After discussing how we approximately sample from Gibbs measures, we illustrate our main results on toy examples and a Bayesian classification task.

6.1 Approximately sampling from the Gibbs measure↩︎

There is no known exact sampling algorithm to draw points from Gibbs measures of the form \[\begin{align} \label{xvpniard} \mathrm{d}\mathbb{P}_{n, \beta_n}^{V}(y_1, \dots, y_n) = \frac{1}{Z_{n, \beta_n}^{V}}\exp\left(-\beta_n H_n(y_1, \dots, y_n)\right)\,\mathrm{d}y_1\,\dots\mathrm{d}y_n, \end{align}\tag{34}\] with \[\begin{align} \label{rditmzbg} H_n(y_1, \dots, y_n) = \frac{1}{2 n^2}\sum_{i \neq j} K(y_i, y_j) + \frac{1}{n}\sum_{i = 1}^{n}V(y_i), \end{align}\tag{35}\] and \[\begin{align} \label{esqybrgd} V(z) = V^{\pi}(z) := - U_{K}^{\pi}(z) + \Phi(z), \quad z\in\mathbb{R}^d. \end{align}\tag{36}\] We thus rely on MCMC to produce points \(y_{1}^{T}, \dots, y_{n}^{T}\) that are approximately distributed according to 4 after \(T\) iterations of an MCMC kernel. More precisely, we use Metropolis-Adjusted Langevin (MALA; [2]) updates targeting 4 . Formally, these are Metropolis–Hastings updates with Gaussian proposal given by \[\begin{align} \label{e:mala95updates} y_{1:n} \,|\, y_{1:n}^{(t)} \sim \mathcal{N}\left(y_{1:n}^{(t)}-\alpha \beta_n \nabla H_n(y_{1:n}^{(t)}),\, 2\alpha I_{dn}\right), \end{align}\tag{37}\] where \(y_{1:n}\) is short for \(y_1, \dots, y_n\) and \(\alpha = \alpha_0 \beta_{n}^{-1}\), with \(\alpha_0\) tuned so as to reach an acceptance proportion close to \(50\%\). This choice comes from the fact that we want to avoid mixing times to depend on \(\beta_n\) in warm start cases, following heuristics from [43]; see [20] for more details. Regarding the computational complexity, each of the \(T\) MALA iterations requires evaluating \(H_n\) and its gradient. Each of the \(T\) iterations thus costs \(\mathcal{O}(n^2)\) for the target 4 , and \(\mathcal{O}(n^2 + n M_n)\) for the quenched Gibbs measure of Definition 2.5, where \(M_n\) is the number of atoms of \(\nu_n\) used to approximate \(\pi\) in 6 .

We consider a toy example where \(d=3\) and \(n=500\). The target distribution \(\pi\) is a centered Gaussian distribution with covariance \(\sigma^2 I_d\), truncated at \(\lvert x \rvert > 5 \sigma\), and with \(\sigma = 0.5\). We consider the quenched Gibbs measure with \(\beta_n = n^2\), and interaction kernel the regularized Riesz kernel \(K_{s, \epsilon}\) given by \[\begin{align} \label{bwlsnjkz} K_{s, \epsilon}(x, y) = \frac{1}{\left(\lvert x -y \rvert^2 + \epsilon^2\right)^{s/2}}, \end{align}\tag{38}\] with \(s = d-2\) and \(\epsilon = 0.1\). As for the background, we set \(\nu_n = \frac{1}{M_n}\sum_{i = 1}^{M_n}\delta_{x_i}\) where \(x_1, \dots, x_{M_n}\) is a realization of the history of a Metropolis–Hastings chain targeting \(\pi\), with \(M_n = 1\,000\).

In 1 (a), we show the first two spatial coordinates of the \(50\,000\)-th iterate of a MALA chain approximating the Gibbs measure. For comparison, in 1 (b), we show the first two spatial coordinates of an approximate sample \(y_1, \dots, y_n\) from \(\pi\), obtained as the history \(y_1, \dots, y_n\) of a Metropolis–Hastings chain targeting \(\pi\) with Gaussian proposal \(N(0, \alpha I_d)\), where \(\alpha\) is tuned to reach \(50\%\) acceptance, and we removed \(5\,000\) burn-in iterations.

We see in 1 that both methods stay confined in the regions of large probability under the target \(\pi\), as expected. Moreover, the approximated Gibbs samples in 1 (a) are more diverse than a Metropolis-Hastings sample, let alone because about half the proposed points in the latter are rejected, explaining the many superimposed disks in 1 (b). On top of these rejections causing autocorrelation in the MCMC chain, we expect the interaction term in 4 to favor diverse configurations. However, this promotion of diversity is balanced by 1) the use of an MCMC chain as a confining background, 2) the fact that we truncate the interaction kernel, 3) the constraint that a significant proportion of the points should remain in the bulk of the Gaussian target, as well as 4) the fact that we project a 3-dimensional sample onto a plane. This explains why repulsiveness is less visually striking in Figure 1 (a) than in 1 (c), where we plot for reference an approximate Gibbs sample when the target \(\pi\) is the 2D uniform measure on \(B(0, 1)\), with interaction kernel \(-\frac{1}{2}\log(\lvert x - y\rvert^2)\), exact potential \(U^{\pi}\) and \(T = 5\,000\) MALA iterations.

In 2, we plot the first two coordinates in the same setup as 1 (a), with the same realization of the background, but this time with an inverse temperature of \(\beta_n = n^3\). In 2 (a) and 2 (b), we see that points spread a lot in the domain of the target, even in regions of low probability under \(\pi\). In light of 2 (c), where outer regions of low probability under \(\pi\) now contain a matching fraction of the points, one can understand the overspread in Figures 2 (a) and 2 (b) as a consequence of MALA needing more iterations to yield an accurate approximation to its Gibbs target. This supports the intuitive claim that higher \(\beta_n\) comes with a worse mixing time for MALA: the lower the temperature, the more difficult it is to approximately sample a Gibbs measure. There thus seems to be a tradeoff in choosing a large \(\beta_n\) for better statistical guarantees, and keeping \(\beta_n\) low to preserve tractable approximate sampling procedures. There are to our knowledge no theoretical results that enable to quantify this tradeoff, and one could hope that Markov kernels specifically designed for 4 could yield better performance, e.g. smaller confidence intervals for linear statistics. While those are important questions, we do not address them in this paper. In practice, we recommend \(\beta_n = n^2\), which already yields improved statistical guarantees over independent and MCMC quadratures, and has remained computationally tractable in our experiments.

6.2 Worst-case integration error and variance↩︎

We consider the same setting where \(\pi\) is a three-dimensional Gaussian distribution \(\mathcal{N}(0, \sigma^2 I_{3})\), truncated at \(\lvert x \rvert > 5 \sigma\) with \(\sigma = 0.5\). We estimate the worst-case integration error of various quadratures \(\mu_n := \frac{1}{n}\sum_{i = 1}^{n}\delta_{y_i}\) as \(n\) grows.

We first evaluate \(\mathrm{E}_{\mathcal{H}_K}(\mu_n) = \sqrt{I_{K}(\mu_n - \pi)}\) when \(K = K_{s, \epsilon}\) again with \(s = d-2\) and \(\epsilon = 0.1\). In particular, we compute \[I_{K}(\mu_n - \pi) - I_{K}(\pi) = I_{K}(\mu_n)-2I_{K}(\mu_n, \pi),\] using a long MH chain targeting \(\pi\) to approximate the integral with respect to \(\pi\), with size \(M = 90,000\) and \(10,000\) burn-in iterations. We show in 3 (a) how \(\mathrm{E}_{\mathcal{H}_K}(\mu_n)\) decays when \(n\) grows for i) \(\mu_{n}^{\mathrm{MCMC}}\), where the \((y_i)\) is the history of a MH chain targeting \(\pi\) after \(5000\) burn in iterations, ii) \(\mu_{n}^{\mathrm{Q}}\), where the \((y_i)\) are an approximate sample from the quenched Gibbs measure with \(T = 10\,000\) iterations, temperature \(\beta_n = n^2\), kernel \(K = K_{s, \epsilon}\) and \(\nu_n = \frac{1}{M_n}\sum_{i = 1}^{M_n} \delta_{x_i}\) where \(M_n = 1,000\) and the \(x_i\) are the history of a MH chain targeting \(\pi\) again, iii) \(\mu_{n}^{\mathrm{KT}}\) where the \((y_i)\) are a coreset of size \(n\) of an \(n^2\) MCMC chain \(x_1, \dots x_{n^2}\) targeting \(\pi\). In the latter method, the coreset if obtained via the kernel thinning algorithm of [13], [44], who proved that the \((y_i)_{i = 1, \dots, n}\) preserve the same sample quality as the original \((x_i)_{i = 1, \dots, n^2}\) with high probability. We compute the \(90\%\)-quantile of the worst-case error, over \(100\) independent realizations of the quadrature for each method and each cardinality \(n\), with fixed background for ii) and iii).

As expected, we see that both kernel thinning and our method outperform MCMC in terms of worst-case integration error, at least up to \(n = 350\) particles. For kernel thinning, we used the vanilla version presented in [13] which requires to store all the \(n^4\) pairwise interactions for the underlying MCMC background of size \(n^2\), which is why we could not run the experiment for more than \(n = 350\) due to storage issues. Note however that recent refinements of this algorithm might help circumvent this issue [45].

Beyond \(n = 350\), we observe that the energy increases after for the quenched Gibbs measure, becoming similar to the MCMC error, which we attribute again to the fact that \(T = 10\,000\) MALA iterations are not enough to approximate the Gibbs measure for that large a number of particles.

We now provide further experiments that investigate the impact of the background, as well as a conjectured faster decaying variance for linear statistics.

In 4, we run the same experiment as in 6.2, where we vary the size \(M_n\) of the MCMC background \((x_i)\) to see whether there is an impact on the worst-case integration error. We also add to the comparison a background \(\nu_n\) obtained via reweighting the realization of a Coulomb gas \(x_1, \dots, x_{n}\) as in Definition (2.8), with quadratic potential \(V\) so that \(\mu_V\) is the uniform measure on \(B(0, 5\sigma)\). We use \(T = 10\,000\) iterations of MALA to approximately sample that Coulomb background.

When the background is drawn using MCMC, we see a slight improvement between a background of size \(n\) and a background of (larger) fixed size \(1000\) or \(5000\). When the background is drawn from a first well spread Coulomb gas with reweighting, we observe a small loss in terms of energy, especially for small values of \(n\). We conjecture that this is due to the quality of the estimation when using weights to approximate the potential with a small number of points \(n\).

Then, we look at the effect of increasing the temperature schedule \(\beta_n\) in 5 (a) for MCMC background with size \(M_n = 1000\), and in 5 (b) for a reweighted Coulomb background \(x_1, \dots, x_{n}\) with uniform equilibrium measure. Again, we clearly see a downgrade in terms of energy for lower temperature schedules which is again due to the fact that \(T = 10\,000\) are not enough to sample from the Gibbs measure in those regimes as we saw in 2. Running the same experiments with \(T = 200\,000\) MALA iterations for instance is very costly, and one could argue that it is not reasonable for a three-dimensional truncated Gaussian target.

We now look at the variance of linear statistics with respect to \(\mu_n\) for the three different methods (i-iii) listed above. We conjecture that Gibbs measures can significantly improve on MCMC in that criterion, although we defer the corresponding theoretical investigation to further work. We consider a simple test function in \(\mathcal{H}_{K_{s, \epsilon}}\) given by \(f(x) = \sum_{i = 1}^{10} a_i K_{s, \epsilon}(x, z_i)\) where the \((z_i)\) are chosen uniformly at random on \([-1, 1]^{d}\) and the \(a_i\) are chosen uniformly at random on \([-1, 1]\). We compute in 6 the empirical estimate of the variance \(\mathrm{Var}\left[\frac{1}{n}\sum_{i = 1}^{n} f(y_i)\right]\) over \(100\) independent realizations of \(y_1, \dots, y_n\) for each method and each value of \(n\), with MCMC background of size \(M_n = 1\,000\) and \(T = 10\,000\) MALA iterations for the quenched Gibbs measure.

While kernel thinning clearly exhibits a faster rate of convergence than MCMC, and the variance of the Gibbs estimator can –when tuned well– fall between that of MCMC and kernel thinning, it is hard to conclude anything for the rates corresponding to Gibbs measures. This is again related to the difficulty of tuning the MALA sampler as the Gibbs target changes.

6.3 Bayesian classification↩︎

In this section, we give a practical statistical illustration that a quenched Gibbs quadrature provides better coverage of simultaneous confidence intervals for a given computational budget, compared to classical MCMC.

First, we define the target integration task starting from the foundational expected utility principle; this is classical in Bayesian statistics, we include it for completeness. Consider a training data set \((z_j, t_j)_{j = 1}^{N}\) where \(z_j \in \mathcal{Z} = \mathbb{R}^{d}\) encodes the features of the j-th individual and \(t_j \in \mathcal{T} = \{1, \dots, C\}\) is the associated class label. An action is the choice of a classifier \(g : \mathcal{Z}^{N} \times \mathcal{T}^{N} \times \mathcal{Z} \rightarrow \mathcal{T}\). For a new point \((z, t)\) with unknown class \(t\), if one considers a joint probability measure \(p\) on \(z_1, t_1, \dots, z_N, t_N, z, t\) and the so-called 0 – 1 utility function, the expected utility of the action \(g\) is given by \[\label{e:seu} -\int 1_{t\neq g(z_{1:N}, t_{1:N}, z)} \mathrm{d} p(z_{1:N}, t_{1:N}, z, t),\tag{39}\] where \(z_{1:N}\) is short for \((z_1, \dots, z_N)\). The Bayes action or classifier is the function \(g^\star\) which maximizes the expected utility 39 . It is easy to see that it is given by \[\label{e:bayes95action} g^\star: z_{1:N}, t_{1:N}, z \mapsto \arg\max_{k=1}^C p(t=k \vert z_{1:N}, t_{1:N}, z).\tag{40}\] In other words, the Bayes classifier predicts the class \(k\) with the highest conditional probability under \(p\). Making the classical assumption that \(p\) factorizes as \[\begin{align} \label{e:factorization} p(z_{1:N}, t_{1:N}, z, t) = \int p(y) \cdot \prod_{j=1}^N p(t_j \vert z_j, y) p(z_j) \cdot p(t\vert z, y) p(z) \mathrm{d}y, \end{align}\tag{41}\] where we introduced a parameter \(y \in\mathcal{Y}\) and a reference measure \(\mathrm{d}y\), the conditional in 40 writes \[\begin{align} \label{e:posterior95predictive} p(t = k |z_{1:N}, t_{1:N}, z) = \int p(t = k | z, y) \pi(y) \mathrm{d}y, \end{align}\tag{42}\] where \[\pi(y) \propto \prod_{j = 1}^{N} p(t_j | z_j, y) p(y)\] is the density of the posterior probability measure on the parameter \(y\). Using a Monte Carlo method to replace \(\pi\) by an approximation \(\mu_n = \frac{1}{n} \sum_{i = 1}^{n} \delta_{y_i}\), we approximate the Bayes action by \[\label{e:mc95action} g^{\mathrm{MC}}_n: z_{1:N}, t_{1:N}, z \mapsto \arg\max_{k=1}^C \int p(t=k \vert z, y) \mathrm{d}\mu_n(y).\tag{43}\]

One way to guarantee that the expected utility 39 of \(g^{\mathrm{MC}}_n\) is close to the one of \(g^{\star}\) is to ensure that the \(C\) integrals in the argmax of 43 simultaneously approximate the \(C\) integrals in the argmax of 40 . Furthermore, if the dataset on which we want to perform classification contains \(M\) test points instead of only one and we want to ensure that we approximate all predictions of the Bayes classifier using a single Monte Carlo run \(y_1, \dots, y_n\), then there are \(M\times C\) integrals to simultaneously approximate. In that case, Theorem (2.7) and Theorem (2.10) show that the asymptotic confidence intervals obtained with a single run \(y_{1:n}^{\mathrm{Q}}\) from the quenched Gibbs measure are tighter than the ones obtained for a single run \(y_{1:n}^{\mathrm{MCMC}}\) from MCMC, depending on which smoothness one wants to assume on \(p(t = k |z, y)\). We empirically illustrate in the following that it is also the case at finite \(n\) and for a MCMC-approximated sample of the quenched Gibbs measure, using proposal 37 .

To specify our setting, we consider binary logistic regression, where \(C = 2\), \(\mathcal{Z} = \mathbb{R}^{2}\), \(\mathcal{T} = \{0, 1\}\), \(d =3\), and we assume \(p\) to factorize as in 41 with \[\label{e:logistic95regression} p(t = 1 | z, y) = 1 - p(t = 0 | z, y) \propto \exp\left(y^{\top}[z, 1]\right),\tag{44}\] where \([z, 1]\) stands for \(x\) concatenated with \(1\) to include an intercept. To keep the posterior \(\pi\) compactly-supported, we consider the prior \(p(y)\) to be a three-dimensional Gaussian distribution with covariance \(\sigma^2 I_3\) truncated at \(\lvert x \rvert > 5 \sigma\), with \(\sigma = 0.5\) as before, so that \(\pi\) is compactly supported. We mimick a toy dataset from the documention of the blackjax⁴ library: the features \(z_{1:N}\) are drawn from a balanced mixture of two isotropic Gaussian distributions, and the class \(t_{1:N}\) corresponds to the component membership. We consider \(N = 50\) individuals for the training set, and \(M = 10\) test individuals to classify. Computing the Bayes classifier for these \(10\) individuals then requires to compute \(M (C-1) = 10\) integrals with respect to \(\pi\). In order to be able to compare with a reference value for the Bayes classifier, we first approximate the conditional \(p(t = 1 | z_{1:N}, t_{1:N}, z)\) in 42 , for each \(z\) in the test set, using a long Metropolis–Hastings chain targeting \(\pi\), with \(10^5\) iterations. Then, we compare the results to the Monte Carlo approximation \[\label{e:bayes95action95mc} \hat{p}_{n}(t = 1 | z_{1:N}, t_{1:N}, z) = \frac{1}{n} \sum_{i = 1}^{n} p(t = 1 | z, y_i),\tag{45}\] when \(y_1, \dots, y_n\) is i) the history of a Metropolis–Hastings chain targeting \(\pi\) after \(500\) burn-in iterations, ii) an approximate sample from the quenched Gibbs measure, after \(T = 10\,000\) MALA iterations, with regularized Riesz kernel with \(s = d -2\) and \(\epsilon = 0.1\), at inverse temperature \(\beta_n = n^2\), and with background \(\nu_n\) of size \(M_n = 1000\) given by an iid subsample of the history of a Metropolis–Hastings chain. We build \(50\) independent samples for each method, and compute for each value of \(\delta\) the proportions \(\mathfrak{p}_{\delta, \mathrm{MCMC}}\) and \(\mathfrak{p}_{\delta, \mathrm{Gibbs}}\) of realizations of samples for which the \(M\) estimators 45 simultaneously fall in the \(M\) confidence regions \[[\hat{p}_{n}(t = 1 | z_{1:N}, t_{1:N}, z) - \delta\,;\, \hat{p}_{n}(t = 1 | z_{1:N}, t_{1:N}, z) + \delta]\] corresponding to the \(M\) values \(z\) in the test set. We plot \(\mathfrak{p}_{\delta, \mathrm{MCMC}}\) and \(\mathfrak{p}_{\delta, \mathrm{Gibbs}}\) for different values of \(\delta\) in 3 (b) when \(n = 100\) for ii) and for different values of \(n\) for i). We use Gaussian confidence intervals at confidence level \(95\%\) with Bonferroni correction for those \(5\) values of \(\delta\).

As expected, the Gibbs measure yields tighter simultaneous confidence intervals. In particular, \(100\) particles drawn from the quenched Gibbs measure exhibit a coverage quality that is at least the same the history of a Metropolis–Hastings chain of size \(10\,000\).

7 Proof of auxiliary lemmas↩︎

7.1 Preliminaries↩︎

We will repeatedly use a particular integral representation for \(g\) and \(K_{\zeta}\). First, a quick computation shows that for any \(c > 0\) and \(\alpha > 0\), \[\Gamma(\alpha) := \int_{0}^{+\infty} t^{\alpha -1}e^{-t}\,\mathrm{d}t = c^{\alpha} \int_{0}^{+\infty} t^{\alpha -1}e^{-ct}\,\mathrm{d}t.\] choosing \(\alpha = \frac{d-2}{2}\), this yields \[\begin{align} g(z, y) &= \frac{1}{\Gamma\left(\frac{d-2}{2}\right)}\int_{0}^{+\infty} t^{\frac{d-2}{2}-1}\exp\left(- t \lvert z - y\rvert^2\right)\,\mathrm{d}t,\tag{46}\\ K_{\zeta}(z, y) &= \frac{1}{\Gamma\left(\frac{d-2}{2}\right)}\int_{0}^{+\infty}t^{\frac{d-2}{2}-1}\exp\left(-n^{-2\zeta} t\right)\exp\left(- t \lvert z - y\rvert^2\right)\,\mathrm{d}t.\tag{47} \end{align}\]

It will be often convenient to work in Fourier space to bound the energies. In fact [4] states that for bounded translation invariant kernels \(K(x, y) = \Psi(\lVert x - y\rVert)\), the energy \(I_g\) of a signed measure can be seen as the \(L^{2}\)-norm of the characteristic function of that measure, w.r.t. to a spectral measure which \(\Psi\) is the Fourier transform of. We give explicit expressions for this spectral measure in the following lemma for \(K_{\zeta}\) and \(g\), which proves the formula at the same time even though \(g\) is not bounded.

We give estimates on \(\Lambda_{\zeta}\) in the next lemma, which we will use in the proof of 7.

We will also use the following boundedness property of Riesz potentials.

7.2 Auxiliary lemmas for 6↩︎

7.3 Auxiliary proofs for 7↩︎

References↩︎

---

1. Corresponding author: rouault.martin@protonmail.com↩︎

2. e.g. by Cauchy-Schwarz in \(\mathcal{H}_{K}\).↩︎

3. When \(K=g\), it has been called a Coulomb gas, one-component plasma, or jellium.↩︎

4. https://blackjax-devs.github.io/sampling-book/models/logistic_regression.html↩︎