# A Kernel Measure of Dissimilarity between $M$ Distributions

Zhen Huang
Department of Statistics, Columbia University
e-mail<EMAIL_ADDRESS>
and
Bodhisattva Sen
Department of Statistics, Columbia University
e-mail<EMAIL_ADDRESS>Supported by NSF grant DMS-2015376

###### Abstract

Given $M\geq 2$ distributions defined on a general measurable space, we
introduce a nonparametric (kernel) measure of multi-sample dissimilarity (KMD)
— a parameter that quantifies the difference between the $M$ distributions.
The population KMD, which takes values between 0 and 1, is 0 if and only if
all the $M$ distributions are the same, and 1 if and only if all the
distributions are mutually singular. Moreover, KMD possesses many properties
commonly associated with $f$-divergences such as the data processing
inequality and invariance under bijective transformations. The sample estimate
of KMD, based on independent observations from the $M$ distributions, can be
computed in near linear time (up to logarithmic factors) using $k$-nearest
neighbor graphs (for $k\geq 1$ fixed). We develop an easily implementable test
for the equality of $M$ distributions based on the sample KMD that is
consistent against all alternatives where at least two distributions are not
equal. We prove central limit theorems for the sample KMD, and provide a
complete characterization of the asymptotic power of the test, as well as its
detection threshold. The usefulness of our measure is demonstrated via real
and synthetic data examples; our method is also implemented in an R package.

Keywords: Asymptotic power behavior, detection threshold, $k$-nearest neighbor
graph, multi-distribution $f$-divergence, nonparametric test for equality of
distributions

##  1 Introduction

Suppose that $\mathcal{Z}$ is a general measurable space, and we have $M$
distributions $P_{1},\ldots,P_{M}$ on $\mathcal{Z}$ from which we observe
independent samples. A natural statistical question to ask here is: “Are these
$M$ distribution the same? If they are not the same, then how different are
they?”. In this paper, we answer these questions by proposing a nonparametric
measure that quantifies the differences between the multiple samples. To
rephrase this in statistical language, we define a measure
$\eta\equiv\eta(P_{1},\ldots,P_{M})$ such that:

  * (i)

$\eta$ is a deterministic number between $[0,1]$;

  * (ii)

$\eta=0$ if and only if $P_{1}=\ldots=P_{M}$ (i.e., all the $M$ distributions
are the same), and

  * (iii)

$\eta=1$ if and only if the $M$ distributions are mutually singular, i.e.,
there exist disjoint measurable sets $\\{A_{i}\\}_{i=1}^{M}$ such that
$P_{i}(A_{i})=1$, for $i=1,\ldots,M$. Thus, $\eta=1$ quantifies that the $M$
distributions are very different.

Moreover, any value between 0 and 1 of $\eta$ would convey an idea of how
different these $M$ distributions are. We call our proposed $\eta$ as a kernel
measure of multi-sample dissimilarity (KMD) as its definition involves a
positive semi-definite kernel matrix.

While there is a rich literature on the multi-sample testing problem:

${\rm H}_{0}:P_{1}=\ldots=P_{M}\quad\qquad{\rm against}\quad\qquad{\rm
H}_{1}:P_{i}\neq P_{j},\ {\rm for\ some\ }1\leq i<j\leq M,$ (1)

(see Section 1.1 for a detailed review) most of these tests do not quantify to
what extent these distributions are different, when the null is violated.
Moreover, many popular distances or similarities — such as the KL-divergence,
Hellinger distance, and other $f$-divergences [59, 38, 7] — that quantify the
difference between two distributions can be difficult to estimate when
$\mathcal{Z}$ is not a Euclidean space. In modern statistical applications, it
is quite often necessary to compare more than two distributions, and they can
be defined on a general measurable space. Two motivating applications in this
direction are:

Example 1 (Multivariate functional data): In speech recognition, data may be
viewed as functional inputs.

Figure 1: Four instances of spoken Arabic digits, each of which contains time
series of 13 Mel Frequency Cepstrum Coefficients (MFCCs).

For example, Figure 1 shows four instances of spoken Arabic digits from the
ArabicDigits data set in Gorecki et al. [32], each of which contains time
series of 13 Mel Frequency Cepstrum Coefficients (MFCCs). Hence, in this
scenario, $\mathcal{Z}$ is the space of all 13-dimensional time series. A
natural problem here is to distinguish the 10 spoken digits, which can be seen
as $M=10$ different distributions on $\mathcal{Z}$. Knowing how different
these $M$ distributions are on a scale between 0 to 1 provides a sense of how
well any machine learning algorithm can perform in distinguishing these $M$
digits. One can also study how the digits spoken may differ by gender. For
example, how males speak the number 1 should be different from how males speak
the number 2, but this difference should be smaller than that between how
females speak 1 and how males speak 2. A statistical measure that quantifies
the extent of these differences can be obtained by our procedure.

Example 2 (Distribution over documents): In sentiment analysis [48, 65], the
goal is to identify the sentiment (e.g., positive, negative, neutral) of a
document. These documents may each contain hundreds of words, and are
typically of different lengths. For example, Maas et al. [48] considered 2000
movie reviews, each of which is either positive or negative. If there are in
total $M$ sentiments, then there are naturally $M$ distributions in the space
of documents — each distribution corresponding to a different sentiment. Let
$V$ be the vocabulary, i.e., the set of all possible words. Then, in this
case, $\mathcal{Z}=\\{(v_{1},\ldots,v_{n}):v_{i}\in V,n\in\mathbb{N}\\}$ with
each element in $\mathcal{Z}$ denoting a generic document. Given a specific
data set with $n_{i}$ documents, conveying the $i$-th sentiment, for
$i=1,\ldots,M$, a data analyst may want to know how different the $M$
sentiments are in this data set. Our procedure can quantify this difference,
so as to suggest whether it is possible to effectively classify the different
sentiments.

Often in practice, we have a distance between two objects in $\mathcal{Z}$
which has already been demonstrated to be useful in the domain area. For the
first example above, a useful metric is the dynamic time warping (DTW)
distance between multi-dimensional time series [9, 31], a distance which could
account for differences in speaking rates between speakers and can be computed
even when the two time series are discretized into different lengths. For the
second example above, one can use the Jaccard distance or other distances
between documents [65]. In these situations, applying usual Euclidean methods
may require additional embedding efforts (see e.g., [35] on document
embedding). In contrast, it will be seen that our method can directly provide
an easy, practical, and interpretable way to quantify the difference between
multiple distributions, as long as a metric is available on the space
$\mathcal{Z}$.

In this paper, given $n_{i}$ independent observations from $P_{i}$,
$i=1,\ldots,M$, we propose and study an empirical estimator $\hat{\eta}$ of
$\eta$, constructed using geometric graphs [21, 10] (e.g., the $k$-nearest
neighbor ($k$-NN) graph for $k\geq 1$). The main contributions of the paper,
and some important properties of $\eta$ and $\hat{\eta}$ are summarized below:

  1. 1.

We propose a nonparametric measure of multi-sample dissimilarity
$\eta=\eta(P_{1},\ldots,P_{M})\in[0,1]$ that satisfies properties (i)–(iii)
mentioned above (see Theorem 1). Moreover, any value of $\eta$, between 0 and
1, conveys an idea about how different these distributions are. For example,
in a large class of location and scale distributions, $\eta$ increases as the
“difference” between the parameters gets larger (see Proposition 1). Our
$\eta$ also satisfies many desirable properties commonly associated with
$f$-divergences [20], including the data processing inequality and invariance
under bijective transformations (Proposition 2) and joint convexity
(Proposition 3). Indeed, $\eta$ is a member of the multi-distribution
$f$-divergence proposed in García-García and Williamson [30]; however, no
estimation strategy was given in [30].

  2. 2.

We develop an estimator $\hat{\eta}$ of $\eta$, which is consistent (Theorem
2), interpretable, easily implementable and computationally efficient. It can
be computed in linear time (up to logarithmic factors); this is an enormous
reduction from the $O(n^{2})$ complexity of energy and kernel based methods
[70, 33], and other $O(n^{3})$ complexity / NP complete approaches [62, 13,
22, 49, 37]. We further show that $\hat{\eta}$ is a generalization of the two-
sample statistic based on $k$-NN proposed in Schilling [63] and Henze [40]
(see Remark 2). The multi-sample test (1) based on $\hat{\eta}$ is also
consistent against all alternatives for which $P_{i}\neq P_{j}$, for some
$1\leq i<j\leq M$ (Corollary 2).

  3. 3.

The asymptotic distribution of $\hat{\eta}$ is Gaussian, which yields easy-to-
use asymptotic tests. Under the null hypothesis (1), the permutation
distribution111Given the pooled sample without the information on their sample
identities, any permutation of the sample identities is equally likely under
the null hypothesis (1). This conditional distribution is referred to as the
permutation distribution in this paper. and the unconditional distribution of
$\hat{\eta}$ are both asymptotically normal (Theorem 3). Further, when
$\mathcal{Z}$ is a Euclidean space, the asymptotic null distribution is
distribution-free if the common distribution has a Lebesgue density (see
Theorem 4 for a precise statement).

  4. 4.

We further provide in Section 5 a complete characterization of the asymptotic
power of the test for (1) based on $\hat{\eta}$ (see Theorem 5), and its
detection threshold, using the technique of Poissonization (cf., the recent
paper Bhattacharya [11] where the detection threshold for as a class of graph-
based two-sample tests is studied). In particular, we show that, under both
fixed and shrinking alternatives converging to the null, $\hat{\eta}$ has an
asymptotic normal distribution after proper centering; see Appendix B
(Theorems B.1 and B.2).

  5. 5.

In Section 6, we demonstrate the usefulness of the proposed methodology via
real and synthetic data examples. Our method is also implemented in an R
package222See https://cran.r-project.org/package=KMD and
https://github.com/zh2395/KMD..

The outline of the paper is as follows. In Section 2, we formally define
$\eta$ and investigate its properties. Its estimator $\hat{\eta}$ is studied
in Section 3 along with its basic characteristics. In Section 4 we provide the
asymptotic distribution of $\hat{\eta}$, under the null (1). A rigorous study
of the power behavior of the test based on $\hat{\eta}$ and its detection
threshold is given in Section 5. Simulations and real data experiments are
provided in Section 6. All the proofs of our main results, further
discussions, additional numerical experiments, and more results on the
behavior of $\hat{\eta}$ under the alternative are given in Appendices A-D.

###  1.1 Related Works

The $M$-sample testing problem (1) has been extensively studied in the
statistics literature, both in parametric and nonparametric regimes.
Parametric tests (e.g., $t$-test, MANOVA, likelihood ratio tests, Wald tests)
are provably powerful when the underlying model assumptions hold true (see
e.g., [44] and the references therein), but could have poor performance when
the model is misspecified. In comparison, nonparametric methods are usually
powerful against general alternatives, under much more relaxed assumptions on
the data distributions; this will be the framework adopted in this paper.

There is a long history of nonparametric two-sample tests. In the one-
dimensional case, classical well-known distribution-free tests include the
Kolmogorov-Smirnov test [66] and the Wald-Wolfowitz run test [71]. In the
multivariate setting, Friedman and Rafsky [28] proposed generalizations of the
Wald-Wolfowitz test based on the minimum spanning tree of the pooled sample.
Its theoretical properties were further analyzed by Henze and Penrose [41].
Multivariate two-sample tests based on nearest neighbor ideas were proposed in
[63, 40, 36]. Chen and Friedman [17] extended a number of two-sample tests
based on type coincidences in a geometric graph. Bhattacharya [10] proposed a
general asymptotic framework for studying graph-based two-sample tests; also
see [11]. Apart from the above tests based on geometric graphs, there are two-
sample tests based on data depth [46, 73]. In the past decade, energy
statistics [70] and tests based on the maximum mean discrepancy (MMD) [33] in
the kernel literature have also drawn great attention in two-sample testing,
and their equivalence has been established [64]. There are also methods that
achieve finite-sample distribution-freeness by using minimum non-bipartite
matching [62, 49], the shortest Hamiltonian path [13], or multivariate ranks
defined via optimal transport [22, 37].

When moving to general $M$-sample testing, Petrie [57] generalized a number of
graph-theoretic tests to the multi-sample scenario. Recently Mukherjee et al.
[49] generalized the test in [62] based on minimum non-bipartite matching,
retaining the exact distribution-freeness in finite samples. Energy statistic
has also been used in Rizzo and Székely [61] for $M$-sample testing, providing
a nonparametric extension of ANOVA. Liu and Singh [46] proposed an index $Q$
between $[0,1]$ measuring the dissimilarity between two distributions, based
on data depth, with the null being achieved when $Q=\frac{1}{2}$; however $Q$
does not satisfy properties (i)-(iii). Although some $f$-divergences, such as
the Hellinger distance [38] or the total variation distance, satisfy
(i)-(iii), they cannot be easily extended beyond $M=2$.

##  2 KMD: The Population Version

In this section we define the population version of our measure of
dissimilarity $\eta$ between the $M$ distributions $P_{1},\ldots,P_{M}$. Our
definition of $\eta$ involves the use of a reproducing kernel Hilbert space
(RKHS) over the finite discrete space ${\mathcal{S}}:=\\{1,\ldots,M\\}$, and
hence we call our proposal the kernel measure of multi-sample dissimilarity
(KMD). Although the discrete kernel $K(x,y):=I(x=y)$, for
$x,y\in{\mathcal{S}}$, seems to be the most natural choice of the kernel over
the discrete space ${\mathcal{S}}$, our results are applicable to other
kernels as well.

Reproducing Kernel Hilbert Space (RKHS): While there is a general theory of
RKHS on arbitrary spaces, the RKHS over a finite space is much simpler, which
will be introduced in the following. For an introduction to the theory of RKHS
and its applications in statistics we refer the reader to [8, 69]. By a kernel
function $K:{\mathcal{S}}\times{\mathcal{S}}\to\mathbb{R}$ we mean a symmetric
and nonnegative definite function, i.e., the matrix $[K(i,j)]_{i,j=1}^{M}$ is
positive semi-definite. The kernel $K(\cdot,\cdot)$ is said to be
characteristic if for any $(\alpha_{1},\ldots,\alpha_{M})\neq(0,\ldots,0)$
with $\sum_{i=1}^{M}\alpha_{i}=0$, we have
$\sum_{i,j=1}^{M}\alpha_{i}\alpha_{j}K(i,j)>0$. Note that the usual definition
of a characteristic kernel is through the uniqueness of the kernel mean
embedding [67], but it is equivalent to this simpler definition when the space
${\mathcal{S}}$ is finite; see e.g., [67, Section 4.4].

###  2.1 Definition of $\eta$

Suppose that we have $n_{i}$ independent observations
$X_{i1},\ldots,X_{in_{i}}$ from the distribution $P_{i}$ taking values in
$\mathcal{Z}$, for $i=1,\ldots,M$. Denote the pooled sample
$\\{X_{11},\ldots,X_{1n_{1}},\ldots,X_{M1},\ldots,$ $X_{Mn_{M}}\\}$ as
$\\{Z_{1},\ldots,Z_{n}\\}$ with $n=n_{1}+\ldots+n_{M}$, and the corresponding
labels as $\Delta_{1},\ldots,\Delta_{n}$, i.e.,
$\Delta_{j}=i\in\\{1,\ldots,M\\}$ if $Z_{j}$ comes from distribution $P_{i}$.
If $\frac{n_{i}}{n}\approx\pi_{i}\in(0,1)$ such that
$\sum_{i=1}^{M}\pi_{i}=1$, then $\\{(\Delta_{i},Z_{i})\\}_{i=1}^{n}$ can be
“approximately” thought of as an i.i.d. sample of size $n$ from
$(\tilde{\Delta},\tilde{Z})$ whose distribution is specified as follows:

  1. 1.

$\mathbb{P}(\tilde{\Delta}=i)=\pi_{i}$, for $i=1,\ldots,M$.

  2. 2.

Given $\tilde{\Delta}=i$, $\tilde{Z}$ is drawn from distribution $P_{i}$.

The following lemma (proved in Appendix C.1) is a crucial observation that
will allow us to formally define $\eta$.

######  Lemma 1.

The $M$ distributions $P_{1},\ldots,P_{M}$ are the same if and only if
$\tilde{\Delta}\perp\\!\\!\\!\perp\tilde{Z}$, and the $M$ distributions are
mutually singular if and only if $\tilde{\Delta}$ is a function of
$\tilde{Z}$.

The connection between $M$ sample testing and measures of association between
$\tilde{\Delta}$ and $\tilde{Z}$ has been noted before; see e.g., [29, 50].
The above lemma motivates the use of a certain measure of association [16, 4,
21, 42] between $\tilde{\Delta}$ and $\tilde{Z}$ to quantify the dissimilarity
between the $M$ distributions. In particular, we adopt the ideas from the
kernel measure of association (KMAc) proposed in Deb et al. [21]. Let $\mu$ be
the distribution of $(\tilde{Z},\tilde{\Delta})$ as defined above. Here we
assume that $\frac{n_{i}}{n}$ converges to some $\pi_{i}\in(0,1)$ as
$n\to\infty$, such that $\sum_{i=1}^{M}\pi_{i}=1$. Let
$(\tilde{Z}_{1},\tilde{\Delta}_{1})$ and $(\tilde{Z}_{2},\tilde{\Delta}_{2})$
be i.i.d. $\mu$. Also, let
$(\tilde{Z},\tilde{\Delta},\tilde{\Delta}^{\prime})$ be such that
$(\tilde{Z},\tilde{\Delta})\sim\mu$,
$(\tilde{Z},\tilde{\Delta}^{\prime})\sim\mu$, and $\tilde{\Delta}$ and
$\tilde{\Delta}^{\prime}$ are conditionally independent given $\tilde{Z}$.
Define the kernel measure of multi-sample dissimilarity (KMD) as:

$\displaystyle\eta\equiv\eta(P_{1},\ldots,P_{M}):=\frac{\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})\big{]}-\mathbb{E}\big{[}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\big{]}}{\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta})\big{]}-\mathbb{E}\big{[}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\big{]}}.$
(2)

Note that $\eta$, as defined above, also depends on the mixing proportions
$\\{\pi_{i}\\}_{i=1}^{M}$, but for notational simplicity we do not highlight
this dependence here. If we let $\bar{P}:=\sum_{i=1}^{M}\pi_{i}P_{i}$ be the
distribution of $\tilde{Z}$. Given $\tilde{Z}=z$, the conditional probability
of $\tilde{Z}$ coming from $P_{i}$ is $\frac{{\rm d}(\pi_{i}P_{i})}{{\rm
d}\bar{P}}(z)$, where $\frac{{\rm d}(\pi_{i}P_{i})}{{\rm d}\bar{P}}$ denotes
the Radon-Nikodym derivative of $\pi_{i}P_{i}$ with respect to $\bar{P}$.
Hence $\eta$ in (2) has the alternative expression:

$\eta=\frac{\int\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)\frac{{\rm d}P_{i}}{{\rm
d}\bar{P}}(z)\frac{{\rm d}P_{j}}{{\rm d}\bar{P}}(z){\rm
d}\bar{P}(z)-\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)}{\sum_{i=1}^{M}\pi_{i}K(i,i)-\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)}.$
(3)

When $P_{1},\ldots,P_{M-1}$ are absolutely continuous w.r.t. $P_{M}$, $\eta$
reduces to a member of the multi-dimensional $f$-divergences proposed in
García-García and Williamson [30].

