Martin Rouault
email: martin.rouault@univ-lille.fr

Rémi Bardenet

Mylène Maïda

Monte Carlo with kernel-based Gibbs measures: Guarantees for probabilistic herding

Kernel herding belongs to a family of deterministic quadratures that seek to minimize the worst-case integration error over a reproducing kernel Hilbert space (RKHS). In spite of strong experimental support, it has revealed difficult to prove that this worst-case error decreases at a faster rate than the standard square root of the number of quadrature nodes, at least in the usual case where the RKHS is infinite-dimensional. In this theoretical paper, we study a joint probability distribution over quadrature nodes, whose support tends to minimize the same worst-case error as kernel herding. We prove that it does outperform i.i.d. Monte Carlo, in the sense of coming with a tighter concentration inequality on the worst-case integration error. While not improving the rate yet, this demonstrates that the mathematical tools of the study of Gibbs measures can help understand to what extent kernel herding and its variants improve on computationally cheaper methods. Moreover, we provide early experimental evidence that a faster rate of convergence, though not worst-case, is likely.

Introduction

Numerical integration with respect to a possibly unnormalized target distribution π on R d has become routine in computational statistics [START_REF] Robert | The Bayesian choice: from decision-theoretic foundations to computational implementation[END_REF] and probabilistic machine learning [START_REF] Murphy | Probabilistic machine learning: Advanced topics[END_REF]. Monte Carlo algorithms [START_REF] Robert | Monte Carlo statistical methods[END_REF] are randomized algorithms that tackle this task, defining estimators that rely on n evaluations of the target integrand at suitably chosen random points in R d , called nodes. Classical Monte Carlo algorithms, such as Markov Chain Monte Carlo (MCMC), come with a set of probabilistic error controls, such as central limit theorems, that involve errors of magnitude n -1/2 . The popularity of MCMC in practice has justified a continuous research effort to improve on that rate, which is considered too slow when evaluating the integrand is computationally costly. Quasi-Monte Carlo methods (QMC), for instance, rely on smoothness assumptions to obtain worst-case error controls of order 1/n. A common such smoothness assumption is that the target integrands belong to the unit ball of a particular reproducing kernel Hilbert space (RKHS); see e.g. (Dick et al., 2013, Section 3).

