mmdshort = MMD, long = maximum mean discrepancy qmcshort = QMC, long = quasi
Monte Carlo pdfshort = PDF, long = probability density function ksdshort =
KSD, long = kernel Stein discrepancy mcmcshort = MCMC, long = Markov chain
Monte Carlo smcshort = SMC, long = sequential Monte Carlo

# Minimum Kernel Discrepancy Estimators

Chris. J. Oates
Newcastle University, UK

###### Abstract

For two decades, reproducing kernels and their associated discrepancies have
facilitated elegant theoretical analyses in the setting of quasi Monte Carlo.
These same tools are now receiving interest in statistics and related fields,
as criteria that can be used to select an appropriate statistical model for a
given dataset. The focus of this article is on _minimum kernel discrepancy
estimators_ , whose use in statistical applications is reviewed, and a general
theoretical framework for establishing their asymptotic properties is
presented.

##  1 Introduction

The abstract problem tackled in quasi Monte Carlo is to construct an empirical
distribution $P_{n}$ that provides an accurate approximation to a given
probability distribution $P$ of interest. To make this problem well-defined a
_discrepancy_ is required, to precisely quantify the quality of the
approximation $P_{n}$ to $P$. The mathematical tractability of the discrepancy
is an important consideration, both for the ease of algorithmic design and to
enable the mathematical determination of how close any particular
approximation $P_{n}$ is to being optimal. In a seminal paper, Hickernell
(1998) advocated _kernel discrepancies_ in this context, explaining how
reproducing kernels can be used to construct explicit discrepancies with a
range of favourable mathematical properties. Subsequent researchers have
exploited kernel discrepancies to analyse a variety of modern and powerful
quasi Monte Carlo methods (see e.g. Dick and Pillichshammer, 2010), and more
recently also Markov chain Monte Carlo methods (see e.g. Gorham and Mackey,
2017).

The abstract problem tackled in statistics is to construct a probabilistic
model $P$ for a population quantity of interest, given a finite sample from
this population with empirical distribution $P_{n}$. Thus, on the face of it,
the problems tacked in quasi Monte Carlo and in statistics are diametrically
opposed. Nevertheless, an equally valid approach to the statistical problem is
to first select a discrepancy and proceed to select a probabilistic model
$P_{\theta}$ from a parametric set such that the discrepancy between $P_{n}$
and $P_{\theta}$ is minimised. It is well-known that maximum likelihood
estimation asymptotically minimises the Kullback–Leibler divergence between
$P$ and $P_{\theta}$ (Akaike, 1973), and that this approach can be non-robust
when the statistical model is misspecified, meaning that $P$ and $P_{\theta}$
do not coincide for any value $\theta$ in the parameter set. Several robust
alternatives to maximum likelihood have been proposed, such as M-estimation
(Huber, 1964) and the more general minimum distance estimation (Donoho and
Liu, 1988), which remains an active research field (Basu et al., 2011; Pardo,
2018). More interestingly, as far as this article is concerned, _minimum
kernel discrepancy estimation_ has recently been proposed and studied. The use
of kernel discrepancies in the statistical context can be traced back at least
to Song et al. (2008), before receiving considerable attention in the machine
learning community, where minimum kernel discrepancy estimation has been used
to train machine learning models (Dziugaite et al., 2015; Li et al., 2015;
Sutherland et al., 2017), to calibrate computer simulations (Park et al.,
2016; Mitrovic et al., 2016; Dellaporta et al., 2022), and to facilitate
goodness-of-fit testing (Key et al., 2021). It seems, to this author, that
there could be a useful dialogue between the quasi Monte Carlo and the
statistical communities, centring on whether the theoretical insight into
kernel discrepancies gleaned in the former can be used to explain the
remarkable success in machine learning applications achieved by the latter.
This article serves as an introduction to minimum kernel discrepancy
estimation, and along the way some opportunities for dialogue between the
communities will be highlighted.

In Section 2 both kernel discrepancy and of minimum kernel discrepancy
estimators are defined. The use of minimum kernel discrepancy estimators in
statistics is described in Section 3. A concise set of theoretical results for
minimum kernel discrepancy estimation are presented in Section 4. Some open
questions and possible research directions are suggested in Section 5.

##  2 Minimum Kernel Discrepancy Estimators

Section 2.1 introduces the set-up and notation that will be used. Then,
minimum kernel discrepancy estimators are defined in Section 2.2.

###  2.1 Set-Up and Notation

Our setting is a measurable space $\mathcal{X}$. A _kernel_ is a function
$k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$, that is _measurable_ ,
_symmetric_ and _positive semi-definite_ , meaning respectively that

  * •

$k(\cdot,x)$ is measurable for all $x\in\mathcal{X}$

  * •

$k(x,y)=k(y,x)$ for all $x,y\in\mathcal{X}$

  * •

$\sum_{i=1}^{n}\sum_{j=1}^{n}w_{i}w_{j}k(x_{i},x_{j})\geq 0$ for all
$w_{1},\dots,w_{n}\in\mathbb{R}$, $x_{1},\dots,x_{n}\in\mathcal{X}$,
$n\in\mathbb{N}$ .

The _reproducing kernel Hilbert space_ associated to a kernel $k$ is a Hilbert
space $\mathcal{H}(k)$ of real-valued functions on $\mathcal{X}$, such that

  * •

$k(\cdot,x)\in\mathcal{H}(k)$ for all $x\in\mathcal{X}$

  * •

$\langle h,k(\cdot,x)\rangle_{\mathcal{H}(k)}=h(x)$ for all $x\in\mathcal{X}$,
$h\in\mathcal{H}(k)$ .

The latter requirement is known as the _reproducing property_. It can be shown
that $\mathcal{H}(k)$ exists, that $\mathcal{H}(k)$ is uniquely determined,
and that the elements of $\mathcal{H}(k)$ are measurable functions on
$\mathcal{X}$ (see Chapter 4 of Steinwart and Christmann, 2008). Our interest
in reproducing kernel Hilbert spaces derives from their suitability for
constructing a discrepancy between probability distributions on $\mathcal{X}$,
as explained next.

Let $\mathcal{P}_{k}(\mathcal{X})$ be the set of probability distributions $P$
on $\mathcal{X}$ for which $\mathrm{P}:\mathcal{H}(k)\rightarrow\mathbb{R}$,
$\mathrm{P}(h)=\int h\;\mathrm{d}P$, is a bounded linear functional.
Sufficient conditions for $P\in\mathcal{P}_{k}(\mathcal{X})$ are discussed in
Appendix A. If $P\in\mathcal{P}_{k}(\mathcal{X})$ then there exists a Riesz
representer $\mu_{k}(P)\in\mathcal{H}(k)$, called the _kernel mean element_ ,
such that $\mathrm{P}f=\langle f,\mu_{k}(P)\rangle_{\mathcal{H}(k)}$. The
_kernel discrepancy_ on $\mathcal{P}_{k}(\mathcal{X})$ defined as

$\displaystyle
D_{k}:\mathcal{P}_{k}(\mathcal{X})\times\mathcal{P}_{k}(\mathcal{X})$
$\displaystyle\rightarrow[0,\infty)$ $\displaystyle(P,Q)$
$\displaystyle\mapsto\|\mu_{k}(P)-\mu_{k}(Q)\|_{\mathcal{H}(k)},$

being the distance between the Riesz representers $\mu_{k}(P)$ and
$\mu_{k}(Q)$ in $\mathcal{H}(K)$. It is immediate that $D_{k}$ is a pseudo-
metric on $\mathcal{P}_{k}(\mathcal{X})$, meaning that $D_{k}$ satisfies all
the requirements of a metric on $\mathcal{P}_{k}(\mathcal{X})$ except for the
_distinguishability_ requirement, since in general we may have $D_{k}(P,Q)=0$
even if $P$ and $Q$ are distinct. The pseudo-metric $D_{k}$ appears in several
fields; in machine learning it is called the _maximum mean discrepancy_
(Gretton et al., 2012; Muandet et al., 2017), in probability it is called an
_integral probability (pseudo-)metric_ (Müller, 1997), in statistics it is
called a _kernel scoring rule_ (Gneiting and Raftery, 2007; Dawid, 2007), and
in quasi Monte Carlo it is called the _worst case integration error_ (Dick and
Pillichshammer, 2010). The emergence of the kernel discrepancy as a powerful
tool in these fields is due in large part to an explicit formula that enables
it to be computed:

######  Proposition 1.

For $P,Q\in\mathcal{P}_{k}(\mathcal{X})$,

$\displaystyle D_{k}(P,Q)=\sqrt{\iint
k_{Q}(x,y)\mathrm{d}P(x)\mathrm{d}P(y)},$

where

$\displaystyle k_{Q}(x,y):=k(x,y)-\int k(x,y)\mathrm{d}Q(x)-\int
k(x,y)\mathrm{d}Q(y)+\iint k(x,y)\mathrm{d}Q(x)\mathrm{d}Q(y).$ (1)

###### Proof.

Taking $h(x)=k(y,x)$ and using the reproducing property, together with the
definition of $\mu_{k}(P)$ as the Riesz representer of $\mathrm{P}$, leads to
$\mu_{k}(P)(y)=\langle\mu_{k}(P),k(\cdot,y)\rangle_{\mathcal{H}(k)}=\mathrm{P}(k(\cdot,y))=\int
k(x,y)\mathrm{d}P(x)$. Thus the kernel mean element is $\mu_{k}(P)(\cdot)=\int
k(x,\cdot)\mathrm{d}P(x)$, where all such integrals in this article are
understood as a Bochner integrals, for all $P\in\mathcal{P}_{k}(\mathcal{X})$.
In particular,
$\langle\mu_{k}(P),\mu_{k}(Q)\rangle_{\mathcal{H}(k)}=\int\mu_{k}(P)(y)\mathrm{d}Q(y)=\iint
k(x,y)\mathrm{d}P(x)\mathrm{d}Q(y)$ for all
$P,Q\in\mathcal{P}_{k}(\mathcal{X})$. Therefore, recalling that $k$ is
symmetric,

$\displaystyle D_{k}(P,Q)^{2}$
$\displaystyle=\langle\mu_{k}(P),\mu_{k}(P)\rangle_{\mathcal{H}(k)}-2\langle\mu_{k}(P),\mu_{k}(Q)\rangle_{\mathcal{H}(k)}+\langle\mu_{k}(Q),\mu_{k}(Q)\rangle_{\mathcal{H}(k)}$
$\displaystyle=\iint k(x,y)\mathrm{d}P(x)\mathrm{d}P(x)-2\iint
k(x,y)\mathrm{d}P(x)\mathrm{d}Q(y)+\iint k(x,y)\mathrm{d}Q(x)\mathrm{d}Q(y)$
$\displaystyle=\iint k_{Q}(x,y)\mathrm{d}P(x)\mathrm{d}P(y),$

which completes the proof. ∎

It can be verified that the map
$k_{Q}:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ in (1) is a
$Q$-dependent kernel, satisfying $\int k_{Q}(x,y)\mathrm{d}P(x)=0$ for all
$y\in\mathcal{X}$ whenever $Q=P$. This presentation, which assigns different
roles to $P$ and $Q$, is useful in our setting, where the first argument $P$
will be fixed and the second argument $Q$ will be varied.

###  2.2 Minimum Kernel Discrepancy Estimators

Consider data that are assumed to arise as a sequence of independent and
identically distributed random variables $(x_{n})_{n\in\mathbb{N}}$, with
$x_{n}\sim P$ and $P\in\mathcal{P}_{k}(\mathcal{X})$. Let
$P_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$ denote the empirical
distribution of the first $n$ terms in the dataset. A _statistical model_ is a
collection
$\\{P_{\theta}\\}_{\theta\in\Theta}\subset\mathcal{P}_{k}(\mathcal{X})$, where
the _parameter_ $\theta$ takes values in an index set $\Theta$. _Minimum
kernel discrepancy estimators_ , if they exist, are defined as

$\displaystyle\theta_{n}\in\operatorname*{argmin}_{\theta\in\Theta}D_{k}(P_{n},P_{\theta}),$

as illustrated in Figure 1. The properties of such estimators depend on the
choice of kernel $k$ and the statistical model
$\\{P_{\theta}\\}_{\theta\in\Theta}$, but are independent of how the
statistical model is parametrised (in $\theta$). Under regularity conditions,
established in Section 4, the estimators $\theta_{n}$ converge (in an
appropriate sense) to a population minimiser

$\displaystyle\theta_{\star}$
$\displaystyle\in\operatorname*{argmin}_{\theta\in\Theta}D_{k}(P,P_{\theta}),$

if these exist. Since in all cases we take $Q=P_{\theta}$ in (1), we
henceforth adopt the shorthand $k_{\theta}$ for $k_{Q}$. Through selection of
the kernel $k$, one can trade statistical efficiency with robustness to model
misspecification, which is not possible within the classical framework of
maximum likelihood; see Section 4.5. However, important aspects of kernel
selection remains unsolved, for example it is unclear how to retain
statistical efficiency when $\mathcal{X}$ is high-dimensional. Such open
challenges are highlighted in Section 5.

$\\{\mu_{k}(P_{\theta})\\}_{\theta\in\Theta}$$\mu_{k}(P_{n})$$\mu_{k}(P_{\theta_{n}})$$\mu_{k}(P)$$\mu_{k}(P_{\theta_{\star}})$$P_{\theta_{n}}$$P_{n}$$P$$\mathcal{H}(k)$
Figure 1: Minimum kernel discrepancy estimation. The distance, as measured in
$\mathcal{H}(k)$, between the kernel mean element $\mu_{k}(P_{\theta})$
associated to the statistical model $P_{\theta}$ and the kernel mean element
$\mu_{k}(P_{n})$ associated to the dataset is minimised. Under conditions
established in Section 4, a minimum kernel discrepancy estimator $\theta_{n}$
converges to a value $\theta_{\star}$ that minimises the distance to the
kernel mean element $\mu_{k}(P)$ associated to the true data-generating
distribution $P$.

##  3 Applications

This section presents three distinct applications for minimum kernel
discrepancy estimation: the generalised method of moments (Section 3.1),
generative adversarial networks (Section 3.2), and energy-based models
(Section 3.3).

###  3.1 Generalised Method of Moments

Let $\phi:\mathcal{X}\rightarrow\mathbb{R}^{p}$ have components that are
integrable with respect to each $\\{P_{\theta}\\}_{\theta\in\Theta}$, and let
$\phi(P):=\int\phi\mathrm{d}P_{\theta}$ denote a vector of _summary
statistics_ , meaning that this vector summarises some of the salient aspects
of $P_{\theta}$. The _generalised method of moments_ estimator, if it exists,
is defined as

$\displaystyle\theta_{n}\in\operatorname*{argmin}_{\theta\in\Theta}\|\phi(P_{n})-\phi(P_{\theta})\|,$

