# MMD Two-sample Testing in the Presence
of Arbitrarily Missing Data

Yijin Zeng  Niall M. Adams  Dean A. Bodenham

(Department of Mathematics, Imperial College London,
South Kensington Campus, London SW7 2AZ, U.K.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>)

###### Abstract

In many real-world applications, it is common that a proportion of the data
may be missing or only partially observed. We develop a novel two-sample
testing method based on the Maximum Mean Discrepancy (MMD) which accounts for
missing data in both samples, without making assumptions about the missingness
mechanism. Our approach is based on deriving the mathematically precise bounds
of the MMD test statistic after accounting for all possible missing values. To
the best of our knowledge, it is the only two-sample testing method that is
guaranteed to control the Type I error for both univariate and multivariate
data where data may be arbitrarily missing. Simulation results show that our
method has good statistical power, typically for cases where 5% to 10% of the
data are missing. We highlight the value of our approach when the data are
missing not at random, a context in which either ignoring the missing values
or using common imputation methods may not control the Type I error.

##  1 Introduction

Two-sample hypothesis testing is a fundamental statistical method used to
determine if two samples of data are different. Numerous two-sample testing
methods are available, including Student’s $t$-test [31], the Wilcoxon-Mann-
Whitney $U$ test [18], the Kolmogorov-Smirnov test [13], the Energy Distance
[23], and the Maximum Mean Discrepancy (MMD) test [9]. These methods have
proven useful in various fields including medicine [32, 20], finance [17, 12],
psychology [22] and machine learning [9, 10, 28]. Nearly all two-sample
testing methods are designed solely for samples that are fully observed.
However, in many real-world datasets, a subset of univariate values may be
missing or multivariate values may only be partially observed.

When data are missing, common practices are either to ignore all missing
values or impute these values using some imputation scheme. Following this
step, the data is treated as complete, allowing any standard two-sample
testing method to be used. However, except in special cases, such as when the
missing data are missing completely at random [24], these practices are often
invalid as they risk increasing the probability of a Type I error occurring.
Under certain assumptions, such as when the data are missing completely at
random, missing at random, or when the data are missing not at random but the
missingness mechanisms (the samples can be subject to different missingness
mechanisms) can be explicitly specified [27], then certain sophisticated
missing data methods exist, such as the expectation-maximization algorithm and
the multiple imputation method [26]. However, without these assumptions,
relying on these methods is fraught with risk.

Recently a two-sample testing method for data with univariate missing values
was proposed which does not rely on a specific imputation procedure [33]. This
method, based the Wilcoxon-Mann-Whitney test, rejects the null hypothesis
after accounting for all possible values of the missing data by deriving
bounds for the $U$ statistic. This approach avoids ignoring or imputing the
missing data and is shown to control the Type I error without making
assumptions about the missingness mechanisms, while also having good
statistical power when the proportion of missing data is around $10\%$.
However, this approach is restricted to univariate samples and can only detect
location shifts since it is based on the Wilcoxon-Mann-Whitney test.

In this paper, we propose a novel two-sample testing method based on the MMD
test in the presence of missing data that makes no assumption about the
missingness mechanisms, following the approach in [33]. However, our approach
can be used for both univariate and multivariate samples, and can detect any
distributional change since it is based on the MMD test.

Our approach, named MMD-Miss, is based on deriving the bounds of the unbiased
MMD test statistic [9] in the presence of missing data. To do so, we use the
Laplacian kernel, a popular choice of characteristic kernel that enables MMD
to detect any distributional shift [30]. To compute the $p$-value for the
test, we either use Monte Carlo sampling from the permuation distribution of
the observed data or, in the case of high-dimensional multivariate data, we
use the normality approximation recently derived in [8]. We prove that this
approach controls the Type I error, regardless of the values of the missing
data. Numerical simulations are presented, showing that our method has good
testing power when the proportion of missing data is around 5% to 10%, while
controlling the Type I error.

#### Two-sample testing for a small proportion of missing data

In the context of performing two-sample testing with missing data, a common
misunderstanding is that when the proportion of missing data is small, e.g. 5%
as suggested by [26, 11], the testing result will not be skewed by missing
data after either ignoring or imputing these values. While these practices may
be justified under certain conditions, such as when the data are missing
completely at random [24], we provide experiments in Section 5 showing that
these practices can lead to a Type I error asymptotically equals to 1 even
with a small proportion (5%) of missing data, when the data are missing not at
random.

##  2 Background

#### Related Work

Various approaches have been proposed to mitigate the issue of missing data
for two-sample hypothesis testing. When the missing data are missing at
random, rank-based two-sample testing methods have been proposed which control
the Type I error and achieve good testing power [5, 4]. Paired two-sample
testing with missing data are studied in [19, 6], under the assumption that
data are missing completely at random and at least one of the paired samples
is observed. When part of the samples are interval-censored, i.e. the samples
are bounded within an interval, [21] proposed an algorithm based on multiple
imputation [16] by employing approximate Bayesian bootstrap [25]. A method for
detecting a change in the mean of high-dimensional normal data assumes the two
samples share the same missingness mechanism [29]. By considering certain
unobserved samples as the “worst case”, [14] proposes a test method based on
rank called the worst rank test.

#### Missingness Mechanisms.

Missingness mechanisms are typically classified as missing completely at
random (MCAR), missing at random (MAR), and missing not at random (MNAR) [24].
Let $Z_{\text{com}}=(z_{1},\ldots,z_{n})^{T}$ denote the complete data, where
each $z_{i}$ is a $\mathbb{R}^{d}$ real value. The complete data
$Z_{\text{com}}$ can be split into observed parts $Z_{\text{obs}}$ and missing
parts $Z_{\text{miss}}$. Let $R$ be a binary indicator matrix matching the
dimensions of $Z_{\text{com}}$, with elements indicating observed (1) or
missing (0) values. The core idea of [24] is to admit $R$ as a probabilistic
phenomenon.

The mechanism is MCAR if $P\left(R|Z_{\text{com}}\right)=P\left(R\right)$,
meaning $R$ is independent of data values, justifying the practice of ignoring
missing data. It is MAR if
$P\left(R|Z_{\text{com}}\right)=P\left(R|Z_{\text{obs}}\right)$, where $R$
depends only on observed data, making certain imputation methods viable. If
neither condition holds, the mechanism is MNAR. For known mechanisms [24, 14],
specific imputations might be justified. However, with unknown mechanisms,
ignoring or imputing data can lead to invalid results, as shown in our
simulations (Section 5). Notably, our method does not depend on the missing
data values, ensuring valid results regardless of the missingness mechanism.

#### Two-Sample Testing.

A two-sample testing method is used to determine whether two groups of data
are statistically significantly different. Consider two samples of
observations $X=\\{x_{1},\ldots,x_{n_{1}}\\}$ and
$Y=\\{y_{1},\ldots,y_{n_{2}}\\}$, where each observation is in
$\mathbb{R}^{d}$, and suppose that $X$, $Y$ are independent random samples
from distributions $p$, $q$, respectively. We define the null hypothesis
$H_{0}$ as $p=q$ and the alternative hypothesis $H_{1}$ as $p\neq q$. The two-
sample testing method first computes a statistic based on the data, and then a
$p$-value based on the statistic. After comparing the $p$-value to a pre-
specified significance threshold $\alpha$, the null hypothesis $H_{0}$ is
either rejected or fails to be rejected.

A Type I error occurs when the null hypothesis $H_{0}$ is true but the test
incorrectly rejects it. Conversely, a Type II error occurs when $H_{0}$ is
false but the test fails to reject it. A two-sample testing method is usually
derived so that, under $H_{0}$, the Type I error is not larger than the pre-
specified significance threshold $\alpha\in(0,1)$. In this case, we say the
test controls the Type I error. For a given $\alpha$, a two-sample testing
method is preferred if it has a lower probability of making a Type II error,
denoted by $\beta$. The power is defined as $1-\beta$, and a two-sample
testing method is more desirable if it demonstrates higher power. Two-sample
testing methods are usually assessed on their power, given that the Type I
error is controlled.

#### Maximum Mean Discrepancy.

The Maximum Mean Discrepancy (MMD) is a kernel-based measure of distance
between two distributions $p$ and $q$, by comparing their mean embeddings in a
reproducing kernel Hilbert space [1]. More formally, let $\mathcal{H}_{k}$
denote a reproducing kernel Hilbert space with kernel function $k$, the MMD
between $p$ and $q$ is defined as [9]

$\displaystyle\text{MMD}^{2}[\mathcal{H}_{k},p,q]=\sup_{f\in\mathcal{H}_{k}:||f||\leq
1}(\mathbb{E}_{X\sim p}[f(X)]-\mathbb{E}_{Y\sim q}[f(Y)]).$

In the context of two-sample testing, an unbiased estimate of
$\text{MMD}^{2}[\mathcal{H}_{k},p,q]$ is:

$\displaystyle\begin{split}{\text{MMD}}^{2}_{u}(X,Y)&=\tfrac{1}{n_{1}(n_{1}-1)}\sum_{i=1}^{n_{1}}\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{n_{2}}k(x_{i},x_{j})+\tfrac{1}{n_{2}(n_{2}-1)}\sum_{i=1}^{n_{2}}\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{n_{2}}k(y_{i},y_{j})-\tfrac{2}{n_{1}n_{2}}\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}k(x_{i},y_{j}).\end{split}$

The MMD test rejects the null hypothesis when ${\text{MMD}}^{2}_{u}(X,Y)$
exceeds certain threshold, usually chosen to control the probability of Type I
error under a pre-specified number $\alpha\in(0,1)$. It is shown in [9] that
when $k$ is a characteristic kernel [30],
$\text{MMD}^{2}[\mathcal{H}_{k},p,q]=0$ if and only if $p=q$, otherwise
$\text{MMD}^{2}[\mathcal{H}_{k},p,q]>0$. Hence, when a characteristic kernel
is used, the MMD test is able to detect any distribution changes. However, on
finite sample sizes, the choices of kernel have significant impact on the
power of MMD test [2]. Two popular choices of characteristic kernels are
Laplacian and Gaussian kernels, which are defined as

$\displaystyle
k_{L}({x},{y})=\exp(-\beta||x-y||_{1}),~{}~{}k_{G}({x},{y})=\exp(-\gamma||x-y||_{2}),$

respectively, where $\gamma>0$ and $\beta>0$ are hyperparameters.

#### Permutation Test.

The permutation test is a common numerical procedure for approximating the
distribution of ${\text{MMD}}^{2}_{u}(X,Y)$ under the null hypothesis and
deciding the test threshold. It is proved to be able to control the Type I
error of MMD test non-asymptotically [28]. Let
$X=\\{x_{1},\ldots,x_{n_{1}}\\}$ and $Y=\\{y_{1},\ldots,y_{n_{2}}\\}$, where
each sample is $\mathbb{R}^{d}$ real value. The permutation test begins with
uniformly sampling $B$ i.i.d. permutations from $\\{1,\ldots,n_{1}+n_{2}\\}$.
Let us denote the $B$ samplings as $(\sigma^{(1)},\ldots,\sigma^{(B)})$, where
each $\sigma^{(i)}=(\sigma^{(i)}(1),\ldots,\sigma^{(i)}(n_{1}+n_{2}))$,
$i=1,\ldots,B$. Subsequently, define
$z_{1}=x_{1},\ldots,z_{n_{1}}=x_{n_{1}},z_{n_{1}+1}=y_{1},\ldots,z_{n_{1}+n_{2}}=y_{n_{2}}$
and for any $i=1,\ldots,B$, denote

$\displaystyle
X_{\sigma^{(i)}}=\\{z_{\sigma^{(i)}(1)},\ldots,z_{\sigma^{(i)}(n_{1})}\\},Y_{\sigma^{(i)}}=\\{z_{\sigma^{(i)}(n_{1}+1)},\ldots,z_{\sigma^{(i)}(n_{1}+n_{2})}\\}.$

The test threshold of MMD test, using ${\text{MMD}}^{2}_{u}(X,Y)$ as test
statistic with significance level $0<\alpha<1$, is then computed as the
$\left\lceil\alpha(B+1)\right\rceil$-th largest numbers in the set

$\displaystyle\left\\{{\text{MMD}}^{2}_{u}(X,Y)\right\\}\cup\left\\{{\text{MMD}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}}),i=1,\ldots,B\right\\}.$

If ${\text{MMD}}^{2}_{u}(X,Y)$ exceeds this threshold, the MMD test rejects
the null hypothesis. Otherwise, the null hypothesis is not rejected.

#### Normality Approximation.

It is well known that for fixed dimension $d$, the asymptotic distribution of
${\text{MMD}}^{2}_{u}(X,Y)$ under the null hypothesis (i.e. $p=q$) takes the
form of an infinite weighted sum of $\chi^{2}$ random variables with the
weights depending on $p$ [9], which are normally unknown. An estimation of the
asymptotic distribution of ${\text{MMD}}^{2}_{u}(X,Y)$ is establish in [10],
which is proved asymptotically correct regardless of the distribution $p$
under certain conditions. However, empirical experiments show this estimation
is rather conservative.

It was recently shown [8] that the studentized ${\text{MMD}}^{2}_{u}(X,Y)$
given by Equation (3) convergences to normal distribution when all
$n,m,d\to\infty$. This result holds for a wide range of kernels, including
Laplacian and Gaussian kernels. This normality approximation is proved
empirically rather accurate [8] even for relatively small $n,m,d$, e.g.
$n=m=25,d=50$, and the higher approximation accuracy is shown for larger
$n,m,d$.

##  3 Bounding MMD with Missing Data

In this section, we provide bounds for the MMD test statistic in the presence
of missing data. To do so, we use the Laplacian kernel which is a
characteristic kernel [7] and able to detect any change in distribution. In
the following, $k(x,y)$ will be used to denote the Laplacian kernel only.

#### Univariate Samples.

If not all samples in $X$ and $Y$ are observed, the
${\text{MMD}}^{2}_{u}(X,Y)$ statistic cannot be directly computed. However,
for any two given real values $x$ and $y$, the Laplacian kernel
$k(x,y)=\exp\left(-\beta|x-y|\right)$ is bounded by $(0,1]$ since $\beta>0$
and $|x-y|>0$. Hence, a straightforward way to obtain a lower bound of
${\text{MMD}}^{2}_{u}(X,Y)$ is to let $k(x_{i},x_{j})=0$ if at least one
sample of $x_{i}$ and $x_{j}$ is unobserved. Similarly, if at least one sample
of $y_{i}$ and $y_{j}$ is unobserved, let $k(x_{i},x_{j})=0$. On the other
hand, to minimize $-k(x_{i},y_{j})$, we can let $-k(x_{i},y_{j})=-1$ if either
$x_{i}$ or $y_{j}$ is unobserved. The upper bound of
${\text{MMD}}^{2}_{u}(X,Y)$ can be obtained following the similar manner. This
simplistic method, while providing an theoretically correct bound of
${\text{MMD}}^{2}_{u}(X,Y)$, yields bounds that are too conservative for
effective two-sample testing in practice. In the following, we construct
tighter bounds of ${\text{MMD}}^{2}_{u}(X,Y)$ in order to make it useful for
the two-sample testing problem. One starts by decomposing
${\text{MMD}}^{2}_{u}(X,Y)$ into terms which have either none, one or both
arguments of the kernel function missing; see Appendix A.1 for details. Given
this decomposition, one only needs to bound functions of the form in Equation
(1) below.

######  Lemma 1.

Let $x_{1},\ldots,x_{\ell_{1}}$ and $y_{1},\ldots,y_{\ell_{2}}$ be real
values, that are observed. Suppose $a_{1},\ldots,a_{\ell_{1}}$,
$b_{1},\ldots,b_{\ell_{2}}$ and $\beta$ are positive constants. Define

$\displaystyle
T(z)=\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta|x_{i}-z|)-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta|y_{i}-z|)$
(1)

as a function of $z\in\mathbb{R}$. Subsequently, for any given
$z_{0}\in\mathbb{R}$,

$\displaystyle
T(z_{0})\geq\min\\{0,T(x_{1}),\ldots,T(x_{\ell_{1}}),T(y_{1}),\ldots,T(y_{\ell_{2}})\\},$

and

$\displaystyle
T(z_{0})\leq\max\\{0,T(x_{1}),\ldots,T(x_{\ell_{1}}),T(y_{1}),\ldots,T(y_{\ell_{2}})\\}.$

This lemma is proved in Appendix A.2. It provides a linear time algorithm for
computing the bounds of $T(z)$ by computing all the values in the set
$\\{T(z):z\in\\{x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}\\}\\}$ and
identifying the maximum and minimum values.

#### Multivariate Samples.