At another end of the algorithmic spectrum, variational Bayesian methods (VB; [START_REF] Blei | Variational inference: A review for statisticians[END_REF][START_REF] Liu | Stein variational gradient descent: A general purpose bayesian inference algorithm[END_REF] sacrifice some of the error controls to gain in scalability. At its core, VB is the minimization of a dissimilarity measure between a candidate approximation and the target distribution π. Minimizing a relative entropy, for instance, yields algorithms amenable to stochastic gradient techniques [START_REF] Hoffman | Stochastic variational inference[END_REF], yet that usually come with no theoretical guarantee on how well integrals w.r.t. π are approximated.

An intermediate method between Monte Carlo and relative entropy-based VB is the minimization of an integral probability metric (IPM) of the form

ν → I K (ν -π) = K(x, y) d(ν -π) ⊗2 (x, y), (1)
where K is a positive definite kernel, known as the interaction kernel. It is known [START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF] that the square root of I K (ν -π) in ( 1) is the worst-case integration error for integrands in the unit ball of the RKHS defined by K, when approximating π by ν; see [START_REF] Pronzato | Bayesian quadrature, energy minimization, and space-filling design[END_REF] for a recent survey. Loosely speaking, minimizing (1) is thus an attempt at designing efficient algorithms in the vein of VB, yet that come with a control on the integration error like Monte Carlo and QMC. Kernel herding, for instance, which rose to attention in the context of learning Markov random fields [START_REF] Welling | Herding Dynamical Weights to Learn[END_REF][START_REF] Welling | Herding dynamic weights for partially observed random field models[END_REF]Chen & Welling, 2010), has been shown to actually be a conditional gradient descent that greedily minimizes (1) (Chen et al., 2010;[START_REF] Bach | On the Equivalence between Herding and Conditional Gradient Algorithms[END_REF]. One practical limitation of kernel herding is the requirement to evaluate the kernel embedding K(•, x) dπ(x). When the support of the target measure π is all of R d , this requirement can be circumvented by a specific choice of K, namely the Stein kernel [START_REF] Anastasiou | Stein's Method Meets Computational Statistics: A Review of Some Recent Developments[END_REF]. In that case, the IPM (1) coincides with the kernel Stein discrepancy (KSD), and simple gradient descent schemes can yield efficient minimization algorithms [START_REF] Korba | Kernel stein discrepancy descent[END_REF]. To our knowledge, besides the need to evaluate the kernel embedding, the main theoretical limitation of kernel herding and its variants is that, while there is experimental support in favor of an n -1 convergence rate of the worst-case integration error in the RKHS induced by K (Chen et al., 2010;[START_REF] Pronzato | Performance analysis of greedy algorithms for minimising a Maximum Mean Discrepancy[END_REF], there is no result that shows an improvement over the Monte Carlo rate n -1/2 when the RKHS is infinite-dimensional [START_REF] Bach | On the Equivalence between Herding and Conditional Gradient Algorithms[END_REF], see Section 2.2. Our paper provides a step in this direction, for a randomized relaxation 1 of herding.

Gibbs measures are probability distributions that describe systems of interacting particles. By choosing the interaction carefully, one can arrange the corresponding Gibbs measure to favor configurations of points that tend to minimize I K in (1), when a suitable inverse temperature parameter goes to infinity. Gibbs measures have been studied for decades in probability and mathematical physics, with a focus on models that relate to electromagnetism [START_REF] Serfaty | Systems with Coulomb interaction[END_REF]. Classical results include large deviation principles (LDPs; [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF] and concentration inequalities in some cases (Chafaï et al., 2018). Our main contribution is to prove that the Gibbs measure whose points repel each other by an amount given by a bounded kernel K(x, y) satisfies a concentration inequality for the worst-case integration error (1); see Theorem 3.5. Our Corollary 3.6 shows a faster sub-Gaussian decay than the i.i.d or MCMC case, coming from the kernel-dependent repulsion. In other words, our probabilistic relaxation of herding provably outperforms classical Monte Carlo methods, in the limited sense that it requires fewer nodes to reach a given worst-case integration error in any given RKHS. We believe this is important, in the sense that this is a first step in establishing the faster convergence of IPM minimization algorithms. There are many limitations to be studied in future work, however. In particular, we have no exact sampling algorithm for our Gibbs measure, which forces us to use MCMC in our experiments. Moreover, we assume that the target measure has compact support, which prevents using the Stein kernel trick to avoid evaluating the kernel embedding.

The rest of the paper is organized as follows. In Section 2, we quickly survey worst-case controls on the integration error. In Section 3, we introduce a family of Gibbs measures and state our theoretical results. In Section 4, we explain how to approximately sample from such Gibbs measures, and experimentally validate our claims. We discuss perspectives in Section 5. Proofs are deferred to the appendix.

Related work

We survey worst-case guarantees for Monte Carlo and IPM minimization algorithms.

Uniform concentration for Monte Carlo

The introduction of a new Monte Carlo method is typically backed up by a central limit theorem [START_REF] Robert | Monte Carlo statistical methods[END_REF]. In practice, where the number n of quadrature nodes is fixed, one prefers a concentration inequality, to derive a confidence interval for f dπ. While rarely put forward, many applications further require a uniform control over several integrands at a time. For instance, in multi-class classification with 0/1 loss and M classes, determining the Bayesian predictor involves giving a joint confidence region over M -1 integrals. This motivates studying the simultaneous approximation of several integrals by a single set of n Monte Carlo nodes. One way to formalize this problem is by upper bounding the Wasserstein distance

W 1 (µ n , π) = sup ∥f ∥Lip≤1 f d(µ n -π) , (2)
where µ n is the Monte Carlo empirical approximation of the target measure π, and the supremum is taken over all Lipschitz functions of Lipschitz constant less than 1. Sanov's theorem (Dembo & Zeitouni, 2009, Theorem 6.2.10) gives a large deviation principle (LDP), i.e. an asymptotic control on the tails of the random variable (2) when µ n is the empirical measure of i.i.d. draws from π. Non-asymptotic counterparts have been obtained through sub-Gaussian concentration inequalities with speed n, see [START_REF] Bolley | Quantitative concentration inequalities for empirical measures on non-compact spaces[END_REF][START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure[END_REF] for more details.

Theorem 2.1 [START_REF] Bolley | Quantitative concentration inequalities for empirical measures on non-compact spaces[END_REF]. Assume that x 1 , . . . , x n are drawn i.i.d from π. Under a suitable moment assumption on π, for any d ′ > d, there exists n 0 such that, for any n ≥ n 0 and

r > n -1/(d ′ +2) , P [W 1 (µ n , π) > r] ≤ exp -αnr 2 , ( 3
)
where α is a constant depending on π.

Note the regime restriction r > n -1/(d ′ +2) . Analogous LDPs and concentration bounds with the same rate exist for the empirical measure of Markov chains, see (Dembo & Zeitouni, 2009, Chapter 6.5) and [START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure[END_REF]. From a quadrature point of view, results like Theorem 2.1 give simultaneous confidence intervals.

Corollary 2.2. Fix ϵ > 0, and δ ∈ (0, 1). For any n big enough so that n -1/(d ′ +2) ≤ ϵ and exp -αnϵ 2 ≤ δ, with probability higher than 1 -δ,

sup ∥f ∥Lip≤1 f d(µ n -π) ≤ ϵ. ( 4
)
where x 1 , . . . , x n ∼ π are i.i.d. and µ n = 1 n n i=1 δ xi . If we take the supremum over smoother functions, like the functions in the RKHS H K [START_REF] Berlinet | Reproducing kernel Hilbert spaces in probability and statistics[END_REF] of positive definite kernel K, one can further hope to reach inequality (4) for smaller values of n. In that case, i.i.d. sampling from the so-called kernel leverage score yields a concentration like (4) with a faster rate, though not fully explicit (Bach, 2017, Proposition 2). [START_REF] Belhadji | Kernel interpolation with continuous volume sampling[END_REF] give a more explicit rate, but depart from the i.i.d. setting by studying a kernel-dependent joint distribution on the quadrature nodes called volume sampling. In short, under generic assumptions on the kernel, Markov's inequality applied to (Belhadji et al., 2020, Theorem 4) shows that there is a constant C > 0 such that under volume sampling, with probability 1 -Cσ n+1 /ϵ,

sup f ∈H K f dπ - n i=1 w i f (x i ) 2 ≤ ϵ, (5)
where σ n is the n-th eigenvalue of the operator on L 2 (π) with kernel K, and the weights (w i ) are suitably chosen. Because σ n can go to zero arbitrarily fast with n (e.g., exponentially for the Gaussian kernel), ( 5) attains a given confidence and error levels at smaller n than under i.i.d. sampling. Downsides are that (i) there is no exact algorithm yet for volume sampling that does not require to evaluate the eigenvalues and eigenfunctions of the integral operator with kernel K, and (ii) the dependence of w i on all nodes makes it hard to derive, e.g., a central limit theorem.

Variational Bayes and Kernel herding

Convergence guarantees for VB are often formulated in terms of the minimized dissimilarity measure; see [START_REF] Alquier | On the properties of variational approximations of Gibbs posteriors[END_REF][START_REF] Lambert | Variational inference via Wasserstein gradient flows[END_REF] and references therein. For instance, under strong assumptions on the target π and the allowed variational approximation, [START_REF] Lambert | Variational inference via Wasserstein gradient flows[END_REF] have given rates for the convergence to the minimal achievable relative entropy KL(ρ t ||π) between π and the t-th iterate ρ t of an idealized (continuous-time) VB algorithm. For Stein variational gradient descent, [START_REF] Liu | Stein variational gradient descent: A general purpose bayesian inference algorithm[END_REF] and [START_REF] Korba | A non-asymptotic analysis for Stein variational gradient descent[END_REF] proved a decay of KL(•||π) along with non-asymptotic bounds at rate n -1 for the kernel Stein discrepancy (KSD) between an (at most) n-point empirical measure based on the algorithm and the target measure π. Since the KSD is a particular case of (1), this implies a control on a worst-case integration error. Indeed, if K is positive definite, so that it defines an RKHS H K , it can be shown that 1) is the square of the worst-case integration error over the unit ball of H K when replacing π by ν; see e.g. [START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF]. For KSD, the RKHS corresponds to the Stein kernel. Under some assumptions on H K , kernel herding algorithms have been proved to achieve I K (µ nπ) ≤ cn -2 , so that the worst-case quadrature error decreases at rate n -1 (Chen et al., 2010). Other variants of conditional gradient algorithms led to further improvement, up to convergences of the type [START_REF] Bach | On the Equivalence between Herding and Conditional Gradient Algorithms[END_REF]. While those assumptions are reasonable when the dimension of H K is finite, [START_REF] Bach | On the Equivalence between Herding and Conditional Gradient Algorithms[END_REF] have shown that they are never fulfilled in the infinite-dimensional setting. In that case, the only general result is the "slow" rate I K (µ n -π) ≤ cn -1 [START_REF] Bach | On the Equivalence between Herding and Conditional Gradient Algorithms[END_REF]. Proving a faster rate for a variant of herding remains an open problem.

I K (ν -π) in (
I K (µ n -π) ≤ exp (-cn)

LDPs and concentration for Gibbs measures

We informally define a (Gibbs) measure on R d n by

P V n,βn (dx 1 ...dx n ) ∝ e -βnHn(x1,...,xn) dx 1 ...dx n , (6)
where β n > 0 is called inverse temperature, and where

H n (x 1 , ..., x n ) = 1 2n 2 i̸ =j K(x i , x j ) + 1 n n i=1 V (x i ), (7)
with V : R d → R. There are assumptions to be made on K and V to guarantee that (6) defines a bona fide probability distribution; see Section 3. H n in ( 7) can be recognized to be a discrete analogous to I K in (1).

Intuitively, points distributed according to (6) tend to correspond to low pairwise kernel values K(x i , x j ) (we say that they repel by a force given by the kernel), yet stay confined in regions where V is not too large. We also emphasize that the zero temperature limit -informally taking β n = +∞ in (6) -corresponds to finding the deterministic minimizers of H n , which in turn intuitively correspond to minimizer of I K ; see [START_REF] Serfaty | Systems with Coulomb interaction[END_REF] for precise results. We focus in this paper on the so-called low-temperature regime β n /n → +∞ (denoted in the sequel by β n ≫ n), in which one can hope to observe properties of the Gibbs measure (6) that depart from those of i.i.d. sets of n points. Asymptotic properties of (6) as n → ∞ have been studied by [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF], who prove an LDP like Sanov's classical result, but with the rate 1/n replaced by 1/β n , which goes faster to 0 in the regime β n ≫ n. Moreover, the convergence of the empirical measure is now towards the so-called equilibrium measure µ V , which depends in a non-trivial way on V and K. As for non-asymptotic counterparts to this LDP, concentration inequalities have been obtained for some singular2 kernels, known as the Coulomb and Riesz kernels [START_REF] Chafaï | Concentration for Coulomb gases and Coulomb transport inequalities[END_REF][START_REF] García-Zelada | Generalized transport inequalities and concentration bounds for Riesz-type gases[END_REF]. For instance, Chafaï et al. (2018) prove for the Coulomb kernel that whenever r > n -1/d ,

P V n,βn (W 1 (µ n , µ V ) > r) ≤ exp -cβ n r 2 . ( 8
)
The concentration result (8) improves on the i.i.d. concentration in (3). Besides being valid for values of r down to n -1/d , the speed n in the exponential is replaced by β n , which can increase arbitrarily fast, at the price of replacing the target measure by the equilibrium measure of the system3 . After choosing a suitable potential V such that µ V = π, we rephrase this bound as a uniform quadrature guarantee.

Corollary 2.3 (Chafaï et al., 2018). Fix ϵ > 0, and δ ∈ (0, 1). For any n big enough so that n -1/d ≤ ϵ and exp -cβ n ϵ 2 ≤ δ, with probability higher than 1 -δ,

sup ∥f ∥Lip≤1 f d(µ n -π) ≤ ϵ, (9)
where x 1 , . . . , x n ∼ P V n,βn in (6) and µ n = 1 n n i=1 δ xi . As long as β n ≫ n, for a fixed confidence level δ and uniform worst-case error ϵ, as soon as n is big enough so that n -1/d ≤ ϵ, the constraints in Corollary 2.3 are achieved with a smaller number n of points than for i.i.d. samples in Corollary 2.2. Fewer quadrature nodes are required by the Gibbs measure to achieve the same guarantee.

Results like (Chafaï et al., 2018) are motivated by statistical physics and focus on a particular family of singular kernels. The price of singularity is quite long and technical proofs. On the other hand, in machine learning, we typically consider bounded kernels like the Gaussian or Matern kernel. Our main result is a version of (8) that is valid for very general bounded kernels, bringing an improvement over i.i.d. sampling similar to Corollary 2.3, with n -1/d even replaced by n -1/2 . Maybe surprisingly, while our proof follows the lines of (Chafaï et al., 2018), we were able to considerably simplify the more technical arguments. We hope our work helps transfer tools and concepts from the theory of Gibbs measures to the study of IPM-based quadrature.

As a final note on existing work, and to prepare for the discussion in Section 5, we remark that central limit theorems for Gibbs measures like (6) (with speed depending on β n ) are very subtle mathematical results. [START_REF] Leblé | Fluctuations of two dimensional Coulomb gases[END_REF] and [START_REF] Bauerschmidt | The two-dimensional Coulomb plasma: quasi-free approximation and central limit theorem[END_REF] have obtained a CLT only in dimension two so far, and yet only for the (singular) Coulomb kernel. While important steps have been made towards larger dimensions for the Coulomb kernel (Serfaty, 2023), this remains an important and difficult open problem in statistical physics. As a consequence, direct comparison with the n -1/2 rate appearing in the CLTs of MCMC chains is currently out of reach.

Main results

We first rigorously introduce some key notions to understand the limiting behavior of Gibbs measures, like the equilibrium measure. We then introduce our Gibbs measure on quadrature nodes, and state our main result, which features the equilibrium measure. In the last paragraph, we explain how to choose the parameters of the Gibbs measure so that the equilibrium measure is a given target distribution π.

Energies and the equilibrium measure

Let d ≥ 1, K : R d × R d → R ∪ {±∞} and V : R d → R ∪ {+∞}.
For reasons that shall become clear shortly, we call K the interaction kernel, and V the external potential. Assumptions on K and V will be given to make the following definition meaningful. Definition 3.1 (Energies). Whenever they are well-defined, we introduce the following quantities, for signed Borel measures µ, ν on R d . The interaction potential, or kernel embedding, of µ, is defined as

U µ K (z) = K(z, y) dµ(y), z ∈ R d .
The interaction energy between µ and ν is defined as

I K (µ, ν) = K(x, y) dµ(x) dν(y).
When µ = ν, we simply write I K (µ) = I K (µ, µ). Finally, we let4

I V K (µ) = 1 2 {K(x, y) + V (x) + V (y)} dµ(x) dµ(y).
The physics-inspired energy vocabulary is useful to the intuition: in a world where the Coulomb interaction is given by K, U µ K (z) would be the electric potential created at point z by charges distributed according to µ. In the same way, I V K (µ) is the energy of points distributed according to µ, repelling each other according to K, and confined by some external potential V. We henceforth denote by E K (respectively E V K ) the set of finite signed Borel measures with finite interaction energy I K (µ) (respectively, with finite energy I V K (µ)). We will work under the following assumptions on K and V . The first one restricts our class of interaction kernels, in particular insisting that points should repel, but that the interaction cannot be singular.

Assumption 1. K is symmetric, non-negative, continuous, and bounded on the diagonal: there exists some constant C ≥ 0 such that K(x, x) ≤ C < ∞ for all x ∈ R d .

Assumption 1 in particular ensures that I K (µ) is well-defined for any probability measure, with possibly infinite value. Our next assumption excludes pathological cases where I K does not induce a distance on probability distributions.

Assumption 2. K is integrally strictly positive definite (ISDP), i.e. I K (µ) > 0 for any non-zero finite signed measure Borel µ.

Assumption 1 and 2 allow most kernels used in machine learning [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], like the Gaussian or isotropic Matern kernels, as well as truncated singular kernels like

K s,ϵ (x, y) = 1 (|x -y| 2 + ϵ 2 ) s/2 ,
for ϵ > 0 and s > 0 [START_REF] Pronzato | Bayesian quadrature, energy minimization, and space-filling design[END_REF].

Under Assumptions 1 and 2, K is finite on the diagonal and ISDP, so that it is in particular positive definite. We can then consider the RKHS H K induced by the kernel K [START_REF] Berlinet | Reproducing kernel Hilbert spaces in probability and statistics[END_REF]). An easy consequence of the Cauchy-Schwarz inequality in H K is that for any x, y ∈ R d , 0 ≤ K(x, y) ≤ C. In particular, I K (µ, ν) and U µ K (z) from Definition 3.1 are well-defined and finite for all finite signed Borel measures; see [START_REF] Pronzato | Bayesian quadrature, energy minimization, and space-filling design[END_REF] for more details. The following known duality formula then links energy minimization and quadrature guarantees for integrands in the unit ball of H K .

Proposition 3.2 [START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF]). Under Assumptions 1 and 2, for probabilities µ, ν in

E K , let γ K (µ, ν) = sup ∥f ∥ H K ≤1 f d(µ -ν) .
(10)

Then γ K (µ, ν) = (I K (µ -ν)) 1/2 .
We now add an assumption to make sure that V is strong enough a confining term. Together with Assumption 1, this ensures that I V K is well-defined for any probability measure5 , again with possibly infinite value.

Assumption 3. V is lower semi-continuous, finite everywhere and V (x) -→ +∞ when |x| -→ +∞. Moreover, there exists a constant c > 0 such that exp (-cV (x)) dx < ∞.

We are now ready to consider the minimizers of I V K .

Proposition 3.3. Let K satisfy Assumptions 1 and 2, and V satisfy Assumption 3. Then 1. I V K is lower semi-continuous, has compact level sets and I V K (µ) > -∞ for any probability distribution µ on R d .

2. If µ ∈ E V K , then I K (µ) and |V | dµ are finite, and

I V K (µ) = 1 2 I K (µ) + V dµ. 3. I V
K is strictly convex on the convex non-empty set E V K .

4. I V K has a unique minimizer µ V over the set of probability measures on R d , called the equilibrium measure, and the support of µ V is compact.

The proof is given in Appendix A. It is a careful assembly of arguments from [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF] and [START_REF] Pronzato | Bayesian quadrature, energy minimization, and space-filling design[END_REF].

Concentration for the Gibbs measure

We saw in Section 1 that herding-like algorithms rely on finding a configuration of points {x 1 , . . . , x n } that minimizes the interaction energy I K ( 1 n n i=1 δ xi -π). In this paper, we rather consider points that are drawn from a distribution that favors small values for an empirical proxy of this interaction energy. To properly define our distribution, consider first, for x 1 , . . . , x n ∈ R d , the discrete energy

H n (x 1 , ..., x n ) = 1 2n 2 i̸ =j K(x i , x j ) + 1 n n i=1 V (x i ), ( 11
)
which we copy here from (7) for ease of reference. Note the similarity, up to diagonal terms, with

I V K ( 1 n δ xi ), where I V K is defined in (3.1).
Definition 3.4. Let K satisfy Assumptions 1 and 2, and V satisfy Assumption 3. Let β n ≥ 2cn, where c is the constant of Assumption 3. The Gibbs measure P V n,βn is the probability measure on R d n defined by

dP V n,βn (X n ) = 1 Z V n,βn exp (-β n H n (X n )) dx 1 ... dx n ,
where β n is called the inverse temperature, X n is short for (x 1 , . . . , x n ), and Z V n,βn is the normalization constant, which, under our assumptions, is finite and positive [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF](Chafaï et al., , 2018)).