Observe that if all the $M$ distributions $P_{1},\ldots,P_{M}$ are the same,
then $\tilde{\Delta}$ is independent of $\tilde{Z}$ (by Lemma 1), and thus
$(\tilde{\Delta},\tilde{\Delta}^{\prime})$ has the same distribution as
$(\tilde{\Delta}_{1},\tilde{\Delta}_{2})$; so the numerator of $\eta$ in (2)
equals 0. If all the $M$ distributions are mutually singular, then
$\tilde{\Delta}$ is a function of $\tilde{Z}$ (by Lemma 1), and $\eta=1$ (as
$\tilde{\Delta}^{\prime}=\tilde{\Delta}$). The converse is also true, i.e.,
$\eta=0$ (resp. $\eta=1$) also implies $P_{1}=\ldots=P_{M}$ (resp. all the $M$
distributions are mutually singular); thus $\eta$ satisfies properties
(i)–(iii) mentioned at the beginning of the Introduction. We formalize this in
the following result (proved in Appendix C.2).

######  Theorem 1.

Suppose the kernel $K(\cdot,\cdot)$ is characteristic. Then $\eta$, defined in
(2), satisfies properties (i)–(iii) mentioned in the Introduction.

Theorem 1 describes two extreme cases corresponding to $\eta=0$ or 1. The
following proposition further illustrates that any value of $\eta$ between 0
and 1 indeed conveys an idea of how different the distributions are. More
specifically, for common location and scale families, $\eta$ increases as the
“difference” between the distributions becomes “larger”.

######  Proposition 1 (Monotonicity of $\eta$).

Consider $M=2$ and $\mathcal{Z}=\mathbb{R}^{d}$ for $d\geq 1$. Suppose $P_{1}$
is a log-concave333A distribution is called log-concave if its density can be
written as $f(x)=\exp\left(-h(x)\right)$ for some convex function
$h:\mathbb{R}^{d}\to\mathbb{R}\cup\\{+\infty\\}$. Many common probability
distributions are log-concave, such as Gaussian distribution, uniform
distribution over a convex set, and gamma distribution if the shape parameter
is $\geq 1$. distribution with density $f(\cdot)$, and
$\pi_{1},\pi_{2}\in(0,1)$ are fixed.

  1. 1.

(Location family) Suppose $P_{1}$ and $P_{2}$ have densities $f(\cdot)$ and
$f(\cdot-\theta)$ respectively, where $\theta=\lambda h$ for a fixed
$h\in\mathbb{R}^{d}\backslash\\{\mathbf{0}\\}$ and $\lambda\geq 0$. Then
$\eta(P_{1},P_{2})$ is a function of $\lambda$. Moreover, $\eta(P_{1},P_{2})$
monotonically increases from 0 to 1 as $\lambda$ increases from 0 to $\infty$.

  2. 2.

(Scale family) Suppose $P_{1}$ and $P_{2}$ have densities $f(\cdot)$ and
$\lambda f(\lambda\times\cdot)$ respectively, where $\lambda>0$. Suppose
further that $f$ is twice differentiable in the interior of its support. Then
$\eta(P_{1},P_{2})$ monotonically decreases from 1 to 0 as $\lambda$ increases
from 0 to 1, and monotonically increases from 0 to 1 as $\lambda$ grows from 1
to $+\infty$.

Further, the above monotonicity of $\eta$ for scenarios 1 and 2 is strict
(i.e., strictly increasing and strictly decreasing) when $\eta(P_{1},P_{2})<1$
(i.e., when $P_{1}$ and $P_{2}$ are not mutually singular).

In the following, we state a few important properties of $\eta$. These
properties are commonly associated with $f$-divergences [45]. Our first
result, Proposition 2 (proved in Appendix C.3), shows that “processing” the
$M$ distributions makes them “less different”, as measured by $\eta$.

######  Proposition 2 (Data processing inequality and invariance).

Recall $\tilde{\Delta}$ and $\tilde{Z}\in\mathcal{Z}$ defined at the beginning
of Section 2.1. Let $\mathcal{Z}^{\prime}$ be another measurable space, and
$\kappa(\cdot,\cdot)$ be a transition kernel from $\mathcal{Z}$ to
$\mathcal{Z}^{\prime}$, i.e., for any $z\in\mathcal{Z}$, $\kappa(z,\cdot)$
specifies a distribution on $\mathcal{Z}^{\prime}$. Suppose $P_{i}$ is
transitioned to $Q_{i}$ by $\kappa$, i.e., $Q_{i}(B):=\int\kappa(z,B){\rm
d}P_{i}(z)$, for every measurable set $B\subset\mathcal{Z}^{\prime}$, for
$i=1,\ldots,M$. If the mixture proportions $\\{\pi_{i}\\}_{i=1}^{M}$ are held
fixed, then

$\eta(P_{1},\ldots,P_{M})\geq\eta(Q_{1},\ldots,Q_{M}).$

Further, equality holds in the above display if and only if
$\tilde{\Delta}\mid\tilde{Z}\stackrel{{\scriptstyle
d}}{{=}}\tilde{\Delta}\mid\tilde{Z}^{\prime}$ where $\tilde{Z}^{\prime}$ is
obtained by passing $\tilde{Z}$ through the transition kernel $\kappa$. In
particular, $\eta$ is invariant under any measurable bijective transformation.

The above result has an interesting consequence. Suppose that $T(Y_{i})\sim
Q_{i}$ where $Y_{i}\sim P_{i}$ and $T:\mathcal{Z}\to\mathcal{Z}^{\prime}$ is a
bijection. Then Proposition 2 implies that
$\eta(P_{1},\ldots,P_{M})=\eta(Q_{1},\ldots,Q_{M})$, thereby showing that a
bijective transformation of the $M$ distributions does not change our measure
of dissimilarity $\eta$.

Many distance measures, including $f$-divergences, satisfy a convexity
property [45]. The following result (proved in Appendix C.4) shows that $\eta$
is also jointly convex in its inputs, i.e., the $M$ distributions.

######  Proposition 3 (Joint convexity of $\eta$).

Let $P_{1},\ldots,P_{M}$, $Q_{1},\ldots,Q_{M}$ be distributions on
$\mathcal{Z}$, and $\lambda\in[0,1]$. If the mixture proportions
$\\{\pi_{i}\\}_{i=1}^{M}$ are held fixed, then

$\eta\left(\lambda P_{1}+(1-\lambda)Q_{1},\ldots,\lambda
P_{M}+(1-\lambda)Q_{M}\right)\leq\lambda\eta(P_{1},\ldots,P_{M})+(1-\lambda)\eta(Q_{1},\ldots,Q_{M}).$

Combining the above convexity result with Jensen’s inequality yields the
following corollary of Proposition 3.

######  Corollary 1 (Conditioning increases $\eta$).

Let $Y_{i}\sim P_{i}$, for $i=1,\ldots,M$, and $W$ be any random variable.
With the mixture proportions $\\{\pi_{i}\\}_{i=1}^{M}$ held fixed, we have:

$\mathbb{E}_{W}\left[\eta\left(Y_{1}|W,\ldots,Y_{M}|W\right)\right]\geq\eta\left(P_{1},\ldots,P_{M}\right).$

The above corollary has an information theoretic interpretation: If we view
$Y_{i}\sim P_{i}$ as the output distribution after passing $W$ through a
“channel” $p(Y_{i}|W)$ [19], the above relation tells us that the average
“difference” between the corresponding channel transitions is at least the
“difference” between the output distributions.

######  Remark 1 (Connection to Henze-Penrose dissimilarity [41]).

In fact, $\eta$, in the special case when $\mathcal{Z}=\mathbb{R}^{d}$ and
$M=2$, has a close connection to existing measures of dissimilarity between
distributions. Suppose that $P_{1}$ and $P_{2}$ have densities $f$ and $g$
w.r.t. the Lebesgue measure on $\mathbb{R}^{d}$, then, it can be shown that
$\eta=1-\int\frac{f(x)g(x)}{\pi_{1}f(x)+\pi_{2}g(x)}\mathrm{d}x$ (see Appendix
C.5 for a proof). This has a close connection to the Henze-Penrose
dissimilarity [41] defined as
$\delta(f,g):=1-2\pi_{1}\pi_{2}\int\frac{f(x)g(x)}{\pi_{1}f(x)+\pi_{2}g(x)}\mathrm{d}x,$
which belongs to a general class of separation measures between distributions
[34], and arises as the limit of a large class of graph-based two-sample
tests, including the Friedman-Rafsky test [28, 41], nearest-neighbor based
tests [63, 40], and the crossmatch test [2, 49].

##  3 Estimation

In this section the sample version $\hat{\eta}$ of $\eta$ will be introduced.
While the definition of $\eta$ does not require $\mathcal{Z}$ to be a metric
space, to establish some useful properties of $\hat{\eta}$ we will assume that
$\mathcal{Z}$ is a metric space with distance function
$\rho_{\mathcal{Z}}(\cdot,\cdot)$. The notion of a _geometric graph_ [21, 10]
is crucial in our construction. Intuitively, in a geometric graph on a set of
points $z_{1},\ldots,z_{m}\in\mathcal{Z}$ ($m\geq 1$), an edge $(z_{i},z_{j})$
will appear if $z_{i}$ and $z_{j}$ are “close” in distance.

Formally, $\mathcal{G}$ is said to be a geometric graph on $\mathcal{Z}$ if,
given any finite subset $S$ of $\mathcal{Z}$, $\mathcal{G}(S)$ defines a graph
with vertex set $S$ and the corresponding edge set
$\mathcal{E}(\mathcal{G}(S))$. The graph can be both directed or undirected,
and we will restrict ourselves to simple graphs, i.e., graphs without multiple
edges and self loops. Examples of such graphs include minimum spanning trees
(MSTs) and $k$-nearest neighbor ($k$-NN) graphs (where $k\geq 1$), as
described below.

  1. 1.

$k$-NN graph: The directed $k$-NN graph puts an edge from each node $z_{i}$ to
its $k$-NNs among $z_{1},\ldots,z_{i-1},z_{i+1},\ldots,z_{m}$ (so $z_{i}$ is
excluded from the set of its $k$-NNs). Ties are broken at random if they occur
to ensure the out-degree is always $k$. The undirected $k$-NN graph is
obtained by ignoring the direction in the directed $k$-NN graph and removing
multiple edges if they exist.

  2. 2.

MST: An MST is a subset of edges of an edge-weighted undirected graph which
connects all the vertices in the graph with the least possible sum of edge
weights and contains no cycles. For instance, given the set of points
$z_{1},\ldots,z_{m}\in\mathcal{Z}$ one can construct an MST for the complete
graph with vertices as $z_{i}$’s and edge weights being the distance
$\rho_{\mathcal{Z}}(z_{i},z_{j})$ between vertices $z_{i}$ and $z_{j}$.

In practice, a $k$-NN graph is recommended as the primary choice over the MST
for its flexibility and computational convenience (see Remark 3 below).

###  3.1 Definition of $\hat{\eta}$

Recall our $M$ sample problem setting (in Section 2.1) and let
$Z_{1},\ldots,Z_{n}$ be the pooled sample and $\Delta_{1},\ldots,\Delta_{n}$
be the corresponding observation labels. Let
$\mathcal{G}_{n}:=\mathcal{G}(Z_{1},\ldots,Z_{n})$ where $\mathcal{G}$ is a
geometric graph on $\mathcal{Z}$ such that
$(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})$ implies $Z_{i}$ and $Z_{j}$ are
“close”. Let $d_{i}$ be the out-degree of $Z_{i}$ in $\mathcal{G}_{n}$. We
consider the following estimator of $\eta$:

$\displaystyle\hat{\eta}:=\frac{\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})-\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K(\Delta_{i},\Delta_{j})}{\frac{1}{n}\sum_{i=1}^{n}K(\Delta_{i},\Delta_{i})-\frac{1}{n(n-1)}\sum_{i\neq
j}K(\Delta_{i},\Delta_{j})}.$ (4)

The definition of this estimator is intuitive: first,
$\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta})\big{]}$ in the denominator
of $\eta$ (see (2)) is estimated by
$\frac{1}{n}\sum_{i=1}^{n}K(\Delta_{i},\Delta_{i})$; second,
$\mathbb{E}\big{[}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\big{]}$ is
estimated by

$\frac{\sum_{i\neq
j}K(\Delta_{i},\Delta_{j})}{n(n-1)}=\frac{\sum_{i=1}^{M}n_{i}(n_{i}-1)K(i,i)+\sum_{i\neq
j}^{M}n_{i}n_{j}K(i,j)}{n(n-1)}.$ (5)

For estimating the remaining term
$\mathbb{E}\left[\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|\tilde{Z}\big{]}\right]$,
ideally we would want two independent observations
$\tilde{\Delta},\tilde{\Delta}{{}^{\prime}}$ from the conditional distribution
$\tilde{\Delta}|\tilde{Z}=Z_{i}$, so we take one to be $\Delta_{i}$, and the
other to be the label $\Delta_{j}$ of an observation $Z_{j}$ which is ‘close’
to $Z_{i}$. The first term in the numerator of (4) formalizes this intuition
via the geometric graph $\mathcal{G}_{n}$. Our estimator $\hat{\eta}$ also has
a nice interpretation: it can be shown that $\hat{\eta}$ is linearly related
to the leave-one-out cross-validation accuracy of a $k$-NN classifier if the
discrete kernel is used; see Appendix A.3 for the details.

With a $k$-NN graph, $\hat{\eta}$ has near linear computational complexity
$O(kn\log n)$. Note that finding the $k$-NN graph has computational complexity
$O(kn\log n)$. Further, the first term in the numerator of (4) is a sum of
$kn$ terms (as $d_{i}\equiv k$ for all $i$); the second term in the numerator
of (4) can in fact be computed in $O(1)$ time (see (5)).

The following theorem (proved in Appendix C.7) generalizes the above special
case and shows that $\hat{\eta}$ is strongly consistent in estimating $\eta$
under mild assumptions.

######  Theorem 2 (Consistency).

Recall our setting as described at the beginning of Section 2.1 where we
assume that $\frac{n_{i}}{n}\to\pi_{i}\in(0,1)$ as $n\to\infty$, for
$i=1,\ldots,M$. Suppose that
$\tilde{\mathcal{G}}_{n}:=\mathcal{G}(\tilde{Z}_{1},\tilde{Z}_{2},\ldots,\tilde{Z}_{n})$,
where $\tilde{Z}_{1},\tilde{Z}_{2},\ldots,\tilde{Z}_{n}$ are i.i.d. from the
mixture distribution $\sum_{i=1}^{M}\pi_{i}P_{i}$ on $\mathcal{Z}$, satisfies
Assumptions 1–3 (detailed in Appendix A.1). Then,
$\hat{\eta}\overset{a.s.}{\longrightarrow}\eta$, as $n\to\infty$.

Assumptions 1–3 required on the geometric graph $\tilde{\mathcal{G}}_{n}$ for
the above result were made in Deb et al. [21, Theorem 3.1]. For the $k$-NN
graph and the MST, these conditions are satisfied under mild assumptions. For
example, in an Euclidean space, they hold for the $k$-NN graph when
$\|\tilde{Z}_{1}-\tilde{Z}_{2}\|$ has a continuous distribution and
$k=o(\frac{n}{\log n})$; for the MST these are satisfied when $\tilde{Z}_{1}$
has an absolutely continuous distribution [21, Proposition 3.2].

The following remark shows that $\hat{\eta}$ in (4) generalizes, in various
directions, many previously known nearest-neighbor type statistics for two-
sample testing.

######  Remark 2 (Connection to [63] and [40]).

Under a Euclidean setting with $M=2$, if the discrete kernel and the directed
$k$-NN graph (with fixed $k$) are used, then from (4),

$\hat{\eta}=\frac{\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}I(\Delta_{i}=\Delta_{j})-\frac{1}{n(n-1)}\sum\limits_{i=1}^{M}n_{i}(n_{i}-1)}{1-\frac{1}{n(n-1)}\sum\limits_{i=1}^{M}n_{i}(n_{i}-1)}.$

A linearly transformed version of this statistic was first proposed in
Schilling [63] for two-sample testing. Its detailed asymptotic properties were
later analyzed by Henze [40]. Our proposed $\hat{\eta}$ in this paper can be
viewed as a normalized version of this statistic, allowing $M>2$ distributions
defined on general metric spaces besides the Euclidean space, using more
general kernel functions in addition to the discrete kernel, and other
geometric graphs besides the $k$-NN graph. Even for the directed $k$-NN graph,
$k$ is allowed to be unbounded and grow with $n$ in the analysis of our
consistency result (Theorem 2) and the CLT (see Theorem 3 below), instead of
being fixed as in Schilling [63], Henze [40], Bhattacharya [11].

##  4 Asymptotic Behavior of $\hat{\eta}$ under ${\rm H}_{0}$

In this section the asymptotic normality of $\hat{\eta}$ under ${\rm H}_{0}$
will be derived, which will yield a simple asymptotic test for the equality of
the $M$ distributions in (1). We show that under the null, the permutation
distribution of $\hat{\eta}$, i.e., the distribution of $\hat{\eta}$ given the
pooled sample $\\{Z_{1},\ldots,Z_{n}\\}$, is asymptotically normal, and as a
result, the unconditional asymptotic distribution of $\hat{\eta}$ is also
normal.

###  4.1 Asymptotic Permutation Distribution

Let $\mathcal{F}_{n}$ be the $\sigma$-algebra generated by the unordered
pooled sample $\\{Z_{1},\ldots,Z_{n}\\}$ without the labeling information
(note that the number of observations from each distribution $n_{i}$, for
$i=1,\ldots,M$, is known). The following result (proved in Appendix C.8)
states that the permutation distribution, i.e., the conditional distribution
of $\hat{\eta}$ given $\mathcal{F}_{n}$, is asymptotically normal.
Furthermore, the asymptotic variance of $\hat{\eta}$ under ${\rm H}_{0}$ is
distribution-free under suitable conditions, in the sense that it does not
depend on the underlying common distribution $P_{1}=\ldots=P_{M}$.

######  Theorem 3.

Suppose Assumptions 2 and 3 (in Appendix A.1) are satisfied,
$\frac{n_{i}}{n}\to\pi_{i}\in(0,1)$ as $n\to\infty$, for $i=1,\ldots,M$, the
kernel $K(\cdot,\cdot)$ is characteristic, and the vertex degrees of
$\mathcal{G}_{n}$ (recall that $\mathcal{G}_{n}$ is the geometric graph
constructed on $\\{Z_{1},\ldots,Z_{n}\\}$) are bounded above by $t_{n}$ such
that $\frac{t_{n}^{r}}{n}\to 0$ for all $r\in\mathbb{N}$ (e.g., $t_{n}=C(\log
n)^{\gamma}$ for some $C>0,\gamma\geq 0$). Then, under the null hypothesis (1)
that all the $M$ distributions are equal,