being the parameter vector of a statistical model whose summary statistics
agree most closely with the empirical summary statistics $\phi(P_{n})$
calculated using the dataset (Hansen, 1982). This procedure can be cast as
minimum kernel discrepancy estimation using a kernel
$k(x,y)=\langle\phi(x),\phi(y)\rangle$. For this kernel, $\int
k(x,y)\mathrm{d}P_{\theta}(y)=\langle\phi(x),\phi(P_{\theta})\rangle$ and
$\iint
k(x,y)\mathrm{d}P_{\theta}(x)\mathrm{d}P_{\theta}(y)=\langle\phi(P_{\theta}),\phi(P_{\theta})\rangle$.
Then
$k_{\theta}(x,y)=\langle\phi(x),\phi(y)\rangle-\langle\phi(x),\phi(P_{\theta})\rangle-\langle\phi(P_{\theta}),\phi(y)\rangle+\langle\phi(P_{\theta}),\phi(P_{\theta})\rangle$,
so that
$D_{k}(P_{n},P_{\theta})^{2}=\langle\phi(P_{n}),\phi(P_{n})\rangle-\langle\phi(P_{n}),\phi(P_{\theta})\rangle-\langle\phi(P_{\theta}),\phi(P_{n})\rangle+\langle\phi(P_{\theta}),\phi(P_{\theta})\rangle=\|\phi(P_{n})-\phi(P_{\theta})\|^{2}$.
In contrast to minimum kernel discrepancy estimation in general, the
generalised method of moments is mathematically transparent in the special
case where we assume the statistical model $P_{\theta}$ can be parametrised
such that $\theta=\phi(P_{\theta})$ and $\Theta=\mathbb{R}^{p}$. Indeed, in
this case $D_{k}(P_{n},P_{\theta})=\|\theta_{n}-\theta\|$, so the generalised
method of moments estimator is

$\theta_{n}=\frac{1}{n}\sum_{i=1}^{n}\phi(x_{i}),$

which converges to $\theta_{\star}=\phi(P)$ and satisfies the central limit
$\sqrt{n}(\theta_{n}-\theta_{\star})\operatorname*{\stackrel{{\scriptstyle\textup{d}}}{{\longrightarrow}}}N(0,\mathbb{C}_{X\sim
P}[\phi(X)])$ whenever the associated covariance matrix is well-defined. This
will provide a convenient opportunity to verify, in a concrete setting, the
general theoretical results presented in Section 4.3.

Though originating as a statistically efficient and computationally favourable
alternative to maximum likelihood, it has more recently been observed that,
through careful selection of the summary statistics $\phi$, one can usefully
trade-off the efficiency and the robustness of the estimator. The related
areas of _approximate Bayesian computation_ (Beaumont, 2019) and _Bayesian
synthetic likelihood_ (e.g. Frazier and Drovandi, 2021) contain extensive
practical guidance on how summary statistics can be selected. Indeed, the
estimator $\theta_{n}$ can be recognised as a _maximum a posteriori_ estimator
in an approxmate Bayesian computation framework when a flat _a priori_
distribution is employed.

###  3.2 Generative Adversarial Networks

A _generative model_ consists of a measurable space $\mathcal{Y}$, equipped
with a probability measure $\mathbb{P}$, and statistical model of the form

$\displaystyle P_{\theta}=G_{\\#}^{\theta}\;\mathbb{P}$ (2)

where $G^{\theta}:\mathcal{Y}\rightarrow\mathcal{X}$ is measurable for each
$\theta\in\Theta$. Here $G_{\\#}^{\theta}\;\mathbb{P}$ denotes the
_pushforward_ of $\mathbb{P}$ through $G^{\theta}$, meaning that samples
$X\sim P_{\theta}$ are generated by first sampling $Y\sim\mathbb{P}$ and then
setting $X=G^{\theta}(Y)$. In the specific case where $G^{\theta}$ is a (deep)
neural network, we call $P_{\theta}$ a (deep) _neural network generative
model_. It is clear that many rich and interesting statistical models can be
constructed in this manner, using arbitrarily large neural networks with a
corresponding parameter vector $\theta$ that is arbitrarily high-dimensional.
However, conventional statistical wisdom suggests that such models may not be
useful in practice, due to the prohibitively large number of parameters that
will need to be estimated. Maximum likelihood estimation has shown to struggle
when real data are used since, despite its flexibility, the statistical model
will necessarily be misspecified (Theis et al., 2016). Alternative strategies
have nevertheless been developed to estimate the parameters $\theta$ of such
models and – remarkably – these have demonstrated a high level of success
provided one has access to a sufficiently rich training dataset. Here we focus
on a line of research called _generative adversarial networks_ (GANs)
(Goodfellow et al., 2020), and specifically the so-called _integral
probability metric_ GANs. Let $P_{n}$ denote the empirical distribution of the
training dataset. Then a GAN estimator for $\theta$, if one exists, is defined
as

$\displaystyle\theta_{n}\in\operatorname*{argmin}_{\theta\in\Theta}\;\max_{f\in\mathcal{F}}\;\left|\int
f\mathrm{d}P_{n}-\int f\mathrm{d}P_{\theta}\right|,$

where $\mathcal{F}$ is a set of test functions such that $\mathcal{F}\subset
L^{1}(P_{\theta})$ for each $\theta\in\Theta$. The choice of $\mathcal{F}$
determines the properties of the GAN estimator. The GAN nomenclature
originates from the case where the elements of $\mathcal{F}$ are also neural
networks, so that the _witness function_ or _critic_ $f$ for which the maximum
is realised can be interpreted as an _adversarial_ neural network. Alternative
choices of $\mathcal{F}$ lead to so-called _generative moment matching
networks_ (Dziugaite et al., 2015; Li et al., 2015), _MMD_ GANs (Li et al.,
2017), _Wasserstein_ GANs (Arjovsky et al., 2017), _Mc_ GANs (Mroueh et al.,
2017), _Fisher_ GANs (Mroueh and Sercu, 2017), and _Sobolev_ GANs (Mroueh et
al., 2018). Here we focus on generative moment matching networks, which take
$\mathcal{F}$ to be the unit ball in a reproducing kernel Hilbert space
$\mathcal{H}(k)$ of real-valued functions on $\mathcal{X}$; see Figure 2. This
procedure can be cast as minimum kernel discrepancy estimation:

######  Proposition 2.

The GAN estimator $\theta_{n}$ is a minimum kernel discrepancy estimator when
$\mathcal{F}$ is the set of $f\in\mathcal{H}(k)$ such that
$\|f\|_{\mathcal{H}(k)}\leq 1$.

###### Proof.

For $f\in\mathcal{F}$, using Cauchy–Schwarz,

$\displaystyle\left|\int f\mathrm{d}P_{n}-\int f\mathrm{d}P_{\theta}\right|$
$\displaystyle=\left|\langle
f,\mu_{k}(P_{n})-\mu_{k}(P_{\theta})\rangle\right|\leq\|f\|_{\mathcal{H}(k)}\|\mu_{k}(P_{n})-\mu_{k}(P_{\theta})\|_{\mathcal{H}(k)}$

which shows that

$\displaystyle 0\leq\max_{f\in\mathcal{F}}\;\int f\mathrm{d}P_{n}-\int
f\mathrm{d}P_{\theta}\leq D_{k}(P_{n},P_{\theta}).$

If $D_{k}(P_{n},P_{\theta})=0$ then we must have equality, so suppose
$D_{k}(P_{n},P_{\theta})\neq 0$. Then one can verify that
$f=[\mu_{k}(P_{n})-\mu_{k}(P_{\theta})]/D_{k}(P_{n},P_{\theta})$ is the
witness function for which equality is realised. ∎

This class of GANs is appealing in the sense that the optimisation over
$\mathcal{F}$ can be explicitly solved, and the solution is the kernel
discrepancy $D_{k}(P_{n},P_{\theta})$. Although the kernel mean element
$\mu_{k}(P_{\theta})$ is not closed-form in general, it can be consistently
approximated using an empirical approximation to $P_{\theta}$; this is
discussed further in Section 5. The performance of these GANs depends on the
choice of kernel $k$, and methods have been proposed to flexibly _learn_ a
suitable $k$ based on the training dataset (Li et al., 2017). These methods
have demonstrated remarkable performance on machine learning benchmarks, out-
performing other forms of GAN (Bińkowski et al., 2018).

$\Theta$$\theta$$\theta_{n}$$C_{n}^{0.95}$$P_{\theta}$$P_{\theta_{n}}$ Figure
2: Minimum kernel discrepancy estimation for generative models. The parameters
$\theta$ of sophisticated generative models may be estimated by minimising the
kernel discrepancy between the generative model $P_{\theta}$ and the empirical
distribution of the dataset. Confidence sets $C_{n}^{\gamma}$ for
$\theta_{\star}$ can be constructed to ensure the notional coverage of
$100\gamma\%$ is asymptotically correct; see Section 4.

###  3.3 Energy-Based Models

An _energy-based model_ on $\mathcal{X}=\mathbb{R}^{d}$ is a statistical model
$P_{\theta}$ with pdf

$\displaystyle p_{\theta}(x)=\frac{\exp(-E_{\theta}(x))}{Z(\theta)}$

where $E_{\theta}:\mathcal{X}\rightarrow\mathbb{R}$ is a flexibly parametrised
_potential_ such that $x\mapsto\exp(-E_{\theta}(x))$ is integrable and
$Z(\theta):=\int\exp(-E_{\theta}(x))\mathrm{d}x>0$. Energy-based models arise
both in statistics, where they are referred to as an instance of _intractable
likelihood_ (Lyne et al., 2015), and in machine learning, where they provide
an alternative to GANs that can be more convenient when an explicit pdf is
required (LeCun et al., 2007). Similarly to GANs, maximum likelihood
estimation has been shown to struggle in this setting, where, despite its
flexibility, the statistical model will necessarily be misspecified. However,
there is an additional barrier to application of minimum kernel discrepancy
estimation to energy-based models; the challenge of approximating the kernel
mean element $\mu_{k}(P_{\theta})$. Indeed, due to the intractable
normalisation constant $Z(\theta)$, sampling from an energy-based model is not
straight-forward and an algorithm such as mcmc is required to approximately
sample from $P_{\theta}$ for each value of $\theta\in\Theta$. This can entail
a prohibitive computational cost when attempting optimisation over $\Theta$
(Song and Kingma, 2021).

To address this computational challenge, we can employ a class of kernels $k$
introduced in Oates et al. (2017). These kernels $k$ solve the computational
problem by ensuring that $\mu_{k}(P_{\theta})$ is equal to 0 in
$\mathcal{H}(k)$, and thus $\mu_{k}(P_{\theta})$ can be trivially computed. To
achieve this, the kernel $k$ must be allowed to depend on $\theta$, and in
what follows we write $k(x,y;\theta)$ to make this dependence explicit. Before
describing the kernels with this property, we first introduce some notation
and assumptions. Let $\nabla
f=(\partial_{x_{1}}f,\dots,\partial_{x_{d}}f)^{\top}$ denote the gradient of a
differentiable function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ and let
$\nabla\cdot f=\partial_{x_{1}}f_{1}+\dots+\partial_{x_{d}}f_{d}$ denote the
divergence of a differentiable function
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$. The presentation in the remainder
of this section follows Barp et al. (2022):

######  Proposition 3.

Let $g:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ satisfy
$g,\mathcal{A}_{\theta}g\in L^{1}(P_{\theta})$ where

$(\mathcal{A}_{\theta}g)(x):=(\nabla\cdot g)(x)-\langle g(x),(\nabla
E_{\theta})(x)\rangle.$

Then $\int\mathcal{A}_{\theta}g\;\mathrm{d}P_{\theta}=0$.

###### Proof.

Let $\varphi_{n}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a compactly supported
function with $\varphi_{n}(x)=1$ for $\|x\|\leq n$ and
$\|\nabla\varphi_{n}\|<n^{-1}$, for each $n\in\mathbb{N}$. From the divergence
theorem,

$\displaystyle
0=\int(\nabla\cdot(\varphi_{n}p_{\theta}g))(x)\;\mathrm{d}x=\int\langle\nabla\varphi_{n}(x),(p_{\theta}g)(x)\rangle\;\mathrm{d}x+\int\varphi_{n}(x)(\nabla\cdot(p_{\theta}g))(x)\;\mathrm{d}x.$

Since $\varphi_{n}\rightarrow 1$ pointwise and $\nabla\cdot(p_{\theta}g)\in
L^{1}(\mathbb{R}^{d})$, from the dominated convergence theorem

$\displaystyle\int\varphi_{n}(x)(\nabla\cdot(p_{\theta}g))(x)\;\mathrm{d}x\rightarrow\int(\nabla\cdot(p_{\theta}g))(x)\;\mathrm{d}x.$

On the other hand, since $p_{\theta}g\in L^{1}(\mathbb{R}^{d})$,

$\displaystyle\left|\int\langle\nabla\varphi_{n}(x),(p_{\theta}g)(x)\rangle\;\mathrm{d}x\right|\leq\|\nabla\varphi_{n}\|\int\|(p_{\theta}g)(x)\|\;\mathrm{d}x\rightarrow
0.$

Thus we have shown that

$\displaystyle\int\mathcal{A}_{\theta}g\;\mathrm{d}P_{\theta}=\int\frac{1}{p_{\theta}}\nabla\cdot(p_{\theta}g)\;\mathrm{d}P_{\theta}=\int(\nabla\cdot(p_{\theta}g))(x)\;\mathrm{d}x=0,$

as claimed. ∎

######  Proposition 4.

Let $c:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a kernel
with
$c(\cdot,y),\nabla_{y}c(\cdot,y),\mathcal{A}_{\theta}c(\cdot,y)\mathrm{1}_{d},\mathcal{A}_{\theta}\nabla_{y}c(\cdot,y)\in
L^{1}(P_{\theta})$ for each $y\in\mathbb{R}^{d}$. Then

$\displaystyle k(x,y;\theta)$
$\displaystyle=\nabla_{x}\cdot\nabla_{y}c(x,y)-\langle\nabla_{x}E_{\theta}(x),\nabla_{y}c(x,y)\rangle$
$\displaystyle\qquad-\langle\nabla_{y}E_{\theta}(y),\nabla_{x}c(x,y)\rangle+\langle\nabla_{x}E_{\theta}(x),\nabla_{y}E_{\theta}(y)\rangle
c(x,y)$

is a $\theta$-dependent kernel with $\mu_{k}(P_{\theta})=0$ for all
$\theta\in\Theta$.

###### Proof.

For the proof that $k(\cdot,\cdot;\theta)$ is indeed a kernel we refer to Barp
et al. (2022, Theorem 2.6). For the final part, note that

$\displaystyle k(x,y;\theta)=\mathcal{A}_{\theta}g(x),\qquad
g(x)=\nabla_{y}c(x,y)+c(x,y)\nabla_{y}\log p_{\theta}(y)$

where, under our assumptions, $g,\mathcal{A}_{\theta}g\in L^{1}(P_{\theta})$.
Then we may apply Proposition 3 to deduce that $\mu_{k}(P_{\theta})(y)=\int
k(x,y;\theta)\mathrm{d}P_{\theta}(x)=\int\mathcal{A}_{\theta}g(x)\mathrm{d}P_{\theta}(x)=0$
for all $y\in\mathcal{X}$, and thus $\mu_{k}(P_{\theta})=0$, as claimed. ∎

This construction enables minimum kernel discrepancy estimation to be applied
to energy-based models, with positive initial result reported in a range of
statistical applications (Barp et al., 2019; Matsubara et al., 2022a). Indeed,
$k(\cdot,\cdot;\theta)$ can be evaluated without the normalisation constant
$Z(\theta)$, and any kernel for which the conclusion of Proposition 4 holds
satisfies $k_{\theta}(\cdot,\cdot)=k(\cdot,\cdot;\theta)$. The discrepancy
associated to $k(\cdot,\cdot;\theta)$ is called a _kernel Stein discrepancy_
(Liu et al., 2016; Chwialkowski et al., 2016; Gorham and Mackey, 2017); these
kernelise the _score-matching divergence_ of Hyvärinen and Dayan (2005) and
have recently been applied to a variety of other tasks in statistics and
machine learning (Anastasiou et al., 2022).

This completes our short tour of how minimum kernel discrepancy methods can be
used. Next we turn our attention to the asymptotic properties of these
estimators, aiming insofar as possible for a general framework.

##  4 Asymptotic Theory

The aim of this section is to present general asymptotic theory that is
applicable to each of the applications we have just discussed. Conditions for
strong consistency of the minimum kernel discrepancy estimator are established
for a $\theta$-independent kernel in Section 4.1 and extended to the case of a
$\theta$-dependent kernel in Section 4.2. Conditions for asymptotic normality
in the case of a possibly $\theta$-dependent kernel are established in Section
4.3, and the strength of these conditions is examined in Section 4.4. The
presentation is intended to be self-contained; we relate our results to
existing literature in Section 4.5.

To simplify the presentation, it will be assumed throughout that
$P\in\mathcal{P}_{k}(\mathcal{X})$ and
$\\{P_{\theta}\\}_{\theta\in\Theta}\subset\mathcal{P}_{k}(\mathcal{X})$, so
that the kernel discrepancy is well-defined. It will also be assumed that the
statistical model is parameterised in such a way that $\theta_{\star}$ and
$\theta_{n}$ actually exist, in the latter case for $n$ sufficiently large
(but in neither case is uniqueness assumed). Convergence in distribution,
convergence in probability, and almost sure convergence will be respectively
denoted
$\operatorname*{\stackrel{{\scriptstyle\textup{d}}}{{\longrightarrow}}}$,
$\operatorname*{\stackrel{{\scriptstyle\textup{p}}}{{\longrightarrow}}}$ and
$\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}$.

###  4.1 Strong Consistency

In this section we consider a topological space $\Theta$; no additional
structure on $\Theta$ is required until Section 4.2. Our first task is to find
weak sufficient conditions under which a.s. every sequence
$(\theta_{n})_{n\in\mathbb{N}}$ of minimum kernel discrepancy estimators
satisfies

$\displaystyle\\{\text{accumulation points
of}\;(\theta_{n})_{n\in\mathbb{N}}\\}\subseteq\operatorname*{argmin}_{\theta\in\Theta}D_{k}(P,P_{\theta}).$
(3)

That is, there is a probability 1 set on which all accumulation points of all
sequences $\theta_{n}$ of minimum kernel discrepancy estimators are minimisers
$\theta_{\star}$ of the kernel discrepancy between the distribution $P$, from
which the data are sampled, and the statistical model $P_{\theta}$. To begin
with we consider the case where the kernel $k$ does not depend on the
parameter $\theta$; the case of a $\theta$-dependent kernel is treated in
Section 4.2.

######  Lemma 5 (Strong consistency of the kernel mean element).

Assume that

  1. A1

$\mathcal{H}(k)$ is separable.

  2. A2

$\int\sqrt{k(x,x)}\;\mathrm{d}P(x)<\infty$.

Then
$\mu_{k}(P_{n})\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\mu_{k}(P)$.

###### Proof.

The random variables $\xi_{i}=k(\cdot,x_{i})$ are independent, identically
distributed, and satisfy
$\mathbb{E}[\|k(\cdot,x_{i})\|_{\mathcal{H}(k)}]=\int\sqrt{k(x,x)}\;\mathrm{d}P<\infty$
(from A2), so we may appeal to the strong law of large numbers (Theorem 18) in
the separable (from A1) Banach space $\mathcal{B}=\mathcal{H}(k)$, to get

$\displaystyle\mu_{k}(P_{n})=\frac{1}{n}\sum_{i=1}^{n}k(\cdot,x_{i})\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\int
k(\cdot,x)\mathrm{d}P(x)=\mu_{k}(P),$

which establishes the result. ∎

######  Lemma 6 (Uniform convergence in discrepancy).

Assume A1, A2. Then

$\displaystyle\sup_{\theta\in\Theta}|D_{k}(P_{n},P_{\theta})-D_{k}(P,P_{\theta})|$
$\displaystyle\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}0.$
(4)

