Scalable Kernel-Based Distances for Statistical Inference and Integration

Masha Naslidnyk

stat.MLFeb 25, 2026

Scalable Kernel-Based Distances for Statistical Inference and Integration

Authors: Masha Naslidnyk

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

$\mathbb{N}$	The set of natural numbers, $\{0,1,\dots\}$ .
$\mathbb{N}_{\geq a}$	The subset $\{a,a+1,\dots\}\subseteq\mathbb{N}$ , for $a\in\mathbb{N}$ .
$\mathrm{Id}_{d}$	The $d\times d$ identity matrix.
$x_{1:N}\in\mathcal{X}^{N}$	Tuple $(x_{1},x_{2},\dots,x_{N})$ , for $x_{1},\dots,x_{N}\in\mathcal{X}$ .
$x_{1:N}\in\mathbb{R}^{N\times d}$	Vector (for $d=1$ ) or matrix $\begin{bmatrix}x_{1}&x_{2}&\dots&x_{N}\end{bmatrix}^{\top}$ , for $x_{1},\dots,x_{N}\in\mathcal{X}\subseteq\mathbb{R}^{d}$ .
$\\|\cdot\\|_{\mathcal{F}}$ , $\langle\cdot,\cdot\rangle_{\mathcal{F}}$	The norm and inner product of a space $\mathcal{F}$ .
$f(x_{1:N})\in\mathbb{R}^{N}$	Column vector $\begin{bmatrix}f(x_{1})&f(x_{2})&\dots&f(x_{N})\end{bmatrix}^{\top}$ , for $x_{1},\dots,x_{N}\in\mathcal{X}$ and $f:\mathcal{X}\to\mathbb{R}$ .

Table 2.1 : Summary of notation used in the thesis.

Polynomial of degree $q$	$k_{\mathrm{Poly}}(x,x^{\prime})=\tau^{2}(x^{\top}x^{\prime}+c)^{q}$ , for $c\geq 0$ and $q\in\mathbb{N}$ .
Matérn of order $\nu$ , and smoothness $s=\nu+d/2$	$k_{\nu}(x,x^{\prime})=\frac{\tau^{2}}{\Gamma(\nu)2^{\nu-1}}(\frac{\sqrt{2\nu}}{l}\\|x-x^{\prime}\\|_{2})^{\nu}K_{\nu}(\frac{\sqrt{2\nu}}{l}\\|x-x^{\prime}\\|_{2})$ , where $K_{\nu}$ for $\nu>0$ is the modified Bessel function of the second kind.
Gaussian	$k_{\mathrm{Gauss}}(x,x^{\prime})=\tau^{2}\exp(-\\|x-x^{\prime}\\|^{2}_{2}/(2l^{2}))$
Brownian motion	$k_{\mathrm{BM}}(x,x^{\prime})=\tau^{2}\min(x,x^{\prime})$
Fractional Brownian motion	$k_{0,H}(x,x^{\prime})=\tau^{2}\big(x^{2H}+{x^{\prime}}^{2H}-\lvert x-x^{\prime}\rvert^{2H}\big)/2$

Table 2.2 : Example kernels on

\mathcal{X}=\mathbb{R}^{d}

. For Brownian motion and fractional Brownian motion,

\mathcal{X}=[0,\infty)

Figure 3.1 : Estimating the MMD requires approximating the embedding $\mu_{k,\mathbb{P}_{\theta}}$ of the model $\mathbb{P}_{\theta}$ in the RKHS $\mathcal{H}_{k}$ . The classical approach approximates it using $N$ equally-weighted i.i.d. samples from $\mathbb{P}_{\theta}$ , denoted $\mu_{k,\mathbb{P}_{\theta,N}}$ . We show that this estimator can be improved by using optimally-weighted samples, denoted $\mu_{k,\mathbb{P}^{w}_{\theta,N}}$ .

Model	$r$	$d$	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

{}^{2}\penalty\ (\times 10^{-3})

between

\mathbb{P}_{\theta,N}