When x 1 , . . . , x n ∼ dP V n,βn , we will henceforth denote by µ n = 1 n n i=1 δ xi the associated empirical measure. We saw in Section 2 that µ n converges to the equilibrium measure µ V , with a large deviations principle at speed β n . We are able to give non-asymptotic guarantees on this convergence via a concentration inequality, which is the main result of the paper.

Theorem 3.5 (Concentration inequality). Let K satisfy Assumptions 1 and 2, and V satisfy Assumption 3. Further assume that the associated equilibrium measure µ V has finite entropy. Let β n ≥ 2cn, where c is the constant of Assumption 3. Then there exist constants c 0 , c 1 , c 2 > 0, depending on K and V , such that for any n ≥ 2 and for any r > 0,

P V n,βn I K (µ n -µ V ) > r 2 ≤ exp - c 0 2 β n r 2 + n(c 1 + β n n 2 c 2 ) . ( 12
)
We emphasize again that by Proposition 3.2, ( 12) provides a non-asymptotic confidence interval for the worst-case quadrature error in the unit ball of the RKHS H K . Note that the bound is only interesting in the regime β n /n -→ +∞, where the temperature 1/β n goes down quickly enough. A classical choice of temperature scale is β n = βn 2 where β > 0. We can rephrase a bit to get a more explicit sub-Gaussian decay in the bound.

