#  A Kernelised Stein Statistic for Assessing
Implicit Generative Models

Wenkai Xu
Department of Statistics
University of Oxford
<EMAIL_ADDRESS>
&Gesine Reinert
Department of Statistics
University of Oxford
<EMAIL_ADDRESS>

###### Abstract

Synthetic data generation has become a key ingredient for training machine
learning procedures, addressing tasks such as data augmentation, analysing
privacy-sensitive data, or visualising representative samples. Assessing the
quality of such synthetic data generators hence has to be addressed. As (deep)
generative models for synthetic data often do not admit explicit probability
distributions, classical statistical procedures for assessing model goodness-
of-fit may not be applicable. In this paper, we propose a principled procedure
to assess the quality of a synthetic data generator. The procedure is a
kernelised Stein discrepancy (KSD)-type test which is based on a non-
parametric Stein operator for the synthetic data generator of interest. This
operator is estimated from samples which are obtained from the synthetic data
generator and hence can be applied even when the model is only implicit. In
contrast to classical testing, the sample size from the synthetic data
generator can be as large as desired, while the size of the observed data
which the generator aims to emulate is fixed. Experimental results on
synthetic distributions and trained generative models on synthetic and real
datasets illustrate that the method shows improved power performance compared
to existing approaches.

##  1 Introduction

Synthetic data capturing main features of the original dataset are of
particular interest for machine learning methods. The use of original dataset
for machine learning tasks can be problematic or even prohibitive in certain
scenarios, e.g. under authority regularisation on privacy-sensitive
information, training models on small-sample dataset, or calibrating models
with imbalanced groups. High quality synthetic data generation procedures
surpass some of these challenges by creating de-identified data to preserve
privacy and to augment small or imbalance datasets. Training deep generative
models has been widely studied in the recent years (Kingma and Welling, 2013;
Radford et al., 2015; Song and Kingma, 2021) and methods such as those based
on Generative Adversarial Networks (GANs) (Goodfellow et al., 2014) provide
powerful approaches that learn to generate synthetic data which resemble the
original data distributions. However, these deep generative models usually do
not provide theoretical guarantees on the _goodness-of-fit_ to the original
data (Creswell et al., 2018).

To the best of our knowledge, existing mainstream developments for deep
generative models (Song and Ermon, 2020; Li et al., 2017) do not provide a
systematic approach to assess the quality of the synthetic samples. Instead,
heuristic methods are applied, e.g. for image data, the quality of samples are
generally decided via visual comparisons. The training quality has been
studied relying largely on the specific choice of training loss, which does
not directly translate into a measure of sample quality; in the case of the
log-likelihood (Theis et al., 2015). Common quality assessment measures for
implicit generative models, on images for example, include Inception Scores
(IS) (Salimans et al., 2016) and Fréchet Inception Distance (FID) (Heusel et
al., 2017), which are motivated by human inception systems in the visual
cortex and pooling (Wang et al., 2004). Bińkowski et al. (2018) pointed out
issues for IS and FID and developed the Kernel Inception Distance (KID) for
more general datasets. Although these scores can be used for for comparisons,
they do not provide a statistical significance test which would assess whether
a deemed good generative model is “good enough”. A key stumbling block is that
the distribution from which a synthetic method generates samples is not
available; one only ever observes samples from it.

For models in which the density is known explicitly, at least up to a
normalising constant, some assessment methods are available. Gorham and Mackey
(2017) proposed to assess sample quality using discrepancy measures called
kernelised Stein discrepancy (KSD). Schrab et al. (2022) assesses the quality
of generative models on the MNIST image dataset from LeCun et al. (1995) using
an aggregated kernel Stein discrepancy (KSDAgg) test; still an explicit
density is required. The only available implicit goodness-of-fit test, AgraSSt
(Xu and Reinert, 2022), applies only to generators of finite graphs; it is
also of KSD form and makes extensive use of the discrete and finite nature of
the problem. To date, quality assessment procedures of implicit deep
generative models for continuous data remains unresolved. This paper provides
a solution of this problem.

The underlying idea can be sketched as follows. Traditionally, given a set of
$n$ observations, each in $\mathbb{R}^{m}$, one would estimate the
distribution of these observations from the data and then check whether the
synthetic data can be viewed as coming from the data distribution. Here
instead we characterise the distribution which is generated possibly
implicitly from the synthetic data generator, and then test whether the
observed data can be viewed as coming from the synthetic data distribution.
The advantage of this approach is that while the observed sample size $n$ may
be fairly small, the synthetic data distribution can be estimated to any
desirable level of accuracy by generating a large number of samples. Similarly
to the works mentioned in the previous paragraph for goodness-of-fit tests, we
use a KSD approach, based on a Stein operator which characterises the
synthetic data distribution. As the synthetic data generator is usually
implicit, this Stein operator is not available. We show however that it can be
estimated from synthetic data samples to any desired level of accuracy.

#### Our contributions

We introduce a method to assess (deep) generative models, which are often
black-box approaches, when the underlying probability distribution is
continuous, usually in high-dimensions. To this purpose, we develop a non-
parametric Stein operator and the corresponding non-parametric kernel Stein
discrepancies (NP-KSD), based on estimating conditional score functions.
Moreover, we give theoretical guarantees for NP-KSD.

This paper is structured as follows. We start with a review of Stein’s method
and KSD goodness-of-fit tests for explicit models in Section 2 before we
introduce the NP-KSD in Section 3 and analyse the model assessment procedures.
We show results of experiments in Section 4 and conclude with future
directions in Section 5. Theoretical underpinnings, and additional results are
provided in the supplementary material. The code is available at
https://github.com/wenkaixl/npksd.git.

##  2 Stein’s method and kernel Stein discrepancy tests

#### Stein identities, equations, and operators

Stein’s method (Stein, 1972) provides an elegant tool to characterise
distributions via _Stein operators_ , which can be used to assess distances
between probability distributions (Barbour and Chen, 2005; Barbour, 2005;
Barbour et al., 2018). Given a distribution $q$, an operator $\mathcal{A}_{q}$
is called a Stein operator w.r.t. $q$ and Stein class $\mathcal{F}$ if the
following Stein identity holds for any _test function_ $f{\in\mathcal{F}}$:
${\mathbb{E}}_{q}[\mathcal{A}_{q}{f}]=0.$ For a test function $h$ one then
aims to find a function $f=f_{h}{\in\mathcal{F}}$ which solves the Stein
equation

$\mathcal{A}_{q}f(\boldsymbol{x})=h(\boldsymbol{x})-\mathbb{E}_{q}[h(\boldsymbol{x})].$
(1)

Then for any distribution $p$, taking expectations $\mathbb{E}_{p}$ in Eq. 1
assesses the distance $|\mathbb{E}_{p}h-\mathbb{E}_{q}h|$ through
$|\mathbb{E}_{p}\mathcal{A}_{q}f|$, an expression in which randomness enters
only through the distribution $p$.

When the density function $q$ is given explicitly, with smooth support
$\Omega_{q}{\subset\mathbb{R}^{m}}$, is differentiable and vanishes at the
boundary of $\Omega_{q}$, a common choice of Stein operator in the literature
utilises the score-function, see for example Mijoule et al. (2021). The
gradient operator is denoted by $\nabla$ and taken to be a column vector. The
_score function_ of $q$ is defined as $\mathbf{s}_{q}=\nabla\log
q=\frac{\nabla q}{q}$ (with the convention that $\mathbf{s}_{q}\equiv 0$
outside of $\Omega_{q}$). Let $\mathbf{f}=(f_{1},\dots,f_{{m}})^{\top}$ where
$f_{i}:\mathbb{R}^{{m}}\to\mathbb{R},\forall i,$ are differentiable. The
_score-Stein operator_ 111also referred to as Langevin Stein operator (Barp et
al., 2019). is the vector-valued operator acting on (vector-valued) function
$\mathbf{f}$,

$\displaystyle\mathcal{A}_{q}\mathbf{f}(\boldsymbol{x})$
$\displaystyle=\mathbf{f}(\boldsymbol{x})^{\top}\nabla\log
q(\boldsymbol{x})+\nabla\cdot\mathbf{f}(\boldsymbol{x}),$ (2)

and the Stein identity $\mathbb{E}_{q}[\mathcal{A}_{q}f]=0$ holds for
functions $f$ which belong to the so-called canonical Stein class defined in
Mijoule et al. (2021), Definition 3.2. As it requires knowledge of the density
$q$ only via its score function, this Stein operator is particularly useful
for unnormalised densities (Hyvärinen, 2005), appearing e.g. in energy based
models (EBM) (LeCun et al., 2006).

#### Kernel Stein discrepancy

Stein operators can be used to assess discrepancies between two probability
distributions; the Stein discrepancy between probability distribution $p$ and
$q$ (w.r.t. class ${\mathcal{B}}\subset{\mathcal{F}}$) is defined as (Gorham
and Mackey, 2015)

$\operatorname{SD}(p\|q,{{\mathcal{B}}})=\sup_{f\in\mathcal{B}}\\{|\mathbb{E}_{p}[{\mathcal{A}}_{q}f]-\underset{=0}{\underbrace{\mathbb{E}_{p}[\mathcal{A}_{p}f]}}|\\}=\sup_{f\in\mathcal{B}}{|}\mathbb{E}_{p}[{\mathcal{A}}_{q}f]{|}.$
(3)

As the $\sup f$ over a general class $\mathcal{B}$ can be difficult to
compute, taking $\mathcal{B}$ as the unit ball of a reproducing kernel Hilbert
space (RKHS) has been considered, resulting in the kernel Stein discrepancy
(KSD) defined as (Gorham and Mackey, 2017)

$\operatorname{KSD}(p\|q,{{\mathcal{H}}})=\sup_{f\in\mathcal{B}_{1}(\mathcal{H})}{|}\mathbb{E}_{p}[{\mathcal{A}}_{q}f]{|}.$
(4)

Denoting by $k$ the reproducing kernel associated with the RKHS $\mathcal{H}$
over a set $\mathcal{X}$, the reproducing property ensures that $\forall
f\in\mathcal{H}$, $f(\boldsymbol{x})=\langle
f,k(\boldsymbol{x},\cdot)\rangle_{\mathcal{H}},\forall\boldsymbol{x}\in\mathcal{X}$.
Algebraic manipulations yield

$\displaystyle\mathrm{KSD}^{2}(q\|p)={\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[u_{q}(\boldsymbol{x},\tilde{\boldsymbol{x}})],$ (5)

where
$u_{q}(\boldsymbol{x},\tilde{\boldsymbol{x}})=\langle\mathcal{A}_{q}k(\boldsymbol{x},\cdot),\mathcal{A}_{q}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$,
which takes the exact $\sup$ without approximation and does not involve the
(sample) distribution $p$. Then, KSD2 can be estimated through empirical
means, over samples from $p$, e.g. V-statistic (Van der Vaart, 2000) and
U-statistics (Lee, 1990) estimates are

${{\operatorname{KSD}}_{v}^{2}}(q\|p)=\frac{1}{m^{2}}\sum_{i{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},}j}u_{q}(\boldsymbol{x}_{i},\boldsymbol{x}_{j}),\qquad{{\operatorname{KSD}}_{u}^{2}}(q\|p)=\frac{1}{m(m-1)}\sum_{i\neq
j}u_{q}(\boldsymbol{x}_{i}\,\boldsymbol{x}_{j}).$ (6)

KSD has been studied as discrepancy measure between distributions for testing
model goodness-of-fit (Chwialkowski et al., 2016; Liu et al., 2016).

#### KSD testing procedure

Suppose we have observed samples $\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{n}$
from the _unknown_ distribution $p$. To test the null hypothesis
$\mathrm{H_{0}}:p=q$ against the (broad class of) alternative hypothesis
$\mathrm{H_{1}}:p\neq q$, KSD can be empirically estimated via Eq. 6. The null
distribution is usually simulated via the wild-bootstrap procedure
(Chwialkowski et al., 2014). Then if the empirical quantile, i.e. the
proportion of wild bootstrap samples that are larger than
${\operatorname{KSD}_{v}^{2}}(q\|p)$, is smaller than the pre-defined test
level (or significance level) $\alpha$, the null hypothesis is rejected;
otherwise the null hypothesis is not rejected. In this way, a systematic non-
parametric goodness-of-fit testing procedure is obtained, which is applicable
to unnormalised models.

##  3 Non-Parametric kernel Stein discrepancies

The construction of a KSD relies on the knowledge of the density model, up to
normalisation. However, for deep generative models where the density function
is not explicitly known, the computation for Stein operator in Eq. 2, which is
based on an explicit parametric density, is no longer feasible.

While in principle one could estimate the multivariate density function from
synthetic data, density estimation in high dimensions is known to be
problematic, see for example Scott and Sain (2005). Instead, Stein’s method
allows to use a two-step approach: For data in $\mathbb{R}^{m}$, we first pick
a coordinate $i\in[m]:=\\{1,\ldots,m\\}$, and then we characterize the uni-
variate conditional distribution of that coordinate, given the values of the
other coordinates. Using score Stein operators from Ley et al. (2017), this
approach only requires knowledge or estimation of uni-variate conditional
score functions.

