<figure><img src="image_1.png" /></figure>

HAL
open science

# KSD Aggregated Goodness-of-fit Test

Antonin Schrab, Benjamin Guedj, Arthur Gretton

► To cite this version:

Antonin Schrab, Benjamin Guedj, Arthur Gretton. KSD Aggregated Goodness-of-fit Test. 2022. [https://hal-03554423](https://hal-03554423)

HAL Id: hal-03554423

https://inria.hal.science/hal-03554423

Preprint submitted on 3 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# KSD Aggregated Goodness-of-fit Test

Antonin Schrab, Benjamin Guedj, Arthur Gretton

► To cite this version:

Antonin Schrab, Benjamin Guedj, Arthur Gretton. KSD Aggregated Goodness-of-fit Test. Conference on Neural Information Processing Systems, Nov 2022, New Orleans, United States. hal-03554423v2

HAL Id: hal-03554423

https://inria.hal.science/hal-03554423v2

Submitted on 22 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
KSD Aggregated Goodness-of-fit Test

Antonin Schrab<sup>1,2,3</sup> Benjamin Guedj<sup>1,3</sup> Arthur Gretton<sup>2</sup>

Abstract

We investigate properties of goodness-of-fit tests based on the Kernel Stein Discrepancy (KSD). We introduce a strategy to construct a test, called KSDAGG, which aggregates multiple tests with different kernels. KSDAGG avoids splitting the data to perform kernel selection (which leads to a loss in test power), and rather maximises the test power over a collection of kernels. We provide theoretical guarantees on the power of KSDAGG: we show it achieves the smallest uniform separation rate of the collection, up to a logarithmic term. KSDAGG can be computed exactly in practice as it relies either on a parametric bootstrap or on a wild bootstrap to estimate the quantiles and the level corrections. In particular, for the crucial choice of bandwidth of a fixed kernel, it avoids resorting to arbitrary heuristics (such as median or standard deviation) or to data splitting. We find on both synthetic and real-world data that KSDAGG outperforms other state-of-the-art adaptive KSD-based goodness-of-fit testing procedures.

# 1. Introduction

Kernel selection remains a fundamental question in kernel-
based nonparametric hypothesis testing, as it significantly
impacts the test power. Kernel selection has attracted a sig-
nificant interest in the literature, and a number of methods
have been proposed in different settings, such as in the two-
sample, independence and goodness-of-fit testing frame-
works. Those methods include using heuristics (Gretton
et al., 2012a), relying on data splitting (Gretton et al., 2012b;
Sutherland et al., 2017), learning deep kernels (Grathwohl
et al., 2020; Liu et al., 2020), working in the post-selection
inference framework (Yamada et al., 2019; Lim et al., 2019;
2020; Freidling et al., 2021), to name but a few.

In this work, we focus on aggregated tests, which have been
investigated for the two-sample problem by Fromont et al.
(2013), Kim et al. (2020) and Schrab et al. (2021) using
the Maximum Mean Discrepancy (MMD, Gretton et al.,
2012a) and for the independence problem by Albert et al.
(2019) and Kim et al. (2020) using the Hilbert Schmidt
Independence Criterion (HSIC, Gretton et al., 2005). We
extend the use of aggregated tests to the goodness-of-fit
setting, where we are given a model and some samples,
and we are interested in deciding whether the samples have
been drawn from the model. We employ the Kernel Stein
Discrepancy (KSD, Chwialkowski et al., 2016; Liu et al.,
2016) as our test statistic, which is an ideal measure of
distance for this setting: it admits an estimator which can
be computed without requiring samples from the model,
and does not require the model to be normalised. To the
best of our knowledge, ours represents the first aggregation
procedure for the KSD test in the literature.

**Related work.** Fromont et al. (2012; 2013) introduced aggregated tests for the two-sample problem with the equal sample sizes following a Poisson process, using as a test statistic an unscaled version of the MMD. They also provided theoretical results in terms of uniform separation rates using a wild bootstrap. Albert et al. (2019) then proposed an aggregated test for the independence problem using the HSIC, with guarantees using a theoretical quantile, but relying on permutations to obtain the test threshold. Kim et al. (2020) then extended those results to also hold for the estimated quantile, and generalised the two-sample results to hold for the MMD estimator with different sample sizes using a wild bootstrap. All those aforementioned results were proved for the Gaussian kernel only. Schrab et al. (2021) generalised the two-sample results to hold for a wide range of kernels using either a wild bootstrap or a permutation-based procedure, and provided optimality results which hold with fewer restrictions on the class of functions. Our work builds and extends on the above results: we consider the goodness-of-fit framework, where we have samples from only one of the two densities. The main challenges arise from working with the Stein kernel that defines the KSD test statistic: for example, we lose the transition invariant property of the kernel which is crucial to work in the Fourier domain. Finally, Key et al. (2021) addresses the complementary task of Stein test design for a family of models, rather than for a single model.

<footnote n="1">¹Centre for Artificial Intelligence, Department of Com-
puter Science, University College London
²Gatsby Com-
putational Neuroscience Unit, University College London
³Inria, Lille - Nord Europe research centre and Inria Lon-
don Programme. Correspondence to: Antonin Schrab
<a.schrab@ucl.ac.uk>, Benjamin Guedj <b.guedj@ucl.ac.uk>,
Arthur Gretton <arthur.gretton@gmail.com>.</footnote>
# KSD Aggregated Goodness-of-fit Test

Antonin Schrab

Centre for Artificial Intelligence
Gatsby Computational Neuroscience Unit
University College London & Inria London
a.schrab@ucl.ac.uk

Benjamin Guedj

Centre for Artificial Intelligence
University College London & Inria London
b.guedj@ucl.ac.uk

Arthur Gretton

Gatsby Computational Neuroscience Unit
University College London
arthur.gretton@gmail.com

## Abstract

We investigate properties of goodness-of-fit tests based on the Kernel Stein Discrepancy (KSD). We introduce a strategy to construct a test, called KSDAGG, which aggregates multiple tests with different kernels. KSDAGG avoids splitting the data to perform kernel selection (which leads to a loss in test power), and rather maximises the test power over a collection of kernels. We provide non-asymptotic guarantees on the power of KSDAGG: we show it achieves the smallest uniform separation rate of the collection, up to a logarithmic term. For compactly supported densities with bounded model score function, we derive the rate for KSDAGG over restricted Sobolev balls; this rate corresponds to the minimax optimal rate over unrestricted Sobolev balls, up to an iterated logarithmic term. KSDAGG can be computed exactly in practice as it relies either on a parametric bootstrap or on a wild bootstrap to estimate the quantiles and the level corrections. In particular, for the crucial choice of bandwidth of a fixed kernel, it avoids resorting to arbitrary heuristics (such as median or standard deviation) or to data splitting. We find on both synthetic and real-world data that KSDAGG outperforms other state-of-the-art quadratic-time adaptive KSD-based goodness-of-fit testing procedures.

## 1 Introduction

Kernel selection remains a fundamental problem in kernel-based nonparametric hypothesis testing, as it significantly impacts the test power. Kernel selection has attracted a significant interest in the literature, and a number of methods have been proposed in the two-sample, independence and goodness-of-fit testing frameworks, such as using heuristics (Gretton et al., 2012a), relying on data splitting (Gretton et al., 2012b; Sutherland et al., 2017; Kübler et al., 2022), learning deep kernels (Grathwohl et al., 2020; Liu et al., 2020), working in the post-selection inference framework (Yamada et al., 2019; Lim et al., 2019, 2020; Kübler et al., 2020; Freidling et al., 2021), to name but a few.

In this work, we focus on aggregated tests, which have been investigated for the two-sample problem by Fromont et al. (2013), Kim et al. (2022) and Schrab et al. (2021) using the Maximum Mean Discrepancy (MMD; Gretton et al., 2012a), and for the independence problem by Albert et al. (2022) and Kim et al. (2022) using the Hilbert Schmidt Independence Criterion (HSIC; Gretton et al., 2005). We extend the use of aggregated tests to the goodness-of-fit setting, where we are given a model and some samples, and test whether the samples have been drawn from the model. We employ the Kernel Stein Discrepancy (KSD; Chwialkowski et al., 2016; Liu et al., 2016) as our test statistic, which is an ideal measure of distance for this setting: it admits an estimator which can be computed without
**Contributions.** We introduce our adaptive test KSDAGG which aggregates multiple tests with different kernels. This is a general setting which has many applications, such as kernel bandwidth selection. Our aggregated test allows for two numerical methods for estimating the test thresholds: the parametric bootstrap and the wild bootstrap. We conduct a theoretical analysis: we provide a lower bound on the uniform separation rate (Baraud, 2002) of KSDAGG, a condition which guarantees test power. We discuss the implementation of KSDAGG and its use in practice through experiments on synthetic and real-world data, on which we observe that KSDAGG obtains higher power than other KSD-based adaptive state-of-the-art tests. The code is publicly available at https://github.com/antoninschrab/ksdagg-paper/.

**Outline.** Section 2 presents our framework and our notation. Section 3 introduces our algorithm KSDAGG (in Algorithm 1) and contains our main theoretical results, and Section 4 presents numerical experiments to support KSDAGG. We close the paper with avenues for future research in Section 5.

## 2. Notation

We consider the goodness-of-fit problem where given access to a known probability density $p$ (model) and to some i.i.d. $d$-dimensional samples $\mathcal{X}_N := (X_i)_{i=1}^N$ drawn from an unknown density $q$, we want to decide whether $p \neq q$ holds. This can be expressed as a statistical hypothesis testing problem with null hypothesis $H_0: p = q$ and alternative $H_a: p \neq q$.

As a measure of distance between $p$ and $q$, we use the *Kernel Stein Discrepancy (KSD)* introduced by Chwialkowski et al. (2016) and Liu et al. (2016). For a kernel $k$, the KSD is defined as the Maximum Mean Discrepancy (MMD, Gretton et al., 2012a) between $p$ and $q$ using the Stein kernel associated to $k$, that is

$$ 
\begin{aligned}
\text{KSD}_{p,k}^2(q) &:= \text{MMD}_{h_{p,k}}^2(p, q) \\
&:= \mathbb{E}_{q,q}[h_{p,k}(X, Y)] - 2\mathbb{E}_{q,q}[h_{p,k}(X, Y)] \\
&\quad + \mathbb{E}_{p,p}[h_{p,k}(X, Y)] \\
&= \mathbb{E}_{q,q}[h_{p,k}(X, Y)]
\end{aligned}
$$

where the *Stein kernel* $h_{p,k}: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ is defined as

$$ 
\begin{aligned}
h_{p,k}(x, y) &:= (\nabla \log p(x)^{\top} \nabla \log p(y)) k(x, y) \\
&\quad + \nabla \log p(y)^{\top} \nabla_1 k(x, y) \\
&\quad + \nabla \log p(x)^{\top} \nabla_2 k(x, y) \\
&\quad + \sum_{i=1}^{d} \frac{\partial}{\partial x_i \partial y_i} k(x, y)
\end{aligned}
$$

where $\nabla_i$ denotes the gradient with respect to the $i$-th component. The Stein kernel satisfies the *Stein identity*

$$ \mathbb{E}_p[h_{p,k}(X, \cdot)] = 0. $$

A quadratic-time KSD estimator can be computed as the U-statistic (Hoeffding, 1992)

$$ \widehat{\text{KSD}}_{p,k}^2(\mathcal{X}_N) := \frac{1}{N(N-1)} \sum_{1 \le i \ne j \le N} h_{p,k}(X_i, X_j). \quad (1) $$

In this work, the model density $p$ is always known, we do not always explicitly write the dependence of $p$ for all variables as we do for $\text{KSD}_{p,k}^2$ and $h_{p,k}$.

We assume that the kernel $k$ is such that

$$ \text{KSD}_{p,k}^2(q) = \mathbb{E}_{q,q}[h_{p,k}(X, Y)] < \infty \quad (2) $$

and

$$ C_k := \mathbb{E}_{q,q}[h_{p,k}(X, Y)^2] < \infty. \quad (3) $$

We now address the requirements for consistency of the Stein test (Chwialkowski et al., 2016, Theorem 2.2): we assume that the kernel $k$ is $C_0$-universal (Carmeli et al., 2010, Definition 4.1) and that $\mathbb{E}_q \left\| \nabla \left( \log \frac{p(X)}{q(X)} \right) \right\|_2^2 < \infty$.

We use the notations $\mathbb{P}_p$ and $\mathbb{P}_q$ to denote the probability under the model $p$ and under $q$, respectively. Given a kernel $\kappa: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ and a function $f: \mathbb{R}^d \to \mathbb{R}$, we consider the integral transform

$$ (\kappa \diamond f)(y) := \int_{\mathbb{R}^d} \kappa(x, y) f(x) \, dx $$

following the notation of Fromont et al. (2012). When the kernel $\kappa$ is translation invariant, the operation $\diamond$ corresponds to a convolution, however this is not true of the Stein kernel.

## 3. Construction of tests and bounds

We now introduce the single and aggregated KSD tests. We show that these control the probability of type I error as desired, and provide conditions for the control of type II error probability.

### 3.1. Single test

We first construct a KSD test for a fixed kernel $k$ as proposed by Chwialkowski et al. (2016) and Liu et al. (2016). To estimate the test threshold, we can either use a wild bootstrap (Shao, 2010; Leucht & Neumann, 2013; Fromont et al., 2012) or a parametric bootstrap (Stute et al., 1993). Both methods work by simulating sampling values $(\bar{K}_k^1, \dots, \bar{K}_k^{B_1})$ from the (asymptotic) distribution of
requiring samples from the model, and does not require the model to be normalised. To the best of our knowledge, ours represents the first aggregation procedure for the KSD test in the literature.

**Related work.** Fromont et al. (2012, 2013) introduced non-asymptotic aggregated tests for the two-sample problem with sample sizes following a Poisson process, using an unscaled version of the MMD. Using a wild bootstrap, they derived uniform separation rates (Ingster, 1987, 1989, 1993a,b; Baraud, 2002). Albert et al. (2022) then proposed an independence aggregated test using the HSIC, with guarantees using a theoretical quantile, but relying on permutations to obtain the test threshold in practice. Kim et al. (2022) then extended those theoretical results to also hold for the quantile estimated using permutations, however, they did not obtain the desired level dependency in their HSIC bound. All those aforementioned results were proved for the Gaussian kernel only. Schrab et al. (2021) generalised the two-sample results to hold for the usual MMD estimator and for a wide range of kernels using either a wild bootstrap or permutations, and provided optimality results which hold with fewer restrictions. Our work builds and extends on the above non-asymptotic results: we consider the goodness-of-fit framework, where we have samples from only one of the two densities. The main challenges arise from working with the Stein kernel which defines the KSD test statistic: for example, we lose the transition-invariant property of the kernel which is crucial to work in the Fourier domain. Balasubramanian et al. (2021) considered adaptive MMD-based goodness-of-fit tests and obtained their uniform separation rates over Sobolev balls in the asymptotic regime. More generally, Li and Yuan (2019) studied asymptotic adaptive kernel-based tests for the three testing frameworks. Tolstikhin et al. (2016) derived minimax rates for MMD estimators using radial universal kernel (Sriperumbudur et al., 2011). Schrab et al. (2022) extends this work, together with those of Albert et al. (2022) and Schrab et al. (2021), to construct efficient aggregated tests for the three testing frameworks using incomplete U-statistics. They quantify the cost in the minimax rate over Sobolev balls incurred for computational efficiency, and prove minimax optimality of the quadratic-time HSIC permuted aggregated test by improving the bound of Kim et al. (2022). See Appendix C for details.

**Contributions.** We propose a solution to the fundamental kernel selection problem for the widely-used KSD goodness-of-fit tests: we construct an adaptive test KSDAGG which aggregates multiple tests with different kernels. Our contribution is in showing, both theoretically and experimentally, that the aggregation procedure works in this novel setting in which it has never been considered before. We work in the kernel selection framework; this general setting has many applications including the one of kernel bandwidth selection. Our aggregated test allows for two numerical methods for estimating the test thresholds: the wild bootstrap and the parametric bootstrap (a procedure unique to the goodness-of-fit framework). We conduct a theoretical analysis: we derive a general condition which guarantees test power for KSDAGG in terms of its uniform separation rate, with extra assumptions including regularity over restricted Sobolev balls we prove that KSDAGG attains the minimax rate over (unrestricted) Sobolev balls. We discuss the implementation of KSDAGG and experimentally validate our proposed approach on benchmark problems, not only on datasets classically used in the literature but also on original data obtained using state-of-the-art generative models (i.e. Normalizing Flows). We observe, both on synthetic and real-world data, that KSDAGG obtains higher power than other KSD-based adaptive state-of-the-art tests. Contributing to the real-world applications of these goodness-of-fit tests, we provide publicly available code to allow practitioners to employ our method: https://github.com/antoninschrab/ksdagg-paper.

## 2 Notation

We consider the goodness-of-fit problem where given access to a known model probability density $p$ (which can be unnormalised since we actually only need access to its score function $\nabla \log p(\cdot)$) and to some i.i.d. $d$-dimensional samples $\mathbb{X}_N := (X_i)_{i=1}^N$ drawn from an unknown data density $q$, we want to decide whether $p \neq q$ holds. This can be expressed as a statistical hypothesis testing problem with null hypothesis $H_0: p = q$ and alternative $H_a: p \neq q$.

As a measure of distance between $p$ and $q$, we use the *Kernel Stein Discrepancy* (KSD) introduced by Chwialkowski et al. (2016) and Liu et al. (2016). For a kernel $k$, the KSD is the Maximum Mean Discrepancy (MMD; Gretton et al., 2012a) between $p$ and $q$ using the Stein kernel associated to $k$

$$ 
\begin{aligned}
\text{KSD}_{p,k}^2(q) &:= \text{MMD}_{h_{p,k}}^2(p, q) := \mathbb{E}_{q,q}[h_{p,k}(X, Y)] - 2\mathbb{E}_{p,q}[h_{p,k}(X, Y)] + \mathbb{E}_{p,p}[h_{p,k}(X, Y)] \\
&= \mathbb{E}_{q,q}[h_{p,k}(X, Y)]
\end{aligned}
$$
$\widehat{\text{KSD}}_{p,k}^2$ under the null hypothesis and estimating the $(1-\alpha)$-quantile using a Monte Carlo approximation<fnref n="1" />

$$ \hat{q}_{1-\alpha}^k := \inf \left\{ u \in \mathbb{R} : 1 - \alpha \le \frac{1}{B_1} \sum_{b=1}^{B_1} 1(\bar{K}_k^b \le u) \right\} \\ = \bar{K}_k \cdot^{[B_1(1-\alpha)]} $$

where $\bar{K}_k \cdot^{[B_1(1-\alpha)]} \le \dots \le \bar{K}_k \cdot^{B_1}$ are the sorted elements $(\bar{K}_k^1, \dots, \bar{K}_k^{B_1})$. The single test is then defined as

$$ \Delta_{\alpha}^{k}(\mathbb{X}_{N}) := 1(\widehat{\text{KSD}}_{p,k}^{2}(\mathbb{X}_{N}) > \hat{q}_{1-\alpha}^{k}). $$

For the parametric bootstrap, we directly draw new samples $(X'_i)_{i=1}^{N'}$ from the model $p$ and compute the KSD

$$ \bar{K}_k := \frac{1}{N'(N' - 1)} \sum_{1 \le i \ne j \le N'} h_{p,k}(X'_i, X'_j). \quad (4) $$

For the wild bootstrap, we first generate $n$ i.i.d. Rademacher random variables $\epsilon_1, \dots, \epsilon_n$ taking values in $\{-1, 1\}^n$, and then compute

$$ \bar{K}_k := \frac{1}{N(N-1)} \sum_{1 \le i \ne j \le N} \epsilon_i \epsilon_j h_{p,k}(X_i, X_j). \quad (5) $$

Both these processes are then repeated $B_1$ times.

Since it uses samples from the model $p$, the parametric bootstrap results in a test with non-asymptotic level $\alpha$. This comes at the cost of being computationally more expensive and assuming that we are able to sample from $p$ (which may be out of reach in some settings). Conversely, the wild bootstrap has the advantage of not requiring to sample from $p$, which makes it computationally more efficient as only one kernel matrix needs to be computed, but it only achieves the desired level $\alpha$ asymptotically (Shao, 2010; Leucht & Neumann, 2013; Chwialkowski et al., 2016; 2014). Note that we cannot obtain a non-asymptotic level for the wild bootstrap by relying on the result of Romano & Wolf (2005, Lemma 1) as done in the two-sample framework by Fromont et al. (2013) and Schrab et al. (2021). This is because in our case $\bar{K}_k$ and $\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N)$ are not exchangeable variables under the null hypothesis, due to the asymmetry of the KSD statistic with respect to $p$ and $q$.

