# Kernel-Based Tests for
Likelihood-Free Hypothesis Testing

Patrik Róbert Gerber
Department of Mathematics, MIT
Cambridge, MA 02139
<EMAIL_ADDRESS>
&Tianze Jiang
Department of Mathematics, MIT
Cambridge, MA 02139
<EMAIL_ADDRESS>
&Yury Polyanskiy
Department of EECS, MIT
Cambridge, MA 02139
<EMAIL_ADDRESS>
&Rui Sun
Department of Mathematics, MIT
Cambridge, MA 02139
<EMAIL_ADDRESS>
###### Abstract

Given $n$ observations from two balanced classes, consider the task of
labeling an additional $m$ inputs that are known to all belong to _one_ of the
two classes. Special cases of this problem are well-known: with complete
knowledge of class distributions ($n=\infty$) the problem is solved optimally
by the likelihood-ratio test; when $m=1$ it corresponds to binary
classification; and when $m\approx n$ it is equivalent to two-sample testing.
The intermediate settings occur in the field of likelihood-free inference,
where labeled samples are obtained by running forward simulations and the
unlabeled sample is collected experimentally. In recent work it was discovered
that there is a fundamental trade-off between $m$ and $n$: increasing the data
sample $m$ reduces the amount $n$ of training/simulation data needed. In this
work we (a) introduce a generalization where unlabeled samples come from a
mixture of the two classes – a case often encountered in practice; (b) study
the minimax sample complexity for non-parametric classes of densities under
maximum mean discrepancy (MMD) separation; and (c) investigate the empirical
performance of kernels parameterized by neural networks on two tasks:
detection of the Higgs boson and detection of planted DDPM generated images
amidst CIFAR-10 images. For both problems we confirm the existence of the
theoretically predicted asymmetric $m$ vs $n$ trade-off. Our code can be found
at https://github.com/Sr-11/LFI.

##  1 Likelihood-Free Inference

The goal of likelihood-free inference (LFI) [17; 23; 8; 12], also called
simulation-based inference (SBI), is to perform statistical inference in a
setting where the data generating process is a black-box, but can be
simulated. Given the ability to generate samples $X_{\theta}\sim
P_{\theta}^{\otimes n}$ for any parameter $\theta$, and given real-world data
$Z\sim P_{\theta^{\star}}^{\otimes m}$, we want to use our simulations to
learn about the truth $\theta^{\star}$. LFI is particularly relevant in areas
of science where we have precise but complex laws of nature, for which we can
do (stochastic) forward simulations, but can not directly compute the
(distribution) density $P_{\theta}$. The Bayesian community approached the
problem under the name of Approximate Bayesian Computation (ABC) [13; 44; 6].
More recent ML-based methods where regressors and classifiers are used to
summarize data, select regions of interest, approximate likelihoods or
likelihood-ratios [27; 28; 39; 14; 46] have also emerged for this challenge.

Despite empirical advances, the theoretical study of frequentist LFI is still
in its infancy. We focus on the non-parametric and non-asymptotic setting,
which we justify as follows. For applications where tight error control is
critical one might be reluctant to rely on asymptotics. More broadly, the non-
asymptotic regime can uncover new phenomena and provide insights for algorithm
design. Further, parametric models are clearly at odds with the black-box
assumption. Recently, [18] proposed likelihood-free hypothesis testing
($\mathsf{LFHT}$) as a simplified model and found minimax optimal tests for a
range of non-parametric distribution classes, thereby identifying a
fundamental simulation-experimentation trade-off between the number of
simulated observations $n$ and the size of the experimental data sample $m$.
Here we extend [18] to a new setting designed to model experimental setups
more truthfully and derive sample complexity (upper and lower bounds) for
kernel-based tests over non-parametric classes.

While minimax optimal, the algorithms of [18] are impractical as they rely on
discretizing the observations on a regular grid. Thus, both in our theory as
well as experiments we turn to kernel methods which provide an empirically
more powerful set of algorithms that have shown success in non-parametric
testing [19; 20; 29; 21; 36].

Contributions Our contributions are twofold. Theoretically, we introduce
_mixed likelihood-free hypothesis testing_ ($\mathsf{mLFHT}$), which is a
generalization of ($\mathsf{LFHT}$) and provides a better model of
applications such as the search for new physics [11; 35]. We describe a robust
kernel-based test and derive both upper and lower bounds on its minimax sample
complexity over a large non-parametric class of densities, extending results
of [34; 18]. Although the simulation-experimentation ($m$ vs $n$) trade-off
has been proven in the minimax sense (that is, for some worst-case data
distribution), it is not clear whether it actually occurs in real data. Our
second contribution is the empirical confirmation of the existence of an
asymmetric trade-off, cf. Figure 1. To this end we construct state-of-the-art
tests building on ideas of [45; 36] on learning good kernels from the data. We
execute this program in two settings: the Higgs boson discovery [5], and
detecting diffusion [24] generated images planted in the CIFAR-10 [32]
dataset.

Figure 1: $n$ versus $m$ trade-off for the Higgs and CIFAR experiments using
(4). Error probabilities are estimated by normal approximation for Higgs and
simulated for CIFAR.

###  1.1 $\mathsf{LFHT}$ and the Simulation-Experimentation Trade-off

Let $\mathcal{X}$ be some measurable space and suppose we have i.i.d. samples
$X,Y$ each of size $n$ from two unknown distributions $P_{X},P_{Y}$ on
$\mathcal{X}$, as well as a third i.i.d. sample $Z\sim P_{Z}$ of size $m$. In
the context of LFI, we may think of the samples $X,Y$ as being generated by
our simulator, and $Z$ being the data collected in the real world. The problem
we refer to as likelihood-free hypothesis testing is the task of deciding
between the hypotheses

$H_{0}:P_{Z}=P_{X}\quad\text{versus}\quad H_{1}:P_{Z}=P_{Y}.$
($\mathsf{LFHT}$)

This problem originates in [22; 49], where authors study the exponents of
error decay for finite $\mathcal{X}$ and fixed $P_{X},P_{Y}$ as $n\sim
m\to\infty$; more recently [30; 2; 25; 18] it is studied in the non-asymptotic
regime. Assuming that $P_{X},P_{Y}$ belong to a known (non-parametric) class
of distributions $\mathcal{P}$ and are guaranteed to be $\epsilon$-separated
with respect to total variation (TV) distance, [18] characterizes the sample
sizes $n$ and $m$ required for the sum of type-I and type-II errors to be
small, as a function of $\epsilon$ and for several different $\mathcal{P}$’s.
In particular, their results show that (for three settings of $\mathcal{P}$)
$(i)$ testing ($\mathsf{LFHT}$) is possible even when $n$ is not large enough
to accurately estimate $P_{X}$ and $P_{Y}$, and that $(ii)$ to achieve total
error $\alpha$ one can trade off $m$ vs. $n$ along a curve of the form
$\\{\min\\{n,\sqrt{mn}\\}\gtrsim
n_{\mathsf{TS}}(\alpha,\epsilon,\mathcal{P}),m\gtrsim\log(1/\alpha)/\epsilon^{2}\\}$,
where $n_{\mathsf{TS}}$ is the minimax sample complexity of two-sample testing
over $\mathcal{P}$. Here $\gtrsim$ suppresses dependence on constants and
untracked parameters. It is unclear, however, whether predictions drawn from
minimax sample complexities over specified distribution classes can be
observed in real-world data. Without the theory, a natural expectation is that
the error contour $\\{(m,n):\exists\text{ a test with total error
$\leq\alpha$}\\}$ would look similar to that of minimax two-sample testing
with unequal sample size, namely $\\{(m,n):\min\\{n,m\\}\gtrsim
n_{\mathsf{TS}}(\alpha,\epsilon,\mathcal{P})\\}$, i.e. $n$ and $m$ simply need
to be above a certain threshold simultaneously (as is the case for e.g. two-
sample testing over smooth densities [4; 34]). However, from Figure 1 we see
that there is indeed a non-trivial trade-off between $n$ and $m$: the contours
are not parallel to the axes and aren’t symmetric about the line $m=n$. The
importance of Fig. 1 is in demonstrating that said trade-off is not a kink of
a theory that arises due to some esoteric worst-case data distribution, but is
instead a real effect observed in state-of-the-art LFI algorithms ran on
actual data. We also remark that the $n$ used in this plot is the total number
of simulated samples (most of which are used for choosing a neural-network
parameterized kernel) and are not just the $n$ occuring in Theorems 3.2 and
3.3 which apply to a fixed kernel. See Section 4 for details on sample
division.

###  1.2 Mixed Likelihood-Free Hypothesis Testing

A prominent application of likelihood-free inference lies in the field of
particle physics. Scientists run sophisticated experiments in the hope of
finding a new particle or phenomenon. Often said phenomenon can be predicted
from theory, and thus can be simulated, as was the case for the Higgs boson
whose existence was verified after nearly $50$ years at the Large Hadron
Collider (LHC) [9; 3].

Suppose we have $n$ simulations from the _background_ distribution $P_{X}$ and
the _signal_ distribution $P_{Y}$. Further, we also have $m$ (real-world)
datapoints from $P_{Z}=(1-\nu)P_{X}+\nu P_{Y}$, i.e. the observed data is a
mixture between the background and signal distributions with rate parameter
$\nu$. The goal of physicists is to construct confidence intervals for $\nu$,
and a _discovery_ corresponds to a $5\sigma$ confidence interval that excludes
$\nu=0$. We model this problem by testing

$H_{0}:\nu=0\quad\text{versus}\quad H_{1}:\nu\geq\delta$ ($\mathsf{mLFHT}$)

for fixed (usually predicted) $\delta>0$. In particular, a discovery can be
claimed if $H_{0}$ is rejected. In the case of the Higgs boson discovery this
was done by first training a boosted decision tree classifier
$\mathcal{C}:\mathcal{X}\to[0,1]$ to distinguish background from signal
events, and then comparing $\frac{1}{m}\sum_{j}\mathcal{C}(Z_{j})$ to the
expected average under $P_{X}$ and $P_{Y}$ (estimated from simulated samples)
[9; 3]. Finally, the $p$-value and the associated significance is derived
using asymptotic distributions [11; 35].

##  2 The Likelihood-Free Test Statistic

###  2.1 Kernel Embeddings and MMD

Given a set $\mathcal{X}$, we call the function
$K:\mathcal{X}^{2}\to\mathbb{R}$ a kernel if the $n\times n$ matrix with
$ij$’th entry $K(x_{i},x_{j})$ is symmetric positive semidefinite for all
choices of $x_{1},\dots,x_{n}\in\mathcal{X}$ and $n\geq 1$. There is a unique
reproducing kernel Hilbert space (RKHS) $\mathcal{H}_{K}$ associated to $K$.
$\mathcal{H}_{K}$ consists of functions $\mathcal{X}\mapsto\mathbb{R}$ and
satisfies the reproducing property $\langle
K(x,\cdot),f\rangle_{\mathcal{H}_{K}}=f(x)$ for all $f\in\mathcal{H}_{k}$ and
$x\in\mathcal{X}$, in particular $K(x,\cdot)\in\mathcal{H}_{K}$. Given a
probability measure $P$ on $\mathcal{X}$, define its kernel embedding
$\mu_{P}$ as

$\mu_{P}\triangleq\operatorname{\mathbb{E}}_{X\sim
P}K(X,\cdot)=\int_{\mathcal{X}}K(x,\cdot)P(\mathrm{d}x).$ (1)

Given the kernel embeddings of two probability measures $P,Q$, we can measure
their distance in the RKHS by
$\operatorname{MMD}(P,Q)\triangleq\|\mu_{P}-\mu_{Q}\|_{\mathcal{H}_{K}}$,
where MMD stands for maximum mean discrepancy. MMD has a closed form thanks to
the reproducing property and linearity:

$\displaystyle\operatorname{MMD}^{2}(P,Q)$
$\displaystyle=\operatorname{\mathbb{E}}\Big{[}K(X,X^{\prime})+K(Y,Y^{\prime})-2K(X,Y)\Big{]},$

where $(X,X^{\prime},Y,Y^{\prime})\sim P^{\otimes 2}\otimes Q^{\otimes 2}$. In
particular, if $P,Q$ are empirical measures based on observations, we can
evaluate the MMD exactly, which is crucial in practice. Yet another attractive
property of MMD is that (under mild integrability conditions) it is an
integral probability metric (IPM) where the supremum is over the unit ball of
the RKHS $\mathcal{H}_{K}$. See e.g. [43; 38] for references.

Removing the Bias In our proofs we work with the kernel embedding of empirical
measures for which we need to modify the inner product
$\langle\cdot,\cdot\rangle_{\mathcal{H}_{K}}$ (and thus $\operatorname{MMD}$)
slightly by removing the diagonal terms. Namely, given i.i.d. samples $X,Y$ of
size $n,m$ respectively and corresponding empirical measures
$\widehat{P}_{X},\widehat{P}_{Y}$, we define

$\displaystyle\operatorname{MMD}_{u}^{2}(\widehat{P}_{X},\widehat{P}_{Y})\triangleq\sum_{i\neq
j}\frac{K(X_{i},X_{j})}{n(n-1)}+\sum_{i\neq
j}\frac{K(Y_{i},Y_{j})}{m(m-1)}-2\sum_{i,j}\frac{K(X_{i},Y_{j})}{mn}.$ (2)

We also write
$\langle\mu_{\widehat{P}_{X}},\mu_{\widehat{P}_{X}}\rangle_{u,\mathcal{H}_{K}}\triangleq\|\mu_{\widehat{P}_{X}}\|_{u,\mathcal{H}_{K}}^{2}\triangleq\frac{1}{n(n-1)}\sum_{i\neq
j}K(X_{i},X_{j})$ and extend linearly. The $u$ stands for unbiased, since
$\operatorname{\mathbb{E}}\operatorname{MMD}_{u}^{2}(\widehat{P}_{X},\widehat{P}_{Y})=\operatorname{MMD}^{2}(P_{X},P_{Y})\neq\operatorname{\mathbb{E}}\operatorname{MMD}^{2}(\widehat{P}_{X},\widehat{P}_{Y})$.

###  2.2 Test Statistic

Suppose that we have samples $X,Y,Z$ of sizes $n,n,m$ from the probability
measures $P_{X},P_{Y},P_{Z}$. Write $\widehat{P}_{X}$ for the empirical
measure of sample $X$, and analogously for $Y,Z$. The core of our test
statistic for ($\mathsf{mLFHT}$) is the following:

$\displaystyle T(X,Y,Z)$
$\displaystyle\triangleq\langle\mu_{\widehat{P}_{Z}},\mu_{\widehat{P}_{Y}}-\mu_{\widehat{P}_{X}}\rangle_{u,\mathcal{H}_{K}}=\frac{1}{nm}\sum_{i=1}^{n}\sum_{j=1}^{m}\Big{\\{}K(Z_{j},Y_{i})-K(Z_{j},X_{i})\Big{\\}}.$
(3)

Note that $T$ is of the additive form $\frac{1}{m}\sum_{j=1}^{m}f(Z_{i})$
where $f(z)\triangleq\mu_{\widehat{P}_{Y}}(z)-\mu_{\widehat{P}_{X}}(z)$ can be
interpreted as the _witness function_ of [20; 29]. Given some $\pi\in[0,1]$
(taken to be half the predicted signal rate ($\delta/2$) in our proofs), the
output of our test is

$\Psi_{\pi}=\mathbbm{1}\Big{\\{}T(X,Y,Z)\geq\gamma(X,Y,\pi)\Big{\\}},\text{
where
}\gamma(X,Y,\pi)=\pi\operatorname{MMD}^{2}_{u}(\widehat{P}_{X},\widehat{P}_{Y})+T(X,Y,X).$
(4)

The threshold $\gamma$ gives $\Psi_{\pi}$ a natural geometric interpretation:
it checks whether the projection of
$\mu_{\widehat{P}_{Z}}-\mu_{\widehat{P}_{X}}$ onto the vector
$\mu_{\widehat{P}_{Y}}-\mu_{\widehat{P}_{X}}$ falls further than $\pi$ along
the segment joining $\mu_{\widehat{P}_{X}}$ to $\mu_{\widehat{P}_{Y}}$ (up to
deviations due to the missing diagonal terms, see Section 2.1).

Note that when $\delta=1$, the problem ($\mathsf{mLFHT}$) reduces to
($\mathsf{LFHT}$), and $\Psi_{\delta/2}=\Psi_{1/2}=1$ if and only if
$\operatorname{MMD}_{u}(\widehat{P}_{Z},\widehat{P}_{X})-\operatorname{MMD}_{u}(\widehat{P}_{Z},\widehat{P}_{Y})\geq
0$. This statistic has been considered in the past for relative goodness-of-
fit testing [7], and when MMD is replaced by the $L^{2}$-distance we recover
the test statistic studied by [25; 30; 18].

Variance Cancellation It may seem more natural to the reader to threshold the
distance $\operatorname{MMD}_{u}(\widehat{P}_{Z},\widehat{P}_{X})$, resulting
in the test, say,
$\mathbbm{1}\\{\operatorname{MMD}_{u}(\widehat{P}_{Z},\widehat{P}_{X})\geq\operatorname{MMD}_{u}(\widehat{P}_{X},\widehat{P}_{Y})\delta/2\\}$.
The geometric interpretation of this would be similar to the one outlined in
Section 2.2. However, we note a crucial difference: ($\mathsf{LFHT}$) (the
case $\delta=1$) is possible with very little experimental data $m$ due to the
cancellation of variance. More precisely, the statistic
$\operatorname{MMD}^{2}(\widehat{P}_{Z},\widehat{P}_{X})$ contains the term
$\frac{1}{m(m-1)}\sum_{i\neq j}K(Z_{i},Z_{j})$ — whose variance is
prohibitively large and would inflate the $m$ required for reliable testing —
but this can be canceled by subtracting
$\operatorname{MMD}^{2}(\widehat{P}_{Z},\widehat{P}_{Y})$. Our statistic
$T(X,Y,Z)-\gamma(X,Y,\pi)$ simply generalizes this idea to ($\mathsf{mLFHT}$).

##  3 Minimax Rates of Testing

###  3.1 Technical Preliminaries

######  Theorem 3.1 (Hilbert–Schmidt theorem).

Suppose that $K\in L^{2}(\mu\otimes\mu)$ is symmetric. Then there exists a
sequence $(\lambda_{j})_{j\geq 1}\in\ell^{2}$ and an orthonormal basis
$\\{e_{j}\\}_{j\geq 1}$ of $L^{2}(\mu)$ such that $K(x,y)=\sum_{j\geq
1}\lambda_{j}e_{j}(x)e_{j}(y)$ for all $j\geq 1$, where convergence is in
$L^{2}(\mu\otimes\mu)$.

######  Assumption 1.

Throughout the rest of this section we _fix_ a choice of a non-negative
measure $\mu$ and kernel $K$ for which the conditions of Theorem 3.1 hold
(note that $\lambda_{j}\geq 0$ and depend on $\mu$).

Let $H_{0},H_{1}$ be two hypotheses, which we regard abstractly as sets of
probability measures. We say that a test $\Psi$ taking values in $\\{0,1\\}$
has total error probability less than $\alpha$ if

$\sup_{P\in H_{0}}P(\Psi=1)+\sup_{Q\in H_{1}}Q(\Psi=0)\leq\alpha.$ (5)

###  3.2 Upper Bounds on the Minimax Sample Complexity of ($\mathsf{mLFHT}$)

Suppose that we are given three i.i.d. samples $X,Y,Z$ from distributions
$P_{X},P_{Y},P_{Z}$ composed of $n,n,m$ observations respectively. Let
$\nu=\operatorname*{arg\,min}_{\nu^{\prime}\in\mathbb{R}}\operatorname{MMD}(P_{Z},(1-\nu^{\prime})P_{X}+\nu^{\prime}P_{Y})$
and let $\delta>0$. Assume the following:

  1. $(i)$

$P_{X},P_{Y}$ and $P_{Z}$ have $\mu$-densities bounded by $C$,

  2. $(ii)$

$\operatorname{MMD}(P_{X},P_{Y})\geq\epsilon$ holds, and

  3. $(iii)$

$\operatorname{MMD}(P_{Z},(1-\nu)P_{X}+\nu P_{Y})\leq
R\cdot\operatorname{MMD}(P_{X},P_{Y})$.

######  Theorem 3.2.