We denote observed data ${\boldsymbol{z}}_{1},\ldots,{\boldsymbol{z}}_{n}$
with
${\boldsymbol{z}}_{i}=(z_{i}^{(1)},\ldots,z_{i}^{(m)})^{{\top}}\in\mathbb{R}^{m}$;
and denoting the generative model as $G$, we write ${\boldsymbol{X}}\sim G$ to
denote a random $\mathbb{R}^{m}$-valued element from the (often only given
implicitly) distribution which is underlying $G$. Using $G$, we generate $N$
samples denoted by ${\boldsymbol{y}}_{1},\ldots,{\boldsymbol{y}}_{N}$. In our
case, $n$ is fixed and $n\ll N$, allowing $N\rightarrow\infty$ in theoretical
results. The kernel of an RKHS is denoted by $k$ and is assumed to be bounded.
For $\boldsymbol{x}\in\mathbb{R}^{m}$, $x\in\mathbb{R}$ and
$g(\boldsymbol{x}):\mathbb{R}^{m}\rightarrow\mathbb{R}$, we write
$g_{x^{(-i)}}(x):\mathbb{R}\rightarrow\mathbb{R}$ for the uni-variate function
which acts only on the coordinate $i$ and fixes the other coordinates to equal
$x^{(j)},j\neq i$, so that
$g_{x^{(-i)}}(x)=g(x^{(1)},\ldots,x^{(i-1)},x,x^{(i+1)},\ldots,x^{(m)})$.

For $i\in[m]$ let $\mathcal{T}^{(i)}$ denote a Stein operator for the
conditional distribution $Q^{(i)}=Q^{(i)}_{x^{(-i)}}$ with
$\mathbb{E}_{Q^{(i)}_{x^{(-i)}}}g_{x^{(-i)}}(x)=\mathbbm{E}[g_{y^{(-i)}}(Y)|Y^{(j)}=y^{(j)},j\neq
i]$. The proposed Stein operator $\mathcal{A}$ acting on functions
$g:\mathbb{R}^{m}\rightarrow\mathbb{R}$ underlying the non-parametric Stein
operator is

$\mathcal{A}g(x^{{(1)}},\ldots,x^{{(m)}})=\frac{1}{m}{{\sum_{i=1}^{m}\mathcal{T}^{(i)}g_{x^{(-i)}}(x^{(i)})}}.$
(7)

We note that for ${\boldsymbol{X}}\sim q$, the Stein identity
$\mathbb{E}\mathcal{A}g({\boldsymbol{X}})=0$ holds and thus $\mathcal{A}$ is a
Stein operator. The domain of the operator will depend on the conditional
distribution in question. Instead of using the weights $w_{i}=\frac{1}{m}$,
other positive weights which sum to 1 would be possible, but for simplicity we
use equal weights. A more detailed theoretical justification of Eq. 7 is given
in Appendix A.

In what follows we use as Stein operator for a differentiable uni-variate
density $q$ the score operator from Eq. 2, given by

$\displaystyle\mathcal{T}_{q}^{(i)}f(x)=f^{\prime}(x)+f(x)\frac{q^{\prime}(x)}{q(x)}.$
(8)

In Proposition D.1 of Appendix D we shall see that the operator in Eq. 7
equals the score-Stein operator in Eq. 2; in Appendix D an example is also
given. For the development in this paper, Eq. 7 is more convenient as it
relates directly to conditional distributions. Other choices of Stein
operators are discussed for example in Ley et al. (2017); Mijoule et al.
(2021); Xu (2022).

#### Re-sampling Stein operators

The Stein operator Eq. 7 depends on all coordinates $i\ \in[m]$. When $m$ is
large we can estimate this operator via re-sampling with replacement, as
follows. We draw $B$ samples $\\{i_{1},\dots,i_{B}\\}$ with replacement from
$[m]$ such that
$\\{i_{1},\dots,i_{B}\\}\sim\operatorname{Multinom}({B},\\{\frac{1}{m}\\}_{i\in[m]})$.
The re-sampled Stein operator acting on
$f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is

$\mathcal{A}^{B}f({{{{\boldsymbol{z}}}}}):=\frac{1}{B}\sum_{b=1}^{B}\mathcal{A}^{(i_{b})}f({{{{\boldsymbol{z}}}}}).$
(9)

Then we have
$\mathbb{E}\mathcal{A}{{}^{B}}f({\boldsymbol{X}})=\frac{1}{B}\sum_{b=1}^{B}\mathbb{E}\mathcal{A}^{(i_{b})}f({\boldsymbol{X}})=0.$
So $\mathcal{A}{{}^{B}}$ is again a Stein operator.

In practice, when $m$ is large, the stochastic operator in Eq. 9 creates a
computationally efficient way for comparing distributions. A similar re-
sampling strategy for constructing stochastic operators are considered in the
context of Bayesian inference (Gorham et al., 2020), where conditional score
functions, which are given in parametric form, are re-sampled to derive score-
based (or Langevin) Stein operators for posterior distributions. The
conditional distribution has been considered (Wang et al., 2018) and (Zhuo et
al., 2018) in the context of graphical models (Liu and Wang, 2016). In
graphical models, the conditional distribution is simplified to conditioning
on the Markov blanket (Wang et al., 2018), which is a subset of the full
coordinate; however, no random re-sampling is used. Conditional distributions
also apply in message passing, but there, the sequence of updates is ordered.

#### Estimating Stein operators via score matching

Usually the score function $q^{\prime}/q$ in Eq. 8 is not available but needs
to be estimated. An efficient way of estimating the score function is through
score-matching, see for example (Hyvärinen, 2005; Song and Kingma, 2021;
Wenliang et al., 2019). Score matching relies on the following score-matching
(SM) objective (Hyvärinen, 2005),

$J(p\|q)=\mathbb{E}_{p}\left[\left\|\nabla\log p({\boldsymbol{x}})-\nabla\log
q({\boldsymbol{x}})\right\|^{2}\right],$ (10)

which is particularly useful for unnormalised models such as EBMs. Additional
details are included in Appendix E. Often score matching estimators can be
shown to be consistent, see for example Song et al. (2020). Proposition 3.1,
proven in Appendix B, gives theoretical guarantees for the consistency of a
general form of Stein operator estimation, as follows.

######  Proposition 3.1.

Suppose that for $i\in[m]$, $\widehat{s}_{N}^{(i)}$ is a consistent estimator
of the uni-variate score function $s^{(i)}$. Let $\mathcal{T}^{(i)}$ be a
Stein operator for the uni-variate differentiable probability distribution
${Q}^{(i)}$ of the generalised density operator form Eq. 8. Let

$\displaystyle{\widehat{\mathcal{T}}}^{(i)}_{N}g(x)=g^{\prime}(x)+g(x)\widehat{s}_{N}^{(i)}\quad\quad\operatorname{and}\quad\quad\widehat{\mathcal{A}}g=\frac{1}{m}\sum_{i}{\widehat{\mathcal{T}}}^{(i)}_{N}g_{x^{(-i)}}.$

Then ${\widehat{\mathcal{T}}}^{(i)}_{N}$ is a consistent estimator for
${{\mathcal{T}}}^{(i)}$, and $\widehat{\mathcal{A}}$ is a consistent estimator
of $\mathcal{A}.$

#### Non-parametric Stein operators with summary statistics

In practice, the data ${y}^{(-i)}\in\mathbb{R}^{m-1}$ can be high dimensional,
e.g. image pixels, and the observations can be sparse. Thus, estimation of the
conditional distribution can be unstable or exponentially large sample size is
required. Inspired by Xu and Reinert (2021) and Xu and Reinert (2022), we use
low-dimensional measurable non-trivial summary statistics $t$ and the
conditional distribution of the data given $t$ as new target distributions.
Heuristically, if two distributions match, then so do their conditional
distributions. Thus, the conditional distribution $Q{{}^{(i)}}(A)$ is replaced
by $Q_{t}^{(i)}(A)=\mathbb{P}(X^{(i)}\in A|t(x^{(-i)}))$. Setting
$t(x^{(-i)})=x^{(-i)}$ replicates the actual conditional distribution. We
denote the uni-variate score function of $q_{t}(x|t(x^{(-i)}))$ by
${s_{t}^{(i)}(x|t(x^{(-i)}}))$, or by $s_{t}^{(i)}(x)$ when the context is
clear. The summary statistics $t(x^{(-i)})$ can be uni-variate or multi-
variate, and they may attempt to capture useful distributional features. Here
we consider uni-variate summary statistics such as the sample mean.

0: Generator $G$; summary statistics $t(\cdot)$; number of samples $N$ from
$G$; re-sample size $B$

0:

1: Generate samples
$\\{{{\boldsymbol{y}}}_{1},\dots,{{\boldsymbol{y}}}_{N}\\}$ from $G$.

2: Generate coordinate index sample $\\{i_{1},\dots,i_{B}\\}$

3: For ${i_{b}\in[m],l\in[N]}$, estimate $q(z^{(i_{b})}|t(z^{-i_{b}})$ from
samples $\\{{y}^{(i_{b})}_{l},t({y_{l}}^{-i_{b}})\\}_{l\in[N]}$ via the score-
matching objective in Eq. 10.

3: $\widehat{s}_{t,N}^{(i)}(z^{(i)}|t(z^{(-i)})),\forall i\in[m]$.

Algorithm 1 Estimating the conditional probability via summary statistics

The non-parametric Stein operator enables the construction of Stein-based
statistics based on Eq. 7with estimated score functions
$\widehat{s}_{t,N}^{(i)}$ using generated samples from the model $G$, as shown
in Algorithm 1. The re-sampled non-parametric Stein operator is

$\widehat{\mathcal{A}^{B}_{t,N}}g=\frac{1}{B}\sum_{b}{\widehat{\mathcal{T}}}^{(i_{b})}_{t,N}g_{x^{(-i_{b})}}{=\frac{1}{B}\sum_{b}\left(g_{x^{(-i_{b})}}^{\prime}+g_{x^{(-i_{b})}}\widehat{s}_{t,N}^{(i)}\right)}.$

#### Non-parametric kernel Stein discrepancy

With the well-defined non-parametric Stein operator, we define the
corresponding non-parametric Stein discrepancy (NP-KSD) using the Stein
operator in Eq. 9, the Stein discrepancy notion in Eq. 3 and choosing as set
of test functions the unit ball of the RKHS within unit ball RKHS. Similarly
to Eq. 4, we define the NP-KSD with summary statistic $t$ as

$\operatorname{NP-
KSD}_{t}(G\|p)=\sup_{f\in\mathcal{B}_{1}(\mathcal{H})}\mathbb{E}_{p}[\widehat{\mathcal{A}}_{t,N}^{B}f].$
(11)

A similar quadratic form as in Eq. 5 applies to give

$\displaystyle{{\operatorname{NP-
KSD}}}_{t}^{2}(G\|p)={\mathbb{E}}_{{\boldsymbol{x}},{\tilde{\boldsymbol{x}}}\sim
p}[\widehat{u}^{B}_{t,N}({\boldsymbol{x}},{\tilde{\boldsymbol{x}}})],$ (12)

where
$\widehat{u}^{B}_{t,N}({\boldsymbol{x}},{\tilde{\boldsymbol{x}}})=\langle\widehat{\mathcal{A}}^{B}_{t,N}k(\boldsymbol{x},\cdot),\widehat{\mathcal{A}}^{B}_{t,N}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$.
The empirical estimate is

$\displaystyle\widehat{{\operatorname{NP-
KSD}}}_{t}^{2}(G\|p)=\frac{1}{n^{2}}\sum_{i,j\in[n]}[\widehat{u}^{B}_{t,N}({\boldsymbol{z}}_{i},{\boldsymbol{z}}_{j})],$
(13)

where $\mathbb{S}=\\{{\boldsymbol{z}}_{1},\dots,{\boldsymbol{z}}_{n}\\}\sim
p$. Thus, NP-KSD allows the computation between a set of samples and a
generative model, enabling the quality assessment of synthetic data generators
even for implicit models.

The relationship between NP-KSD and KSD is clarified in the following result;
we use the notation $\hat{\bf s}_{t,N}=(\hat{s}_{t,N}(x^{{(i)}}),i\in[m])$.
Here we set

${{\operatorname{KSD}}}_{t}^{2}(q_{t}\|p)={\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle\mathcal{A}_{t}k(\boldsymbol{x},\cdot),\mathcal{A}_{t}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}\quad\mbox{
with
}\quad\mathcal{A}_{t}g(\boldsymbol{x}):=\frac{1}{m}\sum_{i=1}^{m}\mathcal{T}_{q_{t}}^{(i)}g_{x^{(-i)}}(x^{(i)})$
(14)

as in Eq. 7, and following Eq. 8,
$\mathcal{T}_{q_{t}}^{(i)}g_{x^{(-i)}}(x)=g_{x^{(-i)}}^{\prime}(x)+g_{x^{(-i)}}(x)s_{t}^{(i)}(x|t(x^{(-i)})).$
More details about the interpretation of this quantity are given in App. B.1.

######  Theorem 3.2.

Assume that the score function estimator vector
$\hat{{\bf{s}}}_{t,N}{=(\hat{{s}}_{t,N}^{(i)},i=1,\dots,m)^{\top}}$ is
asymptotically normal with mean $0$ and covariance matrix $N^{-1}\Sigma_{s}$.
Then ${{\operatorname{NP-KSD}}}_{t}^{2}(G\|p)$ converges in probability to
${{\operatorname{KSD}}}_{t}^{2}(q_{t}\|p)$ at rate at least
$\min(B^{-\frac{1}{2}},N^{-\frac{1}{2}})$.

The proof of Theorem 3.2, which is found in App. B, also shows that the
distribution ${{\operatorname{NP-
KSD}}}_{t}^{2}(G\|p)-{{\operatorname{KSD}}}_{t}^{2}(q_{t}\|p)$ involves
mixture of normal variables. The assumption of asymptotic normality for score
matching estimators is often satisfied, see for example Song et al. (2020).

#### Model assessment with NP-KSD

Given an implicit generative model $G$ and a set of observed samples
$\mathbb{S}=\\{{{\boldsymbol{z}}}_{1},\ldots,{{\boldsymbol{z}}}_{n}\\}$, we
aim to test the null hypothesis $H_{0}:\mathbb{S}\sim G$ versus the
alternative $H_{1}:\mathbb{S}\not\sim G$. This test assumes that samples
generated from $G$ follows some (unknown) distribution $q$ and $\mathbb{S}$
are generated according to some (unknown) distribution $p$. The null
hypothesis is $H_{0}:p=q$ while the alternative is $H_{1}:p\neq q$. We note
that the observed sample size $n$ is fixed.

#### NP-KSD testing procedures

NP-KSD can be applied for testing the above hypothesis using the testing
procedure outlined in Algorithm 2. In contrast to the KSD testing procedure in
Section 2, the NP-KSD test in Algorithm 2 is a Monte Carlo based test (Xu and
Reinert, 2021, 2022; Schrab et al., 2022) for which the null distribution is
approximated via samples generated from $G$ instead of the wild bootstrap
procedure (Chwialkowski et al., 2014). The reasons for employing the Monte
Carlo testing strategy instead of the wild-bootstrap are 1). The non-
parametric Stein operator depends on the random function $\widehat{s}_{t}$ so
that classical results for V-statistics convergence which assume that the sole
source of randomness is the bootstrap may not apply222A KSD with random Stein
kernel has been briefly discussed in Fernández et al. (2020) when the $h_{q}$
function requires estimation from relevant survival functions.; 2). While the
wild-bootstrap is asymptotically consistent as observed sample size
$n\rightarrow\infty$, it may not necessarily control the type-I error in a
non-asymptotic regime where $n$ is fixed. More details can be found in
Appendix F.