Having discussed control of the probability of type I error of the single test $\Delta_{\alpha}^k$, we now provide a condition on $\|p-q\|_2$ which ensures that the probability of type II error is controlled by some $\beta \in (0, 1)$. The smallest such value of $\|p-q\|_2$, provided that $p-q$ lies in some given class of functions, is called the uniform separation rate (Baraud, 2002).

**Theorem 3.1.** Let $\psi := p - q$ and assume that $\max(\|p\|_{\infty}, \|q\|_{\infty}) \le M$. Consider $C_k$ as defined in Equation (3), $\alpha \in (0, e^{-1})$, $\beta \in (0, 1)$ and $B_1 \in \mathbb{N}$ satisfying $B_1 \ge \frac{3}{\alpha^2}(\ln(\frac{8}{\beta}) + \alpha(1-\alpha))$. There exists a positive constant $C$ such that the condition

$$ \|\psi\|_2^2 \ge \|\psi - h_{p,k} \diamond \psi\|_2^2 + C \log\left(\frac{1}{\alpha}\right) \frac{\sqrt{C_k}}{\beta N} $$

guarantees control over the probability of type II error

$$ \mathbb{P}_q(\Delta_{\alpha}^{k}(\mathbb{X}_N) = 0) \le \beta. $$

*Proof.* We can obtain this result by following the reasoning of Schrab et al. (2021, Theorem 5) for the two-sample case. We highlight the main differences in the proof. Similarly to the result of Schrab et al. (2021, Lemma 2), using Chebyshev's inequality, we obtain the sufficient condition for power

$$ \text{KSD}_{p,k}^2(q) \ge \sqrt{\frac{2}{\beta}\text{var}(\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N)) + \hat{q}_{1-\alpha}^k} $$

We then need to upper bound the variance and quantile terms. One can obtain using the explicit variance formula (Lee, 1990) that

$$ \text{var}(\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N)) \le C \left( \frac{\sigma_1^2}{N} + \frac{\sigma_2^2}{N^2} \right) $$

where

$$ \sigma_1^2 := \mathbb{E}_X [ (\mathbb{E}_{X'} [h_{p,k}(X, X')])^2 ] \le M \|h_{p,k} \diamond \psi\|_2^2 $$

by Stein's identity and

$$ \sigma_2^2 := \mathbb{E}_{X,X'} [h_{p,k}(X, X')^2] = C_k. $$

A similar reasoning to the one of Schrab et al. (2021, Proposition 4) gives us

$$ \hat{q}_{1-\alpha}^k \le C \log \left( \frac{1}{\alpha} \right) \frac{\sqrt{C_k}}{N} $$

provided that the estimated quantile is computed using $B_1 \ge \frac{3}{\alpha^2}(\ln(\frac{8}{\beta}) + \alpha(1-\alpha))$ values. The upper bounds involving the $C_k$ constant can then be combined using the basic inequalities $2\sqrt{xy} \le x+y$ and $\sqrt{x+y} \le \sqrt{x} + \sqrt{y}$, together with the fact that

$$ \begin{aligned} \text{KSD}_{p,k}^2(q) &= \langle \psi, h_{p,k} \diamond \psi \rangle \\ &= \frac{1}{2} (\|\psi\|_2^2 + \|h_{p,k} \diamond \psi\|_2^2 - \|\psi - h_{p,k} \diamond \psi\|_2^2) \end{aligned} $$

to obtain the desired result. A similar 'trick' is also present in the works of Fromont et al. (2013) and Albert et al. (2019). □

<footnote n="1">¹We do not write explicitly the dependence of $\hat{q}_{1-\alpha}^k$ on other variables, but those are implicitly considered when writing probabilistic statements.</footnote>
where the Stein kernel $h_{p,k}: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ is defined as

$$h_{p,k}(x, y) := (\nabla \log p(x)^{\top} \nabla \log p(y)) k(x, y) + \nabla \log p(y)^{\top} \nabla_x k(x, y) \\ + \nabla \log p(x)^{\top} \nabla_y k(x, y) + \sum_{i=1}^{d} \frac{\partial}{\partial x_i \partial y_i} k(x, y)$$

and satisfies the *Stein identity* $\mathbb{E}_p[h_{p,k}(X, \cdot)] = 0$. Additional background details on the KSD are presented in Appendix C. A quadratic-time KSD estimator can be computed as the U-statistic (Hoeffding, 1992)

$$\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) := \frac{1}{N(N-1)} \sum_{1 \le i \ne j \le N} h_{p,k}(X_i, X_j). \quad (1)$$

In this work, the score of the model density $p$ is always known, we do not always explicitly write the dependence on $p$ of all the variables we consider. We assume that the kernel $k$ is such that

$$\text{KSD}_{p,k}^2(q) = \mathbb{E}_{q,q}[h_{p,k}(X, Y)] < \infty \quad \text{and} \quad C_k := \mathbb{E}_{q,q}[h_{p,k}(X, Y)^2] < \infty. \quad (2)$$

We now address the requirements for consistency (i.e. test power converges to 1 as the sample size goes to $\infty$) of the Stein test (Chwialkowski et al., 2016, Theorem 2.2): we assume that the kernel $k$ is

$C_0$-universal (Carmeli et al., 2010, Definition 4.1) and that $\mathbb{E}_q[\|\nabla(\log \frac{p(X)}{q(X)})\|_2^2] < \infty$.

We use the notations $\mathbb{P}_p$ and $\mathbb{P}_q$ to denote the probability under the model distribution $p$ and under the data distribution $q$, respectively. Given a kernel $\kappa: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ and a function $f: \mathbb{R}^d \to \mathbb{R}$ in $L^2(\mathbb{R}^d)$, we consider the integral transform $T_\kappa$ defined as

$$(T_\kappa f)(y) := \int_{\mathbb{R}^d} \kappa(x, y)f(x) \, dx$$

for $y \in \mathbb{R}^d$. When the kernel $\kappa$ is translation-invariant, the integral transform corresponds to a convolution. However, the lack of translation invariance of the Stein kernel introduces new challenging problems. First, working the integral transform of the Stein kernel is more complicated since it does not correspond to a simple convolution. Second, for the expectation of the Stein kernel squared, $C_k := \mathbb{E}_{q,q}[h_{p,k}(X, Y)^2]$, it is not possible to extract the bandwidth parameter $\lambda$ outside of the expectation as it is the case when using the usual kernel directly as for the MMD and HSIC.

# 3 Construction of tests and bounds

We now introduce the single and aggregated KSD tests. We show that these control the probability of type I error as desired, and provide conditions for the control of the probability of type II error.

## 3.1 Single test

We first construct a KSD test for a fixed kernel $k$ as proposed by Chwialkowski et al. (2016) and Liu et al. (2016). To estimate the test threshold, we can either use a wild bootstrap (Shao, 2010; Leucht and Neumann, 2013; Fromont et al., 2012; Chwialkowski et al., 2014) or a parametric bootstrap (Stute et al., 1993). Both methods work by simulating sampling values $(\bar{K}_k^1, \dots, \bar{K}_k^{B_1})$ from the (asymptotic) distribution of $\widehat{\text{KSD}}_{p,k}^2$ under the null hypothesis and estimating the $(1-\alpha)$-quantile for $\alpha \in (0, 1)$ using a Monte Carlo approximation<fnref n="1" />

$$\hat{q}_{1-\alpha}^k := \inf \left\{ u \in \mathbb{R} : 1 - \alpha \le \frac{1}{B_1 + 1} \sum_{b=1}^{B_1+1} 1(\bar{K}_k^b \le u) \right\} = \bar{K}_k^{\alpha(B_1+1)(1-\alpha)} \quad (3)$$

where $\bar{K}_k^{\alpha 1} \le \dots \le \bar{K}_k^{\alpha B_1+1}$ are the sorted elements $(\bar{K}_k^1, \dots, \bar{K}_k^{B_1+1})$ with $\bar{K}_k^{B_1+1} := \widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N)$. The single test is then defined as (a test function outputs 1 when the null is rejected and 0 otherwise)

$$\Delta_{\alpha}^{k}(\mathbb{X}_{N}) := 1(\widehat{\text{KSD}}_{p,k}^{2}(\mathbb{X}_{N}) > \hat{q}_{1-\alpha}^{k}).$$

<footnote n="1">¹We do not write explicitly the dependence of $\hat{q}_{1-\alpha}^k$ on other variables, but those are implicitly considered when writing probabilistic statements.</footnote>
## 3.2. Aggregated test

We can now introduce our aggregated test, which is motivated by the earlier works of Fromont et al. (2012; 2013), Albert et al. (2019), and Schrab et al. (2021).

We compute $\tilde{K}_k^1, \dots, \tilde{K}_k^{B_2}$ further simulated KSD values from the null hypothesis obtained using either a parametric bootstrap or a wild bootstrap as in Equations (4) or (5), respectively. Consider a finite collection of kernels $\mathcal{K}$ satisfying the properties presented in Section 2. We construct an aggregated test $\Delta_{\alpha}^{\mathcal{K}}$, called KSDAGG, which rejects the null hypothesis if one of the single tests $(\Delta_{u_{\alpha} w_k}^k)_{k \in \mathcal{K}}$ rejects it, that is

$$ \Delta_{\alpha}^{\mathcal{K}}(\mathbb{X}_N) := \mathbb{1}(\Delta_{u_{\alpha} w_k}^k(\mathbb{X}_N) = 1 \text{ for some } k \in \mathcal{K}). $$

The levels of the single tests are adjusted to ensure the aggregated test has the prescribed level $\alpha$. This adjustment is performed by introducing positive weights $(w_k)_{k \in \mathcal{K}}$ satisfying $\sum_{k \in \mathcal{K}} w_k \le 1$ and some correction

$$ u_{\alpha} := \sup \left\{ u \in (0, \min_{k \in \mathcal{K}} w_k^{-1}) : \hat{P}_u \le \alpha \right\} \quad (6) $$

where

$$ \hat{P}_u := \frac{1}{B_2} \sum_{b=1}^{B_2} \mathbb{1} \left( \max_{k \in \mathcal{K}} \left( \tilde{K}_k^{\prime b} - \tilde{K}_k^{\prime \cdot [B_1(1-uw_k)]} \right) > 0 \right) $$

is a Monte Carlo approximation of the probability of type I error of our aggregated test with correction $u$

$$ P_u := \mathbb{P}_p \left( \max_{k \in \mathcal{K}} \left( \widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) - \hat{q}_{1-uw_k}^k \right) > 0 \right). $$

In practice, we work with a finite collection of kernels, and thus the extrema over $\mathcal{K}$ are achieved. To compute $u_{\alpha}$, we estimate the supremum in Equation (6) by performing $B_3$ steps of the bisection method, as proposed by Schrab et al. (2021). Detailed pseudocode for KSDAGG is provided in Algorithm 1.

We verify in the next proposition that performing this correction indeed ensures that our aggregated test $\Delta_{\alpha}^{\mathcal{K}}$ has the prescribed level $\alpha$.

**Proposition 3.2.** For $\alpha \in (0, 1)$ and $B_1, B_2, B_3 \in \mathbb{N} \setminus \{0\}$, the aggregated test $\Delta_{\alpha}^{\mathcal{K}}$ satisfies

$$ \mathbb{P}_p(\Delta_{\alpha}^{\mathcal{K}}(\mathbb{X}_N) = 1) \le \alpha $$

non-asymptotically when using a parametric bootstrap and asymptotically when using a wild bootstrap.

*Proof.* This result follows from the proof provided by Schrab et al. (2021, Proposition 8) and the facts that the parametric bootstrap (Stute et al., 1993) and wild bootstrap (Chwialkowski et al., 2016) have non-asymptotic and asymptotic levels $\alpha$, respectively. □

We now provide guarantees for the power of our aggregated test KSDAGG in terms of its uniform separation rate.

### Algorithm 1 KSDAGG

**Inputs:** samples $\mathbb{X}_N = (x_i)_{i=1}^N$, density $p$ or score function $\nabla \log p(\cdot)$, finite kernel collection $\mathcal{K}$, weights $(w_k)_{k \in \mathcal{K}}$, level $\alpha \in (0, e^{-1})$, estimation parameters $B_1, B_2, B_3 \in \mathbb{N}$, parametric or wild bootstrap

**for** $k \in \mathcal{K}$ **do**

Compute $(\tilde{K}_k^b)_{1 \le b \le B_1}$ as in Equations (4) or (5)

Sort in ascending order to obtain $(\tilde{K}_k^{\bullet b})_{1 \le b \le B_1}$

Compute $(\tilde{K}_k^b)_{1 \le b \le B_2}$ as in Equations (4) or (5)

**end for**

$$ u_{\min} = 0 $$

$$ u_{\max} = \min_{k \in \mathcal{K}} w_k^{-1} $$

**for** $t = 1, \dots, B_3$ **do**

$$ u = \frac{1}{2}(u_{\min} + u_{\max}) $$

$$ \hat{P}_u = \frac{1}{B_2} \sum_{b=1}^{B_2} \mathbb{1} \left( \max_{k \in \mathcal{K}} \left( \tilde{K}_k^b - \tilde{K}_k^{\bullet [B_1(1-uw_k)]} \right) > 0 \right) $$

**if** $\hat{P}_u \le \alpha$ **then**

$$ u_{\min} = u $$

**else**

$$ u_{\max} = u $$

**end if**

**end for**

$$ u_{\alpha} = u_{\min} $$

**if** $\max_{k \in \mathcal{K}} (\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) - \tilde{K}_k^{\bullet [B_1(1-u_{\alpha}w_k)]}) > 0$ **then**

**return** 1 (reject $\mathcal{H}_0$)

**else**

**return** 0 (fail to reject $\mathcal{H}_0$)

**end if**

**Theorem 3.3.** Let $\psi := p - q$ denote the difference in densities and assume that $\max(\|p\|_{\infty}, \|q\|_{\infty}) \le M$. Consider a collection of kernels $\mathcal{K}$ with associated positive weights $(w_k)_{k \in \mathcal{K}}$ satisfying $\sum_{k \in \mathcal{K}} w_k \le 1$. Consider $C_k$ as defined in Equation (3), $\alpha \in (0, e^{-1})$, $\beta \in (0, 1)$ and $B_1, B_2, B_3 \in \mathbb{N}$ satisfying $B_1 \ge \frac{3}{\alpha^2}(\ln(\frac{8}{\beta}) + \alpha(1 - \alpha))$, $B_2 \ge \frac{8}{\alpha^2}\ln(\frac{2}{\beta})$ and $B_3 \ge \log_2(\frac{4}{\alpha}\min_{k \in \mathcal{K}} w_k^{-1})$. There exists a positive constant $C$ such that if

$$ \|\psi\|_2^2 \ge \min_{k \in \mathcal{K}} \left( \|\psi - h_{p,k} \diamond \psi\|_2^2 + C \log \left( \frac{1}{\alpha w_k} \right) \frac{\sqrt{C_k}}{\beta N} \right) $$

then the probability of type II error of $\Delta_{\alpha}^{\mathcal{K}}$ is controlled,

$$ \mathbb{P}_q(\Delta_{\alpha}^{\mathcal{K}}(\mathbb{X}_N) = 0) \le \beta. $$