There exists a universal constant $c>0$ such that under assumptions
$(i$-$iii)$ above, $\Psi_{\delta/2}$ defined Section 2.2 tests
($\mathsf{mLFHT}$) at total error $\alpha$ provided

$\displaystyle\min\\{m,n\\}\geq
c\frac{C\|\lambda\|_{\infty}\log(1/\alpha)}{\left(\epsilon\delta/(1+R)\right)^{2}}\qquad\text{
and }\qquad\min\\{n,\sqrt{nm}\\}\geq
c\frac{C\|\lambda\|_{2}\log(1/\alpha)}{\delta\epsilon^{2}}.$

In an appropriate regime of the parameters, the sufficient sample complexity
in Theorem 3.2 exhibits a trade-off of the form
$\min\\{n,\sqrt{mn}\\}\gtrsim\|\lambda\|_{2}\log(1/\alpha)/(\delta\epsilon^{2})$
between the number of simulation samples $n$ and real observations $m$. The
existence of an asymmetric simulation-experimentation trade-off is confirmed
empirically in Figure 1, where we plot the sample sizes against the
probability of error. Note that we plot the _total_ number $n$ of simulation
samples, including those used for _training_ the kernel itself (see Section
4); which ensures that Figure 1 gives a realistic picture of data
requirements.

Note that Theorem 3.2 does _not_ place assumptions on the distributions
$P_{X},P_{Y}$ beyond bounded density, which is different from the usual
results in statistics where prior specification of distribution classes is
crucial and we consider it as a major advantage. In Section 3.5 we specialize
Theorem 3.2 to include such prior(s) and recover some known results.

###  3.3 Lower Bounds on the Minimax Sample Complexity of ($\mathsf{mLFHT}$)

Given $J\geq 2$, let
$\|\lambda\|_{2J}^{2}\triangleq\sum_{j=2}^{J}\lambda_{j}^{2}$. Let
$J^{\star}_{\epsilon}\triangleq\max\Big{\\{}J:\sup\limits_{\eta_{j}=\pm
1}\Big{\|}\sum_{j=2}^{J}\eta_{j}\sqrt{\lambda_{j}}e_{j}\Big{\|}_{\infty}\leq\frac{\|\lambda\|_{2J}}{2\epsilon}\Big{\\}}$.

######  Theorem 3.3 (Lower Bounds for mLFHT).

Suppose that $\int_{\mathcal{X}}K(x,y)\mu(\mathrm{d}x)\equiv\lambda_{1}$,
$\mu(\mathcal{X})=1$ and $\sup_{x\in\mathcal{X}}K(x,x)\leq 1$. There exists a
universal constant $c>0$ such that any test performing ($\mathsf{mLFHT}$) at
total error at most $\alpha$ must use a number $(n,m)$ of observations that
satisfy

$\displaystyle m\geq
c\frac{\log(1/\alpha)}{\epsilon^{2}\delta^{2}}\quad\text{and}\quad n\geq
c\frac{\|\lambda\|_{2J^{\star}_{\epsilon}}\sqrt{\log(1/\alpha)}}{\epsilon^{2}}\quad\text{and}\quad\delta
m+\sqrt{mn}$ $\displaystyle\geq
c\frac{\|\lambda\|_{2J^{\star}_{\epsilon}}\sqrt{\log(1/\alpha)}}{\epsilon^{2}\delta}.$

######  Remark 3.4.

Recall that the eigenvalues $\lambda$ depend on the choice of $\mu$, so that
by choosing a different base measure $\mu$ one can optimize the lower bound.

######  Remark 3.5.

The requirements $\sup_{x\in\mathcal{X}}K(x,x)\leq 1$ and $\mu(\mathcal{X})=1$
are is essentially without loss of generality, as $\mu$ and $K$ can be
rescaled. The condition
$\int_{\mathcal{X}}K(x,y)\mu(\mathrm{d}x)\equiv\lambda_{1}$ implies that the
top eigenfunction $e_{1}$ is equal to a constant or equivalently, that
$K(\mathrm{d}x,y)\mu(\mathrm{d}x)$ defines a Markov kernel (up to a
normalizing constant) for $\mu$-almost every $y\in\mathcal{X}$.

###  3.4 Tightness of Theorems 3.2 and 3.3

Dependence on $\boldsymbol{\|\lambda\|_{2}}$ An apparent weakness of Theorem
3.3 is its reliance on the unknown value $J^{\star}_{\epsilon}$, which depends
on the specifics of the kernel $K$ and base measure $\mu$. Determining it is
potentially highly nontrivial even for simple kernels. Slightly weakening
Theorem 3.3 we obtain the following corollary, which shows that the dependence
on $\|\lambda\|_{2}$ is tight, at least for small $\epsilon$.

######  Corollary 3.6.

Suppose $J\geq 2$ is such that $\sum_{j=2}^{J}\lambda_{j}^{2}\geq
c^{2}\|\lambda\|_{2}^{2}$ for some $c\leq 1$. Then
$\|\lambda\|_{2J^{\star}_{\epsilon}}$ can be replaced by $c\|\lambda\|_{2}$ in
Theorem 3.3 whenever $\epsilon\leq c\|\lambda\|_{2}/(2\sqrt{J-1})$.

Dependence on $\boldsymbol{R}$ and $\boldsymbol{\alpha}$ Due to the general
nature of our lower bound constructions, it is difficult to capture the
dependence on the misspecification parameter $R$. As for the probability of
error $\alpha$, based on related work [16] we expect the gap of size
$\sqrt{\log(1/\alpha)}$ to be a shortcoming of Theorem 3.2 and not the lower
bound. Closing this gap may require a different approach, however, as tests
based on empirical $L^{2}$ distances are known to have a sub-optimal
concentration in related problems [25].

Dependence on $\boldsymbol{\delta}$ The correct dependence on the signal rate
$\delta$ is the most important question left open by our theoretical results.
Any method requiring $n$ larger than a function of $\delta$ irrespective of
$m$ (as in Theorem 3.2) is provably sub-optimal because taking $m\gtrsim
1/(\delta\epsilon)^{2}$ and $n$ large enough to estimate both $P_{X},P_{Y}$ to
within accuracy $\epsilon/10$ always suffices to reach a fixed level of total
error.

###  3.5 Relation to Prior Results

In the following, we connect Theorem 3.2 to some special cases. Specifically,
we note that Theorem 3.2 recovers some known results in the literature [4; 34;
18] that are _minimax optimal_.

Binary Hypothesis Testing Suppose the two distributions $P_{X},P_{Y}$ are
known, we are given $m$ i.i.d. observations $Z_{1},\dots,Z_{m}\sim P_{Z}$ and
our task is to decide between the hypotheses $H_{0}:P_{X}=P_{Z}$ versus
$H_{1}:P_{Y}=P_{Z}$. Then, we may take $n=\infty,R=0,\delta=1$ in Theorem 3.2
to conclude that

$m\geq c\cdot\frac{C\|\lambda\|_{\infty}\log(1/\alpha)}{\epsilon^{2}}$

observations suffice to perform the test at total error $\alpha$.

Two-Sample Testing Suppose we have two i.i.d. samples $X$ and $Y$, both of
size $n$, from unknown distributions $P_{X},P_{Y}$ respectively and our task
is to decide between $H_{0}:P_{X}=P_{Y}$ against
$H_{1}:\operatorname{MMD}(P_{X},P_{Y})\geq\epsilon$. We split our $Y$ sample
in half resulting in $Y^{(1)}$ and $Y^{(2)}$ and form the statistic
$\Psi_{\mathsf{TS}}\triangleq\Psi_{1/2}(X,Y^{(1)},Y^{(2)})-\Psi_{1/2}(Y^{(1)},X,Y^{(2)})$,
where $\Psi_{1/2}$ is defined in Section 2.2. Then
$|\operatorname{\mathbb{E}}\Psi_{\mathsf{TS}}|$ is equal to $0$ under the null
hypothesis and is at least $1-2\alpha_{1}$ under the alternative, where
$\alpha_{1}$ is the target total error probability of $\Psi_{1/2}$. Taking
$\alpha_{1}=5\%$, by repeated sample splitting and majority voting we may
amplify the success probability to $\alpha$ provided

$n\geq c^{\prime}\frac{C\|\lambda\|_{2}\log(1/\alpha)}{\epsilon^{2}},$ (6)

where $c^{\prime}>0$ is universal (see Appendix for details). The upper bound
(6) partly recovers [34, Theorem 3 and 5] where authors show that thresholding
the MMD with Gaussian kernel
$G_{\sigma}(x,y)=\sigma^{-d}\exp(-\|x-y\|^{2}/\sigma^{2})$ achieves the
minimax optimal sample complexity $n\asymp\epsilon^{-(2\beta+d/2)/\beta}$ for
the problem of two-sample testing over the class $\mathcal{P}_{\beta,d}$ of
$d$-dimensional $(\beta,2)$-Sobolev-smooth distributions under
$\epsilon$-$L^{2}$-separation. For this, let $\sigma\asymp\epsilon^{1/\beta}$
which ensures that $\|P-Q\|_{L^{2}}\lesssim\operatorname{MMD}(P,Q)$ over
$P,Q\in\mathcal{P}_{\beta,d}$. Taking e.g.
$\mathrm{d}\mu(x)=\exp(-\|x\|_{2}^{2})\mathrm{d}x$ as the base measure, (6)
recovers the claimed sample complexity since $\|\lambda\|_{2}^{2}=\int
G^{2}_{\sigma}(x,y)\mathrm{d}\mu(x)\mathrm{d}\mu(y)=\mathcal{O}(\sigma^{-d})$
hiding dimension dependent constants. This result requires a bounded density
with respect to a Gaussian.

Likelihood-Free Hypothesis Testing By taking $\alpha\asymp
R\asymp\delta=\Theta(1)$ in Theorem 3.2 we can recover many of the results in
[18, Theorem 1]. When $\mathcal{X}$ is finite, we can take the kernel
$K(x,y)=\sum_{z\in\mathcal{X}}\mathbbm{1}\\{x=y=z\\}$ in Theorem 3.2 to obtain
the results for regular discrete distributions which state that under
$\epsilon$-$L^{1}$-separation the minimax optimal sample complexity is given
by $m\gtrsim
1/\epsilon^{2};\min\\{n,\sqrt{nm}\\}\gtrsim\sqrt{|\mathcal{X}|}/\epsilon^{2}$.
A similar kernel recovers the optimal result for the class of $\beta$-Hölder
smooth densities on the hypercube $[0,1]^{d}$. See Appendix for details.

Curse of Dimensionality Using the Gaussian kernel $G_{\sigma}$ as for two-
sample testing above, one can conclude by Theorem 3.2 that the required number
of samples for ($\mathsf{mLFHT}$) over the class $\mathcal{P}_{\beta,d}$ under
$\epsilon$-$L^{2}$-separation grows at most like
$\Omega\left((\frac{c}{\epsilon})^{\Omega(d)}\frac{1}{\delta^{2}}\right)$ for
some $c>1$, instead of the expected
$\Omega\left((\frac{c}{\epsilon\delta})^{\Omega(d)}\right)$. This may be
interpreted as theoretical support for the success of LFI in practice where
signal and background can be rather different (cf. [5, Figures 2-3]) and the
difficulty of the problem stems from the rate of signal events being small
(i.e. $\epsilon\approx 1$ but $\delta\ll 1$).

##  4 Learning Kernels from Data

Given a _fixed_ kernel $K$, our Theorems 3.2 and 3.3 show that the sample
complexity is heavily dependent on the separation $\epsilon$ under the given
$\operatorname{MMD}$ as well as the spectrum $\lambda=\lambda(\mu,K)$ of the
kernel. Thus, to have good test performance we need to use a kernel $K$ that
is well-adapted to the problem at hand. In practice, however, instead of using
a fixed kernel it would be only natural to use part of the simulation sample
to try to _learn_ a good kernel.

In Sections 4 and 5 we report experimental results after _training_ a kernel
parameterized by a neural network on part of the simulation data. In
particular, due to the dependence between the data and the kernel, Theorem 3.2
doesn’t directly apply. Our main contribution here is confirming the existence
of an asymmetric simulation-experimentation trade-off (cf. Figure 1 and also
Section 1.1). Figure 1 plots the _total_ number $n$ of simulations used,
including those used for training, so as to provide a realistic view of the
amount of data used. The experiments also illustrate that the
(trained-)kernel-based statistic of Section 2.2 achieves state-of-the-art
performance.

###  4.1 Proposed Training Algorithm

Consider splitting the data into three parts:
$(X^{\mathsf{tr}},Y^{\mathsf{tr}})$ is used for training (optimizing) the
kernel; $(X^{\mathsf{ev}},Y^{\mathsf{ev}})$ is used to evaluate our test
statistic at test time; and $(X^{\mathsf{cal}},Y^{\mathsf{cal}})$ is used for
calibrating the distribution of the test statistic under the null hypothesis.
We write $n_{s}=|X^{s}|=|Y^{s}|$ for
$s\in\\{\mathsf{tr},\mathsf{ev},\mathsf{cal}\\}$. Given the training data
$X^{\mathsf{tr}},Y^{\mathsf{tr}}$ with empirical measures
$\widehat{P}_{X^{\mathsf{tr}}},\widehat{P}_{Y^{\mathsf{tr}}}$, we maximize the
objective in

$\displaystyle\widehat{J}(X^{\mathsf{tr}},Y^{\mathsf{tr}};K)=\frac{\operatorname{MMD}_{u}^{2}(\widehat{P}_{X^{\mathsf{tr}}},\widehat{P}_{Y^{\mathsf{tr}}};K)}{\widehat{\sigma}(X^{\mathsf{tr}},Y^{\mathsf{tr}};K)},$
(7)

which was introduced in [45]. Here $\widehat{\sigma}^{2}$ is an estimator of
the variance of
$\operatorname{MMD}^{2}_{u}(\widehat{P}_{X^{\mathsf{tr}}},\widehat{P}_{Y^{\mathsf{tr}}};K)$
which is defined in (2). Intuitively, the objective $J$ aims to separate
$P_{X}$ from $P_{Y}$ while keeping variance low. For a heuristic justification
of its use for ($\mathsf{mLFHT}$) see Appendix.

In Algorithm 1 we describe the training and testing procedure, which produces
unbiased $p$-values for ($\mathsf{mLFHT}$) when there is no misspecification
($R=0$ in Theorem 3.2). During training, we use the Adam optimizer [31] with
stochastic batches.

######  Proposition 4.1.

When there is no misspecification ($R=0$ in Theorem 3.2), Algorithm 1 outputs
an unbiased estimate of the $p$-value that is consistent as
$\min\\{n_{\mathsf{cal}},k\\}\to\infty$.

Algorithm 1 mLFHT with a learned deep kernel

Input: $(X^{\mathsf{tr}}$, $X^{\mathsf{ev}},X^{\mathsf{cal}})$,
$(Y^{\mathsf{tr}},Y^{\mathsf{ev}},Y^{\mathsf{cal}}$); parametrized kernel
$K_{\omega}$; hyperparameters and initialization.

# Phase 1: Kernel training (optimization) on $X^{\mathsf{tr}}$ and
$Y^{\mathsf{tr}}$.

$\omega\leftarrow\arg\max_{\omega}^{\textnormal{Optimizer}}\widehat{J}(X^{\mathsf{tr}},Y^{\mathsf{tr}};K_{\omega})$;
# maximize objective as defined in (7)

# Phase 2: Distributional calibration of test statistic (under null
hypothesis).

for $r=1,2,\dots,k$ do

$Z^{\mathsf{cal},r}\leftarrow$ sample $m$ points without replacement from
$X^{\mathsf{cal}}$;

$\mathit{T}_{r}\leftarrow\frac{1}{n_{\mathsf{ev}}m}\sum_{i,j}\Big{(}K_{\omega}(Z^{\mathsf{cal},r}_{i},Y^{\mathsf{ev}}_{j})-K_{\omega}(Z^{\mathsf{cal},r}_{i},X^{\mathsf{ev}}_{j})\Big{)}$;

end for

# Phase 3: Inference with input $Z$.

$\widehat{T}\leftarrow\frac{1}{n_{\mathsf{ev}}m}\sum_{i,j}\Big{(}K_{\omega}(Z_{i},Y^{\mathsf{ev}}_{j})-K_{\omega}(Z_{i},X^{\mathsf{ev}}_{j})\Big{)}$;

Output:  Estimated $p$-value:
$\frac{1}{k}\sum_{i=1}^{k}\mathbbm{1}\\{\widehat{T}<T_{i}\\}$.

Time complexity Algorithm 1 runs in three separate stages: training,
calibration, and inference. The first two take $O\left(\\#\text{epochs}\cdot
B^{2}+kn_{\mathsf{ev}}m\right)$ total time, where $B$ is the batch size,
whereas Phase 3 takes only $O(n_{\mathsf{ev}}m)$ time, which is generally much
faster especially if $n_{\mathsf{ev}}<\\!\\!<n_{\mathsf{tr}}$.

Sample usage Empirically, data splitting in Algorithm 1 can have non-trivial
effects on performance. Instead of training the kernel on only a fraction of
the data
($\\{X^{\mathsf{tr}},Y^{\mathsf{tr}}\\}\cap\\{X^{\mathsf{ev}},Y^{\mathsf{ev}}\\}=\varnothing$),
we discovered that taking
$\\{X^{\mathsf{ev}},Y^{\mathsf{ev}}\\}\subseteq\\{X^{\mathsf{tr}},Y^{\mathsf{tr}}\\}$
results in more efficient use of data. The stronger condition
$\\{X^{\mathsf{ev}},Y^{\mathsf{ev}}\\}=\\{X^{\mathsf{tr}},Y^{\mathsf{tr}}\\}$
can also be applied; we take $\subseteq$ to reduce time complexity. We do,
however, crucially require $X^{\mathsf{cal}},Y^{\mathsf{cal}}$ in Phase 2 to
be independently sampled (“held-out”) for consistent $p$-value estimation.
Finally, we remark also that splitting this way is only valid in the context
of Algorithm 1. For the test (4) using the data-dependent threshold $\gamma$,
one needs
$\\{X^{\mathsf{tr}},Y^{\mathsf{tr}}\\}\cap\\{X^{\mathsf{ev}},Y^{\mathsf{ev}}\\}=\varnothing$
to estimate $\gamma$.

###  4.2 Classifier-Based Tests and Other Benchmarks

Let $\phi:\mathcal{X}\to[0,1]$ be a classifier, assigning small values to
$P_{X}$ and high values to $P_{Y}$ by minimizing cross-entropy loss of a
classifier net. There are two natural test statistics based on $\phi$:

Scheffé’s Test. The first idea, attributed to Scheffé in folklore [15, Section
6], is to take the statistic
$T(Z)=\frac{1}{m}\sum_{i=1}^{m}\mathbbm{1}\\{\phi(Z_{i})>t\\}$ where $t$ is
some (learn-able) threshold.

Approximate Neyman-Pearson / Logit Methods. If $\phi$ is trained to
perfection, then $\phi(z)=\mathbb{P}(P_{X}|z)$ would be the likelihood and
$\phi(z)/(1-\phi(z))$ would equal precisely the likelihood ratio between
$P_{Y}$ and $P_{X}$ at $z$. This motivates the use of
$T(Z)=\frac{1}{m}\sum_{i=1}^{m}\log(\phi(Z_{i})/(1-\phi(Z_{i})))$. See also
[10].

Let us list the testing procedures that we benchmark against each other in our
experiments.

  1. 1.

MMD-M: The $\operatorname{MMD}$ statistic (3) using $K$ with the mixing
architecture

$\displaystyle
K(x,y)=[(1-\tau)G_{\sigma}(\varphi_{\omega}(x),\varphi_{\omega}(y))+\tau]\cdot
G_{\sigma_{0}}(x+\varphi_{\omega}^{\prime}(x),y+\varphi_{\omega}^{\prime}(y)).$

Here $G_{\sigma}$ is the Gaussian kernel with variance $\sigma^{2}$;
$\varphi_{\omega},\varphi^{\prime}_{\omega^{\prime}}$ are NN’s (with
parameters $\omega,\omega^{\prime}$), and
$\sigma,\sigma_{0},\tau,\omega,\omega^{\prime}$ are trained.

  2. 2.

MMD-G: The $\operatorname{MMD}$ statistic (3) using the Gaussian kernel
architecture $K(x,y)=G_{\sigma}(\varphi_{\omega}(x),\varphi_{\omega}(y))$
where $\varphi_{\omega}$ is the feature mapping parametrized by a trained
network and $\sigma$ is a trainable parameter.

  3. 3.