$\frac{\hat{\eta}}{\sqrt{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}}\big{|}\mathcal{F}_{n}\overset{d}{\to}N(0,1),\quad{\rm
where}$ $\frac{\hat{\eta}}{\sqrt{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}}=\frac{\sqrt{n}\left(\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{{d}_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K({\Delta}_{i},{\Delta}_{j})-\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K(\Delta_{i},\Delta_{j})\right)}{\sqrt{\tilde{a}\left(\tilde{g}_{1}+\tilde{g}_{3}-\frac{2}{n-1}\right)+\tilde{b}\left(\tilde{g}_{2}-2\tilde{g}_{1}-2\tilde{g}_{3}-1+\frac{4}{n-1}\right)+\tilde{c}\left(\tilde{g}_{1}-\tilde{g}_{2}+\tilde{g}_{3}+\frac{n-3}{n-1}\right)}},$

with $\tilde{g}_{1}:=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}$,
$\tilde{g}_{2}:=\frac{1}{n}\sum_{i,j=1}^{n}\frac{T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}}$,
$T^{\mathcal{G}_{n}}(i,j):=\sum_{k=1}^{n}I\\{(Z_{i},Z_{k}),(Z_{j},Z_{k})\in\mathcal{E}(\mathcal{G}_{n})\\}$
being the number of common out-neighbors,
$\tilde{g}_{3}:=\frac{1}{n}{\sum_{i,j:(Z_{i},Z_{j}),(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}}\frac{1}{d_{i}d_{j}}$,
and

$\displaystyle\tilde{a}$
$\displaystyle:=\frac{1}{n(n-1)}{\sum}^{\prime}K^{2}(\Delta_{i},\Delta_{j}),\qquad\quad\tilde{b}\;:=\frac{1}{n(n-1)(n-2)}{\sum}^{\prime}K(\Delta_{i},\Delta_{j})K(\Delta_{i},\Delta_{l}),$
(6) $\displaystyle\tilde{c}$
$\displaystyle:=\frac{1}{n(n-1)(n-2)(n-3)}{\sum}^{\prime}K(\Delta_{i},\Delta_{j})K(\Delta_{l},\Delta_{m}).$

Here $\sum^{\prime}$ means the summation is over distinct indices. Further,
${\rm Var}(\hat{\eta}|\mathcal{F}_{n})=O\left(\frac{1}{n}\right)$. From the
above conditional CLT, the unconditional CLT also follows:
$\frac{\hat{\eta}}{\sqrt{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}}\overset{d}{\to}N(0,1).$

Although the statement of the above theorem is similar to Deb et al. [21,
Theorem 4.1], the proof technique is quite different, since we are dealing
with the conditional distribution given $\mathcal{F}_{n}$ instead of the
unconditional one. The main technical tool used here is a modification of the
CLT in Pham et al. [58], based on moment matching. Such a technique has been
used to show the asymptotics for a variety of graph-based statistics [14, 40,
41, 57]. Our general result (see Theorem C.1 in Appendix C.8), compared to
[58], can deal with unbounded vertex degrees, which may be of independent
interest.

Theorem 3 has the benefit that if we regard the randomness as coming from
random permutations (with the pooled sample fixed), then the variance ${\rm
Var}(\hat{\eta}|\mathcal{F}_{n})$ can be computed exactly. This could lead to
a better approximation of the sampling distribution of $\hat{\eta}$, compared
to the test that uses the limiting value of ${\rm
Var}(\hat{\eta}|\mathcal{F}_{n})$ (see Section 4.2 below). Note that we allow
$k$, the number of nearest neighbors, to be unbounded and grow with $n$
instead of being fixed as in some previous relevant works [63, 40, 57].

######  Corollary 2.

Consider the testing problem (1). Under the assumptions of Theorem 3, if
furthermore Assumption 1 (in Appendix A.1) holds, then the $M$-sample test
with rejection region:

$\frac{\hat{\eta}}{\sqrt{{\rm Var}(\hat{\eta}|\mathcal{F}_{n})}}\geq
z_{1-\alpha},$ (7)

where $z_{\alpha}$ is the $\alpha^{\rm th}$ quantile of the standard normal
distribution, has asymptotic level $\alpha\in(0,1)$ and is consistent against
any alternative where at least two distributions are different.

The above result is a direct consequence of Theorem 3 and the fact that
$\hat{\eta}$ converges in probability to $\eta>0$ under any alternative (by
Theorem 2).

######  Remark 3 (Computational efficiency).

Another appealing property of our method is that it is computationally
efficient and easy to implement. When a Euclidean $k$-NN graph is used, the
computation complexity for $\hat{\eta}$ and $\frac{\hat{\eta}}{\sqrt{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}}$ is $O(kn\log n)$ (see Appendix A.2 for
more detailed implementation schemes; in particular, $\tilde{a}$, $\tilde{b}$,
$\tilde{c}$ can be computed in $O(1)$ time given $n_{1},\ldots,n_{M}$), which
is near linear time when $k$ does not grow too fast.

###  4.2 Asymptotic Distribution-Free Property

We will show in this sub-section that the asymptotic distribution of
$\hat{\eta}$ under ${\rm H}_{0}$ (see (1)) is distribution-free under mild
assumptions. Such a property actually holds true for a variety of graph-based
test statistics [40, 41, 62], provided that the null distribution is
absolutely continuous and the geometric graph $\mathcal{G}$ is “local” in the
sense of a stabilizing graph [56, 5, 11], which includes MST and $k$-NN graphs
(with fixed $k$) [56].

From Theorem 3 above, it is clear that the variance of $\hat{\eta}$ under
${\rm H}_{0}$ involves $\tilde{a}$, $\tilde{b}$, $\tilde{c}$, $\tilde{g}_{1}$,
$\tilde{g}_{2}$, and $\tilde{g}_{3}$. From (6), it is easy to see that
$\tilde{a}$, $\tilde{b}$, $\tilde{c}$ converge to some $a$, $b$, $c$ depending
only on the kernel and the mixture proportions $\\{\pi_{i}\\}_{i=1}^{M}$. It
can also be shown that $\tilde{g}_{1}$, $\tilde{g}_{2}$, $\tilde{g}_{3}$
converge to some $g_{1}$, $g_{2}$, $g_{3}$ respectively. To motivate these
limits, take $\tilde{g}_{1}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}$ as an
example. Under the null, $Z_{1},\ldots,Z_{n}$ are i.i.d. from a common
distribution. In a “local” graph, $d_{i}$, the out-degree of $Z_{i}$, may only
depend on the points near $Z_{i}$. If the common distribution has a density
continuous at $Z_{i}$, then the points near $Z_{i}$ are approximately sampled
from a distribution with constant density. This means that if the geometric
graph is translation and scale invariant (defined formally later), then
$\frac{1}{d_{i}}$ may be similar in distribution to
$\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}$, where
$\mathcal{P}_{1}$ is the homogeneous Poisson process444For a homogeneous
Poisson process on $\mathbb{R}^{d}$ with intensity $\lambda$ (denoted by
$\mathcal{P}_{\lambda}$), the number of points in a set
$A\subset\mathbb{R}^{d}$ follows a ${\rm Poisson}(\int_{A}\lambda{\rm d}x)$
distribution, and the number of points in disjoint sets are independent. of
intensity 1 on $\mathbb{R}^{d}$, $\mathcal{P}_{1}^{\mathbf{0}}$ is
$\mathcal{P}_{1}\cup\\{\mathbf{0}\\}$, and
$d\big{(}\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\big{)}$ is the
out-degree of $\mathbf{0}$ in $\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})$ —
the geometric graph constructed on $\mathcal{P}_{1}^{\mathbf{0}}$. Hence, it
is reasonable to expect that
$\tilde{g}_{1}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}$ converges to
$g_{1}:=\mathbb{E}\left[\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]$,
which is distribution-free.

In the following, we introduce the necessary mathematical concepts in order to
state our formal results. A geometric graph $\mathcal{G}$ on $\mathbb{R}^{d}$
is said to be translation invariant if translation by $x$ induces a graph
isomorphism555That is, for $x_{1},x_{2}\in\mathcal{X}$, $(x_{1},x_{2})$ is an
edge in $\mathcal{G}(\mathcal{X})$ if and only if $(x_{1}+x,x_{2}+x)$ is an
edge in $\mathcal{G}(\mathcal{X}+x)$. from $\mathcal{G}(\mathcal{X})$ to
$\mathcal{G}(x+\mathcal{X})$ for all $x\in\mathbb{R}^{d}$ and all finite
$\mathcal{X}\subset\mathbb{R}^{d}$. Similarly, $\mathcal{G}$ is scale
invariant if scalar multiplication by $a$ induces a graph isomorphism from
$\mathcal{G}(\mathcal{X})$ to $\mathcal{G}(a\mathcal{X})$ for all finite
$\mathcal{X}$ and all $a>0$. For $\lambda\geq 0$, denote by
$\mathcal{P}_{\lambda}$ the homogeneous Poisson process of intensity $\lambda$
in $\mathbb{R}^{d}$, and define
$\mathcal{P}_{\lambda}^{x}:=\mathcal{P}_{\lambda}\cup\\{x\\}$. $\mathcal{G}$
is said to be stabilizing on $\mathcal{P}_{\lambda}$ if, for almost all
realizations $\mathcal{P}_{\lambda}$, there exists a random variable
$R<\infty$ such that the set of edges incident at the origin is not changed by
modifying points outside a ball of radius $R$, i.e.,
$\mathcal{E}\left(\mathbf{0},\mathcal{G}\big{(}(\mathcal{P}_{\lambda}^{\mathbf{0}}\cap
B(\mathbf{0},R))\cup\mathcal{A}\big{)}\right)=\mathcal{E}\left(\mathbf{0},\mathcal{G}\big{(}\mathcal{P}_{\lambda}^{\mathbf{0}}\cap
B(\mathbf{0},R)\big{)}\right)$ for all finite
$\mathcal{A}\subset\mathbb{R}^{d}\backslash B(\mathbf{0},R)$, where
$B(\mathbf{0},R)$ is the closed Euclidean ball of radius $R$ centered at the
origin $\mathbf{0}\in\mathbb{R}^{d}$, and
$\mathcal{E}(x,\mathcal{G}(\mathcal{X}))=\\{(x,y):(x,y){\rm\ is\ an\ edge\ in\
}\mathcal{G}(\mathcal{X})\\}$ is the set of edges of
$\mathcal{G}(\mathcal{X})$ incident666Note that in a directed graph,
$\mathcal{E}(x,\mathcal{G}(\mathcal{X}))$ only includes the edges starting
from $x$. to $x\in\mathcal{X}$. In such a case, the definition of
$\mathcal{G}$ can be extended to the infinite point set
$\mathcal{P}_{\lambda}^{\mathbf{0}}$, with
$\mathcal{E}\left(\mathbf{0},\mathcal{G}(\mathcal{P}_{\lambda}^{\mathbf{0}})\right):=\mathcal{E}\left(\mathbf{0},\mathcal{G}\big{(}\mathcal{P}_{\lambda}^{\mathbf{0}}\cap
B(\mathbf{0},R)\big{)}\right)$. It is known that both MST and $k$-NN graphs
(with fixed $k$) are translation and scale invariant, and stabilizing on
$\mathcal{P}_{\lambda}$ for all $\lambda>0$; see e.g., [56]. The following
result (proved in Appendix C.10) formally states the asymptotic distribution-
free property of $\hat{\eta}$.

######  Theorem 4.

Under the same assumptions as in Theorem 3, if furthermore $\mathcal{G}$ is
translation and scale invariant, and stabilizing on $\mathcal{P}_{\lambda}$
for some $\lambda>0$, then under the null hypothesis (1) where
$P_{1}=\ldots=P_{M}$ is assumed to have a Lebesgue density on
$\mathbb{R}^{d}$, we have:

$\sqrt{n}\hat{\eta}\overset{d}{\to}N\left(0,\sigma^{2}_{\mathcal{G},K,d,\pi}\right),$

where $\sigma^{2}_{\mathcal{G},K,d,\pi}$ is a positive constant not depending
on the common density. More specifically,