*Proof.* Applying the result of Schrab et al. (2021, Theorem 9), we obtain that the uniform separation rate of our aggregated test $\Delta_{\alpha}^{\mathcal{K}}(\mathbb{X}_N)$ is lower bounded by the minimum over the uniform separation rates of the single tests $(\Delta_{\alpha w_k}^k)_{k \in \mathcal{K}}$. Combining this result with Theorem 3.1 concludes the proof. □
For the parametric bootstrap, we directly draw new samples $(X'_i)_{i=1}^N$ from the model distribution $p$ (it might not always be possible to so) and compute the KSD

$$ \bar{K}_k := \frac{1}{N(N-1)} \sum_{1 \le i \ne j \le N} h_{p,k}(X'_i, X'_j). \quad (4) $$

For the wild bootstrap, we first generate $n$ i.i.d. Rademacher random variables $\epsilon_1, \dots, \epsilon_n$, each taking value in $\{-1, 1\}$, and then compute

$$ \bar{K}_k := \frac{1}{N(N-1)} \sum_{1 \le i \ne j \le N} \epsilon_i \epsilon_j h_{p,k}(X_i, X_j). \quad (5) $$

By repeating either procedure $B_1$ times, we obtain the bootstrapped samples $\bar{K}_k^1, \dots, \bar{K}_k^{B_1}$.

Since it uses samples from the model $p$, the parametric bootstrap (Stute et al., 1993) results in a test with non-asymptotic level $\alpha$. This comes at the cost of being computationally more expensive and assuming that we are able to sample from $p$ (which may be out of reach in some settings). Conversely, the wild bootstrap has the advantage of not requiring to sample from $p$, which makes it computationally more efficient as only one kernel matrix needs to be computed, but it only achieves the desired level $\alpha$ asymptotically (Shao, 2010; Leucht and Neumann, 2013; Leucht, 2012; Chwialkowski et al., 2014, 2016) assuming Lipschitz continuity of $h_{p,k}$ (see Appendix D for details). Note that we cannot obtain a non-asymptotic level for the wild bootstrap by relying on the result of Romano and Wolf (2005, Lemma 1) as done in the two-sample framework by Fromont et al. (2013) and Schrab et al. (2021). This is because in our case $\bar{K}_k$ and $\widehat{\text{KSD}}_{p,k}^2(X_N)$ are not exchangeable variables under the null hypothesis, due to the asymmetry of the KSD statistic with respect to $p$ and $q$.

Having discussed control of the probability of type I error of the single test $\Delta_{\alpha}^k$, we now provide a condition on $\|p - q\|_2$ which ensures that the probability of type II error is controlled by some $\beta \in (0, 1)$. The smallest such value of $\|p - q\|_2$, provided that $p - q$ lies in some given class of regular functions, is called the uniform separation rate (Ingster, 1987, 1989, 1993a,b; Baraud, 2002).

**Theorem 3.1.** Suppose the assumptions listed in Appendix A.2 hold, and let $\psi := p - q$. There exists a positive constant $C$ such that the condition

$$ \|\psi\|_2^2 \ge \|\psi - T_{h_{p,k}}\psi\|_2^2 + C \log\left(\frac{1}{\alpha}\right) \frac{\sqrt{C_k}}{\beta N} $$

guarantees control over the probability of type II error, such that $\mathbb{P}_q(\Delta_{\alpha}^k(X_N) = 0) \le \beta$.

Theorem 3.1, which is proved in Appendix I.1, provides a power guaranteeing condition consisting of two terms. The first term $\|\psi - T_{h_{p,k}}\psi\|_2^2$ indicates the size of the effect of the Stein integral transform operator on the difference in densities $\psi := p - q$, it is a measure of distance from the null (where this quantity is zero). The second term $C \log(1/\alpha)(\beta N)^{-1}\sqrt{C_k}$ is obtained from upper bounding the variance of the KSD U-statistic, it depends on the expectation of the squared Stein kernel $C_k := \mathbb{E}_{q,q}[h_{p,k}(X, Y)^2]$. This second term also controls the quantile of the test.

## 3.2 Aggregated test

We can now introduce our aggregated test, which is motivated by the earlier works of Fromont et al. (2012, 2013), Albert et al. (2022), and Schrab et al. (2021) for two-sample and independence testing.

We compute $\tilde{K}_k^1, \dots, \tilde{K}_k^{B_2}$ further KSD values simulated from the null hypothesis obtained using either a parametric bootstrap or a wild bootstrap as in Equations (4) or (5), respectively. We consider a finite collection $\mathcal{K}$ of kernels satisfying the properties presented in Section 2. We construct an aggregated test $\Delta_{\alpha}^{\mathcal{K}}$, called KSDAGG, which rejects the null hypothesis if one of the single tests $(\Delta_{u_{\alpha}w_k}^k)_{k \in \mathcal{K}}$ rejects it, that is

$$ \Delta_{\alpha}^{\mathcal{K}}(\mathbb{X}_N) := \mathbb{1}(\Delta_{u_{\alpha}w_k}^k(\mathbb{X}_N) = 1 \text{ for some } k \in \mathcal{K}). $$

The levels of the single tests are adjusted to ensure the aggregated test has the prescribed level $\alpha$. This adjustment is performed by introducing positive weights $(w_k)_{k \in \mathcal{K}}$ satisfying $\sum_{k \in \mathcal{K}} w_k \le 1$ and some correction

$$ u_{\alpha} := \sup \left\{ u \in (0, \min_{k \in \mathcal{K}} w_k^{-1}) : \hat{P}_u \le \alpha \right\} \quad (6) $$
We observe that the aggregation procedure allows to achieve the smallest uniform separation rate of the single tests $(\Delta_{\alpha}^k)_{k \in \mathcal{K}}$ up to some logarithmic weighting term $\log(1/w_k)$.

### 3.3. Bandwidth selection

A specific application of the setting we have considered is the problem of bandwidth selection for a fixed kernel. Given a kernel $k: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$, the function

$$k_{\lambda}(x, y) := k\left(\frac{x}{\lambda}, \frac{y}{\lambda}\right)$$

is also a kernel for any bandwidth $\lambda > 0$. A common example is the Gaussian kernel, for which we have $k(x, y) = \exp(-\|x - y\|_2^2)$ and

$$k_{\lambda}(x, y) = \exp\left(-\frac{\|x - y\|_2^2}{\lambda^2}\right).$$

As shown by Gorham & Mackey (2017), a more appropriate kernel for goodness-of-fit testing using the KSD is the IMQ (inverse multiquadric) kernel, which is defined with $k(x, y) = (1 + \|x - y\|_2^2)^{-\beta_k}$ for a fixed parameter $\beta_k \in (0, 1)$ as

$$\begin{aligned} k_{\lambda}(x, y) &= \left(1 + \frac{\|x - y\|_2^2}{\lambda^2}\right)^{-\beta_k} \\ &= \lambda^{2\beta_k} \left(\lambda^2 + \|x - y\|_2^2\right)^{-\beta_k} \\ &\propto \left(\lambda^2 + \|x - y\|_2^2\right)^{-\beta_k} \end{aligned} \quad (7)$$

which is the well-known form of the IMQ kernel with parameters $\lambda > 0$ and $\beta_k \in (0, 1)$. Note that it is justified to consider the kernel up to a multiplicative constant because our single and aggregated tests are invariant under this kernel transformation.

In practice, as suggested by Gretton et al. (2012a) the bandwidth is often set to a heuristic such as the median or the standard deviation of the $L^2$-distances between the samples $(X_i)_{i=1}^N$, however, these are arbitrary choices with no theoretical guarantees. Another common approach proposed by Gretton et al. (2012b); Liu et al. (2020) is to resort to data splitting in order to select a bandwidth on held-out data, by maximising for a proxy for asymptotic power (see Section 4.1 for details). Both methods were originally proposed for the two-sample problem, but extend straightforwardly to the goodness-of-fit setting.

By considering a kernel collection $\mathcal{K}_{\Lambda} = \{k_{\lambda} : \lambda \in \Lambda\}$ for a collection of bandwidths $\Lambda$, we can use our aggregated test KSDAGG to test multiple bandwidths using all the data and without resorting to arbitrary heuristics. We now obtain an expression for the uniform separation rate of $\Delta_{\alpha}^{\mathcal{K}_{\Lambda}}$ in terms of the bandwidth $\lambda$.

**Corollary 3.4.** Consider $\alpha \in (0, e^{-1})$, $\beta \in (0, 1)$ and $B_1, B_2, B_3 \in \mathbb{N}$ satisfying the conditions of Theorem 3.3 and assume that $\max(\|p\|_{\infty}, \|q\|_{\infty}) \le M$. Given a collection $\Lambda$ of positive bandwidths with associated positive weights $(w_{\lambda})_{\lambda \in \Lambda}$ satisfying $\sum_{\lambda \in \Lambda} w_{\lambda} \le 1$, we consider $\mathcal{K}_{\Lambda} = \{k_{\lambda} : \lambda \in \Lambda\}$. There exists a positive constant $C$ such that the condition

$$\|\psi\|_2^2 \ge \min_{\lambda \in \Lambda} \left( \|\psi - h_{p,k_{\lambda}} \diamond \psi\|_2^2 + C \log \left( \frac{1}{\alpha w_{\lambda}} \right) \frac{\sqrt{C_{k_{\lambda}}}}{\beta N} \right)$$

ensures control over the probability of type II error of the aggregated test

$$\mathbb{P}_q(\Delta_{\alpha}^{\mathcal{K}_{\Lambda}}(\mathbb{X}_N) = 0) \le \beta.$$

*Proof.* We simply apply Theorem 3.3 to the collection of kernels $\mathcal{K}_{\Lambda}$. $\square$

We do not impose any restrictions on $\psi := p - q$ such as assuming it belongs to a specific regularity class. For this reason, our result holds more generally but the dependence on $\lambda$ in the terms $\|\psi - h_{p,k_{\lambda}} \diamond \psi\|_2^2$ and $-\log(\alpha w_{\lambda}) (\beta N)^{-1} \sqrt{C_{k_{\lambda}}}$ is not explicit. For a particular regularity class, one can obtain a uniform separation rate $N^{-R}$ for some $R > 0$ by choosing appropriate collections of bandwidths and weights (depending on $N$) such that the two terms have matching orders of $N$.

## 4. Implementation and experiments

We consider three different experiments based on a Gamma one-dimensional distribution, a Bernoulli Gaussian Restricted Boltzmann Machine, and a Normalizing Flow for the MNIST dataset. We compare our proposed aggregated test KSDAGG against three alternatives: the KSD test which uses the median bandwidth, a test which splits the data to select an ‘optimal’ bandwidth according to a proxy for asymptotic test power, and a test which uses extra data for bandwidth selection. The ‘extra data’ test is designed simply to provide a best case oracle for the ‘asymptotic test power’ bandwidth selection procedure, and cannot be used in practice (i.e. any extra samples from $q$ would normally be incorporated into the sample being tested). In order to ensure that our tests always have correct levels for all bandwidth values, dimensions and sample sizes, we use the parametric bootstrap in our experiments.

### 4.1. Alternative bandwidth selection approaches

Gretton et al. (2012a) propose to use the median heuristic as kernel bandwidth, consisting in the median of the $L^2$-distances between the samples given by

$$\lambda_{\text{med}} := \text{median}\{\|x_i - x_j\|_2 : 0 \le i < j \le N\}.$$
Algorithm 1 KSDAGG

**Inputs:** samples $\mathbb{X}_N = (x_i)_{i=1}^N$, density $p$ or score $\nabla \log p(\cdot)$, finite kernel collection $\mathcal{K}$, weights $(w_k)_{k \in \mathcal{K}}$, level $\alpha \in (0, e^{-1})$, estimation parameters $B_1, B_2, B_3 \in \mathbb{N}$, parametric or wild bootstrap

**Output:** 0 (fail to reject $\mathcal{H}_0$) or 1 (reject $\mathcal{H}_0$)

**Algorithm:**

**for** $k \in \mathcal{K}$ **do**

compute $\bar{K}_k^{B_1+1} := \widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N)$ as in Equation (1)

compute $(\bar{K}_k^b)_{1 \le b \le B_1}$ as in Equations (4) or (5)

sort $(\bar{K}_k^b)_{1 \le b \le B_1+1}$ in ascending order to obtain $(\bar{K}_k^{\bullet b})_{1 \le b \le B_1+1}$

compute $(\tilde{\bar{K}}_k^b)_{1 \le b \le B_2}$ as in Equations (4) or (5)

$$u_{\min} = 0, \quad u_{\max} = \min_{k \in \mathcal{K}} w_k^{-1}$$

**for** $t = 1, \dots, B_3$ **do**

$$u = \frac{1}{2}(u_{\min} + u_{\max}), \quad \hat{P}_u = \frac{1}{B_2} \sum_{b=1}^{B_2} 1 \left( \max_{k \in \mathcal{K}} \left( \tilde{\bar{K}}_k^b - \bar{K}_k^{\bullet \left\lceil (B_1+1)(1-uw_k) \right\rceil} \right) > 0 \right)$$

if $\hat{P}_u \le \alpha$ then $u_{\min} = u$ else $u_{\max} = u$

$$u_{\alpha} = u_{\min}$$

if $\max_{k \in \mathcal{K}} (\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) - \bar{K}_k^{\bullet \left\lceil (B_1+1)(1-u_{\alpha}w_k) \right\rceil}) > 0$ then return 1 else return 0

**Time complexity:** $O(|\mathcal{K}|(B_1 + B_2) N^2)$

where

$$\hat{P}_u := \frac{1}{B_2} \sum_{b=1}^{B_2} 1 \left( \max_{k \in \mathcal{K}} \left( \tilde{\bar{K}}_k^b - \bar{K}_k^{\bullet \left\lceil (B_1+1)(1-uw_k) \right\rceil} \right) > 0 \right)$$

is a Monte Carlo approximation of the type I error probability of the aggregated test with correction $u$

$$P_u := \mathbb{P}_p \left( \max_{k \in \mathcal{K}} \left( \widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) - \hat{q}_{1-uw_k}^k \right) > 0 \right).$$

To compute $u_\alpha$, we estimate the supremum in Equation (6) by performing $B_3$ steps of the bisection method, the theoretical results account for this extra approximation. Detailed pseudocode for KSDAGG is provided in Algorithm 1, and details regarding our aggregation procedure are provided in Appendix G. We discuss potential limitations of KSDAgg in Appendix H.

We verify in the next proposition that performing this correction indeed ensures that our aggregated test $\Delta_\alpha^\mathcal{K}$ has the prescribed level $\alpha$.

**Proposition 3.2.** For $\alpha \in (0, 1)$ and a collection of kernels $\mathcal{K}$, the aggregated test $\Delta_\alpha^\mathcal{K}$ satisfies

$$\mathbb{P}_p(\Delta_\alpha^\mathcal{K}(\mathbb{X}_N) = 1) \le \alpha$$

asymptotically using a wild bootstrap (with Lipschitz continuity of $h_{p,k}$ required) and non-asymptotically using a parametric bootstrap.

The proof of Proposition 3.2 is presented in Appendix I.2. Details about the asymptotic result for the wild bootstrap case are reported in Appendix D. We now provide guarantees for the power of our aggregated test KSDAGG in terms of its uniform separation rate.

**Theorem 3.3.** Suppose the assumptions listed in Appendix A.3 hold, and let $\psi := p - q$. There exists a positive constant $C$ such that if

$$\|\psi\|_2^2 \ge \min_{k \in \mathcal{K}} \left( \|\psi - T_{h_{p,k}} \psi\|_2^2 + C \log \left( \frac{1}{\alpha w_k} \right) \frac{\sqrt{C_k}}{\beta N} \right)$$

then the probability of type II error of $\Delta_\alpha^\mathcal{K}$ is controlled by $\beta$, that is, $\mathbb{P}_q(\Delta_\alpha^\mathcal{K}(\mathbb{X}_N) = 0) \le \beta$.

The proof Theorem 3.3 can be found in Appendix I.3, it relies on upper bounding the probability of the intersection of some events by the minimum of the probabilities of each event, and on applying Theorem 3.1 after having verified that its assumptions are satisfied for the tests with adjusted levels. We observe that the aggregation procedure allows to achieve the smallest uniform separation rate of the single tests $(\Delta_\alpha^k)_{k \in \mathcal{K}}$ up to some logarithmic weighting term $\log(1/w_k)$.
Gretton et al. (2012b) first proposed, for the two-sample problem using a linear-time MMD estimator, to split the data and to use half of it to select an ‘optimal’ bandwidth which maximises a proxy for asymptotic power. This procedure was then extended to quadratic-time estimators and to the goodness-of-fit framework by Jitkrittum et al. (2017), Sutherland et al. (2017) and Liu et al. (2020). These strategies rely on the asymptotic normality of the test statistic under $H_a$. In our setting, the asymptotic power proxy to maximise is the ratio

$$ \widehat{KSD}_{p,k}^2(\mathbb{X}_N) / \widehat{\sigma}_{H_a} \qquad (8) $$

where $\widehat{\sigma}_{H_a}^2$ is a regularised positive estimator of the asymptotic variance of $\widehat{KSD}_{p,k}^2$ under $H_a$ (Liu et al., 2020, Equation 5) which can be computed as

$$ \frac{4}{N^3} \sum_{i=1}^{N} \left( \sum_{j=1}^{N} H_{i,j} \right)^2 - \frac{4}{N^4} \left( \sum_{i=1}^{N} \sum_{j=1}^{N} H_{i,j} \right)^2 $$

with $H_{i,j} := h_{p,k}(X_i, X_j)$. In our experiments, we also consider a test which has access to $N$ extra samples drawn from $q$ to select an ‘optimal’ bandwidth to run the KSD test on the original $N$ samples $\mathbb{X}_N$. This test is interesting to compare to because it uses an ‘optimal’ bandwidth without being detrimental to power.

## 4.2. Experimental details

In our experiments, we use collections of bandwidths of the form

$$ \Lambda(\ell_{-}, \ell_{+}) := \{2^{i}\lambda_{\text{med}} : i = \ell_{-}, \dots, \ell_{+}\} $$

for the median bandwidth $\lambda_{\text{med}}$ and integers $\ell_{-} < \ell_{+}$ with uniform weights

$$ w_{\lambda} := \frac{1}{\ell_{+} - \ell_{-} + 1}. $$

For the tests which split the data, we select the bandwidth out of the collection $\Lambda(\ell_{-}, \ell_{+})$ which maximises the power proxy in Equation (8). All our experiments are run with level $\alpha = 0.05$ using the IMQ kernel defined in Equation (7) with parameter $\beta_k = 0.5$. We use a parametric bootstrap with $B_1 = B_2 = 500$ simulated KSD values to compute the adjusted test thresholds and $B_3 = 50$ steps of bisection method to estimate the correction $u_\alpha$ in Equation (6). To estimate the probability of rejecting the null hypothesis, we average the test outputs across 200 repetitions.

## 4.3. Gamma distribution

For our first experiment, we consider a one-dimensional Gamma distribution with shape parameter 5 and scale

parameter 5 as the model $p$. We draw 500 samples from a Gamma distribution with the same scale parameter 5 and with a shifted shape parameter $5 + s$ for $s \in \{0, 0.1, 0.2, 0.3, 0.4\}$. We consider the collection of bandwidths $\Lambda(0, 10)$.

<figure><img src="image_2.png" /><figcaption>Figure 1. Gamma distribution experiment.</figcaption></figure>

The results we obtained are presented in Figure 1. We observe that all tests have the prescribed level 0.05 under the null hypothesis, corresponding to the case $s = 0$. As the shift parameter $s$ increases, the two densities $p$ and $q$ become more different and rejection of the null becomes an easier task, thus the test power increases. Our aggregated test KSDAGG achieves the same power as the ‘best case’ bound on the performance of the asymptotic power heuristic, yielded by the splitting test with extra data. The median test obtains only slightly lower power, the closeness in power can be explained by the fact that this one-dimensional problem is an easy one. We note that the normal splitting test has significantly lower power: this is because, even though it uses an ‘optimal’ bandwidth, it is then run on only half the data, which results in a loss of power.

## 4.4. Bernoulli-Gaussian Restricted Boltzmann Machine

As first considered by Liu et al. (2016) for goodness-of-fit testing using the KSD, we consider a Bernoulli-Gaussian Restricted Boltzmann Machine. It is a graphical model with a binary hidden variable $h \in \{-1, 1\}^{d_h}$ and a continuous observable variable $x \in \mathbb{R}^d$. Those variables have joint density

$$ p(x, h) = \frac{1}{Z} \exp \left( \frac{1}{2} x^{\top} B h + b^{\top} x + c^{\top} h - \frac{1}{2} \|x\|_{2}^{2} \right) $$
### 3.3 Bandwidth selection

A specific application of the setting we have considered is the problem of bandwidth selection for a fixed kernel. Given a kernel $k: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$, the function

$$k_{\lambda}(x, y) := k\left(\frac{x}{\lambda}, \frac{y}{\lambda}\right)$$

is also a kernel for any bandwidth $\lambda > 0$. A common example is the Gaussian kernel, for which we have $k(x, y) = \exp(-\|x - y\|_2^2)$ and $k_{\lambda}(x, y) = \exp(-\|x - y\|_2^2/\lambda^2)$. As shown by Gorham and Mackey (2017), a more appropriate kernel for goodness-of-fit testing using the KSD is the IMQ (inverse multiquadric) kernel, which is defined with $k(x, y) = (1 + \|x - y\|_2^2)^{-\beta_k}$ for a fixed parameter $\beta_k \in (0, 1)$ as

$$k_{\lambda}(x, y) = \left(1 + \frac{\|x - y\|_2^2}{\lambda^2}\right)^{-\beta_k} = \lambda^{2\beta_k} \left(\lambda^2 + \|x - y\|_2^2\right)^{-\beta_k} \propto \left(\lambda^2 + \|x - y\|_2^2\right)^{-\beta_k} \quad (7)$$

which is the well-known form of the IMQ kernel with parameters $\lambda > 0$ and $\beta_k \in (0, 1)$. Note that it is justified to consider the kernel up to a multiplicative constant because the single and aggregated tests are invariant under this kernel transformation.

In practice, as suggested by Gretton et al. (2012a), the bandwidth is often set to a heuristic such as the median or the standard deviation of the $L^2$-distances between the samples $(X_i)_{i=1}^N$, however, these are arbitrary choices with no theoretical guarantees. Another common approach proposed by Gretton et al. (2012b) for the linear-time setting, and extended to the quadratic-time setting by Liu et al. (2020), is to resort to data splitting in order to select a bandwidth on held-out data, by maximising for a proxy for asymptotic power (see Section 4.1 for details). Both methods were originally proposed for the two-sample problem, but extend straightforwardly to the goodness-of-fit setting.

By considering a kernel collection $K_{\Lambda} = \{k_{\lambda} : \lambda \in \Lambda\}$ for a collection of bandwidths $\Lambda$, we can use our aggregated test KSDAGG to test multiple bandwidths using all the data and without resorting to arbitrary heuristics. We now obtain an expression for the uniform separation rate of $\Delta_{\alpha}^{K_{\Lambda}}$ in terms of the bandwidths $\lambda \in \Lambda$.

**Corollary 3.4.** Suppose the assumptions listed in Appendix A.3 hold for $K = K_{\Lambda} = \{k_{\lambda} : \lambda \in \Lambda\}$, and let $\psi := p - q$. There exists a positive constant $C$ such that the condition

$$\|\psi\|_2^2 \ge \min_{\lambda \in \Lambda} \left( \|\psi - T_{h_{p,k_{\lambda}}} \psi\|_2^2 + C \log \left( \frac{1}{\alpha w_{\lambda}} \right) \frac{\sqrt{C_{k_{\lambda}}}}{\beta N} \right)$$

ensures control over the probability of type II error of the aggregated test $P_q(\Delta_{\alpha}^{K_{\Lambda}}(\mathbb{X}_N) = 0) \le \beta$.