We now consider the case where $X=\\{x_{1},\ldots,x_{n_{1}}\\}$,
$Y=\\{y_{1},\ldots,y_{n_{2}}\\}$ and $X,Y\in\mathbb{R}^{d}$ with $d>1$. For
each sample $x_{i}$ in $X$, denote $x_{i}=(x_{i}(1),\ldots,x_{i}(d))$, where
$x_{i}(l)\in\mathbb{R}$ is the value the $l$th component of $x_{i}$.
Similarly, for each sample $y_{i}$ in $Y$, denote
$y_{i}=(y_{i}(1),\ldots,y_{i}(d))$. We now provide notation for which
components of a $d$-dimensional observation are missing and define what it
means to impute such a value.

######  Definition 1.

For a value $z=(z(1),\dots,z(d))\in\mathbb{R}^{d}$, let $U_{z}$ denote the set
of components of $z$ that are missing; in other words $\\{z(j)\,:\,j\in
U_{z}\\}$ are missing.

######  Definition 2.

For a value $z\in\mathbb{R}^{d}$, that has missing components, $z^{\ast}$ is
called an _imputation_ of $z$ if $z^{\ast}\in\mathbb{R}^{d}$ is fully observed
and $z^{\ast}(j)=z(j)$ for all $j\in\\{1,\dots,d\\}\setminus U_{z}$.

We can now extend Lemma 1 to the $d$-dimensional case:

######  Lemma 2.

Let $x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}\in\mathbb{R}^{d}$ be
values that are fully observed. Suppose $a_{1},\ldots,a_{\ell_{1}}$,
$b_{1},\ldots,b_{\ell_{2}},\beta$ are positive constants. For
$z=(z(1),\ldots,z(d))\in\mathbb{R}^{d}$ with missing components, define

$\displaystyle T(\\{z(j)\,:\,j\in
U_{z}\\})=\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta\sum_{j\in
U_{z}}|x_{i}(j)-z(j)|\right)-\sum_{i=1}^{\ell_{2}}b_{i}\exp\left(-\beta\sum_{j\in
U_{z}}|y_{i}(j)-z(j)|\right)$

as a function of the unobserved components of $z$ and let

$\displaystyle\mathcal{X}=\\{T(\\{z(j)\,:\,j\in U_{z}\\})\
:\,z(i)\in\\{x_{1}(i),\ldots,x_{\ell_{1}}(i),y_{1}(i),\ldots,y_{\ell_{2}}(i)\\},i\in
U_{z}\\}.$

Then, for any possible imputation $z^{\ast}$ of $z$,

$\displaystyle\min\\{0,\min\mathcal{X}\\}\leq T(\\{z^{\ast}(j)\,:\,j\in
U_{z}\\})\leq\max\\{0,\max\mathcal{X}\\}$

Lemma 2 is proved in Appendix A.3. It shows that in order to compute the
maximum and minimum values of $T(\\{z(j)\,:\,j\in U_{z}\\})$ for
$z\in\mathbb{R}^{d}$, we only need to check the imputations of $z$ where its
missing components are imputed using the components of
$x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}$. However, computing
$T(\\{z(j)\,:\,j\in U_{z}\\})$ for all possible imputations using Lemma 2 is
$(\ell_{1}+\ell_{2})^{|U_{x_{i}}|}$, which is exponential in the number of
unobserved components of $z$ and impractical to compute. To address this
computational challenge, we propose to further bound $\min\mathcal{X}$ and
$\max\mathcal{X}$, using the following lemma:

######  Lemma 3.

Following the notation and definitions in Lemma 2, denote

$\displaystyle\tilde{x}_{i}(j)$
$\displaystyle=\max\\{|x_{i}(j)-x_{1}(j)|,\ldots,|x_{i}(j)-x_{\ell_{1}}(j)|,|x_{i}(j)-y_{1}(j)|,\ldots,|x_{i}(j)-y_{\ell_{2}}(j)|\\}$

for any $i\in\\{1,\ldots,\ell_{1}\\},j\in U_{z}$; denote

$\displaystyle\tilde{y}_{i}(j)$
$\displaystyle=\max\\{|y_{i}(j)-\ell_{1}(j)|,\ldots,|y_{i}(j)-x_{n_{1}}(j)|,|y_{i}(j)-y_{1}(j)|,\ldots,|y_{i}(j)-y_{\ell_{2}}(j)|\\}$

for any $i\in\\{1,\ldots,\ell_{2}\\},j\in U_{z}$. Subsequently,

$\displaystyle\min\mathcal{X}\geq\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta\sum_{j\in
U_{z}}\tilde{x}_{i}(j)\right)-\sum_{i=1}^{\ell_{2}}b_{i},~{}\max\mathcal{X}\leq\sum_{i=1}^{\ell_{1}}a_{i}-\sum_{i=1}^{\ell_{2}}b_{i}\exp\left(-\beta\sum_{j\in
U_{z}}\tilde{y}_{i}(j)\right).$

This lemma is proved in Appendix A.4. It can be seen that for each
$i\in\\{1,\ldots,\ell_{1}\\},j\in U_{z}$,
$|\tilde{x}_{i}(j)|=\ell_{1}+\ell_{2}$ according to its definition. Hence, to
bound $\min\mathcal{X}$, the computation times
$|U_{z}|\ell_{1}(\ell_{1}+\ell_{2})$ are needed, which is of order
${O}(dn_{1}(n_{1}+n_{2}))$. Similarly, to bound $\max\mathcal{X}$, the order
of computation complexity is ${O}(dn_{2}(n_{1}+n_{2}))$. This computational
cost reduction, compared with Lemma 2, which grows exponentially with
$|U_{z}|$, makes bounding ${\text{MMD}}^{2}_{u}(X,Y)$ possible in practice.

#### Bounds for MMD statistic

The bounds derived above do not depend on the missing data mechanisms.
Combining Lemmas 1, 2 and 3, we have our main result:

######  Theorem 1.

For data $x_{1},\ldots,x_{n_{1}}\in\mathbb{R}^{d}$ and
$y_{1},\ldots,y_{n_{2}}\in\mathbb{R}^{d}$ with missing values, bounds for the
MMD statistic with the Laplacian kernel can be computed in $O(n_{1}+n_{2})$
when $d=1$ and $O(d(n_{1}+n_{2})^{2})$ when $d>1$.

This theorem is proved by a series of results in Appendix A.5 and A.6.

##  4 Two-Sample Testing in the Presence of Arbitrarily Missing Data

In this section, we discuss methods for employing bounds derived in Section 3
to develop valid two-sample testing methods. Besides computing the test
statistic, we also need a method for computing a $p$-value from the test
statistic.

#### Bounding $p$-value using permutations.

The permutation test, as discussed in Section 2, is a numerical procedure that
can be used to a $p$-value for any statistic based on the observed data.
Instead of computing the threshold as in Section 2, an equivalent approach to
compute the $p$-value of the MMD test is

$\displaystyle
p=\frac{1}{B+1}\left(1+\sum_{i=1}^{B}\text{I}({\text{MMD}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\text{MMD}}^{2}_{u}(X,Y))\right).$
(2)

The null hypothesis is rejected if $p$ is smaller or equal to the pre-
specified significance level $\alpha$.

Using the bounds of MMD in the presence of missing data, we proceed by
providing bounds of the $p$-value, using the following theorem:

######  Theorem 2.

Suppose
$X=\\{x_{1},\ldots,x_{n_{1}}\\},Y=\\{y_{1},\ldots,y_{n_{2}}\\}\in\mathbb{R}^{d}$.
Let $(\sigma^{(1)},\ldots,\sigma^{(B)})$ be $B$ i.i.d. random permutations of
$\\{1,\ldots,n_{1}+n_{2}\\}$ and denoted as
$\sigma^{(i)}=(\sigma^{(i)}(1),\ldots,\sigma^{(i)}(n_{1}+n_{2}))$,
$i=1,\ldots,B$. Subsequently, let
$z_{1}=x_{1},\ldots,z_{n_{1}}=x_{n_{1}},z_{n_{1}+1}=y_{1},\ldots,z_{n_{1}+n_{2}}=y_{n_{2}}$
and for any $i=1,\ldots,B$, denote

$\displaystyle
X_{\sigma^{(i)}}=\\{z_{\sigma^{(i)}(1)},\ldots,z_{\sigma^{(i)}(n_{1})}\\},Y_{\sigma^{(i)}}=\\{z_{\sigma^{(i)}(n_{1}+1)},\ldots,z_{\sigma^{(i)}(n_{1}+n_{2})}\\},$

and define $p$ according to Equation (2). Suppose further

$\displaystyle{\underline{\text{MMD}}}^{2}_{u}(X,Y)\leq{\text{MMD}}^{2}_{u}(X,Y),~{}{\text{MMD}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\leq{\overline{\text{MMD}}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}}),~{}i=1,\ldots,B.$

Define

$\displaystyle\overline{p}=\frac{1}{B+1}\left(1+\sum_{i=1}^{B}\text{I}({\overline{\text{MMD}}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\underline{\text{MMD}}}^{2}_{u}(X,Y))\right).$

Then, we must have $\overline{p}\geq p$.

The proof of this theorem is included in Appendix A.7. In the presence of
missing data, ${\text{MMD}}^{2}_{u}(X,Y)$ cannot be computed directly.
However, as shown in Section 3, the lower and upper bounds of
${\text{MMD}}^{2}_{u}(X,Y)$ can be computed when the Laplacian kernel is used.
Then, Theorem 2 can be applied to compute $\overline{p}$, a value that is
bounded by $p$, defined in Equation (2). It is proved in [28] that, under the
null hypothesis, computing the $p$-value using permutations controls the Type
I error. Therefore, using the approach described in Theorem 2 will also
control the Type I error, since the bounds were derived without any
assumptions about the missing data mechanisms.

#### Bounding $p$-value using normality approximation.

In the above, we have proposed a testing procedure based on permutation test
in the presence of missing data. Notably, this procedure provides a valid test
procedure for any given $n_{1},n_{2},d$. When both $n_{1},n_{2},d$ are large
enough, an alternative method to compute a $p$-value is to use the normality
approximation [8]. It is suggested in [8] that this approximation is effective
when $n_{1},n_{2}\geq 25,d\geq 50$.

To more formally describe the normality approximation proposed in [8], we will
follow their notation in this section. Let us denote the sample sizes for $X$
and $Y$ as $n,m$, respectively, rather than $n_{1},n_{2}$. It is proved in [8]
the studentized ${\text{MMD}}^{2}_{u}(X,Y)$, taking the form as

$\displaystyle\frac{{\text{MMD}}^{2}_{u}(X,Y)}{\sqrt{c_{n,m}{\mathcal{V}}_{n,m}^{k*}}},$
(3)

where $c_{n,m}=2/(n(n-1))+4/(nm)+2/(m(m-1))$, and $\mathcal{V}_{n,m}^{k*}$ is
defined in Proposition 10 in [8], convergence to standard normal distribution
under the null (see Theorem 16 in [8]). In order to use this result, it is
necessary to compute the estimation of variance of
${\text{MMD}}^{2}_{u}(X,Y)$, i.e. $c_{n,m}\mathcal{V}_{n,m}^{k*}$, which can
not be computed directly with missing data. To overcome this problem, we
propose to first bound $\mathcal{V}_{n,m}^{k*}$, leading to the following
theorem:

######  Theorem 3.

A lower bound for the studentized MMD statistic in Equation 3 can be computed,
and this provides an upper bound for the $p$-value computed via the normality
approximation, when not all data are fully missing.

This theorem is proved in Appendix A.8, and provides an approach for using the
normality approximation to compute the $p$-value, when data are missing.

##  5 Experiments

In this section, we investigate the Type I error and statistical power of MMD-
Miss proposed in Section 4 and compare it with three common missing data
approaches: case deletion, mean imputation and hot deck imputation. These
other three approaches are described in further detail in Appendix B.2.

#### Asymptotic Type I error for a small proportion of missing data (5%).

The first experiment investigates the Type I error of the proposed methods and
common missing data approaches, with increasing sample sizes for a given
proportion (5%) of missing data, where data are missing not at random (MNAR).
To assess the Type I error, random samples $X=\\{x_{1},\ldots,x_{n_{1}}\\}$,
$Y=\\{y_{1},\ldots,y_{n_{2}}\\}$ are independently generated from a
$d$-dimensional normal distribution with mean vector $\mu=(0,\ldots,0)^{T}$
and covariance matrix equal to the identity matrix. Subsequently, $5\%$ of
samples will be selected in both $X$ and $Y$ to be labeled as missing (for
$d=1$), or incomplete (for $d>1$). When an observation is labeled as
incomplete, $30\%$ of its components will be missing values.

In this experiment, $d\in\\{1,10,50\\}$ will be considered. When the dimension
$d=1$, only observations $x_{i}\in X$ with $x_{i}<0$ will be randomly selected
to possibly be missing, while only observations $y_{i}\in Y$ with $y_{i}>0$
will be randomly selected to possibly be missing; in this case the data are
missing not at random (MNAR), i.e. the data are informatively missing. For
higher dimensions $d\in\\{10,50\\}$, the observations in $X$ that will
possibly be partially missing are those with
$\sum_{l=1}^{d}x_{i}(l)/\sqrt{d}<-0.8$, and then for each chosen sample,
$30\%$ of its components with values smaller than
$\text{median}{(x_{i}(1),\ldots,x_{i}(d))}$ will be randomly selected to be
missing. On the other hand, only observations in $Y$ with
$\sum_{l=1}^{d}y_{i}(l)/\sqrt{d}>0.8$ will possibly be partially missing, and
again for each chosen observation, $30\%$ its components with values larger
than $\text{median}{(y_{i}(1),\ldots,y_{i}(d))}$ will be randomly selected to
be missing. We use these missingness mechanisms, since if only the fully
observed observations are taken into account, the two samples will appear to
be different.

Since using the normality approximation with MMD-Miss may not be suitable when
$d\in\\{1,10\\}$, we only assess its performance when $d=50$, the dimension
suggested by [8] for when the approximation will be effective. Figure 1 shows
that the Type I error of all common missing data approaches increases with
increasing sample sizes. In particular, when the sample sizes equals to
$5000$, the Type I error approaches to 1. In contrast, MMD-Miss controls the
Type I error. The significance level is $\alpha=0.05$.

Figure 1: The Type I error of the proposed method and common missing data
methods when data are missing not at random. (Left): $d=1$; (Middle): $d=10$;
(Right): $d=50$. MMD-Miss only uses the normality approximation when $d=50$.
For all figures, the significance level is $\alpha=0.05$, and $5\%$ of the
data are missing or partially observed. The plotted values show the average
times of the null hypothesis is rejected over 100 repetitions. The error bars
represent one standard error of the mean.

#### Type I error and power for univariate data.

The second experiment considers the Type I error and power of MMD-Miss and
common missing data approaches for univariate samples where data are missing
not at random (MNAR). To assess the Type I error, the data
$X=\\{x_{1},\ldots,x_{n_{1}}\\}$ and $Y=\\{y_{1},\ldots,y_{n_{2}}\\}$ are
independently sampled from a standard normal distribution $N(0,1)$. To assess
the statistical power, the $X$ data are sampled independently from $N(0,1)$
the $Y$ data are sampled independently from either (a) $N(1,1)$ or (b)
$N(1.5,1)$. Then, a proportion $s\in\\{0,0.01,\ldots,0.20\\}$ of samples are
selected in $X$ and $Y$ to be missing. The missingness mechanisms are the same
as for Figure 1 when $d=1$; only observations $x_{i}\in X$ with $x_{i}<0$ and
observations $y_{i}\in Y$ with $y_{i}>0$ will be randomly selected to possibly
be missing.

Figure 2: The Type I error and power of the proposed method and common missing
data methods for univariate samples when data are missing not at random.
(Left): $N(0,1)$ vs $N(0,1)$; (Middle): $N(0,1)$ vs $N(1,1)$; (Right):
$N(0,1)$ vs $N(1.5,1)$. A significance level $\alpha=0.05$ and sample sizes
$n_{1}=n_{2}=500$ are used. The plotted values show the average times of the
null hypothesis is rejected over 100 repetitions. The error bars represent one
standard error of the mean.

Figure 2 shows that the Type I error is not controlled by the case deletion,
hot deck imputation or mean imputation methods. When the proportion of missing
data is $s=5\%$, the Type I error is above $25\%$ for each of these three
methods, and when the $s=10\%$ the Type I error of these three methods is
beyond $75\%$. On the other hand, the proposed method MMD-Miss controls the
Type I error, regardless of the proportion of missing data. While all the
common missing data methods have good statistical power, the power of the
proposed method decreases significantly when there are more than $5\%$ missing
data for alternative $N(0,1)$ vs $N(1,1)$, and when there are more than $10\%$
missing data for alternative $N(0,1)$ vs $N(1.5,1)$. This experiment
demonstrates that MMD-Miss is useful when the proportion of missing data is in
the range $5\%$ to $10\%$, while the three common missing data approaches fail
to control the Type I error.