###### Proof.

Let $f(\theta)=D_{k}(P,P_{\theta})$ and
$f_{n}(\theta)=D_{k}(P_{n},P_{\theta})$. From the reverse triangle inequality
and Lemma 5, we have that

$\displaystyle\sup_{\theta\in\Theta}|f_{n}(\theta)-f(\theta)|$
$\displaystyle=\sup_{\theta\in\Theta}\left|\|\mu_{k}(P_{n})-\mu_{k}(P_{\theta})\|_{\mathcal{H}(k)}-\|\mu_{k}(P)-\mu_{k}(P_{\theta})\|_{\mathcal{H}(k)}\right|$
$\displaystyle\leq\|\mu_{k}(P_{n})-\mu_{k}(P)\|_{\mathcal{H}(k)}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}0,$
(5)

which establishes the result. ∎

######  Lemma 7 (Discrepancy is minimised).

Assume A1, A2. Then

$D_{k}(P,P_{\theta_{n}})\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\min_{\theta\in\Theta}D_{k}(P,P_{\theta}).$

###### Proof.

Let $f$ and $f_{n}$ be as in the proof of Lemma 6. Let $\epsilon>0$. Pick one
element $\theta_{\star}\in\operatorname*{argmin}_{\theta\in\Theta}f(\theta)$.
Then a.s. there exists $n_{0}$ such that for all $n>n_{0}$ we have that

$f(\theta_{\star})\stackrel{{\scriptstyle\text{defn of
}\theta_{\star}}}{{\leq}}f(\theta_{n})\stackrel{{\scriptstyle\text{\lx@cref{creftypecap~refnum}{lem:
uniform}}}}{{<}}f_{n}(\theta_{n})+\frac{\epsilon}{2}\stackrel{{\scriptstyle\text{defn
of
}\theta_{n}}}{{\leq}}f_{n}(\theta_{\star})+\frac{\epsilon}{2}\stackrel{{\scriptstyle\text{\lx@cref{creftypecap~refnum}{lem:
uniform}}}}{{<}}f(\theta_{\star})+\epsilon.$

Since $\epsilon>0$ was arbitrary, it follows that
$f(\theta_{n})\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\min_{\theta\in\Theta}f(\theta)$.
∎

Using an additional assumption that there exists a unique and sufficiently
distinct global minimum at the parameter level, we obtain a basic strong
consistency result at the parameter level. The proof of the following result
is immediate from Lemma 7:

######  Theorem 8 (Strong consistency with unique minimum).

Assume A1, A2. If $\theta_{\star}\in\Theta$ is such that for every open
neighbourhood $\theta_{\star}\in N(\theta_{\star})\subset\Theta$ we have

$\displaystyle D_{k}(P,P_{\theta_{\star}})<\inf_{\theta\in\Theta\setminus
N(\theta_{\star})}D_{k}(P,P_{\theta}),$ (6)

then
$\theta_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\theta_{\star}$.

In practice the uniqueness of a minimum is often difficult or impossible to
check, since in general the statistical model will be misspecified. This
motivates a more realistic analysis, which we state in two parts. The
_closure_ of a subset $S$ of a topological space $\Theta$ will be denoted
$\text{cl}(S)$, and we recall that a subset $C$ of a topological space
$\Theta$ is _sequentially compact_ if every sequence in $C$ has an
accumulation point in $C$.

######  Theorem 9 (Strong consistency with several minima, part I).

Assume A1, A2, and assume that $\theta\mapsto D_{k}(P,P_{\theta})$ is a
continuous function on $\Theta$. Then any accumulation point of
$(\theta_{n})_{n\in\mathbb{N}}$ is a.s. an element of
$\min_{\theta\in\Theta}D_{k}(P,P_{\theta})$.

###### Proof.

Let $f$ be as in the proof of Lemma 6. Let $\epsilon>0$ and let
$\theta_{\star}$ be an accumulation point, meaning that there exists a
subsequence with $\theta_{n_{m}}\rightarrow\theta_{\star}$ as
$m\rightarrow\infty$. Then a.s. there exists $m_{0}$ such that for all
$m>m_{0}$ we have

$f(\theta_{\star})\stackrel{{\scriptstyle f\text{ cts at
}\theta_{\star}}}{{<}}f(\theta_{n_{m}})+\epsilon.$

Taking the limit $m\rightarrow\infty$ on both sides and using Lemma 7, we have
that a.s.

$f(\theta_{\star})<\min_{\theta\in\Theta}f(\theta)+\epsilon.$

Since $\epsilon>0$ was arbitrary, it follows that a.s.
$\theta_{\star}\in\min_{\theta\in\Theta}f(\theta)$. ∎

######  Theorem 10 (Strong consistency with several minima, part II).

Assume A1, A2, and assume there is an $\epsilon>0$ for which

$C:=\text{cl}\left\\{\theta\in\Theta:D_{k}(P,P_{\theta})<\epsilon+\inf_{\vartheta\in\Theta}D_{k}(P,P_{\vartheta})\right\\}$

is sequentially compact. Then the sequence $(\theta_{n})_{n\in\mathbb{N}}$ has
at least one accumulation point.

###### Proof.

Let $f$ and $f_{n}$ be as in the proof of Lemma 6. From Lemma 6, there a.s.
exists $n_{0}$ such that for all $n>n_{0}$ we have
$\sup_{\theta\in\Theta}|f_{n}(\theta)-f(\theta)|<\epsilon/2$. In particular,
there a.s. exists $n_{0}$ such that $(\theta_{n})_{n>n_{0}}\subset C$. Since
the sequence $(\theta_{n})_{n\in\mathbb{N}}$ is a.s. eventually contained in
$C$, and $C$ is sequentially compact, there a.s. exists an accumulation point.
∎

For $\Theta=\mathbb{R}^{d}$, the topological condition in Theorem 10 is
satisfied when $C$ is bounded.

###  4.2 Extension to $\theta$-Dependent Kernel

To accommodate also the $\theta$-dependent kernels $k(\cdot,\cdot;\theta)$
from Section 3.3, here we provide an alternative to Lemma 6.

######  Lemma 11 (Uniform convergence in discrepancy; $\theta$-dependent
kernel).

Assume A1 for each $\theta\in\Theta$. Assume $\Theta\subset\mathbb{R}^{p}$ is
open, convex and bounded, and that

  1. A1

$\int\sqrt{k_{\theta}(x,x)}\;\mathrm{d}P(x)<\infty$, for all $\theta\in\Theta$

  2. A2

$\int\sup_{\theta\in\Theta}\left\|\partial_{\theta}k_{\theta}(x,x)\right\|\;\mathrm{d}P(x)<\infty$

  3. A3

$\iint\sup_{\theta\in\Theta}\left\|\partial_{\theta}k_{\theta}(x,y)\right\|\;\mathrm{d}P(x)\mathrm{d}P(y)<\infty$
.

Then

$\displaystyle\sup_{\theta\in\Theta}|D_{\theta}(P_{n},P_{\theta})-D_{\theta}(P,P_{\theta})|$
$\displaystyle\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}0.$
(7)

###### Proof.

The aim is to establish the conditions for the uniform law of large numbers in
Theorem 21. Let $g(\theta)=D_{\theta}(P,P_{\theta})^{2}$ and
$g_{n}(\theta)=D_{\theta}(P_{n},P_{\theta})^{2}$. From (1) and the fact that
$k(\cdot,\cdot;\theta)$ is positive semi-definite,

$\displaystyle k(x,x;\theta)$ $\displaystyle\leq k_{\theta}(x,x)+2\int
k(x,y;\theta)\mathrm{d}P_{\theta}(y).$

Taking the square root and integrating,

$\displaystyle\int\sqrt{k(x,x;\theta)}\mathrm{d}P(x)$
$\displaystyle\leq\int\sqrt{k_{\theta}(x,x)+2\int
k(x,y;\theta)\mathrm{d}P_{\theta}(y)}\;\mathrm{d}P(x)$
$\displaystyle\leq\int\sqrt{k_{\theta}(x,x)}\mathrm{d}P(x)+\int\sqrt{2\int
k(x,y;\theta)\mathrm{d}P_{\theta}(y)}\;\mathrm{d}P(x)$ (8)
$\displaystyle\leq\int\sqrt{k_{\theta}(x,x)}\mathrm{d}P(x)+\sqrt{\int 2\int
k(x,y;\theta)\mathrm{d}P_{\theta}(y)\;\mathrm{d}P(x)},$ (9)

where (8) used the fact that $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$, and (9) used
Jensen’s inequality. The first integral in (9) exists from A1 and all other
integrals exist as a consequence of the standing assumption that
$P,P_{\theta}\in\mathcal{P}_{k(\cdot,\cdot;\theta)}(\mathcal{X})$. Thus
$\int\sqrt{k(x,x;\theta)}\;\mathrm{d}P(x)<\infty$. From this fact and A1, we
can apply Lemma 5 to deduce that
$g_{n}(\theta)\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}g(\theta)$
for each $\theta\in\Theta$. Thus the first condition in Theorem 21 is
satisfied.

For the second condition in Theorem 21, we start by noting it is implicit in
A2, A3 that $\partial_{\theta}k_{\theta}$ exists in $\Theta$ and thus $g_{n}$
is differentiable in $\Theta$. Let $\theta,\vartheta\in\Theta$. From the mean
value theorem, and the fact that $\Theta$ is open and convex,

$\displaystyle g_{n}(\theta)-g_{n}(\vartheta)$
$\displaystyle=\langle\theta-\vartheta,\partial_{\phi}g_{n}(\phi)\rangle$

for some $\phi\in\\{t\theta+(1-t)\vartheta:0\leq t\leq 1\\}\subset\Theta$. In
particular,

$\displaystyle\left|g_{n}(\theta)-g_{n}(\vartheta)\right|$
$\displaystyle\leq\|\theta-\vartheta\|\;\sup_{\phi\in\Theta}\|\partial_{\phi}g_{n}(\phi)\|\leq\|\theta-\vartheta\|\;\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\underbrace{\sup_{\phi\in\Theta}\|\partial_{\phi}k_{\phi}(x_{i},x_{j})\|}_{=:v(x_{i},x_{j})}.$