Corollary 3.6. Under the assumptions of Theorem 3.5, let further β n ≫ n. Then there exist constants u 0 , u 1 > 0 such that for any n ≥ 2 and for any

r ≥ u 0 max n -1/2 , (β n /n) -1/2 , ( 13
)
P V n,βn I K (µ n -µ V ) > r 2 ≤ exp -u 1 β n r 2 . ( 14
)
In particular, when β n ≥ vn 2 for some constant v > 0, Condition (13) simply becomes r > u 0 n -1/2 . The proof of Corollary 3.6 is straightforward from Theorem 3.5, itself proved in details in Appendix A.

We thus recover the known dimension-independent decay in n -1/2 of the worst-case quadrature error as proved for deterministic herding [START_REF] Bach | On the Equivalence between Herding and Conditional Gradient Algorithms[END_REF], though for our probabilistic relaxation only with very large probability, and towards the equilibrium measure. The meaning of very large is that under i.i.d. sampling, analogous results such as Equation (3) feature n instead of β n in the right-hand side of ( 14), and β n /n → +∞. This fast-increasing coverage probability of our confidence interval is a trace of the repulsion introduced in the Gibbs measure.

We now explain in concrete terms how Corollary 3.6 implies that Monte Carlo integration with P V n,βn and with respect to a target distribution π outperforms crude Monte Carlo.