#### Type I error and power for multivariate data.

The third experiment compares the Type I error and power of MMD-Miss and the
three common missing data approaches for multivariate observations. For a
value $a\in\mathbb{R}$, let $\mu_{a,d}=(a,a,\dots,a)^{T}$ be a mean vector
with $d$ components all equal to $a$. To assess the Type I error, observations
$X=\\{x_{1},\ldots,x_{n_{1}}\\}$ and $Y=\\{y_{1},\ldots,y_{n_{2}}\\}$ are
independently sampled from the $d$-dimensional normal distribution with zero
mean vector $\mu_{0,d}=(0,0,\dots,0)^{T}$ and covariance matrix
$\Sigma=I_{d}$, the identity matrix. To assess the statistical power, the $X$
data are independently sampled from $N(\mu_{0,d},I_{d})$ and either (a)
$N(\mu_{1,d},I_{d})$ or (b) $N(\mu_{1.5,d},I_{d})$. Then, a proportion
$s\in[0,0.2]$ of observations from both $X$ and $Y$ will be randomly selected
to be partially missing. The missingness mechanisms follows the same approach
used for Figure 1 when $d>1$, with each partially missing observation having
$30\%$ missing components. The results displayed in Figure 3 show a similar
pattern to the results for univariate case: (i) the three common missing data
approaches cannot control the Type I error although they all demonstrate good
power, while (ii) MMD-Miss controls the Type I error, while also have good
power when proportion of missing data is around $5\%$ to $10\%$. Notably, in
the case where $d=50$, MMD-Miss with the normality approximation has better
statistical power than when it uses permutations to compute the $p$-value. For
the alternative $N((1.5,\ldots,1.5)^{T},I_{d})$, MMD-Miss with the normality
approximation has good power when $15\%$ of the data have missing values.

Figure 3: The Type I error and power of MMD-Miss and the three common missing
data approaches for multivariate samples when data are missing not at random.
when the dimension of the data is either $d=10$ or $d=50$. When $d=50$, the
normality approximation can be used with MMD-Miss. A significance level
$\alpha=0.05$, and sample sizes $n_{1}=n_{2}=500$ are used. The plotted values
show the average times of the null hypothesis is rejected over 100
repetitions. The error bars represent one standard error of the mean.

#### MNIST dataset.

We evaluate the performance of MMD-Miss on real-world data using MNIST images
[15], with examples shown in Figure 4 in Appendix B.3. Each image in the MNIST
dataset has dimensions of $28\times 28$ pixels and is labeled from $0$ to $9$.
For our analysis, the pixel values of each image are scaled between $0$ and
$1$. To assess Type I error, datasets $X$ and $Y$ are generated by randomly
sampling with replacement from MNIST training set images labeled as $3$. For
evaluating power, $X$ is generated by randomly sampling with replacement from
images labeled as $0$, while $Y$ continues to be sampled from images labeled
as $3$. A proportion $s\in[0,0.25]$ of the samples in $Y$ are then randomly
selected and labeled as incomplete if there are more than 85 non-zero pixels
in the region defined by rows 1 to 14 and columns 8 to 21 (i.e. a sub-block in
the upper half of each image), which will be marked as missing. In other
words, images with more non-zero pixels in the specified region are more
likely to be partially observed. Examples of these incompletely observed
images are shown in Figure 5 in Appendix B.3 Table 1 shows that the three
common missing data approaches cannot control the Type I error for this task.
On the other hand, MMD-Miss controls the Type I error while enjoying good
power. For MMD-Miss based on normality approximation, the power is good except
for the case when $25\%$ images from $Y$ are missing.

Table 1: Comparison of Type I error and power on MNIST dataset. NA stands for MMD-Miss with the normality approximation; PT for MMD-Miss with a permutation test; CD stands for case deletion; MI stands for mean imputation; HD stands for hot deck imputation. A significance level $\alpha=0.05$, and sample sizes $n_{1}=n_{2}=500$ are used. The values in the table are the average times the null hypothesis is rejected over 100 repetitions. Proportion of Missing Data ($s$) | Type I Error | Power
---|---|---
MMD-Miss | Common | MMD-Miss | Common
NA | PT | CD | MI | HD | NA | PT | CD | MI | HD
0.00 | 0.07 | 0.07 | 0.06 | 0.07 | 0.07 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
0.05 | 0.00 | 0.00 | 0.08 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
0.10 | 0.00 | 0.00 | 0.21 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
0.15 | 0.00 | 0.00 | 0.64 | 1.00 | 1.00 | 1.00 | 0.00 | 1.00 | 1.00 | 1.00
0.20 | 0.00 | 0.00 | 0.91 | 1.00 | 1.00 | 0.96 | 0.00 | 1.00 | 1.00 | 1.00
0.25 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00

##  6 Conclusion

The proposed method MMD-Miss can perform two-sample testing on both univariate
and multivariate data with missing values. It has good statistical power,
typically when $5\%$ to $10\%$ of the data are missing, while controlling the
Type I error.

#### Limitations

First, a limitation of MMD-Miss is that it is restricted to the Laplacian
kernel for computational reasons. However since the Laplacian kernel is
characteristic, MMD with this kernel can detect any distributional change.
Moreover, MMD-Miss will be effective for other kernels if appropriate bounds
can be derived and computed efficiently. Second, MMD-Miss is only suitable
when up to $20\%$ of the values are not fully observed, and typically performs
best when $5\%$ to $10\%$ of the data are missing. However, using imputation
methods when even $5\%$ of the data are informatively missing leads to the
Type I error being out of control.

##  7 Acknowledgements

Yijin Zeng is funded by a Roth Studentship from the Department of Mathematics,
Imperial College London and the EPSRC CDT in Statistics and Machine Learning.

## Overview of Appendix

Appendix A. Mathematical details of developing bounds of
${\text{MMD}}^{2}_{u}(X,Y)$ in Section 3 and constructing valid two-sample
testing method in Section 4 are given.

Appendix B. Details of experiments in Section 5 are provided.

Appendix C. Additional experiments are given for further investigating the
power of the proposed methods.

##  Appendix A Mathematical Details

###  A.1 Decomposition of MMD

Without loss of generality, let us assume $x_{1},\ldots,x_{m_{1}}$,
$y_{1},\ldots,y_{m_{2}}$ are samples that are not observed, or not fully
observed. Then, ${\text{MMD}}^{2}_{u}(X,Y)$ can be decomposed, using the
following lemma:

######  Lemma 4.

Suppose $X=\\{x_{1},\ldots,x_{n_{1}}\\}$, $Y=\\{y_{1},\ldots,y_{n_{2}}\\}$ are
subsets of $\mathbb{R}^{d}$, where $d\geq 1$. Suppose that
${\text{MMD}}^{2}_{u}(X,Y)$ is unbiased MMD test statistic defined as

$\displaystyle\begin{split}{\text{MMD}}^{2}_{u}(X,Y)&=\tfrac{1}{n_{1}(n_{1}-1)}\sum_{i=1}^{n_{1}}\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{n_{2}}k(x_{i},x_{j})+\tfrac{1}{n_{2}(n_{2}-1)}\sum_{i=1}^{n_{2}}\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{n_{2}}k(y_{i},y_{j})-\tfrac{2}{n_{1}n_{2}}\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}k(x_{i},y_{j}),\end{split}$

with $k$ denoting the Laplacian kernel. Then, for any two positive integers
$m_{1}$, $m_{2}$ such that $m_{1}~{}\leq~{}n_{1},~{}m_{2}~{}\leq~{}n_{2}$,
${\text{MMD}}^{2}_{u}(X,Y)$ can be rewritten as:

$\displaystyle{\text{MMD}}^{2}_{u}(X,Y)=A_{1}+A_{2}+A_{3}+A_{4},$

where $c_{1}=\frac{2}{n_{1}(n_{1}-1)}$, $c_{2}=\frac{2}{n_{2}(n_{2}-1)}$,
$c_{3}=\frac{2}{n_{1}n_{2}}$, and

$\displaystyle A_{1}$
$\displaystyle=c_{1}\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k(x_{i},x_{j})+c_{2}\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j}),$
$\displaystyle A_{2}$
$\displaystyle=c_{1}\sum_{i=m_{1}+1}^{n_{1}-1}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j})+c_{2}\sum_{i=m_{2}+1}^{n_{2}-1}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j}),$
$\displaystyle A_{3}$
$\displaystyle=c_{1}\sum_{i=1}^{m_{1}}\sum_{j=m_{1}+1}^{n_{1}}k(x_{i},x_{j})-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j}),$
$\displaystyle A_{4}$
$\displaystyle=c_{2}\sum_{i=1}^{m_{2}}\sum_{j=m_{2}+1}^{n_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j}).$

###### Proof.

To start, for any $x$, $y$ of $\mathbb{R}^{d}$ real values, the Laplacian
kernel is defined as follows

$\displaystyle k({x},{y})=\exp(-\beta||x-y||_{1}).$

It can be seen that $k$ is symmetric, i.e. $k(x,y)=k(y,x)$. Subsequently,
according to Lemma 2 in [3], the ${\text{MMD}}^{2}_{u}(X,Y)$ can then be
rewritten as

$\displaystyle{\text{MMD}}^{2}_{u}(X,Y)$
$\displaystyle:=\frac{2}{n_{1}(n_{1}-1)}\sum_{i=1}^{n_{1}}\sum_{j=1}^{i-1}k(x_{i},x_{j})+\frac{2}{n_{2}(n_{2}-1)}\sum_{i=1}^{n_{2}}\sum_{j=1}^{i-1}k(y_{i},y_{j})$
$\displaystyle-\frac{2}{n_{1}n_{2}}\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}k(x_{i},y_{j}).$

Notice that

$\displaystyle\sum_{i=1}^{n_{1}}\sum_{j=1}^{i-1}k(x_{i},x_{j})$
$\displaystyle=\sum_{i=1}^{n_{1}-1}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j}),$
$\displaystyle\sum_{i=1}^{n_{2}}\sum_{j=1}^{i-1}k(y_{i},y_{j})$
$\displaystyle=\sum_{i=1}^{n_{2}-1}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j}).$

Thus, we further have

$\displaystyle{\text{MMD}}^{2}_{u}(X,Y)$
$\displaystyle:=\frac{2}{n_{1}(n_{1}-1)}\sum_{i=1}^{n_{1}-1}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j})+\frac{2}{n_{2}(n_{2}-1)}\sum_{i=1}^{n_{2}-1}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j})$
$\displaystyle-\frac{2}{n_{1}n_{2}}\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}k(x_{i},y_{j}).$

Let $m_{1},m_{2}$ be any two positive integers such that $m_{1}\leq
n_{1},m_{2}\leq n_{2}$. Denote $c_{1}=\frac{2}{n_{1}(n_{1}-1)}$. Subsequently,

$\displaystyle\frac{2}{n_{1}(n_{1}-1)}\sum_{i=1}^{n_{1}-1}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j})$
$\displaystyle=c_{1}\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j})+c_{1}\sum_{i=m_{1}+1}^{n_{1}-1}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j})$
$\displaystyle=c_{1}\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k(x_{i},x_{j})+c_{1}\sum_{i=1}^{m_{1}}\sum_{j=m_{1}+1}^{n_{1}}k(x_{i},x_{j})$
$\displaystyle+c_{1}\sum_{i=m_{1}+1}^{n_{1}-1}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j}).$

Also, denote $c_{2}=\frac{2}{n_{2}(n_{2}-1)}$, it follows

$\displaystyle\frac{2}{n_{2}(n_{2}-1)}\sum_{i=1}^{n_{2}-1}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j})$
$\displaystyle=c_{2}\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j})+c_{2}\sum_{i=m_{2}+1}^{n_{2}-1}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j})$
$\displaystyle=c_{2}\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k(y_{i},y_{j})+c_{2}\sum_{i=1}^{m_{2}}\sum_{j=m_{2}+1}^{n_{2}}k(y_{i},y_{j})$
$\displaystyle+c_{2}\sum_{i=m_{2}+1}^{n_{2}-1}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j}).$

Denote $c_{3}=\frac{2}{n_{1}n_{2}}$, then

$\displaystyle\frac{2}{n_{1}n_{2}}\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}k(x_{i},y_{j})$
$\displaystyle=c_{3}\sum_{i=1}^{m_{1}}\sum_{j=1}^{n_{2}}k(x_{i},y_{j})+c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=1}^{n_{2}}k(x_{i},y_{j})$
$\displaystyle=c_{3}\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j})+c_{3}\sum_{i=1}^{m_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j})$
$\displaystyle+c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j})+c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j}).$

Combine the above together, ${\text{MMD}}^{2}_{u}(X,Y)$ can be rewritten as

$\displaystyle{\text{MMD}}^{2}_{u}(X,Y)$
$\displaystyle=c_{1}\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k(x_{i},x_{j})+c_{1}\sum_{i=1}^{m_{1}}\sum_{j=m_{1}+1}^{n_{1}}k(x_{i},x_{j})$
$\displaystyle+c_{1}\sum_{i=m_{1}+1}^{n_{1}-1}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j})+c_{2}\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k(y_{i},y_{j})$
$\displaystyle+c_{2}\sum_{i=1}^{m_{2}}\sum_{j=m_{2}+1}^{n_{2}}k(y_{i},y_{j})+c_{2}\sum_{i=m_{2}+1}^{n_{2}-1}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j})$
$\displaystyle-
c_{3}\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j})-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j})$
$\displaystyle-
c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j})-c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j}).$

By rearranging the above equation, we conclude

$\displaystyle{\text{MMD}}^{2}_{u}(X,Y)$
$\displaystyle=c_{1}\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k(x_{i},x_{j})+c_{2}\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j})$
$\displaystyle+c_{1}\sum_{i=m_{1}+1}^{n_{1}-1}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j})+c_{2}\sum_{i=m_{2}+1}^{n_{2}-1}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j})$
$\displaystyle+c_{1}\sum_{i=1}^{m_{1}}\sum_{j=m_{1}+1}^{n_{1}}k(x_{i},x_{j})-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j})$
$\displaystyle-
c_{2}\sum_{i=1}^{m_{2}}\sum_{j=m_{2}+1}^{n_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j}),$

which proves our result. ∎

This lemma divides the ${\text{MMD}}^{2}_{u}(X,Y)$ into four parts: the first
part $A_{1}$ includes the unobserved samples in $X$ and $Y$ only; on the
contrary, the second part $A_{2}$ includes observed samples in $X$ and $Y$
only; the third and the fourth terms $A_{3}$ and $A_{4}$ involve mixed terms
where one group (either $X$ or $Y$) contains only unobserved samples while the
corresponding samples in the other group are all observed.

In order to bound ${\text{MMD}}^{2}_{u}(X,Y)$, we propose to bound the four
terms $A_{1},A_{2},A_{3}$ and $A_{4}$ separately. Given that $A_{2}$ includes
the observed samples only, we focus on bounding $A_{1},A_{2}$ and $A_{4}$. For
the first term $A_{1}$, where all samples are unobserved, it can be seen that
$1\geq k(x,y)>0$ for any $x,y\in\mathbb{R}$. Hence, we have
${m_{1}(m_{1}-1)}/{2}\geq\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k(x_{i},x_{j})>0$,
${m_{2}(m_{2}-1)}/{2}\geq\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k(y_{i},y_{j})>0$,
and $m_{1}m_{2}\geq\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j})>0$,
which then according to the definition of $A_{1}$ in Lemma 4, follows

$\displaystyle\frac{m_{1}(m_{1}-1)}{n_{1}(n_{1}-1)}+\frac{m_{2}(m_{2}-1)}{n_{2}(n_{2}-1)}>A_{1}>-\frac{2}{n_{1}n_{2}}m_{1}m_{2}.$
(4)

The primary challenge of bounding ${\text{MMD}}^{2}_{u}(X,Y)$ lies in
providing tight bounds for terms $A_{3}$ and $A_{4}$, where both observed and
unobserved samples are presented. We focus on the term $A_{3}$, with the
understanding that a similar approach can be applied to the term $A_{4}$. To
start, we notice that the term $A_{3}$ can be rewritten as

$\displaystyle A_{3}$
$\displaystyle=\sum_{i=1}^{m_{1}}\left(c_{1}\sum_{j=m_{1}+1}^{n_{1}}k(x_{i},x_{j})-c_{3}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j})\right),$
(5)