Now, since we have a V-statistic, we use A2 and A3 to verify that

$\displaystyle\iint|v(x,y)|\mathrm{d}P(x)\mathrm{d}P(y)$
$\displaystyle=\iint\sup_{\theta\in\Theta}\|\partial_{\theta}k_{\theta}(x,y)\|\mathrm{d}P(x)\mathrm{d}P(y)<\infty$
$\displaystyle\int|v(x,x)|\mathrm{d}P(x)$
$\displaystyle=\int\sup_{\theta\in\Theta}\|\partial_{\theta}k_{\theta}(x,x)\|\mathrm{d}P(x)<\infty$

and thus from Theorem 19 we have that

$\displaystyle B_{n}$
$\displaystyle:=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\sup_{\theta\in\Theta}\|\partial_{\theta}k_{\theta}(x_{i},x_{j})\|\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\iint\sup_{\theta\in\Theta}\|\partial_{\theta}k_{\theta}(x,y)\|\mathrm{d}P(x)\mathrm{d}P(y)<\infty.$

Since in addition $\Theta$ is assumed to be bounded, the conditions of the
uniform law of large numbers in Theorem 21 have now been established. Thus
Theorem 21 gives that
$\sup_{\theta\in\Theta}|g_{n}(\theta)-g(\theta)|\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}0$,
as claimed. ∎

The arguments of Lemma 7 and Theorems 8, 9 and 10 all go through with A1, A2
and A3 in place of A2. The assumption that $\Theta$ is bounded may not be as
strong as it first appears, since one is free to decide how the statistical
model $P_{\theta}$ is parametrised (we do not assume $\Theta$ is compact).

###  4.3 Asymptotic Normality

The last asymptotic result we present is to show that fluctuations of the
minimum kernel discrepancy estimator are asymptotically normal. For a matrix
$M\in\mathbb{R}^{p\times p}$, we write $M\succ 0$ to indicate that $M$ is
positive definite; i.e. for all $0\neq c\in\mathbb{R}^{p}$ we have
$c^{\top}Mc>0$.

######  Theorem 12 (Asymptotic normality).

Suppose
$\theta_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\theta_{\star}$.
Let there exist an open set $S\subseteq\Theta\subseteq\mathbb{R}^{p}$ such
that $\theta_{\star}\in S$ and the following hold:

  1. A1

$\iint\sup_{\theta\in
S}\left\|\partial_{\theta}k_{\theta}(x,y)\right\|\mathrm{d}P(x)\mathrm{d}P(y)<\infty$

  2. A2

the functions
$\\{\theta\mapsto\partial_{\theta}k_{\theta}(x,y):x,y\in\mathcal{X}\\}$ are
differentiable on $S$

  3. A3

the functions
$\\{\theta\mapsto\partial_{\theta}^{2}k_{\theta}(x,y):x,y\in\mathcal{X}\\}$
are uniformly continuous at $\theta_{\star}$

  4. A4

$\iint\|\partial_{\theta}k_{\theta}(x,y)\|^{2}\mathrm{d}P(x)\mathrm{d}P(y)|_{\theta=\theta_{\star}}<\infty$

  5. A5

$\iint\|\partial_{\theta}^{2}k_{\theta}(x,x)\|\mathrm{d}P(x)|_{\theta=\theta_{\star}}<\infty$

  6. A6

$\iint\|\partial_{\theta}^{2}k_{\theta}(x,y)\|\mathrm{d}P(x)\mathrm{d}P(y)|_{\theta=\theta_{\star}}<\infty$

  7. A7

$\Gamma:=\frac{1}{2}\left.\iint\partial_{\theta}^{2}k_{\theta}(x,y)\mathrm{d}P(x)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}\succ
0$.

Then
$\sqrt{n}(\theta_{n}-\theta_{\star})\operatorname*{\stackrel{{\scriptstyle\textup{d}}}{{\longrightarrow}}}N(0,\Gamma^{-1}\Sigma\Gamma^{-\top})$,
where $\Sigma:=\mathbb{C}_{X\sim
P}\left[\int\partial_{\theta}k_{\theta}(X,y)\mathrm{d}P(y)\right]|_{\theta=\theta_{\star}}$.

###### Proof.

For a function $\theta\mapsto h(\theta)$, let
$h^{\prime}(\theta)=\partial_{\theta}h(\theta)$, and let
$h^{\prime\prime}(\theta)=\partial_{\theta}^{2}h(\theta)$. Let $g$ and $g_{n}$
be as in the proof of Lemma 11. From A2 there is a convex open set
$N(\theta_{\star})\subset\Theta$ with $\theta_{\star}\in N(\theta_{\star})$ on
which the functions $\theta\mapsto\partial_{\theta}k_{\theta}(x,y)$ are
differentiable, and thus the map $\theta\mapsto\partial_{\theta}g_{n}(\theta)$
is differentiable on $N(\theta_{\star})$. Since we are assuming
$\theta_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\theta_{\star}$,
there a.s. exists $n_{0}\in\mathbb{N}$ such that for all $n\geq n_{0}$ we have
$\theta_{n}\in N(\theta_{\star})$. Thus for each $n\geq n_{0}$ we can apply
the mean value theorem on the convex open set $N(\theta_{\star})$ to establish
that there is a $\tilde{\theta}_{n}=t_{n}\theta_{n}+(1-t_{n})\theta_{\star}$,
$t_{n}\in[0,1]$, for which

$\displaystyle
0=g_{n}^{\prime}(\theta_{n})=g_{n}^{\prime}(\theta_{\star})+g_{n}^{\prime\prime}(\tilde{\theta}_{n})\times(\theta_{n}-\theta_{\star}).$

If $g_{n}^{\prime\prime}(\tilde{\theta}_{n})$ is non-singular (as will be
established below by showing convergence to a matrix $2\Gamma$, assumed to be
positive definite due to A7), then we can re-express this as

$\displaystyle\sqrt{n}(\theta_{n}-\theta_{\star})=-[g_{n}^{\prime\prime}(\tilde{\theta}_{n})]^{-1}[\sqrt{n}g_{n}^{\prime}(\theta_{\star})]$
(10)

The stated result follows from Slutsky’s lemma if we can show both that
$g_{n}^{\prime\prime}(\tilde{\theta}_{n})\operatorname*{\stackrel{{\scriptstyle\textup{p}}}{{\longrightarrow}}}2\Gamma$
and
$\sqrt{n}g_{n}^{\prime}(\theta_{\star})\operatorname*{\stackrel{{\scriptstyle\textup{d}}}{{\longrightarrow}}}N(0,4\Sigma)$.
Both of these will now be established.

##### Establishing
$g_{n}^{\prime\prime}(\tilde{\theta}_{n}){\operatorname*{\stackrel{{\scriptstyle\textup{p}}}{{\longrightarrow}}}}2\Gamma$.

First we make the trivial algebraic observation that

$\displaystyle g_{n}^{\prime\prime}(\tilde{\theta}_{n})$
$\displaystyle=g_{n}^{\prime\prime}(\theta_{\star})+[g_{n}^{\prime\prime}(\tilde{\theta}_{n})-g_{n}^{\prime\prime}(\theta_{\star})].$
(11)

From A5 and A6, for $r,s\in\\{1,\dots,p\\}$ we have that

$\displaystyle\left.\int|\partial_{\theta_{r}}\partial_{\theta_{s}}k_{\theta}(x,x)|\mathrm{d}P(x)\right|_{\theta=\theta_{\star}},\quad\left.\iint|\partial_{\theta_{r}}\partial_{\theta_{s}}k_{\theta}(x,y)|\mathrm{d}P(x)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}<\infty$

and thus Theorem 19 on the consistency of V-statistics can be applied,
establishing that

$\displaystyle[g_{n}^{\prime\prime}(\theta_{\star})]_{r,s}$
$\displaystyle=\left.\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\partial_{\theta_{r}}\partial_{\theta_{s}}k_{\theta}(x_{i},x_{j})\right|_{\theta=\theta_{\star}}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\left.\iint\partial_{\theta_{r}}\partial_{\theta_{s}}k_{\theta}(x,y)\mathrm{d}P(x)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}=2\Gamma_{r,s},$

and in particular
$g_{n}^{\prime\prime}(\theta_{\star})\operatorname*{\stackrel{{\scriptstyle\textup{p}}}{{\longrightarrow}}}2\Gamma$.
Next consider the second term in (11). For any $\epsilon>0$, from A3 there
exists a bounded open set $N_{\epsilon}(\theta_{\star})\subset\Theta$ such
that $\theta_{\star}\in N_{\epsilon}(\theta_{\star})$ and
$\|\partial_{\theta}^{2}k_{\theta}(x,y)-\partial_{\theta}^{2}k_{\theta}(x,y)|_{\theta=\theta_{\star}}\|<\epsilon$
for all $\theta\in N_{\epsilon}(\theta_{\star})$ and all $x,y\in\mathcal{X}$.
Since we are assuming
$\theta_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\theta_{\star}$,
this means
$\tilde{\theta}_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\theta_{\star}$,
and thus there a.s. exists $n_{\epsilon}\in\mathbb{N}$ such that for all
$n\geq n_{\epsilon}$ we have $\tilde{\theta}_{n}\in
N_{\epsilon}(\theta_{\star})$. Thus, for all $n\geq n_{\epsilon}$,

$\displaystyle\left\|g_{n}^{\prime\prime}(\tilde{\theta}_{n})-g_{n}^{\prime\prime}(\theta_{\star})\right\|$
$\displaystyle=\left\|\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\partial_{\theta}^{2}k_{\theta}(x_{i},x_{j})|_{\theta=\tilde{\theta}_{n}}-\partial_{\theta}^{2}k_{\theta}(x_{i},x_{j})|_{\theta=\theta_{\star}}\right\|$
$\displaystyle\leq\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\left\|\partial_{\theta}^{2}k_{\theta}(x_{i},x_{j})|_{\theta=\tilde{\theta}_{n}}-\partial_{\theta}^{2}k_{\theta}(x_{i},x_{j})|_{\theta=\theta_{\star}}\right\|<\epsilon$

which demonstrates that the second term in (11) a.s. vanishes.

##### Establishing
$\sqrt{n}g_{n}^{\prime}(\theta_{\star}){\operatorname*{\stackrel{{\scriptstyle\textup{d}}}{{\longrightarrow}}}}N(0,4\Sigma)$.

Let $0\neq c\in\mathbb{R}^{p}$. From A1, the function
$h:S\times(\mathcal{X}\times\mathcal{X})\rightarrow\mathbb{R}^{p}$ defined as
$h(\theta,(x,y))=k_{\theta}(x,y)$ has $\partial_{\theta}h$ locally uniformly
integrably bounded in the sense of Definition 22, and thus the conditions of
Lemma 23 hold in a neighbourhood of $\theta_{\star}$, permitting us to
interchange $\partial_{\theta}$ and the integral with respect to $P\times P$.
This, together with the fact that $\theta_{\star}$ is a minimiser of
$g(\theta)$, shows that

$\displaystyle\left.\iint
c^{\top}\partial_{\theta}k_{\theta}(x,y)\mathrm{d}P(x)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}$
$\displaystyle=c^{\top}\partial_{\theta}\left.\iint
k_{\theta}(x,y)\mathrm{d}P(x)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}=c^{\top}g^{\prime}(\theta_{\star})=0.$
(12)

Further, from A4 and Cauchy–Schwarz,

$\displaystyle\left.\iint(c^{\top}\partial_{\theta}k_{\theta}(x,y))^{2}\mathrm{d}P(x)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}\leq\|c\|^{2}\left.\iint\|\partial_{\theta}k_{\theta}(x,y)\|^{2}\mathrm{d}P(x)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}<\infty.$

Thus we can appeal to asymptotic normality of V-statistics, using Theorem 20,
to obtain that

$\displaystyle
c^{\top}\sqrt{n}g_{n}^{\prime}(\theta_{\star})=\sqrt{n}\left(\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}c^{\top}\partial_{\theta}k_{\theta}(x_{i},x_{j})|_{\theta=\theta_{\star}}-\underbrace{0}_{\text{from
\eqref{eq: gprins is
0}}}\right)\operatorname*{\stackrel{{\scriptstyle\textup{d}}}{{\longrightarrow}}}N(0,4\sigma_{c}^{2})$

where

$\displaystyle\sigma_{c}^{2}$ $\displaystyle=\mathbb{V}_{X\sim
P}\left[\left.\int
c^{\top}\partial_{\theta}k_{\theta}(X,y)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}\right]=c^{\top}\Sigma
c.$

Since $c\neq 0$ was arbitrary, it follows that
$\sqrt{n}g_{n}^{\prime}(\theta_{\star})=\sqrt{n}[g_{n}^{\prime}(\theta_{\star})-g^{\prime}(\theta_{\star})]\operatorname*{\stackrel{{\scriptstyle\textup{d}}}{{\longrightarrow}}}N(0,4\Sigma)$,
as required. ∎

The assumptions of Theorem 12 will be discussed in Section 4.4, but we briefly
remark that A3 implies A1 and that, due to the crude way it appears in the
proof, A3 can likely be weakened. The conclusion of Theorem 12 agrees with our
discussion of the genralised method of moments in Section 3.1; in this case
$\partial_{\theta}k_{\theta}(x,y)=-\phi(x)-\phi(y)+2\theta$ and
$\partial_{\theta}^{2}k_{\theta}(x,y)=2I$, so $\Gamma=I$ and
$\Sigma=\mathbb{C}_{X\sim P}[\phi(X)]$, and Theorem 12 gives
$\sqrt{n}(\theta_{n}-\theta_{\star})\operatorname*{\stackrel{{\scriptstyle\textup{d}}}{{\longrightarrow}}}N(0,\mathbb{C}_{X\sim
P}[\phi(X)])$, in agreement. Finally, we note that Theorem 12 gives rise to a
consistent _sandwich_ estimator $\Gamma_{n}^{-1}\Sigma_{n}\Gamma_{n}^{-1}$ for
the asymptotic covariance $\Gamma^{-1}\Sigma\Gamma^{-1}$ of the minimum kernel
discrepancy estimator, where $\Gamma_{n}$ and $\Sigma_{n}$ are obtained by
substituting $\theta_{n}$ _in lieu_ of $\theta_{\star}$ into the definitions
of $\Gamma$ and $\Sigma$; see Freedman (2006). The resulting $100\gamma\%$
confidence sets for $\theta_{\star}$ are illustrated in Figure 2.

###  4.4 Discussion of Assumptions for the Exponential Family Model

To better understand A1-A7, we consider the consequences of these assumptions
in the setting of canonical exponential family statistical model

$p_{\theta}(x)=\exp(\langle\theta,t(x)\rangle-a(\theta)+b(x))$

