Quenched large deviations for Monte Carlo integration with Coulomb gases

Kernel herding belongs to a family of deterministic quadratures that seek to minimize the maximum mean discrepancy (MMD), that is, the worst-case integration error over a reproducing kernel Hilbert space (RKHS). These MMD minimization procedures come with strong experimental support, but comparatively less theoretical footing. In particular, apart from recent progress in distribution compression, little has been proved in favor of an improvement of MMD minimization over classical Monte Carlo quadrature when the RKHS is infinite-dimensional. In this paper, we study a joint probability distribution over quadrature nodes, a tailored Gibbs distribution, whose support intuitively tends to concentrate around MMD minimizers as a temperature parameter is decreased. Our main contribution is to prove that drawing integration nodes from our distribution does outperform i.i.d Monte Carlo. While our bounds on the worst-case integration error feature the same rate as i.i.d. Monte Carlo, we do obtain a tighter concentration inequality as the temperature parameter decreases. This means smaller confidence intervals as the number of quadrature nodes increases. While arguably a first step, our results demonstrate that the mathematical toolbox developed around Gibbs measures can help understand to what extent kernel herding and its variants improve on computationally cheaper methods. There remains the issue of sampling from our Gibbs distribution. In our numerical experiments, we demonstrate that a simple MCMC chain already yields approximate samples that lead to improved confidence intervals around the target integrals, as supported by our theoretical results.

The standard Monte Carlo estimator $\widehat{I}_N^{\mathrm{MC}}$ of $\int fdω$ relies on independent samples from $ω$ and has variance of order $1/N$. Replacing the samples with a determinantal point process (DPP), a repulsive distribution, makes the estimator consistent, with variance rates that depend on how the DPP is adapted to $f$ and $ω$. We examine two existing DPP-based estimators: one by Bardenet & Hardy (2020) with a rate of $\mathcal{O}(N^{-(1+1/d)})$ for smooth $f$, but relying on a fixed DPP. The other, by Ermakov & Zolotukhin (1960), is unbiased with rate of order $1/N$, like Monte Carlo, but its DPP is tailored to $f$. We revisit these estimators, generalize them to continuous settings, and provide sampling algorithms.

Unadjusted samplers such as unadjusted Hamiltonian Monte Carlo and underdamped Langevin are well-known to be biased. Metropolis--Hastings adjustment has been conventionally incorporated into Hamiltonian Monte Carlo to eliminate the bias. However, this adjustment can significantly increase the iteration complexity due to the small step size required for reasonable Metropolis acceptance rates. In this work, we extend the \emph{delocalization of bias} phenomenon, previously established for the overdamped Langevin algorithm, to these two unadjusted algorithms. We show that to control the $W_2$ bias of any $K$-dimensional marginal of a high-dimensional distribution, $O(\sqrt{K})$ integration steps suffice up to $\log d$ terms, assuming either weak or sparse interactions among variables. The discrete-time integrators here introduce technical difficulties beyond those of the overdamped setting, which we address through a broadly applicable matrix-polynomial framework that characterizes their propagators. Our result for the underdamped Langevin algorithm is valid for all large friction parameters, implying that the Leimkuhler-Matthews integrator for the overdamped Langevin dynamics also exhibits delocalization of bias.