Application to guarantees for probabilistic herding

Let d ≥ 1 and π be a probability measure on R d , which we assume to be our target. To apply Theorem 3.5, we shall work under the following assumption.

Assumption 4. The support S π ⊂ R d of π is compact, and π has finite entropy, in the sense that π has a density π ′ w.r.t. Lebesgue, and that -log π ′ (x)dπ(x) < ∞.

The following proposition shows that, for a given kernel K, we can choose V so that µ V = π, assuming prior computation of the kernel embedding U π K (z) for all z ∈ R d .

Proposition 3.7. Let K satisfying Assumptions 1 and 2 be fixed. Let π satisfy Assumption 4. In particular, there exists R > 0 such that S π ⊂ B(0, R), where B(0, R) is the closed Euclidean ball. Let Φ : R d → R be any continuous, nonnegative function such that Φ = 0 on ∂B(0, R), Φ(z) → +∞ as |z| → +∞ and {|x|>R} e -Φ(x) dx < ∞.

Then, setting V π (z) = -U π K (z) when z ∈ B(0, R) and V π (z) = -U π K (z) + Φ(z) otherwise, V π satisfies Assumption 3 and µ V π = π.
This is a standard result, which relies on the so-called Euler-Lagrange characterization of the equilibrium measure. We give a proof in Appendix A, which is inspired by Corollary 1.4 of [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF], who treat the more difficult case of singular interactions. A classical choice of Φ is Φ(z) = |z| 2 -R 2 . In machine learning terms, Proposition 3.7 says that a suitably penalized kernel embedding is a good choice of confining potential. Of course, this choice of V requires the ability to evaluate the kernel embedding U π K , and we fall back here onto a standard limitation in the herding literature (Chen et al., 2010;[START_REF] Bach | On the Equivalence between Herding and Conditional Gradient Algorithms[END_REF].

With π now our equilibrium measure, Corollary 3.6 implies a uniform quadrature guarantees with respect to π.

Corollary 3.8. Let K satisfying Assumptions 1 and 2 be fixed. Let π satisfy Assumption 4. Set V = V π as in Proposition 3.7, assume that β n ≫ n, and let x 1 , . . . , x n ∼ P V n,βn . Fix further ϵ > 0, and δ ∈ (0, 1). For any n big enough so that n -1/2 ≤ ϵ and exp -cβ n ϵ 2 ≤ δ, with probability larger than 1 -δ,

sup ∥f ∥ H K ≤1 f d(µ n -π) ≤ ϵ. (15)
.

The proof of Corollary 3.8 is a direct application of Corollary 3.6 and Proposition 3.7. Compared to Corollary 2.2, for a fixed uniform worst-case integration error ϵ and confidence level δ, fewer points are required under P V n,βn than under i.i.d samples from π. In particular, β n replaces n in the constraint that links δ and ϵ.

In the (admittedly limited) sense of Corollary 3.8, our probabilistic relaxation of kernel herding provably outperforms i.i.d. sampling, under generic assumptions on the underlying RKHS. In contrast, for deterministic IPM minimization algorithms, we recall that better guarantees than i.i.d sampling were only proved so far in the case where H K is finite-dimensional, which is quite restrictive in that it forbids most classical kernels, such as the Gaussian.

Experiments

We explain how we approximately sample from P V n,βn , and then perform two toy experiments: one to illustrate Theorem 3.5, and one to assess whether our Gibbs measure might come with a better convergence rate for single integrals.

Approximately sampling from P V n,βn

There is no known algorithm to sample from (6) for a generic choice of K and V , so we resort to MCMC, namely the Metropolis-adjusted Langevin algorithm (MALA; [START_REF] Robert | Monte Carlo statistical methods[END_REF]). While [START_REF] Chafaï | Simulating Coulomb gases and log-gases with hybrid Monte Carlo algorithms[END_REF] report having to tame the gradients in their experiments on singular kernels, vanilla MALA has been in our experience enough to get a good approximation to P V n,βn for (smooth) bounded kernels such as the Gaussian and the truncated Riesz kernel. To wit, MALA uses a Metropolis-Hastings Markov kernel with proposal y|y t ∼ N (y t -αβ n ∇H n (y t ), 2αI dn ) ( 16)

where y ∈ (R d ) n , and α is a user-tuned step size parameter. We first show how points sampled from P V n,βn look like in dimension d = 2, for the quadratic potential V : x → |x| 2 /2, and for different scalings of the inverse temperature β n . We consider the truncated logarithmic kernel K(x, y) = -log |x -y| 2 + ϵ 2 , where ϵ is a truncation parameter that we set to ϵ = 10 -2 , and set the number of particles n to 1000. For ϵ = 0, the equilibrium measure is known to be uniform on the unit disk, and we expect it to be close to uniform in the truncated case as well. In -2 each panel of Figure 1, we show the state of the MALA chain after T = 5000 iterations, with step size6 α = α 0 β -1 n , and α 0 is manually tuned at the beginning of each run so that acceptance reaches 50%. We observe that the three empirical measures indeed approximate the uniform distribution on the disk, with more regular spacings as the inverse temperature grows. Figure 1b already shows a more regular arrangement of the points than under i.i.d. draws from the uniform distribution, while the lattice-like structure of Figure 1c is a manifestation of what physicists call crystallization (Serfaty, 2023): the Gibbs measure is concentrated around minimizers of the energy.

0 2 -1 0 1 (a) βn = n 3/2 -2 -1 0 1 2 -1 0 1 (b) βn = n 2 -2 -1 0 1 2 -1 0 1 (c) βn = n 3

Comparing worst-case errors

In this experiment, we take for π the uniform measure on the unit ball of R d , with d = 3. We compare, for various values n of the number of quadrature nodes, the worst-case integration error I K of (i) the empirical measure µ MCMC n of an MH chain of length n targeting π, with an isotropic Gaussian proposal with variance 0.05I d , and of (ii) the empirical measure µ n of an approximate sample of P V n,βn . For the latter, we run MALA for T = 5 000 iterations, with step size tuned as in Section 4.1.

We consider two interaction kernels, the Gaussian kernel

K 1 (x, y) = exp -|x -y| 2 /2
and the truncated Riesz kernel For K ∈ {K 1 , K 2 }, we set V to V π as in Proposition 3.7, with U π K (•) replaced by a Monte Carlo approximation M -1 M i=1 K(•, z i ), where (z 1 , . . . , z M ) are a sample from an MH chain of length M = 1 000 targeting π, independent from any other sample. This approximation is usual in kernel herding experiments [START_REF] Pronzato | Bayesian quadrature, energy minimization, and space-filling design[END_REF].

K 2 (x, y) = |x -y| 2 + 0.1 2 -(d-2)/2 .
Figure 2b and 2a show the results for K 1 and K 2 , respectively. For each value of n, we plot an independent approximation of I K (µ -π) for µ ∈ {µ MCMC n , µ n } the empirical measure of either the baseline MH chain or of the T -th iterate of our MALA chain targeting P V n,βn . The approximation results from writing

I K (µ -π) = I K (µ) -2I K (µ, π) + I K (π), (17)
in which the first term of ( 17) can be computed exactly when µ has finite support, but the other terms require independent MH samples targeting π, here of length 10, 000. We observe that I K (• -π) decays at the same rate under the approximation of P V n,βn as for MCMC samples. This was expected, since our concentration bound (3.6) recovers the n -1 rate for the energy I K (µ n -π) under P V n,βn , the improvement rather being on the sub-Gaussian decay. We see nonetheless that the approximated energy (and hence, the worst-case integration error) is always smaller by about a factor 3 under P V n,βn , which is also expected since P V n,βn favors small values of I K by definition.

Comparing variances for a single integrand

We know from classical CLT arguments that when x 1 , . . . , x n are drawn from an MCMC chain targeting π, the variance of n -1 n i=1 f (x i ) scales like n -1 as n grows, under appropriate assumptions on f . While this is not at all implied by our Corollary 3.6, analogies with the statistical physics literature make us expect a CLT to hold for P V n,βn , at rate β n -1/2 , at least for some temperature schedules (β n ) and smooth enough integrands. Such a result would imply asymptotic confidence intervals for single integrands of width decreasing like β -1/2 n , a faster decay than standard Monte Carlo. To assess whether this expectation is reasonable, we consider the same setting as in Section 4.2: π is uniform on the unit ball in d = 3, the kernel is the truncated Riesz kernel K 2 , and the integrand is f : x -→ K(x, 0), which naturally belongs to H K . For each value of n in Figure 2c, we run 100 independent MH chains of length n targeting π, and plot the empirical variance of the Monte Carlo estimator of f dπ. Similarly, we run 100 independent MALA chains targeting P V n,βn with α tuned as in Section 4.2, for T = 5000 iterations each, and plot the 100 empirical variances obtained from each T -th sample. Again, U π K2 is approximated through long MH chains.

We observe in Figure 2c that the variance is noticeably smaller under the Gibbs measure, for both temperature schedules. Moreover, the rate of decay appears faster, at least in the "usual" temperature schedule β n = n 2 . It is hard to be more quantitative, as the fact we use MALA, and with a fixed number of iterations across all values of n, may impact the convergence we see here. Still, the experiment supports our belief that a fast CLT holds for P V n,βn .

Discussion

Using a Gibbs measure that favors nodes that repel according to a kernel K, we improved on the nonasymptotic worst-case integration guarantees of crude Monte Carlo, through a concentration inequality with a fairly easy proof, at least compared to the classical results in statistical physics that inspired us. When the RKHS of K is infinite-dimensional, such an improvement has yet to be proved for deterministic IPM minimization. A strong argument in favor of a Gibbs measure would be a CLT with a fast rate. We show experimental evidence that supports this expectation. Limitations of our approach that deserve further inquiry are the impact of using an approximation of the kernel embedding U π K and an MCMC sampler, here MALA. Integrating a tractable approximation of U π K without loss on the convergence speed would be an important improvement. Simultaneously, understanding how the accuracy of the MALA approximation relates to n and β n would help find the right trade-off between statistical accuracy and computational cost.

Finally, compared to the bound (5) for volume sampling, our concentration bound features β n , but not the eigenvalues of the kernel operator. In that sense, our confidence intervals are likely to be looser in an RKHS with fast-decaying kernel eigenvalues. Yet, establishing a CLT for the estimator in (5), where the weights in the estimator depend on all quadrature nodes, promises to be particularly hard. Moreover, if both distributions are approximately sampled through MCMC, one evaluation of our Gibbs density is quadratic in n, while it is cubic for volume sampling. A careful experimental comparison at fixed budget would thus be interesting.

Figure 1 :

 1 Figure 1: Three independent approximate samples of P V n,βn corresponding to different temperature schedules.

Figure 2 :

 2 Figure 2: Energy and variance comparisons between approximate samples from P V n,βn and MCMC samples.

Actually, the original herding algorithm was inspired by the zero-temperature limit of a physical particle system(Welling,

2009), so our relaxation is a return to the roots of sorts.

By singular, we mean that K(x, x) = +∞ for all x ∈ R d .

Throughout the section, we neglect sampling costs. While a fast-growing βn implies better theoretical guarantees, the price of (approximately) sampling from P V n,βn intuitively increases with βn, introducing a trade-off in practice; see Section

Note that with our convention I 0 K = 1 2 I K .

Unlike I K , we shall only evaluate I V K on probability measures.

Having α decrease at least as β -1 n intuitively avoids the distance between two consecutive MALA states to grow with n; see the MALA proposal (16).

Acknowledgments

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).