Here we note that any test which is based on the summary statistic $t$ will
only be able to test for a distribution up to equivalence of their
distributions with respect to the summary statistic $t$; two distributions $P$
and $Q$ are equivalent w.r.t. the summary statistics $t$ if
$P({\boldsymbol{X}}|t({\boldsymbol{X}}))=Q({\boldsymbol{X}}|t({\boldsymbol{X}}))$.
Thus the null hypothesis for the NP-KSD test is that the distribution is
equivalent to $P$ with respect to $t$. Hence, the null hypothesis specifies
the conditional distribution, not the unconditional distribution.

0: Observed sample set
$\mathbb{S}=\\{{{\boldsymbol{z}}}_{1},\ldots,{{\boldsymbol{z}}}_{n}\\}$;
generator $G$ and generated sample size $N$; estimation statistics $t$; RKHS
kernel $K$; re-sampling size $B$; bootstrap sample size $b$; confidence level
$\alpha$;

1: Estimate $\widehat{s}(z^{(i)}|t(z^{(-i)}))$ based on Algorithm 1.

2: Uniformly generate re-sampling index $\\{i_{1},\dots,i_{B}\\}$ from $[m]$,
with replacement.

3: Compute $\tau=\widehat{\operatorname{NP-
KSD}}^{2}(\widehat{s}_{t};\mathbb{S})$ in Eq. (13).

4: Simulate
$\mathbb{S}_{i}=\\{{{\boldsymbol{y}}}^{\prime}_{1},\dots,{{\boldsymbol{y}}}^{\prime}_{n}\\}$
for ${i}\in[b]$ from $G$.

5: Compute $\tau_{i}=\widehat{\operatorname{NP-
KSD}}^{2}(\widehat{s}_{t};\mathbb{S}_{i})$ in again with index re-sampling.

6: Estimate the empirical (1- $\alpha$) quantile $\gamma_{1-\alpha}$ via
$\\{\tau_{1},\dots,\tau_{b}\\}$.

6: Reject the null hypothesis if $\tau>\gamma_{1-\alpha}$; otherwise do not
reject.

Algorithm 2 Assessment procedures for implicit generative models

#### Related works

To assess whether an implicit generative models can generate samples that are
_significantly_ good for the desired data model, several hypothesis testing
procedures have been studied. Jitkrittum et al. (2018) has proposed kernel-
based test statistics, Relative Unbiased Mean Embedding (Rel-UME) test and
Relative Finite-Set Stein Discrepancy (Rel-FSSD) test for relative model
goodness-of-fit, i.e. whether model S is a better fit than model R. While Rel-
UME is applicable for implicit generative models, Rel-FSSD still requires
explicit knowledge of the unnormalised density. The idea for assessing sample
quality for implicit generative models is through addressing two-sample
problem, where samples generated from the implicit model are compared with the
observed data. In this sense, maximum-mean-discrepancy (MMD) may also apply
for assessing sample qualities for the implicit models. With efficient choice
of (deep) kernel, Liu et al. (2020) applied MMD tests to assess the
distributional difference for image data, e.g. MNIST (LeCun et al., 1998) v.s.
digits image trained via deep convolutional GAN (DCGAN) (Radford et al.,
2015); CIFAR10 (Krizhevsky, 2009) v.s. CIFAR10.1 (Recht et al., 2019).
However, as the distribution is represented via samples, the two-sample based
assessment suffers from limited probabilistic information from the implicit
model and low estimation accuracy when the sample size for observed data is
small.

##  4 Experiments

###  4.1 Baseline and competing approaches

We illustrate the proposed NP-KSD testing procedure with different choice of
summary statistics. We denote by NP-KSD the version which uses the estimation
of the conditional score, i.e. $t(x^{(-i)})=x^{(-i)}$; by NP-KSD_mean the
version which uses conditioning on the mean statistics, i.e.
$t(x^{(-i)})=\frac{1}{m-1}\sum_{j\neq i}x^{(j)}$; and by NP-KSD_G the version
which fits a Gaussian model as conditional density333NP-KSD_G for non-Gaussian
densities is generally mis-specified. We deliberately check this case to
assess the robustness of the NP-KSD procedure under model mis-specification..

Two-sample testing methods can be useful for model assessment, where the
observed sample set is tested against sample set generated from the model. In
our setting where $n\ll N$, we consider a consistent non-asymptotic MMD-based
test, MMDAgg (Schrab et al., 2021), as our competing approach; see Appendix F
for more details. For synthetic distributions where the null models have
explicit densities, we include the KSD goodness-of-fit testing procedure in
Section 2 as the baseline. Gaussian kernels are used and the median heuristic
(Gretton et al., 2007) is applied for bandwidth selection. As a caveat, in
view of (Gorham and Mackey, 2015), when the kernel decays more rapidly than
the score function grows, then identifiability of $q_{t}$ through a KSD method
may not be guaranteed. Details while MMD is not included in this list are
found in Appendix F.

###  4.2 Experiments on synthetic distributions

(a) GVD: $n=100$

(b) GVD:$\sigma_{per}=0.4$

(c) MoG: $n=200$

(d) MoG: $m=40$

Figure 1: Rejection rates of the synthetic distributions: test level
$\alpha=0.05$; $100$ trials per round of experiment; $10$ rounds of experiment
are taken for average and standard deviation; bootstrap sample size $b=500$;
$m=3$ for (a) and (b); $m=6$ for (c); $n=100$, $\sigma_{per}=0.5$ for (d).

#### Gaussian Variance Difference (GVD)

We first consider a standard synthetic setting, studied in Jitkrittum et al.
(2017), in which the null distribution is multivariate Gaussian with mean zero
and identity covariance matrix. The alternative is set to perturb the the
diagonal terms of the covariance matrix, i.e. the variances, all by the same
amount.

The rejection rate against the variances perturbation is shown in Figure 1(a).
From the result, we see that all the tests presented have controlled type-I
error. For all the tests the power increases with increased perturbation. NP-
KSD and NP-KSD_mean outperform the MMDAgg approach. Using the mean statistics,
NP-KSD_mean is having slightly higher power than KSD. The mis-specified NP-
KSD_G has lower power, but is still competitive to MMDAgg.

The test power against the sample size $N$ generated from the null model is
shown in Figure 1(b). The generated samples are used as another sample set for
the MMDAgg two-sample procedure, while used for estimating the conditional
score for NP-KSD-based methods. As the generated sample size increases, the
power of MMDAgg increases more slowly than that of the NP-KSD-based methods,
which achieve maximum test power in the presented setting. The NP-KSD-based
tests tend to have lower variability of the test power, indicating more
reliable testing procedures than MMDAgg.

#### Mixture of Gaussian (MoG)

Next, we consider as a more difficult problem that the null model is a two-
component mixture of two independent Gaussians. Both Gaussian components have
identity covariance matrix. The alternative is set to perturb the covariance
between adjacent coordinates.

The rejection rate against this perturbation of covariance terms are presented
in Figure 1(c). The results show consistent type I error. The NP-KSD and NP-
KSD_mean tests have better test power compared to KSD and MMDAgg, although NP-
KSD has slightly higher variance. Among the ${\operatorname{NP-KSD}}$ tests,
the smallest variability is achieved by NP-KSD_mean. For the test with $m=40$,
we also vary the re-sample size $B$. As shown in Figure 1(d), while the
variability of the average test power also increased slightly. From the
result, we also see that for $B=20=m/2$ the test power is already competive
compared to $B=40$. Additional experimental results including computational
runtime and training generative models for synthetic distributions are
included in Appendix C.

###  4.3 Applications to deep generative models

For real-world applications, we assess models trained from well-studied
generative modelling procedures, including a Generative Adversarial Network
(GAN) (Goodfellow et al., 2014) with multilayer perceptron (MLP), a Deep
Convolutional Generative Adversarial Network (DCGAN) (Radford et al., 2015),
and a Variational Autoencoder (VAE) (Kingma and Welling, 2013). We also
consider a Noise Conditional Score Network (NCSN) (Song and Ermon, 2020),
which is a score-based generative modelling approach, where the score
functions are learned (Song and Ermon, 2019) to performed annealed Langevin
dynamics for sample generation. We also denote Real as the scheme that
generates samples randomly from the training data, which essentially acts as a
generator of the null distribution.

#### MNIST Dataset

This dataset contains $28\times 28$ grey-scale images of handwritten digits
(LeCun et al., 1998)444
https://pytorch.org/vision/main/generated/torchvision.datasets.MNIST.html . It
consist of $60,000$ training samples and $10,000$ test samples. Deep
generative models in Table 1 are trained using the training samples. We assess
the quality of these trained generative models by testing against the true
observed MNIST samples (from the test set). Samples from both distributions
are visually illustrated in Figure 3 in Appendix C.

| GAN_MLP | DCGAN | VAE | NCSN | Real
---|---|---|---|---|---
NP-KSD | 1.00 | 0.92 | 1.00 | 1.00 | 0.03
NP-KSD_m | 1.00 | 1.00 | 1.00 | 1.00 | 0.01
MMDAgg | 1.00 | 0.73 | 0.93 | 1.00 | 0.06
Table 1: Rejection rate for MNIST generative models.

$600$ samples are generated from the generative models and $100$ samples are
used for the test; test level $\alpha=0.05$. From Table 1, we see that all the
deep generative models have high rejection rate, showing that the trained
models are not good enough. Testing with the Real scheme has controlled type-I
error. Thus, NP-KSD detects that the “real” data are a true sample set from
the underlying dataset.

#### CIFAR10 Dataset

This dataset contains $32\times 32$ RGB coloured images (Krizhevsky, 2009)555
https://pytorch.org/vision/stable/generated/torchvision.datasets.CIFAR10.html
. It consist of $50,000$ training samples and $10,000$ test samples. Deep
generative models in Table 2 are trained using the training samples and test
samples are randomly drawn from the test set.

| DCGAN | NCSN | CIFAR10.1 | Real
---|---|---|---|---
NP-KSD | 0.68 | 0.73 | 0.92 | 0.06
NP-KSD_m | 0.74 | 0.81 | 0.96 | 0.02
MMDAgg | 0.48 | 0.57 | 0.83 | 0.07
Table 2: Rejection rate for CIFAR10 generative models.

Samples are illustrated in Figure 4 in Appendix C. We also compare with the
CIFAR10.1 dataset(Recht et al., 2018)666
https://github.com/modestyachts/CIFAR-10.1/tree/master/datasets, which is
created to differ from CIFAR10 to investigate generalisation power for
training classifiers. $800$ samples are generated from the generative models
and $200$ samples are used for the test; test level $\alpha=0.05$. Table 2
shows higher rejection rates for NP-KSD tests compared to MMDAgg, echoing the
results for synthetic distributions. The trained DCGAN generates samples with
lower rejection rate in the CIFAR10 dataset than in the CIFAR10.1 dataset. We
also see that the score-based NCSN has higher rejection rate than the non-
score-based DCGAN, despite NP-KSD being a score-based test. The distribution
difference between CIFAR10 and CIFAR10.1 can be well-distinguished from the
tests. Testing with the Real scheme again has controlled type-I error.

##  5 Conclusion and future directions

Synthetic data are in high demand, for example for training ML procedures;
quality is important. Synthetic data which miss important features in the data
can lead to erroneous conclusions, which in the case of medical applications
could be fatal, and in the case of loan applications for example could be
detrimental to personal or business development. NP-KSD provides a method for
assessing synthetic data generators which comes with theoretical guarantees.
Our experiments on synthetic data have shown that NP-KSD achieves good test
power and controlled type-I error. On real data, NP-KSD detects samples from
the true dataset. That none of the classical deep learning methods used in
this paper has a satisfactory rejection rate indicates scope for further
developments in synthetic data generation.