where $x\in\mathcal{X}=\mathbb{R}^{d}$, whose parameter
$\theta\in\Theta\subset\mathbb{R}^{p}$ is estimated using either the
generalised method of moments (Section 4.4.1) or using a kernel Stein
discrepancy (Section 4.4.2).

####  4.4.1 Estimation via Generalised Method of Moments

Consider the generalised method of moments from Section 3.1, which uses a
finite rank kernel of the form $k(x,y)=\langle\phi(x),\phi(y)\rangle$. To
simplify the following discussion we focus on the well-specified case where
there is a unique $\theta_{\star}\in\Theta=\mathbb{R}^{p}$ with
$P=P_{\theta_{\star}}$ and seek conditions under which the generalised moment
matching estimator is consistent, with fluctuations that are asymptotically
normal. In this context there are essentially two important conditions that
must be satisfied:

  1. B1

$\phi$ is uniformly continuous and bounded

  2. B2

$M_{i,j}:=\left.\int\phi_{j}(x)[t_{i}(x)-\partial_{\theta_{i}}a(\theta)]\mathrm{d}P(x)\right|_{\theta=\theta_{\star}}$
has full row rank.

The first part, B1, is satisfied by (for example) the _Fourier features_
popularised in Rahimi and Recht (2007), while B2 ensures that the set of
features is rich enough to identify the correct model. In addition to these
two main conditions, we require sufficient regularity that, in a bounded open
neighbourhood $S$ of $\theta_{\star}$:

  1. B1

$\theta\mapsto\int\phi\;\mathrm{d}P_{\theta}$ is continuous at each $\theta\in
S$

  2. B2

$\theta\mapsto\int t\;\mathrm{d}P_{\theta}$ and $\int
tt^{\top}\;\mathrm{d}P_{\theta}$ are continuous at $\theta_{\star}$

  3. B3

$\sup_{\theta\in S}\int\|t\|^{2}\;\mathrm{d}P_{\theta}<\infty$

  4. B4

for each $\theta\in S$, there exists $\gamma_{\theta}>0$ such that
$\int(1+\|t\|^{2})\exp(\gamma_{\theta}\|t\|)\mathrm{d}P_{\theta}<\infty$

  5. B5

$\partial_{\theta}^{2}a(\theta)$ exists for each $\theta\in S$ and is
continuous at $\theta_{\star}$,

which are all reasonably mild. Collectively B1-B5 imply A1-A2 and A1-A7, so
that the generalised moment matching estimator is both consistent and
asymptotically normal. Full derivations are reserved for Section C.1.

####  4.4.2 Estimation via Kernel Stein Discrepancy

Consider the minimum kernel Stein discrepancy estimator from Section 3.3. Here
we suppose that $\Theta$ is an open, convex and bounded subset of
$\mathbb{R}^{p}$ but, in contrast to the last section, we do not assume here
the statistical model is well-specified. For
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, the convention $[\nabla
f]_{i,j}=\partial_{x_{i}}f_{j}(x)$ will be used. In this application, there
are again two main conditions to be satisfied:

  1. C1

$c(x,y),\nabla_{x}c(x,y),\nabla_{y}c(x,y)$ are bounded over
$x,y\in\mathbb{R}^{d}$

  2. C2

$\nabla E,\nabla t,\nabla b$ and $x\mapsto\nabla_{x}\cdot\nabla_{y}c(x,y)$,
$y\in\mathbb{R}^{d}$, are in $L^{1}(P_{\theta})$ for each $\theta\in\Theta$

  3. C3

$\iint c(x,y)[\nabla t(x)]^{\top}[\nabla
t(y)]\;\mathrm{d}P(x)\mathrm{d}P(y)\succ 0$

The requirements C1-C2 ensure the integrability conditions of Proposition 4
are satisfied, while C3 is analogous to B2, ensuring that the optimal
statistical model can be identified. In addition, the following regularity is
required:

  1. C1

$(x,y)\mapsto
c(x,y),\nabla_{x}c(x,y),\nabla_{y}c(x,y),\nabla_{x}\cdot\nabla_{y}c(x,y)$ are
continuous

  2. C2

$\nabla t$ and $\nabla b$ are continuous

  3. C3

$\int\|\nabla t\|^{4}\;\mathrm{d}P$, $\int\|\nabla t\|^{3}\|\nabla
b\|\;\mathrm{d}P$, $\int\|\nabla t\|^{2}\|\nabla b\|^{2}\;\mathrm{d}P$,
$\int\|\nabla b\|\;\mathrm{d}P<\infty$,

which are again reasonably mild. Collectively C1-C3 imply A1 and A1-A7, so
that the minimum kernel Stein discrepancy estimator is both consistent and
asymptotically normal. Full derivations are reserved for Section C.2

###  4.5 Related Work

The theoretical properties of minimum kernel discrepancy estimators have been
studied by a number of authors since 2015, but our analysis is arguably the
most general to have appeared. Indeed, as we will discuss, previous work
focused on either the case where $k$ does not depend on $\theta$, or on the
case where $k$ has a specific form of dependence on $\theta$.

In the case where $k$ does not depend on $\theta$, the original work of
Dziugaite et al. (2015) focused on obtaining non-asymptotic concentration
inequalities at the level of the kernel discrepancy. Although not discussed in
that work, these guarantees can be extended to guarantees at the level of the
parameter under further assumptions similar to those in Theorem 8. These non-
asymptotic results were extended to the regression context, where data
$x_{1},\dots,x_{n}$ are non-independent, in Alquier and Gerber (2020);
Chérief-Abdellatif and Alquier (2022). Both a non-asymptotic analysis and an
asymptotic analysis were presented in Briol et al. (2019), where strong
consistency was established under assumptions that include (1) a unique
$\theta_{\star}$, and (2) for all $n$, there exists $\epsilon_{n}>0$ such that
$\\{\theta\in\Theta:D_{k}(P_{n},P_{\theta})<D_{k}(P,P_{\theta_{\star}})+\epsilon_{n}\\}$
is a.s. bounded. The latter paper is comprehensive and includes also a
discussion of the important issues of estimator robustness, the geometry
induced by kernel discrepancy, and computational aspects of minimum kernel
discrepancy estimation (a point we return to in Section 5). A Bayesian
interpretation of minimum discrepancy estimation was proposed in Chérief-
Abdellatif and Alquier (2020).

The case where $k$ depends on $\theta$ was first studied in Barp et al.
(2019). There the authors focus on asymptotic analysis and a specific form of
$k$ derived from Stein’s method, the canonical form of which we discussed in
Section 3.3. Strong consistency of the minimum kernel Stein discrepancy
estimator was established under assumptions that include (1) the statistical
model is well-specified (i.e. $P=P_{\theta_{\star}}$), (2) $\theta_{\star}$ is
unique, (3) $\theta\mapsto P_{\theta}$ is injective, and either (4a) $\Theta$
is compact, or (4b) the maps $\theta\mapsto D_{k}(P_{n},P_{\theta})$ are a.s.
convex. Sufficient conditions for asymptotic normality are also provided, and
the important issues of estimator robustness, the geometry induced by kernel
discrepancy, and computational aspects of minimum kernel Stein discrepancy
estimation were discussed. Closest in spirit to our analysis is Matsubara et
al. (2022b), who endowed the minimum kernel Stein discrepancy estimator with a
Bayesian interpetation. The present analysis aims to be more concise, and also
more general in the senses that: (1) uniqueness of $\theta_{\star}$ is not
assumed, and (2) existence of third derivatives
$\partial_{\theta}^{3}k_{\theta}$ is not assumed. The aforementioned authors
also discussed estimator robustness and computational aspects of minimum
kernel Stein discrepancy estimation in detail.

Naturally, many alternatives to kernel discrepancy exist and have been
explored in the parameter inference context; to limit scope these have not
been discussed. However, we note that kernel discrepancies form a large class,
some of which can be viewed as the limit of entropy-regularised optimal
transport (Genevay et al., 2018), some of which offer control over smoothed
Wasserstein distances (Nietert et al., 2021), and some of which metrise weak
convergence on compact domains (Simon-Gabriel et al., 2020).

##  5 Open Research Directions

This paper was intended to serve both as a self-contained exposition of
minimum kernel discrepancy estimation and as an invitation to contribute to
this nascent research field. Accordingly, in this section we present a non-
exhaustive list of open research questions that could be addressed:

  1. 1.

Computing the kernel mean element. An important issue which we have not yet
discussed is how to compute the kernel mean element $\mu_{k}(P_{\theta})$,
given that for general choices of $k$ and $P_{\theta}$ there will not be a
closed-form solution to this integral. In the machine learning community, a
Monte Carlo approximation of $P_{\theta}$ is typically used, often in tandem
with stochastic gradient-based optimisation over $\theta\in\Theta$. However,
in some situations one might expect to do better: In the context of generative
models (2), one has $\mu_{k}(P_{\theta})=\int
k(G^{\theta}(x),\cdot)\mathrm{d}\mathbb{P}(x)$ where $\mathbb{P}$ is a user-
specified reference measure, typically the uniform masure on $[0,1]^{d}$. This
raises the possibility of using quasi Monte Carlo methods to approximate the
kernel mean element. A preliminary investigation in Niu et al. (2021)
established the $O(n^{-1+\epsilon})$ rate, but it remains an open problem to
obtain explicit conditions on $G^{\theta}$ and $k$ that ensure higher-order
quasi Monte Carlo convergence rates hold. Related, although the kernel mean
element $\mu_{k}(P_{n})$ is explicit, it is sometimes prohibitive to work with
a very large dataset. In these circumstances, methods have been proposed to
select a representative subset of data such that the kernel mean element
associated to this subset is close to that of the full datset (Cortes and
Scott, 2016; Teymur et al., 2021). However, the interaction between
computational approximations and the performance of the minimum kernel
discrepancy estimator have yet to be studied, beyond the Monte Carlo case in
Briol et al. (2019).

  2. 2.

Selecting a kernel. Beyond the basic requirement that the estimator is
consistent, further guidence on the application-specific selection of the
kernel $k$ is needed. The choice of $k$ engenders an efficiency-robustness
trade-off, so one must speculate about the extent to which a statistical model
is likely to be misspecified. Further, the efficiency of minimum kernel
discrepancy estimation in high-dimensional domains $\mathcal{X}$ is yet to be
explored. It is typical to encounter $P$ whose mass is concentrated around a
subset of $\mathcal{X}$, often a low-dimensional linear subspace or a sub-
manifold. In such circumstances it would be desirable to select a kernel $k$
that is appropriate adapted to this low-dimensional set, to improve the
dimension dependence of the estimator. Insight may be gained by considering
the case where $\mathcal{X}$ is infinite-dimensional; kernel discrepancies are
well-defined in this context (Wynne and Duncan, 2022; Wynne et al., 2022) but
minimum kernel discrepancy estimation is yet to be explored. On the other
hand, the case where $\mathcal{X}$ is discrete raises a different set of
challenges, since selection of a natural kernel in the discrete context can be
difficult. Both settings are further complicated by the fact that minimum
kernel discrepancy estimation is not invariant to how the data are
represented, i.e. the choice of $\mathcal{X}$ itself, and both theoretical and
practical guidence is needed.

  3. 3.

Insight into machine learning. Progress in machine learning is rapid, and in
2021 it was demonstrated that GANs are out-performed by so-called _diffusion
models_ on a range of image-generation tasks (Dhariwal and Nichol, 2021). Such
applications are characterised by large training datasets, so $P_{n}$ is
approximately equal to $P$, and high-dimensional $\mathcal{X}$, so there is a
loose analogy to the setting where quasi Monte Carlo is studied. One might
speculate on whether insights from quasi Monte Carlo might help to explain the
limited performance of GANs on such high-dimensional imaging tasks, or even
suggest ways in which GANs might be improved.

  4. 4.

Probability in Hilbert spaces. The strong consistency of the minimum kernel
discrepancy estimator was established, for a $\theta$-independent kernel $k$,
as a consequence of the strong law of large numbers in $\mathcal{H}(k)$. For a
$\theta$-independent kernel it may also be possible also to establish
asymptotic normality from the Hilbert space version of the central limit
theorem. This approach would arguably be more elegant than the one we
presented, avoiding the need to deal with V-statistics. However, it is not
clear to this author whether such a strategy can be applied in the case of a
$\theta$-dependent kernel.

  5. 5.

Confidence sets. An estimate $\theta_{n}$ for a parameter of interest
$\theta_{\star}$ is only useful if the precision of $\theta_{n}$ can be
estimated. Theorem 12 suggests the asymptotically valid $100\gamma\%$
confidence set

$C_{n}^{\gamma}:=\\{\theta\in\Theta:\|\Sigma_{n}^{-1/2}\Gamma_{n}(\theta-\theta_{n})\|^{2}\leq
F_{\chi_{p}^{2}}(\gamma)\\},$

where for simplicity here we assume that $\Theta=\mathbb{R}^{p}$ and
$\Sigma\succ 0$, and we use $F_{\cdot}$ notation to denote the cumulative
probability distribution, here of a chi-squared distribution with $p$ degrees
of freedom, $\chi_{p}^{2}$. However, at finite sample sizes $n$, the coverage
of these confidence sets could be far from the notional $100\gamma\%$ due to
the variance of the plug-in estimators $\Sigma_{n}$ and $\Gamma_{n}$ for
$\Sigma$ and $\Gamma$. This problem is exacerbated when $\theta$ is high-
dimensional, since these matrices are of dimension $p\times p$ where
$p=\text{dim}(\Theta)$. A variety of strategies for regularised estimation
exist and could be explored in this context (Lam, 2020).

##### Acknowledgements

The author wishes to thank Alessandro Barp, François-Xavier Briol, Andrew
Duncan, Jeremias Knoblauch, Lester Mackey and Takuo Matsubara for comments on
an earlier version of this manuscript, and financial support from EP/W019590/1
and EP/N510129/1.

## Appendices

Appendix A discusses sufficient conditions for
$P\in\mathcal{P}_{k}(\mathcal{X})$. Appendix B contains auxiliary results from
probability and calculus that we used. Appendix C contains full details for
the case of the canonical exponential family statistical model from Section
4.4 of the main text.

##  Appendix A Kernel Mean Embedding

In this section, $P$ is a distribution on a measurable space $\mathcal{X}$,
and $k$ is a measurable and symmetric positive definite kernel with
reproducing kernel Hilbert space denoted $\mathcal{H}(k)$. The central issue
considered here is when $\mathrm{P}(h)=\int h\;\mathrm{d}P$ is a bounded
linear functional on $\mathcal{H}(k)$. Our presentation here follows Barp et
al. (2022, Appendix C).

######  Definition 13 (Scalarly integrable).

A map $\Phi:\mathcal{X}\rightarrow\mathcal{H}(k)$ is said to be _scalarly_
$P$-_integrable_ if $\\{x\mapsto\langle
h,\Phi(x)\rangle_{\mathcal{H}(k)}:h\in\mathcal{H}(k)\\}\subset L^{1}(P)$.

The argument used in the following proof can be traced back to Dunford (1937),
with our account based on Schwabik and Ye (2005, Lemma 2.1.1).