where $c_{1}=\frac{2}{n_{1}(n_{1}-1)}$ and $c_{3}=\frac{2}{n_{1}n_{2}}$. Here,
we introduce $T_{1}(z)$, defined as

$\displaystyle
T_{1}(z)=c_{1}\sum_{j=m_{1}+1}^{n_{1}}k(z,x_{j})-c_{3}\sum_{j=m_{2}+1}^{n_{2}}k(z,y_{j}).$

Thus, $A_{3}$ can be further expressed as:

$\displaystyle A_{3}=\sum_{i=1}^{m_{1}}T_{1}(x_{i}),$ (6)

i.e. $A_{3}$ is the summation of the function values $T_{1}(z)$ applied to all
unobserved samples in $X$. This expression implies that in order to bound
$A_{3}$, it is sufficient to study bounds of function $T_{1}(z)$. We proceed
by examining $T_{1}(z)$’s behavior, using the following lemma:

###  A.2 Proof of Lemma 1

######  Lemma 1.

Let $x_{1},\ldots,x_{\ell_{1}}$ and $y_{1},\ldots,y_{\ell_{2}}$ be univariate
real values, that are observed. Suppose $a_{1},\ldots,a_{\ell_{1}}$,
$b_{1},\ldots,b_{\ell_{2}}$ and $\beta$ are positive constants. Define

$\displaystyle
T(z)=\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta|x_{i}-z|)-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta|y_{i}-z|)$

as a function of $z\in\mathbb{R}$. Subsequently, for any given
$z_{0}\in\mathbb{R}$,

$\displaystyle
T(z_{0})\geq\min\\{0,T(x_{1}),\ldots,T(x_{\ell_{1}}),T(y_{1}),\ldots,T(y_{\ell_{2}})\\}.$
(7)

On the other hand,

$\displaystyle
T(z_{0})\leq\max\\{0,T(x_{1}),\ldots,T(x_{\ell_{1}}),T(y_{1}),\ldots,T(y_{\ell_{2}})\\}.$
(8)

###### Proof.

We will first prove that inequality (7) holds. Let
$z_{1}=\min\\{x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}\\}$. When
$z<z_{1}$, it follows

$\displaystyle T(z)$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta(x_{i}-z))-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta(y_{i}-z))$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta(x_{i}-z_{1}+z_{1}-z))-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta(y_{i}-z_{1}+z_{1}-z))$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta(x_{i}-z_{1}))\exp(-\beta(z_{1}-z))-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta(y_{i}-z_{1}))\exp(-\beta(z_{1}-z))$
$\displaystyle=\exp(-\beta(z_{1}-z))\left\\{\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta(x_{i}-z_{1}))-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta(y_{i}-z_{1}))\right\\}$
$\displaystyle=\exp(-\beta(z_{1}-z))T(z_{1}).$

Notice that, since $\beta>0$ and $z<z_{1}$,

$\displaystyle 0<\exp(-\beta(z_{1}-z))<1.$

Hence, if $T(z_{1})\geq 0$,

$\displaystyle 0\leq\exp(-\beta(z_{1}-z))T(z_{1})\leq T(z_{1});$
$\displaystyle\Rightarrow$ $\displaystyle 0\leq T(z)$

if $T(z_{1})<0$,

$\displaystyle 0>\exp(-\beta(z_{1}-z))T(z_{1})>T(z_{1});$
$\displaystyle\Rightarrow$ $\displaystyle T(z)>T(z_{1}).$

Thus, for any given $z_{0}\in\mathbb{R}$, if $z_{0}\in(-\infty,z_{1})$, we
have

$\displaystyle T(z_{0})\geq\min\\{0,T(z_{1})\\},$

which proves (7) when $z_{0}\in(-\infty,z_{1})$.

Similarly, let
$z_{2}=\max\\{x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}\\}$. When
$z>z_{2}$, it follows

$\displaystyle T(z)$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta(z-x_{i}))-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta(z-y_{i}))$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta(z-z_{2}+z_{2}-x_{i}))-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta(z-z_{2}+z_{2}-y_{i}))$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta(z_{2}-x_{i}))\exp(-\beta(z-z_{2}))-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta(z_{2}-y_{i}))\exp(-\beta(z-z_{2}))$
$\displaystyle=\exp(-\beta(z-z_{2}))\left\\{\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta(z_{2}-x_{i}))-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta(z_{2}-y_{i}))\right\\}$
$\displaystyle=\exp(-\beta(z-z_{2}))T(z_{2}).$

Notice that, since $\beta>0$ and $z>z_{2}$,

$\displaystyle 0<\exp(-\beta(z-z_{2}))<1.$

Hence, if $T(z_{2})\geq 0$,

$\displaystyle 0\leq\exp(-\beta(z-z_{2}))T(z_{2})\leq T(z_{2});$
$\displaystyle\Rightarrow$ $\displaystyle 0\leq T(z)$

if $T(z_{1})<0$,

$\displaystyle 0>\exp(-\beta(z-z_{2}))T(z_{2})>T(z_{2});$
$\displaystyle\Rightarrow$ $\displaystyle T(z)>T(z_{2}).$

Thus, for any given $z_{0}\in\mathbb{R}$, if $z_{0}\in(z_{2},\infty)$, we have

$\displaystyle T(z_{0})>\min\\{0,T(z_{2})\\},$

which proves (7) when $z_{0}\in(z_{2},\infty)$.

Suppose $z_{1}<z_{0}<z_{2}$ and
$z_{0}\notin\\{x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}\\}$. Then,
there must be at least one number in
$\\{x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}\\}$ smaller than $z$,
and at least one number larger than $z$. Subsequently, denote $z_{3}$ as the
maximum number in $\\{x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}\\}$
smaller than $z_{0}$; denote $z_{4}$ as the minimum number in
$\\{x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}\\}$ larger than
$z_{0}$.

Notice that

$\displaystyle T(z)$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta|x_{i}-z|)-\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta|y_{i}-z|)$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}\leq
z_{3})\exp(-\beta|x_{i}-z|)+\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}>z_{3})\exp(-\beta|x_{i}-z|)$
$\displaystyle-\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}\leq
z_{3})\exp(-\beta|y_{i}-z|)-\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}>z_{3})\exp(-\beta|y_{i}-z|).$

When $z_{3}<z<z_{4}$, according to the definition of $z_{3}$ and $z_{4}$, for
any $x_{i},y_{i}\leq z_{3}$, it follows $z>x_{i},y_{i}$. On the other hand,
for for any $x_{i},y_{i}>z_{3}$, it follows $x_{i},y_{i}\geq z_{4}$, which
gives $z<x_{i},y_{i}$.

Hence, when $z_{3}<z<z_{4}$,

$\displaystyle T(z)$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}\leq
z_{3})\exp(-\beta(z-x_{i}))+\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}>z_{3})\exp(-\beta(x_{i}-z))$
$\displaystyle-\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}\leq
z_{3})\exp(-\beta(z-y_{i}))-\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}>z_{3})\exp(-\beta(y_{i}-z)).$

Further,

$\displaystyle T(z)$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}\leq
z_{3})\exp(-\beta(z-z_{3}+z_{3}-x_{i}))+\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}>z_{3})\exp(-\beta(x_{i}-z_{4}+z_{4}-z))$
$\displaystyle-\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}\leq
z_{3})\exp(-\beta(z-z_{3}+z_{3}-y_{i}))-\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}>z_{3})\exp(-\beta(y_{i}-z_{4}+z_{4}-z)).$
$\displaystyle=\exp(-\beta(z-z_{3}))\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}\leq
z_{3})\exp(-\beta(z_{3}-x_{i}))$
$\displaystyle+\exp(-\beta(z_{4}-z))\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}>z_{3})\exp(-\beta(x_{i}-z_{4}))$
$\displaystyle-\exp(-\beta(z-z_{3}))\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}\leq
z_{3})\exp(-\beta(z_{3}-y_{i}))$
$\displaystyle-\exp(-\beta(z_{4}-z))\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}>z_{3})\exp(-\beta(y_{i}-z_{4})).$
$\displaystyle=\exp(-\beta(z-z_{3}))\left\\{\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}\leq
z_{3})\exp(-\beta(z_{3}-x_{i}))-\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}\leq
z_{3})\exp(-\beta(z_{3}-y_{i}))\right\\}$
$\displaystyle+\exp(-\beta(z_{4}-z))\left\\{\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}>z_{3})\exp(-\beta(x_{i}-z_{4}))-\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}>z_{3})\exp(-\beta(y_{i}-z_{4}))\right\\}.$

For notation ease, let us denote

$\displaystyle A:=\left\\{\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}\leq
z_{3})\exp(-\beta(z_{3}-x_{i}))-\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}\leq
z_{3})\exp(-\beta(z_{3}-y_{i}))\right\\},$ $\displaystyle
B:=\left\\{\sum_{i=1}^{\ell_{1}}a_{i}\text{I}(x_{i}>z_{3})\exp(-\beta(x_{i}-z_{4}))-\sum_{i=1}^{\ell_{2}}b_{i}\text{I}(y_{i}>z_{3})\exp(-\beta(y_{i}-z_{4}))\right\\}.$

Thus,

$\displaystyle T(z)=A\exp(-\beta(z-z_{3}))+B\exp(-\beta(z_{4}-z)).$

Notice that if $B=0$, for any $z_{3}<z<z_{4}$, $T(z)=A\exp(-\beta(z-z_{3}))$
is a monotonic increasing function in $(z_{3},z_{4})$ when $A<0$, or a
monotonic decreasing function in $(z_{3},z_{4})$ when $A>0$, or a constant
function in $(z_{3},z_{4})$ when $A=0$. Thus, we always have
$\min\\{T(z_{3}),T(z_{4})\\}\leq T(z_{0})$ for $z_{0}\in(z_{3},z_{4})$, which
proves inequality (7) directly.

If, however, $B\neq 0$, let $T(z)=0$, we have

$\displaystyle A\exp(-\beta(z-z_{3}))=-B\exp(-\beta(z_{4}-z))$
$\displaystyle\Rightarrow$
$\displaystyle-\frac{A}{B}=\frac{\exp(-\beta(z_{4}-z))}{\exp(-\beta(z-z_{3}))}.$

When $-\frac{A}{B}\leq 0$, since
$\frac{\exp(-\beta(z_{4}-z))}{\exp(-\beta(z-z_{3}))}>0$, then $T(z)$ cannot
take 0 in $(z_{3},z_{4})$. If however, $-\frac{A}{B}>0$, we further have

$\displaystyle\log\left(-\frac{A}{B}\right)=-\beta(z_{4}-z)-\\{-\beta(z-z_{3})\\}$
$\displaystyle\Rightarrow$ $\displaystyle\log\left(-\frac{A}{B}\right)=-\beta
z_{4}+\beta z+\beta z-\beta z_{3}$ $\displaystyle\Rightarrow$ $\displaystyle
z=\frac{1}{\beta}\log\left(-\frac{A}{B}\right)+z_{4}+z_{3}.$

That is, $T(z)=0$ when
$z=\frac{1}{\beta}\log\left(-\frac{A}{B}\right)+z_{4}+z_{3}\in(z_{3},z_{4})$.
Overall, we have shown $T(z)$ takes 0 at most once in $(z_{3},z_{4})$ when
$B\neq 0$. This result will be used subsequently for proving our final
conclusion.

Now, taking derivative of $T(z)$ of $z$,

$\displaystyle\frac{\partial T(z)}{\partial z}$ $\displaystyle=-\beta
A\exp(-\beta(z-z_{3}))+\beta B\exp(-\beta(z_{4}-z)).$

Further, taking derivative of $\frac{\partial T(z)}{\partial z}$ of $z$,

$\displaystyle\frac{\partial^{2}T(z)}{\partial^{2}z}$
$\displaystyle=\beta^{2}A\exp(-\beta(z-z_{3}))+\beta^{2}B\exp(-\beta(z_{4}-z))$
$\displaystyle=\beta^{2}T(z).$

We are going to prove $T(z_{0})\geq\min\\{0,T{(z_{3})},T{(z_{4})}\\}$ when
$B\neq 0$ using contradiction. Let us assume that

$\displaystyle T(z_{0})<\min\\{0,T{(z_{3})},T{(z_{4})}\\}.$ (9)

Then, since $\beta>0$, it must have

$\displaystyle\left.\frac{\partial^{2}T(z)}{\partial^{2}z}\right|_{z_{0}}=\beta^{2}T(z_{0})<0.$

For $\left.\frac{\partial T(z)}{\partial z}\right|_{z_{0}}$, it is either
$(i):\left.\frac{\partial T(z)}{\partial z}\right|_{z_{0}}<0$ or
$(ii):\left.\frac{\partial T(z)}{\partial z}\right|_{z_{0}}\geq 0$. We are
going to discuss the two cases separately.

Suppose $(i):\left.\frac{\partial T(z)}{\partial z}\right|_{z_{0}}<0$. If for
any $z\in[z_{0},z_{4})$, $\frac{\partial^{2}T(z)}{\partial^{2}z}<0$. Then,
$\frac{\partial T(z)}{\partial z}$ is a monotonic decreasing function of $z$
in $[z_{0},z_{4})$. Hence, for any $z\in[z_{0},z_{4})$,

$\displaystyle\frac{\partial T(z)}{\partial z}<\left.\frac{\partial
T(z)}{\partial z}\right|_{z_{0}}<0$

Thus, $T(z)$ is a monotonic decreasing function in $z\in[z_{0},z_{4})$, giving
us

$\displaystyle T(z_{0})>T(z_{4}),$

which is contradicted to the assumption (9).

If, however, for any $z\in[z_{0},z_{4})$,
$\frac{\partial^{2}T(z)}{\partial^{2}z}<0$ does not holds, recall that we have
shown $T(z)$ take 0 at most once in $(z_{3},z_{4})$ when $B\neq 0$. This then
implies there is only one point $t\in[z_{0},z_{4})$ such that
$\left.\frac{\partial^{2}T(z)}{\partial^{2}z}\right|_{t}=0$. Since
$\frac{\partial^{2}T(z)}{\partial^{2}z}$ is a continuous function and
$\left.\frac{\partial^{2}T(z)}{\partial^{2}z}\right|_{z_{0}}<0$, we have

$\displaystyle\frac{\partial^{2}T(z)}{\partial^{2}z}<0,~{}z\in[z_{0},t).$

Thus, $\frac{\partial T(z)}{\partial z}$ is a monotonic decreasing function in
$z\in[z_{0},t)$, giving us

$\displaystyle\frac{\partial T(z)}{\partial z}<\left.\frac{\partial
T(z)}{\partial z}\right|_{z_{0}}<0,~{}z\in[z_{0},t).$ (10)

However, notice that we also have

$\displaystyle\left.\frac{\partial^{2}T(z)}{\partial^{2}z}\right|_{t}=0>\left.\frac{\partial^{2}T(z)}{\partial^{2}z}\right|_{z_{0}},$
$\displaystyle\Rightarrow$ $\displaystyle T(t)=0>T(z_{0})$
$\displaystyle\Rightarrow$ $\displaystyle\exists
t^{\prime}\in(z_{0},t)~{}\text{such that}\left.\frac{\partial T(z)}{\partial
z}\right|_{t^{\prime}}=\frac{T(t)-T(z_{0})}{t-z_{0}}>0,~{}~{}~{}(\text{mean
value theorem})$

which contradicts (10).

Suppose $(ii):\left.\frac{\partial T(z)}{\partial z}\right|_{z_{0}}\geq 0$.
The prove method is similar as to when $(i)$ holds.

If for any $z\in(z_{3},z_{0}]$, $\frac{\partial^{2}T(z)}{\partial^{2}z}<0$.
Then, $\frac{\partial T(z)}{\partial z}$ is a monotonic decreasing function of
$z$ in $(z_{3},z_{0}]$. Hence, for any $z\in(z_{3},z_{0}]$,

$\displaystyle\frac{\partial T(z)}{\partial z}>\left.\frac{\partial
T(z)}{\partial z}\right|_{z_{0}}\geq 0.$

Thus, $T(z)$ is a monotonic increasing function in $z\in(z_{3},z_{0}]$, giving

$\displaystyle T(z_{3})<T(z_{0}),$

which contradicts to the assumption (9).

If, however, for any $z\in(z_{3},z_{0}]$,
$\frac{\partial^{2}T(z)}{\partial^{2}z}<0$ does not holds, recall that we have
shown $T(z)$ taken 0 at most once in $(z_{3},z_{4})$. This then implies there
is only one point $t\in(z_{3},z_{0}]$ such that
$\left.\frac{\partial^{2}T(z)}{\partial^{2}z}\right|_{t}=0$. Since
$\frac{\partial^{2}T(z)}{\partial^{2}z}$ is a continuous function and
$\left.\frac{\partial^{2}T(z)}{\partial^{2}z}\right|_{z_{0}}<0$, we have