$\sigma^{2}_{\mathcal{G},K,d,\pi}:=\frac{{a}\left({g}_{1}+{g}_{3}\right)+{b}\left({g}_{2}-2{g}_{1}-2{g}_{3}-1\right)+{c}\left({g}_{1}-{g}_{2}+{g}_{3}+1\right)}{\left(\sum_{i=1}^{M}\pi_{i}K(i,i)-\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)\right)^{2}}>0,\quad{\rm
where}$ (8) $\displaystyle a$
$\displaystyle:=\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K^{2}(i,j),\qquad
b:=\sum_{i,j,l=1}^{M}\pi_{i}\pi_{j}\pi_{l}K(i,j)K(i,l),\qquad
c:=\left(\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)\right)^{2},$ $\displaystyle
g_{1}$
$\displaystyle:=\mathbb{E}\left[\frac{1}{d\big{(}\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\big{)}}\right],\quad\
g_{2}:=\mathbb{E}\left[\sum_{y\neq
z:(y,\mathbf{0}),(z,\mathbf{0})\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{1}{d\left(y,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)d\left(z,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right]+g_{1},$
$\displaystyle g_{3}$
$\displaystyle:=\mathbb{E}\left[\sum_{y:(y,\mathbf{0}),(\mathbf{0},y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{1}{d\left(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)d\left(y,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right].$

The proof of the above result proceeds by showing that the permutation
variance in Theorem 3 converges to the distribution-free limiting variance in
Theorem 4. Compared to Theorem 4, Theorem 3 is more general in the sense that
it does not require a Euclidean space. Even in a Euclidean space, Theorem 3 is
based on the exact permutation variance which may yield a more accurate
approximation of the sampling distribution of $\hat{\eta}$ than Theorem 4,
especially when the dimension is high, where the convergence of
$\tilde{g}_{i}$, for $i=1,2,3$, may be slower; see Appendix D.4 for empirical
evidence.

##  5 Asymptotic Power and Detection Threshold

We have already seen that $\hat{\eta}$ is asymptotically normal under the
null. A natural question to ask now is: “what is the asymptotic distribution
of $\hat{\eta}$ under alternatives?”. It turns out that under a fixed
alternative or certain shrinking alternatives converging to the null,
$\hat{\eta}$ is also asymptotically normal, if it is properly centered. Using
these CLTs, we can provide a complete characterization of the asymptotic power
of the test (7). We will focus on $\mathcal{Z}=\mathbb{R}^{d}$ in this
section, and assume that $P_{i}$ has a density $f_{i}$ w.r.t. the Lebesgue
measure.

The asymptotic distribution under alternatives and the detection threshold of
graph-based statistics were not available until the recent work by
Bhattacharya [11], where the analysis was carried out under a Poissonized
setting — instead of assuming $\frac{n_{i}}{n}\to\pi_{i}$ as $n\to\infty$, it
is assumed that $n_{i}\sim{\rm Poisson}(N_{i})$, and
$\frac{N_{i}}{N}\to\pi_{i}>0$ as $N\to\infty$ where $N:=\sum_{i=1}^{M}N_{i}$.
In the Poissonized framework, the pooled data $\\{Z_{1},\ldots,Z_{n}\\}$
follows a non-homogeneous Poisson point process. The spatial independence of
this Poisson process facilitates the computation of the variance of
$\hat{\eta}$ under alternatives, and helps us establish CLTs of $\hat{\eta}$
using Stein’s method [18]. We consider this Poissonized setting in this
section.

CLT under a general fixed alternative. Under ${\rm H}_{0}$, $\hat{\eta}$ is
unbiased in estimating $\eta\equiv 0$, so it is reasonable to conjecture that
$\sqrt{n}(\hat{\eta}-\eta)$ would converge to a normal distribution (as shown
in Theorem 3). However, this is not the case under a fixed alternative, where
the bias of $\hat{\eta}$ can be large and dominating. Although $\hat{\eta}$
converges to $\eta$ under any fixed alternative (Theorem 2), the distance
between $\hat{\eta}$ and $\eta$ under such an alternative is in general of
order $n^{-1/d}$ [21, Corollary 5.1], having the same order as the distance
between a data point and its nearest neighbor. Hence in such a case, it is
unrealistic to expect that $\sqrt{n}(\hat{\eta}-\eta)$ would be asymptotically
normal. Instead, we show in Theorem B.1 in Appendix B that $\hat{\eta}$ is
asymptotically normal after appropriate centering.

CLT under shrinking alternatives. Given that our level $\alpha$ test (7) is
consistent against all fixed alternatives, it is natural to study its power
behavior under shrinking alternatives converging to the null. For this
purpose, a CLT under shrinking alternatives is needed, which is shown in
Theorem B.2 in Appendix B under the setting where the $i$-th distribution has
a Lebesgue density $f_{i}^{N}$ converging Lebesgue almost everywhere to some
density $f$ as $N\to\infty$, for $i=1,\ldots,M$.

Asymptotic power and detection threshold. Here we answer the following
question: “Can we characterize the exact limiting power of our test (7) under
a sequence of shrinking alternatives?”. To answer the above question, we
consider $M=2$ and study the power of the test along a parametric sub-model
[11]: suppose $P_{i}$ has Lebesgue density $f(\cdot|\theta_{i})$ on
$\mathbb{R}^{d}$ belonging to the parametric family
$\\{f(\cdot|\theta):\theta\in\Theta\subset\mathbb{R}^{p}\\}$ ($p\geq 1$), for
$i=1,2$.

Theorem 5 below (proved in Appendix C.14) describes the asymptotic power of
our test (7), whose proof depends on the CLT under shrinking alternatives. As
a consequence of Theorem 5, we will also be able to answer the closely related
question: “At what rate should $\theta_{2}$ converge to $\theta_{1}$ (assumed
fixed) so that our test (7) would be powerless (i.e., power $\leq\alpha$) if
the convergence is faster than the rate, and would have asymptotic power 1 if
the convergence is slower than the rate, as the sample size increases?”. In
such a situation, the order of $|\theta_{2}-\theta_{1}|$ as a function of the
sample size $N$ is called the detection threshold of the test [11].

######  Theorem 5.

Suppose $\\{\mathbb{P}_{\theta}\\}_{\theta\in\Theta}$ is a parametric family
of distributions with a convex parameter space $\Theta\subset\mathbb{R}^{p}$;
we further assume that $\mathbb{P}_{\theta}$ has a Lebesgue density
$f(\cdot|\theta)$. Suppose that $M=2$ and the discrete kernel $K$ and the
$k$-NN graph (for a fixed $k\geq 1$) are used in defining $\hat{\eta}$. Let
$P_{1}$ and $P_{2}$ have densities $f(\cdot|\theta_{1})$ and
$f(\cdot|\theta_{2})$ respectively, and
$\theta_{2}=\theta_{1}+\varepsilon_{N}$, where $\varepsilon_{N}\to 0$ as
$N\to\infty$. Let $\sigma^{2}_{\mathcal{G},K,d,\pi}$ be as defined in Theorem
4, and ${\rm H}_{x}f(x|\theta_{1})$ be the Hessian of $f(x|\theta_{1})$ (taken
w.r.t. $x$). Also, define, for $h\in\mathbb{R}^{p}$,

$a_{k,\theta_{1}}(h):=-\frac{(1-2\pi_{2})C_{k,2}}{4\,k\,d\,\sigma_{\mathcal{G},K,d,\pi}}\int
h^{\top}\nabla_{\theta_{1}}\left(\frac{{\rm tr}({\rm
H}_{x}f(x|\theta_{1}))}{f(x|\theta_{1})}\right)f^{\frac{d-2}{d}}(x|\theta_{1}){\rm
d}x,$ (9)
$b_{k,\theta_{1}}(h):=\frac{\pi_{1}\pi_{2}}{\sigma_{\mathcal{G},K,d,\pi}}\mathbb{E}_{X\sim
f(\cdot|\theta_{1})}\left[\frac{h^{\top}\nabla_{\theta_{1}}f(X|\theta_{1})}{f(X|\theta_{1})}\right]^{2},\;{\rm
with}\;\;C_{k,2}:=\mathbb{E}\Bigg{[}\sum_{x\in\mathcal{P}_{1}^{\mathbf{0}}:(\mathbf{0},x)\in\mathcal{E}(\mathcal{G}_{k{\rm
NN}}(\mathcal{P}_{1}^{\mathbf{0}}))}\|x\|^{2}\Bigg{]}.$

Under suitable assumptions on (i) the smoothness of the parametric family,
(ii) conditions such that the $k$-NN graph is nicely behaved (see Appendix B.3
for the detailed list of assumptions), and (iii)
$\sqrt{N}\left(\frac{N_{i}}{N}-\pi_{i}\right)\to 0$, as $N\to\infty$, for
$i=1,\ldots,M$, we have the following result.

  * 1.

If $d\leq 8$, then the following hold:

    * (a)

when $\|N^{\frac{1}{4}}\varepsilon_{N}\|\to 0$: the limiting power of the test
(7) is $\alpha$.

    * (b)

when $N^{\frac{1}{4}}\varepsilon_{N}\to h\in\mathbb{R}^{d}$: if $d\leq 7$, the
limiting power of (7) is $\Phi\left(z_{\alpha}+b_{k,\theta_{1}}(h)\right)$; if
$d=8$ the limiting power is
$\Phi\left(z_{\alpha}+a_{k,\theta_{1}}(h)+b_{k,\theta_{1}}(h)\right)$.

    * (c)

when $\|N^{\frac{1}{4}}\varepsilon_{N}\|\to\infty$: the limiting power of the
test (7) is 1.

  * 2.

If $d\geq 9$, then the following hold:

    * (a)

when $\|N^{\frac{1}{2}-\frac{2}{d}}\varepsilon_{N}\|\to 0$: the limiting power
of the test (7) is $\alpha$.

    * (b)

when $N^{\frac{1}{2}-\frac{2}{d}}\varepsilon_{N}\to h\in\mathbb{R}^{d}$: the
limiting power of the test (7) is
$\Phi\left(z_{\alpha}+a_{k,\theta_{1}}(h)\right)$.

    * (c)

when $\|N^{\frac{1}{2}-\frac{2}{d}}\varepsilon_{N}\|\to\infty$ such that
$\|N^{\frac{2}{d}}\varepsilon_{N}\|\to 0$: then depending on whether

$N^{\frac{1}{2}-\frac{2}{d}}(1-2\pi_{2})\int\varepsilon_{N}^{\top}\nabla_{\theta_{1}}\left(\frac{{\rm
tr}({\rm
H}_{x}f(x|\theta_{1}))}{f(x|\theta_{1})}\right)f^{\frac{d-2}{d}}(x|\theta_{1}){\rm
d}x\to\left\\{\begin{aligned} &\infty,\\\ &-\infty,\\\ \end{aligned}\right.$
(10)

the limiting power of the test (7) is 0 or 1, respectively.

    * (d)

when $N^{\frac{2}{d}}\varepsilon_{N}\to h\in\mathbb{R}^{d}$: the limiting
power of the test (7) is 0 or 1, depending on whether
$a_{k,\theta_{1}}(h)+b_{k,\theta_{1}}(h)$ is negative or positive,
respectively.

    * (e)

$N^{\frac{2}{d}}\varepsilon_{N}\to\infty$: the limiting power of the test (7)
is 1.

The above result shows that the detection threshold of $\hat{\eta}$ exhibits a
“$d=8$” phenomenon (see Bhattacharya [11]): when $d\leq 8$, the detection
threshold is $N^{-\frac{1}{4}}$; while when $d\geq 9$, the detection threshold
is somewhere between $N^{-\frac{2}{d}}$ and $N^{-\frac{1}{2}+\frac{2}{d}}$,
depending on the direction of $\varepsilon_{N}$ and the sign of
$a_{k,\theta_{1}}$; see (9) and (10). If $\varepsilon_{N}=\delta_{N}h$ for
$\delta_{N}>0$ and some fixed $h\in\mathbb{R}^{d}$, and
$a_{k,\theta_{1}}(h)>0$, then the detection threshold is $\delta_{N}\sim
N^{-\frac{1}{2}+\frac{2}{d}}$; on the other hand, if $a_{k,\theta_{1}}(h)<0$,
then the detection threshold is $\delta_{N}\sim N^{-\frac{2}{d}}$. When
$a_{k,\theta_{1}}(h)=0$, the precise location of the detection threshold has
to be determined on a case by case basis (see [11]).

The detection threshold for some particular choices of $P_{1}$ and $P_{2}$ are
already known: for distinguishing two truncated normal distributions with
different location parameters the detection threshold is $N^{-\frac{1}{4}}$
for all $d$, while for distinguishing two truncated normals with different
scale parameters the detection threshold can attain both $N^{-\frac{2}{d}}$
and $N^{-\frac{1}{2}+\frac{2}{d}}$ for $d\geq 9$ depending on the sign of
$a_{k,\theta_{1}}(h)$ (see [11, Section 4.2]). Thus in certain cases, this
nonparametric detection threshold can be very close to the parametric rate
$N^{-\frac{1}{2}}$.

In Appendix D.2, we empirically illustrate that the same detection threshold
also holds for the non-Poissonized setting, i.e., the original setting of
Section 2 where the $n_{i}$’s are nonrandom constants instead of
$n_{i}\sim{\rm Poisson}(N_{i})$.

Our results crucially use the general framework established in [11]. However,
we fill in some gaps in the original proof of [11, Theorem 4.2], where only
the CLT under a fixed alternative was shown; the CLT under shrinking
alternatives was not explicitly formulated, but was assumed to hold instead.
Our CLTs (Theorems B.1 and B.2) are also applicable to $M\geq 2$
distributions.

##  6 Numerical Studies

In this section, the finite-sample performance of our methods will be
investigated on both real and synthetic data. We compare the power behavior of
the tests based on $\hat{\eta}$ to other competing methods, based on synthetic
data. We also analyze many real data sets and demonstrate the usefulness of
KMD. Further simulation experiments are relegated to Appendix D.

###  6.1 Power Study on Synthetic Data

The empirical power behavior of some special cases of our test statistic has
been partly investigated in Schilling [63] and Petrie [57] as $\hat{\eta}$ is
equivalent to the statistics in these papers when $M=2$ with $k$-NN graphs and
the discrete kernel. Thus, in this subsection, we focus on the case $M>2$. We
consider the following settings:

  1. 1.

Normal location: $P_{1}=N(\mathbf{0},I_{d})$,
$P_{2}=N(0.1\cdot\mathbf{1},I_{d})$, $P_{3}=N(0.2\cdot\mathbf{1},I_{d})$, as
$d$ varies.

  2. 2.

Normal scale: $P_{1}=N(\mathbf{0},I_{d})$, $P_{2}=N(\mathbf{0},1.5\cdot
I_{d})$, $P_{3}=N(\mathbf{0},2\cdot I_{d})$, as $d$ varies.

  3. 3.

$t$-distribution location: $P_{1}$ (on $\mathbb{R}^{16}$) has each coordinate
drawn i.i.d. from $t(1)$ with noncentrality parameter $\delta$; $P_{2}=P_{3}$
has each coordinate drawn i.i.d. from $t(1)$.

  4. 4.

U-shaped scale: $P_{1}=P_{2}$ is a “U-shaped” distribution which is a mixture
of

$\scriptsize N\left(\begin{pmatrix}0\\\ 0\end{pmatrix},\begin{pmatrix}2&0\\\
0&\frac{1}{8}\end{pmatrix}\right),\quad N\left(\begin{pmatrix}-3\\\
1\end{pmatrix},\begin{pmatrix}\frac{1}{2}&-\frac{1}{3}\\\
-\frac{1}{3}&\frac{1}{2}\end{pmatrix}\right),\quad N\left(\begin{pmatrix}3\\\
1\end{pmatrix},\begin{pmatrix}\frac{1}{2}&\frac{1}{3}\\\
\frac{1}{3}&\frac{1}{2}\end{pmatrix}\right)$

with mixing weights $\frac{1}{2},\frac{1}{4},\frac{1}{4}$. $P_{3}$ is a scalar
multiple (with scale $\in[1,1.4]$) of $P_{1}$.

  5. 5.

S-shaped rotation: $Y_{1}\sim P_{1}\equiv P_{2}$ is a “S-shaped” distribution
which is a mixture of

$\scriptsize N\left(\begin{pmatrix}-\frac{9}{2}\\\
-\frac{1}{2}\end{pmatrix},\begin{pmatrix}\frac{3}{2}&-\sqrt{\frac{3}{8}}\\\
-\sqrt{\frac{3}{8}}&1\end{pmatrix}\right),N\left(\begin{pmatrix}0\\\
-\frac{1}{2}\end{pmatrix},\begin{pmatrix}\frac{3}{2}&\sqrt{\frac{3}{8}}\\\
\sqrt{\frac{3}{8}}&1\end{pmatrix}\right),N\left(\begin{pmatrix}\frac{9}{2}\\\
1\end{pmatrix},\begin{pmatrix}\frac{3}{2}&-\sqrt{\frac{3}{8}}\\\
-\sqrt{\frac{3}{8}}&1\end{pmatrix}\right).$

with weights $\frac{1}{3},\frac{1}{3},\frac{1}{3}$. $P_{3}$ is obtained by
multiplying $Y_{1}$ by $\tiny\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}$, for $\theta\in[0,0.1\pi]$.

  6. 6.

Spherically symmetric: $P_{1},P_{2},P_{3}$ are spherically symmetric
distributions777$Y_{1}\sim P_{1}$ has a spherically symmetric distribution
with radial distribution $L_{1}$ if $Y_{1}=L_{1}\times U$, where $U$ is the
uniform distribution over the unit sphere $\mathcal{S}^{d-1}$ in
$\mathbb{R}^{d}$. on $\mathbb{R}^{d}$ with different radial densities. We
assume that $P_{1},P_{2}$ and $P_{3}$ have Uniform[0,1], ${\rm
Beta}(1-\alpha,1+\alpha)$, and ${\rm Beta}(1+\alpha,1-\alpha)$ radial
densities respectively (as $\alpha$ varies).

The first two settings are classical, which were also considered in Mukherjee
et al. [49].

Figure 2: Empirical power curves for different methods under the various
setings.

The rest of the examples go beyond normality, considering general
distributions exhibiting different kinds of dissimilarities. The U-shaped and
S-shaped distributions are taken from Hallin et al. [37]. We set $n_{i}=100$,
for $i=1,2,3$. All tests in the following (except MANOVA() built in R) use 500
random permutations to compute their $p$-values. The null hypothesis is
rejected if the observed $p$-value is less than 0.05, and the power is
estimated by averaging over 1000 replications. We compare our methods with
other parametric and nonparametric tests: ‘Pillai’ and ‘Roy’ are the
parametric tests implemented by MANOVA() in R. ‘MCM’, ‘MMCM’ are the
nonparametric test statistics proposed by Mukherjee et al. [49]. ‘KMD’ is our
method with $k$-NN graph (see (4)). ‘DISCO’ [61] is a nonparametric
generalization of ANOVA using the energy statistics between samples,
implemented as a permutation test in the R package energy [60] with default
index parameter $\alpha=1$.

It can be seen from the top left plot of Figure 2 that DISCO achieves similar
power as the parametric method ‘Pillai’ in the normal location problem. For
this setting, the power curve of our test ‘KMD’ is sandwiched between these
two methods and the nonparametric methods MCM and MMCM (which have the
attractive property of being fully distribution-free in finite samples but
also at a cost of a lower power compared to other methods). Roy’s largest root
test is known to be powerful in this normal location example [23, Section
11.11.1]. When the dimension is larger than the sample size, the error
variance matrix is singular and the parametric tests like MANOVA() cannot be
applied.

For the normal scale problem in Figure 2 we illustrate that the choice of an
appropriate kernel for our method can lead to improved performance; see
‘KMD(K2)’. For this setting, the observations coming from the distributions
with the larger scale parameter lie on the outer layers of the data cloud.
When the dimension $d$ is large, although the observations generated with the
smallest scale parameter almost always have their nearest neighbors coming
from the same distribution, the points in the outer layers may easily find
themselves ‘closer’ to points in the inner layer (generated from the
distribution with the smaller scales). This can lead to low power for the
‘KMD’ method with the discrete kernel. However, if it is a priori known that
$P_{1}$ may lie in the inner layer, then we can assign $K(1,1)$ a larger
value, which in turn would imply that $\hat{\eta}$ would be large as it would
be dominated by the contributions from the observations from $P_{1}$ (which
have their nearest neighbors coming from the same distribution). Figure 2
shows that when using the kernel $K_{2}$, defined as $K_{2}(1,1)=10$,
$K_{2}(2,2)=K_{2}(3,3)=1$, and $K_{2}(i,j)=0$, for $i\neq j$, the power of the
method is substantially increased.

Note that for the normal scale problem, both Roy’s and Pillai’s tests have low
power. But our method still works due to its generality, and as it is
consistent against any alternative. For the other settings considered here,
our KMD exhibits very good performance, achieving highest power among the
competing procedures.

######  Remark 4 (Choice of $k$ for the $k$-NN graph).

When testing equality of $M$ distributions, we chose $k=0.10n$ (in our
simulations), for samples with up to a few hundred observations (as advocated
in Petrie [57]). Note that when $\eta=0$, $\hat{\eta}$ is unbiased and a
larger $k$ reduces variance of $\hat{\eta}$. While, for estimating $\eta$
using its empirical version $\hat{\eta}$, a much smaller $k$ is recommended:
often $k=1$ would work the best. A larger $k$ often produces a smaller
$\hat{\eta}$ (e.g., $\hat{\eta}=0$ in the extreme case $k=n-1$), inducing a
large ‘bias’ especially when $\eta$ is large; see Appendix D.3 for empirical
evidence.

###  6.2 Measuring Multi-sample Dissimilarity with Real Data

In the following we analyse a few real data examples. Here we use a 1-NN graph
for constructing $\hat{\eta}$ and for testing. We first consider data sets
from the UCI Machine Learning Repository [26]. For each of these multi-sample
data sets, every observation/instance has a label indicating which sample it
came from. However, for these data sets, instances with different labels are
logically or physically different from each other, and have already been shown
to be different in Petrie [57].

To illustrate the usefulness of $\eta$ (and $\hat{\eta}$) as a measure of
dissimilarity between distributions, we construct various $M$ sample scenarios
by adapting these real data sets suitably, as described below. For each data
set, we compare distributions in the three settings: (i) the distributions to
be compared have the same label (so the null hypothesis (1) may hold here);
(ii) two distributions, out of the 3 distributions considered, have the same
label; and (iii) all the three distributions to be compared have distinct
labels. It is natural to expect that $\hat{\eta}$ would increase as we move
from scenario (i) to (iii); see Table 1 for our results.

Amazon commerce reviews [47]: Each observation/instance here is a
10000-dimensional vector representing usage of digits, punctuation, words,
sentence length, word frequencies, etc., of reviews from one of 50 extremely
active customers. The data set has in total $50$ labels — corresponding to the
50 customers — with 30 instances per label. For each instance, most variables
equal zero while the remaining are integers ranging from one to a few dozen.
As explained above we consider comparing: (i) the first 15 reviews from
reviewer 4 with her last 15 reviews; (ii) the first 15 reviews from reviewer
4, the last 15 reviews from reviewer 4, and all 30 reviews from reviewer 8;
(iii) the first 15 reviews from reviewer 4, the first 15 reviews from reviewer
12, and all 30 reviews from reviewer 8.

Semeion Handwritten Digits [68]: In this example we have 10 labels —
corresponding to the 10 digits — and each instance/observation is a
256-dimensional vector representing a written digit, with each coordinate
being 0 or 1 depending on its underlying grayscale value. Consider comparing:
(i) the first 81 instances of digit 6 with the last 80 instances of the same
digit; (ii) the first 81 instances of digit 6, the last 80 instances of digit
6, and all instances of digit 8; and (iii) the first 81 instances of digit 6,
the first 80 instances of digit 7, and all instances of digit 8.

ISOLET [24]: Each instance is a 617-dimensional vector representing a spoken
letter. 150 speakers spoke each letter in the English alphabet twice, so we
have 52 observations from each speaker. Consider comparing: (i) the 7-th
letter, i.e., ‘g’, spoken by the first 75 speakers with ‘g’ spoken by the last
75 speakers; (ii) ‘g’ spoken by the first 75 speakers, ‘g’ spoken by the last
75 speakers, and ‘t’ spoken by all speakers; and (iii) ‘g’ spoken by the first
75 speakers, ‘n’ spoken by the first 75 speakers, and ‘t’ spoken by all
speakers.

LRS [26]: Each instance in this data set is a 93-dimensional vector which
describes the fluxes from astronomical objects. There are in total 10 labels.
Consider comparing: (i) the first 96 instances with label 2 and the last 177
instances with label 2; (ii) the first 96 instances with label 2, the last 177
instances with label 2, and all the instances with label 1 (NA removed); and
(iii) all 96 instances with label 4, the last 177 instances with label 2, and
all the instances with label 1 (NA removed).

Table 1 gives the value of $\hat{\eta}$ for the 3 different settings for each
data set discussed above (for simplicity we use the standard Euclidean
distance to construct the $k$-NN graphs in these examples). It can be seen
from the table that the test for (1) based on $\hat{\eta}$ (and implemented
via a permutation test using 500 random permutations) is rejected (at level
0.05) for both settings (ii) and (iii) for all data sets. Moreover, our
measure $\hat{\eta}$ suggests that, for each of the data sets, the
distributions in setting (iii) are “more different” than those in setting
(ii), in the sense that $\hat{\eta}$ is closer to 1 for setting (iii). Note
that for the data sets ISOLET and LRS, the permutation test based on
$\hat{\eta}$ also reject setting (i) where the samples have the same label.
For the ISOLET data, this suggests that the way the first 75 speakers spoke
the 7-th letter is different from the way the last 75 speakers spoke the same
letter. A similar conclusion can be drawn for the LRS data set. If we randomly
select half of the instances with label 2 in LRS data set to form the first
sample, instead of the first half, then the $p$-value will typically not be
significant. Note that in the ISOLET data set, although all the hypothesis
tests yield the same $p$-value of $1/501$, $\hat{\eta}$ gives a more
meaningful summary of the dissimilarities between the distributions and
produces values that progressively increase from settings (i) to (iii).

Table 1: $\hat{\eta}$ and $p$-values (over 500 random permutations) for real data sets. Settings | (i) Same label | (ii) Mixture | (iii) Different labels
---|---|---|---
(Illustration) | $1$ $2$ | $1$ $2$ $3$ |  $1$ $2$ $3$
Data set | $\hat{\eta}$ |  $p$-value | $\hat{\eta}$ |  $p$-value | $\hat{\eta}$ |  $p$-value
Amazon | 0.033 | 0.513 | 0.423 | 0.004 | 0.607 | 0.002
Semeion | 0.099 | 0.160 | 0.604 | 0.002 | 0.975 | 0.002
ISOLET | 0.382 | 0.002 | 0.638 | 0.002 | 0.886 | 0.002
LRS | 0.168 | 0.014 | 0.508 | 0.002 | 0.868 | 0.002

Next we consider the two non-Euclidean examples introduced in Section 1.

Speech recognition: The ArabicDigits data set consists of a total of 8800
instances with 10 labels (numbers 0-9) from native Arabic speakers. We use the
popular DTW distance (as mentioned in the Introduction) between two multi-
dimensional time series [9, 25, 31] to compute the 1-NN graph. The KMD between
the 10 spoken digits has an incredibly high value of 0.9976, suggesting the 10
distributions have highly disjoint supports. This may be the reason why a
simple nearest neighbor classifier could achieve 99.8% prediction accuracy
[31]. We also investigate, using $\hat{\eta}$, the validity of the following
intuitive supposition: how males speak the number ‘1’ should be different from
how they speak ‘2’, but the difference should be smaller than that between how
males speak ‘1’ and how females speak ‘2’. To make the problem harder, we only
work with the 13th MFCC (instead of the 13-dimensional time-series), and the
corresponding two KMD estimates are 0.5687 and 0.6867 respectively, agreeing
with our intuition.

Sentiment analysis: The movie review data set from Selivanov et al. [65]
consists of 5000 movie reviews, each with a binary sentiment label
corresponding to whether the review is positive or negative. We process the
data by transferring texts to lower case, removing non-alphanumeric symbols,
collapsing multiple spaces, and using vocabulary-based vectorization with the
vocabulary pruned to have minimum number of occurrences over all documents
being 5, and maximum proportion of documents containing a term in the
vocabulary being 10%. This results in a 5000 $\times$ 12644 sparse document-
term matrix [65], and Jaccard distance [43] is used as the metric between
documents. The data processing is implemented using the R package text2vec
[65]. The estimated KMD between positive reviews and negative reviews is
0.445, suggesting some overlap in the supports of the two distributions.
Though a simple linear support vector machine classifier based on the
document-term matrix could provide 82% baseline accuracy (evaluated from
10-fold cross-validation), many later refined approaches, some of which were
specifically designed for sentiment analysis on such data sets, cannot improve
this accuracy to 90% [48]. This agrees with the observed value of KMD which is
much lower than in the previous example.

##  Appendix A Appendix

This Appendix will be organized as follows: In Appendix A, we provide some
general discussions that were deferred from the main paper. An analysis of the
asymptotic behavior of $\hat{\eta}$ under alternatives is given in Appendix B.
We show that both under a fixed alternative and under shrinking alternatives
converging to the null, $\hat{\eta}$ has an asymptotic normal distribution
after proper centering (see Theorems B.1 and B.2). Using the CLT under
shrinking alternatives, we provide a complete characterization of the local
power of our method for $\mathcal{Z}=\mathbb{R}^{d}$ (see Theorem B.3). In
particular, this provides the detection threshold of the test based on
$\hat{\eta}$. All proofs are given in Appendix C. Further simulation results
that were mentioned in the main paper are given in Appendix D.

###  A.1 Assumptions on the Geometric Graph

Let $\tilde{Z}_{1},\tilde{Z}_{2}\ldots$ be i.i.d. from the mixture
$\sum_{i=1}^{M}\pi_{i}P_{i}$ on $\mathcal{Z}$. Let
$\tilde{\mathcal{G}}_{n}:=\mathcal{G}(\tilde{Z}_{1},\ldots,\tilde{Z}_{n})$ be
the geometric graph with vertex set
$\\{\tilde{Z}_{1},\ldots,\tilde{Z}_{n}\\}\subset\mathcal{Z}$. Let
$\mathcal{E}(\mathcal{\tilde{G}}_{n})$ denote the set of (directed/undirected)
edges of $\tilde{\mathcal{G}}_{n}$, i.e.,
$(\tilde{Z}_{i},\tilde{Z}_{j})\in\mathcal{E}(\mathcal{\tilde{G}}_{n})$ if and
only if there is an edge from $\tilde{Z}_{i}$ to $\tilde{Z}_{j}$ in
$\tilde{\mathcal{G}}_{n}$, and $\tilde{d}_{i}$ denotes the out-degree of
$\tilde{Z}_{i}$ in $\tilde{\mathcal{G}}_{n}$. To be specific,
$\tilde{d}_{i}:=\sum_{j:(\tilde{Z}_{i},\tilde{Z}_{j})\in\mathcal{E}(\mathcal{\tilde{G}}_{n})}1$.
We assume the following conditions on $\tilde{\mathcal{G}}_{n}$ (as in Deb et
al. [21]):

######  Assumption 1.

Given the graph $\tilde{\mathcal{G}}_{n}$, let $N(1),\ldots,N(n)$ be
independent random variables where $N(i)$ is a uniformly sampled index from
among the (out-)neighbors of $\tilde{Z}_{i}$ in $\tilde{\mathcal{G}}_{n}$.
Then

$\rho_{\mathcal{Z}}(\tilde{Z}_{1},\tilde{Z}_{N(1)})\overset{p}{\to}0\quad{\rm
as\ }n\to\infty.$

######  Assumption 2.

Assume that there exists a deterministic positive sequence $r_{n}\geq 1$ (may
or may not be bounded), such that almost surely:

$\min_{1\leq i\leq n}\tilde{d}_{i}\geq r_{n}.$

Let $\tilde{\mathcal{G}}_{n,i}$ denote the graph obtained from
$\tilde{\mathcal{G}}_{n}$ by replacing $\tilde{Z}_{i}$ with an i.i.d. random
element $\tilde{Z}_{i}^{\prime}$. Assume that there exists a deterministic
positive sequence $q_{n}$ (may or may not be bounded), such that

$\max_{1\leq i\leq
n}\max\\{|\mathcal{E}(\tilde{\mathcal{G}}_{n})\setminus\mathcal{E}(\tilde{\mathcal{G}}_{n,i})|,|\mathcal{E}(\tilde{\mathcal{G}}_{n,i})\setminus\mathcal{E}(\tilde{\mathcal{G}}_{n})|\\}\leq
q_{n}\ \ a.s.\qquad{\rm and}\qquad\frac{q_{n}}{r_{n}}=O(1).$

######  Assumption 3.

There exists a deterministic sequence $\\{t_{n}\\}_{n\geq 1}$ (may or may not
be bounded) such that the vertex degree (including both in- and out-degrees
for directed graphs, i.e., $\\#\\{j:(\tilde{Z}_{i},\tilde{Z}_{j})\ {\rm or}\
(\tilde{Z}_{j},\tilde{Z}_{i})\in\mathcal{E}(\tilde{\mathcal{G}}_{n})\\}$) of
every point $\tilde{Z}_{i}$ (for $i=1,\ldots,n$) is bounded by $t_{n}$, and
$\frac{t_{n}}{r_{n}}=O(1)$.

Assumption 1 formalizes our intuition that the presence of an edge between two
points implies that the two points are close. Assumption 2 states that the
graph is ‘local’ in the sense that replacing one random point will not change
too many edges. Assumption 3 requires that the degree of each vertex should be
of the same order. See Deb et al. [21, Section 3] for a detailed discussion on
these assumptions.

###  A.2 Computational Complexity and Implementation Schemes

When a Euclidean $k$-NN graph is used, the computation complexity of
$\hat{\eta}$ and $\frac{\hat{\eta}}{\sqrt{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}}$ is $O(kn\log n)$, which is near linear,
when $k$ is bounded. This is due to the fact that the Euclidean $k$-NN graph
can be computed in $O(kn\log n)$ time (for example, using the k-d tree; see
[6]).

The computation of the $k$-NN graph is well implemented in many computational
softwares [51, 3]. The computation of $\tilde{a},\tilde{b},\tilde{c}$ (see
(6)) and $\tilde{g}_{2}$ is described below. Observe that

$\displaystyle\tilde{a}$
$\displaystyle=\frac{1}{n(n-1)}\left(\sum_{i=1}^{M}n_{i}(n_{i}-1)K^{2}(i,i)+\sum_{i\neq
j}n_{i}n_{j}K^{2}(i,j)\right)$
$\displaystyle=\frac{1}{n(n-1)}\left(\sum_{i,j=1}^{M}n_{i}n_{j}K^{2}(i,j)-\sum_{i=1}^{M}n_{i}K^{2}(i,i)\right)$

which can be computed in $O(M^{2})=O(1)$ time since $M$ is a constant.
Further, note that,

$\displaystyle\tilde{b}$
$\displaystyle=\frac{1}{n(n-1)(n-2)}\Big{(}\sum_{i=1}^{M}K(i,i)K(i,i)n_{i}(n_{i}-1)(n_{i}-2)$
$\displaystyle\qquad+{\sum}^{\prime}K(i,i)K(i,j)n_{i}(n_{i}-1)n_{j}+{\sum}^{\prime}K(i,j)K(i,i)n_{i}n_{j}(n_{i}-1)$
$\displaystyle\qquad+{\sum}^{\prime}K(i,j)K(i,j)n_{i}n_{j}(n_{j}-1)+{\sum}^{\prime}K(i,j)K(i,l)n_{i}n_{j}n_{l}\Big{)}$
$\displaystyle=\frac{1}{n(n-1)(n-2)}\Big{(}\sum_{i,j,l=1}^{M}K(i,j)K(i,l)n_{i}n_{j}n_{l}-\sum_{i=1}^{M}K(i,i)^{2}n_{i}(3n_{i}-2)$
$\displaystyle\qquad-2\
{\sum}^{\prime}K(i,i)K(i,j)n_{i}n_{j}-{\sum}^{\prime}K(i,j)^{2}n_{i}n_{j}\Big{)}$
$\displaystyle=\frac{1}{n(n-1)(n-2)}\Big{(}\sum_{i,j,l=1}^{M}K(i,j)K(i,l)n_{i}n_{j}n_{l}+2\sum_{i=1}^{M}K(i,i)^{2}n_{i}$
$\displaystyle\qquad-2\sum_{i,j=1}^{M}K(i,i)K(i,j)n_{i}n_{j}-\sum_{i,j=1}^{M}K(i,j)^{2}n_{i}n_{j}\Big{)}$
$\displaystyle=\frac{1}{n(n-1)(n-2)}\Big{(}\sum_{i,j,l=1}^{M}K(i,j)K(i,l)n_{i}n_{j}n_{l}+\sum_{i=1}^{M}K(i,i)^{2}n_{i}$
$\displaystyle\qquad-2\sum_{i,j=1}^{M}K(i,i)K(i,j)n_{i}n_{j}-n(n-1)\tilde{a}\Big{)}$
$\displaystyle=\frac{1}{n(n-1)(n-2)}\Big{(}\vec{\textbf{n}}^{\top}[K(i,j)]{\rm
diag}(\vec{\textbf{n}})[K(i,j)]\vec{\textbf{n}}+\sum_{i=1}^{M}K(i,i)^{2}n_{i}$
$\displaystyle\qquad-2\sum_{i,j=1}^{M}K(i,i)K(i,j)n_{i}n_{j}-n(n-1)\tilde{a}\Big{)},$

where $\vec{\textbf{n}}=(n_{1},\ldots,n_{M})^{\top}$, ${\rm
diag}(\vec{\textbf{n}})$ is the diagonal matrix with $\vec{\textbf{n}}$ being
in the diagonal, and $[K(i,j)]$ is the $M\times M$ matrix with the entry in
the $i$-th row and $j$-th column being $K(i,j)$. Hence, $\tilde{b}$ can be
computed in $O(M^{2})=O(1)$ time as well. Similarly,

$\displaystyle\tilde{c}$
$\displaystyle=\frac{\big{(}\sum_{i,j=1}^{M}n_{i}n_{j}K(i,j)-\sum_{i=1}^{M}n_{i}K(i,i)\big{)}^{2}-4n(n-1)(n-2)\tilde{b}-2n(n-1)\tilde{a}}{n(n-1)(n-2)(n-3)}$

can be computed in $O(1)$ time. Recall that,
$\tilde{g}_{2}=\frac{1}{n}\sum_{i,j=1}^{n}\frac{T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}+\frac{1}{n}\sum_{i\neq
j}\frac{T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}}$. To compute $\sum_{i\neq
j}\frac{T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}}$, for each $Z_{i}$, we find its
in-neighbors, i.e.,
$\\{Z_{j}:(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})\\}=\\{Z_{i_{1}},\ldots,Z_{i_{k^{\prime}}}\\}$.
Let
$S_{i}:=\left(\frac{1}{d_{i_{1}}}+\ldots+\frac{1}{d_{i_{k^{\prime}}}}\right)^{2}-\sum_{j=1}^{k^{\prime}}\frac{1}{d_{i_{j}}^{2}}$.
Then $\sum_{i\neq
j}\frac{T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}}=\sum_{i=1}^{n}S_{i}$. The
computational complexity for $\tilde{g}_{2}$ is $O(kn)$ since $k^{\prime}\leq
kC(d)$ for some constant $C(d)$ depending on the dimension $d$ of
$\mathcal{Z}$ [72, Lemma 8.4]. The computations of other terms are straight
forward.

###  A.3 Interpretation of $\hat{\eta}$ as $k$-NN Cross-Validated Accuracy

Here we show that $\hat{\eta}$ is linearly related to the leave-one-out cross-
validation accuracy of a $k$-NN classifier.

A classifier uses existing data $\\{(\Delta_{i},Z_{i})\\}_{i\in I}$, where
$I\subset\\{1,\ldots,n\\}$ is an index set, to make prediction for the label
$\Delta_{j}$ of a new observation $Z_{j}$, $j\notin I$. The leave-one-out
cross-validated accuracy (abbreviated as “accuracy” in the following) of a
possibly random classifier $f:\mathcal{Z}\to\\{1,\ldots,M\\}$ is defined as:

$f\ {\rm
Accuracy}:=\frac{1}{n}\sum_{j=1}^{n}\mathbb{P}_{f}\Big{(}\Delta_{j}=f\left(Z_{j}|\\{(\Delta_{i},Z_{i})\\}_{i\in\\{1,\ldots,n\\}\backslash\\{j\\}}\right)\Big{)},$

where
$f\left(Z_{j}|\\{(\Delta_{i},Z_{i})\\}_{i\in\\{1,\ldots,n\\}\backslash\\{j\\}}\right)$
is the prediction for $\Delta_{j}$ using $Z_{j}$ given by $f$ learnt on data
$\\{(\Delta_{i},Z_{i})\\}_{i\in\\{1,\ldots,n\\}\backslash\\{j\\}}$, and
$\mathbb{P}_{f}$ means averaging over possible randomness in $f$.

Consider the accuracy of: (i) $k$-NN classifier — the label of a data point is
predicted by the label of a random $k$-NN of that data point; (ii) random
guess classifier — the label of a data point is predicted by a random guess
according to the proportion of different labels in the existing data.

When $K(x,y)=I(x=y)$ and a $k$-NN graph is used to define $\hat{\eta}$ via
(4), we have

$\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{k}\sum_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}I(\Delta_{i}=\Delta_{j})=k\text{-}{\rm
NN\ Accuracy}$

and

$\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K(\Delta_{i},\Delta_{j})=\frac{1}{n}\sum_{i=1}^{n}\frac{\sum_{j:j\neq
i}I(\Delta_{i}=\Delta_{j})}{n-1}={\rm Random\ Guess\ Accuracy}.$

Hence $\eta$ in (4) reduces to

$\hat{\eta}=\frac{(k\text{-}{\rm NN\ Accuracy})-({\rm Random\ Guess\
Accuracy})}{1-({\rm Random\ Guess\ Accuracy})}.$ (A.1)

Intuitively, if the $M$ distributions are concentrated on different regions of
the space $\mathcal{Z}$, then the $k$-NN accuracy is close to 1 and thus
$\hat{\eta}$ will be close to 1. On the other hand, if the $M$ distributions
are the same, then the $k$-NN accuracy will be similar to the random guess
accuracy and thus $\hat{\eta}$ will be close to 0. Moreover, $\hat{\eta}$
grows linearly with $k$-NN accuracy, which provides a simple and intuitive
interpretation of our measure $\hat{\eta}$.

##  Appendix B Asymptotic Behavior of $\hat{\eta}$ under Alternatives

In this section, we will first establish CLTs for $\hat{\eta}$ under a fixed
alternative and under certain shrinking alternatives converging to the null.
Using these CLTs, we characterize the asymptotic power of the test based on
$\hat{\eta}$. In particular, this provides the detection threshold of the test
based on $\hat{\eta}$. We will focus on $\mathcal{Z}=\mathbb{R}^{d}$ ($d\geq
1$) in this section, and assume that $P_{i}$ has a density $f_{i}$ w.r.t. the
Lebesgue measure, for all $i=1,\ldots,M$.

The asymptotic distribution of geometric graph-based statistics under
alternatives and their detection thresholds were not available until the
recent work of Bhattacharya [11], where the analysis was carried out under a
Poissonized setting. The spatial independence of the Poisson process
facilitates the computation of the asymptotic variance of $\hat{\eta}$, and
also helps establish the CLTs for $\hat{\eta}$ using Stein’s method [18]. We
will also consider this Poissonized framework in this section. Although de-
Poissonization techniques are available in the literature [52, 54], the de-
Poissonized version of our theorem (Theorem B.1) unfortunately does not reduce
to our original setting in the main paper where the sample size $n_{i}$’s are
nonrandom.

In the Poissonized framework, instead of assuming $\frac{n_{i}}{n}\to\pi_{i}$,
it is assumed that $n_{i}\sim{\rm Poisson}(N_{i})$, for $i=1,\ldots,M$. Write
$N:=\sum_{i=1}^{M}N_{i}$ and assume $\frac{N_{i}}{N}\to\pi_{i}>0$ as
$N\to\infty$. An equivalent characterization is that we first decide to draw
in total $n\sim{\rm Poisson}(N)$ data points, and for each data point, with
probability $\frac{N_{i}}{N}$, we draw an observation from the $i$-th
population, for $i=1,\ldots,M$. In the following, we provide CLTs for
$\hat{\eta}$ under fixed and shrinking alternatives, under the Poissonized
setting when we let $N\to\infty$.

###  B.1 Asymptotic Normality under a General Fixed Alternative

In the following, we briefly sketch how a CLT for $\hat{\eta}$, under a fixed
alternative, can be obtained after proper centering; see Theorem B.1 for the
precise statement. Since the Poisson distribution is tightly concentrated
around its mean: ${\rm Var}\left(\frac{n_{i}}{N_{i}}\right)=\frac{1}{N_{i}}\to
0$ as $N_{i}\to\infty$. Thus, we have
$\frac{n_{i}}{N_{i}}\overset{p}{\to}\mathbb{E}\left[\frac{n_{i}}{N_{i}}\right]=1$.
Hence the denominator of $\hat{\eta}$ in (4) converges in probability to
$\sum_{i=1}^{M}\pi_{i}K(i,i)-\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)>0$. So our
analysis focuses on the asymptotic behavior of the numerator in (4):