######  Proposition 14.

If $\Phi:\mathcal{X}\rightarrow\mathcal{H}(k)$ is scalarly $P$-integrable,
then

$\displaystyle\mathrm{P}_{\Phi}:\mathcal{H}(k)$ $\displaystyle\rightarrow$
$\displaystyle\mathbb{R}$ $\displaystyle h$ $\displaystyle\mapsto$
$\displaystyle\int\langle h,\Phi(x)\rangle_{\mathcal{H}(k)}\;\mathrm{d}P(x)$

is a bounded linear functional.

###### Proof.

First we claim that the graph of the linear map

$\displaystyle T:\mathcal{H}(k)$ $\displaystyle\rightarrow$ $\displaystyle
L^{1}(P)$ $\displaystyle h$ $\displaystyle\mapsto$
$\displaystyle(x\mapsto\langle h,\Phi(x)\rangle_{\mathcal{H}(k)}).$

is closed. To see this, let $h_{n}\rightarrow h$ in $\mathcal{H}(k)$ and
suppose that $T(h_{n})\rightarrow g$ in $L^{1}(P)$. The claim is that $g$ and
$T(h)$ are equal in $L^{1}(P)$. Since every sequence converging in $L^{1}(P)$
has an almost surely converging subsequence, there is a subsequence
$(h_{n_{i}})_{i\in\mathbb{N}}$ such that

$\langle h_{n_{i}},\Phi(x)\rangle_{\mathcal{H}(k)}\rightarrow g(x)$

for $P$-almost all $x\in\mathcal{X}$. Since $\langle
h_{n},\Phi(x)\rangle_{\mathcal{H}(k)}\rightarrow\langle
h,\Phi(x)\rangle_{\mathcal{H}(k)}$ for all $x\in\mathcal{X}$, it follows that
$g(x)=\langle h,\Phi(x)\rangle_{\mathcal{H}(k)}=T(h)(x)$ for $P$-almost all
$x\in\mathcal{X}$. Thus $g$ and $T(h)$ are equal in $L^{1}(P)$ and the graph
of $T$ is indeed closed. The conditions of the closed graph theorem, which
states that a linear map $T$ between Banach spaces is bounded if and only if
its graph is closed, have now been verified. Thus $T$ is bounded and we denote

$\|T\|:=\sup_{\|h\|_{\mathcal{H}(k)}\leq 1}\int|T(h)|\;\mathrm{d}P<\infty.$

This allows us to conclude that $\mathrm{P}_{\Phi}$ is a bounded linear
functional, since

$\displaystyle|\mathrm{P}_{\Phi}(h)|=\left|\int\langle
h,\Phi(x)\rangle_{\mathcal{H}(k)}\;\mathrm{d}P(x)\right|\leq\int|\langle
h,\Phi(x)\rangle_{\mathcal{H}(k)}|\;\mathrm{d}P(x)=\int|T(h)|\;\mathrm{d}P\leq\|T\|\|h\|_{\mathcal{H}(k)},$

as required. ∎

######  Corollary 15 (Characterisation of $\mathcal{P}_{k}(\mathcal{X})$).

$P\in\mathcal{P}_{k}(\mathcal{X})$ if and only if $\mathcal{H}(k)\subset
L^{1}(P)$.

###### Proof.

Take $\Phi(x)=k(\cdot,x)$ to be the _canonical feature map_ , so that $k$
being measurable implies the scalar functions $x\mapsto\langle
h,\Phi(x)\rangle=h(x)$ are measurable and $\mathrm{P}=\mathrm{P}_{\Phi}$. If
$\mathcal{H}(k)\subset L^{1}(P)$ then $x\mapsto k(\cdot,x)$ is scalarly
$P$-integrable, and Proposition 14 shows that $\mathrm{P}$ is a bounded linear
functional. Conversely, if $\mathcal{H}(k)\not\subset L^{1}(P)$ then it is
clear that $\mathrm{P}$ is not bounded and thus
$P\notin\mathcal{P}_{k}(\mathcal{X})$. ∎

The question then reduces to when $\mathcal{H}(k)\subset L^{1}(P)$.

######  Proposition 16.

If $\int\sqrt{k(x,x)}\mathrm{d}P(x)<\infty$ then $\mathcal{H}(k)\subset
L^{1}(P)$.

###### Proof.

For $h\in\mathcal{H}(k)$, from the reproducing property and Cauchy–Schwarz,

$\displaystyle\int|h(x)|\;\mathrm{d}P(x)=\int|\langle
h,k(\cdot,x)\rangle|\;\mathrm{d}P(x)\leq\|h\|_{\mathcal{H}(k)}\int\|k(\cdot,x)\|_{\mathcal{H}(k)}\;\mathrm{d}P(x),$

where the reproducing property again yields
$\|k(\cdot,x)\|_{\mathcal{H}(k)}=\sqrt{k(x,x)}$, as required. ∎

One can weaken the above integrability condition under mild assumptions on $k$
and $\mathcal{X}$:

######  Proposition 17.

If $k$ is continuous, $\mathcal{X}$ is separable, and
$\iint|k(x,y)|\;\mathrm{d}P(x)\mathrm{d}P(y)<\infty$, then
$\mathcal{H}(k)\subset L^{1}(P)$.

###### Proof.

Since $k$ is measurable and real-valued, $k$ is _strongly measurable_ in the
sense of Carmeli et al. (2006, Section 3.1). Further, since $k$ is strongly
measurable and $\iint|k(x,y)|\;\mathrm{d}P(x)\mathrm{d}P(y)<\infty$, then $k$
is _$\infty$ -bounded_ in the sense of Carmeli et al. (2006, Definition 4.1);
see Carmeli et al. (2006, Corollary 4.3). Since $k$ is continuous and
$\mathcal{X}$ is separable, it follows that $\mathcal{H}(k)$ is separable
(Carmeli et al., 2006, Corollary 5.2). Since $\mathcal{H}(k)$ is separable,
$k$ being $\infty$-bounded is equivalent to $\mathcal{H}(k)\subset L^{1}(P)$
(Carmeli et al., 2006, Proposition 4.4). ∎

##  Appendix B Auxiliary Results

The following auxiliary results were used for the proofs in the main text:

######  Theorem 18 (Strong law of large numbers in a Banach space).

Let $\mathcal{B}$ be a separable Banach space and let
$(\zeta_{n})_{n\in\mathbb{N}}$ be a sequence of independent and identically
distributed random variables taking values in $\mathcal{B}$ such that
$\mathbb{E}[\|\zeta_{n}\|_{\mathcal{B}}]<\infty$ and
$\mathbb{E}[\zeta_{n}]=0$. Then

$\frac{1}{n}\sum_{i=1}^{n}\zeta_{i}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}0.$

###### Proof.

See Ledoux and Talagrand (1991, Corollary 7.10). ∎

######  Theorem 19 (Strong law of larger numbers for V-statistics).

Let $P$ be a distribution on a measurable space $\mathcal{X}$, and let
$(x_{n})_{n\in\mathbb{N}}$ be a sequence of independent draws from $P$. Let
$v:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ be a measurable
symmetric function with $\iint|v(x,y)|\mathrm{d}P(x)\mathrm{d}P(y)<\infty$ and
$\int|v(x,x)|\mathrm{d}P(x)<\infty$. Then

$\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}v(x_{i},x_{j})\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\iint
v(x,y)\mathrm{d}P(x)\mathrm{d}P(y).$

###### Proof.

Let

$\displaystyle
V_{n}:=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}v(x_{i},x_{j}),\qquad
U_{n}:=\frac{1}{n(n-1)}\sum_{i=1}^{n}\sum_{j\neq i}v(x_{i},x_{j}),\qquad
R_{n}:=\sum_{i=1}^{n}v(x_{i},x_{i}).$

From the usual strong law of large numbers, we have that
$n^{-1}R_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\int
v(x,x)\;\mathrm{d}P(x)<\infty$ and thus
$n^{-2}R_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}0$.
From the strong law of large numbers for U-statistics, due to Hoeffding
(1961), we have that
$U_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\iint
v(x,y)\mathrm{d}P(x)\mathrm{d}P(y)$. Thus

$V_{n}=\frac{n-1}{n}U_{n}+\frac{1}{n^{2}}R_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}1\times\iint
v(x,y)\mathrm{d}P(x)\mathrm{d}P(y)+0,$

which completes the argument. ∎

######  Theorem 20 (Asymptotic normality of V-statistics).

Let $P$ be a distribution on a measurable space $\mathcal{X}$, and let
$(x_{n})_{n\in\mathbb{N}}$ be a sequence of independent draws from $P$. Let
$v:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ be a measurable
symmetric function with $\iint v(x,y)^{2}\mathrm{d}P(x)\mathrm{d}P(y)<\infty$.
Let $\sigma^{2}=\mathbb{V}_{X\sim P}[\int v(X,y)\mathrm{d}P(y)]$. Then

$\sqrt{n}\left(\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}v(x_{i},x_{j})-\iint
v(x,y)\mathrm{d}P(x)\mathrm{d}P(y)\right)\operatorname*{\stackrel{{\scriptstyle\textup{d}}}{{\longrightarrow}}}\mathcal{N}(0,4\sigma^{2}),$

where we interpret $\mathcal{N}(0,0)$ as a point mass $\delta_{0}$ at
$0\in\mathcal{X}$.

###### Proof.

This proof uses the same notation as the proof of Theorem 19. From Hoeffding
(1948), since $\iint v(x,y)^{2}\mathrm{d}P(x)\mathrm{d}P(y)<\infty$, we have
asymptotic normality of the U-statistic

$\sqrt{n}\left(U_{n}-\iint
v(x,y)\mathrm{d}P(x)\mathrm{d}P(y)\right)\operatorname*{\stackrel{{\scriptstyle\textup{d}}}{{\longrightarrow}}}\mathcal{N}(0,4\sigma^{2})$

(see also Section 5.5.1 of Serfling (2009)). From Theorem 1 of Bonner and
Kirschner (1977), since $\iint v(x,y)\mathrm{d}P(x)\mathrm{d}P(y)<\infty$, we
have that
$\sqrt{n}(V_{n}-U_{n})\operatorname*{\stackrel{{\scriptstyle\textup{p}}}{{\longrightarrow}}}0$
(see also Section 5.7.3 of Serfling (2009)). From Slutsky’s theorem, we can
combine these facts to arrive at the required result. ∎

######  Theorem 21 (Uniform law of large numbers).

Let $\Theta\subset\mathbb{R}^{p}$ be bounded. Let
$f:\Theta\rightarrow\mathbb{R}$ be fixed and consider a sequence
$(f_{n})_{n\in\mathbb{N}}$ of stochastic functions
$f_{n}:\Theta\rightarrow\mathbb{R}$. Suppose that

  * •

$f_{n}(\theta)\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}f(\theta)$
for all $\theta\in\Theta$

  * •

$|f_{n}(\theta)-f_{n}(\vartheta)|\leq B_{n}\|\theta-\vartheta\|$ for all
$\theta,\vartheta\in\Theta$, where $B_{n}$ does not depend on
$\theta,\vartheta$ and a.s. $\limsup B_{n}<\infty$.

Then
$\sup_{\theta\in\Theta}|f_{n}(\theta)-f(\theta)|\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}0$.

###### Proof.

This follows as a combination of Theorem 21.8 and Theorem 21.10 from Davidson
(1994). ∎

######  Definition 22.

Let $P$ be a distribution on a measurable space $\mathcal{X}$. A function
$f:\Theta\times\mathcal{X}\rightarrow\mathbb{R}^{p}$ is _locally uniformly
integrably bounded_ with respect to $P$ if, for every $\theta\in\Theta$, there
is a non-negative function $b_{\theta}:\mathcal{X}\rightarrow\mathbb{R}$ that
is integrable with respect to $P$, and an open neighbourhood $U_{\theta}$ of
$\theta$, such that $\|f(\vartheta,x)\|\leq b_{\theta}(x)$ for all
$\vartheta\in U_{\theta}$.

######  Lemma 23 (Differentiate under the Lebesgue integral).

Let $P$ be a distribution on a measurable space $\mathcal{X}$. Let
$\Theta\subseteq\mathbb{R}^{p}$ be an open set. Let
$h:\Theta\times\mathcal{X}\rightarrow\mathbb{R}$ be such that
$\partial_{\theta}h$ is locally uniformly integrably bounded with respect to
$P$. Then

$\partial_{\theta}\int
h(\theta,x)\mathrm{d}P(x)=\int\partial_{\theta}h(\theta,x)\mathrm{d}P(x).$

###### Proof.

See, for example, Aliprantis and Burkinshaw (1998, Theorem 24.5, p.193),
Billingsley (1979, Theorem 16.8, pp.181-182). ∎

##  Appendix C Assumptions for the Exponential Family Model

This appendix contains full details on how our theoretical assumptions were
verified for for the canonical exponential family statistical model in Section
4.4.

###  C.1 Estimation via Generalised Method of Moments

In this section we verify that collectively B1-B5 imply A1-A2 and A1-A7, so
that the generalised moment matching estimator is both consistent and
asymptotically normal. The matrix norm $\|M\|^{2}=\sum M_{i,j}^{2}$ will be
used. From B1, we have $|k(x,y)|\leq k_{\max}$ for all $x,y\in\mathbb{R}^{p}$
where $k_{\max}$ is a positive constant.

Verify A1: The space $\mathcal{H}(k)$ is separable whenever $k$ is continuous
and $\mathcal{X}$ is separable (which is the case for
$\mathcal{X}=\mathbb{R}^{d}$), and thus B1 implies A1.

Verify A2: B1 immediately implies A2.

Verify A1 and A4: At each $\theta\in S$, we claim the functions $x\mapsto
k(x,y)\partial_{\theta}p_{\theta}(x)$ and $x\mapsto\int
k(x,y)\partial_{\theta}p_{\theta}(x)\mathrm{d}P_{\vartheta}(y)$ are locally
uniformly integrably bounded with respect to the Lebesgue measure on
$\mathbb{R}^{p}$. Indeed, for the first of these functions we have a bound

$\displaystyle\int\sup_{\vartheta\in
S_{\gamma,\theta}}\|k(x,y)\partial_{\vartheta}p_{\vartheta}(x)\|\;\mathrm{d}x$
$\displaystyle\leq k_{\max}\int\sup_{\vartheta\in
S_{\gamma,\theta}}\|t(x)-\partial_{\vartheta}a(\vartheta)\|p_{\vartheta}(x)\;\mathrm{d}x$
$\displaystyle=k_{\max}\int\sup_{\vartheta\in
S_{\gamma,\theta}}\|t(x)-\partial_{\vartheta}a(\vartheta)\|\frac{p_{\vartheta}(x)}{p_{\theta}(x)}\;\mathrm{d}P_{\theta}(x),$

where
$S_{\gamma,\theta}=S\cap\\{\vartheta\in\mathbb{R}^{p}:\|\vartheta-\theta\|<\gamma_{\theta}\\}$.
An identical bound holds also for the second function considered. Now,