A Proofs

A.1 Proof of Proposition 3.3

The proof is a simple application of known results.

1. This is a consequence of the first point of Theorem 1.1 in [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF]; 2. This is given by the second point of Lemma 2.2 in [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF]; 3. This is a consequence of the ISDP assumption on the kernel, using Lemma 3.1 in [START_REF] Pronzato | Bayesian quadrature, energy minimization, and space-filling design[END_REF]. To see that E K is non-empty, we can simply consider Dirac measures δ x , since K is bounded and V is finite everywhere.

4. This is a consequence of points 1 and 3 using the arguments of Section 4.1 of [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF].

A.2 Proof of Proposition 3.7

The proof of Proposition 3.7 relies on the so-called Euler-Lagrange equations, which we recall here.

Lemma A.1. Let K be a kernel satisfying Assumptions 1 and 2, and V be an external potential satisfying Assumption 3. Set C V = I K (µ V ) + V dµ V . Then µ V has compact support, and a probability measure ν satisfies ν = µ V if and only if ν has compact support and there exists a constant C such that

In that case, C = C V .

Proof. Let us first check that µ V satisfies this characterization. We already know from Proposition 3.3 that µ V exists, is unique, and has compact support. We can use the same procedure as [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF], Theorem 1.2, proof of item 5) : considering the directional derivative of I V K and using the fact that µ V is the minimizer, their equation (4.5) yields that, for any probability distribution ν ∈ E V K ,