$\displaystyle\frac{\partial^{2}T(z)}{\partial^{2}z}<0,~{}z\in(t,z_{0}].$

Thus, $\frac{\partial T(z)}{\partial z}$ is a monotonic decreasing function in
$z\in(t,z_{0}]$, giving us

$\displaystyle\frac{\partial T(z)}{\partial z}>\left.\frac{\partial
T(z)}{\partial z}\right|_{z_{0}}\geq 0,~{}z\in(t,z_{0}].$ (11)

However, notice that we also have

$\displaystyle\left.\frac{\partial^{2}T(z)}{\partial^{2}z}\right|_{t}=0>\left.\frac{\partial^{2}T(z)}{\partial^{2}z}\right|_{z_{0}},$
$\displaystyle\Rightarrow$ $\displaystyle T(t)=0>T(z_{0})$
$\displaystyle\Rightarrow$ $\displaystyle\exists
t^{\prime}\in(t,z_{0})~{}\text{such that}\left.\frac{\partial T(z)}{\partial
z}\right|_{t^{\prime}}=\frac{T(t)-T(z_{0})}{t-z_{0}}<0.~{}~{}~{}(\text{mean
value theorem})$

which contradicts (11). Hence, we finish our prove for inequality (7).

We now prove inequality (8). Notice that

$\displaystyle-T(z)=\sum_{i=1}^{\ell_{2}}b_{i}\exp(-\beta|y_{i}-z|)-\sum_{i=1}^{\ell_{1}}a_{i}\exp(-\beta|x_{i}-z|).$

Subsequently, using the result of inequality (7), we have

$\displaystyle-T(z_{0})\geq\min\\{0,-T(x_{1}),\ldots,-T(x_{\ell_{1}}),-T(y_{1}),\ldots,-T(y_{\ell_{2}})\\}$
$\displaystyle\Rightarrow$ $\displaystyle
T(z_{0})\leq-\min\\{0,-T(x_{1}),\ldots,-T(x_{\ell_{1}}),-T(y_{1}),\ldots,-T(y_{\ell_{2}})\\}$
$\displaystyle\Rightarrow$ $\displaystyle
T(z_{0})\leq\max\\{0,T(x_{1}),\ldots,T(x_{\ell_{1}}),T(y_{1}),\ldots,T(y_{\ell_{2}})\\},$

which proves inequality (8). ∎

This lemma provides a linear time algorithm for computing the bounds of
$T(z)$, be checking all the function values of $T(z)$ to the set
$\\{x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}\\}$ and comparing all
these function values with 0. Notice that $T_{1}(z)$ is a special form of
$T(z)$ by letting all $a_{1},\ldots,a_{\ell_{1}}$ in Lemma 1 equal to $c_{1}$
and all $b_{1},\ldots,b_{\ell_{2}}$ equal to $c_{2}$, where
$c_{1}=\frac{2}{n_{1}(n_{1}-1)}$, $c_{2}=\frac{2}{n_{2}(n_{2}-1)}$. Hence,
Lemma 1 allows us to compute the bounds of $T_{1}(z)$ efficiently, which then
gives the bounds of $A_{3}$ according to Equation (6). This result, combining
with Equation (4), allows us to conclude the bounds of MMD for univariate
samples, which is proved in Theorem 4.

###  A.3 Proof of Lemma 2

This section proves Lemma 2, which extends Lemma 1 into higher dimensions.

######  Lemma 2.

Let $x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}\in\mathbb{R}^{d}$ be
values that are fully observed. Suppose $a_{1},\ldots,a_{\ell_{1}}$,
$b_{1},\ldots,b_{\ell_{2}},\beta$ are positive constants. For
$z=(z(1),\ldots,z(d))\in\mathbb{R}^{d}$ with missing components, define

$\displaystyle T(\\{z(j)\,:\,j\in
U_{z}\\})=\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta\sum_{j\in
U_{z}}|x_{i}(j)-z(j)|\right)-\sum_{i=1}^{\ell_{2}}b_{i}\exp\left(-\beta\sum_{j\in
U_{z}}|y_{i}(j)-z(j)|\right)$

as a function of the unobserved components of $z$ and let

$\displaystyle\mathcal{X}=\\{T(\\{z(j)\,:\,j\in U_{z}\\})\
:\,z(i)\in\\{x_{1}(i),\ldots,x_{\ell_{1}}(i),y_{1}(i),\ldots,y_{\ell_{2}}(i)\\},i\in
U_{z}\\}.$

Then, for any possible imputation $z^{\ast}$ of $z$,

$\displaystyle T(\\{z^{\ast}(j)\,:\,j\in
U_{z}\\})\geq\min\\{0,\min\mathcal{X}\\},$ (12)

and

$\displaystyle T(\\{z^{\ast}(j)\,:\,j\in
U_{z}\\})\leq\max\\{0,\max\mathcal{X}\\}.$ (13)

###### Proof.

We will only prove inequality (12) and inequality (13) can be proved following
the same method.

Let us prove inequality (12) using mathematical induction. When $|U_{z}|=1$,
this lemma is exactly Lemma 1, which has already been shown. Suppose when
$|U_{z}|=l-1$, where $l$ is any integer such that $1\leq l$, inequality (12)
holds. We are going to show when $|U_{z}|=l$, inequality (12) still holds.

Without loss of generality, let us assume (after relabeling) that
$U_{z}=1,\ldots,l$. Subsequently, for any given
$(z^{*}_{1},\ldots,z^{*}_{l})\in\mathbb{R}^{l}$, it follows

$\displaystyle T(z^{*}_{1},\ldots,z^{*}_{l})$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta\sum_{j=1}^{l}|x_{i}(j)-z^{*}_{j}|\right)-\sum_{i=1}^{\ell_{2}}b_{i}\exp\left(-\beta\sum_{j=1}^{l}|y_{i}(j)-z^{*}_{j}|\right)$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta|x_{i}(1)-z^{*}_{1}|\right)\exp\left(-\beta\sum_{j=2}^{l}|x_{i}(j)-z^{*}_{j}|\right)$
$\displaystyle-\sum_{i=1}^{\ell_{2}}b_{i}\exp\left(-\beta|y_{i}(1)-z^{*}_{1}|\right)\exp\left(-\beta\sum_{j=2}^{l}|y_{i}(j)-z^{*}_{j}|\right).$

Denote

$\displaystyle a_{i}^{\prime}$
$\displaystyle=a_{i}\exp\left(-\beta|x_{i}(1)-z^{*}_{1}|\right),~{}i=1,\ldots,\ell_{1},$
$\displaystyle b_{i}^{\prime}$
$\displaystyle=b_{i}\exp\left(-\beta|y_{i}(1)-z^{*}_{1}|\right),~{}i=1,\ldots,\ell_{2}.$

Consider a new function

$\displaystyle
T^{\prime}(z_{2},\ldots,z_{l})=\sum_{i=1}^{\ell_{1}}a_{i}^{\prime}\exp\left(-\beta\sum_{j=2}^{l}|x_{i}(j)-z_{j}|\right)-\sum_{i=1}^{\ell_{2}}b_{i}^{\prime}\exp\left(-\beta\sum_{j=2}^{l}|y_{i}(j)-z_{j}|\right).$

Let

$\displaystyle\mathcal{X}^{\prime}=\\{T^{\prime}(z_{2},\ldots,z_{l})|z_{j}\in\\{x_{1}(j),\ldots,x_{\ell_{1}}(j),y_{1}(j),\ldots,y_{\ell_{2}}(j)\\},j=2,\ldots,l\\}.$

Then, using the assumption that when $|U_{z}|=l-1$, inequality (12) holds, it
follows that

$\displaystyle T^{\prime}(z^{*}_{2},\ldots,z^{*}_{l})$
$\displaystyle\geq\min\\{0,\min\mathcal{X}^{\prime}\\},$

Notice that

$\displaystyle\begin{split}T^{\prime}(z_{2},\ldots,z_{l})&=\sum_{i=1}^{\ell_{1}}a_{i}^{\prime}\exp\left(-\beta\sum_{j=2}^{l}|x_{i}(j)-z_{j}|\right)-\sum_{i=1}^{\ell_{2}}b_{i}^{\prime}\exp\left(-\beta\sum_{j=2}^{l}|y_{i}(j)-z_{j}|\right)\\\
&=\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta|x_{i}(1)-z^{*}_{1}|\right)\exp\left(-\beta\sum_{j=2}^{l}|x_{i}(j)-z_{j}|\right)\\\
&-\sum_{i=1}^{\ell_{2}}b_{i}\exp\left(-\beta|y_{i}(1)-z^{*}_{1}|\right)\exp\left(-\beta\sum_{j=2}^{l}|y_{i}(j)-z_{j}|\right)\\\
&=T(z^{*}_{1},z_{2},\ldots,z_{l}).\end{split}$ (14)

Hence,

$\displaystyle
T^{\prime}(z^{*}_{2},\ldots,z^{*}_{l})=T(z^{*}_{1},\ldots,z^{*}_{l}).$

We therefore have

$\displaystyle T(z^{*}_{1},\ldots,z^{*}_{l})$
$\displaystyle\geq\min\\{0,\min\chi^{\prime}\\}.$ (15)

Further, let us denote

$\displaystyle(z_{2}^{\prime},\ldots,z_{l}^{\prime})=\operatorname*{arg\,min}_{z_{i}\in\\{x_{1}(i),\ldots,x_{\ell_{1}}(i),y_{1}(i),\ldots,y_{\ell_{2}}(i)\\},i=2,\ldots,l}T^{\prime}(z_{2},\ldots,z_{l}).$
(16)

That is,

$\displaystyle
T^{\prime}(z_{2}^{\prime},\ldots,z_{l}^{\prime})=\min\mathcal{X}^{\prime}.$

Applying (14), it follows

$\displaystyle
T(z^{*}_{1},z_{2}^{\prime},\ldots,z_{l}^{\prime})=T^{\prime}(z_{2}^{\prime},\ldots,z_{l}^{\prime})=\min\mathcal{X}^{\prime}.$
(17)

Denote

$\displaystyle a_{i}^{\prime\prime}$
$\displaystyle=a_{i}\exp\left(-\beta\sum_{j=2}^{l}|x_{i}(j)-z_{j}^{\prime}|\right),~{}i=1,\ldots,\ell_{1},$
$\displaystyle b_{i}^{\prime\prime}$
$\displaystyle=b_{i}\exp\left(-\beta\sum_{j=2}^{l}|y_{i}(j)-z_{j}^{\prime}|\right)~{}i=1,\ldots,\ell_{2}.$

Consider a new function

$\displaystyle
T^{\prime\prime}(z_{1})=\sum_{i=1}^{\ell_{1}}a_{i}^{\prime\prime}\exp\left(-\beta|x_{i}(1)-z_{1}|\right)-\sum_{i=1}^{\ell_{1}}b_{i}^{\prime\prime}\exp\left(-\beta|y_{i}(1)-z_{1}|\right).$

Let

$\displaystyle\mathcal{X}^{\prime\prime}=\\{T^{\prime\prime}(z_{1})|z_{1}\in\\{x_{1}(1),\ldots,x_{\ell_{1}}(1),y_{1}(1),\ldots,y_{\ell_{2}}(1)\\}\\}.$

Then, using the result that when $|U_{z}|=1$, (12) holds, it follows

$\displaystyle T^{\prime\prime}(z^{*}_{1})$
$\displaystyle\geq\min\\{0,\min\mathcal{X}^{\prime\prime}\\}.$ (18)

Notice that

$\displaystyle T^{\prime\prime}(z_{1})$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}^{\prime\prime}\exp\left(-\beta|x_{i}(1)-z_{1}|\right)-\sum_{i=1}^{\ell_{1}}b_{i}^{\prime\prime}\exp\left(-\beta|y_{i}(1)-z_{1}|\right)$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta\sum_{j=2}^{l}|x_{i}(j)-z_{j}^{\prime}|\right)\exp\left(-\beta|x_{i}(1)-z_{1}|\right)$
$\displaystyle-\sum_{i=1}^{\ell_{1}}b_{i}\exp\left(-\beta\sum_{j=2}^{l}|y_{i}(j)-z_{j}^{\prime}|\right)\exp\left(-\beta|y_{i}(1)-z_{1}|\right)$
$\displaystyle=T(z_{1},z_{2}^{\prime},\ldots,z_{l}^{\prime}).$

According to (17),

$\displaystyle T^{\prime\prime}(z^{*}_{1})=\min\mathcal{X}^{\prime}.$

Notice that since (18) holds, we have

$\displaystyle\min\mathcal{X}^{\prime}\geq\min\\{0,\min\mathcal{X}^{\prime\prime}\\},$

put which back into (15), we further have

$\displaystyle T(z^{*}_{1},\ldots,z^{*}_{l})$
$\displaystyle\geq\min\\{0,\min\mathcal{X}^{\prime}\\}\geq\min\\{0,\min\mathcal{X}^{\prime\prime}\\}.$

Notice that

$\displaystyle\min\mathcal{X}^{\prime\prime}$
$\displaystyle=\min\\{T^{\prime\prime}(z_{1})|z_{1}\in\\{x_{1}(1),\ldots,x_{\ell_{1}}(1),y_{1}(1),\ldots,y_{\ell_{2}}(1)\\}\\}$
$\displaystyle=\min\\{T(z_{1},z_{2}^{\prime},\ldots,z_{l}^{\prime})|z_{1}\in\\{x_{1}(1),\ldots,x_{\ell_{1}}(1),y_{1}(1),\ldots,y_{\ell_{2}}(1)\\}\\},$

where

$\displaystyle
z_{i}^{\prime}\in\\{x_{1}(i),\ldots,x_{\ell_{1}}(i),y_{1}(i),\ldots,y_{\ell_{2}}(i)\\},i=2,\ldots,l$

according to its definition in (16). Hence,

$\displaystyle\min\mathcal{X}^{\prime\prime}\geq\min\chi\Rightarrow
T(z^{*}_{1},\ldots,z^{*}_{l})\geq\min\\{0,\min\mathcal{X}\\},$

which proves inequality (12). ∎

###  A.4 Proof of Lemma 3

In Lemma 1, it is proved that in order to compute the bounds of
$T(\\{z(j)\,:\,j\in U_{z}\\})$ for $z\in\mathbb{R}^{d}$, we only need to check
the imputations of $z$ where its missing components are imputed using the
components of $x_{1},\ldots,x_{\ell_{1}},y_{1},\ldots,y_{\ell_{2}}$. However,
computing $T(\\{z(j)\,:\,j\in U_{z}\\})$ for all possible imputations using
Lemma 2 is $(\ell_{1}+\ell_{2})^{|U_{x_{i}}|}$, which is exponential in the
number of unobserved components of $z$ and impractical to compute. To address
this practical computational challenge, we introduce the following lemma:

######  Lemma 3.

Following the notation and definitions in Lemma 2, denote

$\displaystyle\tilde{x}_{i}(j)$
$\displaystyle=\max\\{|x_{i}(j)-x_{1}(j)|,\ldots,|x_{i}(j)-x_{\ell_{1}}(j)|,|x_{i}(j)-y_{1}(j)|,\ldots,|x_{i}(j)-y_{\ell_{2}}(j)|\\}$

for any $i\in\\{1,\ldots,\ell_{1}\\},j\in U_{z}$; denote

$\displaystyle\tilde{y}_{i}(j)$
$\displaystyle=\max\\{|y_{i}(j)-\ell_{1}(j)|,\ldots,|y_{i}(j)-x_{n_{1}}(j)|,|y_{i}(j)-y_{1}(j)|,\ldots,|y_{i}(j)-y_{\ell_{2}}(j)|\\}$

for any $i\in\\{1,\ldots,\ell_{2}\\},j\in U_{z}$. Subsequently,

$\displaystyle\min\mathcal{X}\geq\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta\sum_{j\in
U_{z}}\tilde{x}_{i}(j)\right)-\sum_{i=1}^{\ell_{2}}b_{i},$ (19)
$\displaystyle\max\mathcal{X}\leq\sum_{i=1}^{\ell_{1}}a_{i}-\sum_{i=1}^{\ell_{2}}b_{i}\exp\left(-\beta\sum_{j\in
U_{z}}\tilde{y}_{i}(j)\right).$ (20)

###### Proof.

We will only prove inequality (19) and inequality (20) can be proved
similarly.