Future research will assess alternatives to the computer-intensive Monte Carlo
method for estimating the null distribution, for example adapting wild-
bootstrap procedures. It will explore alternative choices of score estimation
as well as of kernel functions.

Finally, some caution is advised. The choice of summary statistic may have
strong influence on the results and a classification based on NP-KSD may still
miss some features. Erroneous decisions could be reached when training
classifiers. Without scrutiny this could lead to severe consequences for
example in health science applications. Yet NP-KSD is an important step
towards understanding black-box data generating methods and thus understanding
their potential shortcomings.

## References

  * Barbour (2005) AD Barbour.  Multivariate Poisson–binomial approximation using Stein’s method.  In _Stein’s Method And Applications_ , pages 131–142. World Scientific, 2005.
  * Barbour and Chen (2005) AD Barbour and LHY Chen.  An Introduction to Stein‘s method.  _Lecture Notes Series. Institute for Mathematical Sciences. National University of Singapore_ , 4, 2005.
  * Barbour et al. (2018) AD Barbour, Malwina J Luczak, and Aihua Xia.  Multivariate approximation in total variation, ii: Discrete normal approximation.  _The Annals of Probability_ , 46(3):1405–1440, 2018.
  * Barp et al. (2019) Alessandro Barp, Francois-Xavier Briol, Andrew Duncan, Mark Girolami, and Lester Mackey.  Minimum Stein discrepancy estimators.  In _Advances in Neural Information Processing Systems_ , pages 12964–12976, 2019.
  * Bińkowski et al. (2018) Mikołaj Bińkowski, Danica J Sutherland, Michael Arbel, and Arthur Gretton.  Demystifying MMD GANs.  In _International Conference on Learning Representations_ , 2018.
  * Bishop and Nasrabadi (2006) Christopher M Bishop and Nasser M Nasrabadi.  _Pattern Recognition and Machine Learning_.  Springer, 2006.
  * Chwialkowski et al. (2016) Kacper Chwialkowski, Heiko Strathmann, and Arthur Gretton.  A kernel test of goodness of fit.  In _International Conference on Machine Learning_ , pages 2606–2615. PMLR, 2016.
  * Chwialkowski et al. (2014) Kacper P Chwialkowski, Dino Sejdinovic, and Arthur Gretton.  A wild bootstrap for degenerate kernel tests.  In _Advances in Neural Information Processing Systems_ , pages 3608–3616, 2014.
  * Creswell et al. (2018) Antonia Creswell, Tom White, Vincent Dumoulin, Kai Arulkumaran, Biswa Sengupta, and Anil A Bharath.  Generative adversarial networks: An overview.  _IEEE Signal Processing Magazine_ , 35(1):53–65, 2018.
  * Fernández et al. (2020) Tamara Fernández, Wenkai Xu, Marc Ditzhaus, and Arthur Gretton.  A kernel test for quasi-independence.  _Advances in Neural Information Processing Systems_ , 33:15326–15337, 2020.
  * Goodfellow et al. (2014) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  _Advances in neural information processing systems_ , 27, 2014.
  * Gorham and Mackey (2015) Jackson Gorham and Lester Mackey.  Measuring sample quality with Stein’s method.  In _Advances in Neural Information Processing Systems_ , pages 226–234, 2015.
  * Gorham and Mackey (2017) Jackson Gorham and Lester Mackey.  Measuring sample quality with kernels.  In _International Conference on Machine Learning_ , pages 1292–1301. PMLR, 2017.
  * Gorham et al. (2020) Jackson Gorham, Anant Raj, and Lester Mackey.  Stochastic Stein discrepancies.  _Advances in Neural Information Processing Systems_ , 33:17931–17942, 2020.
  * Gretton et al. (2008) A. Gretton, K. Fukumizu, C.H. Teo, L. Song, B. Schölkopf, and A.J. Smola.  A kernel statistical test of independence.  In _Advances in Neural Information Processing Systems_ , pages 585–592, 2008.
  * Gretton et al. (2007) Arthur Gretton, Karsten Borgwardt, Malte Rasch, Bernhard Schölkopf, and Alex J Smola.  A kernel method for the two-sample-problem.  In _Advances in Neural Information Processing Systems_ , pages 513–520, 2007.
  * Gretton et al. (2009) Arthur Gretton, Kenji Fukumizu, Zaid Harchaoui, and Bharath K Sriperumbudur.  A fast, consistent kernel two-sample test.  In _Advances in Neural Information Processing Systems_ , pages 673–681, 2009.
  * Gretton et al. (2012a) Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola.  A kernel two-sample test.  _The Journal of Machine Learning Research_ , 13(1):723–773, 2012a.
  * Gretton et al. (2012b) Arthur Gretton, Dino Sejdinovic, Heiko Strathmann, Sivaraman Balakrishnan, Massimiliano Pontil, Kenji Fukumizu, and Bharath K Sriperumbudur.  Optimal kernel choice for large-scale two-sample tests.  In _Advances in Neural Information Processing Systems_ , pages 1205–1213, 2012b.
  * Heusel et al. (2017) Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.  GANs trained by a two time-scale update rule converge to a local Nash equilibrium.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Holmes and Reinert (2004) Susan Holmes and Gesine Reinert.  Stein’s method for the bootstrap.  In _Stein’s Method_ , volume 46, pages 93–133. Institute of Mathematical Statistics, 2004.
  * Hyvärinen (2005) Aapo Hyvärinen.  Estimation of non-normalized statistical models by score matching.  _Journal of Machine Learning Research_ , 6(Apr):695–709, 2005.
  * Jitkrittum et al. (2016) Wittawat Jitkrittum, Zoltán Szabó, Kacper P Chwialkowski, and Arthur Gretton.  Interpretable distribution features with maximum testing power.  In _Advances in Neural Information Processing Systems_ , pages 181–189, 2016.
  * Jitkrittum et al. (2017) Wittawat Jitkrittum, Wenkai Xu, Zoltán Szabó, Kenji Fukumizu, and Arthur Gretton.  A linear-time kernel goodness-of-fit test.  In _Advances in Neural Information Processing Systems_ , pages 262–271, 2017.
  * Jitkrittum et al. (2018) Wittawat Jitkrittum, Heishiro Kanagawa, Patsorn Sangkloy, James Hays, Bernhard Schölkopf, and Arthur Gretton.  Informative features for model comparison.  In _Advances in Neural Information Processing Systems_ , pages 808–819, 2018.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kingma and Welling (2013) Diederik P Kingma and Max Welling.  Auto-Encoding Variational Bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Krizhevsky (2009) Alex Krizhevsky.  Learning multiple layers of features from tiny images.  _Technical Report, University of Toronto_ , 2009.
  * LeCun et al. (1995) Yann LeCun, Lawrence D Jackel, Léon Bottou, Corinna Cortes, John S Denker, Harris Drucker, Isabelle Guyon, Urs A Müller, Eduard Säckinger, Patrice Simard, and Vladimir Vapnik.  Learning algorithms for classification: A comparison on handwritten digit recognition.  _Neural networks: the statistical mechanics perspective_ , 261(276):2, 1995.
  * LeCun et al. (1998) Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * LeCun et al. (2006) Yann LeCun, Sumit Chopra, Raia Hadsell, M Ranzato, and F Huang.  A tutorial on energy-based learning.  _Predicting Structured Data_ , 1(0), 2006.
  * Lee (1990) A. J. Lee.  _U-Statistics: Theory and Practice_.  CRC Press, 1990.
  * Ley et al. (2017) Christophe Ley, Gesine Reinert, and Yvik Swan.  Stein’s method for comparison of univariate distributions.  _Probability Surveys_ , 14:1–52, 2017.
  * Li et al. (2017) Chun-Liang Li, Wei-Cheng Chang, Yu Cheng, Yiming Yang, and Barnabás Póczos.  MMD GAN: Towards deeper understanding of moment matching network.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Liu et al. (2020) Feng Liu, Wenkai Xu, Jie Lu, Guangquan Zhang, Arthur Gretton, and Danica J Sutherland.  Learning deep kernels for non-parametric two-sample tests.  In _International Conference on Machine Learning_ , pages 6316–6326. PMLR, 2020.
  * Liu and Wang (2016) Qiang Liu and Dilin Wang.  Stein variational gradient descent: A general purpose Bayesian inference algorithm.  In _Advances In Neural Information Processing Systems_ , pages 2378–2386, 2016.
  * Liu et al. (2016) Qiang Liu, Jason Lee, and Michael Jordan.  A kernelized Stein discrepancy for goodness-of-fit tests.  In _International Conference on Machine Learning_ , pages 276–284, 2016.
  * Mijoule et al. (2021) Guillaume Mijoule, Gesine Reinert, and Yvik Swan.  Stein’s density method for multivariate continuous distributions.  _arXiv preprint arXiv:2101.05079_ , 2021.
  * Ouimet (2021) Frédéric Ouimet.  General formulas for the central and non-central moments of the multinomial distribution.  _Stats_ , 4(1):18–27, 2021.
  * Radford et al. (2015) Alec Radford, Luke Metz, and Soumith Chintala.  Unsupervised representation learning with deep convolutional generative adversarial networks.  _arXiv preprint arXiv:1511.06434_ , 2015.
  * Recht et al. (2018) Benjamin Recht, Rebecca Roelofs, Ludwig Schmidt, and Vaishaal Shankar.  Do cifar-10 classifiers generalize to cifar-10?  _arXiv preprint arXiv:1806.00451_ , 2018.
  * Recht et al. (2019) Benjamin Recht, Rebecca Roelofs, Ludwig Schmidt, and Vaishaal Shankar.  Do imagenet classifiers generalize to imagenet?  In _International Conference on Machine Learning_ , pages 5389–5400. PMLR, 2019.
  * Reinert (2005) Gesine Reinert.  Three general approaches to Stein’s method.  _An introduction to Stein’s method_ , 4:183–221, 2005\.
  * Salimans et al. (2016) Tim Salimans, Ian Goodfellow, Wojciech Zaremba, Vicki Cheung, Alec Radford, and Xi Chen.  Improved techniques for training GANs.  _Advances in Neural Information Processing Systems_ , 29, 2016.
  * Schrab et al. (2021) Antonin Schrab, Ilmun Kim, Mélisande Albert, Béatrice Laurent, Benjamin Guedj, and Arthur Gretton.  MMD aggregated two-sample test.  _arXiv preprint arXiv:2110.15073_ , 2021.
  * Schrab et al. (2022) Antonin Schrab, Benjamin Guedj, and Arthur Gretton.  KSD aggregated goodness-of-fit test.  _arXiv preprint arXiv:2202.00824_ , 2022.
  * Scott and Sain (2005) David W Scott and Stephan R Sain.  Multidimensional density estimation.  _Handbook of statistics_ , 24:229–261, 2005.
  * Song and Ermon (2019) Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  In _Advances in Neural Information Processing Systems_ , pages 11895–11907, 2019.
  * Song and Ermon (2020) Yang Song and Stefano Ermon.  Improved techniques for training score-based generative models.  _Advances in Neural Information Processing Systems_ , 33:12438–12448, 2020.
  * Song and Kingma (2021) Yang Song and Diederik P Kingma.  How to train your energy-based models.  _arXiv preprint arXiv:2101.03288_ , 2021.
  * Song et al. (2020) Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon.  Sliced score matching: A scalable approach to density and score estimation.  In _Uncertainty in Artificial Intelligence_ , pages 574–584. PMLR, 2020.
  * Sriperumbudur et al. (2011) Bharath K Sriperumbudur, Kenji Fukumizu, and Gert RG Lanckriet.  Universality, characteristic kernels and RKHS embedding of measures.  _Journal of Machine Learning Research_ , 12(Jul):2389–2410, 2011.
  * Stein (1972) Charles Stein.  A bound for the error in the normal approximation to the distribution of a sum of dependent random variables.  In _Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 2: Probability Theory_. The Regents of the University of California, 1972.
  * Theis et al. (2015) Lucas Theis, Aäron van den Oord, and Matthias Bethge.  A note on the evaluation of generative models.  _arXiv preprint arXiv:1511.01844_ , 2015.
  * Tolstikhin et al. (2017) Ilya Tolstikhin, Bharath K Sriperumbudur, and Krikamol Muandet.  Minimax estimation of kernel mean embeddings.  _The Journal of Machine Learning Research_ , 18(1):3002–3048, 2017.
  * Van der Vaart (2000) Aad W Van der Vaart.  _Asymptotic Statistics_ , volume 3.  Cambridge University Press, 2000.
  * Wang et al. (2018) Dilin Wang, Zhe Zeng, and Qiang Liu.  Stein variational message passing for continuous graphical models.  In _International Conference on Machine Learning_ , pages 5219–5227. PMLR, 2018.
  * Wang et al. (2004) Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli.  Image quality assessment: from error visibility to structural similarity.  _IEEE Transactions on Image Processing_ , 13(4):600–612, 2004.
  * Welling and Teh (2011) Max Welling and Yee W Teh.  Bayesian learning via stochastic gradient langevin dynamics.  In _Proceedings of the 28th international conference on machine learning (ICML-11)_ , pages 681–688. Citeseer, 2011.
  * Wenliang et al. (2019) Li Wenliang, Danica J Sutherland, Heiko Strathmann, and Arthur Gretton.  Learning deep kernels for exponential family densities.  In _International Conference on Machine Learning_ , pages 6737–6746. PMLR, 2019.
  * Xu (2022) Wenkai Xu.  Standardisation-function kernel Stein discrepancy (Sf-KSD): A unifying approach for non-parametric goodness-of-fit testing.  In _International Conference on Artificial Intelligence and Statistics_ , 2022.
  * Xu and Reinert (2021) Wenkai Xu and Gesine Reinert.  A Stein goodness-of-test for exponential random graph models.  In _International Conference on Artificial Intelligence and Statistics_ , pages 415–423. PMLR, 2021.
  * Xu and Reinert (2022) Wenkai Xu and Gesine Reinert.  AgraSSt: Approximate graph Stein statistics for interpretable assessment of implicit graph generators.  _arXiv preprint arXiv:2203.03673_ , 2022.
  * Zhou et al. (2020) Yuhao Zhou, Jiaxin Shi, and Jun Zhu.  Nonparametric score estimators.  In _International Conference on Machine Learning_ , pages 11513–11522. PMLR, 2020.
  * Zhuo et al. (2018) Jingwei Zhuo, Chang Liu, Jiaxin Shi, Jun Zhu, Ning Chen, and Bo Zhang.  Message passing Stein variational gradient descent.  In _International Conference on Machine Learning_ , pages 6018–6027. PMLR, 2018.