MMD-O: The $\operatorname{MMD}$ statistic (3) using the Gaussian Kernel
$K(x,y)=G_{\sigma}(x,y)$ with optimized bandwidth $\sigma$. First proposed in
[7; 36].

  4. 4.

UME: An interpretable model comparison algorithm proposed by [29], which
evaluates the kernel mean embedding on a chosen “witness set”.

  5. 5.

SCHE, LBI: Scheffé’s test and Logit Based Inference methods [10], based on a
binary classifier network $\phi$ trained via cross-entropy, introduced above.

  6. 6.

RFM: Recursive Feature Machines, a recently proposed kernel learning algorithm
by [41].

###  4.3 Additive Statistics and the Thresholding Trick

Given a function $f$ (usually obtained by training) and test data
$Z=(Z_{1},\dots,Z_{m})$, we call a test additive if its output is obtained by
thresholding $T_{f}(Z)\triangleq\frac{1}{m}\sum_{i=1}^{m}f(Z_{i})$. We point
out that all of MMD-M/G/O, SCHE, LBI, UME, RFM are of this form, see the
Appendix for further details. Similarly to [36], we observe that any such
statistic can be realized by our kernel-based approach.

######  Proposition 4.2.

The kernel-based statistic defined in (3) with the kernel $K(x,y)=f(x)f(y)$ is
equal to $T_{f}$ up to a multiplicative and additive constant independent of
$Z$.

Motivated by the Scheffé’s test, instead of directly thresholding the additive
statistic $T_{f}(Z)$, we found empirically that replacing $f$ by
$f_{t}(x)\triangleq\mathbbm{1}\\{f(x)>t\\}$ can yield improved power. We set
$t$ by maximizing an estimate of the significance under the null using a
normal approximation, i.e. by solving
$t_{\mathsf{opt}}\triangleq\operatorname*{arg\,max}_{t}\frac{T_{f_{t}}(Y^{\mathsf{opt}})-T_{f_{t}}(X^{\mathsf{opt}})}{\sqrt{T_{f_{t}}(X^{\mathsf{opt}})(1-T_{f_{t}}(X^{\mathsf{opt}}))}}$,
where $X^{\mathsf{opt}},Y^{\mathsf{opt}}$ satisfy
$\\{X^{\mathsf{opt}},Y^{\mathsf{opt}}\\}\cap(\\{X^{\mathsf{cal}},Y^{\mathsf{cal}}\\}\cup\\{X^{\mathsf{ev}},Y^{\mathsf{ev}}\\})=\varnothing$.
This trick improves the performance of our tester on the Higgs dataset in
Section 5.2 but not for the image detection problem in Section 5.1.

##  5 Experiments

Figure 2: Empirical performance on (8) for the CIFAR detection problem when
$n_{\mathsf{tr}}=1920$. Plots from left to right are as follows. (a) rejection
rate under the alternative if test rejects whenever the estimated $p$-value is
smaller than $5\%$; (b) expected $p$-value [42] under the alternative; (c) the
average of type-I and II error probabilities when thresholded at 0 (different
from (4), see Appendix); and (d) ROC curves for different $m$ using MMD-M and
Algorithm 1. Shaded area shows the standard deviation across $10$ independent
runs. Missing benchmarks (thresholded MMD, MMD-O, LBI, RFM) are weaker; see
Appendix for full plot.

###  5.1 Image Source Detection

Our first empirical study looks at the task of detecting whether images come
from the CIFAR-10 [33] dataset or a SOTA generative model (DDPM) [24; 48].
While source detection is on its own interesting, it turns out that detecting
whether a group of images comes from the generative model versus the real
dataset can be too “easy” (see experiments in [29]). Therefore, we consider a
_mixed alternative_ , where the alternative hypothesis is not simply the
generative model but CIFAR with planted DDPM images. Namely, our $n$ labeled
images come from the following distributions:

$P_{X}=\text{CIFAR},\quad\text{and}\quad
P_{Y}=\pi\cdot\text{DDPM}+(1-\pi)\cdot\text{CIFAR}.$ (8)

The goal is to test whether the $m$ unlabeled observations $Z$ have been
corrupted with $\pi$ or more fraction of DDPM images (versus uncorrupted
CIFAR); this corresponds to ($\mathsf{mLFHT}$) with $\delta=1$. Figure 2 shows
the performance of our approach and the benchmarks when $\pi=1/3$.

Network Architecture With a standard deep CNN, the difference is only at the
final layer: for the kernel-based tests it is a feature output; for
classifiers, we add an extra linear layer to logits.

We see from Figure 2 that our kernel-based test outperforms other benchmarks
at a fixed training set size $n_{tr}$. One potential cause is that
$\operatorname{MMD}$ has an “optimization” subroutine (which it solves in
closed form) as it is an IPM. This additional layer of optimization may lead
to better performance at small sample sizes. The thresholding trick does not
seem to improve power empirically. We omit several benchmarks from this figure
for graphic presentation and they do not exhibit good separating power; see
the Appendix for the complete results. The bottom plot of Figure 1 shows $m$
and $n$ on $\log$-scale against the total probability of error, exhibiting the
simulation-experimentation trade-off.

###  5.2 Higgs-Boson Discovery

Figure 3: Expected significance of discovery on a mixture of 1000 backgrounds
and 100 signals in the Higgs experiment. Shaded area shows the standard
deviation over $10$ independent runs. See Appendix for full plot including
missing benchmarks.

The 2012 announcement of the Higgs boson’s discovery by the ATLAS and CMS
experiments [1; 9] marked a significant milestone in physics. The statistical
problem inherent in the experiment is well-modeled by ($\mathsf{mLFHT}$),
using a signal rate $\delta$ predicted by theory and misspecification
parameter $R=0$ (as was assumed in the original discovery). We consider our
algorithm’s power against past studies in the physics literature [5] as
measured by the _significance of discovery_. We note an important distinction
from Algorithm 1 in this application.

Estimating the Significance In physics, the threshold for claiming a
“discovery” is usually at a significance of $5\sigma$, corresponding to a
$p$-value of $2.87\times 10^{-7}$. Approximately
$n_{\mathsf{cal}}\sim(2.87)^{-1}\times 10^{7}$ samples would be necessary for
Algorithm 1 to reach such a precision. Fortunately the distribution of the
test statistic is approximated by a Gaussian customarily. We adopt this
approach for our experiment hereby assuming that $m$ is large enough for the
CLT to apply. We use the “expected significance of discovery” as our metric
[5] which, for the additive statistic
$T_{f}=\frac{1}{m}\sum_{i=1}^{m}f(Z_{i})$, is given by
$\frac{\delta(T_{f}(Y^{\mathsf{cal}})-T_{f}(X^{\mathsf{cal}}))}{\sqrt{\widehat{\operatorname{var}}(f(X^{\mathsf{cal}}))/m}}$.
If the thresholding trick (Section 4.3) is applied we use the more precise
Binomial tail, in which case the significance is estimated by
$-\Phi^{-1}(\mathbb{P}(\textup{Bin}(m,T_{f_{t_{\mathsf{opt}}}}(X^{\mathsf{cal}}))\geq
T_{f_{t_{\mathsf{opt}}}}(Z)))$, where $\Phi$ is the standard normal CDF.

Architecture and Training The architecture is a 6-layer feedforward net
similar for all tests (kernel-based and classifiers) except for the last
layer. We leave further details to the Appendix.

As can be seen in Figure 3, Scheffé’s test and MMD-M with threshold
$t_{\mathsf{opt}}$ are the best methods, achieving similar performance as the
algorithm of [5]; reaching the significance level of $5\sigma$ on 2.6 million
simulated datapoints and a test sample made up of a mixture of 1000
backgrounds and 100 signals. The top plot of Figure 1 shows $m$ and $n$ on
$\log$-scale against the total probability of error through performing the
test (4), exhibiting the asymmetric simulation-experimentation trade-off.

##  6 Conclusion

In this paper, we introduced ($\mathsf{mLFHT}$) as a theoretical model of
real-world likelihood-free signal detection problems arising in science. We
proposed a kernel-based test statistic and analyzed its minimax sample
complexity, obtaining both upper (Theorem 3.2) and lower bounds (Theorem 3.3)
in terms of multiple problem parameters, and discussed their tightness
(Section 3.4) and connections to prior work (Section 3.5). On the empirical
side, we described a method for training a parametrized kernel and proposed a
consistent $p$-value estimate (Algorithm 1 and Proposition 4.1). We examined
the performance of our method in two experiments and found that parametrized
kernels achieve state-of-the-art performance compared to relevant benchmarks
from the literature. Moreover, we confirmed experimentally the existence of
the asymmetric simulation-experimentation trade-off (Figure 1) which is
suggested by minimax analysis. We defer further special cases of Theorem 3.2,
all relevant proofs and experimental details to the Appendix.

## Acknowledgments and Disclosure of Funding

YP was supported in part by the National Science Foundation under Grant No
CCF-2131115. Research was sponsored by the United States Air Force Research
Laboratory and the Department of the Air Force Artificial Intelligence
Accelerator and was accomplished under Cooperative Agreement Number
FA8750-19-2-1000. The views and conclusions contained in this document are
those of the authors and should not be interpreted as representing the
official policies, either expressed or implied, of the Department of the Air
Force or the U.S. Government. The U.S. Government is authorized to reproduce
and distribute reprints for Government purposes notwithstanding any copyright
notation herein.