Without loss of generality, let us assume (after relabeling)
$U_{z}=1,\ldots,l$, where $l=|U_{z}|$. Subsequently, let

$\displaystyle(z^{*}_{1},\ldots,z^{*}_{l})=\operatorname*{arg\,min}_{\\{z_{i}\in\\{x_{1}(i),\ldots,x_{\ell_{1}}(i),y_{1}(i),\ldots,y_{\ell_{2}}(i)\\},i=1,\ldots,l\\}}T(z_{1},\ldots,z_{l}).$

That is,

$\displaystyle T(z^{*}_{1},\ldots,z^{*}_{l})=\min\mathcal{X}.$

Notice that for any $j\in\\{1,\ldots,l\\}$,

$\displaystyle
z^{*}_{j}\in\\{x_{1}(j),\ldots,x_{\ell_{1}}(j),y_{1}(j),\ldots,y_{\ell_{2}}(j)\\}$

according to its definition. Hence, for any $i\in\\{1,\ldots,\ell_{1}\\}$,

$\displaystyle|z^{*}_{j}-x_{i}(j)|\leq\max\\{|x_{i}(j)-x_{1}(j)|,\ldots,|x_{i}(j)-x_{\ell_{1}}(j)|=\tilde{x}_{i}(j),$

following which

$\displaystyle\sum_{j=1}^{l}\tilde{x}_{i}(j)\geq\sum_{j=1}^{l}|x_{i}(j)-z^{*}_{j}|,~{}i=\\{1,\ldots,\ell_{1}\\}$
$\displaystyle\Rightarrow$
$\displaystyle-\beta\sum_{j=1}^{l}\tilde{x}_{i}(j)\leq-\beta\sum_{j=1}^{l}|x_{i}(j)-z^{*}_{j}|,~{}i=\\{1,\ldots,\ell_{1}\\}.$

Subsequently,

$\displaystyle T(z^{*}_{1},\ldots,z^{*}_{l})$
$\displaystyle=\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta\sum_{j=1}^{l}|x_{i}(j)-z^{*}_{j}|\right)-\sum_{i=1}^{\ell_{2}}b_{i}\exp\left(-\beta\sum_{j=1}^{l}|y_{i}(j)-z^{*}_{j}|\right)$
$\displaystyle\geq\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta\sum_{j=1}^{l}\tilde{x}_{i}(j)\right)-\sum_{i=1}^{\ell_{2}}b_{i}\exp\left(-\beta\sum_{j=1}^{l}|y_{i}(j)-z^{*}_{j}|\right)$
$\displaystyle\geq\sum_{i=1}^{\ell_{1}}a_{i}\exp\left(-\beta\sum_{j=1}^{l}\tilde{x}_{i}(j)\right)-\sum_{i=1}^{\ell_{2}}b_{i},$

which proves inequality (19). ∎

###  A.5 Proof of Theorem 1 when $d=1$

We can now prove our main result Theorem 1. To do so, we divide the proof of
Theorem 1 into univariate case ($d=1$) and multivariate case ($d>1$), which
are provided in Theorem 4 and Theorem 5, respectively.

######  Theorem 4.

Suppose $X=\\{x_{1},\ldots,x_{n_{1}}\\}$ and $Y=\\{y_{1},\ldots,y_{n_{2}}\\}$
are univariate real values. Assume $x_{1},\ldots,x_{m_{1}}$,
$y_{1},\ldots,y_{m_{2}}$ are unobserved. Let $k$ denote the Laplacian kernel
and define

$\displaystyle T_{1}(z)$
$\displaystyle=c_{1}\sum_{j=m_{1}+1}^{n_{1}}k(z,x_{j})-c_{3}\sum_{j=m_{2}+1}^{n_{2}}k(z,y_{j}),$
$\displaystyle T_{2}(z)$
$\displaystyle=c_{2}\sum_{j=m_{2}+1}^{n_{2}}k(z,x_{j})-c_{3}\sum_{j=m_{1}+1}^{n_{1}}k(z,y_{j}),$

where $c_{1}=\frac{2}{n_{1}(n_{1}-1)},c_{2}=\frac{2}{n_{2}(n_{2}-1)}$ and
$c_{3}=\frac{2}{n_{1}n_{2}}$. Further, let

$\displaystyle
S_{1}:=\\{0,T_{1}(x_{m_{1}+1}),\ldots,T_{1}(x_{n_{1}}),T_{1}(y_{m_{2}+1}),\ldots,T_{1}(y_{n_{2}})\\},$
$\displaystyle
S_{2}:=\\{0,T_{2}(x_{m_{2}+1}),\ldots,T_{2}(x_{n_{2}}),T_{2}(y_{m_{1}+1}),\ldots,T_{2}(y_{n_{1}})\\}.$

Then, the ${\text{MMD}}^{2}_{u}(X,Y)$ using Laplacian kernel $k$ is bounded as
follows

$\displaystyle\frac{m_{1}(m_{1}-1)}{n_{1}(n_{1}-1)}+\frac{m_{2}(m_{2}-1)}{n_{2}(n_{2}-1)}+m_{1}\max
S_{1}+m_{2}\max S_{2}+A_{2}>{\text{MMD}}^{2}_{u}(X,Y),$
$\displaystyle{\text{MMD}}^{2}_{u}(X,Y)>A_{2}+m_{1}\min S_{1}+m_{2}\min
S_{2}-\frac{2}{n_{1}n_{2}}m_{1}m_{2},$

where $A_{2}$ is defined in lemma 4.

###### Proof.

To start, in Lemma 4, it is shown that

$\displaystyle{\text{MMD}}^{2}_{u}(X,Y)=A_{1}+A_{2}+A_{3}+A_{4},$

where

$\displaystyle A_{1}$
$\displaystyle=c_{1}\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k(x_{i},x_{j})+c_{2}\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j}),$
$\displaystyle A_{2}$
$\displaystyle=c_{1}\sum_{i=m_{1}+1}^{n_{1}-1}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j})+c_{2}\sum_{i=m_{2}+1}^{n_{2}-1}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j}),$
$\displaystyle A_{3}$
$\displaystyle=c_{1}\sum_{i=1}^{m_{1}}\sum_{j=m_{1}+1}^{n_{1}}k(x_{i},x_{j})-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j}),$
$\displaystyle A_{4}$
$\displaystyle=c_{2}\sum_{i=1}^{m_{2}}\sum_{j=m_{2}+1}^{n_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j}),$

and $c_{1}=\frac{2}{n_{1}(n_{1}-1)},c_{2}=\frac{2}{n_{2}(n_{2}-1)}$ and
$c_{3}=\frac{2}{n_{1}n_{2}}$.

Following Equation (4), it is provided

$\displaystyle\frac{m_{1}(m_{1}-1)}{n_{1}(n_{1}-1)}+\frac{m_{2}(m_{2}-1)}{n_{2}(n_{2}-1)}>A_{1}>-\frac{2}{n_{1}n_{2}}m_{1}m_{2}.$

Further, using Lemma 1, it follows

$\displaystyle m_{1}\max S_{1}\geq A_{3}\geq m_{1}\min S_{1},$

and

$\displaystyle m_{2}\max S_{2}\geq A_{4}\geq m_{2}\min S_{2},$

which concludes our proof. ∎

In order to compute the bounds of ${\text{MMD}}^{2}_{u}(X,Y)$ for univariate
samples with missing data using Theorem 4, we need to compute $\min S_{1},\max
S_{1},\min S_{2}$ and $\max S_{2}$. Since $|S_{1}|=n_{1}-m_{1}$ and
$|S_{2}|=n_{2}-m_{2}$. The computation complexity is of order
$O(n_{1}+n_{2})$.

###  A.6 Proof of Theorem 1 when $d>1$

In Theorem 4, the result of Theorem 1 when $d=1$ is provided. This section
proves the result of Theorem 1 when $d>1$. The final conclusion of this
section is provided in Theorem 5. To obtain this result, we first introduce
the following definition of the incomplete Laplacian kernel, which can be
computed between two incomplete samples.

######  Definition 3.

For any $x,y\in\mathbb{R}^{d}$, let $U_{x,y}\subset\\{1,\ldots,d\\}$ be the
index of dimensions for which either $x$ or $y$ have components that are not
observed. Let $k$ denote the Laplacian kernel with parameter $\beta$. Then,
the incomplete Laplacian kernel is defined as

$\displaystyle
k^{*}(x,y)=\exp\left(-\beta\sum_{i\in\\{1,\ldots,d\\}\setminus\left(U_{x}\cup
U_{y}\right)}|x(i)-y(i)|\right).$

Subsequently, we have the following result.

######  Lemma 5.

Suppose $x$ and $y$ are two samples of $\mathbb{R}^{d}$ real values and assume
not all dimensions of values of $x$ and $y$ are observed. Denote
$[d]=\\{1,\ldots,d\\}$. Let $U_{x}$ be a set includes all unobserved
dimensions in $x$ and $U_{y}$ be a set includes all unobserved dimensions in
$y$. Let $k^{*}$ denote the incomplete Laplacian kernel. Subsequently, it
follows

$\displaystyle k^{*}(x,y)\geq k(x,y)\geq 0.$

###### Proof.

According to the definition of the Laplacian kernel, we have

$\displaystyle k(x,y)=\exp\left(-\beta\sum_{l=1}^{d}|x(l)-y(l)|\right)\geq 0.$

Also, notice that

$\displaystyle\sum_{l\in[d]\setminus\left(U_{x}\cup
U_{y}\right)}|x(l)-y(l)|\leq\sum_{l=1}^{d}|x(l)-y(l)|.$

Then, since $\beta>0$,

$\displaystyle-\beta\sum_{i\in[d]\setminus\left(U_{x}\cup
U_{y}\right)}|x(i)-y(i)|\geq-\beta\sum_{l=1}^{d}|x(l)-y(l)|.$

According to the definition of incomplete Laplacian kernel, we therefore have

$\displaystyle k^{*}(x,y)$
$\displaystyle=\exp\left(-\beta\sum_{l\in[d]\setminus\left(U_{x}\cup
U_{y}\right)}|x(l)-y(l)|\right)$
$\displaystyle\geq\exp\left(-\beta\sum_{l=1}^{d}|x(l)-y(l)|\right)$
$\displaystyle=k(x,y),$

which finishes our proof. ∎

Recall that in Lemma 4, it is shown that ${\text{MMD}}^{2}_{u}(X,Y)$ can be
decomposed into four parts $A_{1},A_{2},A_{3}$ and $A_{4}$ with part $A_{1}$
includes incomplete samples only, $A_{2}$ includes complete samples only, and
$A_{3}$, $A_{4}$ mixed with both complete and incomplete samples.

The Lemma 5 allows us to bound term $A_{1}$ by bounding any two incompletely
observed samples separately. Let $U_{x_{i}}$ be a set including all unobserved
dimensions in $x_{i}$, where $1\leq i\leq m_{1}$. Similarly, let $U_{y_{i}}$
be a set including all unobserved dimensions in $y_{i}$, where $1\leq i\leq
m_{2}$. Subsequently, according to Lemma 5, it can be seen that
$\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k^{*}(x_{i},x_{j})\geq\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k(x_{i},x_{j})>0$,
$\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k^{*}(y_{i},y_{j})\geq\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k(y_{i},y_{j})>0$,
and
$\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k^{*}(x_{i},y_{j})\geq\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j})>0$.
Then, according to the definition of $A_{1}$ in Lemma 4, it follows

$\displaystyle
c_{1}{\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k^{*}(x_{i},x_{j})}+c_{2}{\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k^{*}(y_{i},y_{j})}>A_{1}>-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k^{*}(x_{i},y_{j}),$
(21)

where $c_{1}=\frac{2}{n_{1}(n_{1}-1)},c_{2}=\frac{2}{n_{2}(n_{2}-1)}$ and
$c_{3}=\frac{2}{n_{1}n_{2}}$. We now bound the terms $A_{3}$ and $A_{4}$.
Recall that the term $A_{3}$ can be rewritten as (5) in both univariate and
multivariate cases. For any $x_{i}$ where $1\leq i\leq m_{1}$, let $U_{x_{i}}$
be a set including all unobserved dimensions in $x_{i}$ and let us introduce
the function

$\displaystyle T_{1}(x_{i}(l);l\in
U_{x_{i}})=c_{1}\sum_{j=m_{1}+1}^{n_{1}}k(x_{i},x_{j})-c_{3}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j}),~{}i\in\\{1,\ldots,m_{1}\\}$
(22)

where $c_{1}=\frac{2}{n_{1}(n_{1}-1)}$ and $c_{3}=\frac{2}{n_{1}n_{2}}$. Then,
$A_{3}$ can be computed as

$\displaystyle A_{3}=\sum_{i=1}^{m_{1}}T_{1}(x_{i}(l);l\in U_{x_{i}}).$

Subsequently, by combining Lemma 2 and Lemma 3, we have the results of bounds
of ${\text{MMD}}^{2}_{u}(X,Y)$ when $d>1$.

######  Theorem 5.

Suppose $X=\\{x_{1},\ldots,x_{n_{1}}\\}$ and $Y=\\{y_{1},\ldots,y_{n_{2}}\\}$
are samples of $\mathbb{R}^{d}$ real values. Assume $x_{1},\ldots,x_{m_{1}}$,
$y_{1},\ldots,y_{m_{2}}$ are samples that are not observed completely. Let $k$
denote the Laplacian kernel and $k^{*}$ denote the incomplete Laplacian kernel
defined in Definition 3. Further, denote

$\displaystyle\tilde{x}_{i}(l)$
$\displaystyle=\max\\{|x_{i}(l)-x_{m_{1}+1}(l)|,\ldots,|x_{i}(l)-x_{n_{1}}(l)|,|x_{i}(l)-y_{m_{2}+1}(l)|,\ldots,|x_{i}(l)-y_{n_{2}}(l)|\\}$

for any $i\in\\{m_{1}+1,\ldots,n_{1}\\},l\in\\{1,\ldots,d\\}$; denote

$\displaystyle\tilde{y}_{i}(l)$
$\displaystyle=\max\\{|y_{i}(l)-x_{m_{1}+1}(l)|,\ldots,|y_{i}(l)-x_{n_{1}}(l)|,|y_{i}(l)-y_{m_{2}+1}(l)|,\ldots,|y_{i}(l)-y_{n_{2}}(l)|\\}$

for any $i\in\\{m_{2}+1,\ldots,n_{2}\\},l\in\\{1,\ldots,d\\}$. Let

$\displaystyle C_{1}$
$\displaystyle=\sum_{i=1}^{m_{1}}\max\left\\{0,c_{1}\sum_{j=m_{1}+1}^{n_{1}}k^{*}(x_{i},x_{j})-c_{3}\sum_{j=m_{2}+1}^{n_{2}}k^{*}(x_{i},y_{j})\exp\left(-\beta\sum_{l\in
U_{x_{i}}}\tilde{y}_{j}(l)\right)\right\\},$ $\displaystyle C_{2}$
$\displaystyle=\sum_{i=1}^{m_{1}}\min\left\\{0,c_{1}\sum_{j=m_{1}+1}^{n_{1}}k^{*}(x_{i},x_{j})\exp\left(-\beta\sum_{l\in
U_{x_{i}}}\tilde{x}_{j}(l)\right)-c_{3}\sum_{j=m_{2}+1}^{n_{2}}k^{*}(x_{i},y_{j})\right\\},$
$\displaystyle C_{3}$
$\displaystyle=\sum_{i=1}^{m_{2}}\max\left\\{0,c_{2}\sum_{j=m_{2}+1}^{n_{2}}k^{*}(y_{i},y_{j})-c_{3}\sum_{j=m_{1}+1}^{n_{1}}k^{*}(y_{i},x_{j})\exp\left(-\beta\sum_{l\in
U_{y_{i}}}\tilde{x}_{j}(l)\right)\right\\},$ $\displaystyle C_{4}$
$\displaystyle=\sum_{i=1}^{m_{2}}\min\left\\{0,c_{2}\sum_{j=m_{2}+1}^{n_{2}}k^{*}(y_{i},y_{j})\exp\left(-\beta\sum_{l\in
U_{y_{i}}}\tilde{y}_{j}(l)\right)-c_{3}\sum_{j=m_{1}+1}^{n_{1}}k^{*}(y_{i},x_{j})\right\\},$

where $c_{1}=\frac{2}{n_{1}(n_{1}-1)},c_{2}=\frac{2}{n_{2}(n_{2}-1)}$ and
$c_{3}=\frac{2}{n_{1}n_{2}}$. Subsequently, the ${\text{MMD}}^{2}_{u}(X,Y)$
using Laplacian kernel $k$ is bounded as follows