##  Appendix A Justification of the Stein operator

Here we justify the two-step approach for constructing a Stein operator.

###  A.1 Step 1: A non-parametric Stein operator

Suppose we can estimate the conditional distribution from data. Then we can
create a Markov chain with values in $(\mathbb{R}^{d})^{m}$ as follows.
Starting with $Z_{0}=\\{x_{1},\dots,x_{m}\\}$ with $x_{i}\in\mathbb{R}^{d}$
for $i=1,\ldots,m$ (often we choose $d=1$), we pick an index $I\in[m]$ at
random. If $I=i$ we replace $x_{i}$ by $X_{i}^{\prime}$ drawn from the
conditional distribution of $X_{i}$ given $(X_{j}:j\neq i)$. This gives
$Z_{1}=(x_{1},\ldots,x_{i-1},X_{i}^{\prime},x_{i+1},\ldots,x_{m})$777Denote
$Z_{1}=(x^{(1)},\ldots,x^{(i-1)},{X^{(i)}}^{\prime},x^{(i+1)},\ldots,x^{(m)})\in\mathbb{R}^{m}$
where the superscript $(i)$ is used for coordinate index. The subscript is
used to denote different samples.; see for example Reinert [2005]. To make
this a continuous-time Markov process generator, we wait an
exponential(1)-distributed time before every change.

This generator induces a Stein operator for the target distribution as
follows. Here we take $d=1$ for clarity; the generalisation to other $d$ is
straightforward. Let $f:{\mathbb{R}^{m}}\to\mathbb{R}$ and consider the
expectation w.r.t. the one-step evolution of the Markov chain

$\displaystyle\mathbb{E}_{-i}[f(x^{(1)},\ldots,x^{(i-1)},X^{(i)},x^{(i+1)},\ldots,x^{(m)})]$
$\displaystyle=\int
f(x^{(1)},\ldots,x^{(i-1)},y,x^{(i+1)},\ldots,x^{(m)})\mathbb{P}(X^{(i)}=y|X^{(j)}=x^{(j)},j\neq
i).$

We now consider the index $i$ as the $i$-th coordinate of multivariate random
variables in $\mathbb{R}^{m}$. The conditional expectation here fixing all but
the $i$-th _coordinate_ term only depends on the uni-variate conditional
distribution $Q{{}^{(i)}}$ with $Q{{}^{(i)}}(A)=\mathbb{P}(X^{(i)}\in
A|X^{(j)}=x^{(j)},j\neq i)$. Thus, the Stein operator induced from the Markov
chain has the form

$\mathcal{A}f(z)=\mathcal{A}^{(I)}f(z)$ (15)

where

$\mathcal{A}^{(i)}f({{\boldsymbol{x}}})=\mathbb{E}_{-i}[f(x^{(1)},\ldots,x^{(i-1)},X^{(i)},x^{(i+1)},\ldots,x^{(m)})]-f({{\boldsymbol{x}}}).$
(16)

From the law of total expectation it follows that the Stein identity holds.

###  A.2 Step 2: marginal Stein operators

In Eq. (16), the expectation

$\mathbb{E}_{-i}[f(x^{{(i)}},\ldots,x^{{(i-1)}},X^{{(i)}},x^{{(i+1)}},\ldots,x^{{(m)}})]-f(x^{{(1)}},\ldots,x^{{(m)}})$

can itself be treated via Stein’s method, by finding a Stein operator
$\mathcal{T}^{{(i)}}$ and a function $g$ such that $g=g_{f}$ solves the
$\mathcal{T}^{(i)}$-Stein equation Eq. (1) for $f$;

$\mathcal{T}^{(i)}g(x)=\mathbb{E}_{-i}[f(x^{(1)},\ldots,x^{(i-1)},X^{(i)},x^{(i+1)},\ldots,x^{(m)})]-f(x^{(1)},\ldots,x^{(m)}).$
(17)

Fixing $x_{j},j\neq i$ and setting
$f^{(i)}(x)=f(x^{(1)},\ldots,x^{(i-1)},x,x^{(i+1)},\ldots,x^{(m)})$, we view
$\mathcal{T}^{(i)}$ as a Stein operator for a uni-variate distribution, acting
on functions $g=g_{x^{(-i)}}:\mathbb{R}\rightarrow\mathbb{R}$.

Summarising the approach, the Stein operator $\mathcal{A}$ acting on functions
$f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ underlying the non-parametric Stein
operator is

$\mathcal{A}f(x^{(1)},\ldots,x^{(m)})=\mathcal{T}^{(I)}g_{f,x^{-I}}(x^{(I)})$
(18)

where $I\in[m]$ is a randomly chosen index. In view of Eq. (18) we take
$g:\mathbb{R}^{m}\rightarrow\mathbb{R}$, write
$g_{x^{(-i)}}(x):\mathbb{R}\rightarrow\mathbb{R}$ for the uni-variate function
which acts only on the coordinate $i$ and fixes the other coordinates to equal
$x^{(-i)}$, we as Stein operator (using the same letter $\mathcal{A}$ as
before, which is abuse of notation);

$\mathcal{A}g(x^{(1)},\ldots,x^{(m)})=\mathcal{T}^{(I)}g_{x^{-I}}(x^{(I)}).$

This formulation simplifies Eq. (18) in that we no longer have to consider the
connection between $f$ an $g$. The final step is to note that when we
condition on the random index $I$, again a Stein operator is obtained, as
follows. As

$\mathbbm{E}_{I}[\mathcal{A}g(x^{(1)},\ldots,x^{(m)})]=\frac{1}{m}\sum_{i=1}^{m}\mathcal{T}^{(i)}g_{x^{(-i)}}(x^{(i)}).$
(19)

As $\mathbb{E}[\mathcal{T}^{(i)}g_{X^{(-i)}}(X^{(i)})]=0$, the Stein identity
is satisfied. The operator in Eq. (19) is the Stein operator given in Eq. (7).
The strategy of averaging over all coordinate terms $i\in[m]$ has also studied
in variational inference, via coordinate ascent variational inference (CAVI)
[Bishop and Nasrabadi, 2006] which focuses on latent variable inference.

##  Appendix B Proofs and additional results

Assuming that if $f\in\mathcal{H}$ then $-f\in\mathcal{H}$ we can assume that
the supremum over the expectation is non-negative, and with Eq. 11,

$\displaystyle 0\leq{{\operatorname{NP-KSD}}}_{t}(P\|Q)$ $\displaystyle=$
$\displaystyle\sup_{f\in\mathcal{B}_{1}(\mathcal{H})}\mathbb{E}_{p}[\widehat{\mathcal{A}}_{t,N}^{B}f]$
(20) $\displaystyle=$
$\displaystyle\sup_{f}\\{\mathbb{E}_{p}\mathcal{A}_{t}f+\mathbb{E}_{p}[\widehat{\mathcal{A}}_{t,N}^{B}-\mathcal{A}_{t}]f\\}$
$\displaystyle=$
$\displaystyle\sup_{f}\\{\mathbb{E}_{p}\mathcal{A}_{t}f+\mathbb{E}_{p}[{\widehat{\mathcal{A}}}_{t,N}^{B}-\widehat{\mathcal{A}}_{t,N}]f+\mathbb{E}_{p}[f(\widehat{s}^{(i)}_{t,N}-\log
q_{t}^{\prime})]\\}.$

Here ${\widehat{\mathcal{A}}}_{t,N}$ is the Stein operator using the estimated
conditional score function $\hat{s}_{t,N}$ with the estimation based on $N$
synthetic observations. We now assess the contribution to 20 which stems from
estimating the score function. Note that here we only need to estimate a one-
dimensional score function and hence the pitfalls of score estimation in high
dimensions do not apply. We note however the contribution Zhou et al. [2020]
for a general framework.

Assume that we estimate the uni-variate conditional density ${q}_{t}^{(i)}$
based on $N$ samples. We assume that ${q}_{t}^{(i)}$ is differentiable, and we
denote its score function by

$s_{t}^{(i)}(x)=\frac{({q}_{t}^{(i)})^{\prime}(x)}{{q}_{t}^{(i)}(x)}.$

We next prove an extension of Proposition 3.1.

######  Proposition B.1.

Suppose that for $i\in[m]$, $\widehat{s}_{N}^{(i)}$ is a consistent estimator
of the uni-variate score function $s^{(i)}$. Let $\mathcal{T}^{(i)}$ be a
Stein operator for the uni-variate differentiable probability distribution
${Q}^{(i)}$ of the generalised density operator form Eq. (8). Let

$\displaystyle{\widehat{\mathcal{T}}}^{(i)}_{N}g(x)$ $\displaystyle=$
$\displaystyle g^{\prime}(x)+g(x)\widehat{s}_{N}^{(i)}$
$\displaystyle\widehat{\mathcal{A}}g(\boldsymbol{x})$ $\displaystyle=$
$\displaystyle{\widehat{\mathcal{T}}}^{(I)}_{N}g_{x^{(-I)}}(x^{(I)})\quad\quad\mbox{
and }$ $\displaystyle\widehat{\mathcal{A}}_{N}g(\boldsymbol{x})$
$\displaystyle=$
$\displaystyle{\frac{1}{m}\sum_{i\in[m]}{\widehat{\mathcal{T}}}^{(i)}_{N}g_{x^{(-i)}}(x^{(i)})}.$

Then ${\widehat{\mathcal{T}}}^{(i)}_{N}$ is a consistent estimator for
${{\mathcal{T}}}^{(i)}$, and $\widehat{\mathcal{A}}$ as well as
$\widehat{\mathcal{A}}_{N}$ are consistent estimators of $\mathcal{A}.$

###### Proof.

Take a fixed $x$. As $\widehat{s}_{N}^{(i)}$ is a consistent estimator of
$s^{(i)}$, it holds that for any $\epsilon>0$ and for any $x$ in the range of
$s^{(i)}$,

$\mathbb{P}(|\widehat{s}_{N}^{(i)}{(x)}-{s}^{(i)}(x)|>\epsilon)\rightarrow 0$

as $N\rightarrow\infty$. Here $\omega$ denotes the random element for the
estimation, which is implicit in $\widehat{q}_{N}$. On the set

$A_{\epsilon}=\left\\{|\widehat{s}_{N}^{(i)}{(x)}-{s}^{(i)}(x)|\leq\epsilon\right\\}$

we have that

$|{\widehat{\mathcal{T}}}^{{(i)}}_{N}g(x)-{{\mathcal{T}}}^{(i)}g(x)|\leq\epsilon
f(x).$

For every fixed $x$ this expression tends to 0 as $\epsilon\rightarrow 0$.
Hence consistency of $\widehat{T}_{N}$ follows. The last two assertions follow
immediately from Eq. (18) and Eq. (19). ∎

###  B.1 Asymptotic behaviour of NP-KSD

Here we assess the asymptotic behaviour of NP-KSD2. With $\boldsymbol{s}_{t}$
denoting the conditional score function,

${\operatorname{NP-
KSD}}^{2}_{t}(G\|p)=\mathbb{E}_{\boldsymbol{x},\boldsymbol{x}^{\prime}\sim
p}\left\langle\mathcal{A}_{Q_{t}}k(\boldsymbol{x},\cdot),\mathcal{A}_{Q_{t}}k(\boldsymbol{x}^{\prime},\cdot)\right\rangle_{\mathcal{H}}$

where
$\mathcal{A}_{Q_{t}}k(\boldsymbol{x},\cdot)=\mathcal{A}_{t}k(\boldsymbol{x},\cdot)$
can be written as

$\displaystyle\mathcal{A}_{t}k(\boldsymbol{x},\cdot)$ $\displaystyle=$
$\displaystyle\frac{1}{m}\sum_{{i\in[m]}}\left\\{\mathcal{A}_{\widehat{Q}^{(i)}_{t}}k(\boldsymbol{x},\cdot)+k(\boldsymbol{x},\cdot)(\widehat{s}^{(i)}_{t,{N}}-s^{(i)}_{t})\right\\}$
(21) $\displaystyle=$
$\displaystyle{\frac{1}{m}\sum_{i\in[m]}}\left\\{\frac{\partial}{\partial
x^{(i)}}k(\boldsymbol{x},\cdot)+k(\boldsymbol{x},\cdot){\widehat{s}}^{(i)}_{t,{N}}+k(\boldsymbol{x},\cdot)(\widehat{s}^{(i)}_{t,{N}}-s^{(i)}_{t})\right\\}.$

Recall that KSD2 is given in Eq. (5) by

$\mathrm{KSD}^{2}(q\|p)={\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle\mathcal{A}_{q}k(\boldsymbol{x},\cdot),\mathcal{A}_{q}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}],$