Corollary 3.4 follows from applying Theorem 3.3 to the collection of kernels $K_{\Lambda}$. Our results hold with great generality as we have not imposed any restrictions on $\psi := p-q$ such as assuming it belongs to a specific regularity class. For this reason, the dependence on $\lambda$ in the terms $\|\psi - T_{h_{p,k_{\lambda}}} \psi\|_2^2$ and $\log(1/(\alpha w_{\lambda})) (\beta N)^{-1} \sqrt{C_{k_{\lambda}}}$ is not explicit. In Section 3.4, we characterise this dependence with regularity assumptions on $\psi := p-q$, and derive uniform separation rates in terms of the sample size.

### 3.4 Uniform separation rates over restricted Sobolev balls

In this section, we derive uniform separation rates with stronger assumptions on $p$ and $q$. This provides settings in which the power guaranteeing conditions of Theorems 3.1 and 3.3 are satisfied, and illustrates the interactions between the two terms in those conditions. We make the following assumptions.

* The model density $p$ is strictly positive on its connected and compact support $S \subseteq \mathbb{R}^d$.

* The score function $\nabla \log p(x)$ is continuous and bounded on $S$.

* The support of the density $q$ is a connected and compact subset of $S$.

* The kernel used is a scaled Gaussian kernel $k_{\lambda}(x, y) := \lambda^{2-d} \exp(-\|x - y\|_2^2 / \lambda^2)$.
where Z is an unknown normalizing constant. By marginalising over h, we obtain the density p of x

$$p(x) = \sum_{h \in \{-1, 1\}^{d_h}} p(x, h).$$

We can sample from it using a Gibbs sampler with 2000 burn-in iterations. We use the dimensions $d = 50$ and $d_h = 40$ as considered by Jitkrittum et al. (2017) and Grathwohl et al. (2020). Even though computing $p$ is intractable for large dimension $d_h$, the score function admits a convenient closed form

$$\nabla \log p(x) = b - x + B \frac{\exp(2(B^\top x + c)) - 1}{\exp(2(B^\top x + c)) + 1}.$$

We draw the components of $b$ and $c$ from Gaussian standard distributions and sample Rademacher variables taking values in $\{-1, 1\}$ for the elements of $B$ for the model $p$. We draw 1000 samples from a distribution $q$ which is constructed in a similar way as $p$ but with the difference that some Gaussian noise $N(0, \sigma)$ is injected into the elements of $B$. We consider the standard deviations of the perturbations $\sigma \in \{0, 0.01, 0.02, 0.03\}$. We run our experiments with the collection of bandwidths $\Lambda(-20, 0)$ and provide the results in Figure 2.

<figure><img src="image_3.png" /><figcaption>Figure 2. Bernoulli-Gaussian Restricted Boltzmann Machine experiment.</figcaption></figure>

Again, we observe that our aggregated test KSDAGG matches the power obtained by the test which uses extra data to select an 'optimal' bandwidth. This means that, in this experiment, KSDAGG obtains the same power as the 'best' single test. The difference between KSDAGG and the median heuristic test is now significant, and the splitting test obtains lowest power of the four tests. Again, all tests have well-calibrated levels ($\sigma = 0$) and increasing the noise level $\sigma$ results in more power for all the tests.

<figure><img src="image_4.png" /><figcaption>Figure 3. MNIST and Normalizing Flow digits. The top row corresponds to digits from the MNIST dataset. The digits in the middle and bottom rows have been sampled from the Normalizing Flow.</figcaption></figure>

## 4.5. MNIST Normalizing Flow

Finally, we consider a high-dimensional problem working with images in dimensions $28^2 = 784$. We consider a multi-scale Normalizing Flow (Dinh et al., 2017; Kingma & Dhariwal, 2018) which has been trained on the MNIST dataset (LeCun et al., 1998; 2010), it is a generative model which has a probability density $p$. As observed in Figure 3, some samples produced by the model can look exactly like MNIST digits, while other do not resemble digits. This Normalizing Flow has been trained to 'ideally' produce samples from the MNIST dataset. We are interested in whether or not we can detect the difference in densities. Given some images of digits, are we able to tell if those were sampled from the Normalizing Flow model? We consider the case where the images from $q$ are sampled from the model (level experiment, confirming performance under $H_0$), and the case where the samples from $q$ are drawn from the true MNIST dataset (power experiment). The experiments are run with the collection of bandwidths $\Lambda(-20, 0)$. The results are displayed in Table 1 and Figure 4.

<table><caption>Table 1. Normalizing Flow MNIST level experiment.</caption><thead><tr><th>SAMPLE SIZE</th><th>KSDAGG</th><th>KSD MEDIAN</th><th>KSD SPLIT</th><th>KSD SPLIT EXTRA DATA</th></tr></thead><tbody><tr><td>100</td><td>0.075</td><td>0.05</td><td>0.05</td><td>0.05</td></tr><tr><td>200</td><td>0.05</td><td>0.05</td><td>0.04</td><td>0.05</td></tr><tr><td>300</td><td>0.055</td><td>0.09</td><td>0.04</td><td>0.09</td></tr><tr><td>400</td><td>0.05</td><td>0.06</td><td>0.04</td><td>0.06</td></tr><tr><td>500</td><td>0.04</td><td>0.055</td><td>0.08</td><td>0.055</td></tr></tbody></table>
In particular, any strictly positive twice-differentiable density on $\mathbb{R}^d$ truncated to some $d$-dimensional interval will satisfy the assumptions for the model $p$. Any density truncated to the same $d$-dimensional interval will satisfy the assumption for the density $q$. As an example, one can consider truncated normal densities.

Given some smoothness parameter $s > 0$, radius $R > 0$ and dimension $d \in \mathbb{N} \setminus \{0\}$, the Sobolev ball $S_d^s(R)$ is defined as the function space

$$S_d^s(R) := \left\{ f \in L^1(\mathbb{R}^d) \cap L^2(\mathbb{R}^d) : \int_{\mathbb{R}^d} \|\xi\|_2^{2s} |\hat{f}(\xi)|^2 d\xi \le (2\pi)^d R^2 \right\},$$

where $\hat{f}$ denotes the Fourier transform of $f$. For $s, t, d, R, L$ all strictly positive, we define the restricted Sobolev ball $S_d^{s,t}(R, L)$ as containing all functions $f \in S_d^s(R)$ satisfying

$$\int_{\|\xi\|_2 \le t} |\hat{f}(\xi)|^2 d\xi \le \frac{1}{L} \int_{\mathbb{R}^d} |\hat{f}(\xi)|^2 d\xi. \quad (8)$$

With the Sobolev assumption $p-q \in S_d^s(R)$, the densities $p$ and $q$ can differ at any frequencies. The restricted Sobolev assumption $p-q \in S_d^{s,t}(R, L)$ does not include the case in which the densities $p$ and $q$ differ only at low frequencies due to the additional restriction in Equation (8).

**Theorem 3.5.** Suppose the assumptions listed in Appendix A.4 hold for $K_\Lambda = \{k_\lambda : \lambda \in \Lambda\}$. Let $\psi := p-q$. We show that there exists some $L > 0$ such that if $p-q \in S_d^{s,t}(R, L)$ then (i) & (ii) hold.

(i) Under the assumptions of Theorem 3.1, for the KSD test with bandwidth $\lambda := N^{-2/(4s+d)}$, the condition

$$\|\psi\|_2^2 \ge CN^{-4s/(4s+d)}$$

for some $C > 0$ guarantees control over the probability of type II error $\mathbb{P}_q(\Delta_\lambda^\Lambda(\mathbb{X}_N) = 0) \le \beta$.

(ii) Under the assumptions of Theorem 3.3, for KSDAGG with the collection

$$\Lambda := \left\{ 2^{-\ell} : \ell \in \{1, \dots, \left\lfloor \frac{2}{d} \log_2 \left( \frac{N}{\ln(\ln(N))} \right) \right\rfloor\} \right\}$$

and weights $w_\lambda := 6/\pi^2\ell^2$, we have $\mathbb{P}_q(\Delta_\lambda^\Lambda(\mathbb{X}_N) = 0) \le \beta$ provided that, for some $C > 0$,

$$\|\psi\|_2^2 \ge C \left( \frac{N}{\ln(\ln(N))} \right)^{-4s/(4s+d)}.$$

The proof of Theorem 3.5 is presented in Appendix I.4. The uniform separation rate $N^{-4s/(4s+d)}$ in Theorem 3.5 is known to be optimal in the minimax sense over (unrestricted) Sobolev balls $S_d^s(R)$ (Li and Yuan, 2019; Balasubramanian et al., 2021; Albert et al., 2022; Schrab et al., 2021). However, this rate is for the KSD test with bandwidth depending on the smoothness parameter $s > 0$ which is unknown in practice. Note that even methods which split the data to perform kernel selection are not able to select this bandwidth depending on $s$ and cannot achieve the minimax rate. The aggregation procedure is adaptive to this unknown parameter $s$: crucially, the collection of bandwidths $\Lambda$ does not depend on $s$, and so, the aggregated KSDAGG test of Theorem 3.5 (ii) can be implemented in practice, unlike the single KSD test of Theorem 3.5 (i). The price to pay for this adaptivity is only an iterated logarithmic factor in the minimax rate over Sobolev balls $S_d^s(R)$. The uniform separation rates presented in Theorem 3.5 are over the restricted Sobolev balls $S_d^{s,t}(R, L)$. Whether those rates can also be derived under less restrictive assumptions and in the more general setting of (unrestricted) Sobolev balls $S_d^s(R)$, is a challenging problem, which is left for future work.

# 4 Implementation and experiments

We consider three different experiments based on a Gamma one-dimensional distribution, a Gaussian-Bernoulli Restricted Boltzmann Machine, and a Normalizing Flow for the MNIST dataset. We compare our proposed aggregated test KSDAGG against three alternatives: the KSD test which uses the median bandwidth, a test which splits the data to select an ‘optimal’ bandwidth according to a proxy for asymptotic test power, and a test which uses extra data for bandwidth selection. The ‘extra data’ test is designed simply to provide a best-case scenario for the bandwidth selection procedure which maximises asymptotic test power, but it cannot be used in practice as it has access to extra data compared to the other tests. In order to ensure that our tests always have correct levels for all bandwidth values, dimensions and sample sizes, we use the parametric bootstrap in our experiments.
In Table 1, we observe that the four tests have correct level 0.05 for the five different sample sizes considered (the small fluctuations about the designed test level can be explained by the fact that we are averaging 200 test outputs to estimate these levels). The well-calibrated levels obtained in Table 1 demonstrate the validity of the power results presented in Figure 4.

In Figure 4, we observe that only our aggregated test KSDAGG obtains high power, that is, it is able to detect that MNIST samples are not drawn from the Normalizing Flow. The power of the other tests increases only marginally as the sample size increases. We notice that the test which uses extra data to select an ‘optimal’ bandwidth performs poorly when compared to KSDAGG. This could be explained by the fact that this test uses a proxy for the asymptotic power for bandwidth selection, and that in this high-dimensional setting, the asymptotic regime is not attained with sample sizes below 500.

<figure><img src="image_5.png" /><figcaption>Figure 4. Normalizing Flow MNIST power experiment.</figcaption></figure>

# 5. Discussion

We have introduced KSDAGG, an aggregated goodness-of-fit test based on the Kernel Stein Discrepancy. We have investigated the theoretical properties of this adaptive test. We have shown that it achieves the desired level and have provided conditions to guarantee high power by exhibiting a lower bound on its uniform separation rate. We have observed in our experiments that KSDAGG outperforms alternative state-of-the-art approaches to KSD kernel adaptation for goodness-of-fit testing.

This work covers the problem of KSD adaptivity in the goodness-of-fit framework, this complements the works of Fromont et al. (2013), Kim et al. (2020) and Schrab et al.

(2021) on MMD adaptivity for the two-sample problem and of Albert et al. (2019) and Kim et al. (2020) for HSIC adaptivity for independence testing. A potential future direction of interest could be to tackle the adaptivity problem of the KSD-based linear-time goodness-of-fit test proposed by Jitkrittum et al. (2017). In this setting, the data is split to select feature locations (and the kernel bandwidth), and the KSD test is then run using those adaptive features in data space or in the Fourier domain. A challenging problem would be to obtain those adaptive features using an aggregation procedure which avoids splitting the data.

# References

Albert, M., Laurent, B., Marrel, A., and Meynaoui, A. Adaptive test of independence based on HSIC measures. *To appear in The Annals of Statistics*, arXiv preprint arXiv:1902.06441, 2019.

Baraud, Y. Non-asymptotic minimax rates of testing in signal detection. *Bernoulli*, 1(8(5):577–606), 2002.

Carmeli, C., De Vito, E., Toigo, A., and Umanitá, V. Vector valued reproducing kernel Hilbert spaces and universality. *Analysis and Applications*, 8(01):19–61, 2010.

Chwialkowski, K., Sejdinovic, D., and Gretton, A. A wild bootstrap for degenerate kernel tests. In *Advances in neural information processing systems*, pp. 3608–3616, 2014.

Chwialkowski, K., Strathmann, H., and Gretton, A. A kernel test of goodness of fit. In *International Conference on Machine Learning*, pp. 2606–2615. PMLR, 2016.

Dinh, L., Sohl-Dickstein, J., and Bengio, S. Density estimation using real NVP. In *International Conference on Learning Representations*, 2017.

Freidling, T., Poignard, B., Climente-González, H., and Yamada, M. Post-selection inference with HSIC-Lasso. In *International Conference on Machine Learning*, pp. 3439–3448. PMLR, 2021.

Fromont, M., Laurent, B., Lerasle, M., and Reynaud-Bouret, P. Kernels based tests with non-asymptotic bootstrap approaches for two-sample problems. In *Conference on Learning Theory*, PMLR, 2012.

Fromont, M., Laurent, B., and Reynaud-Bouret, P. The two-sample problem for Poisson processes: Adaptive tests with a nonasymptotic wild bootstrap approach. *The Annals of Statistics*, 41(3):1431–1461, 2013.

Gorham, J. and Mackey, L. Measuring sample quality with kernels. In *International Conference on Machine Learning*, pp. 1292–1301. PMLR, 2017.
## 4.1 Alternative bandwidth selection approaches

Gretton et al. (2012a) proposed to use the median heuristic as kernel bandwidth, it consists in the median of the $L^2$-distances between the samples given by

$$ \lambda_{\text{med}} := \text{median}\{\|x_i - x_j\|_2 : 0 \le i < j \le N\}. $$

Gretton et al. (2012b) first proposed, for the two-sample problem using a linear-time MMD estimator, to split the data and to use half of it to select an ‘optimal’ bandwidth which maximises a proxy for asymptotic power. This procedure was extended to quadratic-time estimators and to the goodness-of-fit framework by Jitkrittum et al. (2017), Sutherland et al. (2017) and Liu et al. (2020). These strategies rely on the asymptotic normality of the test statistic under $\mathcal{H}_a$. In our setting, the asymptotic power proxy to maximise is the ratio

$$ \widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) / \widehat{\sigma}_{\mathcal{H}_a} \qquad (9) $$

where $\widehat{\sigma}_{\mathcal{H}_a}^2$ is a closed-form regularised positive estimator of the asymptotic variance of $\widehat{\text{KSD}}_{p,k}^2$ under $\mathcal{H}_a$ (Liu et al., 2020, Equation 5). In our experiments, we also consider a test which has access to $N$ extra samples drawn from $q$ to select an ‘optimal’ bandwidth to run the KSD test on the original $N$ samples $\mathbb{X}_N$. This test is interesting to compare to because it uses an ‘optimal’ bandwidth without being detrimental to power (as it uses all $N$ samples $\mathbb{X}_N$ to run the test with the selected bandwidth).

## 4.2 Experimental details

Inspired from Theorem 3.5, in each experiments, we use KSDAGG with different collections of bandwidths of the form $\Lambda(\ell_-, \ell_+) := \{2^i \lambda_{\text{med}} : i = \ell_-, \dots, \ell_+\}$ for the median bandwidth $\lambda_{\text{med}}$ and integers $\ell_- < \ell_+$ with uniform weights $w_\lambda := 1/(\ell_+ - \ell_- + 1)$. For the tests which split the data, we select the bandwidth, out of the collection $\Lambda(\ell_-, \ell_+)$, which maximises the power proxy discussed in Section 4.1. All our experiments are run with level $\alpha = 0.05$ using the IMQ kernel defined in Equation (7) with parameter $\beta_k = 0.5$. We use a parametric bootstrap with $B_1 = B_2 = 500$ bootstrapped KSD values to compute the adjusted test thresholds, and $B_3 = 50$ steps of bisection method to estimate the correction $u_\alpha$ in Equation (6). For our aggregated test, we have also designed another collection which is parameter-free (no varying parameters across experiments) and performs as well as the previous variant whose parameters are chosen appropriately for each experiments. The collection is obtained using the maximal inter-sample distance normalised by the dimension (see details in Appendix B). We denote the resulting robust test by KSDAGG\* for which we use either a wild or parametric bootstrap with $B_1 = B_2 = 2000$ and $B_3 = 50$, this is the variant we recommend using in practice. To estimate the probability of rejecting the null hypothesis, we average the test outputs across 200 repetitions. All experiments have been run on an AMD Ryzen Threadripper 3960X 24 Cores 128Gb RAM CPU at 3.8GHz, the runtime is of the order of a couple of hours (significant speedup can be obtained by using parallel computing). We have used the implementation of Jitkrittum et al. (2017) to sample from a Gaussian-Bernoulli Restricted Boltzmann Machine, and Phillip Lippe’s implementation of MNIST Normalizing Flows, both under the MIT license.

## 4.3 Gamma distribution

For our first experiment, we consider a one-dimensional Gamma distribution with shape parameter 5 and scale parameter 5 as the model $p$. For $q$, we draw 500 samples from a Gamma distribution with the same scale parameter 5 and with a shifted shape parameter $5 + s$ for $s \in \{0, 0.1, 0.2, 0.3, 0.4\}$. We consider the collection of bandwidths $\Lambda(0, 10)$.

The results we obtained are presented in Figure 1a. We observe that all tests have the prescribed level 0.05 under the null hypothesis, which corresponds to the case $s = 0$. As the shift parameter $s$ increases, the two densities $p$ and $q$ become more different and rejection of the null becomes an easier task, thus the test power increases. Our aggregated test KSDAGG achieves the same power as the KSD test with an ‘optimal’ bandwidth selected using extra data by maximizing the proxy for asymptotic power. The median test obtains only slightly lower power, this closeness in power can be explained by the fact that this one-dimensional problem is a simple one. All four tests KSDAGG $\Lambda(0, 10)$, KSDAGG\*, KSD median, and KSD split, all achieve very similar power. We note that the normal splitting test has significantly lower power: this is because, even though it uses an ‘optimal’ bandwidth, it is then run on only half the data, which results in a loss of power.
Grathwohl, W., Wang, K.-C., Jacobsen, J.-H., Duvenaud, D., and Zemel, R. Learning the Stein discrepancy for training and evaluating energy-based models without sampling. In *International Conference on Machine Learning*, pp. 3732–3747. PMLR, 2020.

Gretton, A., Bousquet, O., Smola, A., and Schölkopf, B. Measuring statistical dependence with Hilbert-Schmidt norms. In *International Conference on Algorithmic Learning Theory*. Springer, 2005.

Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., and Smola, A. A kernel two-sample test. *Journal of Machine Learning Research*, 13:723–773, 2012a.

Gretton, A., Sejdinovic, D., Strathmann, H., Balakrishnan, S., Pontil, M., Fukumizu, K., and Sriperumbudur, B. K. Optimal kernel choice for large-scale two-sample tests. In *Advances in Neural Information Processing Systems*, volume 1, pp. 1205–1213, 2012b.

Hoeffding, W. A class of statistics with asymptotically normal distribution. In *Breakthroughs in Statistics*, pp. 308–334. Springer, 1992.