$\displaystyle
c_{1}{\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k^{*}(x_{i},x_{j})}+c_{2}{\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k^{*}(y_{i},y_{j})}+A_{2}+C_{1}+C_{3}>{\text{MMD}}^{2}_{u}(X,Y),$
$\displaystyle{\text{MMD}}^{2}_{u}(X,Y)>-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k^{*}(x_{i},y_{j})+A_{2}+C_{2}+C_{4},$

where $A_{2}$ is defined in Lemma 4.

###### Proof.

In Lemma 4, it is proved that

$\displaystyle{\text{MMD}}^{2}_{u}(X,Y)=A_{1}+A_{2}+A_{3}+A_{4},$

where

$\displaystyle A_{1}$
$\displaystyle=c_{1}\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k(x_{i},x_{j})+c_{2}\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j}),$
$\displaystyle A_{2}$
$\displaystyle=c_{1}\sum_{i=m_{1}+1}^{n_{1}-1}\sum_{j=i+1}^{n_{1}}k(x_{i},x_{j})+c_{2}\sum_{i=m_{2}+1}^{n_{2}-1}\sum_{j=i+1}^{n_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j}),$
$\displaystyle A_{3}$
$\displaystyle=c_{1}\sum_{i=1}^{m_{1}}\sum_{j=m_{1}+1}^{n_{1}}k(x_{i},x_{j})-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j}),$
$\displaystyle A_{4}$
$\displaystyle=c_{2}\sum_{i=1}^{m_{2}}\sum_{j=m_{2}+1}^{n_{2}}k(y_{i},y_{j})-c_{3}\sum_{i=m_{1}+1}^{n_{1}}\sum_{j=1}^{m_{2}}k(x_{i},y_{j}),$

and $c_{1}=\frac{2}{n_{1}(n_{1}-1)},c_{2}=\frac{2}{n_{2}(n_{2}-1)}$ and
$c_{3}=\frac{2}{n_{1}n_{2}}$.

According to inequality (21),

$\displaystyle
c_{1}{\sum_{i=1}^{m_{1}}\sum_{j=i+1}^{m_{1}}k^{*}(x_{i},x_{j})}+c_{2}{\sum_{i=1}^{m_{2}}\sum_{j=i+1}^{m_{2}}k^{*}(y_{i},y_{j})}>A_{1}>-c_{3}\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}k^{*}(x_{i},y_{j}).$
(23)

Notice that

$\displaystyle A_{3}$
$\displaystyle=\sum_{i=1}^{m_{1}}\left(c_{1}\sum_{j=m_{1}+1}^{n_{1}}k(x_{i},x_{j})-c_{3}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j})\right).$

Let us introduce the function

$\displaystyle T_{1}(x_{i}(l);l\in
U_{x_{i}})=c_{1}\sum_{j=m_{1}+1}^{n_{1}}k(x_{i},x_{j})-c_{3}\sum_{j=m_{2}+1}^{n_{2}}k(x_{i},y_{j}),~{}i\in\\{1,\ldots,m_{1}\\}.$

Subsequently, $A_{3}$ can be computed as

$\displaystyle A_{3}=\sum_{i=1}^{m_{1}}T_{1}(x_{i}(l);l\in U_{x_{i}}).$ (24)

Denote

$\displaystyle\mathcal{X}_{x_{i}}=\\{T_{1}(x_{i}(l);l\in
U_{x_{i}})|x_{l}\in\\{x_{m_{1}+1}(l),\ldots,x_{n_{1}}(l),y_{m_{2}+1}(l),\ldots,y_{n_{2}}(l)\\},l\in
U_{x_{i}}\\},$

for any $i=1,\ldots,m_{1}$. Let $[d]=\\{1,\ldots,d\\}$. Notice that

$\displaystyle T_{1}(x_{i}(l);l\in U_{x_{i}})$
$\displaystyle=\sum_{j=m_{1}+1}^{n_{1}}a_{i,j}^{\prime}\exp\left(-\beta\sum_{l\in
U_{x_{i}}}|x_{i}(l)-x_{j}(l)|\right)-$
$\displaystyle\sum_{j=m_{2}+1}^{n_{2}}b_{i,j}^{\prime}\exp\left(-\beta\sum_{l\in
U_{x_{i}}}|x_{i}(l)-y_{j}(l)|\right),$

where

$\displaystyle a_{i,j}^{\prime}$
$\displaystyle=c_{1}\exp\left(-\beta\sum_{l\in[d]\setminus
U_{x_{i}}}|x_{i}(l)-x_{j}(l)|\right),$ $\displaystyle b_{i,j}^{\prime}$
$\displaystyle=c_{3}\exp\left(-\beta\sum_{l\in[d]\setminus
U_{x_{i}}}|x_{i}(l)-y_{j}(l)|\right).$

Then, according to Lemma 2, it follows

$\displaystyle\max\\{0,\max\mathcal{X}_{x_{i}}\\}\geq T_{1}(x_{i}(l);l\in
U_{x_{i}})\geq\min\\{0,\min\mathcal{X}_{x_{i}}\\},i=1,\ldots,m_{1}.$ (25)

According to Lemma 3,

$\displaystyle\min\mathcal{X}_{x_{i}}\geq\sum_{j=m_{1}+1}^{n_{1}}a_{i,j}^{\prime}\exp\left(-\beta\sum_{l\in
U_{x_{i}}}\tilde{x}_{i}(l)\right)-\sum_{j=m_{2}+1}^{n_{2}}b_{i,j}^{\prime},$
$\displaystyle\max\mathcal{X}_{x_{i}}\leq\sum_{j=m_{1}+1}^{n_{1}}a_{i,j}^{\prime}-\sum_{j=m_{2}+1}^{n_{2}}b_{i,j}^{\prime}\exp\left(-\beta\sum_{l\in
U_{x_{i}}}\tilde{y}_{j}(l)\right).$

Notice that according to the definition of incomplete kernel in Definition 3,

$\displaystyle a_{i,j}^{\prime}$
$\displaystyle=c_{1}\exp\left(-\beta\sum_{l\in[d]\setminus
U_{x_{i}}}|x_{i}(l)-x_{j}(l)|\right)=c_{1}k^{*}(x_{i},x_{j}),$ $\displaystyle
b_{i,j}^{\prime}$ $\displaystyle=c_{3}\exp\left(-\beta\sum_{l\in[d]\setminus
U_{x_{i}}}|x_{i}(l)-y_{j}(l)|\right)=c_{3}k^{*}(x_{i},y_{j}).$

Hence, for any $i=1,\ldots,m_{1}$,

$\displaystyle\min\mathcal{X}_{x_{i}}\geq\sum_{j=m_{1}+1}^{n_{1}}c_{1}k^{*}(x_{i},x_{j})\exp\left(-\beta\sum_{l\in
U_{x_{i}}}\tilde{x}_{i}(l)\right)-\sum_{j=m_{2}+1}^{n_{2}}c_{3}k^{*}(x_{i},y_{j}),$
$\displaystyle\max\mathcal{X}_{x_{i}}\leq\sum_{j=m_{1}+1}^{n_{1}}c_{1}k^{*}(x_{i},x_{j})-\sum_{j=m_{2}+1}^{n_{2}}c_{3}k^{*}(x_{i},y_{j})\exp\left(-\beta\sum_{l\in
U_{x_{i}}}\tilde{y}_{j}(l)\right).$

According to (25), we then have

$\displaystyle T_{1}(x_{i}(l);l\in
U_{x_{i}})\geq\min\\{0,\min\mathcal{X}_{x_{i}}\\}$
$\displaystyle\geq\min\left\\{0,\sum_{j=m_{1}+1}^{n_{1}}c_{1}k^{*}(x_{i},x_{j})\exp\left(-\beta\sum_{l\in
U_{x_{i}}}\tilde{x}_{i}(l)\right)-\sum_{j=m_{2}+1}^{n_{2}}c_{3}k^{*}(x_{i},y_{j})\right\\},$

and

$\displaystyle T_{1}(x_{i}(l);l\in
U_{x_{i}})\leq\max\\{0,\max\mathcal{X}_{x_{i}}\\}$
$\displaystyle\leq\max\left\\{0,\sum_{j=m_{1}+1}^{n_{1}}c_{1}k^{*}(x_{i},x_{j})-\sum_{j=m_{2}+1}^{n_{2}}c_{3}k^{*}(x_{i},y_{j})\exp\left(-\beta\sum_{l\in
U_{x_{i}}}\tilde{y}_{j}(l)\right)\right\\}.$

Following (24), then

$\displaystyle C_{1}\geq A_{3}\geq C_{2}.$

Using the similar method, it can be prove that

$\displaystyle C_{3}\geq A_{4}\geq C_{4}.$

Combining with (25), we conclude our result. ∎

For computing bounds of ${\text{MMD}}^{2}_{u}(X,Y)$ using Theorem 5, we need
to compute $\tilde{x}_{i}(l)$ for any
$i\in\\{m_{1}+1,\ldots,n_{1}\\},l\in\\{1,\ldots,d\\}$, and compute
$\tilde{y}_{i}(l)$ for any
$i\in\\{m_{2}+1,\ldots,n_{2}\\},l\in\\{1,\ldots,d\\}$. Notice that for any
$\tilde{x}_{i}(l)$, we have

$\displaystyle|\\{|x_{i}(l)-x_{m_{1}+1}(l)|,\ldots,|x_{i}(l)-x_{n_{1}}(l)|,|x_{i}(l)-y_{m_{2}+1}(l)|,\ldots,|x_{i}(l)-y_{n_{2}}(l)|\\}|$
$\displaystyle=n_{1}-m_{1}+n_{2}-m_{2}.$

Hence, computing all $\tilde{x}_{i}(l)$ is of computational order at most
$O(dn_{1}(n_{1}+n_{2}))$. Similarly, computing all $\tilde{y}_{i}(l)$ is of
computational complexity order at most $O(dn_{2}(n_{1}+n_{2}))$. After
computing $\tilde{x}_{i}(l)$ for any
$i\in\\{m_{1}+1,\ldots,n_{1}\\},l\in\\{1,\ldots,d\\}$, and $\tilde{y}_{i}(l)$
for any $i\in\\{m_{2}+1,\ldots,n_{2}\\},l\in\\{1,\ldots,d\\}$,
$C_{1},C_{2},C_{3}$ and $C_{4}$ can be computed with computational complexity
order at most $O((n_{1}+n_{2}))$. Overall, the computational complexity order
of computing bounds of ${\text{MMD}}^{2}_{u}(X,Y)$ using Theorem 5, is
$O((n_{1}+n_{2}))$. This concludes our proof for Theorem 1.T

###  A.7 Proof of Theorem 2

This section provides bounds of $p$-value of MMD test statistic based on the
bounds of MMD derived in Section 3, and the permutation test.

######  Theorem 2.

Suppose $X=\\{x_{1},\ldots,x_{n_{1}}\\}$ and $Y=\\{y_{1},\ldots,y_{n_{2}}\\}$
are samples of $\mathbb{R}^{d}$ real value. Let
$(\sigma^{(1)},\ldots,\sigma^{(B)})$ be $B$ i.i.d. random samplings, each
being random permutation of $\\{1,\ldots,n_{1}+n_{2}\\}$ and denoted as
$\sigma^{(i)}=(\sigma^{(i)}(1),\ldots,\sigma^{(i)}(n_{1}+n_{2}))$,
$i=1,\ldots,B$. Subsequently, let
$z_{1}=x_{1},\ldots,z_{n_{1}}=x_{n_{1}},z_{n_{1}+1}=y_{1},\ldots,z_{n_{1}+n_{2}}=y_{n_{2}}$
and for any $i=1,\ldots,B$, denote

$\displaystyle
X_{\sigma^{(i)}}=\\{z_{\sigma^{(i)}(1)},\ldots,z_{\sigma^{(i)}(n_{1})}\\},Y_{\sigma^{(i)}}=\\{z_{\sigma^{(i)}(n_{1}+1)},\ldots,z_{\sigma^{(i)}(n_{1}+n_{2})}\\},$

and define $p$ according to equation (2). Suppose further

$\displaystyle{\underline{\text{MMD}}}^{2}_{u}(X,Y)\leq{\text{MMD}}^{2}_{u}(X,Y),~{}{\text{MMD}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\leq{\overline{\text{MMD}}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}}),~{}i=1,\ldots,B.$

Define

$\displaystyle\overline{p}=\frac{1}{B+1}\left(1+\sum_{i=1}^{B}\text{I}({\overline{\text{MMD}}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\underline{\text{MMD}}}^{2}_{u}(X,Y))\right).$

Then, we must have $\overline{p}\geq p$.

###### Proof.

Notice that for any $i\in\\{1,\ldots,B\\}$, we have

$\displaystyle\text{I}({\text{MMD}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\text{MMD}}^{2}_{u}(X,Y))=1$
$\displaystyle\Rightarrow$
$\displaystyle{\text{MMD}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\text{MMD}}^{2}_{u}(X,Y).$

Since

$\displaystyle{\text{MMD}}^{2}_{u}(X,Y)\geq{\underline{\text{MMD}}}^{2}_{u}(X,Y),~{}{\overline{\text{MMD}}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\text{MMD}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}}),~{}i=1,\ldots,B,$

it then follows

$\displaystyle{\overline{\text{MMD}}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\text{MMD}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\underline{\text{MMD}}}^{2}_{u}(X,Y),$
$\displaystyle\Rightarrow$
$\displaystyle\text{I}({\overline{\text{MMD}}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\underline{\text{MMD}}}^{2}_{u}(X,Y))=1.$

for any $i\in\\{1,\ldots,B\\}$.

If, for any $i\in\\{1,\ldots,B\\}$,

$\displaystyle\text{I}({\text{MMD}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\text{MMD}}^{2}_{u}(X,Y))=0,$

then we have

$\displaystyle\text{I}({\overline{\text{MMD}}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\underline{\text{MMD}}}^{2}_{u}(X,Y))=0~{}\text{or}~{}1.$

Hence, for any $i\in\\{1,\ldots,B\\}$, we have

$\displaystyle\text{I}({\overline{\text{MMD}}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\underline{\text{MMD}}}^{2}_{u}(X,Y))\geq\text{I}({\text{MMD}}^{2}_{u}(X_{\sigma^{(i)}},Y_{\sigma^{(i)}})\geq{\text{MMD}}^{2}_{u}(X,Y)).$

According to the definitions of $p$ and $\overline{p}$, we then must have
$\overline{p}\geq p$.

∎

###  A.8 Proof of Theorem 3

This section provides bounds of $p$-value of MMD test statistic based on the
bounds of MMD derived in Section 3, and the normality approximation proposed
in [8].

######  Theorem 3.

Following the notation and definitions in Proposition 10 in [8], let $k$
denote the Laplacian kernel. Assume $X_{1},\ldots,X_{m_{1}}$ and
$Y_{1},\ldots,Y_{m_{2}}$ are not observed. For $1\leq s,t\leq N$ let us
further denote