## References

  * Aad et al. [2012] Aad, G., Abajyan, T., Abbott, B., Abdallah, J., Khalek, S. A., Abdelalim, A. A., Aben, R., Abi, B., Abolins, M., AbouZeid, O., et al.  Observation of a new particle in the search for the standard model Higgs boson with the ATLAS detector at the LHC.  _Physics Letters B_ , 716(1):1–29, 2012.
  * Acharya et al. [2012] Acharya, J., Das, H., Jafarpour, A., Orlitsky, A., Pan, S., and Suresh, A.  Competitive classification and closeness testing.  In _Conference on Learning Theory_ , pages 22–1. JMLR Workshop and Conference Proceedings, 2012.
  * Adam-Bourdarios et al. [2015] Adam-Bourdarios, C., Cowan, G., Germain, C., Guyon, I., Kégl, B., and Rousseau, D.  The Higgs boson machine learning challenge.  In _NIPS 2014 workshop on high-energy physics and machine learning_ , pages 19–55. PMLR, 2015.
  * Arias-Castro et al. [2018] Arias-Castro, E., Pelletier, B., and Saligrama, V.  Remember the curse of dimensionality: The case of goodness-of-fit testing in arbitrary dimension.  _Journal of Nonparametric Statistics_ , 30(2):448–471, 2018.
  * Baldi et al. [2014] Baldi, P., Sadowski, P., and Whiteson, D.  Searching for exotic particles in high-energy physics with deep learning.  _Nature communications_ , 5(1):1–9, 2014.
  * Beaumont [2019] Beaumont, M. A.  Approximate bayesian computation.  _Annual review of statistics and its application_ , 6:379–403, 2019.
  * Bounliphone et al. [2016] Bounliphone, W., Belilovsky, E., Blaschko, M. B., Antonoglou, I., and Gretton, A.  A test of relative similarity for model selection in generative models.  In Bengio, Y. and LeCun, Y., editors, _4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings_ , 2016.
  * Brehmer et al. [2020] Brehmer, J., Louppe, G., Pavez, J., and Cranmer, K.  Mining gold from implicit models to improve likelihood-free inference.  _Proceedings of the National Academy of Sciences_ , 117(10):5242–5249, 2020.
  * Chatrchyan et al. [2012] Chatrchyan, S., Khachatryan, V., Sirunyan, A. M., Tumasyan, A., Adam, W., Aguilo, E., Bergauer, T., Dragicevic, M., Erö, J., Fabjan, C., et al.  Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC.  _Physics Letters B_ , 716(1):30–61, 2012.
  * Cheng and Cloninger [2019] Cheng, X. and Cloninger, A.  Classification logit two-sample testing by neural networks.  _CoRR_ , abs/1909.11298, 2019.
  * Cowan et al. [2011] Cowan, G., Cranmer, K., Gross, E., and Vitells, O.  Asymptotic formulae for likelihood-based tests of new physics.  _The European Physical Journal C_ , 71(2):1–19, 2011.
  * Cranmer et al. [2020] Cranmer, K., Brehmer, J., and Louppe, G.  The frontier of simulation-based inference.  _Proceedings of the National Academy of Sciences_ , 117(48):30055–30062, 2020.
  * Csilléry et al. [2010] Csilléry, K., Blum, M. G., Gaggiotti, O. E., and François, O.  Approximate bayesian computation (abc) in practice.  _Trends in ecology & evolution_, 25(7):410–418, 2010.
  * Dalmasso et al. [2020] Dalmasso, N., Izbicki, R., and Lee, A.  Confidence sets and hypothesis testing in a likelihood-free inference setting.  In _International Conference on Machine Learning_ , pages 2323–2334. PMLR, 2020.
  * Devroye and Lugosi [2012] Devroye, L. and Lugosi, G.  _Combinatorial Methods in Density Estimation_.  Springer Series in Statistics. Springer New York, 2012.  ISBN 9781461301257.
  * Diakonikolas et al. [2021] Diakonikolas, I., Gouleakis, T., Kane, D. M., Peebles, J., and Price, E.  Optimal testing of discrete distributions with high probability.  In _Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing_ , pages 542–555, 2021.
  * Diggle and Gratton [1984] Diggle, P. J. and Gratton, R. J.  Monte carlo methods of inference for implicit statistical models.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 46(2):193–212, 1984.
  * Gerber and Polyanskiy [2022] Gerber, P. R. and Polyanskiy, Y.  Likelihood-free hypothesis testing.  _CoRR_ , abs/2211.01126, 2022.  doi: 10.48550/arXiv.2211.01126.
  * Gretton et al. [2009] Gretton, A., Fukumizu, K., Harchaoui, Z., and Sriperumbudur, B. K.  A fast, consistent kernel two-sample test.  _Advances in neural information processing systems_ , 22, 2009.
  * Gretton et al. [2012a] Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., and Smola, A.  A kernel two-sample test.  _The Journal of Machine Learning Research_ , 13(1):723–773, 2012a.
  * Gretton et al. [2012b] 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.  _Advances in neural information processing systems_ , 25, 2012b.
  * Gutman [1989] Gutman, M.  Asymptotically optimal classification for multiple tests with empirically observed statistics.  _IEEE Transactions on Information Theory_ , 35(2):401–408, 1989.
  * Gutmann et al. [2018] Gutmann, M. U., Dutta, R., Kaski, S., and Corander, J.  Likelihood-free inference via classification.  _Statistics and Computing_ , 28:411–425, 2018.
  * Ho et al. [2020] Ho, J., Jain, A., and Abbeel, P.  Denoising diffusion probabilistic models.  _Advances in Neural Information Processing Systems_ , 33:6840–6851, 2020.
  * Huang and Meyn [2012] Huang, D. and Meyn, S.  Classification with high-dimensional sparse samples.  In _2012 IEEE International Symposium on Information Theory Proceedings_ , pages 2586–2590. IEEE, 2012.
  * Ingster [1987] Ingster, Y. I.  Minimax testing of nonparametric hypotheses on a distribution density in the $l_{p}$ metrics.  _Theory of Probability & Its Applications_, 31(2):333–337, 1987.  doi: 10.1137/1131042.
  * Izbicki et al. [2014] Izbicki, R., Lee, A., and Schafer, C.  High-dimensional density ratio estimation with extensions to approximate likelihood computation.  In _Artificial intelligence and statistics_ , pages 420–429. PMLR, 2014.
  * Jiang et al. [2017] Jiang, B., Wu, T.-y., Zheng, C., and Wong, W. H.  Learning summary statistic for approximate bayesian computation via deep neural network.  _Statistica Sinica_ , pages 1595–1618, 2017.
  * Jitkrittum et al. [2018] Jitkrittum, W., Kanagawa, H., Sangkloy, P., Hays, J., Schölkopf, B., and Gretton, A.  Informative features for model comparison.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Kelly et al. [2012] Kelly, B. G., Wagner, A. B., Tularak, T., and Viswanath, P.  Classification of homogeneous data with large alphabets.  _IEEE transactions on information theory_ , 59(2):782–795, 2012.
  * Kingma and Ba [2015] Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  In Bengio, Y. and LeCun, Y., editors, _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.
  * Krizhevsky et al. [2009] Krizhevsky, A., Hinton, G., et al.  Learning multiple layers of features from tiny images.  2009\.
  * Krizhevsky et al. [2020] Krizhevsky, A., Nair, V., and Hinton, G.  The cifar-10 dataset (2014).  _Online: http://www. cs. toronto. edu/kriz/cifar. html_ , 55, 2020\.
  * Li and Yuan [2019] Li, T. and Yuan, M.  On the optimality of gaussian kernel based nonparametric tests against smooth alternatives.  _arXiv preprint arXiv:1909.03302_ , 2019.
  * Lista [2017] Lista, L.  _Statistical methods for data analysis in particle physics_ , volume 941.  Springer, 2017.
  * Liu et al. [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_ , pages 6316–6326. PMLR, 2020.
  * Lopez-Paz and Oquab [2017] Lopez-Paz, D. and Oquab, M.  Revisiting classifier two-sample tests.  In _International Conference on Learning Representations_ , 2017.
  * Muandet et al. [2017] Muandet, K., Fukumizu, K., Sriperumbudur, B., Schölkopf, B., et al.  Kernel mean embedding of distributions: A review and beyond.  _Foundations and Trends® in Machine Learning_ , 10(1-2):1–141, 2017.
  * Papamakarios et al. [2019] Papamakarios, G., Sterratt, D., and Murray, I.  Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 837–848. PMLR, 2019.
  * Polyanskiy and Wu [2023+] Polyanskiy, Y. and Wu, Y.  _Information Theory: From Coding to Learning_.  Cambridge University Press, 2023+.
  * Radhakrishnan et al. [2022] Radhakrishnan, A., Beaglehole, D., Pandit, P., and Belkin, M.  Feature learning in neural networks and kernel machines that recursively learn features.  _arXiv preprint arXiv:2212.13881_ , 2022.
  * Sackrowitz and Samuel-Cahn [1999] Sackrowitz, H. and Samuel-Cahn, E.  P values as random variables-expected p values.  _The American Statistician_ , 53(4):326–331, 1999\.  ISSN 00031305.
  * Schölkopf and Smola [2001] Schölkopf, B. and Smola, A. J.  _Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond_.  The MIT Press, 06 2001.  ISBN 9780262256933.  doi: 10.7551/mitpress/4175.001.0001.
  * Sisson et al. [2018] Sisson, S. A., Fan, Y., and Beaumont, M.  _Handbook of approximate Bayesian computation_.  CRC Press, 2018.
  * Sutherland et al. [2017] Sutherland, D. J., Tung, H., Strathmann, H., De, S., Ramdas, A., Smola, A. J., and Gretton, A.  Generative models and model criticism via optimized maximum mean discrepancy.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net, 2017.
  * Thomas et al. [2022] Thomas, O., Dutta, R., Corander, J., Kaski, S., and Gutmann, M. U.  Likelihood-free inference by ratio estimation.  _Bayesian Analysis_ , 17(1):1–31, 2022.
  * Tsybakov [2008] Tsybakov, A. B.  _Introduction to Nonparametric Estimation_.  Springer Publishing Company, Incorporated, 1st edition, 2008.  ISBN 0387790519.
  * von Platen et al. [2022] von Platen, P., Patil, S., Lozhkov, A., Cuenca, P., Lambert, N., Rasul, K., Davaadorj, M., and Wolf, T.  Diffusers: State-of-the-art diffusion models.  https://github.com/huggingface/diffusers, 2022.
  * Ziv [1988] Ziv, J.  On classification with empirically observed statistics and universal data compression.  _IEEE Transactions on Information Theory_ , 34(2):278–286, 1988.

###### Table of Contents

  1. 1 Likelihood-Free Inference
    1. 1.1 $\mathsf{LFHT}$ and the Simulation-Experimentation Trade-off
    2. 1.2 Mixed Likelihood-Free Hypothesis Testing
  2. 2 The Likelihood-Free Test Statistic
    1. 2.1 Kernel Embeddings and MMD
    2. 2.2 Test Statistic
  3. 3 Minimax Rates of Testing
    1. 3.1 Technical Preliminaries
    2. 3.2 Upper Bounds on the Minimax Sample Complexity of ($\mathsf{mLFHT}$)
    3. 3.3 Lower Bounds on the Minimax Sample Complexity of ($\mathsf{mLFHT}$)
    4. 3.4 Tightness of Theorems 3.2 and 3.3
    5. 3.5 Relation to Prior Results
  4. 4 Learning Kernels from Data
    1. 4.1 Proposed Training Algorithm
    2. 4.2 Classifier-Based Tests and Other Benchmarks
    3. 4.3 Additive Statistics and the Thresholding Trick
  5. 5 Experiments
    1. 5.1 Image Source Detection
    2. 5.2 Higgs-Boson Discovery
  6. 6 Conclusion
  7. A Notation
  8. B Applications of Theorem 3.2
    1. B.1 Bounded Discrete Distributions Under $L^{2}/L^{1}$-Separation
    2. B.2 $\beta$-Hölder Smooth Densities on $[0,1]^{d}$ Under $L^{2}/L^{1}$-Separation
    3. B.3 $(\beta,2)$-Sobolev Smooth Densities on $\mathbb{R}^{d}$ Under $L^{2}$-Separation
  9. C Black-box Boosting of Success Probability
  10. D  Proof of Theorem 3.2
    1. D.1 Notation and Technical Tools
    2. D.2 Mean and Variance Computation
  11. E  Proof of Theorem 3.3
    1. E.1 Information theoretic tools
    2. E.2 Constructing hard instances
      1. E.2.1 Lower Bound on $m$
      2. E.2.2 Lower Bound on $n$
      3. E.2.3 Lower Bound on $m\cdot n$
  12. F Proofs From Section 4
    1. F.1 Heuristic Justification of the Objective (7)
    2. F.2 Proof of Proposition 4.1
    3. F.3 Proof of Proposition 4.2
    4. F.4 Additive Test Statistics
  13. G Application: Diffusion Models vs CIFAR
    1. G.1 Dataset Details
    2. G.2 Experiment Setup and Benchmarks
    3. G.3 Sample Allocation
    4. G.4 Remarks on Results
  14. H Application: Higgs-Boson Detection
    1. H.1 Dataset Details
    2. H.2 Experiment Setup and Training Models
      1. H.2.1 Configuration and Model Architecture
      2. H.2.2 Training
    3. H.3 Evaluating the Performance
      1. H.3.1 Evaluating the p-Value with the Methodology of Algorithm 1
      2. H.3.2 Evaluating the Error of the Test (4)

##  Appendix A Notation

We use $A\gtrsim B,A\lesssim B,A\asymp B$ to denote $A=\Omega(B),B=\Omega(A)$
and $A=\Theta(B)$ respectively, where the hidden constants depend on untracked
parameters multiplicatively.111For example, the first equation in (9) means
that there exists a constant $c$ independent of $\alpha,k,\epsilon,\delta,R$,
such that $\min\\{m,n\\}\geq
c\frac{\log(1/\alpha)(1+R)^{2}}{k\epsilon^{2}\delta^{2}}$.

We write $\operatorname{TV},\operatorname{KL},\chi^{2}$ for total-variation,
KL-divergence and $\chi^{2}$-divergence, respectively. We write
$D(P_{Y|X}\|Q_{Y|X}|P_{X})=\operatorname{\mathbb{E}}_{X\sim
P_{X}}D(P_{Y|X}\|Q_{Y|X})$ as the _conditional divergence_ for any probability
measures $P,Q$ on two variables $X,Y$ and divergence
$D\in\\{\operatorname{TV},\operatorname{KL},\chi^{2}\\}$.

We write $\ell^{p}$ for the usual $\ell^{p}$ sequence space and $L^{p}$ for
the usual $L^{p}$ space with respect to the Lebesgue measure. Both the
$\ell^{p}$ norm and the $L^{p}$ norm are written as $\|\cdot\|_{p}$ if no
ambiguity arises.

For real numbers $a,b\in\mathbb{R}$ we also write $\max\\{a,b\\}$ as $a\vee b$
and $\min\\{a,b\\}$ as $a\wedge b$.

We use $\vec{1}_{d}$ to denote an $d$-dimensional all $1$’s vector.

For an integer $k\in\mathbb{Z}^{+}$, we write $[k]$ as a short notation for
the set $\\{1,2,\dots,k\\}$.

In the proofs of Theorem 3.2 and Theorem 3.3, we use
$\stackrel{{\scriptstyle!}}{{=}}$ for an equality that we are trying to prove.

##  Appendix B Applications of Theorem 3.2

Usually, minimax rates of testing are proven under separation assumptions
using more traditional measures of distance such as $L^{p}$, where
$p\in[1,\infty]$. In this section we show one example of how Theorem 3.2 can
be used to recover known results, and also obtain some novel results under
$L^{2}$-separation and $L^{1}$-separation.

###  B.1 Bounded Discrete Distributions Under $L^{2}/L^{1}$-Separation

Sample Complexity Upper Bounds Let $\mathcal{P}_{\mathsf{Db}}(k,C)$ be the set
of all discrete distributions $P$ supported on $[k]=\\{1,2,\dots,k\\}$
satisfying $\max_{1\leq i\leq k}p(i)\leq C/k$, where $p$ is the probability
mass function of $P$ (here $\sum_{i=1}^{k}p(k)=1$). For distributions
$P_{X},P_{Y},P_{Z}$ we shall write $p_{X},p_{Y},p_{Z}$ as their probability
mass functions, respectively.

Let us apply Theorem 3.2 with underlying space $\mathcal{X}=[k]$ and measure
$\mu=\frac{1}{k}\sum_{i=1}^{k}\delta_{i}$. Take the kernel
$K(x,y)=\mathbbm{1}\\{x=y\\}=\sum_{i=1}^{k}\mathbbm{1}\\{x=y=i\\}$, and note
that for any two distributions $P_{X},P_{Y}$ we have

$\operatorname{MMD}^{2}(P_{X},P_{Y})=\operatorname{\mathbb{E}}\Big{[}K(X,X^{\prime})+K(Y,Y^{\prime})-2K(X,Y)\Big{]}=\sum_{i}|p_{X}(i)-p_{Y}(i)|^{2}$

where $(X,X^{\prime},Y,Y^{\prime})\sim P_{X}^{\otimes 2}\otimes P_{Y}^{\otimes
2}$. So the corresponding MMD is the $\ell^{2}$-distance on probability mass
functions. Note also that
$K=\sum_{i=1}^{k}\frac{1}{k}\left(\sqrt{k}\mathbbm{1}\\{x=i\\}\right)\left(\sqrt{k}\mathbbm{1}\\{y=i\\}\right)$,
where $\left\\{\sqrt{k}\mathbbm{1}\\{x=i\\}\right\\}_{i=1}^{k}$ forms an
orthonormal basis of $L^{2}(\mu)$. So $K$ has only one nonzero eigenvalue,
namely

$\lambda_{1}=\lambda_{2}=\ldots=\lambda_{k}=1/k,$

of multiplicity $k$. Suppose that we observe samples $X,Y,Z$ of size $n,n,m$
from $P_{X},P_{Y},P_{Z}\in\mathcal{P}_{\mathsf{Db}}(k,C)$, where
$\operatorname{MMD}(P_{X},P_{Y})=\sqrt{\sum_{i}|p_{X}(i)-p_{Y}(i)|^{2}}\geq\epsilon$.
Plugging into Theorem 3.2 shows that:

######  Proposition B.1.

For any two $P_{X},P_{Y}\in\mathcal{P}_{\mathsf{Db}}(k,C)$, if the
$\ell^{2}$-distance between $p_{X},p_{Y}$ is at least $\epsilon$, then testing
($\mathsf{mLFHT}$) is possible at total error $\alpha$ using $n$ simulation
samples and $m$ real data samples provided that

$\displaystyle\min\\{m,n\\}$
$\displaystyle\gtrsim\frac{C\|\lambda\|_{\infty}\log(1/\alpha)(1+R)^{2}}{\delta^{2}\epsilon^{2}}\asymp\frac{\log(1/\alpha)(1+R)^{2}}{k\epsilon^{2}\delta^{2}},$
(9) $\displaystyle\min\\{n,\sqrt{mn}\\}$
$\displaystyle\gtrsim\frac{C\|\lambda\|_{2}\sqrt{\log(1/\alpha)}}{\epsilon^{2}\delta}\asymp\frac{\sqrt{\log(1/\alpha)}}{\sqrt{k}\epsilon^{2}\delta}.$

where $R$ is defined as in the assumption $(iii)$ of Item $(iii)$.

We can convert the above results to measure separation with respect to total
variation (recall
$\operatorname{TV}(p,q)=\frac{1}{2}\sum_{i}|p(i)-q(i)|=\frac{1}{2}\|p-q\|_{1}$)
using the AM-QM inequality
$\|p_{X}-p_{Y}\|_{1}\leq\sqrt{k}\|p_{X}-p_{Y}\|_{2}$. Then, taking
$R\asymp\alpha\asymp\delta=\Theta(1)$ recovers the minimax optimal results of
[18] for $\mathsf{LFHT}$ over the class $\mathcal{P}_{\mathsf{Db}}$. Note that
analogous results for two-sample testing follow from the above using the
reduction presented in Section 3.5.

Sample Complexity Lower Bounds Recall the definition of $J^{\star}_{\epsilon}$
and note that $\|\lambda\|_{2,J}^{2}=\frac{\min(J-1,k)}{k^{2}}$ for all $J\geq
2$. By Corollary 3.6 we see that $J^{\star}_{\epsilon}\gtrsim k$ as soon as
$\epsilon\lesssim 1/k$. Thus, for $\epsilon\lesssim 1/k$ the necessity of

$m\gtrsim\frac{\log(1/\alpha)}{k\epsilon^{2}\delta^{2}},\,n\gtrsim\frac{\sqrt{\log(1/\alpha)}}{\sqrt{k}\epsilon^{2}}\,\text{
and
}\,m+\sqrt{mn}\gtrsim\frac{\sqrt{\log(1/\alpha)}}{\sqrt{k}\epsilon^{2}\delta}$
(10)

follows by Theorem 3.3. Here it is crucial to note that when
$\delta=\Theta(1)$, we have

$\displaystyle
m+\sqrt{mn}\gtrsim\frac{\sqrt{\log(1/\alpha)}}{\sqrt{k}\epsilon^{2}}\text{ and
}n\gtrsim\frac{\sqrt{\log(1/\alpha)}}{\sqrt{k}\epsilon^{2}}\iff\sqrt{mn}\gtrsim\frac{\sqrt{\log(1/\alpha)}}{\sqrt{k}\epsilon^{2}}\text{
and }n\gtrsim\frac{\sqrt{\log(1/\alpha)}}{\sqrt{k}\epsilon^{2}}$

and hence the upper bound (9) meets with the lower bound (10) provided
$R\asymp\delta=\Theta(1)$. Once again, setting
$R\asymp\delta\asymp\alpha=\Theta(1)$ we the optimal lower bounds recovering
the results of [18] (in the regime $\epsilon\lesssim 1/k$). In short we can
also recover the following result for $\mathsf{LFHT}$.

######  Proposition B.2.

On the class $P_{\mathsf{Db}}(k,C)$, using $n$ simulation samples and $m$ real
data samples, if

$\displaystyle n\gtrsim\frac{1}{\sqrt{k}\epsilon^{2}},\quad
m\gtrsim\frac{1}{k\epsilon^{2}},\quad\sqrt{mn}\gtrsim\frac{1}{\sqrt{k}\epsilon^{2}},$
(11)

then for any two distributions $P_{X},P_{Y}\in\mathcal{P}_{\mathsf{Db}}(k,C)$
with $\|p_{X}-p_{Y}\|_{2}\geq\epsilon$, testing ($\mathsf{LFHT}$) is possible
with a total error of $1\%$. Conversely, to ensure the existence a test that
can test ($\mathsf{LFHT}$) with a total error of $1\%$ for any
$P_{X},P_{Y}\in\mathcal{P}_{\mathsf{Db}}(k,C)$ with
$\|p_{X}-p_{Y}\|_{2}\geq\epsilon$, the number of observations $(n,m)$ must
satisfy

$\displaystyle n\gtrsim\frac{1}{\sqrt{k}\epsilon^{2}},\quad
m\gtrsim\frac{1}{k\epsilon^{2}},\quad\sqrt{mn}\gtrsim\frac{1}{\sqrt{k}\epsilon^{2}}.$
(12)

The implied constants in (11) and (12) do not depend on $k$ and $\epsilon$,
but may differ.

###  B.2 $\beta$-Hölder Smooth Densities on $[0,1]^{d}$ Under
$L^{2}/L^{1}$-Separation

Sample Complexity Upper Bounds Let $\mathcal{P}_{\mathsf{H}}(\beta,d,C)$ be
the set of all distributions on $[0,1]^{d}$ with $\beta$-Hölder smooth
Lebesgue-density $p$ satisfying $\|p\|_{\mathcal{C}^{\beta}}\leq C$ for some
constant $C>1$, where

$\displaystyle\|p\|_{\mathcal{C}^{\beta}}\triangleq\max_{0\leq|\alpha|\leq\lceil\beta-1\rceil}\|f^{(\alpha)}\|_{\infty}+\sup_{x\neq
y\in[0,1]^{d},|\alpha|=\lceil\beta-1\rceil}\frac{|f^{(\alpha)}(x)-f^{(\alpha)}(y)|}{\|x-y\|^{\beta-\lceil\beta-1\rceil}_{2}},$

where $\lceil\beta-1\rceil$ is the largest integer strictly smaller than
$\beta$ and $|\alpha|=\sum_{i}\alpha_{i}$ is the norm of a multi-index
$\alpha\in\mathbb{N}^{d}$. Abusing notation, we also use
$\mathcal{P}_{\mathsf{H}}(\beta,d,C)$ to denote the set of all corresponding
density functions.

We take $K(x,y)=\sum_{j}\mathbbm{1}\\{x,y\in B_{j}\\}$, where
$\\{B_{j}\\}_{j\in[\kappa]^{d}}$ is the $j$’th cell of the regular grid of
size $\kappa^{d}$ on $[0,1]^{d}$, i.e.,
$B_{j}=[(j-\vec{1}_{d})/\kappa,j/\kappa]$ for $j\in[\kappa]^{d}$. Clearly
there are $\kappa^{d}$ nonzero eigenvalues, each equal to $1$. The following
approximation result is due to Ingster [26], see also [4, Lemma 7.2].

######  Lemma B.3.

Let $f,g\in\mathcal{P}_{\mathsf{H}}(\beta,d,C)$ with
$\|f-g\|_{2}\geq\epsilon$. Then, there exist constants $c,c^{\prime}$
independent of $\epsilon$ such that for any $\kappa\geq c\epsilon^{-1/\beta}$,

$\operatorname{MMD}(f,g)\geq c^{\prime}\|f-g\|_{2}.$

Now, suppose that we have samples $X,Y,Z$ of size $n,n,m$ from
$P_{X},P_{Y},P_{Z}\in\mathcal{P}_{\mathsf{H}}(\beta,d,C)$ with densities
$p_{X},p_{Y},p_{Z}$ such that $\|p_{X}-p_{Y}\|_{2}\geq\epsilon$. Then, Theorem
3.2 combined with Lemma B.3 and the choice $\kappa\asymp\epsilon^{-1/\beta}$
shows that

######  Proposition B.4.

Testing ($\mathsf{mLFHT}$) on $\mathcal{P}_{\mathsf{H}}(\beta,d,C)$ at total
error $\alpha$ using $n$ simulation and $m$ real data samples is possible
provided

$\displaystyle\min\\{m,n\\}$
$\displaystyle\gtrsim\frac{C\|\lambda\|_{\infty}\log(1/\alpha)(1+R)^{2}}{\delta^{2}\epsilon^{2}}\asymp\frac{\log(1/\alpha)(1+R)^{2}}{\delta^{2}\epsilon^{2}},$
$\displaystyle\min\\{n,\sqrt{nm}\\}$
$\displaystyle\gtrsim\frac{C\|\lambda\|_{2}\sqrt{\log(1/\alpha)}}{\epsilon^{2}\delta}\asymp\frac{\sqrt{\log(1/\alpha)}}{\epsilon^{(2\beta+d/2)/\beta}\delta},$

where $\epsilon$ is an $L^{2}$-distance lower bound between $P_{X},P_{Y}$ and
$R$ is defined as in the assumption $(iii)$ of Item $(iii)$.

Setting $R\asymp\alpha\asymp\delta=\Theta(1)$ recovers the optimal results of
[18] for the class $\mathcal{P}_{\mathsf{H}}$. Once again, identical results
under $L^{1}$ separation follow from Jensen’s inequality
$\|\cdot\|_{L^{1}([0,1]^{d})}\leq\|\cdot\|_{L^{2}([0,1]^{d})}$. Note that
analogous results for two-sample testing follow from the above using the
reduction presented in Section 3.5.

Sample Complexity Lower Bounds The kernel defined in the previous paragraph is
not suitable for constructing lower bounds over the class
$\mathcal{P}_{\mathsf{H}}$ because its eigenfunctions do not necessarily lie
in $\mathcal{P}_{\mathsf{H}}$. It would be possible to consider a different
kernel that is more adapted to this problem/class but we do not pursue this
here.

###  B.3 $(\beta,2)$-Sobolev Smooth Densities on $\mathbb{R}^{d}$ Under
$L^{2}$-Separation

Sample Complexity Upper Bounds Let $\mathcal{P}_{\mathsf{S}}(\beta,d,C)$ be
the class of distributions that are supported on $\mathbb{R}^{d}$ and whose
Lebesgue density $p$ satisfies $\|p\|_{\beta,2}\leq C$, where

$\|p\|_{\beta,2}\triangleq\left\|(1+\|\cdot\|)^{\beta}\mathcal{F}[p]\right\|_{2}$
(13)

and $\mathcal{F}$ denotes the Fourier transform. Again, abusing notation, we
write $\mathcal{P}_{\mathsf{S}}(\beta,d,C)$ both as the set of distributions
and the set of density functions.

We take the Gaussian kernel
$G_{\sigma}(x,y)=\sigma^{-d}\exp(-\|x-y\|_{2}^{2}/\sigma^{2})$ on
$\mathcal{X}=\mathbb{R}^{d}$ with base measure
$\mathrm{d}\mu(x)=\exp(-x^{2})\mathrm{d}x$. In [34] the authors showed that
the two-sample test that thresholds the Gaussian MMD with appropriately chosen
variance $\sigma^{2}$ achieves the minimax optimal sample complexity over
$\mathcal{P}_{\mathsf{S}}$, when separation is measured by $L^{2}$. A key
ingredient in their proof is the following inequality.

######  Lemma B.5 ([34, Lemma 5]).

Let $f,g\in\mathcal{P}_{\mathsf{S}}(\beta,d,C)$ with
$\|f-g\|_{2}\geq\epsilon$. Then, there exist constants $c,c^{\prime}$
independent of $\epsilon$ such that for any $\sigma\leq
c\,\epsilon^{1/\beta}$, we have

$\operatorname{MMD}(f,g)\geq c^{\prime}\|f-g\|_{2}.$

Now, suppose that we have samples $X,Y,Z$ of sizes $n,n,m$ from
$P_{X},P_{Y},P_{Z}\in\mathcal{P}_{\mathsf{S}}(\beta,d,C)$ for some constant
$C$ with densities $p_{X},p_{Y},p_{Y}$ satisfying
$\|p_{X}-p_{Y}\|_{2}\geq\epsilon$.

Note that the heat-semigroup is an $L^{2}$-contraction
($\|\lambda\|_{\infty}\leq 1$) and that

$\|\lambda\|_{2}^{2}=\int
G_{\sigma}(x,y)^{2}\mathrm{d}\mu(x)\mathrm{d}\mu(y)\asymp\sigma^{-d}$

up to constants depending on the dimension. Theorem 3.2 combined with Lemma
B.5 and a choice $\sigma\asymp\epsilon^{1/\beta}$ yields the following result.

######  Proposition B.6.

Testing ($\mathsf{mLFHT}$) over the class $\mathcal{P}_{\mathsf{S}}$ with
total error $\alpha$ is possible provided

$\displaystyle\min\\{m,n\\}$
$\displaystyle\gtrsim\frac{C\|\lambda\|_{\infty}\log(1/\alpha)(1+R)^{2}}{\delta^{2}\epsilon^{2}}\asymp\frac{\log(1/\alpha)(1+R)^{2}}{\delta^{2}\epsilon^{2}}$
$\displaystyle\min\\{n,\sqrt{nm}\\}$
$\displaystyle\gtrsim\frac{C\|\lambda\|_{2}\sqrt{\log(1/\alpha)}}{\epsilon^{2}\delta}\asymp\frac{\sqrt{\log(1/\alpha)}}{\epsilon^{(2\beta+d/2)/\beta}\delta},$

where $\epsilon$ is the lower bound on the $L^{2}$-distance between
$P_{X},P_{Y}$ and $R$ is defined as in the assumption $(iii)$ of Item $(iii)$.

Taking $R\asymp\delta\asymp\alpha=\Theta(1)$ above, we obtain new results for
$\mathsf{LFHT}$ and using the reduction from two-sample testing given in
Section 3.5 we partly recover [34, Theorem 5]. Only partly, because the above
requires bounded density with respect to our base measure
$\mathrm{d}\mu(x)=\exp(-x^{2})\mathrm{d}x$.

Sample Complexity Lower Bounds Note that our lower bound Theorem 3.3 doesn’t
apply because the top eigenfunction of the Gaussian kernel is not constant.
Once again, a more careful choice of the base measure (or kernel) might lead
to a more suitable argument for the lower bound. We leave such pursuit as
open.

##  Appendix C Black-box Boosting of Success Probability

In this section we briefly describe how upper bounds on the minimax sample
complexity in the constant error probability regime ($\alpha=\Theta(1)$) can
be used to obtain the dependence $\log(1/\alpha)$ in the small error
probability regime ($\alpha=o(1)$). We will argue abstractly in a way that
applies to the setting of Theorem 3.2.

Suppose that from some distributions $P_{1},P_{2},\dots,P_{k}$ we take samples
$X^{1},X^{2},\dots,X^{k}$ of size $n_{1},n_{2},\dots,n_{k}$ respectively and
are able to decide between two hypotheses $H_{0}$ and $H_{1}$ (fixed but
arbitrary) with total error probability at most $1/3$. Call this test as
$\Psi(X^{1},\dots,X^{k})\in\\{0,1\\}$, so that

$\mathbb{P}(\Psi(X^{1},\dots,X^{k})=0|H_{0})\geq
2/3\qquad\text{and}\qquad\mathbb{P}(\Psi(X^{1},\dots,X^{k})=1|H_{1})\geq 2/3.$

Now, to each an error of $o(1)$, instead, we take
$18n_{1}\log(2/\alpha),\dots,18n_{k}\log(2/\alpha)$ observations from $P_{1}$
through $P_{k}$, and split each sample into $18\log(2/\alpha)$ equal sized
batches $\\{X^{i,j}\\}_{i\in[k],j\in[18\log(2/\alpha)]}$. Here
$18\log(2/\alpha)$ is assumed to be an integer without loss of generality. The
split samples form $18\log(2/\alpha)$ independent binary random variables

$A_{j}\triangleq\Psi(X^{1,j},\dots,X^{k,j})$

for $j=1,2,\dots,18\log(2/\alpha)$. We claim that the majority voting test

$\Psi_{\alpha}(\\{X^{i,j}\\}_{i,j})=\begin{cases}1&\text{if }\bar{A}\geq
1/2\\\ 0&\text{otherwise}\end{cases}$

tests $H_{0}$ against $H_{1}$ with total probability of error at most
$\alpha$, where

$\bar{A}\triangleq\frac{1}{18\log(2/\alpha)}\sum_{j=1}^{18\log(2/\alpha)}A_{j}.$

Indeed, by Hoeffding’s inequality, we have

$\displaystyle\mathbb{P}\left(\left.\bar{A}\geq 1/2\right|H_{0}\right)$
$\displaystyle\leq\alpha/2$ $\displaystyle\mathbb{P}\left(\left.\bar{A}\leq
1/2\right|H_{1}\right)$ $\displaystyle\leq\alpha/2.$

Therefore, in the remainder of our upper bound proofs, we only focus on
achieving a constant probability of error ($\alpha=\Theta(1)$) as the
logarithmic dependence follows by the above.

######  Remark C.1.

As mentioned in the discussion succeeding Corollary 3.6, we do conjecture the
_tight_ dependence in the upper bound to be $\sqrt{\log(\alpha^{-1})}$ instead
of $\log(\alpha^{-1})$ shown by this method.

##  Appendix D  Proof of Theorem 3.2

###  D.1 Notation and Technical Tools

We use the expansion

$K(x,y)=\sum_{\ell}\lambda_{\ell}e_{\ell}(x)e_{\ell}(y)$

extensively, where $\lambda\triangleq(\lambda_{1},\lambda_{2},\dots)$ are
$K$’s eigenvalues (regarded as an integral operator on $L^{2}(\mu)$) in non-
increasing order and $e_{1},e_{2},\dots$ are the corresponding eigenfunctions
forming an orthonormal basis for $L^{2}(\mu)$, and convergence is to be
understood in $L^{2}(\mu)$. We use the notation
$\langle\,\cdot\,\rangle\triangleq\int\cdot\,\mathrm{d}\mu$. For all $u\in
L^{2}(\mu)$ we define

$u_{\ell}\triangleq\langle ue_{\ell}\rangle,\quad
u_{\ell{\ell^{\prime}}}\triangleq\langle
ue_{\ell}e_{\ell^{\prime}}\rangle,\quad\ell=1,2,\dots$

and consequently $u=\sum_{\ell}u_{\ell}e_{\ell}$. We also define

$K[u](\cdot)\triangleq\int
K(t,\cdot)u(t)\mu(\mathrm{d}t)=\sum_{\ell}\lambda_{\ell}u_{\ell}e_{\ell}(\cdot),$

where the second equality follows from the orthonormality of
$\\{e_{\ell}\\}_{\ell=1}^{\infty}.$ Note that the RKHS embedding satisfies
$\mu_{u}\triangleq\int K(x,\cdot)u(x)\mathrm{d}\mu(x)=K[u]$. Now, for $P_{X}$
we write

$x_{\ell}\triangleq(p_{X})_{\ell}=\langle p_{X}e_{\ell}\rangle,\quad
x_{\ell{\ell^{\prime}}}\triangleq(p_{X})_{\ell{\ell^{\prime}}}=\langle
p_{X}e_{\ell}e_{\ell^{\prime}}\rangle,\quad\ell,{\ell^{\prime}}=1,2,\dots$

where $p_{X}$ is the $\mu$-density of $P_{X}$. The similar notations also
apply to $P_{Y},P_{Z}$. The following identities will be very useful in our
proofs.

######  Lemma D.1.

For each identity below, let $f,g,h\in L^{2}(\mu)$ be such that the quantity
is well defined. Then,

$\displaystyle\|\mu_{f}\|_{\mathcal{H}_{K}}^{2}$
$\displaystyle=\sum_{\ell}\lambda_{\ell}f_{\ell}^{2}$ (14)
$\displaystyle\operatorname{MMD}^{2}(f,g)$
$\displaystyle=\sum_{\ell}\lambda_{\ell}(f_{\ell}-g_{\ell})^{2}$ (15)
$\displaystyle\|K[f]\|_{2}^{2}$
$\displaystyle=\sum_{\ell}\lambda_{\ell}^{2}f_{\ell}^{2}$ (16)
$\displaystyle\sum_{\ell}\lambda_{\ell}f_{\ell}g_{\ell}$
$\displaystyle=\langle fK[g]\rangle=\langle K[f]g\rangle$ (17)
$\displaystyle\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}h_{\ell{\ell^{\prime}}}f_{\ell}g_{\ell^{\prime}}$
$\displaystyle=\langle hK[f]K[g]\rangle$ (18)
$\displaystyle\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}g_{\ell{\ell^{\prime}}}f_{\ell{\ell^{\prime}}}$
$\displaystyle=\sum_{\ell}\lambda_{\ell}\langle fe_{\ell}K[ge_{\ell}]\rangle.$
(19)