$\tilde{H}_{n}:=\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(i,j)\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})-\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K(\Delta_{i},\Delta_{j}).$

When the above quantity cannot be well-defined, i.e., $n=0$ or $1$, or $n\leq
k$ when a $k$-NN graph is used, we set $\tilde{H}_{n}=\hat{\eta}=0$. Using
$U$-statistics projection theory [21, Lemma D.4], it can be shown that
$n\tilde{H}_{n}/\sqrt{N}$ has the same asymptotic distribution as $H_{n}$
defined as

$H_{n}:=\frac{1}{\sqrt{N}}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})-\frac{1}{\sqrt{N}}\sum_{i=1}^{n}g_{N}(\Delta_{i}),$

which has a non-degenerate asymptotic normal distribution after centering by
its mean. Here,

$g_{N}(\Delta_{i}):=2\sum_{p=1}^{M}\frac{N_{p}}{N}K(\Delta_{i},p)-\sum_{p,q=1}^{M}\frac{N_{p}N_{q}}{N^{2}}K(p,q).$

Similar to the CLT result under ${\rm H}_{0}$ (see Theorem 3), if we let
$\mathcal{F}_{n}:=\sigma(n,Z_{1},\ldots,Z_{n})$, the $\sigma$-algebra
generated by the unlabelled data and the number of total observations, then
$H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})$ given $\mathcal{F}_{n}$ converges in
distribution to a normal limit $N(0,\kappa_{1}^{2})$ (see Theorem B.1 for the
exact expression for $\kappa_{1}^{2}$). We can also show that the
unconditional distribution of
$\mathbb{E}(H_{n}|\mathcal{F}_{n})-\mathbb{E}(H_{n})$ converges to another
normal limit $N(0,\kappa_{2}^{2})$; see Theorem B.1 for the exact expression
for $\kappa_{2}^{2}$. Finally, a simple argument using characteristic
functions shows that
$H_{n}-\mathbb{E}(H_{n})\overset{d}{\to}N(0,\kappa_{1}^{2}+\kappa_{2}^{2})$.

If ${\rm H}_{0}$ holds, then it can be shown that
$\mathbb{E}(H_{n}|\mathcal{F}_{n})=\mathbb{E}(H_{n})=0$, and $\kappa_{1}^{2}$
is exactly the distribution-free variance derived in Theorem 4, and
$\kappa_{2}^{2}=0$. However, in a general situation,
$\mathbb{E}(H_{n}|\mathcal{F}_{n})$ is a function of
$\\{Z_{1},\ldots,Z_{n}\\}$, which follows a non-homogeneous Poisson
process888For a non-homogeneous Poisson process with intensity function $f$
(denoted by $\mathcal{P}_{f}$), the number of points in a set
$A\subset\mathbb{R}^{d}$ follows a ${\rm Poisson}(\int_{A}f(x){\rm d}x)$
distribution, and the number of points in disjoint sets are independent.
$\mathcal{P}_{N\phi_{N}}$ under the Poissonized setting, with