and

\mathbb{P}_{\theta,M}

over 100 runs for the V-statistic and our optimally-weighted (OW) estimator. Settings:

M=10000

N=256

Panel of Figure 3 — 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 $k$ and $c$ . Increasing the smoothness of $k$ improves the performance. (b) Comparison of V-statistic and OW estimator as a function of dimension. OW performs better for both parameterisations of $\mathbb{U}$ , with the Gaussian giving lowest error. (c) Value of $\theta_{4}$ also impacts the performance of the OW estimator. (d) Error vs. total computation cost for different $M$ . OW performs better than V-statistic for similar cost: $N=M$ for V-statistic, whereas $N=(68,126,200,317)$ for OW.

Cases	i.i.d. V-stat	i.i.d. OW (ours)
$\theta_{4}=0.1$ (null holds)	0.040	0.047
$\theta_{4}=0.5$ (alternative holds)	0.040	0.413

Table 3.2 : Fraction of repeats for which the null was rejected. An ideal test would have

0.05

when the null holds, and

1

otherwise.

Figure 3.3 : ABC posteriors for the wind farm model. Our OW estimator yields posterior samples that are more concentrated around the true

\theta_{0}

than the V-statistic. Settings:

M=100

\theta_{0}=20

Figure 4.1 : Illustration of conditional Bayesian quadrature (CBQ). Stage one fits a GP to $f(x,\theta_{t})$ for each $\theta_{t}\in\{\theta_{1},\dots,\theta_{T}\}$ and integrates to obtain BQ estimates $I_{\mathrm{BQ}}(\theta_{1}),\dots,I_{\mathrm{BQ}}(\theta_{T})$ . Stage two places a GP over $\theta\mapsto I(\theta)$ and fuses these to yield $I_{\mathrm{CBQ}}(\theta)$ with posterior uncertainty shown by shaded regions.

Panel of Figure 4 — Figure 4.2 : Bayesian sensitivity analysis for linear models. Left: RMSE of all methods when $d=2$ and $N=50$ . Middle: RMSE of all methods when $d=2$ and $T=50$ . Right: RMSE of all methods when $N=T=100$ .

Figure 4.3 : Bayesian linear model sensitivity analysis in

d=2

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

f(x)

with Hurst parameter

H=0.2

(smoothness

s+\alpha=0.2

), using the Brownian motion kernel ( 5.3 ) with amplitudes

\tau^{2}=1

(left),

\tau^{2}=\hat{\tau}_{\textup{CV}}^{2}=6.12

given by the LOO-CV estimator (middle), and

\tau^{2}=\hat{\tau}_{\textup{ML}}^{2}=3.7

given by the ML estimator (right). Vertical lines mark the true value of the integral

\int_{0}^{T}f(x)\text{d}x

(solid line) and the BQ posterior mean

I_{\mathrm{BQ}}=\int_{0}^{T}m_{N}(x)\text{d}x

(dashed line); the shaded bands show the 95% BQ credible interval

[I_{\mathrm{BQ}}-1.96\tau\sigma_{\mathrm{BQ}},I_{\mathrm{BQ}}+1.96\tau\sigma_{\mathrm{BQ}}]

. Bottom row: The GP interpolation underlying the top row: an fBm path that is the true integrand

f(x)

(solid line), the training data

x_{1:N},f(x_{1:N})

, and the GP posterior mean

m_{N}(x)

(dashed line); the shaded bands show the 95% GP credible interval

[m_{N}(x)-1.96\tau\sqrt{k_{N}(x,x)},m_{N}(x)+1.96\tau\sqrt{k_{N}(x,x)}]

Figure 5.2 : BQ of an integrated fractional Brownian motion (fBm) integrand

f(x)

with Hurst parameter

H=0.5

(smoothness

s+\alpha=1.5

), using the Brownian motion kernel ( 5.3 ) with amplitudes

\tau^{2}=1

(left),

\tau^{2}=\hat{\tau}_{\textup{CV}}^{2}=0.028