Jitkrittum, W., Xu, W., Szabó, Z., Fukumizu, K., and Gretton, A. A linear-time kernel goodness-of-fit test. In *Advances in Neural Information Processing Systems*, pp. 262–271, 2017.

Key, O., Fernandez, T., Gretton, A., and Briol, F.-X. Composite goodness-of-fit tests with kernels. *arXiv preprint arXiv:2111.10275*, 2021.

Kim, I., Balakrishnan, S., and Wasserman, L. Minimax optimality of permutation tests. *To appear in The Annals of Statistics*, arXiv preprint arXiv:2003.13208, 2020.

Kingma, D. P. and Dhariwal, P. Glow: Generative flow with invertible 1x1 convolutions. In *Advances in Neural Information Processing Systems*, pp. 10236–10245, 2018.

LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. Gradient-based learning applied to document recognition. *Proceedings of the IEEE*, 86(11):2278–2324, 1998.

LeCun, Y., Cortes, C., and Burges, C. MNIST handwritten digit database. AT&T Labs, 2010. URL http://yann.lecun.com/exdb/mnist/.

Lee, J. *U-statistics: Theory and Practice*. Citeseer, 1990.

Leucht, A. and Neumann, M. H. Dependent wild bootstrap for degenerate U- and V-statistics. *Journal of Multivariate Analysis*, 117:257–280, 2013.

Lim, J. N., Yamada, M., Schölkopf, B., and Jitkrittum, W. Kernel stein tests for multiple model comparison. In *Advances in Neural Information Processing Systems*, pp. 2240–2250, 2019.

Lim, J. N., Yamada, M., Jitkrittum, W., Terada, Y., Matsui, S., and Shimodaira, H. More powerful selective kernel tests for feature selection. In *International Conference on Artificial Intelligence and Statistics*, pp. 820–830. PMLR, 2020.

Liu, F., Xu, W., Lu, J., Zhang, G., Gretton, A., and Sutherland, D. J. Learning deep kernels for non-parametric two-sample tests. In *International Conference on Machine Learning*, 2020.

Liu, Q., Lee, J., and Jordan, M. A kernelized Stein discrepancy for goodness-of-fit tests. In *International Conference on Machine Learning*, pp. 276–284. PMLR, 2016.

Romano, J. P. and Wolf, M. Exact and approximate step-down methods for multiple hypothesis testing. *Journal of the American Statistical Association*, 100(469):94–108, 2005.

Schrab, A., Kim, I., Albert, M., Laurent, B., Guedj, B., and Gretton, A. MMD aggregated two-sample test. *arXiv preprint arXiv:2110.15073*, 2021.

Shao, X. The dependent wild bootstrap. *Journal of the American Statistical Association*, 105(489):218–235, 2010.

Stute, W., Manteiga, W. G., and Quindimil, M. P. Bootstrap based goodness-of-fit-tests. *Metrika*, 40(1):243–256, 1993.

Sutherland, D. J., Tung, H.-Y., Strathmann, H., De, S., Ramdas, A., Smola, A., and Gretton, A. Generative models and model criticism via optimized maximum mean discrepancy. In *International Conference on Learning Representations*, 2017.

Yamada, M., Wu, D., Tsai, Y. H., Ohta, H., Salakhutdinov, R., Takeuchi, I., and Fukumizu, K. Post selection inference with incomplete maximum mean discrepancy estimator. In *International Conference on Learning Representations*, 2019.
<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Gamma distribution. (b) Gaussian-Bernoulli Restricted Boltzmann Machine.</figcaption></figure>

## 4.4 Gaussian-Bernoulli Restricted Boltzmann Machine

As first considered by Liu et al. (2016) for goodness-of-fit testing using the KSD, we consider a Gaussian-Bernoulli Restricted Boltzmann Machine. It is a graphical model with a binary hidden variable $h \in \{-1, 1\}^{d_h}$ and a continuous observable variable $x \in \mathbb{R}^d$. Those variables have joint density

$$p(x, h) = \frac{1}{Z} \exp \left( \frac{1}{2} x^\top B h + b^\top x + c^\top h - \frac{1}{2} \|x\|_2^2 \right)$$

where $Z$ is an unknown normalizing constant. By marginalising over $h$, we obtain the density $p$ of $x$

$$p(x) = \sum_{h \in \{-1, 1\}^{d_h}} p(x, h).$$

We can sample from it using a Gibbs sampler with 2000 burn-in iterations. We use the dimensions $d = 50$ and $d_h = 40$ as considered by Jitkrittum et al. (2017) and Grathwohl et al. (2020). Even though computing $p$ is intractable for large dimension $d_h$, the score function admits a convenient closed-form expression

$$\nabla \log p(x) = b - x + B \frac{\exp(2(B^\top x + c)) - 1}{\exp(2(B^\top x + c)) + 1}.$$

We draw the components of $b$ and $c$ from Gaussian standard distributions and sample Rademacher variables taking values in $\{-1, 1\}$ for the elements of $B$ for the model $p$. We draw 1000 samples from a distribution $q$ which is constructed in a similar way as $p$ but with the difference that some Gaussian noise $N(0, \sigma)$ is injected into each of the elements of $B$. For the standard deviations of the perturbations, we consider $\sigma \in \{0, 0.01, 0.02, 0.03\}$. We use the collection $\Lambda(-20, 0)$ (different from Section 4.3), for KSDAGG* we use a wild bootstrap, the results are provided in Figure 1b.

Again, we observe that our aggregated tests KSDAGG $\Lambda(0, 20)$ and KSDAGG* match the power obtained by the test which uses extra data to select an 'optimal' bandwidth. This means that, in this experiment, the aggregated tests obtain the same power as the 'best' single test. The difference with the median heuristic test is significant in this experiment, and the splitting test obtains the lowest power of the four tests. Again, all tests have well-calibrated levels (case $\sigma = 0$) and increasing the noise level $\sigma$ results in more power for all the tests.

## 4.5 MNIST Normalizing Flow

Finally, we consider a high-dimensional problem working with images in dimensions $28^2 = 784$. We consider a multi-scale Normalizing Flow (Dinh et al., 2017; Kingma and Dhariwal, 2018) which has been trained on the MNIST dataset (LeCun et al., 1998, 2010), it is a generative model which has a
<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Digits from the MNIST dataset. (b, c) Digits sampled from the Normalizing Flow.</figcaption></figure>

probability density *p*. As observed in Figure 2, some samples produced by the model can look exactly like MNIST digits, while others do not resemble digits. This Normalizing Flow has been trained to 'ideally' produce samples from the MNIST dataset. We are interested in whether or not we can detect the difference in densities. Given some images of digits, are we able to tell if those were sampled from the Normalizing Flow model? We consider the case where the samples from *q* are drawn from the true MNIST dataset (power experiment), and the case where the images from *q* are sampled from the Normalizing Flow model (level experiment, confirming performance for the power experiment). The experiments are run with the collection of bandwidths $\Lambda(-20, 0)$. The results are displayed in Figure 3a and Figure 3b. We use a parametric bootstrap for KSDAGG*.

In Figure 3a, we observe that the four tests have correct levels (around 0.05) for the five different sample sizes considered (the small fluctuations about the designed test level can be explained by the fact that we are averaging 200 test outputs to estimate these levels). The well-calibrated levels obtained in Figure 3a demonstrate the validity of the power results presented in Figure 3b.

As seen in Figure 3b, only our aggregated tests KSDAGG and KSDAGG* obtain high power; they are able to detect that MNIST samples are not drawn from the Normalizing Flow. The power of the other tests increases only marginally as the sample size increases. The test which uses extra data to select an 'optimal' bandwidth performs poorly when compared to the aggregated tests. This could be explained by the fact that this test selects the bandwidth using a proxy for the asymptotic power, and that in this high-dimensional setting, the asymptotic regime is not attained with sample sizes below 500. See Appendices E and F for details regarding selected bandwidths and for reported runtimes.

<figure><img src="image_4.png" /><figcaption>Figure 3: Normalizing Flow MNIST. (a) Level experiment. (b) Power experiment.</figcaption></figure>

# 5 Acknowledgements

Antonin Schrab acknowledges support from the U.K. Research and Innovation (EP/S021566/1). Benjamin Guedj acknowledges partial support by the U.S. Army Research Laboratory and the U.S. Army Research Office, and by the U.K. Ministry of Defence and the U.K. Engineering and Physical Sciences Research Council (EP/R013616/1), and by the French National Agency for Research (ANR-18-CE40-0016-01 & ANR-18-CE23-0015-02). Arthur Gretton acknowledges support from the Gatsby Charitable Foundation. We also thank Bharath Sriperumbudur for helpful discussions.
References

Albert, M., Laurent, B., Marrel, A., and Meynaoui, A. (2022). Adaptive test of independence based on HSIC measures. *The Annals of Statistics*, 50(2):858–879.

Anastasiou, A., Barp, A., Briol, F.-X., Ebner, B., Gaunt, R. E., Ghaderinezhad, F., Gorham, J., Gretton, A., Ley, C., Liu, Q., et al. (2021). Stein's method meets statistics: A review of some recent developments. *arXiv preprint arXiv:2105.03481*.

Aronszajn, N. (1950). Theory of reproducing kernels. *Transactions of the American Mathematical Society*, 68(3):337–404.

Balasubramanian, K., Li, T., and Yuan, M. (2021). On the optimality of kernel-embedding based goodness-of-fit tests. *Journal of Machine Learning Research*, 22(1).

Baraud, Y. (2002). Non-asymptotic minimax rates of testing in signal detection. *Bernoulli*, 8(85):577–606.

Barp, A., Briol, F.-X., Duncan, A., Girolami, M., and Mackey, L. (2019). Minimum stein discrepancy estimators. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, *Advances in Neural Information Processing Systems*, volume 32. Curran Associates, Inc.

Carmeli, C., De Vito, E., Toigo, A., and Umanitá, V. (2010). Vector valued reproducing kernel Hilbert spaces and universality. *Analysis and Applications*, 8(01):19–61.

Chebyshev, P. L. (1899). *Oeuvres. Commissionaires de l'Académie Impériale des Sciences*, 1.

Chwialkowski, K., Sejdinovic, D., and Gretton, A. (2014). A wild bootstrap for degenerate kernel tests. In *Advances in neural information processing systems*, pages 3608–3616.

Chwialkowski, K., Strathmann, H., and Gretton, A. (2016). A kernel test of goodness of fit. In *International Conference on Machine Learning*, pages 2606–2615. PMLR.

de la Peña, V. H. and Giné, E. (1999). *Decoupling: From Dependence to Independence*. Springer Science & Business Media.

Dinh, L., Sohl-Dickstein, J., and Bengio, S. (2017). Density estimation using real NVP. In *International Conference on Learning Representations*.

Dvoretzky, A., Kiefer, J., and Wolfowitz, J. (1956). Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator. *The Annals of Mathematical Statistics*, pages 642–669.

Fernandez, T. and Gretton, A. (2019). A maximum-mean-discrepancy goodness-of-fit test for censored data. In *Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics*, volume 89 of *Proceedings of Machine Learning Research*, pages 2966–2975. PMLR.

Fernandez, T. and Rivera, N. (2022). A general framework for the analysis of kernel-based tests. *arXiv preprint arXiv:2209.00124*.

Fernandez, T., Rivera, N., Xu, W., and Gretton, A. (2020). Kernelized Stein discrepancy tests of goodness-of-fit for time-to-event data. In *Proceedings of the 37th International Conference on Machine Learning*, volume 119 of *Proceedings of Machine Learning Research*, pages 3112–3122. PMLR.

Freidling, T., Poignard, B., Climente-González, H., and Yamada, M. (2021). Post-selection inference with HSIC-Lasso. In *International Conference on Machine Learning*, pages 3439–3448. PMLR.

Fromont, M., Laurent, B., Lerasle, M., and Reynaud-Bouret, P. (2012). Kernels based tests with non-asymptotic bootstrap approaches for two-sample problems. In *Conference on Learning Theory*, PMLR.
Fromont, M., Laurent, B., and Reynaud-Bouret, P. (2013). The two-sample problem for Poisson processes: Adaptive tests with a nonasymptotic wild bootstrap approach. *The Annals of Statistics*, 41(3):1431–1461.

Gorham, J. and Mackey, L. (2017). Measuring sample quality with kernels. In *International Conference on Machine Learning*, pages 1292–1301. PMLR.

Gorham, J. and Mackey, L. W. (2015). Measuring sample quality with Stein's method. In Cortes, C., Lawrence, N. D., Lee, D. D., Sugiyama, M., and Garnett, R., editors, *Advances in Neural Information Processing Systems 28: Annual Conference on Neural Information Processing Systems 2015*, December 7-12, 2015, Montreal, Quebec, Canada, pages 226–234.

Grathwohl, W., Wang, K.-C., Jacobsen, J.-H., Duvenaud, D., and Zemel, R. (2020). Learning the Stein discrepancy for training and evaluating energy-based models without sampling. In *International Conference on Machine Learning*, pages 3732–3747. PMLR.

Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., and Smola, A. (2012a). A kernel two-sample test. *Journal of Machine Learning Research*, 13:723–773.

Gretton, A., Bousquet, O., Smola, A., and Schölkopf, B. (2005). Measuring statistical dependence with Hilbert-Schmidt norms. In *International Conference on Algorithmic Learning Theory*. Springer.

Gretton, A., Sejdinovic, D., Strathmann, H., Balakrishnan, S., Pontil, M., Fukumizu, K., and Sriperumbudur, B. K. (2012b). Optimal kernel choice for large-scale two-sample tests. In *Advances in Neural Information Processing Systems*, volume 1, pages 1205–1213.

Hoeffding, W. (1992). A class of statistics with asymptotically normal distribution. In *Breakthroughs in Statistics*, pages 308–334. Springer.

Huggins, J. and Mackey, L. (2018). Random feature Stein discrepancies. *Advances in Neural Information Processing Systems*, 31.

Ingster, Y. I. (1987). Minimax testing of nonparametric hypotheses on a distribution density in the $L_p$ metrics. *Theory of Probability & its Applications*, 31(2):333–337.

Ingster, Y. I. (1989). An asymptotically minimax test of the hypothesis of independence. *Journal of Soviet Mathematics*, 1(44):466–476.

Ingster, Y. I. (1993a). Asymptotically minimax hypothesis testing for nonparametric alternatives. *Journal of Soviet Mathematics*, 1(44):466–476.

Ingster, Y. I. (1993b). Minimax testing of the hypothesis of independence for ellipsoids in $l_p$. *Zapiski Nauchnykh Seminarov POMI*, 1(207:77–97).

Jitkrittum, W., Kanagawa, H., and Schölkopf, B. (2020). Testing goodness of fit of conditional density models with kernels. In Peters, J. and Sontag, D., editors, *Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI)*, volume 124 of *Proceedings of Machine Learning Research*, pages 221–230. PMLR.

Jitkrittum, W., Xu, W., Szabó, Z., Fukumizu, K., and Gretton, A. (2017). A linear-time kernel goodness-of-fit test. In *Advances in Neural Information Processing Systems*, pages 262–271.

Kanagawa, H., Jitkrittum, W., Mackey, L., Fukumizu, K., and Gretton, A. (2019). A kernel Stein test for comparing latent variable models. *arXiv:1907.00586 [cs, stat]*. arXiv: 1907.00586.

Key, O., Fernandez, T., Gretton, A., and Briol, F.-X. (2021). Composite goodness-of-fit tests with kernels. *arXiv preprint arXiv:2111.10275*.

Kim, I., Balakrishnan, S., and Wasserman, L. (2022). Minimax optimality of permutation tests. *The Annals of Statistics*, 50(1):225–251.

Kingma, D. P. and Dhariwal, P. (2018). Glow: Generative flow with invertible 1x1 convolutions. In *Advances in Neural Information Processing Systems*, pages 10236–10245.
Korba, A., Aubin-Frankowski, P., Majewski, S., and Ablin, P. (2021). Kernel Stein discrepancy descent. In Meila, M. and Zhang, T., editors, *Proceedings of the 38th International Conference on Machine Learning*, ICML 2021, 18-24 July 2021, volume 139 of *Proceedings of Machine Learning Research*, pages 5719–5730. PMLR.

Kübler, J. M., Jitkrittum, W., Schölkopf, B., and Muandet, K. (2020). Learning kernel tests without data splitting. In *Advances in Neural Information Processing Systems 33*, pages 6245–6255. Curran Associates, Inc.

Kübler, J. M., Jitkrittum, W., Schölkopf, B., and Muandet, K. (2022). A witness two-sample test. In *International Conference on Artificial Intelligence and Statistics*, pages 1403–1419. PMLR.

LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998). Gradient-based learning applied to document recognition. *Proceedings of the IEEE*, 86(11):2278–2324.

LeCun, Y., Cortes, C., and Burges, C. (2010). MNIST handwritten digit database. AT&T Labs.

Leucht, A. (2012). Degenerate $U$- and $V$-statistics under weak dependence: Asymptotic theory and bootstrap consistency. *Bernoulli*, 18(2):552–585.

Leucht, A. and Neumann, M. H. (2013). Dependent wild bootstrap for degenerate $U$- and $V$-statistics. *Journal of Multivariate Analysis*, 117:257–280.

Li, T. and Yuan, M. (2019). On the optimality of gaussian kernel based nonparametric tests against smooth alternatives. *arXiv preprint arXiv:1909.03302*.

Lim, J. N., Yamada, M., Jitkrittum, W., Terada, Y., Matsui, S., and Shimodaira, H. (2020). More powerful selective kernel tests for feature selection. In *International Conference on Artificial Intelligence and Statistics*, pages 820–830. PMLR.

Lim, J. N., Yamada, M., Schölkopf, B., and Jitkrittum, W. (2019). Kernel Stein tests for multiple model comparison. In *Advances in Neural Information Processing Systems*, pages 2240–2250.

Liu, F., Xu, W., Lu, J., Zhang, G., Gretton, A., and Sutherland, D. J. (2020). Learning deep kernels for non-parametric two-sample tests. In *International Conference on Machine Learning*.

Liu, Q., Lee, J., and Jordan, M. (2016). A kernelized Stein discrepancy for goodness-of-fit tests. In *International Conference on Machine Learning*, pages 276–284. PMLR.

Massart, P. (1990). The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality. *The Annals of Probability*, 18(3):1269–1283.

Oates, C. J., Cockayne, J., Briol, F.-X., and Girolami, M. (2019). Convergence rates for a class of estimators based on Stein's method. *Bernoulli*, 25(2):1141–1159.

Oates, C. J. and Girolami, M. A. (2016). Control functionals for quasi-Monte Carlo integration. In Gretton, A. and Robert, C. C., editors, *Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, AISTATS 2016, Cadiz, Spain, May 9-11, 2016*, volume 51 of *JMLR Workshop and Conference Proceedings*, pages 56–65. JMLR.org.

Romano, J. P. and Wolf, M. (2005). Exact and approximate stepdown methods for multiple hypothesis testing. *Journal of the American Statistical Association*, 100(469):94–108.

Schrab, A., Kim, I., Albert, M., Laurent, B., Guedj, B., and Gretton, A. (2021). MMD Aggregated two-sample test. *arXiv preprint arXiv:2110.15073*.

Schrab, A., Kim, I., Guedj, B., and Gretton, A. (2022). Efficient aggregated kernel tests using incomplete $U$-statistics. In *Advances in Neural Information Processing Systems 35: Annual Conference on Neural Information Processing Systems 2022, NeurIPS 2022*.

Shao, X. (2010). The dependent wild bootstrap. *Journal of the American Statistical Association*, 105(489):218–235.

Sriperumbudur, B. K., Fukumizu, K., and Lanckriet, G. R. (2011). Universality, characteristic kernels and RKHS embedding of measures. *Journal of Machine Learning Research*, 12(7).
Stein, C. (1972). 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 6, pages 583–603. University of California Press.

Stein, C., Diaconis, P., Holmes, S., and Reinert, G. (2004). Use of exchangeable pairs in the analysis of simulations. *Lecture Notes-Monograph Series*, pages 1–26.