$\phi_{N}(z):=\sum_{i=1}^{M}\frac{N_{i}}{N}f_{i}(z)\qquad{\rm
and}\qquad\phi(z):=\sum_{i=1}^{M}\pi_{i}f_{i}(z),\quad\mbox{for
}\;z\in\mathcal{Z},$ (B.1)

being the marginal density of $Z$ given the total number of observations, and
its limiting value, respectively. To show the convergence of the variance of
$\mathbb{E}(H_{n}|\mathcal{F}_{n})$ and establish a CLT for
$\mathbb{E}(H_{n}|\mathcal{F}_{n})-\mathbb{E}(H_{n})$ using Stein’s method, we
will need to make another common assumption — the power-law stabilization [55,
54, 11] — on the non-homogeneous Poisson process $\mathcal{P}_{N\phi_{N}}$ (in
addition to the stabilization on the homogeneous Poisson process
$\mathcal{P}_{\lambda}$ defined in Section 4.2), which will be introduced
below.

Recall our setup: $\mathcal{G}$ is a geometric graph on $\mathcal{Z}$ and
$\mathcal{X}$ is a set of points in $\mathcal{Z}$.
$\mathcal{E}\left(x,\mathcal{G}(\mathcal{X})\right)$ is the set of edges in
$\mathcal{G}(\mathcal{X})$ that are incident to $x\in\mathcal{X}$. Let $A$ be
a set with $\phi$-probability 1, e.g., $A={\rm supp}(\phi)$, the support of
$\phi$. Fix $x\in\mathcal{Z}$. A radius of stabilization $R\equiv R(N,x)$ at
$x$ exists if the points in $A\backslash B(x,N^{-1/d}R)$ cannot impact the
neighborhood of $x$ in $\mathcal{G}\left(\mathcal{P}_{N\phi_{N}}^{x}\right)$
where $\mathcal{P}_{N\phi_{N}}^{x}:=\mathcal{P}_{N\phi_{N}}\cup\\{x\\}$, i.e.,

$\mathcal{E}\left(x,\mathcal{G}\left([\mathcal{P}_{N\phi_{N}}^{x}\cap
B(x,N^{-1/d}R)]\cup\mathcal{X}\right)\right)=\mathcal{E}\left(x,\mathcal{G}\left(\mathcal{P}_{N\phi_{N}}^{x}\cap
B(x,N^{-1/d}R)\right)\right),$ (B.2)

for all finite $\mathcal{X}\subset A\backslash B(x,N^{-1/d}R)$. Such a finite
$R$ always exists when $A$ is bounded. For certain graphs such as a $k$-NN
graph, the power-law stabilization states that $R$ is bounded in probability
under mild assumptions, and the tail probability decays sufficiently fast,
uniformly over $N$ and $x$. To define it formally, let

$\tau(t):=\sup_{N\geq 1,x\in A}\mathbb{P}[R(N,x)>t],\qquad\mbox{for}\;\;t>0.$
(B.3)

Then, $\mathcal{G}$ is said to be:

  1. 1.

power-law stabilizing of order $q$ with respect to $\phi_{N}$ if $\sup_{t\geq
1}t^{q}\tau(t)<\infty$,

  2. 2.

exponentially stabilizing with respect to $\phi_{N}$ if
$\limsup_{t\to\infty}t^{-1}\log\tau(t)<0$.

Note that if $\mathcal{G}$ is exponentially stabilizing, then it is power-law
stabilizing for all $q>0$. It is known that a $k$-NN graph is exponentially
stabilizing under suitable conditions (see Proposition B.1 below and Appendix
C.11 for its proof). Other exponentially stabilizing graphs include the
Voronoi and Delaunay graphs [54, 11]. However, there is no similar tail bounds
known for MST, as far as we are aware.

######  Proposition B.1 (Penrose [54, Section 6.3]).

Assume $\phi_{N}$ is bounded away from 0 on its support which is assumed to be
convex. Then the $k$-NN graph (with $k$ fixed, either directed or undirected)
is exponentially stabilizing.

Assuming that the graph $\mathcal{G}$ is power-law stabilizing of sufficient
order w.r.t. our data generating process $\phi_{N}$, the asymptotic normality
of $\hat{\eta}$ can be established, as shown below.

######  Theorem B.1 (CLT under a fixed alternative).

Suppose that we have the Poissonized setting where we have $n_{i}\sim{\rm
Poisson}(N_{i})$ samples from the $i$-th distribution, for $i=1,\ldots,M$.
Suppose the following assumptions hold:

  1. 1.

$\sqrt{N}\left(\frac{N_{i}}{N}-\pi_{i}\right)\to 0$ as $N\to\infty$, for
$i=1,\ldots,M$, where $N=\sum_{i=1}^{M}N_{i}$.

  2. 2.

The $i$-th distribution has a Lebesgue density $f_{i}$ on $\mathbb{R}^{d}$
which is Lebesgue almost everywhere continuous, for $i=1,\ldots,M$.

  3. 3.

$\phi:=\sum_{i=1}^{M}\pi_{i}f_{i}$ is bounded above and has a bounded support.

  4. 4.

$\mathcal{G}$ is translation and scale invariant, stabilizing on
$\mathcal{P}_{\lambda}$ for some $\lambda>0$.

  5. 5.

$\mathcal{G}$ is power-law stabilizing w.r.t. $\phi_{N}$ (see (B.1)) with
order $q_{0}>16d$, with the corresponding full probability set $A$ and radius
of stabilization $R$ as defined in (B.2).

  6. 6.

The degree bound (including in- and out-degrees) for
$\mathcal{G}(\mathcal{P}_{N\phi_{N}})$ is $o_{p}(N^{1/40})$.

Then with $n:=\sum_{i=1}^{M}n_{i}$,

$\frac{1}{\sqrt{n}}\left(n\tilde{H}_{n}-\mathbb{E}[n\tilde{H}_{n}]\right)\overset{d}{\to}N(0,\kappa_{1}^{2}+\kappa_{2}^{2}),$
(B.4)

where

$\displaystyle\kappa_{1}^{2}$ $\displaystyle:=(g_{1}+g_{3})\int{\rm
Var}[K(\tilde{\Delta},\tilde{\Delta}^{\prime})|\tilde{Z}=z]\phi(z){\rm d}z$
$\displaystyle\qquad+(3-2g_{1}-2g_{3}+g_{2})\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),K(\tilde{\Delta},\tilde{\Delta}^{\prime\prime})|\tilde{Z}=z)\phi(z){\rm
d}z$ $\displaystyle\qquad-4\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),g(\tilde{\Delta})|\tilde{Z}=z)\phi(z){\rm
d}z+\mathbb{E}\left[{\rm Var}(g(\tilde{\Delta})|\tilde{Z})\right],$
$\displaystyle\kappa_{2}^{2}$
$\displaystyle:=\int\left(\sum_{i,j=1}^{M}K(i,j)\frac{\pi_{i}f_{i}(z)}{\phi(z)}\frac{\pi_{j}f_{j}(z)}{\phi(z)}-\sum_{i=1}^{M}g(i)\frac{\pi_{i}f_{i}(z)}{\phi(z)}\right)^{2}\phi(z){\rm
d}z$
$\displaystyle\;=\int\left(\mathbb{E}[K(\tilde{\Delta},\tilde{\Delta}^{\prime})|\tilde{Z}=z]-\mathbb{E}[g(\tilde{\Delta})|\tilde{Z}=z]\right)^{2}\phi(z){\rm
d}z.$

Note that in the above displays,
$\tilde{\Delta},\tilde{\Delta}^{\prime},\tilde{\Delta}^{\prime\prime}$ are
drawn independently from the conditional distribution
$\tilde{\Delta}|\tilde{Z}=z$ (recall the definition of
$\tilde{\Delta},\tilde{Z}$ in Section 2.1);
$g(\Delta):=2\sum_{i=1}^{M}\pi_{i}K(\Delta,i)-\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)$,
and $g_{1},g_{2},g_{3}$ are as defined in Theorem 4.

Observe that $\kappa_{1}$ and $\kappa_{2}$ are not distribution-free in
general, unless $\tilde{\Delta}$ is independent of $\tilde{Z}$, i.e., ${\rm
H}_{0}$ holds. The proof of the above result is given in Appendix C.12,
applying similar techniques as in Bhattacharya [11]. Note that Bhattacharya
[11] has applied the results in Penrose [54], but some of the assumptions in
Penrose [54] such as almost everywhere continuity and boundedness of densities
are not explicitly stated in Bhattacharya [11]. Here, due to the presence of a
general kernel and the U-statistic term $\sum_{i\neq
j}K(\Delta_{i},\Delta_{j})$ we cannot directly apply the results in Penrose
[54], and hence we provide a complete proof under the assumptions stated
above.

######  Remark B.1 (Non-degeneracy with characteristic kernel).

Suppose $K$ is characteristic. When $\tilde{\Delta}$ is independent of
$\tilde{Z}$, then $\kappa_{2}^{2}=0$, and $\kappa_{1}^{2}$ reduces to the
distribution-free null variance in Theorem 4, which is strictly positive. When
$\tilde{\Delta}$ is not independent of $\tilde{Z}$, $\kappa_{2}^{2}$ is
strictly positive since the integrand is the fourth power of the maximum mean
discrepancy999 MMD is a distance (provided that the kernel is characteristic)
on the space of probability measures defined through the kernel
$K(\cdot,\cdot)$ [33]. For two probabilities $P$ and $Q$, let independent
samples $\Delta_{1}^{P},\Delta_{2}^{P}$ i.i.d. follow $P$, and
$\Delta_{1}^{Q},\Delta_{2}^{Q}$ i.i.d. follow $Q$. A definition of MMD between
$P,Q$ is given by: ${\rm
MMD}^{2}(P,Q):=\mathbb{E}K(\Delta_{1}^{P},\Delta_{2}^{P})+\mathbb{E}K(\Delta_{1}^{Q},\Delta_{2}^{Q})-2\mathbb{E}K(\Delta_{1}^{P},\Delta_{1}^{Q}).$
(MMD) between the conditional distribution $\tilde{\Delta}|\tilde{Z}=z$ and
the unconditional distribution of $\tilde{\Delta}$ (note that with a
characteristic kernel, the MMD is positive whenever the two distributions are
not equal). Hence the asymptotic variance in (B.4) is always positive.

###  B.2 A CLT under Shrinking Alternatives

In the previous subsection, a CLT was derived under a fixed alternative. In
order to study the detection threshold of our statistic $\hat{\eta}$, a CLT
under shrinking alternatives is needed, which will be studied in this
subsection. Note that Bhattacharya [11] only provided a CLT under a fixed
alternative (cf., Theorem B.1); the CLT under shrinking alternatives was not
explicitly formulated, but was assumed to hold instead. As will be stated in
the following theorem, it turns out that a uniform version of the conditions
for the CLT under fixed alternatives (Theorem B.1) can guarantee a CLT under
shrinking alternatives.

We assume for each $N$, the $M$ distributions have densities
$f_{1}^{N},\ldots,f_{M}^{N}$ on $\mathcal{Z}\subset\mathbb{R}^{d}$
respectively. We write a $N$ at the top right of each $f_{i}$ to emphasize
that the densities change with the sample size $N$.

Assume $f_{i}^{N}(z)\overset{}{\longrightarrow}f(z)$ as $N\to\infty$ for a.e.
$z\in\mathcal{Z}$ (under the Lebesgue measure), for $i=1,\ldots,M$. Let

$\phi^{N}_{N}(z):=\sum_{i=1}^{M}\frac{N_{i}}{N}f_{i}^{N}(z)\qquad{\rm
and}\qquad\phi^{N}(z):=\sum_{i=1}^{M}\pi_{i}f_{i}^{N}(z),\quad{\rm
for}\;\;z\in\mathcal{Z}.$

Then $\phi^{N}_{N}(z)$ is the marginal density of $Z$ given the total number
of observations. To extend the CLT result to shrinking alternatives, we need
to also extend the notion of stabilization for the non-homogeneous Poisson
process $\mathcal{P}_{N\phi^{N}_{N}}$ below. We note that the distribution of
the set of points $\\{Z_{1},\ldots,Z_{n}\\}$ is exactly a non-homogeneous
Poisson process $\mathcal{P}_{N\phi^{N}_{N}}$ [53, Proposition 1.5].

Consider the notation introduced at the start of Appendix B and recall the
Poissonized setting. Let $A^{N}$ be a set with $\phi^{N}$-probability 1. Fix
$x\in\mathcal{Z}$. A radius of stabilization $R\equiv R(N,x)$ at $x$ exists if
the points in $A^{N}\backslash B(x,N^{-1/d}R)$ cannot impact the neighborhood
of $x$ in
$\mathcal{P}_{N\phi^{N}_{N}}^{x}:=\mathcal{P}_{N\phi^{N}_{N}}\cup\\{x\\}$,
i.e.,

$\mathcal{E}\left(x,\mathcal{G}\left([\mathcal{P}_{N\phi^{N}_{N}}^{x}\cap
B(x,N^{-1/d}R)]\cup\mathcal{X}\right)\right)=\mathcal{E}\left(x,\mathcal{G}\left(\mathcal{P}_{N\phi^{N}_{N}}\cap
B(x,N^{-1/d}R)\right)\right),$ (B.5)

for all finite $\mathcal{X}\subset A^{N}\backslash B(x,N^{-1/d}R)$. With

$\tau(t):=\sup_{N\geq 1,x\in A^{N}}\mathbb{P}[R(N,x)>t],\qquad\mbox{for
}\;\;t>0,$

$\mathcal{G}$ is said to be:

  1. 1.

power-law stabilizing of order $q$ with respect to $\phi^{N}_{N}$ if
$\sup_{t\geq 1}t^{q}\tau(t)<\infty$,

  2. 2.

exponentially stabilizing with respect to $\phi^{N}_{N}$ if
$\limsup_{t\to\infty}t^{-1}\log\tau(t)<0$.

Almost the same proof of Proposition B.1 shows that if $\phi^{N}_{N}$ has
convex support and $\phi^{N}_{N}\geq c$ on ${\rm supp}(\phi^{N}_{N})$ for some
$c>0$ and all $N$, then $k$-NN graph is exponentially stabilizing with respect
to $\phi^{N}_{N}$, and thus power-law stabilizing of any order $q>0$. The
following theorem (proved in Appendix C.13) formally states the CLT under
shrinking alternatives.

######  Theorem B.2 (CLT under shrinking alternatives).

Assume the Poissonized setting where we have $n_{i}\sim{\rm Poisson}(N_{i})$
samples from the $i$-th distribution, for $i=1,\ldots,M$. Suppose the
following assumptions hold:

  1. 1.

$\sqrt{N}\left(\frac{N_{i}}{N}-\pi_{i}\right)\to 0$, as $N\to\infty$,
$i=1,\ldots,M$, where $N=\sum_{i=1}^{M}N_{i}$.

  2. 2.

The $i$-th distribution has a Lebesgue density $f_{i}^{N}$ on $\mathbb{R}^{d}$
such that $f_{i}^{N}$ converges pointwise Lebesgue a.e. to a density $f$ as
$N\to\infty$, and $f_{i}^{N}$ is equicontinuous almost everywhere, i.e., for
Lebesgue a.e. $x$, for any $\varepsilon>0$, there exists $\delta>0$ such that
$|f_{i}^{N}(y)-f_{i}^{N}(x)|<\varepsilon$ whenever $|y-x|<\delta$, for
$i=1,\ldots,M$ and all $N\geq 1$.

  3. 3.

$f_{i}^{N}$ is uniformly bounded above and has uniformly bounded support,
i.e., there exists $C>0$ such that $\|f_{i}^{N}\|_{\infty}<C$ for all $N$ and
$i=1,\ldots,M$, and there exists a bounded set that contains ${\rm
supp}(f_{i}^{N})$ for all $N$ and $i=1,\ldots,M$.

  4. 4.

$\mathcal{G}$ is translation and scale invariant, stabilizing on
$\mathcal{P}_{\lambda}$ for some $\lambda>0$.

  5. 5.

$\mathcal{G}$ is power-law stabilizing with respect to $\phi_{N}^{N}$ with
order $q_{0}>d$, with the corresponding full probability sets $A^{N}$ and
radius of stabilization $R$ as defined in (B.5).

  6. 6.

The degree bound (including in- and out-degrees) for
$\mathcal{G}(\mathcal{P}_{N\phi_{N}^{N}})$ is $o_{p}(N^{1/40})$.

Then with $n=\sum_{i=1}^{M}n_{i}$,

$\frac{1}{\sqrt{n}}\left(n\tilde{H}_{n}-\mathbb{E}[n\tilde{H}_{n}]\right)\overset{d}{\to}N(0,\kappa_{1,{\rm
null}}^{2}),$

where

$\kappa_{1,{\rm
null}}^{2}:={a}\left({g}_{1}+{g}_{3}\right)+{b}\left({g}_{2}-2{g}_{1}-2{g}_{3}-1\right)+{c}\left({g}_{1}-{g}_{2}+{g}_{3}+1\right)>0$
(B.6)

is equal to the numerator of the asymptotic null variance given in (8), which
is distribution-free, not depending on $\\{f_{i}^{N}\\}_{i=1}^{M}$ and $f$.

######  Remark B.2 (Conditions 2 and 3 for parametric models).

Suppose that $f_{i}^{N}(\cdot)=f(\cdot|\theta_{i}^{N})$ is parametrized by
$\theta_{i}^{N}$ and the parametric family
$\\{f(\cdot|\theta):\theta\in\Theta\subset\mathbb{R}^{p}\\}$ ($p\geq 1$) with
a common compact support $A$. If $\theta_{i}^{N}$ depends on $N$ and converges
to $\theta_{0}$ in the interior of $\Theta$ as $N\to\infty$, $i=1,\ldots,M$,
and $f(x|\theta)$ is continuous in $(x,\theta)$, then conditions 2 and 3 above
hold, because $f(x|\theta)$ is uniformly continuous in $(x,\theta)\in A\times
B(\theta_{0},\varepsilon)$ (for some $\varepsilon>0$), a compact set.

###  B.3 Asymptotic Power and Detection Threshold

In this subsection, we consider $M=2$ and study the power of the test based on
$\hat{\eta}$ along a parametric sub-model [11]: suppose $P_{i}$ has Lebesgue
density $f(\cdot|\theta_{i})$ on $\mathbb{R}^{d}$ belonging to the parametric
family $\\{f(\cdot|\theta):\theta\in\Theta\subset\mathbb{R}^{p}\\}$ ($p\geq
1$), for $i=1,2$.

Our analysis depends on Theorem B.2, the CLT under shrinking alternatives.
Under the null, $\mathbb{E}[n\tilde{H}_{n}]=0$. Hence, the local power of the
test that rejects the null when

$\frac{\sqrt{n}\tilde{H}_{n}}{\kappa_{1,{\rm null}}}\geq
z_{1-\alpha}\qquad{\rm or\ equivalently\
(asymptotically})\qquad\frac{\sqrt{n}\hat{\eta}}{\sigma_{\mathcal{G},K,d,\pi}}\geq
z_{1-\alpha}$ (B.7)

is determined by the rate at which $\mathbb{E}[n\tilde{H}_{n}]/\sqrt{n}$
converges to 0, since
$\frac{n\tilde{H}_{n}-\mathbb{E}[n\tilde{H}_{n}]}{\sqrt{n}\kappa_{1,{\rm
null}}}\overset{d}{\to}N(0,1)$ (recall that $\sigma_{\mathcal{G},K,d,\pi}$ has
been defined in (8), and $\kappa_{1,{\rm null}}$ is defined in (B.6)). We next
state the detection threshold of the test (B.7), which is of the same form as
in Bhattacharya [11].