where ${\operatorname{KSD}}(q\|p)$ is a deterministic quantity which under
weak assumption vanishes when $p=q$. Moreover,

${{\operatorname{KSD}}}_{t}^{2}(q_{t}\|p)={\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle\mathcal{A}_{t}k(\boldsymbol{x},\cdot),\widehat{\mathcal{A}}_{t}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}.$

Disentangling this expression in general is carried out using Eq. (7).

######  Remark B.2.

For Gaussian kernels $k=k_{G}$ used in this paper, we can exploit its
factorisation;

$\displaystyle k_{G}(\boldsymbol{x},\tilde{\boldsymbol{x}})$ $\displaystyle=$
$\displaystyle\exp\left\\{-\frac{1}{2\sigma^{2}}\sum_{i=1}^{m}\left(x^{(i)}-{\tilde{x}}^{(i)}\right)^{2}\right\\}=\prod_{i=1}^{m}\exp\left\\{-\frac{1}{2\sigma^{2}}\left(x^{(i)}-{\tilde{x}}^{(i)}\right)^{2}\right\\}.$

In this situation, taking
$g_{{\boldsymbol{x}}}(\cdot)=k_{G}(\boldsymbol{x},\cdot)$, with $\cdot$
denoting an element in $\mathbb{R}^{m}$, gives

$g_{{{\boldsymbol{x}}};x^{(-i)}}({\cdot})=\exp\left\\{-\frac{1}{2\sigma^{2}}\sum_{j:j\neq
i}^{m}\left(x^{(j)}-(\cdot)^{(j)}\right)^{2}\right\\}\exp\left\\{-\frac{1}{2\sigma^{2}}\left(x^{(i)}-(\cdot)^{(i)}\right)^{2}\right\\}.$

For the operator ${\mathcal{T}}_{q}^{(i)}$ in Eq. (8) we have

$\displaystyle{\mathcal{T}}_{q}^{(i)}g_{x^{(-i)}}({\cdot})$ $\displaystyle=$
$\displaystyle{\exp\left\\{-\frac{1}{2\sigma^{2}}\sum_{j=1}^{m}\left(x^{(j)}-(\cdot)^{(j)}\right)^{2}\right\\}}\left(\frac{1}{\sigma^{2}}(x^{(i)}-(\cdot)^{(i)})+(\log
q_{t(x^{(-i)})})^{\prime}(x^{(i)})\right).$

Thus, the operator $\mathcal{A}_{t}$ decomposes as

$\displaystyle\frac{1}{m}\sum_{i=1}^{m}{\mathcal{T}}_{q}^{(i)}g_{x^{(-i)}}({\cdot})$
$\displaystyle=$
$\displaystyle\exp\left\\{-\frac{1}{2\sigma^{2}}\sum_{j=1}^{m}\left(x^{(j)}-(\cdot)^{(j)}\right)^{2}\right\\}$
$\displaystyle\frac{1}{m}\sum_{i=1}^{m}\left\\{\frac{1}{\sigma^{2}}\left(x^{(i)}-(\cdot)^{(i)}\right)+(\log
q_{t(x^{(-i)})})^{\prime}(x^{(i)})\right\\}.$

For consistency, in our setting the Stein operators are only applied to the
observations ${\boldsymbol{z}}_{1},\ldots,{\boldsymbol{z}}_{n}$ and hence for
our applications pointwise consistent estimation suffices, in the sense that
that for $i=1,\ldots,m$,
$\widehat{s}_{t,N}^{(i)}=\widehat{s}_{t,N}^{(i)}(x^{(i)})$ is a consistent
estimator of the uni-variate score function
$s_{t}^{(i)}=\\{\log(q(x^{(i)})|t(x^{{(-i)}})\\}^{\prime}$. Score matching
estimators often satisfy not only consistency but also asymptotic normality,
see for example Song et al. [2020]. Such an assumption is required for Theorem
3.2; recall that we use the notation
$\hat{\boldsymbol{s}}_{t,N}=(\hat{s}_{t,N}(x^{(i)}),i\in[m])$. To prove
Theorem 3.2 we re-state it for convenience.

######  Theorem B.3.

Assume that the score function estimator vector $\hat{s}_{t,N}$ is
asymptotically normal with mean $0$ and covariance matrix $N^{-1}\Sigma_{s}$.
Then ${{\operatorname{NP-KSD}}}_{t}^{2}(G\|p)$ converges in probability to
${{\operatorname{KSD}}}_{t}^{2}(q_{t}\|p)$ at rate at least
$\min(B^{-\frac{1}{2}},N^{-\frac{1}{2}})$.

Proof. We have from Eq. (5)

${{\operatorname{NP-
KSD}}}_{t}^{2}(G\|p)={\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle\widehat{\mathcal{A}}^{B}_{t,N}k(\boldsymbol{x},\cdot),\widehat{\mathcal{A}}^{B}_{t,N}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}.$

Expanding this expression, with $\mathcal{A}_{t}$ denoting the score Stein
operator in Eq. 7 for the conditional distribution $q_{t}$,

$\displaystyle{{\operatorname{NP-KSD}}}_{t}^{2}(G\|p)$ $\displaystyle=$
$\displaystyle{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle\mathcal{A}_{t}k(\boldsymbol{x},\cdot),\mathcal{A}_{t}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$
$\displaystyle+{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),\mathcal{A}_{t}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$
$\displaystyle+{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle\mathcal{A}_{t}k(\boldsymbol{x},\cdot),(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$
$\displaystyle+{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),\widehat{(}\mathcal{A}^{B}_{t,N}-\mathcal{A}_{t})k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$
$\displaystyle=$
$\displaystyle{\operatorname{KSD}}^{2}(q_{t}\|p)+2\,{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),\mathcal{A}_{t}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$
$\displaystyle+{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$

where we used the symmetry of the inner product in the last step. Now, for any
function $g$ for which the expression is defined,

$\displaystyle(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})g(\boldsymbol{x})$
$\displaystyle=$
$\displaystyle(\widehat{\mathcal{A}}^{B}_{t,N}-\widehat{\mathcal{A}}_{t,N})g(\boldsymbol{x})+(\widehat{\mathcal{A}}_{t,N}-\mathcal{A}_{t})g(\boldsymbol{x})$
(22)

recalling that ${\widehat{\mathcal{A}}}_{t,N}$ is the Stein operator using the
estimated conditional score function $\hat{s}_{t,N}$ with the estimation based
on $N$ synthetic observations.

To analyse Eq. 22 we first consider
$(\widehat{\mathcal{A}}_{t,N}-\mathcal{A}_{t})g(\boldsymbol{x})$;

$\displaystyle{\widehat{\mathcal{A}}_{t,N}}g(\boldsymbol{x})-\mathcal{A}_{t}g(\boldsymbol{x})=\frac{1}{m}\sum_{i=1}^{m}g(x^{(i)})(\widehat{s}_{t,N}^{(i)}(x^{(i)})-s_{t}^{(i)}(x^{(i)})).$
(23)

We note that it suffices to assume that for $i=1,\ldots,m$,
$\widehat{s}_{t,N}^{(i)}=\widehat{s}_{t,N}^{(i)}(x^{(i)})$ is a consistent
estimator of the uni-variate score function
$s_{t}^{(i)}=\\{\log(q(x^{(i)})|t(x^{(-i)})\\}^{\prime}$. Score matching
estimators often satisfy not only consistency but also asymptotic normality,
see for example Song et al. [2020]. If for $x^{(1)},\ldots,x^{(m)}$ the vector
$\hat{\boldsymbol{s}}_{t,N}=(\hat{s}_{t,N}(x^{(i)}),i\in[m])$ is
asymptotically normal with mean $0$ and covariance matrix $N^{-1}\Sigma_{s}$
then it follows from Eq. 23 that, asymptotically,
$\sqrt{N}(\widehat{\mathcal{A}}_{t,N}g(\boldsymbol{x})-\mathcal{A}_{q_{t}}g(\boldsymbol{x}))$
has a multivariate normal distribution and, in particular,
$(\widehat{\mathcal{A}}_{t,N}g(\boldsymbol{x})-\mathcal{A}_{q_{t}}g(\boldsymbol{x}))$
has fluctuations of the order $N^{-\frac{1}{2}}$.

For the term
$(\widehat{\mathcal{A}}^{B}_{t,N}-\widehat{\mathcal{A}}_{t,N})g(\boldsymbol{x})$
in Eq. 22 we have

$\displaystyle(\widehat{\mathcal{A}}^{B}_{t,N}-\widehat{\mathcal{A}}_{t,N})g(\boldsymbol{x})$
$\displaystyle=$
$\displaystyle\frac{1}{B}\sum_{b=1}^{B}{{\mathcal{T}}_{t,N}}^{(i_{b})}g(x)-{\mathcal{A}_{t,N}}g(\boldsymbol{x})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{m}\left\\{\frac{1}{B}\sum_{b=1}^{B}{{\mathcal{T}}_{t,N}}^{(i_{b})}g(\boldsymbol{x})\mathds{1}(i_{b}=i)-\frac{1}{m}{{\mathcal{T}}_{t,N}}^{(i)}g(\boldsymbol{x})\right\\}.$

Let $k_{i}=\sum_{b=1}^{B}{\mathds{1}}(i_{b}=i)$ the number of times that $i$
is re-sampled. Then $\mathbb{E}(k_{i})=\frac{B}{m}$ and we have

$\displaystyle(\widehat{\mathcal{A}}^{B}_{t,N}-\widehat{\mathcal{A}}_{t,N})g(\boldsymbol{x})$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{m}{\mathcal{T}}_{t,N}^{(i)}g(\boldsymbol{x})\left\\{\frac{1}{B}k_{i}-\frac{1}{m}\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{B}\sum_{i=1}^{m}{\mathcal{T}}_{t,N}^{(i)}g(\boldsymbol{x})\left\\{k_{i}-\mathbb{E}(k_{i})\right\\}.$

This term is known to be approximately mean zero normal with finite variance
$\Sigma(\hat{\boldsymbol{s}}_{t,N};g)$ (which depends on
$\hat{\boldsymbol{s}}_{t,N}$ and $g$) of order $B^{-1}$, see for example
Holmes and Reinert [2004], where an explicit bound on the distance to normal
is provided. This asymptotic normality holds for the operator given the
estimated conditional score function. As the bootstrap samples are drawn
independently of the score function estimator, without conditioning, the
unconditional distribution is a mixture of normal distributions. For an
estimator $\hat{\boldsymbol{s}}_{t,N}$ which is asymptotically normally
distributed, the variances $\Sigma(\hat{\boldsymbol{s}}_{t,N};g)$ will
converge to $\Sigma(\boldsymbol{s}_{t};g)$.

Thus, with Eq. 22,

$\displaystyle{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),\mathcal{A}_{t}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}]$
$\displaystyle=$
$\displaystyle{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}^{B}_{t,N}-\widehat{\mathcal{A}}_{t,N})k(\boldsymbol{x},\cdot),\mathcal{A}_{t}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}]+{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),\mathcal{A}_{t}k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}]$

with the first term is approximately a variance mixture of mean zero normals
tending to $0$ in probability at rate at least $B^{-\frac{1}{2}}$ as
$B\rightarrow\infty$, and the second term approximately a mean zero normal
variable tending to 0 in probability at rate at least $N^{-\frac{1}{2}}$ as
$N\rightarrow\infty$.

It remains to consider

${\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}.$

With Eq. 22 we have

$\displaystyle{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),(\widehat{\mathcal{A}}^{B}_{t,N}-\mathcal{A}_{t})k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$
(27) $\displaystyle=$
$\displaystyle{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}^{B}_{t,N}-\widehat{\mathcal{A}}_{t,N})k(\boldsymbol{x},\cdot),(\widehat{\mathcal{A}}^{B}_{t,N}-\widehat{\mathcal{A}}_{t,N})k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$
$\displaystyle+{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}^{B}_{t,N}-\widehat{\mathcal{A}}_{t,N})k(\boldsymbol{x},\cdot),(\widehat{\mathcal{A}}_{t,N}-\mathcal{A}_{t})k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$
$\displaystyle+{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),(\widehat{\mathcal{A}}^{B}_{t,N}-\widehat{\mathcal{A}}_{t,N})k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}$
$\displaystyle+{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),(\widehat{\mathcal{A}}_{t,N}-\mathcal{A}_{t})k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}.$

In Xu and Reinert [2021], Proposition 2, the following result is shown, using
the notation as above.

Let

$Y=\frac{1}{B^{2}}\sum_{s,t\in{[m]}}(k_{s}k_{t}-\mathbb{E}(k_{s}k_{t}))h_{x}(s,t).$

Assume that $h_{x}$ is bounded and that $Var(Y)$ is non-zero. Then if $Z$ is
mean zero normal with variance $Var(Y)$, there is an explicitly computable
constant $C>0$ such that for all three times continuously differentiable
functions $g$ with bounded derivatives up to order 3,

$|\mathbb{E}[g(Y)]-\mathbb{E}[g(Z)]\leq\frac{C}{B}.$

Moreover, using Equations (17)-(21) from Ouimet [2021], it is easy to see that
$Var(Y)$ is of the order $B^{-1}$. Hence, Term (27) tends to 0 in probability
at rate at least $B^{-1}$. Similarly, using that the bootstrap sampling is
independent of the score function estimation, Terms (27) and (27) tend to 0 in
probability at rate at least $(NB)^{-\frac{1}{2}}$. For Term (27), from Eq.
(23),