Suppose that $f,g$ are probability densities with respect to $\mu$ that are
bounded by $C$. Then

$0\leq\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}g_{\ell{\ell^{\prime}}}f_{\ell{\ell^{\prime}}}\leq
C^{2}\|\lambda\|_{2}^{2}.$ (20)

###### Proof.

We prove each claim, starting with (14). Clearly

$\displaystyle\|\mu_{f}\|^{2}_{\mathcal{H}_{K}}$
$\displaystyle=\|K[f]\|^{2}_{\mathcal{H}_{K}}$ $\displaystyle=\left\|\int
K(x,\cdot)f(x)\mathrm{d}\mu(x)\right\|^{2}_{\mathcal{H}_{K}}$
$\displaystyle=\iint\langle
K(x,\cdot),K(y,\cdot)\rangle_{\mathcal{H}_{K}}f(x)f(y)\mathrm{d}\mu(x)\mathrm{d}\mu(y)$
$\displaystyle=\iint K(x,y)f(x)f(y)\mathrm{d}\mu(x)\mathrm{d}\mu(y)$
$\displaystyle=\sum_{\ell}\lambda_{\ell}f_{\ell}^{2}$

as required. The second claim (15) follows immediately from (14) by
definition. For (16) by orthogonality we have

$\displaystyle\|K[f]\|_{2}^{2}$
$\displaystyle=\|\sum_{\ell}\lambda_{\ell}f_{\ell}e_{\ell}\|_{2}^{2}$
$\displaystyle=\sum_{\ell}\lambda_{\ell}^{2}f_{\ell}^{2}.$

For (17) by the definition of $K[\cdot]$ we have

$\displaystyle\sum_{\ell}\lambda_{\ell}f_{\ell}g_{\ell}$
$\displaystyle=\left\langle\left(\sum_{\ell}\lambda_{\ell}f_{\ell}e_{\ell}\right)g\right\rangle$
$\displaystyle=\langle K[f]g\rangle.$

For (18) we can write

$\displaystyle\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}h_{\ell{\ell^{\prime}}}f_{\ell}g_{\ell^{\prime}}$
$\displaystyle=\sum_{\ell}\lambda_{\ell}f_{\ell}\left\langle\left(\sum_{\ell^{\prime}}\lambda_{\ell^{\prime}}g_{\ell^{\prime}}e_{\ell^{\prime}}\right)he_{\ell}\right\rangle$
$\displaystyle=\sum_{\ell}\lambda_{\ell}f_{\ell}\langle K[g]he_{\ell}\rangle$
$\displaystyle=\langle K[g]hK[f]\rangle.$

Finally, for (19) we have

$\displaystyle\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}f_{\ell{\ell^{\prime}}}g_{\ell{\ell^{\prime}}}$
$\displaystyle=\sum_{\ell}\lambda_{\ell}\left\langle\left(\sum_{\ell^{\prime}}\lambda_{\ell^{\prime}}g_{\ell{\ell^{\prime}}}e_{\ell^{\prime}}\right)fe_{\ell}\right\rangle$
$\displaystyle=\sum_{\ell}\lambda_{\ell}\langle K[ge_{\ell}]fe_{\ell}\rangle.$

Suppose now that $f,g$ are probability densities with respect to $\mu$ that
are bounded by $C>0$. Let $X,Y$ be independent random variables following the
densities $f,g$. Then

$\displaystyle\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}f_{\ell{\ell^{\prime}}}g_{\ell{\ell^{\prime}}}$
$\displaystyle=\operatorname{\mathbb{E}}\left[\left(\sum_{\ell}\lambda_{\ell}e_{\ell}(X)e_{\ell}(Y)\right)^{2}\right]$
$\displaystyle\leq
C^{2}\int_{\mathcal{X}}\int_{\mathcal{X}}\left(\sum_{\ell}\lambda_{\ell}e_{\ell}(x)e_{\ell}(y)\right)^{2}\mathrm{d}\mu(x)\mathrm{d}\mu(y)$
$\displaystyle=C^{2}\|\lambda\|_{2}^{2}$

as claimed, where we used that the $e_{\ell}$ are orthonormal. ∎

###  D.2 Mean and Variance Computation

We take $\pi=\delta/2$. Our statistic reads

$\displaystyle-T(X,Y,Z)+\gamma(X,Y,\pi)$
$\displaystyle=\langle\mu_{\widehat{P}_{Z}}-(\bar{\pi}\mu_{\widehat{P}_{X}}+\pi\mu_{\widehat{P}_{Y}}),\mu_{\widehat{P}_{X}}-\mu_{\widehat{P}_{Y}}\rangle_{u,\mathcal{H}_{K}}$
$\displaystyle=\frac{1}{nm}\underbrace{\sum_{ij}k(X_{i},Z_{j})}_{\mathsf{I}}-\frac{1}{nm}\underbrace{\sum_{ij}k(Y_{i},Z_{j})}_{\mathsf{II}}-\frac{2\bar{\pi}}{n(n-1)}\underbrace{\sum_{i<i^{\prime}}k(X_{i},X_{i^{\prime}})}_{\mathsf{III}}$
$\displaystyle\qquad+\frac{2\pi}{n(n-1)}\underbrace{\sum_{i<i^{\prime}}k(Y_{i},Y_{i^{\prime}})}_{\mathsf{IV}}+\frac{\bar{\pi}-\pi}{n^{2}}\underbrace{\sum_{ij}k(X_{i},Y_{j})}_{\mathsf{V}}.$

Recall that
$\nu=\operatorname*{arg\,min}_{\nu^{\prime}\in\mathbb{R}}\operatorname{MMD}(P_{Z},\bar{\nu}^{\prime}P_{X}+\nu^{\prime}P_{Y})$.
Let us write $z=\bar{\nu}x+\nu y+r$ for $1-\bar{\nu}=\nu$, where the residual
term is denoted as $r\in L^{2}(\mu)$. Let $\mu_{r}=\int
r(t)K(t,\cdot)\mu(\mathrm{d}t)$ be the mean embedding of $r$. Under both
hypotheses we assume that $\|\mu_{r}\|_{\mathcal{H}_{K}}\leq
R\cdot\operatorname{MMD}(P_{X},P_{Y})$, moreover
$\langle\mu_{r},\mu_{P_{Y}}-\mu_{P_{X}}\rangle_{\mathcal{H}_{K}}=0$ by the
definition of $\nu$. We look at each of the $5+\binom{5}{2}=15$ terms of the
variance separately.

$\displaystyle\operatorname{var}(\mathsf{I})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}\Big{\\{}n(n-1)m(z_{\ell{\ell^{\prime}}}-z_{\ell}z_{\ell^{\prime}})x_{\ell}x_{\ell^{\prime}}+nm(m-1)(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})z_{\ell}z_{\ell^{\prime}}$
$\displaystyle\qquad\qquad\qquad\qquad+nm(x_{\ell{\ell^{\prime}}}z_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}}z_{\ell}z_{\ell^{\prime}})\Big{\\}}$
$\displaystyle\operatorname{var}(\mathsf{II})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}\Big{\\{}n(n-1)m(z_{\ell{\ell^{\prime}}}-z_{\ell}z_{\ell^{\prime}})y_{\ell}y_{\ell^{\prime}}+nm(m-1)(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})z_{\ell}z_{\ell^{\prime}}$
$\displaystyle\qquad\qquad\qquad\qquad+nm(y_{\ell{\ell^{\prime}}}z_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}}z_{\ell}z_{\ell^{\prime}})\Big{\\}}$
$\displaystyle\operatorname{var}(\mathsf{III})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}\Big{\\{}\binom{n}{2}(x_{\ell{\ell^{\prime}}}^{2}-x_{\ell}^{2}x_{\ell^{\prime}}^{2})+(\binom{n}{2}^{2}-\binom{n}{2}-\binom{4}{2}\binom{n}{4})(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})x_{\ell}x_{\ell^{\prime}}\Big{\\}}$
$\displaystyle\operatorname{var}(\mathsf{IV})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}\Big{\\{}\binom{n}{2}(y_{\ell{\ell^{\prime}}}^{2}-y_{\ell}^{2}y_{\ell^{\prime}}^{2})+(\binom{n}{2}^{2}-\binom{n}{2}-\binom{4}{2}\binom{n}{4})(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})y_{\ell}y_{\ell^{\prime}}\Big{\\}}$
$\displaystyle\operatorname{var}(\mathsf{V})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}\Big{\\{}n^{2}(n-1)(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})x_{\ell}x_{\ell^{\prime}}+n^{2}(n-1)(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})y_{\ell}y_{\ell^{\prime}}$
$\displaystyle\qquad\qquad\qquad\qquad+n^{2}(x_{\ell{\ell^{\prime}}}y_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}}y_{\ell}y_{\ell^{\prime}})\Big{\\}}$

For the cross terms we obtain

$\displaystyle\operatorname{cov}(\mathsf{I},\mathsf{II})$
$\displaystyle=\sum_{\ell\ell^{\prime}}\lambda_{\ell}\lambda_{\ell^{\prime}}n^{2}m(z_{\ell{\ell^{\prime}}}-z_{\ell}z_{\ell^{\prime}})x_{\ell}y_{\ell^{\prime}}$
$\displaystyle\operatorname{cov}(\mathsf{I},\mathsf{III})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}n(n-1)m(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})z_{\ell}x_{\ell^{\prime}}$
$\displaystyle\operatorname{cov}(\mathsf{I},\mathsf{IV})$ $\displaystyle=0$
$\displaystyle\operatorname{cov}(\mathsf{I},\mathsf{V})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}n^{2}m(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})z_{\ell}y_{\ell^{\prime}}$
$\displaystyle\operatorname{cov}(\mathsf{II},\mathsf{III})$ $\displaystyle=0$
$\displaystyle\operatorname{cov}(\mathsf{II},\mathsf{IV})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}n(n-1)m(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})z_{\ell}y_{\ell^{\prime}}$
$\displaystyle\operatorname{cov}(\mathsf{II},\mathsf{V})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}n^{2}m(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})z_{\ell}x_{\ell^{\prime}}$
$\displaystyle\operatorname{cov}(\mathsf{III},\mathsf{IV})$ $\displaystyle=0$
$\displaystyle\operatorname{cov}(\mathsf{III},\mathsf{V})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}n^{2}(n-1)(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})x_{\ell}y_{\ell^{\prime}}$
$\displaystyle\operatorname{cov}(\mathsf{IV},\mathsf{V})$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}n^{2}(n-1)(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})y_{\ell}x_{\ell^{\prime}}.$

Note that
$\binom{n}{2}^{2}-\binom{n}{2}-\binom{n}{2}\binom{n}{4}=n(n-1)^{2}-n(n-1)$.
Collecting terms, and simplifying, we get the coefficient of the $\frac{1}{n}$
term:

$\displaystyle\operatorname{Coef}\left(\frac{1}{n}\right)$
$\displaystyle=\sum_{\ell,{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}\Bigg{(}\underbrace{(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})z_{\ell}z_{\ell^{\prime}}}_{\operatorname{var}(\mathsf{I})}+\underbrace{(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})z_{\ell}z_{\ell^{\prime}}}_{\operatorname{var}(\mathsf{II})}+\underbrace{4\bar{\pi}^{2}(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})x_{\ell}x_{\ell^{\prime}}}_{\operatorname{var}(\mathsf{III})}$
$\displaystyle\qquad+\underbrace{4\pi^{2}(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})y_{\ell}y_{\ell^{\prime}}}_{\operatorname{var}(\mathsf{IV})}+\underbrace{(\bar{\pi}-\pi)^{2}(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})x_{\ell}x_{\ell^{\prime}}+(\bar{\pi}-\pi)^{2}(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})y_{\ell}y_{\ell^{\prime}}}_{\operatorname{var}(\mathsf{V})}$
$\displaystyle\qquad-\underbrace{4\bar{\pi}(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})z_{\ell}x_{\ell^{\prime}}}_{\operatorname{cov}(\mathsf{I},\mathsf{III})}+\underbrace{2(\bar{\pi}-\pi)(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})z_{\ell}y_{\ell^{\prime}}}_{\operatorname{cov}(\mathsf{I},\mathsf{V})}$
$\displaystyle\qquad-\underbrace{4\pi(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})z_{\ell}y_{\ell^{\prime}}}_{\operatorname{cov}(\mathsf{II},\mathsf{IV})}-\underbrace{2(\bar{\pi}-\pi)(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})z_{\ell}x_{\ell^{\prime}}}_{\operatorname{cov}(\mathsf{II},\mathsf{V})}$
$\displaystyle\qquad-\underbrace{4\bar{\pi}(\bar{\pi}-\pi)(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})x_{\ell}y_{\ell^{\prime}}}_{\operatorname{cov}(\mathsf{III},\mathsf{V})}+\underbrace{4\pi(\bar{\pi}-\pi)(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})y_{\ell}x_{\ell^{\prime}}}_{\operatorname{cov}(\mathsf{IV},\mathsf{V})}\Bigg{)}.$