given by the LOO-CV estimator (middle), and

\tau^{2}=\hat{\tau}_{\textup{ML}}^{2}=0.11

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

s+\alpha

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

s+H

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 $s[f]$ and $\alpha[f]$ in ( 5.22 ). Runs on individual 100 samples from $f$ 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

\mathbb{P}

. Center: Density of the projection of

\mathbb{P}

onto direction

u

on the unit circle, with

\phi_{u}(x)=\langle u,x\rangle

. Right: Different quantiles for all possible directions

u

Figure 6.2 : Illustration of the impact of the slicing direction on KQEs. Suppose

X\sim\mathbb{P}

, the KQEs

\rho_{\mathbb{P}}^{\alpha,u}(x)\coloneqq\rho_{u\#\mathbb{P}}^{\alpha}u(x)

are obtained by considering the

\alpha

-th quantile of

u(X)

. 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	$\mathcal{O}(N^{2})$	$\mathcal{O}(N^{-\nicefrac{{1}}{{2}}})$
OW	$\mathcal{O}(N^{3})$	$\mathcal{O}(N^{-\nicefrac{{s_{c}}}{{r}}})$

Appendix

Table A.1 : Computational and sample complexity rates of the V-statistic and the OW estimator with respect to

N

Figure 24

hyperparameter	value
$\alpha$	0.05	level of the test
$B$	200	number of bootstrap samples
$N$	100	number of samples from the simulator
$M$	500	number of observations in the data
$I$	50	number of initial parameters sampled
$R$	10	number of initial parameters to optimise

Appendix

Table A.2 : Definitions of the hyperparameters.

Variables	Mean	Std	Meaning
$X_{1}$	1000	1.0	Cost of treatment
$X_{2}$	0.1	0.02	Probability of admissions
$X_{3}$	5.2	1.0	Days of hospital
$X_{4}$	400	200	Cost per day
$X_{5}$	0.3	0.1	Utility change if response
$X_{6}$	3.0	0.5	Duration of response

Appendix

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

1

, and a true function that is a linear combination of Matérn kernels of order

\nu_{\mathrm{true}}

. 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

\nu_{\mathrm{true}}

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

\nu_{\mathrm{model}}

, and a true function that is a finite linear combination of Matérn kernels of order

\nu_{\mathrm{true}}

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

0.05

, thus confirming our tests in the main text are valid testing procedures.

Figure D.2 : The experiments in Figure 6.3 repeated for

p=1

. Experimental results comparing our proposed methods with baseline approaches. A higher rejection rate indicates better performance in distinguishing between distributions. Same as for

p=2

, 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

\nu

: down-weighing (

/\setminus

) extremes boosts power on Galaxy MNIST, while up-weighing (

\setminus/

) them helps on CIFAR. Uniform weighting remains a strong default, with optimal

\nu

depending on the dataset. Right: Varying measure

\xi

: using an IQR-scaled Gaussian or uniform default reference measure

\xi

both outperform MMD, indicating potential advantage of a "default"

\xi

over the problem-based

\xi=(\mathbb{P}_{N}+\mathbb{Q}_{N})/2

Figure D.4 : All methods are cost

\mathcal{O}(N\log^{2}N)

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

\mathcal{O}(N\log^{2}N)

Numerical Integration

See all 24 papers

Jul 27, 2026

Minimax Lower Bounds of Kernel Discrepancy Estimation: MMD, HSIC, KSD

Over the past 20 years, kernel discrepancies have been leveraged as a highly powerful tool for quantifying the disagreement of distributions, with numerous successful applications in two-sample, goodness-of-fit, and independence testing, among others. Their fastest estimators are known to converge at a parametric rate---

n^{-1/2}

---under mild conditions. While this rate is known to be minimax optimal on

\mathbb R^d

under strict assumptions with bounded kernels, little is known about its optimality beyond the finite-dimensional Euclidean setting with unbounded kernels. In this work, we prove that the minimax lower bound of estimation of the most popular kernel discrepancies (maximum mean discrepancy, Hilbert-Schmidt independence criterion and kernel Stein discrepancy; MMD, HSIC, KSD) is

n^{-1/2}

on general topological spaces, and under mild assumptions on the kernel; the same rates are shown (as corollaries) to hold for the estimation of the mean embedding and the centered cross-covariance operator. Our results settle the question of optimal estimation of these kernel discrepancies.

Jose Cribeiro-Ramallo, Florian Kalinke, Zoltán Szabó

Chair of Information Systems, Karlsruhe Institute of Technology, Am Fasanengarten 5, 76131 Karlsruhe, Germany · Department of Statistics, London School of Economics, Houghton Street, London, WC2A 2AE, UK

Maximum Mean Discrepancy Minimax Estimation

Ask PDF

May 4, 2026

Measuring Differences between Conditional Distributions using Kernel Embeddings

Comparing conditional distributions is a fundamental challenge in statistics and machine learning, with applications across a wide range of domains. While proposed methods for measuring discrepancies using kernel embeddings of distributions in a reproducing kernel Hilbert space (RKHS) provide powerful non-parametric techniques, the existing literature remains fragmented and lacks a unified theoretical treatment. This paper addresses this gap by establishing a coherent framework for studying kernel-based methods to measure divergence between conditional distributions through what we refer to as conditional maximum mean discrepancy (CMMD). The CMMD consists of a family of metrics which we call levels, with three special cases each using a different type of RKHS embedding: CMMD

_0

(conditional mean operators), CMMD

_1

(conditional mean embeddings), and CMMD

_2

(joint mean embeddings). We additionally introduce a general level

s

CMMD, clarifying the required assumptions, and establishing mathematical connections between the levels through the lens of operator-based smoothing. In addition to reviewing previously proposed estimators, we introduce a novel doubly robust estimator for the CMMD that maintains consistency provided the product of errors of the underlying models reduces sufficiently fast. We provide numerical experiments demonstrating that the CMMD effectively captures complex conditional dependencies for statistical testing.

Peter Moskvichev, Siu Lun Chau, Dino Sejdinovic

School of Mathematical Sciences, Adelaide University · College of Computing and Data Science, Nanyang Technological University

Kernel Mean Embeddings Maximum Mean Discrepancy Kernel Methods

Ask PDF

Jul 22, 2026

Directional Kernel Mean Difference: A Fast Signed Statistic for Univariate Distribution Comparison

We introduce the Directional Kernel Mean Difference (DKMD), a signed statistic for univariate distribution comparison that preserves the direction of distributional shifts. Unlike the squared Maximum Mean Discrepancy (MMD), which discards directional information by squaring the RKHS distance, DKMD integrates the difference of kernel mean embeddings against a fixed odd weighting function. This construction yields three structural properties: antisymmetry, immunity to symmetric distributional differences, and directional monotonicity under stochastic dominance. We derive a data-driven Riemann estimator that ensures asymptotic consistency with the continuous formulation, strictly preserving the theoretical guarantees of the signed statistic in empirical evaluations. To overcome the quadratic computational cost of kernel methods, we develop an

O(N \log N)

prefix--suffix scanning algorithm that exploits the total order of the real line while requiring only

O(N)

memory. Experiments on synthetic benchmarks demonstrate that DKMD correctly isolates directional shifts from symmetric perturbations, remains robust to heavy-tailed outliers that can flip the sign of the mean difference, and scales to millions of samples in seconds.

Shijie Zhong, Jiangfeng Fu

School of Power and Energy Northwestern Polytechnical University Xi’an, Shanxi 710129

Two-Sample Testing Maximum Mean Discrepancy

Ask PDF

Abstract

Numerical Integration

Explore similar work

Minimax Lower Bounds of Kernel Discrepancy Estimation: MMD, HSIC, KSD

Measuring Differences between Conditional Distributions using Kernel Embeddings

Directional Kernel Mean Difference: A Fast Signed Statistic for Univariate Distribution Comparison