$\displaystyle\sup_{\vartheta\in
S_{\gamma,\theta}}\|t(x)-\partial_{\vartheta}a(\vartheta)\|$
$\displaystyle\leq\|t(x)\|+C_{1}$
$\displaystyle\frac{p_{\vartheta}(x)}{p_{\theta}(x)}$
$\displaystyle=\exp\left(\langle\vartheta-\theta,t(x)\rangle-[a(\vartheta)-a(\theta)]\right)\leq\exp\left(C_{2}\|t(x)\|+C_{3}\right)$

where $C_{1}=\sup_{\vartheta\in
S_{\gamma,\theta}}\|\partial_{\vartheta}a(\vartheta)\|$,
$C_{2}=\sup_{\vartheta\in S_{\gamma,\theta}}\|\vartheta-\theta\|$ and
$C_{3}=\sup_{\vartheta\in S_{\gamma,\theta}}|a(\vartheta)-a(\theta)|$. From
B5, we have $C_{1},C_{3}<\infty$, and $C_{2}<\infty$ holds since $S$ is
bounded. Thus

$\displaystyle\int\sup_{\vartheta\in
S_{\gamma,\theta}}\|k(x,y)\partial_{\vartheta}p_{\vartheta}(x)\|\;\mathrm{d}x$
$\displaystyle\leq
k_{\max}\int\left[\|t(x)\|+C_{1}\right]\exp\left(C_{2}\|t(x)\|+C_{3}\right)\;\mathrm{d}P_{\theta}(x)<\infty,$

where the final integral is finite due to B4. Thus we have established the
conditions of Lemma 23 and we may interchange $\partial_{\theta}$ with
integration with respect to $\mathrm{d}x$ and $\mathrm{d}y$:

$\displaystyle\partial_{\theta}k_{\theta}(x,y)$ $\displaystyle=-\int
k(x,y)\partial_{\theta}p_{\theta}(x)\mathrm{d}x-\int
k(x,y)\partial_{\theta}p_{\theta}(y)\mathrm{d}y+\iint
k(x,y)\partial_{\theta}[p_{\theta}(x)p_{\theta}(y)]\mathrm{d}x\mathrm{d}y$
$\displaystyle=-\int
k(x,y)[t(x)-\partial_{\theta}a(\theta)]\mathrm{d}P_{\theta}(x)-\int
k(x,y)[t(y)-\partial_{\theta}a(\theta)]\mathrm{d}P_{\theta}(y)$
$\displaystyle\qquad+\iint
k(x,y)[t(x)+t(y)-2\partial_{\theta}a(\theta)]\mathrm{d}P_{\theta}(x)\mathrm{d}P_{\theta}(y)$

so, using the triangle inequality and Jensen’s inequality,

$\displaystyle\frac{\|\partial_{\theta}k_{\theta}(x,y)\|}{k_{\max}}$
$\displaystyle\leq
4\int\|t(x)\|\mathrm{d}P_{\theta}(x)+4\|\partial_{\theta}a(\theta)\|.$

From B3 we have that $\sup_{\theta\in S}\int\|t\|\mathrm{d}P_{\theta}<\infty$
and, since $S$ is bounded, B5 implies that $\sup_{\theta\in
S}\|\partial_{\theta}a(\theta)\|<\infty$. Thus
$\sup_{x,y\in\mathcal{X}}\sup_{\theta\in
S}\|\partial_{\theta}k_{\theta}(x,y)\|<\infty$ and, in particular, A1 and A4
hold.

Verify A2: An argument analogous to the previous argument shows that B4 is
sufficient to allow a second interchange of $\partial_{\theta}$ with
integration with respect to $\mathrm{d}x$ and $\mathrm{d}y$:

$\displaystyle\partial_{\theta}^{2}k_{\theta}(x,y)$
$\displaystyle=[\partial_{\theta}^{2}a(\theta)]\iint
k(x,y)\mathrm{d}P_{\theta}(x)-\int
k(x,y)[t(x)-\partial_{\theta}a(\theta)][t(x)-\partial_{\theta}a(\theta)]^{\top}\mathrm{d}P_{\theta}(x)$
$\displaystyle\qquad+[\partial_{\theta}^{2}a(\theta)]\iint
k(x,y)\mathrm{d}P_{\theta}(y)-\int
k(x,y)[t(y)-\partial_{\theta}a(\theta)][t(y)-\partial_{\theta}a(\theta)]^{\top}\mathrm{d}P_{\theta}(y)$
$\displaystyle\qquad-2[\partial_{\theta}^{2}a(\theta)]\iint
k(x,y)\mathrm{d}P_{\theta}(x)\mathrm{d}P_{\theta}(y)$ (13)
$\displaystyle\qquad+\iint
k(x,y)[t(x)+t(y)-2\partial_{\theta}a(\theta)][t(x)+t(y)-2\partial_{\theta}a(\theta)]^{\top}\mathrm{d}P_{\theta}(x)\mathrm{d}P_{\theta}(y)$

where the quantities exist on $\theta\in S$ due to B3 and B5, showing that A2
is satisfied.

Verify A3: Since (13) takes the form

$\displaystyle\partial_{\theta}^{2}k_{\theta}(x,y)$ $\displaystyle=\int
k(x,u)f(u,\theta)\;\mathrm{d}u+\int k(y,u)f(u,\theta)\;\mathrm{d}u-\iint
k(u,v)g(u,v,\theta)\;\mathrm{d}u\mathrm{d}v$ $\displaystyle f(u,\theta)$
$\displaystyle:=\left\\{[\partial_{\theta}^{2}a(\theta)]-[t(u)-\partial_{\theta}a(\theta)][t(u)-\partial_{\theta}a(\theta)]^{\top}\right\\}p_{\theta}(u)$
$\displaystyle g(u,v,\theta)$
$\displaystyle=\left\\{2[\partial_{\theta}^{2}a(\theta)]-[t(u)+t(v)-2\partial_{\theta}a(\theta)][t(u)+t(v)-2\partial_{\theta}a(\theta)]^{\top}\right\\}p_{\theta}(u)p_{\theta}(v),$

we have that

$\displaystyle\sup_{x,y\in\mathcal{X}}\frac{\|\partial_{\theta}^{2}k_{\theta}(x,y)-[\partial_{\theta}^{2}k_{\theta}(x,y)|_{\theta=\theta_{\star}}]\|}{k_{\max}}$
$\displaystyle\leq 2\int|f(u,\theta)-f(u,\theta_{\star})|\;\mathrm{d}u$
$\displaystyle\qquad+\iint|g(u,v,\theta)-g(u,v,\theta_{\star})|\;\mathrm{d}u\mathrm{d}v$

Since $\partial_{\theta}^{2}a(\theta)$ is continuous at $\theta_{\star}$ from
B5 and $\int t\;\mathrm{d}P_{\theta}$ and $\int
tt^{\top}\;\mathrm{d}P_{\theta}$ are continuous at $\theta_{\star}$ from B2,
it follows that both integrals on the right hand side vanish as
$\theta\rightarrow\theta_{\star}$, meaning A3 is established.

Verify A5 and A6: Application of the triangle inequality to (13), and
collecting together terms, yields

$\displaystyle\frac{\|\partial_{\theta}^{2}k_{\theta}(x,y)\|}{k_{\max}}$
$\displaystyle\leq
4\int\|t(x)\|^{2}\mathrm{d}P_{\theta}(x)+12\|\partial_{\theta}a(\theta)\|\int\|t(x)\|\mathrm{d}P_{\theta}(x)$
$\displaystyle\qquad+4\|\partial_{\theta}^{2}a(\theta)\|+6\|\partial_{\theta}a(\theta)\|^{2}+2\left(\int\|t(x)\|\mathrm{d}P_{\theta}(x)\right)^{2}.$

Thus B3 and B5 imply
$\sup_{x,y\in\mathcal{X}}\left[\left.\|\partial_{\theta}^{2}k_{\theta}(x,y)\|\right|_{\theta=\theta_{\star}}\right]<\infty$,
so A5 and A6 hold.

Verify A7: Integrating both arguments of (13) with respect to
$P=P_{\theta_{\star}}$, we obtain

$\displaystyle\Gamma$
$\displaystyle=\frac{1}{2}\left.\iint\partial_{\theta}^{2}k_{\theta}(x,y)\mathrm{d}P(x)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}$
$\displaystyle=\iint
k(x,y)[t(x)-\partial_{\theta}a(\theta)][t(y)-\partial_{\theta}a(\theta)]^{\top}\;\mathrm{d}P(x)\mathrm{d}P(y)=MM^{\top}\succ
0,$

where the final equality uses the definition of $M$ and the fact that
$k(x,y)=\langle\phi(x),\phi(y)\rangle$, and the positive definiteness follows
from B2.

Verify
$\theta_{n}{\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}}\theta_{\star}$:
Through the above, we have proven that if
$\theta_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\theta_{\star}$,
then the fluctuations of $\theta_{n}$ are asymptotically normal. However, we
have not yet commented on whether
$\theta_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\theta_{\star}$.
Since we are in the well-specified regime, we have from (13)

$\displaystyle\left.\partial_{\theta}^{2}D_{k}(P,P_{\theta})\right|_{\theta=\theta^{*}}=\left.\iint\partial_{\theta}^{2}k_{\theta}(x,y)\mathrm{d}P(x)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}=2\Gamma\succ
0,$

showing that $\theta_{\star}$ is a local minimum of the kernel discrepancy.
Since we assumed $\theta_{\star}$ is unique, and from B1 the map
$\theta\mapsto D_{k}(P,P_{\theta})$ is continuous, the conditions of Theorem 8
are satisfied, establishing that
$\theta_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\theta_{\star}$.

###  C.2 Estimation via Kernel Stein Discrepancy

In this section we verify that collectively C1-C3 imply A1 and A1-A7, so that
the minimum kernel Stein discrepancy estimator is both consistent and
asymptotically normal. From C1, we have
$|c(x,y)|,\|\nabla_{x}c(x,y)\|,\|\nabla_{y}c(x,y)\|\leq c_{\max}$ for all
$x,y\in\mathbb{R}^{d}$ where $c_{\max}$ is a positive constant. From C2 the
integrability conditions of Proposition 4 are satisfied, so that
$k(\cdot,\cdot;\theta)=k_{\theta}(\cdot,\cdot)$ for each $\theta\in\Theta$.

Verify A1: From C1 and C2, $(x,y)\mapsto k(x,y;\theta)=k_{\theta}(x,y)$ is
continuous, since

$\displaystyle k_{\theta}(x,y)$
$\displaystyle=\nabla_{x}\cdot\nabla_{y}c(x,y)+\langle\nabla_{x}c(x,y),\nabla
t(y)\theta+\nabla b(y)\rangle$
$\displaystyle\qquad+\langle\nabla_{y}c(x,y),\nabla t(x)\theta+\nabla
b(x)\rangle+c(x,y)\langle\nabla t(x)\theta+\nabla b(x),\nabla
t(y)\theta+\nabla b(y)\rangle.$

The space $\mathcal{H}(k)$ is separable whenever $k$ is continuous and
$\mathcal{X}$ is separable (which is the case for
$\mathcal{X}=\mathbb{R}^{d}$), and thus $\mathcal{H}(k)$ is separable and we
have A1.

Verify A2: First we compute

$\displaystyle\partial_{\theta_{i}}k_{\theta}(x,y)$
$\displaystyle=\langle\nabla_{x}c(x,y),[\nabla
t(y)]_{\cdot,i}\rangle+\langle\nabla_{y}c(x,y),[\nabla t(x)]_{\cdot,i}\rangle$
(14) $\displaystyle\qquad+c(x,y)\langle[\nabla t(x)]_{\cdot,i},\nabla
t(y)\theta+\nabla b(y)\rangle+c(x,y)\langle\nabla t(x)\theta+\nabla
b(x),[\nabla t(y)]_{\cdot,i}\rangle,$

which is a linear function of $\theta$, so that we trivially have A2.

Verify A2, A3, A1, A4: The triangle and Cauchy–Schwarz inequalities applied to
(14) yield

$\displaystyle\frac{\|\partial_{\theta_{i}}k_{\theta}(x,y)\|}{c_{\max}}$
$\displaystyle\leq\|\nabla t(x)\|+\|\nabla t(y)\|+\|\nabla
t(x)\|\left(\|\nabla t(y)\|\sup_{\theta\in\Theta}\|\theta\|+\|\nabla
b(y)\|\right)$ $\displaystyle\qquad+\|\nabla t(y)\|\left(\|\nabla
t(x)\|\sup_{\theta\in\Theta}\|\theta\|+\|\nabla b(x)\|\right),$

where $\sup_{\theta\in\Theta}\|\theta\|<\infty$ since we assumed that $\Theta$
was bounded. Thus A2 holds from C3, and both A3 and A1 hold from C3. Squaring,
we see also that A4 holds from C3.

Verify A3, A5, A6: Differentiating (14) again,

$\displaystyle\partial_{\theta_{i}}\partial_{\theta_{j}}k_{\theta}(x,y)$
$\displaystyle=2c(x,y)\langle[\nabla t(x)]_{\cdot,i},[\nabla
t(y)]_{\cdot,j}\rangle$

so A3 holds trivially, since these functions are constant in $\theta$.
Further,

$\displaystyle\frac{\|\partial_{\theta_{i}}\partial_{\theta_{j}}k_{\theta}(x,y)\|}{c_{\max}}$
$\displaystyle\leq 2\|\nabla t(x)\|\|\nabla t(y)\|$

so A5 and A6 hold from C3.

Verify A7: Finally,

$\displaystyle\Gamma=\frac{1}{2}\left.\iint\partial_{\theta}^{2}k_{\theta}(x,y)\;\mathrm{d}P(x)\mathrm{d}P(y)\right|_{\theta=\theta_{\star}}=\iint
c(x,y)[\nabla t(x)]^{\top}[\nabla t(y)]\;\mathrm{d}P(x)\mathrm{d}P(y)\succ 0$

from C3, establishing A7.

Verify
$\theta_{n}{\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}}\theta_{\star}$:
Since we have verified A1, A2 and A3, we have from Lemma 11 and an essentially
identical argument to Lemma 7 that
$D_{k}(P,P_{\theta_{n}})\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}D_{k}(P,P_{\theta_{\star}})$.
Since the squared kernel Stein discrepancy

$\theta\mapsto D_{k}(P,P_{\theta})^{2}=\iint
k_{\theta}(x,y)\;\mathrm{d}P(x)\mathrm{d}P(y)$

is a quadratic with minimum $\theta_{\star}$ and Hessian at $\theta_{\star}$
equal to $2\Gamma\succ 0$, it is necessarily the case that
$\theta_{n}\operatorname*{\stackrel{{\scriptstyle\textup{a.s.}}}{{\longrightarrow}}}\theta_{\star}$
in the $n\rightarrow\infty$ limit.