$\displaystyle{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{\mathcal{A}}_{t,N}-\mathcal{A}_{t})k(\boldsymbol{x},\cdot),(\widehat{\mathcal{A}}_{t,N}-\mathcal{A}_{t})k(\tilde{\boldsymbol{x}},\cdot)\rangle_{\mathcal{H}}]$
$\displaystyle=$
$\displaystyle\frac{1}{m^{2}}\sum_{i=1}^{m}\sum_{j=1}^{m}{\mathbb{E}}_{\boldsymbol{x},\tilde{\boldsymbol{x}}\sim
p}[\langle(\widehat{s}^{{(i)}}_{t,N}(x^{{(i)}})-s^{{(i)}}_{t}(x^{{(i)}}))k(x^{{(i)}},\cdot),(\widehat{s}^{{(j)}}_{t,N}(x^{{(j)}})-s^{{(j)}}_{t}(x^{{(j)}}))k(x_{j},\cdot)\rangle_{\mathcal{H}}].$

If $\hat{\boldsymbol{s}}_{t,N}$ is approximately normal as hypothesised, then
the inner product is approximately a covariance of order $N^{-1}$, and hence
the overall contribution from Term (27) is of order at most $N^{-1}$. This
finishes the proof. $\Box$

##  Appendix C Additional experimental details and results

###  C.1 Additional experiments

NP-KSD | B=5 | B=10 | B=20 | B=40 | (MMDAgg)
---|---|---|---|---|---
Runtime(s) | 4.65 | 6.56 | 8.43 | 10.44 | 5.02
Rejection Rate | 0.24 | 0.40 | 0.51 | 0.55 | 0.27
Table 3: Computational runtime for various re-sample size B: observed sample
size $n=100$; bootstrap size $b=200$; dimension $m=40$. The rejection rate is
used for power comparison; higher rejection rates indicate higher power.

#### Runtime

The computational runtime for each tests are shown in Table 3. MMDAgg runtime
is also shown as a comparison. From the result, we can see that NP-KSD runs
generally slower than permutation-based test, i.e. MMDAgg. This is mainly due
to the learning of conditional score functions and the Monte-Carlo based
bootstrap procedure. As the re-sample size $B$ increase, NP-KSD test requires
longer runtime. However, the rejection rate $B=20$ is approaching to that of
$B=40$ ( similiar observations also shown in 1(d)). MMDAgg generally has
faster computation due to permutation procedure of the test. However, it has
lower test power, which is only comparable to that of $B=5$ at which the
runtime advantage is not that obvious.

#### Training on synthetic distributions

We also train the deep generative models on the synthetic distributions
studied in Section 4 and perform model assessment on the trained models. We
consider the standard Gaussian and Mixture of two-component Gaussian problems.
We train a generative adversarial network with multi-layer perceptron
(GAN_MLP)888DCGAN studied in the main text is particularly useful for the
(high-dimensional) image dataset due to the convolutional neural network (CNN)
layers; DCGAN is not applicable for the problem in $\mathbb{R}^{2}$. and a
variational auto-encoder (VAE) [Kingma and Welling, 2013]. Noise-constrastive
score network NCSN is also trained to learn the score function followed by
annealed Langevin dynamics [Song and Ermon, 2019, 2020]. The trainings are
done via Adam optimiser [Kingma and Ba, 2014] with adaptive learning rate. The
rejection rates are reported in Table 4.

As shown in Table 4, uni-modal Gaussian distribution is easier to be learned
by the generative modelling procedures, as compared to the two-component
Mixture of Gaussian (MoG). As a result, the ${\operatorname{NP-KSD}}$_m
testing procedure shows higher rejection rate on trained MoG generative models
compared to that of Gaussian.999We note that ${\operatorname{NP-KSD}}$ and
${\operatorname{NP-KSD}}$_m with summary statistics taken to be the mean are
equivalent in the two-dimensional problem.. However, as these deep models are
not designed for training and sampling the simple low-dimensional
distribution, it is not surprising the procedure produce samples that not pass
the ${\operatorname{NP-KSD}}$ tests.

Inspired from the settings in Gorham and Mackey [2017], where KSD is used to
measure sample quality, we also apply ${\operatorname{NP-KSD}}$ tests on the
Stochastic Gradient Langevin Dynamics (SGLD) [Welling and Teh, 2011] sampling
procedure studied in Gorham and Mackey [2017]; in Gorham and Mackey [2017],
SGLD is referred to as Stochastic Gradient Fisher Scoring (SGFS). SGLD is
capable of sampling uni-modal distributions, while it can have problems
sampling multi-modal data. The rejection rates shown in Table 4 are slightly
higher than the test level for MoG, while the type-I error is well controlled
for the Gaussian case. Generated samples from SGLD are visualised Figure 2,
illustrating that the SGLD samples look plausible for the Gaussian model, but
less so for the MoG model.

| GAN_MLP | VAE | NCSN | SGLD | Real
---|---|---|---|---|---
Gaussian | 0.36 | 0.61 | 0.25 | 0.06 | 0.03
MoG | 0.78 | 0.92 | 0.45 | 0.12 | 0.04
Table 4: ${\operatorname{NP-KSD}}$_m rejection rate: observed sample size
$n=100$; bootstrap size is $200$. Here a low rejection rate indicates a good
type-1 error. NCSN performs best, among deep generative models, on both tasks
but still has a very high rejection rate. SGLD outperforms the deep generative
models.

Figure 2: Visualisation of samples generated from Stochastic Gradient Langevin
Dynamics (SGLD); top: Gaussian model, bottom: MoG.

###  C.2 Data visualisation

We show samples from the MNIST and CIFAR10 dataset, together with samples from
trained generative models, in Figure 3 and Figure 4, respectively.

(a) Real MNIST samples

(b) NCSN samples  (c) DCGAN samples  (d) GAN samples (e) VAE samples

Figure 3: MNIST samples

(a) Real samples (b) DCGAN samples  (c) CIFAR10.1 samples  (d) NCSN samples

Figure 4: CIFAR10 samples

##  Appendix D Equivalence to the multivariate score-Stein operator

Here we show that the operator in Eq. 7 is equivalent to the corresponding
multivariate score-Stein operator in Eq. 2, when they exist; the difference
being the factor $\frac{1}{m}$. Recall the set-up for score-Stein operators.
Let $q$ with smooth support $\Omega_{q}$ be differentiable. The score function
of $q$ is the function

$\displaystyle\mathbf{s}_{q}=\mathcal{T}_{\nabla,q}1=\nabla\log q=\frac{\nabla
q}{q}$

(with the convention that $\mathbf{s}_{q}\equiv 0$ outside of $\Omega_{q}$).
The score-Stein operator is the vector-valued operator

$\mathcal{A}_{q}=\nabla+\mathbf{s}_{q}\mathrm{I}_{{m}}$ (28)

acting on differentiable functions $g:\mathbb{R}^{m}\to\mathbb{R}$, with
$\mathrm{I}_{m}$ denoting the $m\times m$ identity matrix.

######  Proposition D.1.

When they exist, then the operators in Eq. 7 and in Eq. 2 differ only by a
factor $\frac{1}{m}$.

Proof. Writing $\partial_{i}$ for the derivative in direction $x^{(i)}$, the
score operator acting on differentiable functions
$g:\mathbb{R}^{m}\rightarrow\mathbb{R}$ can be written as

$\mathcal{A}_{p}g(x)=\sum_{i=1}^{m}\left\\{\partial_{i}g(x)+g(x)\partial_{i}(\log
q(x))\right\\}.$ (29)

Now, for $i\in[m],$

$q(x)=q(x^{(i)}|x^{(j)},j\neq i)q(x^{(j)},j\neq i)$

and hence

$\partial_{i}(\log q(x))=\partial_{i}\log q(x^{(i)}|x^{(j)},j\neq i).$

The assertion follows. $\Box$

######  Example D.2 (Bi-variate Gaussian).

Consider $x=(x^{(1)},x^{(2)})^{\top}\in\mathbb{R}^{2}$, i.e. $m=2$ and
$x\sim\mathcal{N}(\mu,\Sigma)$ where
$\mu=(\mu^{(1)},\mu^{(2)})^{\top}\in\mathbb{R}^{2}$,
$\Sigma=\begin{pmatrix}1&\sigma\\\ \sigma&1\end{pmatrix}$. With the
corresponding precision matrix
$\Sigma^{-1}=\frac{1}{1-\sigma^{2}}\begin{pmatrix}1&-\sigma\\\
-\sigma&1\end{pmatrix}$, it is easy to check
$Q^{(1)}(X^{(1)}|X^{(2)}=x^{(2)})\sim\mathcal{N}(\mu^{(1)}+\sigma(x^{(2)}-\mu^{(2)}),1-\sigma^{2})$.
For a bi-variate differentiable test function $g:\mathbb{R}^{2}\to\mathbb{R}$,
applying the Stein operator of the form in Equation 8,