Stute, W., Manteiga, W. G., and Quindimil, M. P. (1993). Bootstrap based goodness-of-fit-tests. *Metrika*, 40(1):243–256.

Sutherland, D. J., Tung, H.-Y., Strathmann, H., De, S., Ramdas, A., Smola, A., and Gretton, A. (2017). Generative models and model criticism via optimized maximum mean discrepancy. In *International Conference on Learning Representations*.

Tolstikhin, I., Sriperumbudur, B. K., and Schölkopf, B. (2016). Minimax estimation of maximum mean discrepancy with radial kernels. *Advances in Neural Information Processing Systems*, 29.

Weckbecker, M., Xu, W., and Reinert, G. (2022). On RKHS choices for assessing graph generators via kernel Stein statistics. *arXiv preprint arXiv:2210.05746*.

Wenliang, L. K. and Kanagawa, H. (2020). Blindness of score-based methods to isolated components and mixing proportions. *arXiv preprint arXiv:2008.10087*.

Wynne, G. and Duncan, A. B. (2022). A kernel two-sample test for functional data. *Journal of Machine Learning Research*, 23(73):1–51.

Wynne, G., Kasprzak, M., and Duncan, A. B. (2022). A spectral representation of kernel stein discrepancy with application to goodness-of-fit tests for measures on infinite dimensional hilbert spaces. *arXiv preprint arXiv:2206.04552*.

Xu, W. (2021). Generalised kernel Stein discrepancy (GKSD): A unifying approach for non-parametric goodness-of-fit testing. *arXiv preprint arXiv:2106.12105*.

Xu, W. (2022). Standardisation-function kernel Stein discrepancy: A unifying view on kernel Stein discrepancy tests for goodness-of-fit. In Camps-Valls, G., Ruiz, F. J. R., and Valera, I., editors, *Proceedings of The 25th International Conference on Artificial Intelligence and Statistics*, volume 151 of *Proceedings of Machine Learning Research*, pages 1575–1597. PMLR.

Xu, W. and Matsuda, T. (2020). A Stein goodness-of-fit test for directional distributions. In Chiappa, S. and Calandra, R., editors, *Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics*, volume 108 of *Proceedings of Machine Learning Research*, pages 320–330. PMLR.

Xu, W. and Matsuda, T. (2021). Interpretable Stein goodness-of-fit tests on Riemannian manifold. In Meila, M. and Zhang, T., editors, *Proceedings of the 38th International Conference on Machine Learning*, volume 139 of *Proceedings of Machine Learning Research*, pages 11502–11513. PMLR.

Xu, W. and Reinert, G. (2021). A Stein goodness-of-test for exponential random graph models. In Banerjee, A. and Fukumizu, K., editors, *Proceedings of The 24th International Conference on Artificial Intelligence and Statistics*, volume 130 of *Proceedings of Machine Learning Research*, pages 415–423. PMLR.

Xu, W. and Reinert, G. (2022a). AgraSSt: Approximate graph Stein statistics for interpretable assessment of implicit graph generators. *arXiv preprint arXiv:2203.03673*.

Xu, W. and Reinert, G. (2022b). A kernelised Stein statistic for assessing implicit generative models. *arXiv preprint arXiv:2206.00149*.

Yamada, M., Wu, D., Tsai, Y. H., Ohta, H., Salakhutdinov, R., Takeuchi, I., and Fukumizu, K. (2019). Post selection inference with incomplete maximum mean discrepancy estimator. In *International Conference on Learning Representations*.

Zhang, M., Key, O., Hayes, P., Barber, D., Paige, B., and Briol, F.-X. (2022). Towards healing the blindness of score matching. *arXiv preprint arXiv:2209.07396*.
# Supplementary material for ‘KSD Aggregated Goodness-of-fit Test’

## A Assumptions

### A.1 Assumptions common to all theoretical results

* kernel $k$ is $C_0$-universal (Carmeli et al., 2010, Definition 4.1)

* KSD$_{p,k}^2(q) = \mathbb{E}_{q,q}[h_{p,k}(X, Y)] < \infty

* $C_k \ge 1$, where $C_k := \mathbb{E}_{q,q}[h_{p,k}(X, Y)^2] < \infty$ as defined in Equation (2)

* $\|q\|_{\infty} \le M$ for some $M > 0$

* $\mathbb{E}_q[\left\|\nabla\left(\log\frac{p(X)}{q(X)}\right)\right\|_2^2] < \infty$

* $\alpha \in (0, e^{-1})$

* $\beta \in (0, 1)$

### A.2 Assumptions for Theorem 3.1

* Assumptions of Appendix A.1

* $B_1 \ge 3(\log(8/\beta) + \alpha(1 - \alpha))/\alpha^2$

* assumption for parametric bootstrap: $n/\sqrt{C_k} \ge \ln(1/\alpha)$

* assumption for wild bootstrap: Lipschitz continuity of $h_{p,k}$

* constant $C$ depends on $M$ and $d$

### A.3 Assumptions for Theorem 3.3 and Corollary 3.4

* Assumptions of Appendix A.1

* $B_1 \ge 12(\max_{k \in \mathcal{K}} w_k^{-2})(\log(8/\beta) + \alpha(1 - \alpha))/\alpha^2$

* $B_2 \ge 8\log(2/\beta)/\alpha^2$

* $B_3 \ge \log_2(4 \min_{k \in \mathcal{K}} w_k^{-1}/\alpha)$

* positive weights $(w_k)_{k \in \mathcal{K}}$ satisfying $\sum_{k \in \mathcal{K}} w_k \le 1$

* assumption for parametric bootstrap: $n/\sqrt{C_k} \ge \ln(1/\alpha w_k)$ for $k \in \mathcal{K}$

* assumption for wild bootstrap: Lipschitz continuity of $h_{p,k}$

* constant $C$ depends on $M$ and $d$

### A.4 Assumptions for Theorem 3.5

* Assumptions of Appendix A.1

* $B_1 \ge 12(\max_{\lambda \in \Lambda} w_{k_\lambda}^{-2})(\log(8/\beta) + \alpha(1 - \alpha))/\alpha^2$

* $B_2 \ge 8\log(2/\beta)/\alpha^2$

* $B_3 \ge \log_2(4 \min_{\lambda \in \Lambda} w_{k_\lambda}^{-1}/\alpha)$

* assumption for parametric bootstrap: $n/\sqrt{C_{k_\lambda}} \ge \ln(1/\alpha w_{k_\lambda})$ for $\lambda \in \Lambda$

* assumption for wild bootstrap: Lipschitz continuity of $h_{p,k_\lambda}$

* constant $C$ is independent of the sample size $N$

* $\|T_{h_{p,k_\lambda}}\psi\|_2 \le \tilde{C}_0 \|T_{\bar{k}_\lambda}\psi\|_2$ for some $\tilde{C}_0 > 0$ with $\bar{k}_\lambda$ as in Equation (16)
## B Collection of bandwidths for KSDAGG*

In this section, we explain how the collection for the robust test KSDAGG* is constructed. First, we compute the maximal inter-sample distance

$$ 
\begin{aligned}
\lambda_M &:= \max \{ \|x_i - x_j\|_2 : 0 \le i < j \le N \}, \\
\lambda_{\max} &:= \max \{ \lambda_M, 2 \}.
\end{aligned}
$$

The collection of B bandwidths is then defined as

$$ \{ d^{-1} \lambda_{\max}^{(i-1)/(B-1)} : i = 1, \dots, B \} $$

where d is the dimension of the samples. This collection is a discretisation of the interval $[d^{-1}, d^{-1}\lambda_{\max}]$. In the limit as the number of bandwidths B goes to infinity, the collection becomes the continuous interval. In practice, we use $B = 10$ bandwidths and we observed that increasing B in all three experiments ($B = 100, 200, \dots$) does not affect the power, that is, by using 10 bandwidths we obtain the same power as if we were to consider the aggregated test with the continuous collection $[d^{-1}, d^{-1}\lambda_{\max}]$ of bandwidths.

## C Background on Kernel Stein Discrepancy

**Background.** Stein's methods (Stein, 1972) have been widely used in the machine learning and statistics communities. At the heart of this field, for the Langevin Stein operator $A_p$ defined as

$$ (A_p f)(x) := f(x)^{\top} \nabla \log p(x) + \mathrm{Tr}(\nabla f(x)), $$

is the fact that smooth densities p and q are equal if and only if

$$ \mathbb{E}_q[(A_p f)(x)] = 0 $$

for smooth functions f vanishing at the boundaries. This is known as Stein's identity (Stein, 1972; Stein et al., 2004), and also holds more generally for other Stein operators. Gorham and Mackey (2015) utilised this identity to define Stein discrepancies as

$$ \sup_{f \in \mathcal{F}} \mathbb{E}_q [f(x)^{\top} \nabla \log p(x) + \mathrm{Tr}(\nabla f(x))] $$

for some space of smooth functions $\mathcal{F}$ satisfying Stein's identity. As proposed by Oates and Girolami (2016), choosing $\mathcal{F}$ to be the unit ball in a Reproducing Kernel Hilbert Space (Aronszajn, 1950, RKHS) denoted $\mathcal{H}$, that is $\mathcal{F} := \{f : \|f\|_{\mathcal{H}} \le 1\}$, we obtain the Kernel Stein Discrepancy (KSD; Liu et al., 2016; Chwialkowski et al., 2016). The KSD can also be expressed in terms of the Stein kernel as presented in Section 2.

**Goodness-of-fit related work.** Stein methods have been used to construct goodness-of-fit tests in various types of data: directional data (Xu and Matsuda, 2020), time-to-event data (Fernandez et al., 2020), data on Riemannian manifolds (Xu and Matsuda, 2021), conditional data (Jitkrittum et al., 2020), graph data (Weckbecker et al., 2022; Xu and Reinert, 2021, 2022a), functional data (Wynne et al., 2022; Wynne and Duncan, 2022), and generative data (Xu and Reinert, 2022b). We also point out the works of Fernandez and Gretton (2019) on an MMD-based goodness-of-fit test for censored data, of Key et al. (2021) on composite goodness-of-fit tests using either the MMD or KSD, of Jitkrittum et al. (2017) on a linear-time KSD goodness-of-fit test, of Gorham and Mackey (2017) on measuring sample quality with kernels and on KSD dominating weak convergence for some kernels, of Huggins and Mackey (2018) on random feature Stein discrepancies, of Korba et al. (2021) on the KSD Wasserstein gradient flow, of Oates et al. (2019) on convergence rates for a class of estimators based on Stein's method, and of Barp et al. (2019) on minimum Stein discrepancy estimators. Lim et al. (2019) and Kanagawa et al. (2019) use KSD tests for multiple model comparison and for comparing latent variable models, respectively. It is worth noting those score-based methods suffer from blindness to isolated components and mixing proportions (Wenliang and Kanagawa, 2020; Zhang et al., 2022). Xu (2021, 2022) present unified approaches for goodness-of-fit testing, and Fernandez and Rivera (2022) introduce a general framework for analysing kernel-based tests working directly with random functionals on Hilbert spaces. Finally, Anastasiou et al. (2021) provide a review of recent developments based on Stein's methods.
## D Asymptotic level for KSD test with wild bootstrap

The proof of asymptotic level of the KSD test using a wild bootstrap was proved by Chwialkowski et al. (2016, Proposition 3.2) relying on the result of Leucht (2012, Theorem 2.1) and of Chwialkowski et al. (2014, Lemma 5). Their results for V-statistics can be extended to the use of U-statistics. They proved this result for the more general non-i.i.d. case, while we need only the result for the simpler i.i.d. case. They proved that the difference between true quantiles and the wild bootstrap quantiles converges to zero in probability under the null hypothesis with the following dependence on N:

$$ \sup_{t \in \mathbb{R}} \left| \mathbb{P}(NB_N > t | \mathcal{X}_N) - \mathbb{P}(NK_N > t) \right| $$

converges to 0 in probability, where $K_N$ is the KSD estimator of Equation (1), and $B_N$ is the wild bootstrap KSD of Equation (5). Lipschitz continuity of the Stein kernel $h_{p,k}$ (defined in Section 2) is required for the result to hold.

## E Runtimes comparison of all KSD tests considered

The time complexity of KSDAGG is $O(|K|(B_1 + B_2)N^2)$ as provided in Algorithm 1. It grows linearly with the number of kernels $|K|$, quadratically with the sample size $N$, and linearly with the number of bootstrap samples $B_1 + B_2$.

<table><caption>Table 1: Runtimes (averaged over 10 runs and reported in seconds) for the Normalizing Flow MNIST experiment presented in Section 4.5 using a wild bootstrap.</caption><thead><tr><th>SAMPLE SIZE</th><th>KSDAGG</th><th>KSD MEDIAN</th><th>KSD SPLIT</th><th>KSD SPLIT EXTRA DATA</th></tr></thead><tbody><tr><td>100</td><td>0.037</td><td>0.005</td><td>0.022</td><td>0.023</td></tr><tr><td>200</td><td>0.084</td><td>0.010</td><td>0.064</td><td>0.070</td></tr><tr><td>300</td><td>0.162</td><td>0.020</td><td>0.132</td><td>0.145</td></tr><tr><td>400</td><td>0.276</td><td>0.034</td><td>0.230</td><td>0.253</td></tr><tr><td>500</td><td>0.421</td><td>0.051</td><td>0.359</td><td>0.395</td></tr></tbody></table>

We report in Table 1 the runtimes of KSDAGG, KSD median, KSD split, and KSD split extra data, all using a wild bootstrap to estimate the quantiles in the MNIST Normalizing Flow experimental setting considered in Section 4.5.

Since we consider a collection of bandwidths $\Lambda$, the time complexity of KSDAGG is $O(|\Lambda|(B_1 + B_2)N^2)$ and the one for median KSD is $O(B_1N^2)$ where, in the setting considered, we have $|\Lambda| = 21$ and $B_1 = B_2 = 500$. While KSDAGG takes roughly 10 times longer to run than KSD median, we could have expected a larger difference looking at the time complexities. This can be explained by the fact that there are two major time-consuming steps: (i) computing the kernel matrices, and (ii) computing the wild bootstrap samples. While (i) has complexity $O(N^2)$ and (ii) complexity $O(BN^2 + NB^2)$, the constant for step (i) is much larger than the one for step (ii) since (i) requires computing the matrix of pairwise distances while (ii) only consists in computing some matrix multiplications. Note that to compute the $|\Lambda|$ kernel matrices for KSDAGG, we need to compute the matrix of pairwise distances only once.

When splitting the data, the computationally expensive step is to select the bandwidth. All the $|\Lambda|$ kernel matrices need to be computed, as for KSDAGG, which is the expensive step (i). The KSD split test runs only slightly faster than the KSD split extra test; it runs faster than KSDAGG but their runtimes are of the same order of magnitude.

## F Details on MNIST Normalizing Flow experiment

We provide details about the results reported in Figure 3b of the MNIST Normalizing Flow experiment presented in Section 4.5. In that case, the median bandwidth $\lambda_{med}$ is on average 2437. The collection of bandwidths considered for KSDAGG is $\Lambda(-20, 0) := \{2^i \lambda_{med} : i = -20, \dots, 0\}$. When
KSDAGG rejects the null hypothesis, the smallest bandwidth $2^{-20}\lambda_{\text{med}} \approx 2^{-20} \cdot 2437 \approx 0.002$ (among others) rejects the single test with adjusted level. The bandwidth selected by KSD split extra data is the largest bandwidth $2^0\lambda_{\text{med}} \approx 2437$. The bandwidth selection is performed by maximizing the criterion in Equation (9). Sutherland et al. (2017) showed that this is equivalent to maximizing asymptotic power in the case of the MMD; the same result holds straightforwardly for the KSD due to similar asymptotic properties. The criterion only maximizes asymptotic power and has no guarantee when using limited data. In this high-dimensional setting ($d = 784$) with sample sizes smaller than 500, the asymptotic regime is not attained, and the criterion used for bandwidth selection does not maximize the power in this non-asymptotic setting. So, even though split extra has access to some extra data, it does not have an accurate criterion to select the bandwidth and ends up selecting the largest bandwidth, which is not well-adapted to this problem. This explains the low power obtained by KSD split extra data.

## G Details on our aggregation procedure

**More powerful than Bonferroni correction.** The Bonferroni correction for multiple testing corresponds to using the adjusted level $\alpha/|K|$ for each of the $|K|$ tests. With the correction used for KSDAGG with uniform weights, the adjusted level for each test is $u_{\alpha}/|K|$ with $u_{\alpha}$ defined in Equation (6). It can be shown that $u_{\alpha} \ge \alpha$ (Albert et al., 2022, Lemma 4), which means that KSDAGG will always reject the null when the test with Bonferroni correction would reject it. When using uniform weights, the multiple testing correction we use is guaranteed to always lead to a test as least as powerful as the one using the Bonferroni correction. The proof that the Bonferroni correction guarantees control of the probability of type I error uses a loose union bound argument, the method we use aims to tighten this loose upper bound. We present edge-case examples to illustrate the strengths of our multiple testing correction. First, assume that the $|K|$ events are all disjoints, then the union bound is tight and both methods yield the adjusted levels $\alpha/|K|$. Second, assume that all events are the same (e.g. same kernels), then the Bonferroni correction still yields adjusted levels $\alpha/|K|$, while our multiple testing strategy can detect that there is nothing to correct for and provide 'adjusted' levels $\alpha$. We note that, given some fixed weights and level, there does not exist a single $u_{\alpha}$ associated to them, it will depend on what the events are.

**Choice of weights.** Without any prior knowledge (which is often the case in practice), we recommend using uniform weights since we do not expect any particular kernels/bandwidths to be better-suited than others. If the user has some prior knowledge of which kernels/bandwidths would be better for the task considered, then those should be given higher weights. We allow for weights whose sum is strictly smaller than 1 simply for the convenience of being able to add a new kernel/bandwidth with a new weight without changing the previous weights (as in Theorem 3.5 with weights $6/\pi^2\ell^2$ for $\ell \in \mathbb{N} \setminus \{0\}$). Multiplying all the weights by a constant simply results in dividing the correction $u_{\alpha}$ defined in Equation (6) by the same constant. This means that the product $u_{\alpha}w_{\lambda}$ remains the same, and hence the definition of the aggregated test is not affected by this change. For simplicity, in practice, we use weights whose sum is equal to 1.

**Interpretability.** When KSDAGG rejects the null hypothesis, we can check which specific kernels reject the adjusted tests: this provides the user with information about the kernels which are well-adapted to the problem considered. The 'best' selection of kernels is naturally returned as a side-effect of the test (without requiring data splitting). This contributes to the interpretability of KSDAGG, for instance in the case where different kernels prioritise different features.

**Extension to uncountable collections of kernels.** Our KSDAGG test aggregates over a finite collection of kernels. While even for large collections of kernels KSDAGG retains high power (due our multiple testing correction), it would be theoretically interesting to be able to consider an uncountable/continuous collection of kernels (parametrised by the kernel bandwidth on the positive real line for example). Optimizing the kernel parameters continuously can be done by selecting the parameters which maximise a proxy for test power as in Equation (9), but to the best of our knowledge, this has currently never been done without data splitting (which usually result in a loss of power). The extension to the case of uncountable collections remains a very challenging problem, and it would be of great theoretical interest to derive such a method in the future. However, as shown empirically in Appendix B, with a discretisation using only 10 bandwidths, our aggregated test KSDAGG\* already obtains the same power as the test using the continuous interval (in limit of the discretisation).
**Extension to KSD tests with adaptive features.** Our proposed test, KSDAGG, provides a solution to the problem of KSD adaptivity in the goodness-of-fit framework without requiring data splitting. A potential future direction of interest could be to tackle the adaptivity problem of the KSD-based linear-time goodness-of-fit test proposed by Jitkrittum et al. (2017). In this setting, the data is split to select feature locations (and the kernel bandwidth), the KSD test is then run using those adaptive features. A challenging problem would be to obtain those adaptive features using an aggregation procedure which avoids splitting the data.