After expanding $z_{\ell}$ as $z_{\ell}=\bar{\nu}x_{\ell}+\nu
y_{\ell}+r_{\ell}$, we split the calculation into multiple parts to simplify
it. First, we focus on terms that are multiplied by
$(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})$ and do not contain
$r_{\ell}$ or $r_{\ell^{\prime}}$. Using Lemma D.1 extensively and the fact
that $\bar{\pi}=1-\pi,\bar{\nu}=1-\nu$, we find that the sum of these terms
equals

$\displaystyle\bar{\nu}^{2}\langle xK[x]^{2}\rangle+\nu^{2}\langle
xK[y]^{2}\rangle+2\bar{\nu}\nu\langle xK[x]K[y]\rangle-\bar{\nu}^{2}\langle
xK[x]\rangle^{2}-\nu^{2}\langle xK[y]\rangle^{2}-2\bar{\nu}\nu\langle
xK[x]\rangle\langle xK[y]\rangle$ $\displaystyle+4\bar{\pi}^{2}\langle
xK[x]^{2}\rangle-4\bar{\pi}^{2}\langle
xK[x]\rangle^{2}+(\bar{\pi}-\pi)^{2}\langle
xK[y]^{2}\rangle-(\bar{\pi}-\pi)^{2}\langle xK[y]\rangle^{2}$
$\displaystyle-4\bar{\pi}\bar{\nu}\langle
xK[x]^{2}\rangle-4\bar{\pi}\nu\langle
xK[x]K[y]\rangle+4\bar{\pi}\bar{\nu}\langle
xK[x]\rangle^{2}+4\bar{\pi}\nu\langle xK[x]\rangle\langle xK[y]\rangle$
$\displaystyle+2(\bar{\pi}-\pi)\bar{\nu}\langle
xK[x]K[y]\rangle+2(\bar{\pi}-\pi)\nu\langle
xK[y]^{2}\rangle-2(\bar{\pi}-\pi)\bar{\nu}\langle xK[x]\rangle\langle
xK[y]\rangle-2(\bar{\pi}-\pi)\nu\langle xK[y]\rangle^{2}$
$\displaystyle-4\bar{\pi}(\bar{\pi}-\pi)\langle
xK[x]K[y]\rangle+4\bar{\pi}(\bar{\pi}-\pi)\langle xK[x]\rangle\langle
xK[y]\rangle$ $\displaystyle=$
$\displaystyle(\bar{\nu}-2\bar{\pi})^{2}\Big{(}\langle
xK[x-y]^{2}\rangle-\langle xK[x-y]\rangle^{2}\Big{)}$ $\displaystyle\leq$
$\displaystyle\,C\,\|\lambda\|_{\infty}\operatorname{MMD}^{2}(P_{X},P_{Y}).$

Similarly, the terms involving
$(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})$ but not $r_{\ell}$ or
$r_{\ell^{\prime}}$ sum up to the quantity

$\displaystyle(\nu-2\pi)^{2}\Big{(}\langle yK[x-y]^{2}\rangle-\langle
yK[x-y]\rangle^{2}\Big{)}$ $\displaystyle\leq
C\|\lambda\|_{\infty}\operatorname{MMD}^{2}(P_{X},P_{Y}).$

Next, collecting the terms involving both
$(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})$ and $r_{\ell}$ or
$r_{\ell^{\prime}}$ we get

$\displaystyle 2\bar{\nu}\langle xK[r]K[x]\rangle+2\nu\langle
xK[r]K[y]\rangle+\langle xK[r]^{2}\rangle-2\bar{\nu}\langle
xK[x]\rangle\langle xK[r]\rangle-2\nu\langle xK[y]\rangle\langle
xK[r]\rangle-\langle xK[r]\rangle^{2}$ $\displaystyle-4\bar{\pi}\langle
xK[x]K[r]\rangle+4\bar{\pi}\langle xK[x]\rangle\langle xK[r]\rangle$
$\displaystyle+2(\bar{\pi}-\pi)\langle
xK[y]K[r]\rangle-2(\bar{\pi}-\pi)\langle xK[y]\rangle\langle xK[r]\rangle$
$\displaystyle=$ $\displaystyle\,2(\bar{\nu}-2\bar{\pi})\Big{(}\langle
xK[r]K[x-y]\rangle-\langle xK[r]\rangle\langle xK[x-y]\rangle\Big{)}+\langle
xK[r]^{2}\rangle-\langle xK[r]\rangle^{2}$ $\displaystyle\lesssim$
$\displaystyle\,C\,\|\lambda\|_{\infty}(R+R^{2})\operatorname{MMD}^{2}(P_{X},P_{Y}).$

Finally, collecting the terms involving both
$(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})$ and $r_{\ell}$ or
$r_{\ell^{\prime}}$ we get

$\displaystyle 2(\nu-2\pi)\Big{(}\langle yK[r]K[y-x]\rangle-\langle
yK[r]\rangle\langle yK[y-x]\rangle\Big{)}+\langle yK[r]^{2}\rangle-\langle
yK[r]\rangle^{2}$ $\displaystyle\lesssim$
$\displaystyle\,C\|\lambda\|_{\infty}(R+R^{2})\operatorname{MMD}^{2}(P_{X},P_{Y}).$

Similarly we get

$\displaystyle\operatorname{Coef}\left(\frac{1}{m}\right)$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}\Bigg{(}\underbrace{(z_{\ell{\ell^{\prime}}}-z_{\ell}z_{\ell^{\prime}})x_{\ell}x_{\ell^{\prime}}}_{\operatorname{var}(\mathsf{I})}+\underbrace{(z_{\ell{\ell^{\prime}}}-z_{\ell}z_{\ell^{\prime}})y_{\ell}y_{\ell^{\prime}}}_{\operatorname{var}(\mathsf{I})}+\underbrace{2(z_{\ell{\ell^{\prime}}}-z_{\ell}z_{\ell^{\prime}})x_{\ell}y_{\ell^{\prime}}}_{\operatorname{cov}(\mathsf{I},\mathsf{II})}\Bigg{)}$
$\displaystyle=\langle zK[x-y]^{2}\rangle-\langle zK[x-y]\rangle^{2}$
$\displaystyle\lesssim
C\|\lambda\|_{\infty}\operatorname{MMD}^{2}(P_{X},P_{Y}).$

The remaining coefficients don’t rely on subtle cancellations, and simple
bounds yield

$\displaystyle\operatorname{Coef}\left(\frac{1}{n(n-1)}\right)$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}\Bigg{(}\underbrace{4\bar{\pi}^{2}\left(\frac{1}{2}(x_{\ell{\ell^{\prime}}}^{2}-x_{\ell}^{2}x_{\ell^{\prime}}^{2})-(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})x_{\ell}x_{\ell^{\prime}}\right)}_{\operatorname{var}(\mathsf{III})}$
$\displaystyle\qquad\qquad+\underbrace{4\pi^{2}\left(\frac{1}{2}(y_{\ell{\ell^{\prime}}}^{2}-y_{\ell}^{2}y_{\ell^{\prime}}^{2})-(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})y_{\ell}y_{\ell^{\prime}}\right)}_{\operatorname{var}(\mathsf{IV})}\Bigg{)}$
$\displaystyle\lesssim C^{2}\|\lambda\|_{2}^{2}$
$\displaystyle\operatorname{Coef}\left(\frac{1}{nm}\right)$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}\Bigg{(}\underbrace{-(z_{\ell{\ell^{\prime}}}-z_{\ell}z_{\ell^{\prime}})x_{\ell}x_{\ell^{\prime}}-(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})z_{\ell}z_{\ell^{\prime}}+(x_{\ell{\ell^{\prime}}}z_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}}z_{\ell}z_{\ell^{\prime}})}_{\operatorname{var}(\mathsf{I})}$
$\displaystyle\qquad-\underbrace{(z_{\ell{\ell^{\prime}}}-z_{\ell}z_{\ell^{\prime}})y_{\ell}y_{\ell^{\prime}}-(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})z_{\ell}z_{\ell^{\prime}}+(y_{\ell{\ell^{\prime}}}z_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}}z_{\ell}z_{\ell^{\prime}})}_{\operatorname{var}(\mathsf{I})}\Bigg{)}$
$\displaystyle\lesssim C^{2}\|\lambda\|_{2}^{2}$
$\displaystyle\operatorname{Coef}\left(\frac{1}{n^{2}}\right)$
$\displaystyle=\sum_{\ell{\ell^{\prime}}}\lambda_{\ell}\lambda_{\ell^{\prime}}\Bigg{(}\underbrace{(\bar{\pi}-\pi)\left(-(y_{\ell{\ell^{\prime}}}-y_{\ell}y_{\ell^{\prime}})x_{\ell}x_{\ell^{\prime}}-(x_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}})y_{\ell}y_{\ell^{\prime}}+(x_{\ell{\ell^{\prime}}}y_{\ell{\ell^{\prime}}}-x_{\ell}x_{\ell^{\prime}}y_{\ell}y_{\ell^{\prime}})\right)}_{\operatorname{var}(\mathsf{V})}\Bigg{)}$
$\displaystyle\lesssim C^{2}\|\lambda\|_{2}^{2}.$

Summarizing, we’ve found that

$\displaystyle\operatorname{var}(T(X,Y,Z)-\gamma(X,Y,\pi))$
$\displaystyle\lesssim\left(\frac{1}{n}+\frac{1}{m}\right)C\|\lambda\|_{\infty}(1+R^{2})\operatorname{MMD}^{2}(P_{X},P_{Y})$
(21)
$\displaystyle\qquad+\left(\frac{1}{n^{2}}+\frac{1}{nm}\right)C^{2}\|\lambda\|_{2}^{2}.$

Using that
$\langle\mu_{r},\mu_{P_{Y}}-\mu_{P_{X}}\rangle_{\mathcal{H}_{K}}=0$, we
compute the expectation to be

$\operatorname{\mathbb{E}}\left[-T(X,Y,Z)+\gamma(X,Y,\pi)\right]=(\pi-\nu)\operatorname{MMD}^{2}(P_{X},P_{Y}).$

Taking $\pi\triangleq\delta/2$ and applying Chebyshev’s inequality shows that
there exists a universal constant $c>0$, such that the testing problem is
possible at constant error probability (say $\alpha=5\%$), provided that the
sample sizes $m,n$ satisfy the following inequalities:

$\displaystyle\min\\{m,n\\}$ $\displaystyle\geq
c\frac{C\|\lambda\|_{\infty}(1+R^{2})}{\delta^{2}\epsilon^{2}}$
$\displaystyle\min\\{n,\sqrt{nm}\\}$ $\displaystyle\geq
c\frac{C\|\lambda\|_{2}}{\delta\epsilon^{2}}.$

By repeated sample splitting and majority voting (see Appendix C), we can
boost the success probability of this test to the desired level $1-\alpha$ by
incurring a multiplicative $\Theta(\log(1/\alpha))$ factor on the sample sizes
$n,m$, which yields the desired result.

##  Appendix E  Proof of Theorem 3.3

###  E.1 Information theoretic tools

Our lower bounds rely on the method of two fuzzy hypotheses [47]. Given a
measurable space $\mathcal{S}$, let $\mathcal{M}(\mathcal{S})$ denote the set
of all probability measures on $\mathcal{S}$. We call subsets
$H\subseteq\mathcal{M}(\mathcal{S})$ hypotheses. The following is the main
technical result that our proofs rely on.

######  Lemma E.1.

Take hypotheses $H_{0},H_{1}\subseteq\mathcal{M}(\mathcal{S})$ and
$P_{0},P_{1}\in\mathcal{M}(\mathcal{S})$ random with $\mathbb{P}(P_{i}\in
H_{i})=1$. Then

$\displaystyle\inf\limits_{\psi}\max_{i=0,1}\sup\limits_{P\in H_{i}}P(\psi\neq
i)\geq\frac{1}{2}\left(1-\operatorname{TV}(\mathbb{E}P_{0},\mathbb{E}P_{1})\right),$

where the infimum is over all tests $\psi:\mathcal{X}\to\\{0,1\\}$.

###### Proof.

For any $\psi$

$\displaystyle\max_{i=0,1}\sup\limits_{\mathbb{P}_{i}\in
H_{i}}\mathbb{P}_{i}(\psi\neq i)$
$\displaystyle\geq\frac{1}{2}\sup\limits_{\mathbb{P}_{i}\in
H_{i}}(\mathbb{P}_{0}(\psi=1)+\mathbb{P}_{1}(\psi=0))$
$\displaystyle\geq\frac{1}{2}\operatorname{\mathbb{E}}\Big{[}P_{0}(\psi=1)+P_{1}(\psi=0)\Big{]}.$

Optimizing over $\psi$ we get that the RHS above is equal to
$\frac{1}{2}(1-\operatorname{TV}(\operatorname{\mathbb{E}}P_{0},\operatorname{\mathbb{E}}P_{1}))$
as required. ∎

Therefore, to prove a lower bound on the minimax sample complexity of testing
with total error probability $\alpha$, we just need to construct two random
measures $P_{i}\in H_{i}$ such that
$1-\operatorname{TV}(\operatorname{\mathbb{E}}P_{0},\operatorname{\mathbb{E}}P_{1})=\Omega(\alpha)$.
In our proofs we also use the following standard results on $f$-divergences.

######  Lemma E.2 ([40, Section 7]).

For any probability distributions $P,Q$ the inequalities

$\displaystyle
1-\operatorname{TV}(P,Q)\geq\frac{1}{2}\exp(-\operatorname{KL}(P\|Q))\geq\frac{1}{2}\frac{1}{1+\chi^{2}(P\|Q)}$

hold.

######  Lemma E.3 (Chain rule for $\chi^{2}$-divergence).

Let $P_{X,Y},Q_{X,Y}$ be probability measures such that the marginals on $X$
are equal ($P_{X}=Q_{X}$). Then

$\chi^{2}(P_{X,Y}\|Q_{X,Y})=\chi^{2}(P_{Y|X}\|Q_{Y|X}|P_{X}).$

###### Proof.

Let $P_{X,Y},Q_{X,Y}$ have densities $p,q$ with respect to some $\mu$. Then,
by some abuse of notation, we have

$\displaystyle\chi^{2}(P_{X,Y}\|Q_{X,Y})$
$\displaystyle=-1+\int\frac{p(x,y)^{2}}{q(x,y)}\mathrm{d}\mu(x,y)$
$\displaystyle=-1+\int\frac{p(y|x)^{2}p(x)}{q(y|x)}\mathrm{d}\mu(x,y)$
$\displaystyle=\int
p(x)\int\left(\frac{p(y|x)^{2}}{q(y|x)}-1\right)\mathrm{d}\mu(y,x)$
$\displaystyle=\chi^{2}(P_{Y|X}\|Q_{Y|X}|P_{X}).$

∎

###  E.2 Constructing hard instances

Recall that in the statement of Theorem 3.3, we assume that
$\mu(\mathcal{X})=1$, $\sup_{x\in\mathcal{X}}K(x,x)\leq 1$ and $\int
K(x,y)\mu(\mathrm{d}x)\equiv\lambda_{1}$. Let $f_{0}\equiv 1$ and for each
$\eta\in\\{\pm 1\\}^{\mathbb{N}}$ define

$f_{\eta}=1+\epsilon\underbrace{\sum_{j\geq
2}\rho_{j}\eta_{j}e_{j}}_{\triangleq g_{\eta}}$ (22)

where $\\{\rho_{j}\\}_{j\geq 2}$ is chosen as $\rho_{j}=\mathbbm{1}\\{2\leq
j\leq J\\}\sqrt{\lambda_{j}}/\|\lambda\|_{2,J}$, where we define
$\|\lambda\|_{2,J}=\sqrt{\sum_{2\leq j\leq J}\lambda_{j}^{2}}$ for some $J\geq
2$. Notice that $\int f_{\eta}(x)\mu(\mathrm{d}x)=\mu(\mathcal{X})=1$ due to
orthogonality of the eigenfunctions. Assume from here on that $J$ is chosen so
that for all $\eta$ we have $f_{\eta}(x)\geq 1/2$ for all $x\in\mathcal{X}$.
This makes $f_{\eta}$ into a valid probability density with respect to the
base measure $\mu$. Before continuing, we prove the following Lemma, which
gives a lower bound on the maximal $J$ for which $f_{\eta}\geq 1/2$ for all
$\eta$.

######  Lemma E.4.

$J\leq J^{\star}_{\epsilon}$ holds provided
$2\epsilon\sqrt{J-1}\leq\|\lambda\|_{2J}$.

###### Proof of Lemma E.4.

Notice that

$\|e_{j}\|_{\infty}=\sup_{x\in\mathcal{X}}\langle
K(x,\cdot),e_{j}\rangle_{\mathcal{H}}\leq\sup_{x\in\mathcal{X}}\|K(x,\cdot)\|_{\mathcal{H}}\|e_{j}\|_{\mathcal{H}}\leq\frac{1}{\sqrt{\lambda}_{j}},$
(23)

where we use $\|K(x,\cdot)\|_{\mathcal{H}}=\sqrt{K(x,x)}$. We have

$\displaystyle\|g_{\eta}\|_{\infty}$ $\displaystyle=\epsilon\|\sum_{j\geq
2}\rho_{j}\eta_{j}e_{j}\|_{\infty}=\epsilon\sup\limits_{x\in\mathcal{X}}\langle
K(x,\cdot),\sum_{j\geq 2}\rho_{j}\eta_{j}e_{j}\rangle_{\mathcal{H}}$
$\displaystyle\leq\epsilon\|\sum_{j\geq
2}\rho_{j}\eta_{j}e_{j}\|_{\mathcal{H}}=\epsilon\sqrt{\sum_{j\geq
2}\rho_{j}^{2}/\lambda_{j}}=\frac{\epsilon\sqrt{J-1}}{\|\lambda\|_{2,J}},$

and the result follows. ∎

Note that Lemma E.4 immediately gives us a proof of Corollary 3.6.

###### Proof of Corollary 3.6.

Suppose that $J$ is such that $\sum_{j=2}^{J}\lambda_{j}^{2}\geq
c^{2}\|\lambda\|_{2}^{2}$. Then, by Lemma E.4, if
$\epsilon\leq\|\lambda\|_{2J}/(2\sqrt{J-1})$ then $J\leq
J_{\epsilon}^{\star}$. By assumption, this is implied by the inequality
$\epsilon\leq c\|\lambda\|_{2}/(2\sqrt{J-1})$, and the result follows. ∎

Continuing with our proof, note that by construction we have

$\operatorname{MMD}^{2}(f_{0},f_{\eta})=\sum_{j\geq
2}\lambda_{j}\rho_{j}^{2}=\epsilon^{2},\quad\forall\eta\in\\{\pm
1\\}^{\mathbb{N}}.$ (24)

####  E.2.1 Lower Bound on $m$

Again, we apply Lemma E.1 with the new (deterministic) construction

$P_{0}=f_{0}^{\otimes n}\otimes f^{\otimes
n}_{\mathbbm{1}}\otimes(1+\delta\epsilon g_{\mathbbm{1}})^{\otimes m}\qquad
P_{1}=f_{0}^{\otimes n}\otimes f^{\otimes n}_{\mathbbm{1}}\otimes
f_{0}^{\otimes m},$ (25)

where we write $f_{\mathbbm{1}}=f_{(1,1,\dots)}$ and similarly for
$g_{\mathbbm{1}}$. By the data-processing inequality for $\chi^{2}$-divergence
(also by Lemma E.3), we may drop the first $2n$ coordinates and obtain

$\displaystyle\chi^{2}(\operatorname{\mathbb{E}}P_{0},\operatorname{\mathbb{E}}P_{1})$
$\displaystyle=\chi^{2}((1+\delta\epsilon g_{\mathbbm{1}})^{\otimes
m}\|f_{0}^{\otimes m})$ $\displaystyle=(1+\delta^{2}\epsilon^{2})^{m}-1$
$\displaystyle\leq\exp(\delta^{2}\epsilon^{2}m)-1.$

By Lemma E.2 we

