Source-linked AI summary
Stationary MMD Points
Zonghao Chen, Toni Karvonen, Heishiro Kanagawa, François-Xavier Briol, Chris. J. Oates
TL;DR
Selecting finite points for numerical integration by minimizing MMD is difficult because the objective is non-convex and existing algorithms cannot provably recover its true minimum. This paper analyzes stationary MMD points and noisy MMD gradient flow, showing faster-than-MMD integration error for RKHS integrands and proving that the flow can compute stationary points.
Problem
Finite-point numerical integration requires careful node selection, but MMD non-convexity leaves existing algorithms unable to provably recover true minimum-MMD points.
Source: §1 · p. 1Method
Results
For a wide class of RKHS integrands, stationary MMD points achieve numerical integration error that vanishes faster than the MMD worst-case bound.
Source: §1 · p. 2Takeaways & Limitations
Stationary points integrate exactly over a large linear subspace of the RKHS, and that subspace expands as the particle count increases.
Source: §1 · p. 2The convergence analysis assumes controlled noise injection satisfying gradient dominance in expectation, with smooth-kernel assumptions including Gaussian, qualifying Matérn, and inverse multiquadratic kernels.
Source: §4 · p. 6
Abstract
from arXiv · showhide
Approximation of a target probability distribution using a finite set of points is a problem of fundamental importance in numerical integration. Several authors have proposed to select points by minimising a maximum mean discrepancy (MMD), but the non-convexity of this objective typically precludes global minimisation. Instead, we consider the concept of \emph{stationary points of the MMD} which, in contrast to points globally minimising the MMD, can be accurately computed. Our main contributions are two-fold and theoretical in nature. We first prove the (perhaps surprising) result that, for integrands in the associated reproducing kernel Hilbert space, the numerical integration error of stationary MMD points vanishes \emph{faster} than the MMD. Motivated by this \emph{super-convergence} property, we consider MMD gradient flows as a practical strategy for computing stationary points of the MMD. We then prove that MMD gradient flow can indeed compute stationary MMD points, based on a refined convergence analysis that establishes a novel non-asymptotic finite-particle error bound.
1. Introduction
The paper studies selecting finite point sets to approximate a target distribution for numerical integration, using MMD as a measure of discrepancy. Because non-convexity makes global MMD minimizers difficult to recover, it analyzes stationary MMD points and how they can be computed.
1. Introduction
- Finite point sets approximate a target distribution for numerical integration, including posterior expectations in Bayesian inference and representative-point selection for coresets.The paper considers uniformly weighted empirical measures and MMD-based approximation.Source: §1 · p. 1
- MMD connects point-set discrepancy to integration error for functions in the kernel's RKHS, making it a useful criterion for selecting numerical integration nodes.The paper also notes MMD's practicality and that kernel choice can induce a wide range of topologies.Source: §1 · p. 1
- The non-convex MMD objective leaves a gap between theory for minimum-MMD points and algorithms that often find only local minima.The paper asks why these locally obtained point sets can nevertheless perform well in numerical integration.Source: §1 · pp. 1–2
- Stationary MMD points satisfy convergence of MMD to zero as n grows and zero MMD gradient with respect to the points at fixed n.Minimum-MMD points also satisfy these properties, but stationary points form a strict relaxation that can be obtained by a specific noisy gradient descent scheme.Source: §1 · p. 2
- Stationary MMD points achieve faster-than-MMD numerical integration error for a wide class of RKHS functions.The paper interprets this as super-convergence: stationary points integrate exactly over a large linear subspace of the RKHS, and this subspace expands with n.Source: §1 · p. 2
- The second contribution proves that an interacting particle system based on MMD gradient flow can compute stationary MMD points when simulated until stationarity.The analysis establishes a non-asymptotic convergence rate for time-discretized, finite-particle flow.Source: §1 · p. 2
2. Related Work
Related work studies point-set approximation through MMD, energy objectives, quasi-Monte Carlo, compression methods, and gradient flows. The paper positions stationary MMD points as a computable alternative with theoretical guarantees and applicability when no candidate pool is available.
Related work
- Uniformly weighted empirical measures are the review's focus; weighted alternatives are discussed separately, and uniform weights are motivated by simplicity, numerical stability, and robustness to assumption violations.Source: §2 · p. 2
- Minimum energy points minimize an energy functional, and that functional coincides up to an additive constant with squared MMD when the interaction function is a kernel.The cited work on support points gives an integration-error rate, while general algorithms are not guaranteed to find minimum-energy configurations.Source: §2 · pp. 2–3
- QMC can converge quickly and be explicitly computed for common distributions, but does not currently address general target distributions on Rd.The discussion distinguishes traditional star discrepancy from the more general MMD perspective.Source: §2 · p. 3
- Prior MMD-minimization methods include kernel herding, greedy candidate selection, and kernel thinning, but their properties and candidate-pool constraints differ.Kernel herding and greedy MMD minimization select from large candidate sets; kernel thinning selects representative points from a dataset.Source: §2 · p. 3
- MMD gradient flow has been used in generative modelling and numerical integration, but formal convergence analysis for time-discretized finite-particle flow was missing.The paper identifies non-convexity in the Wasserstein metric as the primary challenge and addresses the finite-particle setting.Source: §2 · p. 3
- Unlike kernel thinning, stationary MMD points are not restricted to a fixed candidate pool, while the paper claims strong theoretical guarantees relative to greedy optimization methods.The paper also contrasts its method with QMC, which does not currently solve the general-target setting described there.Source: §2 · p. 3
3. Numerical Integration with Stationary MMD Points
Stationary MMD points exactly integrate a kernel-gradient subspace, and its increasing approximation capacity yields integration error that vanishes faster than MMD.
3. Numerical Integration with Stationary MMD Points
- Stationary MMD points integrate every function in Fn := span({1} ∪ GXn) exactly.Here GXn is the span of kernel partial derivatives evaluated at the selected points.Source: §3 · p. 4
- Under Assumption 1, connected support, and a C0-universal kernel, stationary MMD points achieve o(MMD(µ, µn)) error for f=c+h with h∈HX non-constant on X.The result uses the seminorm |f|n tending to zero as n increases.Source: §3 · p. 5
- The proof decomposes f into an exactly integrated component in Fn and a residual whose seminorm vanishes as Fn increasingly approximates relevant functions in H.This approximation-capacity result assumes connected support and a C0-universal kernel.Source: §3 · pp. 4–5
- The super-convergence result shows that non-uniform weights are not required to attain super-convergence for numerical integration.Source: §3 · p. 5
- Stationary MMD points need not be unique, even up to permutation, but the integration guarantees apply to any configuration reached by the particle system.Source: §3 · p. 5
4. Computing Stationary MMD Points
Noisy MMD particle descent is presented as a practical way to compute stationary MMD points, unlike global MMD minimisation, which can be obstructed by non-convexity. Under kernel and noise conditions, the analysis bounds finite-particle convergence and shows that the method can reach stationarity.
4. Computing Stationary MMD Points
- Noisy MMD particle descent perturbs MMD gradient-flow updates with injected noise and is shown to compute stationary MMD points.The update is a discretised Wasserstein gradient flow on MMD; the noise scale is β_t.Source: §4 · p. 5
- Assumption 5 controls noise so noisy descent satisfies gradient dominance in expectation, and it can be checked on the current realised trajectory.Selecting β_t from candidates using an empirical check adds per-iteration computational cost.Source: §4 · p. 6
- Theorem 4.1 bounds MMD error by optimisation and finite-sample estimation terms, with larger noise accelerating optimisation while increasing update variance.The optimisation term decays exponentially under the gradient-dominance condition; the estimation term scales as 1/√n.Source: §4 · p. 6
- With β_t tending to zero, Theorem 4.2 establishes that noisy descent can compute stationary MMD points in the joint n, T → ∞ limit.Theorem 4.1 assumes kernel regularity, Assumption 5, and a sufficiently small step size.Source: §4 · pp. 6–7
- The method has O(n^2T) computational cost, and the stated iteration choice gives O(n^2) cost up to logarithmic factors and oP(n^-1/2) integration error.The comparison describes stationary MMD points as more general than QMC, while sharing computational complexity with approaches that minimise MMD or energy functionals.Source: §4 · p. 7
- With a good initial set satisfying S(1) = O(n^-1(log n)^b), finite-particle descent yields MMD convergence O(n^-1(log n)^b) and single-function integration error o(n^-1(log n)^b).The passage attributes the integration improvement to stationarity enforced by perturbing the initial point set.Source: §4 · p. 7
5. Experiments
The experiments test noisy MMD particle descent on a Gaussian mixture and two OpenML datasets. Stationary MMD points show strong integration performance, with results varying by integrand and setting.
5.1. Mixture of Gaussians
- The noisy particle descent experiment verifies exact integration for the tested functions in Fn at stationarity, although good integration performance occurs in fewer iterations than full convergence requires.The reported setup uses a 10-component, two-dimensional Gaussian mixture, a Gaussian kernel, and noise level βt = t^−1/2; Figure 4 indicates Assumption 5 holds for most iterations except the first few.Source: §5.1 · pp. 7–8, Figure 4 · p. 9
- On the Gaussian mixture, stationary MMD points outperform every baseline for f1 and match support points for f2 while outperforming the other baselines.For n = 100, their computation takes about 20 seconds, versus more than 3000 seconds for kernel thinning; estimated rates are −1.39 for f1 integration error and −1.35 for MMD error.Source: §5.1 · p. 8
- For Gaussian mixtures with d = 10 and d = 50, stationary MMD points consistently outperform i.i.d. sampling across all tested sample sizes.The results are summarized in Table 1, which is not among the supplied passages.Source: §5.1 · p. 8
5.2. OpenML Datasets
- Stationary MMD points outperform all baselines for f1 and MMD on House8L and Elevators, with kernel thinning the closest but computationally expensive baseline.For f2, the authors do not report results because f2 is not in the RKHS.Source: §5.2 · p. 9
6. Discussion
The discussion links stationary MMD points to the practical gap between theoretical minimum-MMD performance and what can be computed. It presents super-convergence and a convergence analysis for noisy particle descent as key advances.
6. Discussion
- For integrands in the associated RKHS, stationary MMD points have numerical integration error that vanishes faster than the MMD, a property termed super-convergence.The paper offers this result as an explanation for previously observed strong empirical performance.Source: §6 · p. 9
- The analysis proves that noisy MMD particle descent can compute stationary MMD points, strengthening the convergence analysis for the method.Source: §6 · p. 9
- The authors identify MMD-based coreset construction in large-scale deep learning as a potential application for stationary MMD point sets.They describe this as a promising avenue for future work.Source: §6 · p. 9
Appendices
The appendices review related approximation methods and provide technical results supporting the paper’s theory. These include RKHS restriction properties and density results used in the analysis of stationary MMD points.
Appendices
- The appendices outline the organization of supplementary material: proofs are collected in Section B and reproducibility details in Section C.Source: Appendices · p. 14
A. Further Related Work
- The related-work appendix surveys kernel cubature, DPPs, information-based complexity, and other weighted integration approaches.It notes that optimal kernel-cubature weights can be costly and unstable, while DPPs have limitations in producing uniformly weighted approximations for general targets.Source: A. Further Related Work · pp. 14–15
B.1. Preliminary Results
- The appendix supplies preliminary arguments used in the proofs, including a relation between kernel cubature error and MMD for weighted measures.It states that under weak conditions the cubature error is o(MMD).Source: A. Further Related Work · p. 14, §B.1 · p. 15
B.1.1. INTERCHANGE INTEGRALS AND DERIVATIVES
- The appendix states that differentiation and integration can be interchanged under Assumption 1, yielding the gradient of the kernel integral.Source: §B.1.1 · p. 16
B.1.2. ON THE SUPPORT X AND THE RKHS SUBSET HX
- The subspace HX is the orthogonal complement of the null space of restriction to X, and restriction is an isometry from HX to the RKHS on X.The appendix also explains that the main super-convergence results exclude functions constant on X because their integrals are trivial; for some RKHSs, such as the Gaussian RKHS on sets with non-empty interior, no such constant function exists.Source: §B.1.2 · p. 16
B.1.3. DENSENESS OF STATIONARY MMD POINTS
- When empirical measures converge to μ in MMD under Assumption 3, their points eventually approach every point in the support X.More precisely, for each x in X and ε > 0, sufficiently large point sets contain a point within 2ε of x.Source: §B.1.3 · pp. 16–17
B.2. Statement and Proof of Theorem B.5
The appendix establishes density properties for derivative-kernel spans that support the main theoretical argument. It first shows density in HX, then approximates that span using derivatives at stationary points.
B.2. Statement and Proof of Theorem B.5
- Under Assumptions 1 and 2, the span GX of kernel derivatives at points in X is dense in HX for functions nonconstant on X.Source: §B.2 · p. 17
- The density proof uses connected paths in X: a function with zero gradient throughout X must be constant there, and the stated function class then forces it to be zero.Source: §B.2 · p. 17
B.3. Statement and Proof of Theorem B.6
- Under Assumptions 1–3, spans of kernel derivatives at stationary points approximate any element of GX arbitrarily well for sufficiently large n.Source: §B.3 · p. 17
B.4. Proof of Theorem 4.1
The proof of Theorem 4.1 combines expected descent of MMD with concentration of the noisy particle updates. A union bound and discrete Grönwall argument yield the stated finite-particle bound, which also supports a logarithmic-iteration corollary under specified initialization and step-size conditions.
B.4. Proof of Theorem 4.1
- Theorem 4.1 is proved from an expected descent lemma and a high-probability bound keeping noisy particle updates close to their expected values.The descent lemma gives monotone expected MMD decrease; the concentration result controls deviations after Gaussian noise injection.Source: §B.4 · p. 18
- A union bound over iterations and the discrete Grönwall lemma combine the per-step estimates into a finite-time bound.The proof then chooses τ = δ + log T to obtain the stated probability level.Source: §B.4 · p. 19
- With initial MMD O(n^-1(log n)^b), γβ_t = n^-1, and Assumption 5, the corollary gives a bound after T = O((log n)^b) iterations.The supplied corollary states these conditions and iteration count for the final MMD bound.Source: §B.4 · p. 20
B.4.1. PROOF OF THEOREM B.7
- Under the stated step-size condition, the descent argument establishes E(1) − E(0) ≤ −γβ²κd²E(0).The proof identifies 256γ²d²κ ≤ 1 as the condition used in this step.Source: §B.4.1 · p. 23
B.4.2. PROOF OF THEOREM B.8
The proof of Theorem B.8 bounds deviations caused by the injected Gaussian noise. It uses concentration for independent, bounded, zero-mean Hilbert-space-valued random variables, then combines the estimate with triangle and Jensen inequalities.
B.4.2. PROOF OF THEOREM B.8
- Conditioned on the current particles, the update randomness comes from injected Gaussian noise, so the probability bound is over that noise.The proof treats the particles as fixed when analyzing the update.Source: §B.4.2 · p. 23
- The noise-induced feature-map deviations are independent, zero-mean Hilbert-space-valued variables with norm bounded by 2√κ.The bound uses the kernel condition sup_x k(x, x) ≤ κ.Source: §B.4.2 · p. 23
- Hilbert-space Bernstein concentration gives a high-probability bound, which the proof completes using triangle and Jensen inequalities.For τ > 1, the cited probability level is at least 1 − exp(−τ).Source: §B.4.2 · p. 24
B.4.3. AUXILIARY RESULTS
The auxiliary results provide regularity bounds used in the convergence proof, including a Lipschitz bound for Φ. They also state a discrete Grönwall lemma for controlling recursive inequalities.
B.4.3. AUXILIARY RESULTS
- The auxiliary inequalities include a Lipschitz bound for Φ: changing its two arguments changes it by at most 2dκ times the sum of their distances.The stated bound is ∥Φ(z, w) − Φ(z′, w′)∥ ≤ 2dκ(∥z − z′∥ + ∥w − w′∥).Source: §B.4.3 · p. 24
- The discrete Grönwall lemma bounds a recursively defined a_T from the per-step factors λ_t and additions b_t.The lemma states a_T ≤ ∏_{t=1}^{T−1}(1 + λ_t) times the corresponding accumulated contributions.Source: §B.4.3 · p. 25
C. Additional Experimental Details
The additional experiments specify implementation choices for mixture and OpenML settings, including iteration counts, step sizes, kernels, and baseline configurations. In the Matérn-kernel ablation, results are reported as consistent with the Gaussian-kernel experiments, with a House8L advantage for stationary MMD points.
C. Additional Experimental Details
- For the mixture experiment, stationary MMD points use noisy particle descent with 10,000–300,000 iterations depending on n, fixed γ = 1.0, and kernel lengthscale ℓ = 1.The same lengthscale is used for all methods to avoid tuning bias; the authors note that problem-specific tuning may further improve performance.Source: C. Additional Experimental Details · p. 27
- For House8L and Elevators, the target is the normalized dataset’s empirical measure, and stationary points use 100,000 iterations with dataset-specific step sizes.The step sizes are γ = 1.0 for House8L and γ = 0.3 for Elevators; QMC is not applicable because the dataset domains are unknown.Source: C. Additional Experimental Details · p. 27
C.1. Ablation study with Mat´ern- 3 2 kernel
- With the Matérn-3/2 kernel at ℓ = 1, the ablation results are consistent with the Gaussian-kernel results, and stationary MMD points outperform the Gaussian-kernel version on House8L.The authors leave a more comprehensive investigation of kernel choices for future work.Source: §C.1 · p. 27
- Support points are unchanged in the kernel ablation because that baseline does not involve a choice of kernel.The Matérn-3/2 kernel satisfies the assumptions required for stationary MMD points.Source: §C.1 · p. 27
C.2. More Experiments on Kernel Thinning
The experiments compare centered and uncentered Gaussian kernels and vary the oversampling parameter for KT. Centering improves KT, while increasing oversampling improves performance at added sample-access and computational cost.
C.2. More Experiments on Kernel Thinning
- The earlier KT implementation estimated MMD from samples, unlike kernel herding and stationary MMD points, which use closed-form kernel mean embeddings.Source: §C.2 · pp. 27–28
- In the mixture-of-Gaussians setting, both methods improve as g increases across 0, 2, 4, and 6, but larger g requires 2gN samples and raises computational cost.Source: Figure 6 · p. 28, §C.2 · p. 28
- The centered kernel subtracts the two µ-averaged kernel terms and adds their joint expectation, giving it zero integral under µ.Source: §C.2 · p. 28
- Centered Gaussian kernels greatly improve KT relative to uncentered kernels across the datasets compared in Figure 7.Source: §C.2 · p. 28, Figure 7 · p. 29