## H Adversarially constructed settings in which KSDAGG might fail

**Median bandwidth test.** In settings (possibly adversarially constructed) where the median bandwidth is the ‘best’ bandwidth, the KSD test with median bandwidth could in theory be more competitive than KSDAGG since by considering a large collection of bandwidths we are not only considering the ‘best’ median bandwidth, but also ‘worse’ bandwidths. However, in practice, we cannot know in advance which bandwidth would perform well, and KSDAGG retains power even for large collections of bandwidths (21 bandwidths considered in MNIST Normalizing Flow experiment in Section 4.5).

**Data splitting test.** Another setting (also possibly adversarially constructed) could be one in which the ‘best’ bandwidth lies in between two bandwidths of our collection which are ‘worse’ bandwidths. In that case, a test which uses data splitting to select an ‘optimal’ bandwidth would be able to select it, however, one must bear in mind the loss of power due to data splitting. In such a setting, the power of KSDAGG could then be improved by considering a finer collection of bandwidths. However, this issue did not arise in our experiments where we found that the aggregation procedure outperformed the competing approaches.

## I Proofs

### I.1 Proof of Theorem 3.1

Note that

$$ 
\begin{aligned}
\mathbb{P}_q(\Delta_{\alpha}^k(\mathbb{X}_N) = 0) &= \mathbb{P}_q(\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) \le \hat{q}_{1-\alpha}^k) \\
&= \mathbb{P}_q(\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) - \text{KSD}_{p,k}^2(q) \le \hat{q}_{1-\alpha}^k - \text{KSD}_{p,k}^2(q))
\end{aligned}
$$

where

$$ \text{KSD}_{p,k}^2(q) = \mathbb{E}_q[\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N)] $$

since $\widehat{\text{KSD}}_{p,k}^2$ is an unbiased estimator. By Chebyshev’s inequality (Chebyshev, 1899), we know that

$$ 
\begin{aligned}
\beta &\ge \mathbb{P}_q \left( \left| \text{KSD}_{p,k}^2(q) - \widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) \right| \ge \sqrt{\frac{\text{var}_q(\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N))}{\beta}} \right) \\
&\ge \mathbb{P}_q \left( \text{KSD}_{p,k}^2(q) - \widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) \ge \sqrt{\frac{\text{var}_q(\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N))}{\beta}} \right) \\
&= \mathbb{P}_q \left( \widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) - \text{KSD}_{p,k}^2(q) \le -\sqrt{\frac{\text{var}_q(\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N))}{\beta}} \right).
\end{aligned}
$$
We deduce that $\mathbb{P}_q(\Delta_{\alpha}^k(\mathbb{X}_N) = 0) \le \beta$ if

$$ \hat{q}_{1-\alpha}^k - \text{KSD}_{p,k}^2(q) \le -\sqrt{\frac{\text{var}_q(\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N))}{\beta}} \\ \text{KSD}_{p,k}^2(q) \ge \hat{q}_{1-\alpha}^k + \sqrt{\frac{\text{var}_q(\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N))}{\beta}}. \quad (10) $$

The condition in Equation (10), which controls the probability of type I error of $\Delta_{\alpha}^k$, needs to hold with high probability since $\hat{q}_{1-\alpha}^k$ depends on the randomness of, either the new samples drawn from $p$ for the parametric bootstrap in Equation (4), or of the Rademacher random variables for the wild bootstrap in Equation (5).

We want to derive a condition in terms of $\|p-q\|_2^2$ rather than in terms of $\text{KSD}_{p,k}^2(q)$ as in Equation (10). For this, using Stein's identity, we obtain

$$ \begin{aligned} \langle \psi, T_{h_{p,k}} \psi \rangle_2 &= \int_{\mathbb{R}^d} \psi(y) (T_{h_{p,k}} \psi)(y) \, dy \\ &= \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} h_{p,k}(x, y) \psi(x) \psi(y) \, dx \, dy \\ &= \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} h_{p,k}(x, y) (p(x) - q(x))(p(y) - q(y)) \, dx \, dy \\ &= \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} h_{p,k}(x, y) q(x) q(y) \, dx \, dy \\ &= \mathbb{E}_{q,q}[h_{p,k}(X, Y)] \\ &= \text{KSD}_{p,k}^2(q) \end{aligned} $$

which gives

$$ \text{KSD}_{p,k}^2(q) = \langle \psi, T_{h_{p,k}} \psi \rangle = \frac{1}{2} (\|\psi\|_2^2 + \|T_{h_{p,k}} \psi\|_2^2 - \|\psi - T_{h_{p,k}} \psi\|_2^2). $$

To guarantee $\mathbb{P}_q(\Delta_{\alpha}^k(\mathbb{X}_N) = 0) \le \beta$, an equivalent condition to the one presented in Equation (10) is then

$$ \|\psi\|_2^2 \ge \|\psi - T_{h_{p,k}} \psi\|_2^2 - \|T_{h_{p,k}} \psi\|_2^2 + 2\hat{q}_{1-\alpha}^k + 2\sqrt{\frac{\text{var}_q(\widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N))}{\beta}} \quad (11) $$

which needs to hold with high probability over the randomness in $\hat{q}_{1-\alpha}^k$. We now upper bound the quantile and variance terms in Equation (11) to obtain a new condition ensuring control of the type II error.

To bound the quantile term $\hat{q}_{1-\alpha}^k$ using $B_1$ wild bootstrapped statistics, using the Dvoretzky-Kiefer-Wolfowitz inequality (Dvoretzky et al., 1956; Massart, 1990) as done by Schrab et al. (2021, Appendix E.4), it is sufficient to bound the quantile $\hat{q}_{1-\alpha}^{k,\infty}$ using wild bootstrapped statistics without finite approximation provided that $B_1 \ge \frac{3}{\alpha^2}(\log(\frac{8}{\beta}) + \alpha(1-\alpha))$. Using the concentration bound for i.i.d. Rademacher chaos of de la Peña and Giné (1999, Corollary 3.2.6) as presented by Kim et al. (2022, Equation 39), there exists some constant $C > 0$ such that

$$ \mathbb{P}_e \left( \left| \frac{1}{N(N-1)} \sum_{1 \le i \ne j \le N} \epsilon_i \epsilon_j h_{p,k}(X_i, X_j) \right| \ge t \mid \mathbb{X}_N \right) \le 2 \exp \left( - \frac{CtN(N-1)}{\sqrt{\sum_{1 \le i \ne j \le N} h_{p,k}(X_i, X_j)^2}} \right) $$

which gives

$$ \hat{q}_{1-\alpha}^{k,\infty} \le C \frac{\log(2/\alpha)}{N(N-1)} \sqrt{\sum_{1 \le i \ne j \le N} h_{p,k}(X_i, X_j)^2} $$
for some different constant $C > 0$. Using Markov's inequality, for any $\delta \in (0, 1)$, with probability at least $1 - \delta$ we have

$$ \sum_{1 \le i \ne j \le N} h_{p,k}(X_i, X_j)^2 \le \frac{1}{\delta} N(N-1) \mathbb{E}_{q,q} [h_{p,k}(X, Y)^2] = \frac{C_k N(N-1)}{\delta} $$

which gives

$$ \hat{q}_{1-\alpha}^{k, \infty} \le C \frac{\log(2/\alpha)}{N(N-1)} \sqrt{\frac{C_k N(N-1)}{\delta}} = \frac{C}{\sqrt{\delta}} \frac{\log(2/\alpha)}{\sqrt{N(N-1)}} \sqrt{C_k} \le C \log\left(\frac{1}{\alpha}\right) \frac{\sqrt{C_k}}{N} $$

since $\alpha \in (0, e^{-1})$, where in the last inequality $C > 0$ is a different constant depending on $\delta$. Using Dvoretzky-Kiefer-Wolfowitz inequality (Dvoretzky et al., 1956) as explained above, we deduce that

$$ 2\hat{q}_{1-\alpha}^{k} \le C \log\left(\frac{1}{\alpha}\right) \frac{\sqrt{C_k}}{N} \qquad (12) $$

with arbitrarily high probability, for some constant $C > 0$. The same quantile bound also holds for the parametric bootstrap provided that $n/\sqrt{C_k} \ge \ln(1/\alpha)$ as derived by Albert et al. (2022, Appendix C.4.1). Their reasoning holds in our setting by replacing $1/\sqrt{\lambda_1 \cdots \lambda_p \mu_1 \cdots \mu_q}$ by $\sqrt{C_k}$, this is justified as for a kernel $k_\lambda$ with bandwidths $\lambda_1, \dots, \lambda_d$ we have $\mathbb{E}_{q,q}[k_\lambda(X, Y)^2] \le C/\lambda_1 \cdots \lambda_d$.

We now bound the variance term in Equation (10). Using the result of Albert et al. (2022, Equation 6), there exists $C > 0$ such that

$$ \mathrm{var}_q(\widehat{\mathrm{KSD}}_{p,k}^2(\mathbb{X}_N)) \le C\left(\frac{\sigma_1^2}{N} + \frac{\sigma_2^2}{N^2}\right) $$

where

$$ \sigma_2^2 := \mathbb{E}_{q,q}[h_{p,k}(X, Y)^2] = C_k $$

and

$$ \begin{aligned} \sigma_1^2 &:= \mathbb{E}_{Y \sim q} \left[ \left( \mathbb{E}_{X \sim q} [h_{p,k}(X, Y)] \right)^2 \right] \\ &= \mathbb{E}_{Y \sim q} \left[ \left( (T_{h_{p,k}} \psi)(Y) \right)^2 \right] \\ &= \int_{\mathbb{R}^d} \left( (T_{h_{p,k}} \psi)(y) \right)^2 q(y) \, dy \\ &\le \|q\|_{\infty} \int_{\mathbb{R}^d} \left( (T_{h_{p,k}} \psi)(y) \right)^2 \, dy \\ &\le M \|T_{h_{p,k}} \psi\|_2^2 \end{aligned} $$

since

$$ \begin{aligned} (T_{h_{p,k}} \psi)(y) &:= \int_{\mathbb{R}^d} h_{p,k}(x, y) \psi(x) \, dx \\ &= \int_{\mathbb{R}^d} h_{p,k}(x, y) p(x) \, dx - \int_{\mathbb{R}^d} h_{p,k}(x, y) q(x) \, dx \\ &= 0 - \int_{\mathbb{R}^d} h_{p,k}(x, y) q(x) \, dx \\ &= -\mathbb{E}_{X \sim q} [h_{p,k}(X, y)] \end{aligned} $$

by Stein's identity. We deduce that

$$ \mathrm{var}_q(\widehat{\mathrm{KSD}}_{p,k}^2(\mathbb{X}_N)) \le C \left( \frac{\|T_{h_{p,k}} \psi\|_2^2}{N} + \frac{C_k}{N^2} \right) $$

for some different constant $C > 0$ depending on $M$ and $d$.
Using the classical inequalities $\sqrt{x+y} \le \sqrt{x} + \sqrt{y}$ and $2\sqrt{xy} \le x+y$ for $x, y > 0$, which are also considered in the works of Fromont et al. (2013), Albert et al. (2022) and Schrab et al. (2021), we obtain

$$ 
\begin{aligned}
2\sqrt{\frac{\mathrm{var}_q(\widehat{\mathrm{KSD}}_{p,k}^2(\mathbb{X}_N))}{\beta}} &\le 2\sqrt{C \frac{\|T_{h_{p,k}}\psi\|_2^2}{\beta N}} + C \frac{C_k}{\beta N^2} \\
&\le 2\sqrt{\|T_{h_{p,k}}\psi\|_2^2 \frac{C}{\beta N}} + 2\sqrt{C \frac{C_k}{\beta N^2}} \\
&\le \|T_{h_{p,k}}\psi\|_2^2 + \frac{C}{\beta N} + 2\sqrt{C} \frac{\sqrt{C_k}}{\sqrt{\beta N}} \\
&\le \|T_{h_{p,k}}\psi\|_2^2 + (C + 2\sqrt{C}) \frac{\sqrt{C_k}}{\beta N} \\
&\le \|T_{h_{p,k}}\psi\|_2^2 + C \log\left(\frac{1}{\alpha}\right) \frac{\sqrt{C_k}}{\beta N}
\end{aligned}
\quad (13) $$

since $\alpha \in (0, e^{-1})$, $\beta \in (0, 1)$ and $C_k \ge 1$ by assumption, where the constant $C > 0$ is different on the last line. By considering the largest of the two constants in the quantile and variance bounds of Equations (12) and (13) multiplied by two, we obtain

$$ 2\hat{q}_{1-\alpha}^k + 2\sqrt{\frac{\mathrm{var}_q(\widehat{\mathrm{KSD}}_{p,k}^2(\mathbb{X}_N))}{\beta}} \le \|T_{h_{p,k}}\psi\|_2^2 + C \log\left(\frac{1}{\alpha}\right) \frac{\sqrt{C_k}}{\beta N}. $$

By combining this bound with the condition in Equation (11), we get that $\mathbb{P}_q(\Delta_{\alpha}^k(\mathbb{X}_N) = 0) \le \beta$ if

$$ \|\psi\|_2^2 \ge \|\psi - T_{h_{p,k}}\psi\|_2^2 - \|T_{h_{p,k}}\psi\|_2^2 + \|T_{h_{p,k}}\psi\|_2^2 + C \log\left(\frac{1}{\alpha}\right) \frac{\sqrt{C_k}}{\beta N}, $$

$$ \|\psi\|_2^2 \ge \|\psi - T_{h_{p,k}}\psi\|_2^2 + C \log\left(\frac{1}{\alpha}\right) \frac{\sqrt{C_k}}{\beta N}, $$

which concludes the proof.

## I.2 Proof of Proposition 3.2

Recall that the correction term in Equation (6) is defined as

$$ u_{\alpha} := \sup \left\{ u \in (0, \min_{k \in \mathcal{K}} w_k^{-1}) : \frac{1}{B_2} \sum_{b=1}^{B_2} 1 \left( \max_{k \in \mathcal{K}} (\tilde{K}_k^b - \hat{q}_{1-uw_k}^k) > 0 \right) \le \alpha \right\} $$

where $\hat{q}_{1-uw_k}^k = \tilde{K}_k^{\alpha \cdot [(B_1+1)(1-uw_k)]}$, as defined in Equation (3). Hence, we have

$$ \frac{1}{B_2} \sum_{b=1}^{B_2} 1 \left( \max_{k \in \mathcal{K}} (\tilde{K}_k^b - \hat{q}_{1-u_{\alpha}w_k}^k) > 0 \right) \le \alpha. \quad (14) $$

Recall that our estimator is $\widehat{\mathrm{KSD}}_{p,k}^2(\mathbb{X}_N)$ where $\mathbb{X}_N := (X_i)_{i=1}^N$ are drawn from $q$. Recall also from Equation (4), that for the parametric bootstrap each element $\tilde{K}_k^b$ is computed by drawing new samples $(X'_i)_{i=1}^{N'}$ from the model $p$ and computing $\widehat{\mathrm{KSD}}_{p,k}^2((X'_i)_{i=1}^{N'})$. Hence, at every sample size, under the null hypothesis $\mathcal{H}_0: p = q$ the quantities $\tilde{K}_k^b$ and $\widehat{\mathrm{KSD}}_{p,k}^2(\mathbb{X}_N)$ are identically distributed. For the wild bootstrap, as defined in Equation (5), Chwialkowski et al. (2014, 2016) show that, under $\mathcal{H}_0$, $\tilde{K}_k^b$ and $\widehat{\mathrm{KSD}}_{p,k}^2(\mathbb{X}_N)$ have the same asymptotic distribution (details are presented in Appendix D).
Therefore, by taking the expectation in Equation (14), we obtain

$$
\begin{align*}
\alpha &\ge \mathbb{P}_p \left( \max_{k \in \mathcal{K}} \left( \widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) - \hat{q}_{1-u_\alpha w_k}^k \right) > 0 \right) \\
&= \mathbb{P}_p \left( \widehat{\text{KSD}}_{p,k}^2(\mathbb{X}_N) > \hat{q}_{1-u_\alpha w_k}^k \text{ for some } k \in \mathcal{K} \right) \\
&= \mathbb{P}_p \left( \Delta_{u_\alpha w_k}^k(\mathbb{X}_N) = 1 \text{ for some } k \in \mathcal{K} \right) \\
&= \mathbb{P}_p \left( \Delta_{\alpha}^K(\mathbb{X}_N) = 1 \right)
\end{align*}
$$

which holds non-asymptotically using the parametric bootstrap, and asymptotically using the wild bootstrap. Note that this is different from the two-sample case where using a wild bootstrap yields well-calibrated non-asymptotic levels (Schrab et al., 2021, Proposition 8).

The same reasoning holds by replacing $u_\alpha$ by any value $u \in (0, u_\alpha)$. In particular, it holds for the lower bound of the interval obtained by performing the bisection method to approximate the supremum in the definition of $u_\alpha$. This lower bound is the one we use in practice, as shown in Algorithm 1 with the step 'u_α = u_min'. We have proved that the test with correction term $u_\alpha$ approximated with a bisection method also achieves the desired level $\alpha$.

I.3 Proof of Theorem 3.3

As explained in Appendix I.2 and utilised in Algorithm 1, we use the lower bound $\hat{u}_\alpha$ of the interval obtained by performing the bisection method to approximate the supremum in the definition of $u_\alpha$ in Equation (6). The assumptions $B_2 \ge \frac{8}{\alpha^2} \log(\frac{2}{\beta})$ and $B_3 \ge \log_2(\frac{4}{\alpha} \min_{k \in \mathcal{K}} w_k^{-1})$ ensure that $\hat{u}_\alpha \ge \alpha/2$ as shown by Schrab et al. (2021, Appendix E.9).

The probability of type II error of $\Delta_{\alpha}^{\mathcal{K}}$ is

$$
\begin{align*}
\mathbb{P}_q(\Delta_{\alpha}^{\mathcal{K}}(\mathbb{X}_N) = 0) &= \mathbb{P}_q(\Delta_{\hat{u}_{\alpha} w_k}^{k}(\mathbb{X}_N) = 0 \text{ for all } k \in \mathcal{K}) \\
&\le \mathbb{P}_q(\Delta_{\hat{u}_{\alpha} w_k}^{k}(\mathbb{X}_N) = 0 \text{ for some } k \in \mathcal{K}) \\
&\le \mathbb{P}_q(\Delta_{\alpha w_k/2}^{k}(\mathbb{X}_N) = 0 \text{ for some } k \in \mathcal{K}).
\end{align*}
$$

To guarantee $\mathbb{P}_q(\Delta_{\alpha}^{\mathcal{K}}(\mathbb{X}_N) = 0) \le \beta$, it is then sufficient to guarantee $\mathbb{P}_q(\Delta_{\alpha w_k^{*}/2}^{k}(\mathbb{X}_N) = 0) \le \beta$ for some $k^* \in \mathcal{K}$ to be specified shortly in Equation (15). By assumption, we have

$$
B_1 \geq \left( \max_{k \in \mathcal{K}} w_k^{-2} \right) \frac{12}{\alpha^2} \left( \log \left( \frac{8}{\beta} \right) + \alpha(1-\alpha) \right).
$$

In order to apply Theorem 3.1 to $\Delta_{\alpha_k}^k$ with $\alpha_k := \alpha w_k/2$ for $k \in \mathcal{K}$, we need to ensure that the condition on $B_1$ of Theorem 3.1 is satisfied for all $k \in \mathcal{K}$, that is

$$
B_1 \geq \frac{3}{\alpha_k^2} \left( \log \left( \frac{8}{\beta} \right) + \alpha_k (1 - \alpha_k) \right)
$$

for all $k \in \mathcal{K}$. Since $0 < \alpha_k < \alpha < e^{-1}$, we have $\alpha(1-\alpha) \ge \alpha_k(1-\alpha_k)$ for $k \in \mathcal{K}$. We also have