######  Theorem B.3 (Detection threshold).

Suppose $\\{\mathbb{P}_{\theta}\\}_{\theta\in\Theta}$ is a parametric family
of distributions with an open convex parameter space
$\Theta\subset\mathbb{R}^{p}$ (for $p\geq 1$) and Lebesgue densities
$\\{f(\cdot|\theta)\\}_{\theta\in\Theta}$. Assume:

  1. 1.

For all $\theta\in\Theta$, $f(\cdot|\theta)$ has a compact and convex support
$S\subset\mathbb{R}^{d}$, with a nonempty interior, not depending on $\theta$.

  2. 2.

The conditions in Theorem B.2 (that guarantee the CLT under shrinking
alternatives) hold.

  3. 3.

For all $\theta\in\Theta$, $f(\cdot|\theta)$ and
$\nabla_{\theta}f(\cdot|\theta)$ are three times continuously differentiable
in the interior of $S$, and the Fisher information matrix $\mathbb{E}_{X\sim
f(\cdot|\theta)}\left[\left(\frac{\nabla_{\theta}f(X|\theta)}{f(X|\theta)}\right)\left(\frac{\nabla_{\theta}f(X|\theta)}{f(X|\theta)}\right)^{\top}\right]$
is positive definite for all $\theta\in\Theta$.

  4. 4.

For all $x\in S$, $f(x|\cdot)$ is three times continuously differentiable in
$\Theta$.

Suppose that $M=2$ and the discrete kernel $K$ and the $k$-NN graph (for a
fixed $k\geq 1$) are used in defining $\hat{\eta}$. Let $P_{1}$ and $P_{2}$
have densities $f(\cdot|\theta_{1})$ and $f(\cdot|\theta_{2})$ respectively,
and $\theta_{2}=\theta_{1}+\varepsilon_{N}$, where $\varepsilon_{N}\to 0$ as
$N\to\infty$. Let $\kappa_{1,{\rm null}}$ be as defined in Theorem B.2, and
${\rm H}_{x}f(x|\theta_{1})$ be the Hessian of $f(x|\theta_{1})$ taken w.r.t.
$x$. Also, define, for $h\in\mathbb{R}^{p}$,

$a_{k,\theta_{1}}(h):=-\frac{\pi_{1}\pi_{2}(1-2\pi_{2})C_{k,2}}{2kd\kappa_{1,{\rm
null}}}\int_{S}h^{\top}\nabla_{\theta_{1}}\left(\frac{{\rm tr}({\rm
H}_{x}f(x|\theta_{1}))}{f(x|\theta_{1})}\right)f^{\frac{d-2}{d}}(x|\theta_{1}){\rm
d}x,$ (B.8)
$b_{k,\theta_{1}}(h):=\frac{\pi_{1}\pi_{2}}{\sigma_{\mathcal{G},K,d,\pi}}\mathbb{E}_{X\sim
f(\cdot|\theta_{1})}\Bigg{[}\frac{h^{\top}\nabla_{\theta_{1}}f(X|\theta_{1})}{f(X|\theta_{1})}\Bigg{]}^{2},\
{\rm with}\
C_{k,2}:=\mathbb{E}\Bigg{[}\sum_{x\in\mathcal{P}_{1}^{\mathbf{0}}:(\mathbf{0},x)\in\mathcal{E}(\mathcal{G}_{k{\rm
NN}}(\mathcal{P}_{1}^{\mathbf{0}}))}\|x\|^{2}\Bigg{]}.$

  * 1.

If the dimension $d\leq 8$, then the following hold:

    * (a)

$\|N^{\frac{1}{4}}\varepsilon_{N}\|\to 0$: The limiting power of the test
(B.7) is $\alpha$.

    * (b)

$N^{\frac{1}{4}}\varepsilon_{N}\to h$:

If $d\leq 7$, limiting power of the test (B.7) is
$\Phi\left(z_{\alpha}+b_{k,\theta_{1}}(h)\right)$.

If $d=8$ and the limiting power is
$\Phi\left(z_{\alpha}+a_{k,\theta_{1}}(h)+b_{k,\theta_{1}}(h)\right)$.

    * (c)

$\|N^{\frac{1}{4}}\varepsilon_{N}\|\to\infty$: The limiting power of the test
(B.7) is 1.

  * 2.

If the dimension $d\geq 9$, then the following hold:

    * (a)

$\|N^{\frac{1}{2}-\frac{2}{d}}\varepsilon_{N}\|\to 0$: The limiting power of
the test (B.7) is $\alpha$.

    * (b)

$N^{\frac{1}{2}-\frac{2}{d}}\varepsilon_{N}\to h$: The limiting power of the
test (B.7) is $\Phi\left(z_{\alpha}+a_{k,\theta_{1}}(h)\right)$.

    * (c)

$\|N^{\frac{1}{2}-\frac{2}{d}}\varepsilon_{N}\|\to\infty$ such that
$\|N^{\frac{2}{d}}\varepsilon_{N}\|\to 0$: Then depending on whether

$N^{\frac{1}{2}-\frac{2}{d}}(1-2\pi_{2})\int_{S}\varepsilon_{N}^{\top}\nabla_{\theta_{1}}\left(\frac{{\rm
tr}({\rm
H}_{x}f(x|\theta_{1}))}{f(x|\theta_{1})}\right)f^{\frac{d-2}{d}}(x|\theta_{1}){\rm
d}x\to\left\\{\begin{aligned} &\infty,\\\ &-\infty,\\\ \end{aligned}\right.$
(B.9)

the limiting power of the test (B.7) is 0 or 1, respectively.

    * (d)

$N^{\frac{2}{d}}\varepsilon_{N}\to h$: The limiting power of the test (B.7) is
0 or 1, depending on whether $a_{k,\theta_{1}}(h)+b_{k,\theta_{1}}(h)$ is
negative or positive, respectively.

    * (e)

$N^{\frac{2}{d}}\varepsilon_{N}\to\infty$: The limiting power of the test
(B.7) is 1.

A proof of the above result is given in Appendix C.14. In particular, the
above result shows that the detection threshold of $\hat{\eta}$ exhibits a
“$d=8$” phenomenon (see Bhattacharya [11]): When $d\leq 8$, the detection
threshold is $N^{-\frac{1}{4}}$; while when $d\geq 9$, the detection threshold
is somewhere between $N^{-\frac{2}{d}}$ and $N^{-\frac{1}{2}+\frac{2}{d}}$,
depending on the direction of $\varepsilon_{N}$ and the sign of
$a_{k,\theta_{1}}$; see (B.8) and (B.9). If $\varepsilon_{N}=\delta_{N}h$ for
$\delta_{N}>0$ and some fixed $h\in\mathbb{R}^{d}$, and
$a_{k,\theta_{1}}(h)>0$, then the detection threshold is $\delta_{N}\sim
N^{-\frac{1}{2}+\frac{2}{d}}$; on the other hand, if $a_{k,\theta_{1}}(h)<0$,
then the detection threshold is $\delta_{N}\sim N^{-\frac{2}{d}}$. When
$a_{k,\theta_{1}}(h)=0$, the precise location of the detection threshold has
to be determined on a case by case basis (see [11]). In Appendix D.2, we
empirically show that the same detection threshold also holds for the non-
Poissonized setting, i.e., the original setting in Section 2 of our paper
where the $n_{i}$’s are nonrandom constants (instead of $n_{i}\sim{\rm
Poisson}(N_{i})$).

##  Appendix C Proofs of the Main Results

###  C.1 Proof of Lemma 1

If $P_{1}=\ldots=P_{M}=:P$, then $\tilde{Z}$ follows $P$ regardless of
$\Delta$, and so $\tilde{\Delta}\perp\\!\\!\\!\perp\tilde{Z}$. Conversely, if
$\tilde{\Delta}\perp\\!\\!\\!\perp\tilde{Z}$, then the conditional
distribution $\tilde{Z}|(\tilde{\Delta}=i)$ should be the same for all
$i=1,\ldots,M$, which implies $P_{1}=\ldots=P_{M}$. Hence the $M$
distributions are the same if and only if
$\tilde{\Delta}\perp\\!\\!\\!\perp\tilde{Z}$.

If there exist disjoint measurable sets $\\{A_{i}\\}_{i=1}^{M}$ such that
$P_{i}(A_{i})=1$, then we can almost surely determine the label of
$\tilde{Z}$, by finding the $A_{i}$ that $\tilde{Z}$ belongs to, i.e.,
$\tilde{\Delta}=\sum_{i=1}^{M}iI(\tilde{Z}\in A_{i})$ can be written as a
measurable function of $\tilde{Z}$. Conversely, if $\tilde{\Delta}$ is a
measurable function of $\tilde{Z}$, then the
$\mathcal{E}(\tilde{Z})$-measurable101010$\mathcal{E}(\tilde{Z})$ is the
smallest $\sigma$-algebra such that $\tilde{Z}$ is measurable. set $A_{i}$
defined as $\tilde{\Delta}^{-1}(i)$ satisfies
$P_{i}(A_{i})=\mathbb{P}(\tilde{Z}\in A_{i}|\tilde{\Delta}=i)=1$. Hence there
exist disjoint measurable sets $\\{A_{i}\\}_{i=1}^{M}$ such that
$P_{i}(A_{i})=1$ if and only $\tilde{\Delta}$ is a measurable function of
$\tilde{Z}$. ∎

###  C.2 Proof of Theorem 1

Note that $\eta$ is the kernel measure of association between $\tilde{Z}$ and
$\tilde{\Delta}$ proposed in [21, Theorem 2.1]. Given the kernel
$K(\cdot,\cdot)$ is characteristic, it satisfies:

  * (i)

$\eta$ is a deterministic number between $[0,1]$;

  * (ii)

$\eta=0$ if and only if $\tilde{\Delta}\perp\\!\\!\\!\perp\tilde{Z}$, and

  * (iii)

$\eta=1$ if and only if $\tilde{\Delta}$ is a measurable function of
$\tilde{Z}$.

Hence the results follow from Lemma 1. ∎

###  C.3 Proof of Proposition 2

Define
$(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}},\tilde{Z},\tilde{Z}^{\prime})$ as
follows:

  1. 1.

Draw $(\tilde{\Delta},\tilde{\Delta}^{\prime},\tilde{Z})$ according to the
distribution mentioned in Section 2.1 (see (2)), i.e.,

    1. (a)

Draw $\tilde{\Delta}\sim\sum_{i=1}^{M}\pi_{i}\delta_{i}$, where $\delta_{i}$
is the Dirac measure at $i$.

    2. (b)

Draw $\tilde{Z}\mid(\tilde{\Delta}=i)\sim P_{i}$.

    3. (c)

Draw $\tilde{\Delta}^{\prime}$ independently from the conditional distribution
of $\tilde{\Delta}\mid\tilde{Z}$, given the value of $\tilde{Z}$ from (b).

  2. 2.

Draw $\tilde{Z}^{\prime}$ from the distribution $\kappa(\tilde{Z},\cdot)$.

Then $\tilde{Z}^{\prime}\sim\sum_{i=1}^{M}\pi_{i}Q_{i}$, and
$\eta(Q_{1},\ldots,Q_{M})$ can be defined from
$(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}},\tilde{Z}^{\prime})$, as
$\eta(P_{1},\ldots,P_{M})$ can be defined from
$(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}},\tilde{Z})$.

The denominator and the second term of the numerator of $\eta$ in (2) are
unchanged when $P_{1},\ldots,P_{M}$ are transitioned to $Q_{1},\ldots,Q_{M}$
as these only depend on the mixture proportion $\\{\pi_{i}\\}_{i=1}^{M}$ which
is unchanged. Hence it suffices to consider the first term in the numerator
and show:

$\mathbb{E}\left[\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|\tilde{Z}\big{]}\right]\geq\mathbb{E}\left[\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|\tilde{Z}^{\prime}\big{]}\right],$
(C.1)

where on the left-hand side, $\tilde{\Delta}$ and
$\tilde{\Delta}{{}^{\prime}}$ are i.i.d. drawn from the conditional
distribution of $\tilde{\Delta}\mid\tilde{Z}$, and on the right-hand side,
$\tilde{\Delta}$ and $\tilde{\Delta}{{}^{\prime}}$ are i.i.d. drawn from the
conditional distribution of $\tilde{\Delta}\mid\tilde{Z}^{\prime}$.

We first consider the case where $\tilde{\Delta}$ is not a deterministic
function of $\tilde{Z}^{\prime}$. Consider the kernel partial correlation [42]
between $\tilde{\Delta}$ and $\tilde{Z}$ given $\tilde{Z}^{\prime}$, defined
by

$\rho^{2}(\tilde{\Delta},\tilde{Z}|\tilde{Z}^{\prime})=\frac{\mathbb{E}\left[\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|(\tilde{Z},\tilde{Z}^{\prime})\big{]}\right]-\mathbb{E}\left[\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|\tilde{Z}^{\prime}\big{]}\right]}{\mathbb{E}\left[K(\tilde{\Delta},\tilde{\Delta})\right]-\mathbb{E}\left[\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|\tilde{Z}^{\prime}\big{]}\right]},$

where in
$\mathbb{E}\left[\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|(\tilde{Z},\tilde{Z}^{\prime})\big{]}\right]$,
$\tilde{\Delta}$ and $\tilde{\Delta}{{}^{\prime}}$ are i.i.d. drawn from the
conditional distribution of $\tilde{\Delta}\mid(\tilde{Z},\tilde{Z}^{\prime})$
and then averaged over $(\tilde{Z},\tilde{Z}^{\prime})$. Since
$\tilde{Z}^{\prime}$ is drawn from $\kappa(\tilde{Z},\cdot)$,
$\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|(\tilde{Z},\tilde{Z}^{\prime})\big{]}=\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|\tilde{Z}\big{]}$,
and therefore

$\rho^{2}(\tilde{\Delta},\tilde{Z}|\tilde{Z}^{\prime})=\frac{\mathbb{E}\left[\mathbb{E}\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|\tilde{Z}\big{]}\right]-\mathbb{E}\left[\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|\tilde{Z}^{\prime}\big{]}\right]}{\mathbb{E}\left[K(\tilde{\Delta},\tilde{\Delta})\right]-\mathbb{E}\left[\big{[}K(\tilde{\Delta},\tilde{\Delta}{{}^{\prime}})|\tilde{Z}^{\prime}\big{]}\right]}.$

It is shown in Huang et al. [42, Lemma 2] that the denominator of
$\rho^{2}(\tilde{\Delta},\tilde{Z}|\tilde{Z}^{\prime})$ is positive, and the
numerator of $\rho^{2}(\tilde{\Delta},\tilde{Z}|\tilde{Z}^{\prime})$ is
nonnegative. Hence inequality (C.1) holds. By Huang et al. [42, Theorem 1],
$\rho^{2}(\tilde{\Delta},\tilde{Z}|\tilde{Z}^{\prime})=0$ if and only if
$\tilde{\Delta}$ is conditionally independent of $\tilde{Z}$ given
$\tilde{Z}^{\prime}$, which is the same as saying
$\tilde{\Delta}\mid(\tilde{Z},\tilde{Z}^{\prime})\stackrel{{\scriptstyle
d}}{{=}}\tilde{\Delta}\mid\tilde{Z}^{\prime}$. Since $\tilde{Z}^{\prime}$ is
obtained by passing $\tilde{Z}$ through $\kappa$, this is further equivalent
to $\tilde{\Delta}\mid\tilde{Z}\stackrel{{\scriptstyle
d}}{{=}}\tilde{\Delta}\mid\tilde{Z}^{\prime}$.

If $\tilde{\Delta}$ is a deterministic function of $\tilde{Z}^{\prime}$, then
$\tilde{\Delta}$ must also be a deterministic function of $\tilde{Z}$, and the
inequality (C.1) holds with equality, and
$\tilde{\Delta}\mid\tilde{Z}\stackrel{{\scriptstyle
d}}{{=}}\tilde{\Delta}\mid\tilde{Z}^{\prime}$ holds as well.

In the case of a deterministic data processing, a function
$T:\mathcal{Z}\to\mathcal{Z}^{\prime}$ is applied to all the data points from
the $M$ distributions. In such a case, $\kappa(z,\cdot)=\delta_{T(z)}$, the
Dirac measure at $T(z)$. If $T$ is a bijection, then
$\tilde{\Delta}\mid\tilde{Z}\stackrel{{\scriptstyle d}}{{=}}\tilde{\Delta}\mid
T(\tilde{Z})=\tilde{\Delta}\mid\tilde{Z}^{\prime}$, so the equality in (C.1)
is attained, showing that a bijective transformation of the $M$ distributions
does not change our measure of dissimilarity $\eta$. ∎

###  C.4 Proof of Proposition 3

From (3), it suffices to show the convexity for the first term in the
numerator. Let $\mu=\sum_{i=1}^{M}(P_{i}+Q_{i})$ be a dominating measure. The
first term in the numerator of (3) can be re-written as:

$\int\sum_{i,j=1}^{M}K(i,j)\frac{\frac{{\rm d}(\pi_{i}P_{i})}{{\rm
d}\mu}}{\frac{{\rm d}\bar{P}}{{\rm d}\mu}}\frac{\frac{{\rm
d}(\pi_{j}P_{j})}{{\rm d}\mu}}{\frac{{\rm d}\bar{P}}{{\rm d}\mu}}\frac{{\rm
d}\bar{P}}{{\rm d}\mu}{\rm d}\mu(z)=\int\sum_{i,j=1}^{M}K(i,j)\frac{\frac{{\rm
d}(\pi_{i}P_{i})}{{\rm d}\mu}\frac{{\rm d}(\pi_{j}P_{j})}{{\rm
d}\mu}}{\frac{{\rm d}(\pi_{1}P_{1})}{{\rm d}\mu}+\ldots+\frac{{\rm
d}(\pi_{M}P_{M})}{{\rm d}\mu}}{\rm d}\mu(z).$

Note that $f:\mathbb{R}^{M}\to\mathbb{R}$,
$x\mapsto\frac{x^{\top}[K(i,j)]_{i,j=1}^{M}x}{1^{\top}x}$ is a convex function
on $1^{\top}x>0$ [15, Exercise 3.23]. Hence,

$\displaystyle\quad\int\sum_{i,j=1}^{M}K(i,j)\frac{\frac{{\rm
d}(\pi_{i}(\lambda P_{i}+(1-\lambda)Q_{i}))}{{\rm d}\mu}\frac{{\rm
d}(\pi_{j}(\lambda P_{j}+(1-\lambda)Q_{j}))}{{\rm d}\mu}}{\frac{{\rm
d}(\pi_{1}(\lambda P_{1}+(1-\lambda)Q_{1}))}{{\rm d}\mu}+\ldots+\frac{{\rm
d}(\pi_{M}(\lambda P_{M}+(1-\lambda)Q_{M}))}{{\rm d}\mu}}{\rm d}\mu(z)$
$\displaystyle=\int f\left(\frac{{\rm d}(\pi_{1}(\lambda
P_{1}+(1-\lambda)Q_{1}))}{{\rm d}\mu}(z),\ldots,\frac{{\rm d}(\pi_{M}(\lambda
P_{M}+(1-\lambda)Q_{M}))}{{\rm d}\mu}(z)\right){\rm d}\mu(z)$
$\displaystyle\leq\int\left[\lambda f\left(\frac{{\rm d}(\pi_{1}P_{1})}{{\rm
d}\mu},\ldots,\frac{{\rm d}(\pi_{M}P_{M})}{{\rm
d}\mu}\right)+(1-\lambda)f\left(\frac{{\rm d}(\pi_{1}Q_{1})}{{\rm
d}\mu},\ldots,\frac{{\rm d}(\pi_{M}Q_{M})}{{\rm d}\mu}\right)\right]{\rm
d}\mu(z)$ $\displaystyle=\lambda\int\sum_{i,j=1}^{M}K(i,j)\frac{\frac{{\rm
d}(\pi_{i}P_{i})}{{\rm d}\mu}\frac{{\rm d}(\pi_{j}P_{j})}{{\rm
d}\mu}}{\frac{{\rm d}(\pi_{1}P_{1})}{{\rm d}\mu}+\ldots+\frac{{\rm
d}(\pi_{M}P_{M})}{{\rm d}\mu}}{\rm d}\mu(z)$
$\displaystyle\quad+(1-\lambda)\int\sum_{i,j=1}^{M}K(i,j)\frac{\frac{{\rm
d}(\pi_{i}Q_{i})}{{\rm d}\mu}\frac{{\rm d}(\pi_{j}Q_{j})}{{\rm
d}\mu}}{\frac{{\rm d}(\pi_{1}Q_{1})}{{\rm d}\mu}+\ldots+\frac{{\rm
d}(\pi_{M}Q_{M})}{{\rm d}\mu}}{\rm d}\mu(z),$