$1-\operatorname{TV}(\operatorname{\mathbb{E}}P_{0},\operatorname{\mathbb{E}}P_{1})\gtrsim\frac{1}{\chi^{2}(\operatorname{\mathbb{E}}P_{0},\operatorname{\mathbb{E}}P_{1})-1}\geq\exp(-\delta^{2}\epsilon^{2}m)\stackrel{{\scriptstyle!}}{{=}}\Omega(\alpha).$

The lower bound $m\gtrsim\log(1/\alpha)/(\delta\epsilon)^{2}$ now follows
readily.

####  E.2.2 Lower Bound on $n$

Once again, we apply Lemma E.1 to the new construction

$P_{0}=f_{0}^{\otimes n}\otimes f_{\eta}^{\otimes n}\otimes f_{0}^{\otimes
m},\qquad P_{1}=f_{\eta}^{\otimes n}\otimes f_{0}^{\otimes n}\otimes
f_{0}^{\otimes m},$ (26)

where we put a uniform prior on $\eta\in\\{\pm 1\\}^{\mathbb{N}}$ as before.
Using the subadditivity of total variation under products, we compute

$\displaystyle\operatorname{TV}(\operatorname{\mathbb{E}}P_{0},\operatorname{\mathbb{E}}P_{1})$
$\displaystyle=\operatorname{TV}(f_{0}^{\otimes
n}\otimes\operatorname{\mathbb{E}}f_{\eta}^{\otimes
n},\operatorname{\mathbb{E}}[f_{\eta}^{\otimes n}]\otimes f_{0}^{\otimes n})$
$\displaystyle\leq
2\operatorname{TV}(\operatorname{\mathbb{E}}f_{\eta}^{\otimes
n},f_{0}^{\otimes n}).$

Just as in Section E.2.3 we upper bound by the $\chi^{2}$-divergence to get

$\displaystyle\chi^{2}(\operatorname{\mathbb{E}}f_{\eta}^{\otimes
n}\|f_{0}^{\otimes n})$
$\displaystyle=-1+\operatorname{\mathbb{E}}_{\eta\eta^{\prime}}\int\prod\limits_{i=1}^{n}(f_{\eta}(x_{i})f_{\eta^{\prime}}(x_{i}))\mu(\mathrm{d}x_{1})\dots\mu(\mathrm{d}x_{n})$
$\displaystyle\leq-1+\operatorname{\mathbb{E}}\exp(n\epsilon^{2}\sum_{j\geq
2}\rho_{j}^{2}\eta_{j}\eta^{\prime}_{j})$
$\displaystyle=-1+\prod\limits_{j\geq 2}\cosh(n\epsilon^{2}\rho_{j}^{2})$
$\displaystyle\leq-1+\exp(n^{2}\epsilon^{4}\sum_{j\geq 2}\rho_{j}^{4})$
$\displaystyle=-1+\exp(n^{2}\epsilon^{4}/\|\lambda\|_{2,J}^{2}).$

Again, by Lemma E.2 we obtain

$1-\operatorname{TV}(\operatorname{\mathbb{E}}P_{0},\operatorname{\mathbb{E}}P_{1})\gtrsim\frac{1}{\chi^{2}(\operatorname{\mathbb{E}}P_{0}\|\operatorname{\mathbb{E}}P_{1})-1}\geq\exp(-n^{2}\epsilon^{4}/\|\lambda\|_{2,J}^{2})\stackrel{{\scriptstyle!}}{{=}}\Omega(\alpha).$

The lower bound $n\gtrsim\sqrt{\log(1/\alpha)}\|\lambda\|_{2,J}/\epsilon^{2}$
now follows readily.

####  E.2.3 Lower Bound on $m\cdot n$

We take a uniform prior on $\eta$ and consider the random measures

$P_{0}=f_{0}^{\otimes n}\otimes f_{\eta}^{\otimes
n}\otimes((1-\delta)f_{0}+\delta f_{\eta})^{\otimes m}\qquad\text{and}\qquad
P_{1}=f_{0}^{\otimes n}\otimes f_{\eta}^{\otimes n}\otimes f_{0}^{\otimes m}.$
(27)

Our goal is to apply Lemma E.1 to $P_{0},P_{1}$. Notice that
$(1-\delta)f_{0}+\delta f_{\eta}=1+\delta\epsilon g_{\eta}$. Let us write
$X,Y,Z$ for the marginals first $n$, second $n$ and last $m$ coordinates of
$P_{0}$ and $P_{1}$. By the data processing inequality and the chain rule
Lemma E.3 we have

$\displaystyle\chi^{2}(\operatorname{\mathbb{E}}P_{0}\|\operatorname{\mathbb{E}}P_{1})$
$\displaystyle=\chi^{2}((\operatorname{\mathbb{E}}P_{0})_{Y,Z}\|(\operatorname{\mathbb{E}}P_{1})_{Y,Z})$
$\displaystyle=\chi^{2}((\operatorname{\mathbb{E}}P_{0})_{Z|Y}\|(\operatorname{\mathbb{E}}P_{1})_{Z|Y}|(\operatorname{\mathbb{E}}P_{0})_{Y})$
$\displaystyle=\operatorname{\mathbb{E}}\chi^{2}\left(\operatorname{\mathbb{E}}\left[\left.(1+\delta\epsilon
g_{\eta})^{\otimes m}\right|Y\right]\|f_{0}^{\otimes m}\right)=:(\dagger).$

Notice that the expectation inside the $\chi^{2}$-divergence is with respect
to $\eta$ given the variables $Y$, or in other words, over the posterior of
$\eta$ with uniform prior given $n$ observations from the density $1+\epsilon
g_{\eta}=f_{\eta}$. The outer expectation is over $Y$. Given $Y$, let $\eta$
and $\eta^{\prime}$ be i.i.d. from said posterior. We get the bound

$\displaystyle(\dagger)+1$
$\displaystyle\leq\operatorname{\mathbb{E}}\int\prod\limits_{i=1}^{m}(1+\delta\epsilon
g_{\eta}(x_{i}))(1+\delta\epsilon
g_{\eta^{\prime}}(x_{i}))\mu(\mathrm{d}x_{i})$
$\displaystyle=\operatorname{\mathbb{E}}(1+\delta^{2}\epsilon^{2}\sum_{j\geq
2}\rho_{j}^{2}\eta_{j}\eta^{\prime}_{j})^{m}$
$\displaystyle\leq\operatorname{\mathbb{E}}\exp(\delta^{2}\epsilon^{2}m\sum_{j\geq
2}\rho_{j}^{2}\eta_{j}\eta^{\prime}_{j}).$

Define the collections of variables $\eta_{-j}=\\{\eta_{j}\\}_{j\geq
2}\setminus\\{\eta_{j}\\}$ and $\eta^{\prime}_{-j}$ similarly. We shall prove
the following claim:

$\operatorname{\mathbb{E}}\left[\left.\exp(\delta^{2}\epsilon^{2}m\rho_{j}^{2}\eta_{j}\eta^{\prime}_{j})\right|\eta_{-j}\eta^{\prime}_{-j}\right]\leq\exp(c\delta^{2}\epsilon^{4}(\delta^{2}m^{2}+mn)\rho_{j}^{4})$
(28)

for some universal constant $c>0$. Assuming that (28) holds, by induction we
can show that

$\displaystyle(\dagger)+1$
$\displaystyle\leq\exp(c\delta^{2}(\delta^{2}m^{2}+mn)\epsilon^{4}\sum_{j\geq
2}\rho_{j}^{4})$
$\displaystyle=\exp(c\delta^{2}(\delta^{2}m^{2}+mn)\epsilon^{4}/\|\lambda\|_{2,J}^{2}).$

Thus, if
$mn+\delta^{2}m^{2}=o\left(\|\lambda\|_{2,J}^{2}/(\delta^{2}\epsilon^{4})\right)$
then testing is impossible.

We now prove (28). Since the variable $\eta^{\prime}_{j}\eta^{\prime}_{j}$ is
either $1$ or $-1$, we have

$\displaystyle\operatorname{\mathbb{E}}\left[\left.\exp(\delta^{2}\epsilon^{2}m\rho_{j}^{2}\eta_{j}\eta^{\prime}_{j})\right|\eta_{-j}\eta^{\prime}_{-j}\right]$
$\displaystyle=(e^{\delta^{2}\epsilon^{2}m\rho_{j}^{2}}-e^{-\delta^{2}\epsilon^{2}m\rho_{j}^{2}})\cdot\mathbb{P}(\eta_{j}\eta^{\prime}_{j}=1|\eta_{-j}\eta^{\prime}_{-j})+e^{-\delta^{2}\epsilon^{2}m\rho_{j}^{2}}.$

Let us write $\eta_{\pm 1,j}$ for the vector of signs equal to $\eta$ but
whose $j$’th coordinate is $\pm 1$ respectively. Looking at the probability
above, and using the independence of $\eta,\eta^{\prime}$ given $Y$, we have

$\displaystyle\mathbb{P}(\eta_{j}\eta^{\prime}_{j}=1|Y,\eta_{-j},\eta^{\prime}_{-j})$
$\displaystyle=\mathbb{P}(\eta_{j}=1|Y,\eta_{-j})^{2}+\mathbb{P}(\eta_{j}=-1|Y,\eta_{-j})^{2}$
$\displaystyle=\frac{1}{4}\frac{(f_{\eta_{1j}}^{\otimes n}(Y))^{2}+(f^{\otimes
n}_{\eta_{-1j}}(Y))^{2}}{\left(\frac{1}{2}f^{\otimes
n}_{\eta_{1j}}(Y)+\frac{1}{2}f^{\otimes n}_{\eta_{-1j}}(Y)\right)^{2}}.$

Taking the expectation
$\operatorname{\mathbb{E}}[\,\cdot\,|\eta_{-j},\eta^{\prime}_{-j}]$ and using
the HM-AM inequality
$(\frac{1}{2}(x+y))^{-1}\leq\frac{1}{2}(\frac{1}{x}+\frac{1}{y})$ valid for
all $x,y>0$ gives

$\displaystyle\mathbb{P}(\eta_{j}\eta^{\prime}_{j}=1|\eta_{-j},\eta^{\prime}_{-j})$
$\displaystyle=\frac{1}{4}\int\frac{(\prod_{i=1}^{n}f_{\eta_{1j}}(x_{i}))^{2}+(\prod_{i=1}^{n}f_{\eta_{-1j}}(x_{i}))^{2}}{\frac{1}{2}\prod_{i=1}^{n}f_{\eta_{1j}}(x_{i})+\frac{1}{2}\prod_{i=1}^{n}f_{\eta_{-1j}}(x_{i})}\mu(\mathrm{d}x_{1})\dots\mu(\mathrm{d}x_{n})$
$\displaystyle\leq\frac{1}{4}+\frac{1}{8}\int\left(\frac{(\prod_{i=1}^{n}f_{\eta_{1j}}(x_{i}))^{2}}{\prod_{i=1}^{n}f_{\eta_{-1j}}(x_{i})}+\frac{(\prod_{i=1}^{n}f_{\eta_{-1j}}(x_{i}))^{2}}{\prod_{i=1}^{n}f_{\eta_{1j}}(x_{i})}\right)\mu(\mathrm{d}x_{1})\dots\mu(\mathrm{d}x_{n})=(\star).$

Note that $f_{\eta_{1j}}=f_{\eta_{-1j}}+2\epsilon\rho_{j}e_{j}$. Using the
lower bound $f_{\eta_{\pm 1j}}(x)\geq\frac{1}{2}$ for all $x\in\mathcal{X}$
and the inequality $1+x\leq\exp(x)$, we get

$\displaystyle(\star)$
$\displaystyle\leq\frac{1}{4}+\frac{1}{8}\left[\left(1+\int\frac{4\epsilon^{2}\rho_{j}^{2}e_{j}^{2}(x)}{f_{\eta_{-1j}}(x)}\mu(\mathrm{d}x)\right)^{n}+\left(1+\int\frac{4\epsilon^{2}\rho_{j}^{2}e_{j}^{2}(x)}{f_{\eta_{1j}}(x)}\mu(\mathrm{d}x)\right)^{n}\right]$
$\displaystyle\leq\frac{1}{4}(1+e^{8\epsilon^{2}n\rho_{j}^{2}}).$

Recall that $(\star)$ is a probability so $(\star)\leq 1$, and we obtain

$\displaystyle(\star)\leq\frac{1}{4}(1+e^{8\epsilon^{2}n\rho_{j}^{2}\wedge\ln
3}).$

Putting it together and applying Lemma E.5 we get

$\displaystyle\text{LHS of }\eqref{eqn:lower peeling claim}$
$\displaystyle\leq(e^{\delta^{2}\epsilon^{2}m\rho_{j}^{2}}-e^{-\delta^{2}\epsilon^{2}m\rho_{j}^{2}})\frac{1}{4}(1+e^{8\epsilon^{2}n\rho_{j}^{2}\wedge\ln
3})+e^{-\delta^{2}\epsilon^{2}m\rho_{j}^{2}}$ $\displaystyle\leq
e^{c\delta^{2}\epsilon^{4}\rho_{j}^{4}(\delta^{2}m^{2}+mn)}$

for universal $c=16>0$. Thus, by Lemma E.2 we obtain

$\displaystyle
1-\operatorname{TV}(\operatorname{\mathbb{E}}P_{0},\operatorname{\mathbb{E}}P_{1})\gtrsim\frac{1}{\chi^{2}(\operatorname{\mathbb{E}}P_{0},\operatorname{\mathbb{E}}P_{1})+1}\geq\exp(-c\delta^{2}\epsilon^{4}(\delta^{2}m^{2}+mn)/\|\lambda\|_{2,J}^{2})\stackrel{{\scriptstyle!}}{{=}}\Omega(\alpha).$

The necessity of

$mn+\delta^{2}m^{2}\gtrsim\frac{\log(1/\alpha)\|\lambda\|_{2,J}^{2}}{\delta^{2}\epsilon^{4}}$

follows immediately. In fact, this was even stronger than stated in Theorem
3.3 with
$mn+m^{2}\gtrsim\log(1/\alpha)\|\lambda\|_{2,J}^{2}/(\delta^{2}\epsilon^{4})$.222We
have $mn+m^{2}\leq(\sqrt{mn}+m)^{2}\leq 2(mn+m^{2})$, so
$\sqrt{mn}+m\asymp\sqrt{mn+m^{2}}$.

######  Lemma E.5.

For $a,b\geq 0$, the following inequality holds:

$\frac{1}{4}(e^{a}-e^{-a})(1+e^{b\wedge\ln 3})+e^{-a}\leq e^{2(ab+a^{2})}.$

###### Proof.

If $b\geq\ln 3$ or $a\geq 1$ we have:

$\displaystyle\text{LHS}\leq\frac{1}{4}(e^{a}-e^{-a})(1+e^{\ln
3})+e^{-a}=e^{a}\leq e^{\frac{b}{\ln 3}a+a^{2}}.$

If $b<\ln 3$ and $a<1$, we have

$e^{b}\leq 1+\frac{2}{\ln 3}b\leq 1+2b,\quad\frac{e^{a}+e^{-a}}{2}\leq
e^{a^{2}},\quad\frac{e^{a}-e^{-a}}{2}\leq\frac{e-e^{-1}}{2}a\leq 2a,$

and then

$\displaystyle\frac{1}{4}(e^{a}-e^{-a})(1+e^{b})+e^{-a}$
$\displaystyle=\frac{1}{2}(e^{a}+e^{-a})+\frac{e^{b}-1}{4}(e^{a}-e^{-a})$
$\displaystyle\leq e^{a^{2}}+2ab$ $\displaystyle\leq e^{a^{2}}(1+2ab)$
$\displaystyle\leq e^{a^{2}+2ab}$

The result follows from $\ln 3>1$. ∎

##  Appendix F Proofs From Section 4

###  F.1 Heuristic Justification of the Objective (7)

As usual, let $X,Y,Z$ denotes samples of sizes $n,n,m$ from
$P_{X},P_{Y},P_{Z}$ respectively. Let us give a heuristic justification for
using the training objective defined in (7) for the purpose of obtaining a
kernel for $\mathsf{LFHT}/\mathsf{mLFHT}$. Note that originally it was
proposed as a training objective for kernels to be used in two sample testing.
Recall that our test for $\mathsf{LFHT}$ can be written as

$\displaystyle\Psi_{1/2}(X,Y,Z)=\mathbbm{1}\Big{\\{}T_{\mathsf{LF}}\geq
0\Big{\\}}$

where

$\displaystyle T_{\mathsf{LF}}$
$\displaystyle=\operatorname{MMD}_{u}^{2}(\widehat{P}_{Z},\widehat{P}_{Y};K)-\operatorname{MMD}_{u}^{2}(\widehat{P}_{Z},\widehat{P}_{X};K),$

Heuristically, to maximize the power of ($\mathsf{mLFHT}$), we would like to
maximize the following population quantity

$J_{\mathsf{LF}}\triangleq\frac{\operatorname{\mathbb{E}}_{0}[T_{\mathsf{LF}}]-\operatorname{\mathbb{E}}_{1}[T_{\mathsf{LF}}]}{\sqrt{\operatorname{var}_{0}(T_{\mathsf{LF}})}}$

where

$\displaystyle\operatorname{\mathbb{E}}_{0}[T_{\mathsf{LF}}]$
$\displaystyle=\operatorname{\mathbb{E}}_{X,Y,Z}[T_{\mathsf{LF}}|{P_{Z}=P_{X}}]=+\operatorname{MMD}^{2}(P_{X},P_{Y};K),$
$\displaystyle\operatorname{\mathbb{E}}_{1}[T_{\mathsf{LF}}]$
$\displaystyle=\operatorname{\mathbb{E}}_{X,Y,Z}[T_{\mathsf{LF}}|{P_{Z}=P_{Y}}]=-\operatorname{MMD}^{2}(P_{X},P_{Y};K).$

Let $T_{\mathsf{TS}}=\operatorname{MMD}_{u}(\hat{P}_{X},\hat{P}_{Y})$ be the
usual statistic that is thresholded for two-sample testing. Then, a
computation analogous to that in Section D.2 show (cf. (21)) that

$\displaystyle\operatorname{var}_{0}(T_{\mathsf{LF}})$
$\displaystyle\approx\frac{A(K)}{n}+\frac{A(K)}{m}+\frac{B(K)}{n^{2}}+\frac{B(K)}{mn},$
$\displaystyle\operatorname{var}_{0}(T_{\mathsf{TS}})$
$\displaystyle\approx\frac{A(K)}{n}+\frac{B(K)}{n^{2}}$

for some $A(K)$ and $B(K)$. Therefore, we have approximately

$J_{\mathsf{LF}}\approx\frac{2\operatorname{MMD}^{2}(P_{X},P_{Y};K)}{\sqrt{1+\frac{n}{m}}\sqrt{\operatorname{var}_{0}(T_{\mathsf{TS}})}}\approx
2\sqrt{\frac{m}{m+n}}\widehat{J}(X,Y;K)$

which only differs from our optimization objective defined in (7) by a
constant factor.

Second, notice that
$\frac{\operatorname{MMD}(P_{X},P_{Y};K)}{\sqrt{\operatorname{var}(T_{\mathsf{TS}})}}$
depends only on $P_{X}-P_{Y}$ and that $((1-\delta)P_{X}+\delta
P_{Y})-P_{X}\propto P_{Y}-P_{X}$, therefore it is sensible to use (7) as our
training objective for is also sensible for ($\mathsf{mLFHT}$), and we don’t
even need to observe the sample $Z$.

###  F.2 Proof of Proposition 4.1

###### Proof.

In this proof we regard
$\mathcal{D}\triangleq(X^{\mathsf{tr}},X^{\mathsf{ev}},Y^{\mathsf{tr}},Y^{\mathsf{ev}})$
and the parameters of the kernel $\omega$ as fixed. Recall that we are looking
at the problem $\mathsf{mLFHT}$ with a misspecification parameter $R=0$ (see
Theorem 3.2). Given a test set $\\{z_{i}\\}_{i\in[m]}$, our test statistic is
$T(\\{z_{i}\\}_{i\in[m]})=\frac{1}{m}\sum_{i=1}^{m}f(z_{i})$ where

$f(z_{i})=\frac{1}{n_{\mathsf{ev}}}\sum_{j=1}^{n_{\mathsf{ev}}}\Big{(}K_{\omega}(z_{i},Y_{j}^{\mathsf{ev}})-K_{\omega}(z_{i},X_{j}^{\mathsf{ev}})\Big{)}.$