$\displaystyle\mathcal{A}g(x^{(1)},x^{(2)})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left\\{{\mathcal{T}}^{(1)}g_{x^{(2)}}(x^{(1)})+{\mathcal{T}}^{(2)}g_{x^{(1)}}(x^{(2)})\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left\\{(g_{x^{(2)}})^{\prime}(x^{(1)})-\frac{x^{(1)}-\mu^{(1)}-\sigma(x^{(2)}-\mu^{(2)})}{1-\sigma^{2}}g_{x^{(2)}}(x^{(1)})\right.$
$\displaystyle\left.+(g_{x^{(1)}})^{\prime}(x^{(2)})-\frac{x^{(2)}-\mu^{(2)}-\sigma(x^{(1)}-\mu^{(1)})}{1-\sigma^{2}}g_{x^{(1)}}(x^{(2)})\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left\\{\partial_{1}g(x^{(1)},x^{(2)})-\frac{x^{(1)}-\mu^{(1)}-\sigma(x^{(2)}-\mu^{(2)})}{1-\sigma^{2}}g_{x^{(2)}}(x^{(1)})\right.$
$\displaystyle\left.+\partial_{2}g(x^{(1)},x^{(2)})-\frac{x^{(2)}-\mu^{(2)}-\sigma(x^{(1)}-\mu^{(1)})}{1-\sigma^{2}}g_{x^{(1)}}(x^{(2)})\right\\}$
$\displaystyle=$ $\displaystyle\frac{1}{2}\left\\{\nabla\times
g(x^{(1)},x^{(2)})-\Sigma^{-1}(x^{(1)}-\mu^{(1)},x^{(2)}-\mu^{(2)})^{T}g(x^{(1)},x^{(2)})\right\\}$

where $\partial_{i}$ denotes the derivative with respect to $x^{(i)}$. Thus,
we recover the score operator given in Eq. 2.

##  Appendix E Energy-based models and score matching

Energy-based models (EBMs) [LeCun et al., 2006] have been used in machine
learning contexts for modelling and learning deep generative models. An EBM is
essentially a Gibbs measure with energy function $E(x)$,

$q(x)=\frac{1}{Z}\exp\\{-E(x)\\},$ (30)

where $Z$ is the (generally) intractable normalisation constant (or partition
function). In particular, learning and training complicated EBM has been
studied in machine learning [Song and Kingma, 2021]. One of the most popular
and relatively stable training objective is the score-matching (SM) objective
given in Eq. (10) [Hyvärinen, 2005],

$J(p\|q)=\mathbb{E}_{p}\left[\left(\log p(x)^{\prime}-\log
q(x)^{\prime}\right)^{2}\right],$

which is particularly useful for the unnormalised models such as EBMs.

For an EBM, the SM objective only requires computing $\nabla E(x)$ and
$\nabla\cdot\nabla E(x)$ (or $\Delta E(x)$), which is independent of the
partition function $Z$. Moreover, by learning the SM objective, we can obtain
$\nabla\log q(x)$ directly, to construct the approximate Stein operator.

##  Appendix F More on kernel-based hypothesis tests

###  F.1 Maximum-mean-discrepancy tests

#### Maximum-mean-discrepancy

(MMD) has been introduced as a kernel-based method to tackle two-sample
problems [Gretton et al., 2007], utilising the rich representation of the
functions in a reproducing kernel Hilbert space (RKHS) via a kernel mean
embedding. Let $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ be the kernel
associated with RKHS $\mathcal{H}$. The kernel _mean embedding_ of a
distribution $p$ induced by $k$ is defined as

$\mu_{p}:=\mathbb{E}_{x\sim p}[k(x,\cdot)]\in\mathcal{H},$ (31)

whenever $\mu_{p}$ exist. The kernel mean embedding in Eq.31 can be estimated
empirically from independent and identically distributed (i.i.d.) samples.
Given $x_{1},\dots,x_{n}\sim p$:

$\widehat{\mu}_{p}:=\frac{1}{n}\sum_{i=1}^{n}k(x_{i},\cdot)$ (32)

replacing $p$ by its empirical counterpart
$\widehat{p}=\frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$ where
$\mathcal{\delta}_{x_{i}}$ denotes the Dirac measure at $x_{i}\in\mathcal{X}$.
For i.i.d. samples, the empirical mean embedding $\widehat{\mu}_{p}$ is a
$\sqrt{n}$-consistent estimator for $\mu_{p}$ in RKHS norm [Tolstikhin et al.,
2017], and with $n$ denoting the number of samples,
$\|\mu_{p}-\widehat{\mu}_{p}\|_{\mathcal{H}}=O_{p}(n^{-\frac{1}{2}})$. When
the sample size $n$ is small, the estimation error may not be negligible.

The MMD between two distributions $p$ and $q$ is defined as

$\displaystyle\operatorname{MMD}(p\|q;\mathcal{H})$
$\displaystyle=\sup_{\|{f}\|_{\mathcal{H}}\leq 1}\mathbb{E}_{x\sim
p}[f(x)]-\mathbb{E}_{{\tilde{x}}\sim q}[f({\tilde{x}})]$
$\displaystyle=\sup_{\|{f}\|_{\mathcal{H}}\leq 1}\left\langle
f,\mu_{p}-\mu_{q}\right\rangle_{\mathcal{H}}=\|\mu_{p}-\mu_{q}\|_{\mathcal{H}}.$
(33)

One desirable property for MMD is to be able to distinguish distributions in
the sense that $\operatorname{MMD}(p\|q;\mathcal{H})=0\Longleftrightarrow
p=q$101010Note that MMD is symmetric with respect to $p$, $q$, while KSD is
not symmetric with respect to $p$, $q$.. This property can be achieved via
characteristic kernels [Sriperumbudur et al., 2011]. It is often more
convenient to work with the squared version of MMD:

$\displaystyle\operatorname{MMD}^{2}(p\|q;\mathcal{H})$
$\displaystyle=\|\mu_{p}-\mu_{q}\|^{2}_{\mathcal{H}}=\left\langle\mu_{p},\mu_{p}\right\rangle+\left\langle\mu_{q},\mu_{q}\right\rangle-2\left\langle\mu_{p},\mu_{q}\right\rangle$
$\displaystyle=\mathbb{E}_{x,\tilde{x}\sim
p}k(x,\tilde{x})+\mathbb{E}_{y,\tilde{y}\sim
q}k(y,\tilde{y})-\mathbb{E}_{x\sim p,y\sim q}k(x,y).$ (34)

Given two sets of i.i.d. samples
$\mathbb{S}_{p}=\\{x_{1},\dots,x_{n}\\}\overset{i.i.d.}{\sim}p$ and
$\mathbb{S}_{q}=\\{y_{1},\dots,y_{l}\\}\overset{i.i.d.}{\sim}q$, an unbiased
estimator of Eq.34, based on the empirical estimate of kernel mean embedding
in Eq.32, is given by

$\displaystyle{{\operatorname{MMD}}^{2}_{u}}(\mathbb{S}_{p}\|\mathbb{S}_{q};\mathcal{H})=\frac{1}{\small
n(n-1)}\sum_{i\neq
i^{\prime}}k(x_{i},x_{i^{\prime}})+\frac{1}{l(l-1)}\sum_{j\neq
j^{\prime}}k(y_{j},y_{j^{\prime}})-\frac{2}{nl}\sum_{ij}k(x_{i},y_{j}).$ (35)

A two-sample test (or two-sample problem) aims to test the null hypothesis
$H_{0}:p=q$ against the alternative hypothesis $H_{1}:p\neq q$. It has been
shown that the asymptotic distribution of $n$-scaled statistic
$n\cdot{\operatorname{MMD}}_{u}^{2}(\mathbb{S}_{p}\|\mathbb{S}_{q};\mathcal{H})$
under the null ($p=q$) is that of an infinite weighted sum of
$\chi^{2}$-distribution [Gretton et al., 2012a, Theorem 12], while under the
alternative ($p\neq q$), the $\sqrt{n}$-scaled statistic
$\sqrt{n}\cdot{\operatorname{MMD}}_{u}^{2}(\mathbb{S}_{p}\|\mathbb{S}_{q};\mathcal{H})$
is asymptotically normally distributed with the mean centered at
${\operatorname{MMD}}(p\|q;\mathcal{H})>0$. Thus,
$n\cdot{\operatorname{MMD}}_{u}^{2}(\mathbb{S}_{p}\|\mathbb{S}_{q};\mathcal{H})$
is taken as a test statistic to be compared against the _rejection threshold_.
If the test statistic exceeds the rejection threshold, the empirical
estimation of the MMD statistic is thought to exhibit significant departure
from the null hypothesis so that $H_{0}$ is rejected.

As the null distribution is given by an infinite weighted sum of $\chi^{2}$
random variables which does not have a closed form expression, the null
distribution can be simulated via a permutation procedure [Gretton et al.,
2008]: Combine and order two sets of samples as $z_{i}=x_{i},i\in[n]$ and
$z_{j}=y_{j-n},n+1\leq j\leq n+l$. Let $\mu:[n+l]\to[n+l]$ be a permutation,
and write $z^{\mu}=\\{z_{\mu_{(1)}},\dots z_{\mu_{(n+l)}}\\}$. Then $z^{\mu}$
is re-split into $\mathbb{S}_{p}^{\mu}=\\{z_{i}\\}_{1\leq i\leq n}$ and
$\mathbb{S}_{q}^{\mu}=\\{z_{j}\\}_{n+1\leq j\leq n+l}$. The permuted MMD is
computed via Eq. 35 as

${\operatorname{MMD}}^{2}_{u}(z^{\mu})={{\operatorname{MMD}}^{2}_{u}}(\mathbb{S}^{\mu}_{p}\|\mathbb{S}^{\mu}_{q};\mathcal{H}).$
(36)

For $\mu_{1},\dots,\mu_{B}$, we obtain
${\operatorname{MMD}}^{2}_{u}(z^{\mu_{1}}),\dots,{\operatorname{MMD}}^{2}_{u}(z^{\mu_{B}})$
and use these values to compute the empirical quantile of the test statistics
${{\operatorname{MMD}}^{2}_{u}}(\mathbb{S}_{p}\|\mathbb{S}_{q};\mathcal{H})$.

To test whether the implicit generative model can generate samples following
the same distribution as the observed sample, it is natural to consider the
two-sample problem described above, which tests whether two sets of samples
are from the same distribution. In the model assessment context, one set of
samples (of size $N$) are generated from the implicit model, while the other
set of samples (of size $n$) are observed.

The MMD test often assumes that the sample sizes $n$ and $l$ are equal; the
asymptotic results including consistency are valid under the regime that
$n,l\rightarrow\infty$ [Gretton et al., 2009, 2012b, Jitkrittum et al., 2016];
also the relative model comparisons in Jitkrittum et al. [2018] only
considered the cases $n=l$. In our setting, the sample size $l$ is usually
denoted by $N$. For our model assessment problem setting, when $n$ is fixed
and $N\to\infty$ is allowed to be asymptotically large, i.e. $n\ll N$, we find
that the type-I error may not controlled. Hence it is not always the case that
MMD is able to pick up the distributional difference between two sets of
samples under the null hypothesis. A simple experiment in Table 5 shows an
example in which the type-I error is not controlled when $N$ is increasing.
Hence, MMD is  not used as comparison for NP-KSD.

The high rejection rate of the MMD statistic, i.e. the high type-I error as
$N$ increase, shown in Table 5 can be heuristically explained as follows. Let
$\\{x_{1},\dots,x_{n}\\},\\{\tilde{x}_{1},\dots,\tilde{x}_{N}\\}\overset{i.i.d.}{\sim}p$
where the two sets of samples are generated from the same distribution. Let
$\widehat{\mu}_{p,n}=\frac{1}{n}\sum_{i\in[n]}k(x_{i},\cdot)$, and
$\widehat{\mu}_{p,N}=\frac{1}{N}\sum_{j\in[N]}k(\tilde{x}_{j},\cdot)$. The
empirical MMD between $\widehat{\mu}_{p,n}$ and $\widehat{\mu}_{p,N}$ can be
seen as

$\|\widehat{\mu}_{p,n}-\widehat{\mu}_{p,N}\|^{2}_{\mathcal{H}}=\|(\widehat{\mu}_{p,n}-\mu_{p})-(\widehat{\mu}_{p,N}-\mu_{p})\|^{2}_{\mathcal{H}},$
(37)

where MMD aims to detect the asymptotic equality of
$(\widehat{\mu}_{p,n}-\mu_{p})$ and $(\widehat{\mu}_{p,N}-\mu_{p})$. When $n$
is small and fixed, and $n\ll N$, the difference is non-trivial and a rich-
enough kernel is able detect this difference, leading to MMD rejecting the
null hypothesis although it is true.

Sample size N | 20 | 50 | 100 | 200 | 500 | 1000
---|---|---|---|---|---|---
MMD | 0.08 | 0.06 | 0.36 | 0.9 | 1.00 | 1.00
MMDAgg | 0.06 | 0.07 | 0.02 | 0.03 | 0.02 | 0.05
KSD | 0.07 | 0.04 | 0.04 | 0.02 | 0.08 | 0.06
Table 5: Type-I error with increasing sample size $N$. $H_{0}$ is the standard
Gaussian with $m=3$; $n=50$; $\alpha=0.05$; $n_{sim}=500$, $100$ trials for
rejection rate. Bold values show the uncontrolled type-I error.

#### MMDAgg, a non-asymptotic MMD-based test

Recently, Schrab et al. [2021] proposed an aggregated MMD test that can
incorporate the setting $n\neq N$ as long as there exists a constant $C>0$
such that $n\leq N\leq cn$. Under this condition, MMDAgg is a consistent non-
asymptotic test with controlled type-I error, see Table 5, which we use as
competitor method in the main text. In Table 5 KSD is included as desired
method when the underlying null distribution is known.

The MMDAgg test statistic is computed by aggregating a set of MMD statistic
based on different choices of kernel, e.g. Gaussian bandwidth. MMDAgg takes
into account a number of choices of bandwidth $\lambda\in\Lambda$ where
$\Lambda$ is a finite set. Let $\widehat{M}_{\lambda}$ denote the empirical
MMD using a kernel with bandwidth $\lambda$. Each $\lambda$ can be weighted,
via $w_{\lambda}$, where $\sum_{\lambda\in\Lambda}w_{\lambda}=1$. In Schrab et
al. [2021] (as well as our implemented experiments), uniform weights are
chosen; $w_{\lambda}\equiv w=\frac{1}{|\Lambda|}$. Denote by $B_{1}$ the
number of samples used to simulate the null distribution for quantile
estimation 111111This is the same as notion $B$ in the main text as well as in
Eq. 36. and denote by $B_{2}$ the number of simulated samples used to estimate
the empirical rejection probability. Define
${\widehat{q}}_{\lambda,1-\alpha}^{B_{1}}(z^{B_{1}})$ as the conditional
empirical $(1-\alpha)$-quantile when MMD uses a kernel with bandwidth
$\lambda$, estimated from the permutation procedure with $B_{1}$ permutations
using Eq.36. Then for a fixed test level $\alpha$, u is estimated via the bi-
section method such that

$\mathbb{P}\left(\max_{\lambda\in\Lambda}(\widehat{M}_{\lambda}-{\widehat{q}}_{\lambda,uw}^{B_{1}}(z^{B_{1}}))>0\right){\leq\alpha}.$
(38)

We reject $H_{0}$ if for any $\lambda\in\Lambda$ and with the estimated
$\widehat{u}$, we have that $\widehat{M}_{\lambda}$ exceeds the rejection
probability in Eq. 38; otherwise we do not reject $H_{0}$. In this way, MMDAgg
does not only achieve the desired non-asymptotic type-I error but is able to
explore a wide range of kernels in order to produce stronger test power.

###  F.2 Wild-bootstrap on KSD testing procedures

The wild bootstrap procedure [Chwialkowski et al., 2014] simulates the null
distribution via so-called wild-bootstrap samples. For KSD, Chwialkowski et
al. [2014] has shown weak asymptotic convergence to the null distribution with
deterministic and bounded kernels. For the NP-KSD test statistic, wild-
bootstrap samples do not necessarily converge to the null distribution, due to
the estimation difference $(\widehat{s}^{t}-s)$, creating a random Stein
kernel for NP-KSD. Perhaps therefore unsurprisingly, the wild bootstrap
procedure NP-KSD does not control the type-I error correctly. Instead, we
consider a Monte Carlo procedure to simulate the null distribution of NP-KSD.
While Monte Carlo estimation is more computationally intensive than wild-
bootstrap, it is an accurate method by design.

Figure 5 illustrates this point. Figure 5(a) shows samples from a Gaussian
distribution. The true density is plotted in red. Two score matching density
estimates, SM1 and SM2, are calculated; SM1 presents a good fit whereas SM2 is
a less accurate estimate. For KSD, which is applicable when the underlying
null distribution is known, Figure 5(b) shows that the Monte Carlo
distribution and the wild-bootstrap distribution are close and reach the same
conclusion for the KSD test statistic. Using the well-fitting SM1 score
density estimate, Figure 5(c) gives the Monte Carlo distribution and the wild-
bootstrap distribution. The wild-bootstrap distribution is close to the wild-
bootstrap distribution for the KSD. In contrast, it differs considerably from
the Monte Carlo distribution and would reject the null hypothesis although it
is true. Figure 5(d) shows a similar result for the not so well fitting
estimate SM2. The wild-bootstrap distribution is now more spread out but the
observed test statistic is still in the tail of this distribution, whereas it
is close to the center of the Monte Carlo distribution. In the synthetic
experiments for MoG in the main text, the model misspecification NP-KSD_G
falls under this setting. These plots illustrate that using wild-bootstrap
samples in this scenario could lead to erroneous conclusions. Hence we use
Monte Carlo samples.

(a) Samples and fitted densities

(b) Simulated null distributions from KSD

(c) Simulated null distributions from fitted SM1

(d) Simulated null distributions from fit SM2

Figure 5: Visualisation for NP-KSD and KSD testing procedures. For KSD, the
wild-bootstrap distribution roughly agrees with the Monte Carlo distribution,
whereas for NP-KSD, the wild-bootstrap distribution deviates strongly from the
Monte Carlo distribution, indicating a danger of reaching an erroneous
conclusion when using wild-bootstrap samples in this scenario.