and the joint convexity of $\eta$ follows. ∎

###  C.5 Proof of a Special Case of $\eta$

Suppose $\mathcal{Z}=\mathbb{R}^{d}$, $M=2$, and that $P_{1}$ and $P_{2}$ have
densities $f$ and $g$ w.r.t. the Lebesgue measure on $\mathbb{R}^{d}$. For
$z\in\mathcal{Z}$,

$\mathbb{P}(\tilde{\Delta}=2|Z=z)=\frac{\pi_{2}g(z)}{\pi_{1}f(z)+\pi_{2}g(z)}.$

Further, $\eta$ can be written as
$\eta=1-\frac{\mathbb{E}[\rho(\tilde{\Delta},\tilde{\Delta}^{\prime})]}{\mathbb{E}[\rho(\tilde{\Delta}_{1},\tilde{\Delta}_{2})]}$
(see [21, Equation 2.1]), where
$\rho(s,s^{\prime}):=\|K(s,\cdot)-K(s^{\prime},\cdot)\|_{\mathcal{H}}^{2}$,111111$\mathcal{H}$
is the RKHS induced by kernel $K$ on ${\mathcal{S}}:=\\{1,\ldots,M\\}$
equipped with inner product $\langle
K(s,\cdot),K(s^{\prime},\cdot)\rangle=K(s,s^{\prime})$, for
$s,s^{\prime}\in{\mathcal{S}}$. and

$\displaystyle\mathbb{E}[\rho(\tilde{\Delta},\tilde{\Delta}^{\prime})]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\mathbb{E}[\rho(\tilde{\Delta},\tilde{\Delta}^{\prime})|Z]\right]$
$\displaystyle=$
$\displaystyle\int\left\\{2\rho(1,2)\frac{\pi_{1}\pi_{2}f(x)g(x)}{[\pi_{1}f(x)+\pi_{2}g(x)]^{2}}\right\\}[\pi_{1}f(x)+\pi_{2}g(x)]\mathrm{d}x$
$\displaystyle=$ $\displaystyle
2\rho(1,2)\int\frac{\pi_{1}\pi_{2}f(x)g(x)}{\pi_{1}f(x)+\pi_{2}g(x)}\mathrm{d}x,$

and

$\mathbb{E}\big{[}\rho(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\big{]}=2\rho(1,2)\,\pi_{1}\pi_{2}.$

Therefore, $\eta=1-\int\frac{f(x)g(x)}{\pi_{1}f(x)+\pi_{2}g(x)}\mathrm{d}x.$
This completes the argument. ∎

###  C.6 Proof of Proposition 1

We will use the expression in Appendix C.5 of $\eta$, i.e.,
$\eta=1-\int\frac{f(x)g(x)}{\pi_{1}f(x)+\pi_{2}g(x)}\mathrm{d}x.$. First
consider the location family. Here, we have:

$\eta=1-\int\frac{1}{\pi_{1}\frac{1}{f(x)}+\pi_{2}\frac{1}{f(x-\theta)}}\mathrm{d}x.$

Since $f(x)$ is log-concave, $\frac{1}{f(x)}$ is convex as $\exp(\cdot)$ is
convex and increasing. Write $g:=\frac{1}{f}$. Then

$\displaystyle\eta$
$\displaystyle=1-\int\frac{1}{\pi_{1}g(x)+\pi_{2}g(x-\theta)}\mathrm{d}x$
$\displaystyle=1-\int\frac{1}{\pi_{1}g(x+\pi_{2}\theta)+\pi_{2}g(x-\pi_{1}\theta)}\mathrm{d}x$
$\displaystyle=1-\int\frac{1}{\pi_{1}g(x+\pi_{2}\lambda
h)+\pi_{2}g(x-\pi_{1}\lambda h)}\mathrm{d}x.$

It suffices to show that $\pi_{1}g(x+\pi_{2}\lambda
h)+\pi_{2}g(x-\pi_{1}\lambda h)$ is increasing in $\lambda$, i.e., for
$0\leq\lambda\leq\Lambda$,

$\pi_{1}g(x+\pi_{2}\lambda h)+\pi_{2}g(x-\pi_{1}\lambda
h)\leq\pi_{1}g(x+\pi_{2}\Lambda h)+\pi_{2}g(x-\pi_{1}\Lambda h).$ (C.2)

Write $x-\pi_{1}\lambda h$ and $x+\pi_{2}\lambda h$ as convex combinations of
$x-\pi_{1}\Lambda h$ and $x+\pi_{2}\Lambda h$, i.e.,

$\displaystyle x-\pi_{1}\lambda h$
$\displaystyle=\frac{\pi_{1}(\Lambda-\lambda)}{\Lambda}(x+\pi_{2}\Lambda
h)+\frac{\pi_{1}\lambda+\pi_{2}\Lambda}{\Lambda}(x-\pi_{1}\Lambda h),$
$\displaystyle x+\pi_{2}\lambda h$
$\displaystyle=\frac{\pi_{1}\Lambda+\pi_{2}\lambda}{\Lambda}(x+\pi_{2}\Lambda
h)+\frac{\pi_{2}(\Lambda-\lambda)}{\Lambda}(x-\pi_{1}\Lambda h).$

By the convexity of $g$, we have:

$\displaystyle\pi_{1}g(x+\pi_{2}\lambda h)+\pi_{2}g(x-\pi_{1}\lambda h)$
$\displaystyle\leq\pi_{1}\left(\frac{\pi_{1}\Lambda+\pi_{2}\lambda}{\Lambda}g(x+\pi_{2}\Lambda
h)+\frac{\pi_{2}(\Lambda-\lambda)}{\Lambda}g(x-\pi_{1}\Lambda h)\right)$
$\displaystyle\qquad+\pi_{2}\left(\frac{\pi_{1}(\Lambda-\lambda)}{\Lambda}g(x+\pi_{2}\Lambda
h)+\frac{\pi_{1}\lambda+\pi_{2}\Lambda}{\Lambda}g(x-\pi_{1}\Lambda h)\right)$
$\displaystyle=\pi_{1}g(x+\pi_{2}\Lambda h)+\pi_{2}g(x-\pi_{1}\Lambda h).$

Hence the monotonicity is proved. When $\lambda=0$, $P_{1}$ and $P_{2}$ are
equal, so $\eta=0$. Since $f(\cdot)$ is a density,
$\int_{\mathbb{R}}f(x-\lambda h)\mathrm{d}\lambda<+\infty$ for Lebesgue almost
every $x$, which implies the convex function $\pi_{2}g(x-\lambda
h)=\pi_{2}\frac{1}{f(x-\lambda h)}\to+\infty$ as $\lambda\to+\infty$ for
Lebesgue almost every $x$. Hence, by dominated convergence,

$\eta=1-\int\frac{1}{\pi_{1}g(x)+\pi_{2}g(x-\lambda h)}\mathrm{d}x\to
1\quad\mathrm{as}\quad\lambda\to+\infty.$

To show that the monotonicity is strict, we can first suppose that $f(x)=0$ on
the boundary of its support since the boundary of a convex set has 0 Lebesgue
measure, so it will not change the value of $\eta$. Since $\eta<1$,
$\pi_{1}g(x+\pi_{2}\lambda h)+\pi_{2}g(x-\pi_{1}\lambda h)$ is not $+\infty$
for Lebesgue almost every $x$. Hence we can find an $x$ such that the left-
hand side of (C.2) is finite. We want to further find an $x$ such that the
inequality in (C.2) is strict. Note that if equality holds in (C.2), then by
the convexity, $g(\cdot)$ must be a linear function on the line segment
$[x-\pi_{1}\Lambda h,x+\pi_{2}\Lambda h]$. Since $g$ is log-convex, it must be
a constant on $[x-\pi_{1}\Lambda h,x+\pi_{2}\Lambda h]$. Note that
$\lim_{x\to\infty}g(x)=\lim_{x\to\infty}\frac{1}{f(x)}=\infty$. Hence we can
move $x$ along the line passing through $x-\pi_{1}\Lambda h$ and
$x+\pi_{2}\Lambda h$ until inequality (C.2) becomes strict (and the left-hand
side is still finite). Now, since a convex function is continuous in the
interior of its domain, (C.2) with strict inequality can still hold in a
neighborhood of $x$, so that after integration we obtain
$\eta(P_{1},P_{2}(\lambda))<\eta(P_{1},P_{2}(\Lambda))$.

Now consider the case of a scale family. Note that

$\eta=1-\int\frac{1}{\pi_{1}\frac{1}{f(x)}+\pi_{2}\frac{1}{\lambda f(\lambda
x)}}\mathrm{d}x.$

By a change of variable $y=\lambda^{\pi_{2}}x$ we get

$\eta=1-\int\frac{1}{\pi_{1}\lambda^{\pi_{2}}g\left(\frac{y}{\lambda^{\pi_{2}}}\right)+\pi_{2}\frac{1}{\lambda^{\pi_{1}}}g(\lambda^{\pi_{1}}y)}\mathrm{d}y,$

where $g(x)=\frac{1}{f(x)}$ is still a convex function as before.

We will show that
$\phi(\lambda):=\pi_{1}\lambda^{\pi_{2}}g\left(\frac{y}{\lambda^{\pi_{2}}}\right)+\pi_{2}\frac{1}{\lambda^{\pi_{1}}}g(\lambda^{\pi_{1}}y)$
is decreasing on $(0,1]$, and increasing on $[1,\infty)$. We can suppose both
$\frac{y}{\lambda^{\pi_{2}}}$ and $\lambda^{\pi_{1}}y$ lie in the interior of
${\rm supp}(f)$, because if any of them is not in the interior of ${\rm
supp}(f)$, then if $\lambda\geq 1$, $\phi(\Lambda)=+\infty$ for any
$\Lambda>\lambda$ by the convexity of $g$, and similarly if $\lambda<1$, then
$\phi(\delta)=+\infty$ for any $0<\delta<\lambda$. Given that both the points
are in the interior of ${\rm supp}(f)$, we can take the derivative of $\phi$
to obtain

$\displaystyle\phi^{\prime}(\lambda)$
$\displaystyle=\pi_{1}\left(\pi_{2}\lambda^{\pi_{2}-1}g\left(\frac{y}{\lambda^{\pi_{2}}}\right)-\lambda^{\pi_{2}}\nabla
g\left(\frac{y}{\lambda^{\pi_{2}}}\right)^{\top}\frac{\pi_{2}y}{\lambda^{\pi_{2}+1}}\right)$
$\displaystyle\qquad\qquad+\pi_{2}\left(-\frac{\pi_{1}}{\lambda^{\pi_{1}+1}}g(\lambda^{\pi_{1}}y)+\frac{1}{\lambda^{\pi_{1}}}\nabla
g(\lambda^{\pi_{1}}y)^{\top}\pi_{1}\lambda^{\pi_{1}-1}y\right)$
$\displaystyle=\frac{\pi_{1}\pi_{2}}{\lambda}\left(\lambda^{\pi_{2}}g\left(\frac{y}{\lambda^{\pi_{2}}}\right)-\nabla
g\left(\frac{y}{\lambda^{\pi_{2}}}\right)^{\top}y-\frac{1}{\lambda^{\pi_{1}}}g(\lambda^{\pi_{1}}y)+\nabla
g(\lambda^{\pi_{1}}y)^{\top}y\right).$

Note that $\phi^{\prime}(1)=0$. Let
$\psi(\lambda):=\lambda^{\pi_{2}}g\left(\frac{y}{\lambda^{\pi_{2}}}\right)-\nabla
g\left(\frac{y}{\lambda^{\pi_{2}}}\right)^{\top}y-\frac{1}{\lambda^{\pi_{1}}}g(\lambda^{\pi_{1}}y)+\nabla
g(\lambda^{\pi_{1}}y)^{\top}y$. It suffices to show that
$\psi^{\prime}(\lambda)>0$. Observe that

$\displaystyle\psi^{\prime}(\lambda)$
$\displaystyle=\pi_{2}\lambda^{\pi_{2}-1}g\left(\frac{y}{\lambda^{\pi_{2}}}\right)-\lambda^{\pi_{2}}\nabla
g\left(\frac{y}{\lambda^{\pi_{2}}}\right)^{\top}\frac{\pi_{2}y}{\lambda^{\pi_{2}+1}}+y^{\top}\nabla^{2}g\left(\frac{y}{\lambda^{\pi_{2}}}\right)y\cdot\frac{\pi_{2}}{\lambda^{\pi_{2}+1}}$
$\displaystyle\qquad\qquad+\frac{\pi_{1}}{\lambda^{\pi_{1}+1}}g(\lambda^{\pi_{1}}y)-\frac{1}{\lambda^{\pi_{1}}}\nabla
g(\lambda^{\pi_{1}}y)^{\top}\pi_{1}\lambda^{\pi_{1}-1}y+\pi_{1}\lambda^{\pi_{1}-1}y^{\top}\nabla^{2}g(\lambda^{\pi_{1}}y)y$
$\displaystyle=\pi_{2}\lambda^{\pi_{2}-1}L\left(\frac{y}{\lambda^{\pi_{2}}}\right)+\frac{\pi_{1}}{\lambda^{\pi_{1}+1}}L\left(\lambda^{\pi_{1}}y\right),$

where $L(x):=g(x)-\nabla g(x)^{\top}x+x^{\top}\nabla^{2}g(x)x$. Since $f(x)$
is log-concave, $g(x)=\frac{1}{f(x)}=\exp(h(x))$ for some convex function
$h(x)$. Hence

$\displaystyle L(x)$ $\displaystyle=e^{h(x)}-e^{h(x)}\nabla
h(x)^{\top}x+x^{\top}\left(e^{h(x)}\nabla h(x)\nabla
h(x)^{\top}+e^{h(x)}\nabla^{2}h(x)\right)x$
$\displaystyle=e^{h(x)}\left(1-\nabla h(x)^{\top}x+(\nabla
h(x)^{\top}x)^{2}\right)+e^{h(x)}x^{\top}\nabla^{2}h(x)x$
$\displaystyle=\frac{3}{4}e^{h(x)}+e^{h(x)}\left(\frac{1}{2}-\nabla
h(x)^{\top}x\right)^{2}+e^{h(x)}x^{\top}\nabla^{2}h(x)x>0.$

Hence $\psi^{\prime}(\lambda)>0$ and the monotonicity is proved.

When $\lambda=1$, $P_{1}=P_{2}$, and so $\eta=0$. As $\lambda\to 0^{+}$, no
matter $f(0)=0$ or $f(0)>0$, we have
$\pi_{1}\frac{1}{f(x)}+\pi_{2}\frac{1}{\lambda f(\lambda x)}\to+\infty$; so
dominated convergence implies

$\eta=1-\int\frac{1}{\pi_{1}\frac{1}{f(x)}+\pi_{2}\frac{1}{\lambda f(\lambda
x)}}\mathrm{d}x\to 1.$

As $\lambda\to+\infty$, by a change of variable $y=\lambda x$ and again using
dominated convergence theorem, we get

$\eta=1-\int\frac{1}{\pi_{1}\frac{\lambda}{f\left(\frac{y}{\lambda}\right)}+\pi_{2}\frac{1}{f(y)}}\mathrm{d}y\to
1.$

This completes the proof of the result. ∎

###  C.7 Proof of Theorem 2

Here we provide a convenient proof of Theorem 2 using the consistency result
in Deb et al. [21] and a coupling of our data to a process where
$\Delta_{i}$’s are i.i.d. Multinoulli($\pi$) with
$\pi=(\pi_{1},\ldots,\pi_{M})$.

Suppose $X^{(m)}_{1},X^{(m)}_{2}\ldots\overset{i.i.d.}{\sim}P_{m}$,
$m=1,\ldots,M$. Let $\tilde{\Delta}_{i}\overset{i.i.d.}{\sim}{\rm
Multinoulli}(\pi)$. If $\tilde{\Delta}_{i}=m$, we draw $\tilde{Z}_{i}$
according to $P_{m}$. More specifically, for $i=1,\ldots,n$,

$\tilde{Z}_{i}:=\sum_{m=1}^{M}I(\tilde{\Delta}_{i}=m)X^{(m)}_{\sum_{j=1}^{i}I(\tilde{\Delta}_{j}=m)}.$

Let $\tilde{n}_{m}:=\sum_{i=1}^{n}I(\tilde{\Delta}_{i}=m)$, and
$\tilde{\mathcal{G}}_{n}$ be the graph constructed on the pooled sample
$(\tilde{Z}_{1},\ldots,\tilde{Z}_{n})$, which is
$X^{(1)}_{1},\ldots,X^{(1)}_{\tilde{n}_{1}},X^{(2)}_{1},\ldots,X^{(2)}_{\tilde{n}_{2}},\ldots,X_{1}^{(M)},\ldots,X^{(M)}_{\tilde{n}_{M}}$,
after proper permutation. Define:

$\tilde{\eta}:=\frac{\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{\tilde{d}_{i}}\sum\limits_{j:(\tilde{Z}_{i},\tilde{Z}_{j})\in\tilde{\mathcal{G}}_{n}}K(\tilde{\Delta}_{i},\tilde{\Delta}_{j})-\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K(\tilde{\Delta}_{i},\tilde{\Delta}_{j})}{\frac{1}{n}\sum_{i=1}^{n}K(\tilde{\Delta}_{i},\tilde{\Delta}_{i})-\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K(\tilde{\Delta}_{i},\tilde{\Delta}_{j})},$

where $\tilde{d}_{i}$ is the out-degree of vertex $i$ in the graph
$\tilde{\mathcal{G}}_{n}$.

From the proof of Theorem 1 we see that $\eta$ is the kernel measure of
association [21] between $\tilde{Z}$ and $\tilde{\Delta}$. Here $\tilde{\eta}$
is the empirical estimator of $\eta$ proposed in [21]. We can check the
conditions required for the consistency of $\tilde{\eta}$ (i.e.,
$\tilde{\eta}\stackrel{{\scriptstyle p}}{{\to}}\eta$) as follows: Note that
any function $f\in\mathcal{H}$121212$\mathcal{H}$ is the RKHS induced by
kernel $K$ on ${\mathcal{S}}:=\\{1,\ldots,M\\}$ equipped with inner product
$\langle K(s,\cdot),K(s^{\prime},\cdot)\rangle=K(s,s^{\prime})$, for
$s,s^{\prime}\in{\mathcal{S}}$. can be written as
$f(s)=f(1)I(s=1)+\ldots+f(M)I(s=M)$, for $s\in{\mathcal{S}}$, which implies
$\mathcal{H}$ is finite-dimensional and hence separable. As any kernel on the
finite set ${\mathcal{S}}$ is bounded, we have $\mathbb{E}_{\tilde{Z}\sim