In Phase 3 of Algorithm 1, we observe the value
$\widehat{T}=T(Z)=\frac{1}{m}\sum_{i=1}^{m}f(Z_{i})$ and reject the null
hypothesis for large values of $\widehat{T}$. Thus, the $p$-value is defined
as

$\displaystyle p=p(Z,\mathcal{D})\triangleq\mathbb{P}_{\widetilde{Z}\sim
P_{X}^{\otimes m}}(T(\widetilde{Z})>\widehat{T}).$

Phase 2 of our Algorithm 1 produces random variables $T_{1},\dots,T_{k}$ that
all have the distribution of $T(\\{\widetilde{Z}_{i}\\}_{i\in[m]})$, so that
$\mathbbm{1}\\{T_{r}\geq\widehat{T}\\}$ ($r=1,\dots,k$) are unbiased estimates
of the $p$-value. However, the $T_{i}$ are not independent, because they
sample from the finite collection of calibration samples $X^{\mathsf{cal}}$.
However, as $n_{\mathsf{cal}}\to\infty$ the covariances between
$T_{r_{1}},T_{r_{2}}$ for $r_{1}\neq r_{2}$ tend to zero, and we obtain a
consistent estimate of $p$. ∎

###  F.3 Proof of Proposition 4.2

###### Proof.

The test statistic $T(X,Y,Z)$ in (3) is given by

$T(X,Y,Z)=\frac{1}{m}\sum_{i=1}^{m}f_{K}(Z_{i})$

where

$f_{K}(z)=\mu_{\widehat{P}_{Y}}(z)-\mu_{\widehat{P}_{X}}(z).$

This simplifies to (consider $K(x,y)=f(x)f(y)$)

$f_{K}(z)=\left(\frac{1}{n}\sum_{j=1}^{n}f(Y_{j})-\frac{1}{n}\sum_{j=1}^{n}f(X_{j})\right)f(z)=C(X,Y)f(z).$

where $C(X,Y)$ does not depend on $z$. Therefore, for any witness function
$f$, we obtain the desired additive test. ∎

###  F.4 Additive Test Statistics

In this section we prove accordingly that the test statistics of all of
MMD-M/G/O, SCHE, LBI, UME, RFM are of the form
$T_{f}(Z)=\frac{1}{m}\sum_{i=1}^{m}f(Z_{i})$ (where $f$ might depends on
$X,Y$). The test is to compare $T_{f}(Z)$ with some threshold $\gamma(X,Y)$.

Note that in the setting of Algorithm 1, the $X$ and $Y$ here correspond to
$X^{\mathsf{ev}}$ and $Y^{\mathsf{ev}}$.

MMD-M/G/O As described in (3) we have

$\displaystyle
T_{f}(Z)=\frac{1}{m}\sum_{i=1}^{m}\left(\frac{1}{n}\sum_{j=1}^{n}\left(K(Z_{i},Y_{j})-K(Z_{i},X_{j})\right)\right).$

SCHE As described in Section 4.2 we have

$\displaystyle
T_{f}(Z)=\frac{1}{m}\sum_{i=1}^{m}\mathbbm{1}\\{\phi(Z_{i})>t\\}.$

LBI As described in Section 4.2 we have

$\displaystyle
T_{f}(Z)=\frac{1}{m}\sum_{i=1}^{m}\log\left(\frac{\phi(Z_{i})}{1-\phi(Z_{i})}\right).$

UME As described in [29], the UME statistic evaluates the squared witness
function at ${J_{q}}$ test locations
$W=\\{w_{k}\\}_{k=1}^{J_{q}}\subset\mathcal{X}$. Formally for any two
distributions $P,Q$ we define

$\displaystyle
U^{2}(P,Q)=\lVert\mu_{Q}-\mu_{P}\rVert_{L^{2}(W)}^{2}=\frac{1}{{J_{q}}}\sum_{k=1}^{J_{q}}(\mu_{Q}(w_{k})-\mu_{P}(w_{k}))^{2}.$

However, we note a crucial difference that their result only considers the
case of $n=m$, and their proposed estimator for $U^{2}(P_{Z},P_{X})$ can not
be naturally extended to the case of $n\neq m$. Here we generalize it to
$m\neq n$ where we (conveniently) use a biased estimate of their distance.
Given samples $X,Y,Z$ and a set of witness locations $W$, the test statistic
is a (biased yet) consistent estimator of
$U^{2}(P_{Z},P_{Y})-U^{2}(P_{Z},P_{X})$. Let
$\psi_{W}(z)=\frac{1}{\sqrt{J_{q}}}(K(z,w_{1}),\dots,K(z,w_{J_{q}}))\in\mathbb{R}^{|W|}$
be the “feature function,” then:

$\displaystyle\widehat{U}^{2}(Z,X)$
$\displaystyle=\left\|\frac{1}{m}\sum_{i=1}^{m}\psi_{W}\left(Z_{i}\right)-\frac{1}{n}\sum_{j=1}^{n}\psi_{W}\left(X_{i}\right)\right\|^{2}_{2}$
$\displaystyle=\left\|\frac{1}{m}\sum_{i=1}^{m}\psi_{W}\left(Z_{i}\right)\right\|_{2}^{2}+\left\|\frac{1}{n}\sum_{j=1}^{n}\psi_{W}\left(X_{i}\right)\right\|_{2}^{2}-\frac{2}{mn}\sum_{1\leq
i\leq m,1\leq j\leq
n}\langle\psi_{W}\left(Z_{i}\right),\psi_{W}\left(X_{j}\right)\rangle$

Here $\braket{\cdot,\cdot}$ denotes the usual inner product. Therefore, the
difference between distances is

$\displaystyle\widehat{U}^{2}(Z,Y)-\widehat{U}^{2}(Z,X)=\frac{1}{m}\sum_{i=1}^{m}\left\langle\psi_{W}\left(Z_{i}\right),\frac{2}{n}\sum_{j=1}^{n}(\psi_{W}\left(X_{j}\right)-\psi_{W}\left(Y_{j})\right)\right\rangle+F(X,Y)$

where $F$ is sum function based only on $X,Y$. This is clearly an additive
statistic for $Z$.

RFM Algorithm 1 in [41] describes a method for learning a kernel from data
given a binary classification task. For convenience lets concatenate the data
to $X^{\text{RFM}}=(X,Y)\in\mathbb{R}^{2n\times d}$ and labels
$y^{\text{RFM}}=(\vec{0}_{n},\vec{1}_{n})\in\mathbb{R}^{1\times 2n}$. Given a
learned kernel $K$, we write the Gram matrix as
$(K(X^{\text{RFM}},X^{\text{RFM}}))_{i,j}=K(X^{\text{RFM}}_{i},X^{\text{RFM}}_{j})$
($1\leq i,j\leq 2n$). Let $K(X^{\text{RFM}},z)$ be a column vector with
components $K(X^{\text{RFM}}_{i},z)$ ($1\leq i\leq 2n$). The classifier is
then defined as

$\displaystyle f^{\text{RFM}}(z)=y^{\text{RFM}}\cdot
K(X^{\text{RFM}},X^{\text{RFM}})^{-1}\cdot K(X^{\text{RFM}},z).$ (29)

Though in [41] the kernel learned from RFM is used to construct a classifier
as in Equation 29, since RFM is a feature learning method, we also apply the
RFM kernel to our MMD test, namely

$\displaystyle f^{\text{RFM to
MMD}}(z)=\frac{1}{n}\sum_{j=1}^{n}\left(K(z,Y_{j})-K(z,X_{j})\right).$

##  Appendix G Application: Diffusion Models vs CIFAR

We defer a more fine-grained detail to our code submission, which includes
executable programs (with PyTorch) once the data-generating script from DDPM
has been run (see README in the ./codes/CIFAR folder).

###  G.1 Dataset Details

We use the CIFAR-10 dataset available online at
https://www.cs.toronto.edu/~kriz/cifar.html, which contains 50000 colored
images of size $32\times 32$ with 10 classes. For the diffusion generated
images, we use the SOTA Hugging Face model (DDPM) that can be found at
https://huggingface.co/google/ddpm-CIFAR-10-32. We generated 10000 artificial
images for our experiments. The code can be found at our code supplements.

For dataset balancing, we randomly shuffled the CIFAR-10 dataset and used
10000 images as data in our code. Most of our experiments are conducted with
the null $P_{X}$ as CIFAR images, and the alternate as
$P_{Y}=\frac{2}{3}\cdot\text{CIFAR}+\frac{1}{3}\cdot\text{DDPM}$. To this end,
we matched 20000 images from CIFAR to belong to the alternate hypothesis, and
the remaining 30000 images to stay in the null hypothesis. For the alternate
dataset, we simply sample without replacement from the $20000+10000$ mixture.
This sampled distribution is _almost_ the same as mixing (so long as the
sample bank is large enough compared to the acquired data, so that each item
in the alternate has close to $1/3$ probability of being in DDPM, which is
indeed the case).

Figure 4: Data visualization for CIFAR-10 (left) vs DDPM diffusion generated
images (right)

###  G.2 Experiment Setup and Benchmarks

We use a standard deep Conv-net [37], which has been employed for SOTA GAN
discriminator tasks in similar settings. It has four convolutional layers and
one fully connected layer outputting the feature space of size $(300,1)$. For
SCHE and LBI, we simply added a linear layer of $(300,2)$ after applying ReLU
to the 300-dimensional layer and used the cross-entropy loss to train the
network. Note that this is equivalent to first fixing the feature space and
then performing logistic regression to the feature space. For kernels, we add
extra trainable parameters after the $300$-d feature output.

For the $\operatorname{MMD}$-based tests, we simply train the kernel on the
neural net and evaluate our objective. For UME, we used a slightly generalized
version of the original statistic in [29] which allows for comparison on
randomly selected witness locations in the null hypothesis with $m\neq n$ (see
Section F.4). The kernel is trained using our heuristic (see (7) and Section
F.1), with MMD replaced by UME. The formula for UME variance can be found in
[29]. For RFM, we use Algorithm 1 in [41] to learn a kernel on (stochastic
batched) samples, and then use our MMD test on the trained kernel.

We use 80 training epochs for most of our code from the CNN architecture (for
classifiers, this is well after interpolating the training data and roughly
when validation loss stops decreasing), and a batch size of 32 which has a
slight empirical benefit compared to larger batch sizes. The learning rates
are tuned separately in MMD methods for optimality, whereas for classifiers
they follow the discriminator’s original setting from [37]. In Phase 2 of
Algorithm 1, we choose $k=1000$ for the desired precision while not
compromising runtime. For each task, we run $10$ independent models and report
their performances as the mean and standard deviation of those $10$ runs as
estimates. We refer to a full set of hyper-parameters in our code
implementation.

Our code is implemented in Python 3.7 (PyTorch 1.1) and was ran on an NVIDIA
RTX 3080 GPU equipped with a standard torch library and dataset extensions.
Our code setup for feature extraction is similar to that of [36]. For
benchmark implementations, our code follows from the original code templated
provided by the cited papers.

Figure 5: Relevant plots following the setting in Figure 2 (in the main text)
of fixing $n_{\mathsf{tr}}=1920$ and varying sample size $m$ in the x-axis for
the comparison with missing benchmarks. Errorbars are projected showing
standard deviation across 10 runs. We replaced part (d) in Figure 2 (in the
main text) to a sanity check in our FPR when thresholded at $\alpha=0.05$.

###  G.3 Sample Allocation

We make a comment on why (4) is _different_ from just thresholding
$\widehat{\operatorname{MMD}^{2}}(Z,Y^{\mathsf{tr}})-\widehat{\operatorname{MMD}^{2}}(Z,X^{\mathsf{tr}})$
at 0, which was what we did in part (c) of Figure 2 (and hence the difference
along the curve of MMD-M vs Figure 1). Our theory assumes that the samples are
i.i.d. conditioned on the kernel being chosen already. However, in the
experiments, the kernel is dependent on the training data. Therefore, to
evaluate the MMD estimate (between experimentations), one needs extra data
that does not intersect with training.

In fact, it can be experimentally shown by comparing Figure 1 and Figure 2(c)
that doing so (while reducing the sample complexity on $n_{\mathsf{ev}}$)
hurts performance. Indeed, we found out that when
$X^{\mathsf{ev}},Y^{\mathsf{ev}}$ are non-intersecting with training,
performance is (almost) always better at a cost of hurting the overall sample
complexity of $n$.

###  G.4 Remarks on Results

Figure 5 lists all of our benchmarks in the setting of Figure 2 (in the main
text) on missing benchmarks, where the last figure is replaced by the false
positive rate at thresholding at $\alpha=0.05$ to verify our results. As
mentioned in the main text, our MMD-M method consistently outperforms other
benchmarks on both the expected $p$-value (of alternate) and rejection rate at
$\alpha=0.05$, while all of our tests observe a empirical false positive rate
close to $\alpha=0.05\%$ (Part (b)), showing the consistency of methods.

##  Appendix H Application: Higgs-Boson Detection

###  H.1 Dataset Details

We use the Higgs dataset available online at
http://archive.ics.uci.edu/ml/datasets/HIGGS, produced using Monte Carlo
simulations [5]. The dataset is nearly balanced, containing $5,829,122$ signal
instances and $5,170,877$ background instances. Each instance is a
$28$-dimensional vector, consisting of $28$ features. The first $21$ features
are kinematic properties measured by the detectors in the accelerator, such as
momentum and energy. The last $7$ properties are _invariant masses_ , derived
from the first $21$ features.

Figure 6: This figure visualizes the distribution of the $26$th feature, the
invariant mass $m_{Wbb}$. The red and black lines are the histograms of the
original dataset. We employ MMD-M as a classifier, trained and evaluated using
$n_{\mathsf{tr}}=1.3\times 10^{6}$ and
$n_{\mathsf{ev}}=n_{\mathsf{opt}}=2\times 10^{4}$ through Section H.3.2. The
blue(green) line represents all instances $z$’s whose “witness scores”
$f(z;X^{\mathsf{ev}},Y^{\mathsf{ev}})$’s are larger(smaller) than
$t_{\mathsf{opt}}$.

###  H.2 Experiment Setup and Training Models

The modified Algorithm 1 is shown in Section H.3.2 and Section H.3.2. Compared
with Section H.3.2, we implement the thresholding trick (Section 4.3) in
Section H.3.2.

####  H.2.1 Configuration and Model Architecture

We implement all methods in Python 3.9 and PyTorch 1.13 and run them on an
NVIDIA Quadro RTX 8000 GPU.

For all classifier-based methods in this study (SCHE and LBI), we adopt the
same architecture as previously proposed in [5]. The classifiers are six-layer
neural networks with 300 hidden units in each layer, all employing the tanh
activation function. For SCHE, the output layer is a single sigmoid unit and
we utilize the binary cross-entropy loss for training. For LBI, the output
layer is a linear unit and we utilize the binary cross entropy loss combined
with a logit function (which is more numerically stable than simply using a
sigmoid layer followed by a cross entropy loss).

For all MMD-based methods (MMD-M, MMD-G, MMD-O, and UME), the networks
$\varphi$ and $\varphi^{\prime}$ are both six-layer neural networks with 300
ReLU units in each layer. The feature space, which is the output of the neural
network $\varphi$, is set to be $100$-dimensional. Here UME has the same
kernel architecture as MMD-M, and the number of test locations is set to be
$J_{q}=4096$. For RFM, we adopt the same architecture as in [41], where the
kernel is $K_{M}(x,y)=\exp(-\gamma(x-y)^{T}M(x-y))$ with a constant $\gamma$
and a learnable positive semi-definite matrix $M$. We set $\gamma\equiv 1$.

The neural networks are initialized using the default setting in PyTorch, and
the bandwidths $\sigma,\sigma^{\prime}$ are initialized using the _median
heuristic_ [21]. The parameter $\tau$ is initially set to $0.5$. For UME, the
witness locations $W$ are initially randomly sampled from the training set.
For RFM, the initial $M$ equals the median bandwidth times an identity matrix.

####  H.2.2 Training

The size of our training set, denoted as $n_{\mathsf{tr}}$, varies from
$1.0\times 10^{2}$ to $1.6\times 10^{6}$. For a given $n_{\mathsf{tr}}$, we
select the first $n_{\mathsf{tr}}$ datapoints from each class of the Higgs
dataset to form $X^{\mathsf{tr}}$ and $Y^{\mathsf{tr}}$, i.e.,
$|X^{\mathsf{tr}}|=|Y^{\mathsf{tr}}|=n_{\mathsf{tr}}$. Subsequently, we
randomly select
$n_{\mathsf{validation}}=\min(\sqrt{10n_{\mathsf{tr}}},0.1n_{\mathsf{tr}})$
points from each of $X_{\mathsf{tr}},Y_{\mathsf{tr}}$ to constitute the
validation set, while the remainder of $X_{\mathsf{tr}},Y_{\mathsf{tr}}$ are
used for running gradient descent. The optimizer is set to be a minibatch SGD,
with a batch size of $1024$, a learning rate of $0.001$, and a momentum of
$0.99$. Training is halted once the validation loss stops to decrease for $10$
epochs, then we choose the checkpoint (saved for each epoch) with the smallest
validation loss thus far as our trained model. Beyond the general setting
above, in RFM a batch size of $1024$ doesn’t work well and instead we use a
batch size of $20,000$.

###  H.3 Evaluating the Performance

####  H.3.1 Evaluating the p-Value with the Methodology of Algorithm 1

We call the “witness score” of an instance $z\in\mathcal{X}$ as

$f(z;X^{\mathsf{ev}},Y^{\mathsf{ev}})=\frac{1}{n_{\mathsf{cal}}}\sum_{i=1}^{n_{\mathsf{cal}}}\left(k(z,Y_{i}^{\mathsf{ev}})-k(z,X_{i}^{\mathsf{ev}})\right).$
(30)

For a vector of instances $Z=(Z_{1},\ldots,Z_{m})$, we write

$f(Z;X^{\mathsf{ev}},Y^{\mathsf{ev}})=(f(Z_{1};X^{\mathsf{ev}},Y^{\mathsf{ev}}),\ldots,f(Z_{m};X^{\mathsf{ev}},Y^{\mathsf{ev}})).$

The testing procedure is summarized in Phases 2, 3 and 4 in Section H.3.2 and
Section H.3.2. In the Higgs experiment, we utilize the Gaussian approximation
method to determine the p-values when the witness function $f$ is not
thresholded, which allows us to reach very small p-values and errors under
limited computational resource. In cases where the score function $f$ is
thresholded by a value $t$, using the Binomial distribution as in Section
H.3.2 is more precise and also fast enough.

Given a trained kernel $K$ trained on $X^{\mathsf{tr}}$ and $Y^{\mathsf{tr}}$,
we set $X^{\mathsf{ev}}=X^{\mathsf{tr}}$ and
$Y^{\mathsf{ev}}=Y^{\mathsf{tr}}$, and accordingly
$n_{\mathsf{ev}}=n_{\mathsf{tr}}$. This results in a more efficient use of
data (since we reuse $X^{\mathsf{tr}},Y^{\mathsf{tr}}$ also as
$X^{\mathsf{ev}},Y^{\mathsf{ev}}$). Then, out of the untouched portion of the
data, we randomly choose $n_{\mathsf{cal}}=20,000$ datapoints from both
classes to populate $X^{\mathsf{cal}}$ and $Y^{\mathsf{cal}}$, i.e.,
$|X^{\mathsf{cal}}|=|Y^{\mathsf{cal}}|=n_{\mathsf{cal}}=20,000$. In addition
to the general setting above, for RFM, we need to solve a
$2n_{\mathsf{ev}}$-dimensional linear equation during inference, which arises
from the inverse matrix in Equation 29 (solving
$K(X^{\text{RFM}},X^{\text{RFM}})\boldsymbol{u}=(y^{\text{RFM}})^{T}$ for
$\boldsymbol{u}\in\mathbb{R}^{2n_{\mathsf{ev}}}$). So we set
$n_{\mathsf{ev}}=\min(n_{\mathsf{tr}},10,000)$ that
$X_{\mathsf{ev}},Y_{\mathsf{ev}}$ are randomly sampled from the training set.

In order to compare different benchmarks, we evaluate the expected
significance of discovery on a mixture of 1000 backgrounds and 100 signals.
For each benchmark and each $n_{\mathsf{tr}}$, we train $10$ independent
models. Then for each trained model we proceed through the Phases 2, 3 (and 4)