Since Dirac measures δ z have finite interaction energy I V K (δ z ), we get Point (i) by taking ν = δ z . The second point is obtained exactly as in the second part of the proof of item 5 of Theorem 1.2 of [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF]. Finally, the converse implication can be similarly obtained along the lines of the proof of item 6 of Theorem 1.2 in [START_REF] Chafaï | First-order global asymptotics for confined particles with singular pair repulsion[END_REF], using the strict convexity of the energy functional in Proposition 3.3.

We are now ready to prove Proposition 3.7, by checking that the Euler-Lagrange equations are satisfied for π and V π , and that V π satisfies the assumptions of Proposition 3.3.

First note that since the kernel K is nonnegative and bounded on the diagonal by assumption, Cauchy-Schwarz in the RKHS H K implies that K(x, y) ≤ C for all x, y ∈ R d . As a consequence, since K is further assumed to be continuous, Lebesgue's dominated convergence theorem yields that z -→ U π K (z) is finite everywhere and continuous. In particular, V π is continuous.

Moreover, the bound on

Finally, the integrability assumption of Assumption 3 is satisfied by our assumption on Φ and because U π K is bounded. Hence, V π satisfies Assumption 3 and the equilibrium measure µ V π is well-defined. We conclude upon noting that, since Φ ≥ 0, Lemma A.1 yields that µ V π = π.