$\displaystyle\overline{a}_{s,t}^{k}=\left\\{\begin{array}[]{ll}k^{*}(X_{s},X_{t})&1\leq
s,t\leq n,s\neq t\\\ k^{*}(X_{s},Y_{t-n})&1\leq s\leq n<t\leq N,s\neq t\\\
k^{*}(X_{t},Y_{s-n})&1\leq t\leq n<s\leq N,\\\
k^{*}(Y_{s-n},Y_{t-n})&n+1<s,t\leq
N,\end{array}\right.,\underline{a}_{s,t}^{k}=\left\\{\begin{array}[]{ll}0&1\leq
s\leq m_{1},1\leq t\leq N,\\\ 0&1\leq s\leq N,n+1\leq t\leq n+m_{2},\\\
1&s=t\\\ a_{s,t}^{k}&\mbox{Otherwise}.\end{array}\right.$

Let
$\overline{A}_{s,t}^{k*}=\overline{a}_{s,t}^{k}-\underline{a}_{\cdot,t}^{k}-\underline{a}_{s,\cdot}^{k}+\overline{a}_{\cdot,\cdot}^{k}$,
where

$\displaystyle\underline{a}_{\cdot,t}^{k}=\frac{1}{N-2}\sum_{i=1}^{N}\underline{a}_{i,t}^{k},~{}\underline{a}_{s,\cdot}^{k}=\frac{1}{N-2}\sum_{j=1}^{N}\underline{a}_{s,j}^{k},~{}\overline{a}_{\cdot,\cdot}^{k}=\frac{1}{(N-1)(N-2)}\sum_{i,j=1}^{N}\overline{a}_{i,j}^{k}.$

Define

$\displaystyle\overline{\mathcal{V}}_{n,m}^{k*}=\frac{1}{N(N-3)}\sum_{s\neq
t}\left(\overline{A}_{s,t}^{k*}\right)^{2}-\frac{1}{(N-1)(N-3)}.$

Then, we must have
$\overline{\mathcal{V}}_{n,m}^{k*}\geq{\mathcal{V}}_{n,m}^{k*},$ Further,
suppose

$\displaystyle
0\leq{\underline{\text{MMD}}}^{2}_{u}(X,Y)\leq{\text{MMD}}^{2}_{u}(X,Y),$

we then have

$\displaystyle\frac{{\underline{\text{MMD}}}^{2}_{u}(X,Y)}{\sqrt{c_{n,m}\overline{\mathcal{V}}_{n,m}^{k*}}}\leq\frac{{\text{MMD}}^{2}_{u}(X,Y)}{\sqrt{c_{n,m}{\mathcal{V}}_{n,m}^{k*}}}.$

###### Proof.

To start, recall that it is proved in Lemma 5 that

$\displaystyle k^{*}(x,y)\geq k(x,y)>0,$

where $x,y$ are potentially incompletely observed $\mathbb{R}^{d}$ samples.
Hence, for any $1\leq s,t\leq N$, we have

$\displaystyle\overline{a}_{s,t}^{k}\geq{a}_{s,t}^{k}\geq\underline{a}_{\cdot,t}^{k},$

which then gives

$\displaystyle\overline{A}_{s,t}^{k*}\geq{A}_{s,t}^{k*},$

where ${A}_{s,t}^{k*}$ is defined in Proposition 10 in [8] as

$\displaystyle{A}_{s,t}^{k*}:={a}_{s,t}^{k}-{a}_{\cdot,t}^{k}-{a}_{s,\cdot}^{k}+{a}_{\cdot,\cdot}^{k}.$

Subsequently,

$\displaystyle\overline{\mathcal{V}}_{n,m}^{k*}\geq{\mathcal{V}}_{n,m}^{k*},$

where ${\mathcal{V}}_{n,m}^{k*}$ is defined in Proposition 10 in [8] as

$\displaystyle{\mathcal{V}}_{n,m}^{k*}=\frac{1}{N(N-3)}\sum_{s\neq
t}\left({A}_{s,t}^{k*}\right)^{2}-\frac{1}{(N-1)(N-3)}.$

Further, suppose

$\displaystyle
0\leq{\underline{\text{MMD}}}^{2}_{u}(X,Y)\leq{\text{MMD}}^{2}_{u}(X,Y).$

Then, we must have

$\displaystyle\frac{{\underline{\text{MMD}}}^{2}_{u}(X,Y)}{\sqrt{c_{n,m}\overline{\mathcal{V}}_{n,m}^{k*}}}\leq\frac{{\text{MMD}}^{2}_{u}(X,Y)}{\sqrt{c_{n,m}{\mathcal{V}}_{n,m}^{k*}}},$

which completes our proof. ∎

In [8], the $p$-value is defined as
$p=1-\Phi\left(\frac{{\text{MMD}}^{2}_{u}(X,Y)}{\sqrt{c_{n,m}{\mathcal{V}}_{n,m}^{k*}}}\right)$,
where $\Phi$ denotes the cumulative density function of the standard normal
distribution. In the presence of missing data, we can define
$\overline{p}=1-\Phi\left(\frac{{\underline{\text{MMD}}}^{2}_{u}(X,Y)}{\sqrt{c_{n,m}\overline{\mathcal{V}}_{n,m}^{k*}}}\right)$.
Then, according to Theorem 3, we must have $\overline{p}\geq p$ when
$0\leq{\underline{\text{MMD}}}^{2}_{u}(X,Y)$. When
${\underline{\text{MMD}}}^{2}_{u}(X,Y)<0$, it suggests that the
${\text{MMD}}^{2}_{u}(X,Y)$ could potentially be smaller than 0. Hence, we
then define $\underline{p}=1$ in this case, i.e. the null hypothesis will not
be rejected if ${\underline{\text{MMD}}}^{2}_{u}(X,Y)<0$.

##  Appendix B Simulation Details

###  B.1 Computing Environment and Resources

The experiments were run on an high performance computing cluster with 325
compute nodes, each equipped with 2x AMD EPYC 7742 processors (128 cores, 1TB
RAM per node). Each job utilized 32 cores and 64 GB of memory, with compute
times varying from a couple of hours to 2-3 days, based on the samples sizes
and dimensions of data of the job.

###  B.2 Simulation Settings

#### Median Heuristic.

For all the experiments in the section 5, Laplacian kernel is used and median
heuristic [9] will be used to compute the parameter $\beta$ of the kernel.
This is a popular method of computing $\beta$, and it generally works well
[9]. For case deletion and the proposed methods, only completely observed
samples will be used for computing $\beta$. For mean and hot deck imputation,
the $\beta$ will be computed using samples after imputations.

#### Case deletion.

For univariate data, case deletion uses only the observed data for testing.
For multivariate data, only data that are completely observed are used.

#### Mean imputation.

In univariate data, mean imputation replaces missing values with the mean of
the observed samples in the same group, $X$ or $Y$. For multivariate data, the
missing components of an incomplete sample are imputed with the mean of the
observed components of that same sample.

#### Hot deck imputation.

In univariate data, hot deck imputation imputes missing values with randomly
selected observed values (with replacement) from the same group, $X$ or $Y$.
For multivariate data, the missing components of an incomplete sample are
imputed with randomly selected observed components (with replacement) from
that same sample.

#### Implementations of common missing data approaches.

After performing case deletion, mean imputation, or hot deck imputation,
standard MMD two-sample testing with a Laplacian kernel is conducted, with
parameter $\beta$ computed using the median heuristic. Permutation tests with
$B=100$ are used to compute the $p$-value.

#### Implementations of proposed methods.

In Theorem 2, it is proved that the $p$-value of ${\text{MMD}}^{2}_{u}(X,Y)$,
denoted as $p$, can be bounded by $\overline{p}$, which can be computed using
the bounds of ${\text{MMD}}^{2}_{u}(X,Y)$ in the presence of missing data. The
bounds of ${\text{MMD}}^{2}_{u}(X,Y)$ are computed according to Theorem 4 when
$d=1$ and according to Theorem 5 when $d>1$. Permutation times $B=100$ will be
used. The null hypothesis is rejected when $\overline{p}$ is larger or equal
than the significance level $\alpha$. In Section 5, we denote this proposed
method based on permutations as MMD-Miss: Perm.

When the sample sizes $n_{1},n_{2}\geq 100,d\geq 50$, the normality
approximation proposed in [8] may be suitable for computing the $p$-value of
the MMD test. In Theorem 3, we provide the upper bound the $p$-value, denoted
as $\overline{p}$, by first bounding the value of estimation of variance and
combining with the lower bounds of ${\text{MMD}}^{2}_{u}(X,Y)$ provided in
Section 3. The null hypothesis will then be rejected when $\overline{p}$ is
larger or equal than the significance level $\alpha$. In Section 5, we denote
this proposed method based on normality approximation as MMD-Miss: Normality.

#### Implementations for generating missing data.

The missingness mechanisms for each experiment are detailed in Section 5. If
there are not enough samples in $X$ or $Y$ satisfying the conditions to be
missing, randomly missing will be introduced. For example, in Figure 1,
$s\in\\{0,\ldots,0.2\\}$ proportion of samples in $X$ will be selected and
labeled as missing and the missingness mechanism is that only sample smaller
than 0 in $X$ will be randomly selected. If the proportion of samples smaller
than 0 in $X$, denoted as $a$, is smaller than a given $s$, then all samples
in $X$ smaller than 0 will be selected and $\lfloor(s-a)n_{1}\rfloor$ number
of unselected samples in $X$ will be randomly selected.

###  B.3 Examples of MNIST dataset.

In this section, we present examples of MNIST images used in Section 5.

Figure 4: Examples of MNIST images [15]. Each image has dimensions of
$d=28\times 28$. The labels for the images, ordered from left to right and top
to bottom, are 5, 0, 4, 1, 9, 2, 1, 3, 1, 4, 3, 5, 3, 6, 1, and 7. Figure 5:
Examples of incompletely observed MNIST images, where white pixels represent
unobserved values. Each image has dimensions of $d=28\times 28$ and is labeled
as 3. The first column shows fully observed images; the second column shows
images with upper half pixels missing (rows 1 to 14 and columns 8 to 21); the
third column shows images after mean imputation; and the fourth column shows
images after hot deck imputation. Each row presents different variations of
the same original image.

##  Appendix C Additional Experiments

In this section, additional experiments are provided for investigating the
power of proposed method. Different sample sizes and alternatives are
considered compared to those in Section 5.

###  C.1 Univariate Data

The missingness mechanism for producing the following two figures is the same
as the mechanisms used to produce Figure 2.

Figure 6: The power of MMD-Miss and the three common missing data approaches
for univariate samples when data are missing not at random (MNAR). (Left):
Sample size $n_{1}=n_{2}=50$; (Middle): Sample size $n_{1}=n_{2}=100$;
(Right): Sample size $n_{1}=n_{2}=200$. For all figures, significance level
$\alpha=0.05$, and alternative hypothesis $N(0,1)$ vs $N(1,1)$ are used. The
plot values are the average times of rejecting the null hypothesis over 100
repetitions. Figure 7: The power of MMD-Miss and the three common missing data
approaches for univariate samples when data are missing not at random (MNAR).
(Left): Sample size $n_{1}=n_{2}=50$; (Middle): Sample size $n_{1}=n_{2}=100$;
(Right): Sample size $n_{1}=n_{2}=200$. For all figures, significance level
$\alpha=0.05$, and alternative hypothesis $N(0,1)$ vs $N(0,3)$ are used. The
plot values are the average times of rejecting the null hypothesis over 100
repetitions.

###  C.2 Multivariate Data

The missingness mechanism for producing the following two figures is the same
as the mechanisms used to produce Figure 3.

Figure 8: The power of MMD-Miss and the three common missing data approaches
for multivariate samples $(d=50)$ when data are missing not at random (MNAR).
(Left): Sample size $n_{1}=n_{2}=50$; (Middle): Sample size $n_{1}=n_{2}=100$;
(Right): Sample size $n_{1}=n_{2}=200$. For all figures, significance level
$\alpha=0.05$, and alternative hypothesis $N((0,\ldots,0)^{T},I_{50})$ vs
$N((1,\ldots,1)^{T},I_{50})$ are used. The plot values are the average times
of rejecting the null hypothesis over 100 repetitions. Figure 9: The power of
MMD-Miss and the three common missing data approaches for multivariate samples
$(d=50)$ when data are missing not at random (MNAR). (Left): Sample size
$n_{1}=n_{2}=50$; (Middle): Sample size $n_{1}=n_{2}=100$; (Right): Sample
size $n_{1}=n_{2}=200$. For all figures, significance level $\alpha=0.05$, and
alternative hypothesis $N((0,\ldots,0)^{T},I_{50})$ vs
$N((0,\ldots,0)^{T},9I_{50})$ are used. The plot values are the average times
of rejecting the null hypothesis over 100 repetitions.

## References

  * [1] N. Aronszajn.  Theory of reproducing kernels.  Transactions of the American mathematical society, 68(3):337–404, 1950.
  * [2] F. Biggs, A. Schrab, and A. Gretton.  MMD-fuse: Learning and combining kernels for two-sample testing without data splitting.  Advances in Neural Information Processing Systems, 36, 2024.
  * [3] D. A. Bodenham and Y. Kawahara.  euMMD: Efficiently computing the mmd two-sample test statistic for univariate data.  Statistics and Computing, 33(5):110, 2023.
  * [4] S. X. Chen, J. Qin, and C. Y. Tang.  Mann-Whitney test with adjustments to pretreatment variables for missing values and observational study.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 75(1):81–102, 2013.
  * [5] Y. K. Cheung.  Exact two-sample inference with missing data.  Biometrics, 61(2):524–531, 2005.
  * [6] Y. Fong, Y. Huang, M. P. Lemos, and M. J. McElrath.  Rank-based two-sample tests for paired data with missing values.  Biostatistics, 19(3):281–294, 2018.
  * [7] K. Fukumizu, A. Gretton, B. Schölkopf, and B. K. Sriperumbudur.  Characteristic kernels on groups and semigroups.  Advances in neural information processing systems, 21, 2008.
  * [8] H. Gao and X. Shao.  Two sample testing in high dimension via maximum mean discrepancy.  Journal of Machine Learning Research, 24(304):1–33, 2023.
  * [9] A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola.  A kernel two-sample test.  The Journal of Machine Learning Research, 13(1):723–773, 2012.
  * [10] Arthur Gretton, Kenji Fukumizu, Zaid Harchaoui, and Bharath K Sriperumbudur.  A fast, consistent kernel two-sample test.  Advances in neural information processing systems, 22, 2009.
  * [11] M. W. Heymans and J. W. R. Twisk.  Handling missing data in clinical research.  Journal of clinical epidemiology, 151:185–188, 2022.
  * [12] J. D. Jobson and B. M. Korkie.  Performance hypothesis testing with the sharpe and treynor measures.  Journal of Finance, pages 889–908, 1981.
  * [13] A. N. Kolmogorov.  Sulla determinazione empirica di una legge didistribuzione.  Giorn Dell’inst Ital Degli Att, 4:89–91, 1933.
  * [14] J. M. Lachin.  Worst-rank score analysis with informatively missing observations in clinical trials.  Controlled clinical trials, 20(5):408–422, 1999.
  * [15] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [16] R. J. A. Little and D. B. Rubin.  Statistical analysis with missing data, volume 793.  John Wiley & Sons, 2019.
  * [17] Ž. Lukić and B. Milošević.  A novel two-sample test within the space of symmetric positive definite matrix distributions and its application in finance.  Annals of the Institute of Statistical Mathematics, pages 1–24, 2024.
  * [18] H. B. Mann and D. R. Whitney.  On a test of whether one of two random variables is stochastically larger than the other.  The Annals of Mathematical Statistics, pages 50–60, 1947.
  * [19] P. Martínez-Camblor, N. Corral, and J. Maria de la Hera.  Hypothesis test for paired samples in the presence of missing data.  Journal of Applied Statistics, 40(1):76–87, 2013.
  * [20] P. C. O’Brien and T. R. Fleming.  A multiple testing procedure for clinical trials.  Biometrics, pages 549–556, 1979.
  * [21] W. Pan.  A two-sample test with interval censored data via multiple imputation.  Statistics in Medicine, 19(1):1–11, 2000.
  * [22] R. Pfister and M. Janczyk.  Confidence intervals for two sample means: Calculation, interpretation, and a few simple rules.  Advances in Cognitive Psychology, 9(2):74, 2013.
  * [23] M. L. Rizzo and G. J. Székely.  Energy distance.  wiley interdisciplinary reviews: Computational statistics, 8(1):27–38, 2016.
  * [24] D. B. Rubin.  Inference and missing data.  Biometrika, 63(3):581–592, 1976.
  * [25] D. B. Rubin and N. Schenker.  Multiple imputation for interval estimation from simple random samples with ignorable nonresponse.  Journal of the American statistical Association, 81(394):366–374, 1986.
  * [26] J. L. Schafer.  Multiple imputation: a primer.  Statistical methods in medical research, 8(1):3–15, 1999.
  * [27] J. L Schafer and J. W. Graham.  Missing data: our view of the state of the art.  Psychological methods, 7(2):147, 2002.
  * [28] A. Schrab, I. Kim, M. Albert, B. Laurent, B. Guedj, and A. Gretton.  MMD aggregated two-sample test.  Journal of Machine Learning Research, 24(194):1–81, 2023.
  * [29] N. Shutoh, M. Kusumi, W. Morinaga, S. Yamada, and T. Seo.  Testing equality of mean vectors in two sample problem with missing data.  Communications in Statistics—Simulation and Computation®, 39(3):487–500, 2010.
  * [30] B. K. Sriperumbudur, K. Fukumizu, and G. R. G. Lanckriet.  Universality, characteristic kernels and rkhs embedding of measures.  Journal of Machine Learning Research, 12(7), 2011.
  * [31] Student.  The probable error of a mean.  Biometrika, pages 1–25, 1908.
  * [32] A. A. Tsiatis, M. Davidian, M. Zhang, and X. Lu.  Covariate adjustment for two-sample treatment comparisons in randomized clinical trials: a principled yet flexible approach.  Statistics in medicine, 27(23):4658–4677, 2008.
  * [33] Y. Zeng, N. M. Adams, and D. A. Bodenham.  On two-sample testing for data with arbitrarily missing values.  arXiv preprint arXiv:2403.15327, 2024.