## References

  * Akaike (1973) H. Akaike.  Information theory and an extension of the likelihood principle.  In _Proceedings of the Second International Symposium of Information Theory_ , 1973.
  * Aliprantis and Burkinshaw (1998) C. D. Aliprantis and O. Burkinshaw.  _Principles of Real Analysis_.  Academic Press, 1998.
  * Alquier and Gerber (2020) P. Alquier and M. Gerber.  Universal robust regression via maximum mean discrepancy.  _arXiv:2006.00840_ , 2020.
  * Anastasiou et al. (2022) A. Anastasiou, A. Barp, F.-X. Briol, B. Ebner, R. E. Gaunt, F. Ghaderinezhad, J. Gorham, A. Gretton, C. Ley, Q. Liu, L. Mackey, C. J. Oates, G. Reinert, and Y. Swan.  Stein’s method meets statistics: A review of some recent developments.  _Statistical Science_ , 2022.  To appear.
  * Arjovsky et al. (2017) M. Arjovsky, S. Chintala, and L. Bottou.  Wasserstein generative adversarial networks.  In _Proceedings of the 34th International Conference on Machine Learning_ , 2017.
  * Barp et al. (2019) A. Barp, F.-X. Briol, A. Duncan, M. Girolami, and L. Mackey.  Minimum Stein discrepancy estimators.  In _Proceedings of the 33rd Conference on Neural Information Processing Systems_ , 2019.
  * Barp et al. (2022) A. Barp, C.-J. Simon-Gabriel, M. Girolami, and L. Mackey.  Targeted separation and convergence with kernel discrepancies.  _arXiv:2209.12835_ , 2022.
  * Basu et al. (2011) A. Basu, H. Shioya, and C. Park.  _Statistical Inference: The Minimum Distance Approach_.  CRC Press, 2011.
  * Beaumont (2019) M. A. Beaumont.  Approximate bayesian computation.  _Annual Review of Statistics and its Application_ , 6:379–403, 2019.
  * Billingsley (1979) P. Billingsley.  _Probability and Measure_.  John Wiley and Sons, 1979.
  * Bińkowski et al. (2018) M. Bińkowski, D. J. Sutherland, M. Arbel, and A. Gretton.  Demystifying MMD GANs.  In _Proceedings of the 6th International Conference on Learning Representations_ , 2018.
  * Bonner and Kirschner (1977) N. Bonner and H.-P. Kirschner.  Note on conditions for weak convergence of von Mises’ differentiable statistical functions.  _The Annals of Statistics_ , 5(2):405–407, 1977\.
  * Briol et al. (2019) F.-X. Briol, A. Barp, A. B. Duncan, and M. Girolami.  Statistical inference for generative models with maximum mean discrepancy.  _arXiv:1906.05944_ , 2019.
  * Carmeli et al. (2006) C. Carmeli, E. De Vito, and A. Toigo.  Vector valued reproducing kernel Hilbert spaces of integrable functions and Mercer theorem.  _Analysis and Applications_ , 4(04):377–408, 2006\.
  * Chérief-Abdellatif and Alquier (2020) B.-E. Chérief-Abdellatif and P. Alquier.  MMD-Bayes: Robust Bayesian estimation via maximum mean discrepancy.  In _Symposium on Advances in Approximate Bayesian Inference_ , pages 1–21. PMLR, 2020.
  * Chérief-Abdellatif and Alquier (2022) B.-E. Chérief-Abdellatif and P. Alquier.  Finite sample properties of parametric MMD estimation: Robustness to misspecification and dependence.  _Bernoulli_ , 28(1):181–213, 2022.
  * Chwialkowski et al. (2016) K. Chwialkowski, H. Strathmann, and A. Gretton.  A kernel test of goodness of fit.  In _Proceedings of the 33rd International Conference on Machine Learning_ , 2016.
  * Cortes and Scott (2016) E. C. Cortes and C. Scott.  Sparse approximation of a kernel mean.  _IEEE Transactions on Signal Processing_ , 65(5):1310–1323, 2016.
  * Davidson (1994) J. Davidson.  _Stochastic Limit Theory: An Introduction for Econometricians_.  OUP Oxford, 1994.
  * Dawid (2007) A. P. Dawid.  The geometry of proper scoring rules.  _Annals of the Institute of Statistical Mathematics_ , 59(1):77–93, 2007.
  * Dellaporta et al. (2022) C. Dellaporta, J. Knoblauch, T. Damoulas, and F.-X. Briol.  Robust Bayesian inference for simulator-based models via the MMD posterior bootstrap.  In _Proceedings of the 25th International Conference on Artificial Intelligence and Statistics_ , 2022.
  * Dhariwal and Nichol (2021) P. Dhariwal and A. Nichol.  Diffusion models beat GANs on image synthesis.  In _Proceedings of the 35th Conference on Neural Information Processing Systems_ , 2021.
  * Dick and Pillichshammer (2010) J. Dick and F. Pillichshammer.  _Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo integration_.  Cambridge University Press, 2010.
  * Donoho and Liu (1988) D. L. Donoho and R. C. Liu.  The “automatic” robustness of minimum distance functionals.  _The Annals of Statistics_ , 16(2):552–586, 1988\.
  * Dunford (1937) N. Dunford.  Integration of vector-valued functions.  _Bulletin of the American Mathematical Society_ , page 43, 1937.
  * Dziugaite et al. (2015) G. K. Dziugaite, D. M. Roy, and Z. Ghahramani.  Training generative neural networks via maximum mean discrepancy optimization.  In _Proceedings of the 31st Conference on Uncertainty in Artificial Intelligence_ , 2015.
  * Frazier and Drovandi (2021) D. T. Frazier and C. Drovandi.  Robust approximate Bayesian inference with synthetic likelihood.  _Journal of Computational and Graphical Statistics_ , 30(4):958–976, 2021.
  * Freedman (2006) D. A. Freedman.  On the so-called “Huber sandwich estimator” and “robust standard errors”.  _The American Statistician_ , 60(4):299–302, 2006\.
  * Genevay et al. (2018) A. Genevay, G. Peyré, and M. Cuturi.  Learning generative models with sinkhorn divergences.  In _Proceedings of the 21st International Conference on Artificial Intelligence and Statistics_ , 2018.
  * Gneiting and Raftery (2007) T. Gneiting and A. E. Raftery.  Strictly proper scoring rules, prediction, and estimation.  _Journal of the American statistical Association_ , 102(477):359–378, 2007.
  * Goodfellow et al. (2020) I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.  Generative adversarial networks.  _Communications of the ACM_ , 63(11):139–144, 2020.
  * Gorham and Mackey (2017) J. Gorham and L. Mackey.  Measuring sample quality with kernels.  In _Proceedings of the 34th International Conference on Machine Learning_ , 2017.
  * Gretton et al. (2012) A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola.  A kernel two-sample test.  _The Journal of Machine Learning Research_ , 13(1):723–773, 2012.
  * Hansen (1982) L. P. Hansen.  Large sample properties of generalized method of moments estimators.  _Econometrica_ , pages 1029–1054, 1982.
  * Hickernell (1998) F. Hickernell.  A generalized discrepancy and quadrature error bound.  _Mathematics of Computation_ , 67(221):299–322, 1998.
  * Hoeffding (1948) W. Hoeffding.  A class of statistics with asymptotically normal distribution.  _The Annals of Mathematical Statistics_ , 19(3):293–325, 1948.
  * Hoeffding (1961) W. Hoeffding.  The strong law of large numbers for ${U}$-statistics.  Technical report, North Carolina State University. Dept. of Statistics, 1961.
  * Huber (1964) P. J. Huber.  Robust estimation of a location parameter.  _The Annals of Mathematical Statistics_ , pages 73–101, 1964.
  * Hyvärinen and Dayan (2005) A. Hyvärinen and P. Dayan.  Estimation of non-normalized statistical models by score matching.  _Journal of Machine Learning Research_ , 6(4), 2005.
  * Key et al. (2021) O. Key, T. Fernandez, A. Gretton, and F.-X. Briol.  Composite goodness-of-fit tests with kernels.  _arXiv:2111.10275_ , 2021.
  * Lam (2020) C. Lam.  High-dimensional covariance matrix estimation.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 12(2):e1485, 2020.
  * LeCun et al. (2007) Y. LeCun, S. Chopra, R. Hadsell, M. Ranzato, and F. Huang.  A tutorial on energy-based learning.  In B. Schölkopf, A. J. Smola, B. Taskar, and S. Vishwanathan, editors, _Predicting Structured Data_. 2007.
  * Ledoux and Talagrand (1991) M. Ledoux and M. Talagrand.  _Probability in Banach Spaces: Isoperimetry and Processes_.  Springer Science & Business Media, 1991.
  * Li et al. (2017) C.-L. Li, W.-C. Chang, Y. Cheng, Y. Yang, and B. Póczos.  MMD GAN: Towards deeper understanding of moment matching network.  In _Proceedings of the 31st Conference on Neural Information Processing Systems_ , 2017.
  * Li et al. (2015) Y. Li, K. Swersky, and R. Zemel.  Generative moment matching networks.  In _Proceedings of the 32nd International Conference on Machine Learning_ , 2015.
  * Liu et al. (2016) Q. Liu, J. Lee, and M. Jordan.  A kernelized Stein discrepancy for goodness-of-fit tests.  In _Proceedings of the 33rd International Conference on Machine Learning_ , 2016.
  * Lyne et al. (2015) A.-M. Lyne, M. Girolami, Y. Atchadé, H. Strathmann, and D. Simpson.  On Russian roulette estimates for Bayesian inference with doubly-intractable likelihoods.  _Statistical Science_ , 30(4):443–467, 2015.
  * Matsubara et al. (2022a) T. Matsubara, J. Knoblauch, F.-X. Briol, and C. J. Oates.  Robust generalised bayesian inference for intractable likelihoods.  _Journal of the Royal Statistical Society, Series B_ , 84(3):997–1022, 2022a.
  * Matsubara et al. (2022b) T. Matsubara, J. Knoblauch, F.-X. Briol, and C. J. Oates.  Robust generalised Bayesian inference for intractable likelihoods.  _Journal of the Royal Statistical Society: Series B_ , 84(3):997–1022, 2022b.
  * Mitrovic et al. (2016) J. Mitrovic, D. Sejdinovic, and Y.-W. Teh.  DR-ABC: Approximate Bayesian computation with kernel-based distribution regression.  In _Proceedings of the 33rd International Conference on Machine Learning_ , 2016.
  * Mroueh and Sercu (2017) Y. Mroueh and T. Sercu.  Fisher GAN.  In _Proceedings of the 31st Conference on Neural Information Processing Systems_ , 2017.
  * Mroueh et al. (2017) Y. Mroueh, T. Sercu, and V. Goel.  McGAN: Mean and covariance feature matching GAN.  In _Proceedings of the 34th International Conference on Machine Learning_ , 2017.
  * Mroueh et al. (2018) Y. Mroueh, C.-L. Li, T. Sercu, A. Raj, and Y. Cheng.  Sobolev GAN.  In _Proceedings of the 6th International Conference on Learning Representations_ , 2018.
  * Muandet et al. (2017) K. Muandet, K. Fukumizu, B. Sriperumbudur, and B. Schölkopf.  Kernel mean embedding of distributions: A review and beyond.  _Foundations and Trends® in Machine Learning_ , 10(1-2):1–141, 2017.
  * Müller (1997) A. Müller.  Integral probability metrics and their generating classes of functions.  _Advances in Applied Probability_ , 29(2):429–443, 1997.
  * Nietert et al. (2021) S. Nietert, Z. Goldfeld, and K. Kato.  Smooth $p$-wasserstein distance: Structure, empirical approximation, and statistical applications.  In _Proceedings of the 38th International Conference on Machine Learning_ , 2021.
  * Niu et al. (2021) Z. Niu, J. Meier, and F.-X. Briol.  Discrepancy-based inference for intractable generative models using quasi-Monte Carlo.  _arXiv:2106.11561_ , 2021.
  * Oates et al. (2017) C. J. Oates, M. Girolami, and N. Chopin.  Control functionals for Monte Carlo integration.  _Journal of the Royal Statistical Society, Series B_ , 79:695–718, 2017.
  * Pardo (2018) L. Pardo.  _Statistical Inference Based on Divergence Measures_.  Chapman and Hall/CRC, 2018.
  * Park et al. (2016) M. Park, W. Jitkrittum, and D. Sejdinovic.  K2-ABC: Approximate Bayesian computation with kernel embeddings.  In _Proceedings of the 18th International Conference on Artificial Intelligence and Statistics_ , 2016.
  * Rahimi and Recht (2007) A. Rahimi and B. Recht.  Random features for large-scale kernel machines.  In _Proceedings of the 21st Conference on Neural Information Processing Systems_ , 2007.
  * Schwabik and Ye (2005) S. Schwabik and G. Ye.  _Topics in Banach Space Integration_.  World Scientific, 2005.
  * Serfling (2009) R. J. Serfling.  _Approximation Theorems of Mathematical Statistics_.  John Wiley & Sons, 2009.
  * Simon-Gabriel et al. (2020) C.-J. Simon-Gabriel, A. Barp, and L. Mackey.  Metrizing weak convergence with maximum mean discrepancies.  _arXiv:2006.09268_ , 2020.
  * Song et al. (2008) L. Song, X. Zhang, A. Smola, A. Gretton, and B. Schölkopf.  Tailoring density estimation via reproducing kernel moment matching.  In _Proceedings of the 25th International Conference on Machine Learning_ , 2008.
  * Song and Kingma (2021) Y. Song and D. P. Kingma.  How to train your energy-based models.  _arXiv:2101.03288_ , 2021.
  * Steinwart and Christmann (2008) I. Steinwart and A. Christmann.  _Support Vector Machines_.  Springer Science & Business Media, 2008.
  * Sutherland et al. (2017) D. J. Sutherland, H.-Y. Tung, H. Strathmann, S. De, A. Ramdas, A. J. Smola, and A. Gretton.  Generative models and model criticism via optimized maximum mean discrepancy.  In _Proceedings of the 5th International Conference on Learning Representations_ , 2017.
  * Teymur et al. (2021) O. Teymur, J. Gorham, M. Riabiz, and C. J. Oates.  Optimal quantisation of probability measures using maximum mean discrepancy.  In _Proceedings of the 24th International Conference on Artificial Intelligence and Statistics_ , 2021.
  * Theis et al. (2016) L. Theis, A. van den Oord, and M. Bethge.  A note on the evaluation of generative models.  In _Proceedings of the 4th International Conference on Learning Representations_ , 2016.
  * Wynne and Duncan (2022) G. Wynne and A. B. Duncan.  A kernel two-sample test for functional data.  _Journal of Machine Learning Research_ , 23(73):1–51, 2022.
  * Wynne et al. (2022) G. Wynne, M. Kasprzak, and A. B. Duncan.  A spectral representation of kernel Stein discrepancy with application to goodness-of-fit tests for measures on infinite dimensional Hilbert spaces.  _arXiv:2206.04552_ , 2022.