stat.MLFeb 25, 2026
SaveScalable Kernel-Based Distances for Statistical Inference and Integration
Organizations: Department of Computer Science University College London
Abstract
Representing, comparing, and measuring the distance between probability distributions is a key task in computational statistics and machine learning. The choice of representation and the associated distance determine properties of the methods in which they are used: for example, certain distances can allow one to encode robustness or smoothness of the problem. Kernel methods offer flexible and rich Hilbert space representations of distributions that allow the modeller to enforce properties through the choice of kernel, and estimate associated distances at efficient nonparametric rates. In particular, the maximum mean discrepancy (MMD), a kernel-based distance constructed by comparing Hilbert space mean functions, has received significant attention due to its computational tractability and is favoured by practitioners. In this thesis, we conduct a thorough study of kernel-based distances with a focus on efficient computation, with core contributions in Chapters 3 to 6. Part I of the thesis is focused on the MMD, specifically on improved MMD estimation. In Chapter 3 we propose a theoretically sound, improved estimator for MMD in simulation-based inference. Then, in Chapter 4, we propose an MMD-based estimator for conditional expectations, a ubiquitous task in statistical computation. Closing Part I, in Chapter 5 we study the problem of calibration when MMD is applied to the task of integration. In Part II, motivated by the recent developments in kernel embeddings beyond the mean, we introduce a family of novel kernel-based discrepancies: kernel quantile discrepancies. These address some of the pitfalls of MMD, and are shown through both theoretical results and an empirical study to offer a competitive alternative to MMD and its fast approximations. We conclude with a discussion on broader lessons and future work emerging from the thesis.
Figures & tables
| The set of natural numbers, . | |
| The subset , for . | |
| The identity matrix. | |
| Tuple , for . | |
| Vector (for ) or matrix , for . | |
| , | The norm and inner product of a space . |
| Column vector , for and . |
Table 2.1 : Summary of notation used in the thesis.
| Polynomial of degree | , for and . |
| Matérn of order , and smoothness | , where for is the modified Bessel function of the second kind. |
| Gaussian | |
| Brownian motion | |
| Fractional Brownian motion |
Table 2.2 : Example kernels on . For Brownian motion and fractional Brownian motion, .
Figure 3.1 : Estimating the MMD requires approximating the embedding of the model in the RKHS . The classical approach approximates it using equally-weighted i.i.d. samples from , denoted . We show that this estimator can be improved by using optimally-weighted samples, denoted .
| Model | References | i.i.d. V-stat | i.i.d. OW (ours) | RQMC V-stat | RQMC OW (ours) | ||
|---|---|---|---|---|---|---|---|
| g-and-k | 1 | 1 | [ Bharti et al., 2022b , Niu et al., 2023 ] | 2.25 (1.52) | 0.086 (0.049) | 0.060 (0.037) | 0.059 (0.037) |
| Two moons | 2 | 2 | [ Lueckmann et al., 2021 , Wiqvist et al., 2021 ] | 2.36 (1.94) | 0.057 (0.054) | 0.056 (0.044) | 0.055 (0.044) |
| Bivariate Beta | 5 | 2 | [ Nguyen et al., 2020 , Niu et al., 2023 ] | 2.13 (1.17) | 0.555 (0.227) | 0.222 (0.111) | 0.193 (0.088) |
| MA(2) | 12 | 10 | [ Marin et al., 2011 , Nguyen et al., 2020 ] | 2.42 (0.796) | 0.705 (0.107) | 0.381 (0.054) | 0.322 (0.052) |
| M/G/1 queue | 10 | 5 | [ Pacchiardi et al., 2024 , Jiang, 2018 ] | 2.52 (1.19) | 1.71 (0.568) | 0.595 (0.134) | 0.646 (0.202) |
| Lotka-Volterra | 600 | 2 | [ Briol et al., 2019b , Wiqvist et al., 2021 ] | 2.13 (1.10) | 2.04 (0.956) | 1.44 (0.955) | 1.42 (0.942) |
Table 3.1 : Average and standard deviation (in parentheses) of estimated MMD between and over 100 runs for the V-statistic and our optimally-weighted (OW) estimator. Settings: , .
Figure 3.2 : Error in estimating MMD 2 for the multivariate g-and-k distribution. (a) Error of our OW estimator for different choices of and . Increasing the smoothness of improves the performance. (b) Comparison of V-statistic and OW estimator as a function of dimension. OW performs better for both parameterisations of , with the Gaussian giving lowest error. (c) Value of also impacts the performance of the OW estimator. (d) Error vs. total computation cost for different . OW performs better than V-statistic for similar cost: for V-statistic, whereas for OW.
| Cases | i.i.d. V-stat | i.i.d. OW (ours) |
|---|---|---|
| (null holds) | 0.040 | 0.047 |
| (alternative holds) | 0.040 | 0.413 |
Table 3.2 : Fraction of repeats for which the null was rejected. An ideal test would have when the null holds, and otherwise.
Figure 3.3 : ABC posteriors for the wind farm model. Our OW estimator yields posterior samples that are more concentrated around the true than the V-statistic. Settings: , .
Figure 4.1 : Illustration of conditional Bayesian quadrature (CBQ). Stage one fits a GP to for each and integrates to obtain BQ estimates . Stage two places a GP over and fuses these to yield with posterior uncertainty shown by shaded regions.
Figure 4.2 : Bayesian sensitivity analysis for linear models. Left: RMSE of all methods when and . Middle: RMSE of all methods when and . Right: RMSE of all methods when .
Figure 4.3 : Bayesian linear model sensitivity analysis in .
Figure 4.4 : Bayesian sensitivity analysis for SIR Model Option pricing in mathematical finance. Left: RMSE of all methods for the SIR example with . Middle: The computational cost (in wall clock time) for CBQ ( ) and for obtaining one single numerical solution from SIR under different discretisation step sizes. In practice, the process of obtaining samples from SIR equations is repeated times. Right: RMSE of all methods for the finance example with .
Figure 4.5 : Uncertainty decision making in health economics. We study RMSE for different estimators of EVPPI.
Figure 5.1 : Top row: BQ of a fractional Brownian motion (fBm) integrand with Hurst parameter (smoothness ), using the Brownian motion kernel ( 5.3 ) with amplitudes (left), given by the LOO-CV estimator (middle), and given by the ML estimator (right). Vertical lines mark the true value of the integral (solid line) and the BQ posterior mean (dashed line); the shaded bands show the 95% BQ credible interval . Bottom row: The GP interpolation underlying the top row: an fBm path that is the true integrand (solid line), the training data , and the GP posterior mean (dashed line); the shaded bands show the 95% GP credible interval .
Figure 5.2 : BQ of an integrated fractional Brownian motion (fBm) integrand with Hurst parameter (smoothness ), using the Brownian motion kernel ( 5.3 ) with amplitudes (left), given by the LOO-CV estimator (middle), and given by the ML estimator (right). For the explanation of the figures, see the caption of Figure 5.1 .
Figure 5.3 : Rates of decay for the ML, CV and ICV estimators from Theorems 14 , 15 and 20 . Observe that the CV estimator’s range of adaptation to the smoothness is wider than the ML estimator’s, and the ICV estimator’s range of adaptation is wider than that for both the CV and ML estimators.
Figure 5.4 : Expected decay rates for the ML, CV and ICV estimators from Theorems 18 , 19 and 21 . Observe that the CV estimator’s range of adaptation to the smoothness is wider than the ML estimator’s, and the ICV estimator’s range of adaptation is wider than that for both the CV and ML estimators.
Figure 5.5 : Asymptotics of CV estimators for functions of varying smoothness as quantified by and in ( 5.22 ). Runs on individual 100 samples from are in grey, means and confidence intervals (of two standard deviations) are in black.
Figure 5.6 : Asymptotics of CV estimator compared to asymptotics of ML estimators, for once differentiable functions.
Figure 6.1 : Illustration of bivariate quantiles. Left: Bivariate distribution . Center: Density of the projection of onto direction on the unit circle, with . Right: Different quantiles for all possible directions .
Figure 6.2 : Illustration of the impact of the slicing direction on KQEs. Suppose , the KQEs are obtained by considering the -th quantile of . Clearly, these quantiles might vary significantly depending on the slicing direction used.
Figure 6.3 : Experimental results comparing our proposed methods with baseline approaches. Methods represented by dotted lines exhibit quadratic complexity for a single computation of the test statistic, while the remaining methods achieve near-linear or linear computational efficiency. A higher rejection rate indicates better performance in distinguishing between distributions. Overall, quadratic-time quantile-based estimators perform comparably to quadratic-time MMD estimators, while near-linear time quantile-based estimators often outperform their MMD-based counterparts.
Figure 6.4 : Comparing the time (in seconds) required to complete the CIFAR-10 vs. CIFAR-10.1 experiment, plotted on a logarithmic scale. A shorter time indicates a faster algorithm. These results align with our complexity analysis.
Appendix figures & tables10 assets
Supplementary material from the paper’s appendix.
Appendix
| Cost | Error | |
|---|---|---|
| V-statistic | ||
| OW |
Table A.1 : Computational and sample complexity rates of the V-statistic and the OW estimator with respect to .
| hyperparameter | value | |
|---|---|---|
| 0.05 | level of the test | |
| 200 | number of bootstrap samples | |
| 100 | number of samples from the simulator | |
| 500 | number of observations in the data | |
| 50 | number of initial parameters sampled | |
| 10 | number of initial parameters to optimise |
Table A.2 : Definitions of the hyperparameters.
| Variables | Mean | Std | Meaning |
|---|---|---|---|
| 1000 | 1.0 | Cost of treatment | |
| 0.1 | 0.02 | Probability of admissions | |
| 5.2 | 1.0 | Days of hospital | |
| 400 | 200 | Cost per day | |
| 0.3 | 0.1 | Utility change if response | |
| 3.0 | 0.5 | Duration of response |
Table B.1 : Variables in the health economics experiment.
Figure C.1 : Rates of decay for the ML and CV estimators for the Matérn kernel of order , and a true function that is a linear combination of Matérn kernels of order . The ML rate is given in Karvonen et al. [2020, Equation 5.2] . The CV rate is empirically observed in Figure C.2 . Observe that the CV estimator’s range of adaptation to the smoothness is wider than the ML estimator’s.
Figure C.2 : Asymptotics of CV estimator compared to asymptotics of the ML estimator, for the Matérn kernel of order , and a true function that is a finite linear combination of Matérn kernels of order .
Figure D.1 : Type I control results for our experiment on CIFAR-10 vs. CIFAR-10.1. We see all methods control their Type I error around or below the specified Type I error rate , thus confirming our tests in the main text are valid testing procedures.
Figure D.2 : The experiments in Figure 6.3 repeated for . Experimental results comparing our proposed methods with baseline approaches. A higher rejection rate indicates better performance in distinguishing between distributions. Same as for , quadratic-time quantile-based estimators perform comparably to quadratic-time MMD estimators, while near-linear time quantile-based estimators often outperform their MMD-based counterparts.
Figure D.3 : Gaussian KQD test power under different weighting measures. Left, middle: Varying measure : down-weighing ( ) extremes boosts power on Galaxy MNIST, while up-weighing ( ) them helps on CIFAR. Uniform weighting remains a strong default, with optimal depending on the dataset. Right: Varying measure : using an IQR-scaled Gaussian or uniform default reference measure both outperform MMD, indicating potential advantage of a "default" over the problem-based .
Figure D.4 : All methods are cost unless specified otherwise. Left: Gaussian KQD compared with sliced Wasserstein with uniform or data-driven directions, on the power decay problem. Sliced Wasserstein fall well below KQD, consistent with their equivalence to KQD using the less expressive linear kernel. Middle: Comparison with alternative approximate KME methods, at matching cost. ME matches MMD-multi power, while Nyström-MMD suffers high Type II error. Right: Comparison with Median-of-Means (MOM) KME approximation, at matching cost. MOM is primarily a robustness-enforcing method, not a cheap-approximation method, and does not perform well at set cost of .