$$
\left( \max_{k \in \mathcal{K}} w_k^{-2} \right) \frac{12}{\alpha^2} \geq 3 \left( \frac{2}{\alpha w_k} \right)^2 = \frac{3}{\alpha_k^2}
$$

for all $k \in \mathcal{K}$. We deduce that

$$
\begin{align*}
B_1 &\ge \left( \max_{k \in \mathcal{K}} w_k^{-2} \right) \frac{12}{\alpha^2} \left( \log \left( \frac{8}{\beta} \right) + \alpha(1-\alpha) \right) \\
&\ge \frac{3}{\alpha_k^2} \left( \log \left( \frac{8}{\beta} \right) + \alpha_k(1-\alpha_k) \right)
\end{align*}
$$

for all $k \in \mathcal{K}$, and so, applying Theorem 3.1 to $\Delta_{\alpha_k}^k$ for $k \in \mathcal{K}$ is justified. We obtain that
$$ \mathbb{P}_q \left( \Delta_{\alpha w_k/2}^k (\mathbb{X}_N) = 0 \right) \le \beta \text{ if} $$

$$
\|\psi\|_2^2 \ge \|\psi - T_{h_{p,k}}\psi\|_2^2 + C \log \left(\frac{2}{\alpha w_k}\right) \frac{\sqrt{C_k}}{\beta N},
$$
or, with a different constant $C > 0$, if

$$ \|\psi\|_2^2 \ge \|\psi - T_{h_{p,k}}\psi\|_2^2 + C \log\left(\frac{1}{\alpha w_k}\right) \frac{\sqrt{C_k}}{\beta N} $$

since $\log(\frac{2}{\alpha w_k}) \le (\log(2) + 1) \log(\frac{1}{\alpha w_k})$ as $\alpha \in (0, e^{-1})$ and $w_k \in (0, 1)$. Now, let

$$ k^* := \underset{k \in \mathcal{K}}{\operatorname{argmin}} \left( \|\psi - T_{h_{p,k}}\psi\|_2^2 + C \log\left(\frac{1}{\alpha w_k}\right) \frac{\sqrt{C_k}}{\beta N} \right). \quad (15) $$

Finally, we have $\mathbb{P}_q(\Delta_{\alpha}^{\mathcal{K}}(\mathbb{X}_N) = 0) \le \beta$ if $\mathbb{P}_q(\Delta_{\alpha w_{k^*}/2}^{k^*}(\mathbb{X}_N) = 0) \le \beta$, that is, if

$$ \begin{aligned} \|\psi\|_2^2 &\ge \|\psi - T_{h_{p,k^*}}\psi\|_2^2 + C \log\left(\frac{1}{\alpha w_{k^*}}\right) \frac{\sqrt{C_{k^*}}}{\beta N} \\ &= \min_{k \in \mathcal{K}} \left( \|\psi - T_{h_{p,k}}\psi\|_2^2 + C \log\left(\frac{1}{\alpha w_k}\right) \frac{\sqrt{C_k}}{\beta N} \right), \end{aligned} $$

as desired.

## I.4 Proof of Theorem 3.5

Recall that we suppose that the following assumptions hold.

* The model density $p$ is strictly positive on its connected and compact support $S \subseteq \mathbb{R}^d$.

* The score function $\nabla \log p(x)$ is continuous and bounded on $S$.

* The support of the density $q$ is a connected and compact subset of $S$.

* The kernel used is a scaled Gaussian kernel $k_{\lambda}(x, y) := \lambda^{2-d} \exp(-\|x - y\|_2^2 / \lambda^2)$.

We introduce some notation. For some $c \in \mathbb{R}$, we write $c := (c, \dots, c) \in \mathbb{R}^d$. We write $(a_1, \dots, a_d) \le (b_1, \dots, b_d)$ when $a_i \le b_i$ for $i = 1, \dots, d$. We use $C, C'$ to denote some generic constants which may change on different lines.

Note that, by properties of compactness on $\mathbb{R}^d$, there exists some $a > 0$ such that $S \subseteq [-\frac{a}{2}, \frac{a}{2}]^d$. Since the score function $\nabla \log p(x)$ is continuous and bounded on the compact set $S$, there exists some $c_1 > 0$ such that

$$ |(\nabla \log p(x))_i| \le c_1 \quad \text{for } i = 1, \dots, d \text{ and for all } x \in S. $$

We work with a scaled Gaussian kernel with bandwidth $\lambda \le c_2$, defined as

$$ k_{\lambda}(x, y) := \lambda^{2-d} \exp\left(-\frac{\|x - y\|_2^2}{\lambda^2}\right) $$

which satisfies

$$ \pm \nabla_x k_{\lambda}(x, y) = \pm 2\lambda^{-d}(y-x) \exp\left(-\frac{\|x-y\|_2^2}{\lambda^2}\right) \le 2\lambda^{-d} \mathbf{a} \exp\left(-\frac{\|x-y\|_2^2}{\lambda^2}\right), $$

$$ \pm \nabla_y k_{\lambda}(x, y) = \pm 2\lambda^{-d}(x-y) \exp\left(-\frac{\|x-y\|_2^2}{\lambda^2}\right) \le 2\lambda^{-d} \mathbf{a} \exp\left(-\frac{\|x-y\|_2^2}{\lambda^2}\right), $$

$$ \left| \sum_{i=1}^{d} \frac{\partial}{\partial x_i \partial y_i} k_{\lambda}(x, y) \right| = \left| 2d\lambda^{-d} - 4\lambda^{-2-d} \|x - y\|_2^2 \right| \exp \left( -\frac{\|x - y\|_2^2}{\lambda^2} \right) $$

$$ \le (2d\lambda^{-d} + 4da^2\lambda^{-2-d}) \exp \left( -\frac{\|x - y\|_2^2}{\lambda^2} \right). $$
The Stein kernel associated to $k_\lambda$ with $\lambda c_3 \ge a$ satisfies

$$
\begin{aligned}
|h_{p,k_\lambda}(x,y)| & \le |(\nabla \log p(x)^\top \nabla \log p(y)) k_\lambda(x,y)| + |\nabla \log p(y)^\top \nabla_x k_\lambda(x,y)| \\
& \quad + |\nabla \log p(x)^\top \nabla_y k_\lambda(x,y)| + \left| \sum_{i=1}^d \frac{\partial}{\partial x_i \partial y_i} k_\lambda(x,y) \right| \\
& \le 2 (dc_1^2 c_2^2 \lambda^{-d} + 2dc_1 a \lambda^{-d} + 3dc_3^2 \lambda^{-d}) \exp \left( -\frac{\|x-y\|_2^2}{\lambda^2} \right) \\
& \le C_0 \pi^{-d/2} \lambda^{-d} \exp \left( -\frac{\|x-y\|_2^2}{\lambda^2} \right)
\end{aligned}
$$

where the constant is $C_0 := 2\pi^{d/2}(dc_1^2c_2^2 + 2dc_1a + 3dc_3^2)$. The Stein kernel can be upper bounded by a scaled Gaussian kernel on $S$. Writing

$$
\bar{k}_{\lambda}(x, y) := \pi^{-d/2} \lambda^{-d} \exp \left( - \frac{\|x - y\|_{2}^{2}}{\lambda^{2}} \right),
$$

which is of the form considered by Schrab et al. (2021, Section 3.1) with equal bandwidths for each
dimension, we have shown that

$$
|h_{p,k_{\lambda}}(x,y)| \le C_0 \bar{k}_{\lambda}(x,y) \quad \text{for all } x,y \in S. \tag{16}
$$

Writing $\psi := p - q$, recall from Theorem 3.1 that a sufficient condition for control of the probability of type II error by $\beta$ is

$$
\|\psi\|_2^2 \ge \|\psi - T_{h_{p,k_\lambda}} \psi\|_2^2 + C \log \left(\frac{1}{\alpha}\right) \frac{\sqrt{\mathbb{E}_{q,q}[h_{p,k_\lambda}(X,Y)^2]}}{\beta N}.
$$

Using the upper bound

$$
\mathbb{E}_{q,q}[h_{p,k_{\lambda}}(X, Y)^2] \le C_0^2 \mathbb{E}_{q,q}[\bar{k}_{\lambda}(X, Y)^2] \le C_0^2 \frac{M}{\lambda^d}
$$

where the last inequality holds as in Schrab et al. (2021, Equation (22)), we obtain the sufficient
condition

$$
\|\psi\|_2^2 \ge \|\psi - T_{h_{p,k_\lambda}} \psi\|_2^2 + C \frac{\log(1/\alpha)}{\beta N \lambda^{d/2}} \quad (17)
$$

where we have absorbed the term $C_0\sqrt{M}$ in the constant $C$.

We would like to upper bound the term $\|\psi - T_{h_{p,k_\lambda}} \psi\|_2^2$ by $\|\psi - T_{\bar{k}_\lambda} \psi\|_2^2$ but this is not possible using Equation (16). Instead, we can use the triangle inequality to get

$$
\begin{aligned}
\| \psi - T_{h_{p,k_\lambda}} \psi \|_2^2 &\le \| \psi - T_{\bar{k}_\lambda} \psi \|_2^2 + \| T_{\bar{k}_\lambda} \psi - T_{h_{p,k_\lambda}} \psi \|_2^2 \\
&\le \| \psi - T_{\bar{k}_\lambda} \psi \|_2^2 + \left( \| T_{\bar{k}_\lambda} \psi \|_2 + \| T_{h_{p,k_\lambda}} \psi \|_2 \right)^2 \\
&\le \| \psi - T_{\bar{k}_\lambda} \psi \|_2^2 + 2 \left( \| T_{\bar{k}_\lambda} \psi \|_2^2 + \| T_{h_{p,k_\lambda}} \psi \|_2^2 \right) \\
&\le \| \psi - T_{\bar{k}_\lambda} \psi \|_2^2 + 2(\tilde{C}_0^2 + 1) \| T_{\bar{k}_\lambda} \psi \|_2^2
\end{aligned}
$$

since

$$
\| T_{h_{p,k_{\lambda}}} \psi \|_2 \le \tilde{C}_0 \| T_{\bar{k}_{\lambda}} \psi \|_2
$$

for some constant $\tilde{C}_0 > 0$ by assumption (see Appendix A.4).

Recall that we assume that $\psi \in S_d^{s,t}(R, L)$ for some $L$ to be determined, in particular $\psi \in S_d^s(R)$.

For the term $\|\psi - T_{\bar{k}_{\lambda}} \psi\|_2^2$, we use the fact that $\psi \in S_d^s(R)$. The term $\|\psi - T_{\bar{k}_{\lambda}} \psi\|_2^2$ can then be
upper bounded exactly as done by Schrab et al. (2021, Appendix E.6) with the difference that we
choose $\tilde{t} > 0$ such that $S_1 < 0.5$ rather than $S_1 < 1$ ($S_1$ is defined in Schrab et al. (2021, Appendix E.6)). Following their reasoning, since $\psi \in S_d^s(R)$, we obtain that there exists some $S_1 \in (0, 0.5)$ and some constant $C > 0$ (depending on $d$, $s$ and $R$) such that

$$ \|\psi - T_{\bar{k}_{\lambda}} \psi\|_{2}^{2} \leq S_{1}^{2} \|\psi\|_{2}^{2} + C \lambda^{2s} \leq \frac{1}{4} \|\psi\|_{2}^{2} + C \lambda^{2s}. $$

To upper bound the term $\|T_{\bar{k}_{\lambda}} \psi\|_2^2$, we modify the reasoning of Schrab et al. (2021, Appendix E.6) used for the term $\|\psi - T_{\bar{k}_{\lambda}} \psi\|_2^2$, and utilise the restricted Sobolev ball regularity assumption. First, note that since $\bar{k}_{\lambda}$ is translation-invariant, $T_{\bar{k}_{\lambda}}$ is a convolution as

$$ (T_{\bar{k}_{\lambda}} \psi)(y) = \int_{S} \bar{k}_{\lambda}(x, y) \psi(x) \, dx = \int_{S} \varphi_{\lambda}(x - y) \psi(x) \, dx = (\psi * \varphi_{\lambda})(y) $$

where

$$ \varphi_{\lambda}(u) := \prod_{i=1}^{d} \lambda^{-1} K\left(\frac{u_i}{\lambda}\right), \quad \text{and} \quad K(u) := \pi^{-1/2} \exp(-u^2). $$

By properties of Fourier transforms, we have $\hat{\varphi}_{\lambda}(\xi) = \prod_{i=1}^{d} \hat{K}(\lambda \xi_i)$ for $\xi \in \mathbb{R}^d$. Note that

$$ \left| \prod_{i=1}^{d} \hat{K}(\xi_i) \right| \leq \prod_{i=1}^{d} \int_{\mathbb{R}} |K(x)e^{-ix\xi_i}| \, dx = \prod_{i=1}^{d} \int_{\mathbb{R}} |K(x)| \, dx = 1. $$

For some $L$ to be determined, by assumption, we have $\psi \in S_d^{s,t}(R, L)$, so

$$ \int_{\|\xi\|_2 \leq t} |\hat{\psi}(\xi)|^2 d\xi \leq \frac{1}{L} \int_{\mathbb{R}^d} |\hat{\psi}(\xi)|^2 d\xi = \frac{\|\hat{\psi}\|_2^2}{L}. $$

For $s > 0$, define

$$ T_s := \sup_{\|\xi\|_2 > t} \frac{\left| \prod_{i=1}^{d} \hat{K}(\xi_i) \right|}{\|\xi\|_2^s} \leq \frac{1}{ts} < \infty. $$

Then, using Plancherel's Theorem as in Schrab et al. (2021, Appendix E.6), we obtain

$$ \begin{aligned} & (2\pi)^d \|T_{\bar{k}_{\lambda}} \psi\|_2^2 \\ &= (2\pi)^d \|\psi * \varphi_{\lambda}\|_2^2 \\ &= \|\hat{\varphi}_{\lambda} \hat{\psi}\|_2^2 \\ &= \int_{\mathbb{R}^d} |\hat{\varphi}_{\lambda}(\xi)|^2 |\hat{\psi}(\xi)|^2 d\xi \\ &= \int_{\mathbb{R}^d} \left| \prod_{i=1}^{d} \hat{K}(\lambda \xi_i) \right|^2 |\hat{\psi}(\xi)|^2 d\xi \\ &= \int_{\|\xi\|_2 \leq t} \left| \prod_{i=1}^{d} \hat{K}(\lambda \xi_i) \right|^2 |\hat{\psi}(\xi)|^2 d\xi + \int_{\|\xi\|_2 > t} \left| \prod_{i=1}^{d} \hat{K}(\lambda \xi_i) \right|^2 |\hat{\psi}(\xi)|^2 d\xi \\ &\leq \int_{\|\xi\|_2 \leq t} |\hat{\psi}(\xi)|^2 d\xi + T_s^2 \int_{\|\xi\|_2 > t} \|\lambda \xi\|_2^{2s} |\hat{\psi}(\xi)|^2 d\xi \\ &\leq \frac{1}{L} \|\hat{\psi}\|_2^2 + \lambda^{2s} T_s^2 \int_{\mathbb{R}^d} \|\xi\|_2^{2s} |\hat{\psi}(\xi)|^2 d\xi \\ &\leq \frac{(2\pi)^d}{L} \|\psi\|_2^2 + \lambda^{2s} T_s^2 (2\pi)^d R^2. \end{aligned} $$
Combining those results, the upper bound on $||\psi - T_{h_{p,k_\lambda}}\psi||_2^2$ becomes

$$ \begin{aligned} ||\psi - T_{h_{p,k_\lambda}}\psi||_2^2 &\le ||\psi - T_{\bar{k}_\lambda}\psi||_2^2 + 2(\tilde{C}_0^2 + 1) ||T_{\bar{k}_\lambda}\psi||_2^2 \\ &\le \left(\frac{1}{4}||\psi||_2^2 + C\lambda^{2s}\right) + \left(\frac{2(\tilde{C}_0^2 + 1)}{L}||\psi||_2^2 + C'\lambda^{2s}\right) \\ &= \left(\frac{1}{4}||\psi||_2^2 + C\lambda^{2s}\right) + \left(\frac{1}{4}||\psi||_2^2 + C'\lambda^{2s}\right) \\ &\le \frac{1}{2}||\psi||_2^2 + C\lambda^{2s}, \end{aligned} $$

where we let $L := 8(\tilde{C}_0^2 + 1)$. The condition in Equation (17) then becomes

$$ ||\psi||_2^2 \ge \frac{1}{2} ||\psi||_2^2 + C\lambda^{2s} + C' \frac{\log(1/\alpha)}{\beta N \lambda^{d/2}} $$

which gives

$$ ||\psi||_2^2 \ge C \left( \lambda^{2s} + \frac{\log(1/\alpha)}{\beta N \lambda^{d/2}} \right). \quad (18) $$

Having proved the power guaranteeing condition in Equation (18), we can then derive rates for the KSD and KSDAGG tests as similarly done by Schrab et al. (2021). Set $\lambda := N^{-2/(4s+d)}$ to get the condition

$$ ||\psi||_2^2 \ge C(N^{-4s/(4s+d)} + N^{-1}N^{d/(4s+d)}) = CN^{-4s/(4s+d)} $$

which is the minimax rate over (unrestricted) Sobolev balls $S_d^s(R)$ (Li and Yuan, 2019; Balasubramanian et al., 2021; Albert et al., 2022; Schrab et al., 2021).

In practice, the smoothness parameter of the Sobolev ball is not known, so we cannot set $\lambda := N^{-2/(4s+d)}$ (i.e. it cannot be implemented). Instead, we can use our aggregated test KSDAGG. Adapting the proof of Theorem 3.3 with the single test power condition Equation (18), we obtain the aggregated test power condition

$$ ||\psi||_2^2 \ge C \min_{\lambda \in \Lambda} \left( \lambda^{2s} + \log \left( \frac{1}{\alpha w_{\lambda}} \right) \frac{1}{\beta N \lambda^{d/2}} \right). \quad (19) $$

Similarly to Schrab et al. (2021, Corollary 10), consider

$$ \Lambda := \left\{ 2^{-\ell} : \ell \in \{1, \dots, \left\lfloor \frac{2}{d} \log_2 \left( \frac{N}{\ln(\ln(N))} \right) \right\rfloor\} \right\}, \quad w_{\lambda} := \frac{6}{\pi^2 \ell^2}. $$

Let

$$ \ell_* := \left\lfloor \frac{2}{4s+d} \log_2 \left( \frac{N}{\ln(\ln(N))} \right) \right\rfloor \le \left\lfloor \frac{2}{d} \log_2 \left( \frac{N}{\ln(\ln(N))} \right) \right\rfloor. $$

The bandwidth $\lambda_* := 2^{-\ell_*} \in \Lambda$ satisfies

$$ \ln\left(\frac{1}{w_{\lambda_*}}\right) \le C \ln(\ell_*) \le C \ln(\ln(N)) $$

as $w_{\lambda_*} := 6\pi^{-2}(\ell_*)^{-2}$, and

$$ \frac{1}{2} \left( \frac{N}{\ln(\ln(N))} \right)^{-2/(4s+d)} \le \lambda_* \le \left( \frac{N}{\ln(\ln(N))} \right)^{-2/(4s+d)}. $$

We get

$$ \lambda_*^{2s} \le \left( \frac{N}{\ln(\ln(N))} \right)^{-4s/(4s+d)} $$

and

$$ \log\left(\frac{1}{\alpha w_{\lambda_*}}\right) \frac{1}{\beta N \lambda_*^{d/2}} \le C \left( \frac{N}{\ln(\ln(N))} \right)^{-1} \left( \frac{N}{\ln(\ln(N))} \right)^{d/(4s+d)} \le C \left( \frac{N}{\ln(\ln(N))} \right)^{-4s/(4s+d)}. $$

The KSDAGG power condition of Equation (19) then becomes

$$ ||\psi||_2^2 \ge C \left( \frac{N}{\ln(\ln(N))} \right)^{-4s/(4s+d)}, $$

which concludes the proof.