A.3 Proof of Theorem 3.5

The proof will follow the one of (Chafaï et al., 2018), with some notable simplifications. We first compute a lower bound on the partition function Z V n,βn , which generalizes the one of (Chafaï et al., 2018) to bounded kernels.

Proposition A.2. Let K be a kernel satisfying Assumptions Assumption 1 and Assumption 2, and V be an external potential satisfying Assumption 3. Assume that the associated equilibrium measure µ V has finite entropy, i.e. S(µ V ) = -log µ ′ V dµ V < ∞, where µ ′ V is the density of µ V w.r.t. the Lebesgue measure. Then for n ≥ 2, we have

Proof. The idea is to rephrase a bit the discrete energy H n and use Jensen's inequality. We let X n = (x 1 , . . . , x n ) for brevity, and we start by writing

where

Using the definition I V K (µ V ) = 1 2 I K (µ V ) + V dµ V we get the result. To bound I K (µ n -µ V ), we shall further use the following lemma, which is again inspired by (Chafaï et al., 2018).

Lemma A.3. Let K and V be a kernel and an external potential satisfying Assumptions 1 to 3. Let µ be any probability measure of finite energy, I V K (µ) < +∞. Then

Proof. We write

With the notation of Lemma A.1, we know that

Using Fubini, we thus see that

Plugging this into (18) yields the result.

We are now ready to prove Theorem 3.5. Recall that we work under the assumption β n ≥ 2cn where c is the constant of Assumption 3. Consider a Borel set A ⊂ R d n . For brevity, recall that for (x k ) k≤n ∈ A, we write µ n = 1 n n i=1 δ xi and X n = (x 1 , . . . , x n ).

Let η > 0 and X n ∈ A. The key idea is to split

The first part will be compared with the energy I V K (µ n ), while the second part will be kept to ensure integrability. Remembering that

we write, by definition,

Using Proposition A.2, we continue

As noted in the proof of Lemma A.1, 0 ≤ K ≤ C, so that the last integral in ( 23) is easily bounded,

Now we choose a particular value for η, namely η = cn/β n where c is the constant of Assumption 3. Further let C 2 = log R d exp (-cV (x)) dx, which is finite by assumption. Setting c 1 = cI V K (µ V ) + C 2 -S(µ V ) and c 2 = 1 2 C -1 2 I K (µ V ), (23) yields

Note that c 1 and c 2 are indeed finite by definition of µ V , since E V K is non-empty. The comparison inequality of Lemma A.3 applied to (24), with

yields Theorem 3.5, upon noting that β n -cn ≥ β n /2 by assumption.