In recent years, there has been a growing body of literature aimed at investigating the power of MMD-based tests and enhancing their performance. For example, the work of Li and Yuan, (2019); Balasubramanian et al., (2021) demonstrated that MMD tests equipped with a fine-tuned kernel can achieve minimax optimality with respect to the $L_{2}$ separation in an asymptotic sense. To establish a similar but non-asymptotic guarantee, Schrab et al., (2023) introduced a MMD aggregated test calibrated by using either permutations or a wild bootstrap. It is also worth noting that the minimax optimality of MMD two-sample tests has been established for separations other than the $L_{2}$ distance, such as MMD distance (Kim and Schrab,, 2023), and Hellinger distance (Hagrass et al.,, 2022). In addition to these works, several other MMD-based minimax tests have been proposed using techniques such as aggregation (Fromont et al.,, 2013; Chatterjee and Bhattacharya,, 2023; Biggs et al.,, 2023) and studentization (Kim and Ramdas,, 2024; Shekhar et al.,, 2023). Despite significant recent advancements made in this field, the quadratic time complexity of these methods remains a barrier in large-scale applications, which highlights the need for more efficient yet powerful testing approaches.

To address the computational concern of quadratic-time MMD tests, several time-efficient approaches have emerged, which leverage subsampled estimation techniques, such as linear-time statistics (Gretton et al., 2012a, ; Gretton et al., 2012b, ), block-based statistics (Zaremba et al.,, 2013) and incomplete U-statistics (Yamada et al.,, 2019; Schrab et al.,, 2022). However, in terms of power, these methods are either sub-optimal or ultimately require quadratic time complexity to achieve optimality (Domingo-Enrich et al.,, 2023, Proposition 2). Other advancements in accelerating two-sample tests have involved techniques, such as Nyström approximations (Chatalic et al.,, 2022), analytic mean embeddings and smoothed characteristic functions (Chwialkowski et al.,, 2015; Jitkrittum et al.,, 2016), deep linear kernels (Kirchler et al.,, 2020), as well as random Fourier features (Zhao and Meng,, 2015). These tests can also run in sub-quadratic time, while their theoretical guarantees on power remain largely unknown. We also mention the recent method using kernel thinning (Dwivedi and Mackey,, 2021; Domingo-Enrich et al.,, 2023), which achieves the same MMD separation rate as the quadratic-time test but with sub-quadratic running time. However, this guarantee is valid under specific distributional assumptions that differ from those we consider. Moreover, their result focuses solely on alternatives that deviate from the null in terms of the MMD metric.

With this context in mind, we revisit the RFF-MMD test (Zhao and Meng,, 2015) and delve into its time-power trade-off concerning the number of random features. Despite an extensive body of literature on random features for kernel approximation, prior work has mainly focused on the estimation quality of kernel approximation (Rahimi and Recht,, 2007; Zhao and Meng,, 2015; Sriperumbudur and Szabo,, 2015; Sutherland and Schneider,, 2015; Yao et al.,, 2023), and a theoretical guarantee on the power of the RFF-MMD test has not been explored. In this work, we seek to bridge this gap by thoroughly analyzing the trade-off between computation time and statistical power in the context of the RFF-MMD test.

Having reviewed the prior work, we now summarize the key contributions of this paper.

• •

  Inconsistency result for RFF-MMD (Section 3). We first investigate the setting where the number of random Fourier features is fixed, and demonstrate that the RFF-MMD test fails to achieve pointwise consistency (Theorem 3 and Corollary 4). Concretely, we prove that there exist infinitely many pairs of distinct distributions for which the power of the RFF-MMD test using a fixed number of random Fourier features is almost equal to the size even asymptotically.

• •

  Sufficient conditions for consistency (Section 4). Our previous negative result clearly indicates that increasing the number of random Fourier features is necessary to achieve pointwise consistency. In Theorem 5, we show that it is indeed sufficient to increase the number of Fourier features to infinity to achieve pointwise consistency, even at an arbitrarily slow rate.

• •

  Time-power trade-off (Section 4). As mentioned before, there exists a clear trade-off between computational efficiency and statistical power in terms of the number of random Fourier features. To balance this trade-off, we adopt the non-asymptotic minimax testing framework and analyze how changes in the number of random Fourier features impact both computational efficiency and separation rates in terms of the $L_{2}$ metric (Theorem 6) and the MMD metric (Theorem 7).

• •

  Achieving optimality in sub-quadratic time (Section 4). We firmly demonstrate in Theorem 6 that it is possible to achieve the minimax separation rate against $L_{2}$ alternatives in sub-quadratic time when the underlying distributions are sufficiently smooth. Similarly, we establish in Proposition 8 that a parametric separation rate against MMD alternatives can be achieved in linear time for certain classes of distributions including Gaussian distributions.

Our theoretical results are validated through simulation studies under various scenarios and the code that reproduces our numerical results can be found at https://github.com/ikjunchoi/rff-mmd.

Let $\mathcal{X}_{n_{1}}:=\{X_{i}\}^{n_{1}}_{i=1}$ be ${n_{1}}$ i.i.d. random samples from the distribution $P_{X}$, and $\mathcal{Y}_{n_{2}}:=\{Y_{j}\}^{n_{2}}_{j=1}$ be ${n_{2}}$ i.i.d. random samples from the distribution $P_{Y}$ where ${n_{1}},{n_{2}}\geq 2$. Based on these mutually independent samples, the problem of two-sample testing is concerned with determining whether $P_{X}$ and $P_{Y}$ agree or not. More formally, let $\mathcal{P}$ be a class of all possible pairs of distributions on some generic space $\mathbb{S}$, and consider two disjoint subsets in $\mathcal{P}$, namely $\mathcal{P}_{0}:=\left\{(P_{X},P_{Y})\in\mathcal{P}\,|\,P_{X}=P_{Y}\right\}$ and $\mathcal{P}_{1}:=\left\{(P_{X},P_{Y})\in\mathcal{P}\,|\,P_{X}\neq P_{Y}\right\}$. Then, the null hypothesis $H_{0}$ and the alternative hypothesis $H_{1}$ of two-sample testing can be formulated as follows:

$$H_{0}:(P_{X},P_{Y})\in\mathcal{P}_{0}\quad\text{vs.}\quad H_{1}:(P_{X},P_{Y})\in\mathcal{P}_{1}.$$

In order to decide whether to reject $H_{0}$ or not, we devise a test function $\Delta_{{n_{1}},{n_{2}}}:(\mathbb{S}^{n_{1}},\,\mathbb{S}^{{n_{2}}})\rightarrow\{0,1\}$, and reject the null hypothesis if and only if $\Delta_{{n_{1}},{n_{2}}}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}})=1$. This decision-making process naturally leads to two types of errors, which we would like to minimize. The first error, called the type I error, occurs by rejecting the null hypothesis despite being true. Conversely, the second error, called the type II error, arises when the null hypothesis is accepted despite being false. One common approach to design an ideal test is to first bound the probability of the type I error uniformly over $\mathcal{P}_{0}$ as

$$\sup_{(P_{X},P_{Y})\in\mathcal{P}_{0}}\mathbb{P}_{X\times Y}(\Delta_{{n_{1},n_{2}}}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}})=1)\leq\alpha,\quad\text{for a given {level} $\alpha\in(0,1),$}$$

where $\mathbb{P}_{X\times Y}$ denotes the probability operator over $\mathcal{X}_{n_{1}}\overset{\mathrm{i.i.d.}}{\sim}P_{X}$ and $\mathcal{Y}_{n_{2}}\overset{\mathrm{i.i.d.}}{\sim}P_{Y}$. We say that such a test is a level-$\alpha$ test. Next, our focus shifts to controlling the type II error. Given a fixed pair $(P_{X},P_{Y})$ in $\mathcal{P}_{1}$ and a level-$\alpha$ test $\Delta^{\alpha}_{{n_{1},n_{2}}}$, suppose that the probability of the type II error is upper bounded by some constant $\beta\in(0,1)$. Equivalently, the probability of correctly rejecting the null, referred to as the power, is lower bounded by $1-\beta$. Ideally, we expect that the power of the test $\Delta^{\alpha}_{{n_{1},n_{2}}}$ against any fixed alternative $(P_{X},P_{Y})\in\mathcal{P}_{1}$ converges to one as we increase the data size ${n_{1}}$ and ${n_{2}}$. More formally, we desire a test $\Delta^{\alpha}_{{n_{1},n_{2}}}$ to be pointwise consistent, satisfying

$\displaystyle\lim_{{n_{1},n_{2}}\rightarrow\infty}\mathbb{P}_{X\times Y}\bigl{(}\Delta^{\alpha}_{{n_{1},n_{2}}}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}})=1\bigr{)}=1,\quad\text{for any fixed $(P_{X},P_{Y})\in\mathcal{P}_{1}$.}$

(1)

A stronger notion of the power is uniform consistency, guaranteeing that the power converges to one uniformly over a class of alternative distributions. See Section 4 for a discussion. For simplicity, in the rest of this paper, we consider $\mathbb{S}$ to be the $d$-dimensional Euclidean space denoted as $\mathbb{R}^{d}$.

As an example of integral probability metrics, the MMD measures the discrepancy between two distributions in a nonparametric manner. Specifically, given a reproducing kernel Hilbert space (RKHS) $\mathcal{H}_{k}$ equipped with a positive definite kernel $k$, the MMD between $P_{X}$ and $P_{Y}$ is defined as

$\displaystyle\mathrm{MMD}(P_{X},P_{Y};\mathcal{H}_{k}):=\sup_{f\in\mathcal{H}_{k}:\|f\|_{\mathcal{H}_{k}}\leq 1}\bigl{|}\mathbb{E}_{X}\left[f(X)\right]-\mathbb{E}_{Y}\left[f(Y)\right]\bigr{|}.$

It can also be represented as the RKHS distance between two mean embeddings of $P_{X}$ and $P_{Y}$, i.e., $\mathrm{MMD}(P_{X},P_{Y};\mathcal{H}_{k})=\|\mu_{X}-\mu_{Y}\|_{\mathcal{H}_{k}}$ where $\mu_{X}(\cdot):=\mathbb{E}_{X}[k(X,\cdot)]$ and $\mu_{Y}(\cdot):=\mathbb{E}_{Y}[k(Y,\cdot)]$. For a characteristic kernel $k$, the mean embedding of the kernel is injective (Sriperumbudur et al.,, 2010), which means that $\mathrm{MMD}(P_{X},P_{Y};\mathcal{H}_{k})=0$ if and only if $P_{X}=P_{Y}$. Among several ways to estimate the MMD, one straightforward way is to substitute the population mean embeddings $\mu_{X}$ and $\mu_{Y}$ with the empirical counterparts $\hat{\mu}_{X}(\cdot)=\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}k(X_{i},\cdot)$ and $\hat{\mu}_{Y}(\cdot)=\frac{1}{{n_{2}}}\sum_{i=1}^{n_{2}}k(Y_{i},\cdot)$. This plug-in approach results in a biased quadratic-time estimator of the squared MMD, also referred to as the V-statistic, given as

	$\displaystyle\widehat{\mathrm{MMD}}_{b}^{2}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\mathcal{H}_{k})$	$\displaystyle=\left\|\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}k(X_{i},\cdot)-\frac{1}{{n_{2}}}\sum_{i=1}^{n_{2}}k(Y_{i},\cdot)\right\|_{\mathcal{H}_{k}}^{2}$		(2)
		$\displaystyle=\frac{1}{{n_{1}^{2}}}\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{1}}k\left(X_{i},X_{j}\right)+\frac{1}{{n_{2}^{2}}}\sum_{i=1}^{n_{2}}\sum_{j=1}^{n_{2}}k\left(Y_{i},Y_{j}\right)-\frac{2}{{n_{1}}{n_{2}}}\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}k\left(X_{i},Y_{j}\right).$

Denoting ${N}:={n_{1}}+{n_{2}}$, this plug-in estimator requires a quadratic-time cost of $O({N}^{2}d)$ in terms of the sample size $N$ as it involves evaluating pairwise kernel similarities between samples. Another common approach to estimate $\mathrm{MMD}^{2}(P_{X},P_{Y};\mathcal{H}_{k})$ is using the U-statistic (e.g., Gretton et al., 2012a, , Lemma 6), which is given as

	$\displaystyle\widehat{\mathrm{MMD}}_{u}^{2}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\mathcal{H}_{k})\,=$	$\displaystyle\frac{1}{{n_{1}}({n_{1}}-1)}\sum_{1\leq i\neq j\leq{n_{1}}}k\left(X_{i},X_{j}\right)+\frac{1}{{n_{2}}({n_{2}}-1)}\sum_{1\leq i\neq j\leq{n_{2}}}k\left(Y_{i},Y_{j}\right)$
	$\displaystyle-$	$\displaystyle\frac{2}{{n_{1}}{n_{2}}}\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}k\left(X_{i},Y_{j}\right).$

This estimator is an unbiased estimator of the squared MMD and also requires quadratic-time computational costs.

Numerous approaches have been introduced to mitigate the computational cost of quadratic-time statistics often at the cost of sacrificing power performance. As reviewed in Section 1.1, some notable approaches include incomplete U-statistics (Gretton et al., 2012a, ; Zaremba et al.,, 2013; Schrab et al.,, 2022), Nyström approximations (Chatalic et al.,, 2022), kernel thinning (Dwivedi and Mackey,, 2021; Domingo-Enrich et al.,, 2023) and random Fourier features (Rahimi and Recht,, 2007; Zhao and Meng,, 2015). This work focuses on the method utilizing random Fourier features and investigates the effect of the number of random features on the power of a test. At the heart of this method is Bochner’s theorem (Lemma 9), which offers a means to approximate the kernel using a low-dimensional feature mapping $\psi_{\omega}$, satisfying $k(x,y)\approx\langle\psi_{\omega}(x),\psi_{\omega}(y)\rangle$. If a bounded continuous positive definite kernel $k$ is translation invariant on $\mathbb{R}^{d}$, that is, $k(x,y)=\kappa(x-y)$, Bochner’s theorem guarantees the existence of a nonnegative Borel measure $\Lambda$. It can be shown that $\Lambda$ is the inverse Fourier transform of $\kappa$ and meets

$$k(x,y)=\int_{\mathbb{R}^{d}}e^{\sqrt{-1}\omega^{\top}(x-y)}d\Lambda(\omega)\stackrel{{\scriptstyle(\dagger)}}{{=}}\int_{\mathbb{R}^{d}}\cos\bigl{(}{\omega^{\top}(x-y)}\bigr{)}d\Lambda(\omega),$$

where the equality $(\dagger)$ is obtained by the fact that $\kappa$ is both real and symmetric. Without loss of generality, we assume that $\Lambda$ is a probability measure, allowing the last integral to be expressed as $\mathbb{E}_{\omega\sim\Lambda}[\langle\psi_{\omega}(x),\psi_{\omega}(y)\rangle]$ with $\psi_{\omega}(x):=[\cos(\omega^{\top}x),\sin(\omega^{\top}x)]^{\top}$. If not, we instead work with the scaled versions of $\Lambda$ and $\psi_{\omega}$, given as $\Lambda^{\prime}:=\kappa^{-1}(0)\Lambda$ and $\psi_{\omega}^{\prime}(\cdot):=[\sqrt{\kappa(0)}\cos(\omega^{\top}\cdot),\sqrt{\kappa(0)}\sin(\omega^{\top}\cdot)]^{\top}$. In this case, $k(x,y)$ can be represented as $\mathbb{E}_{\omega\sim\Lambda^{\prime}}[\langle\psi_{\omega}^{\prime}(x),\psi_{\omega}^{\prime}(y)\rangle]$.

Now, by drawing a sequence of i.i.d. $R$ random frequencies $\boldsymbol{\omega}_{R}:=\{\omega_{r}\}^{R}_{r=1}$ from $\Lambda$, we construct an unbiased estimator of $k(x,y)$ defined as an inner product of random feature maps:

$$\hat{k}(x,y):=\frac{1}{R}\sum\limits_{r=1}^{R}\langle\psi_{\omega_{r}}(x),\psi_{\omega_{r}}(y)\rangle=\langle\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(x),\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(y)\rangle,$$

(3)

where $\psi_{\omega_{r}}(x)=[\cos(\omega_{r}^{\top}x),\sin(\omega_{r}^{\top}x)]^{\top}$ and $\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(x)=\frac{1}{\sqrt{R}}[\psi_{\omega_{1}}(x)^{\top},~{}\ldots~{},\psi_{\omega_{R}}(x)^{\top}]^{\top}\in\mathbb{R}^{2R}$. Let us define the vector in $\mathbb{R}^{2R}$ representing the difference in sample means of random feature maps as follows:

$$T(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R}):=\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(X_{i})-\frac{1}{{n_{2}}}\sum_{j=1}^{n_{2}}\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(Y_{j}).$$

Also, denote the quadratic form of $T:=T(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})$ as $V:=T^{\top}T$. When we replace the kernel $k$ in Equation (2) with the estimated $\hat{k}$, we obtain a RFF-MMD estimator of $\mathrm{MMD}^{2}$ that can run with a time complexity of $O({N}Rd)$:

$$\text{r}\widehat{\mathrm{MMD}}_{b}^{2}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R}):=V(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})=\Bigg{\|}\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(X_{i})-\frac{1}{{n_{2}}}\sum_{j=1}^{n_{2}}\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(Y_{j})\Bigg{\|}_{\mathbb{R}^{2R}}^{2}.$$

(4)

Notably, this estimator can be computed in linear time in terms of the pooled sample size ${N}$, and this computational benefit has motivated the prior work, such as Zhao and Meng, (2015), Sutherland and Schneider, (2015) and Cevid et al., (2022), that consider RFF-MMD statistics.

One may also consider an unbiased RFF-MMD statistic, given as

	$\displaystyle\text{r}\widehat{\mathrm{MMD}}_{u}^{2}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})$	$\displaystyle:=\frac{1}{{n_{1}}({n_{1}}-1)}\sum_{1\leq i\neq j\leq{n_{1}}}\langle\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(X_{i}),\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(X_{j})\rangle$		(5)
		$\displaystyle\quad+\frac{1}{{n_{2}}({n_{2}}-1)}\sum_{1\leq i\neq j\leq{n_{2}}}\langle\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(Y_{i}),\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(Y_{j})\rangle$
		$\displaystyle\quad-\frac{2}{{n_{1}}{n_{2}}}\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}\langle\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(X_{i}),\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(Y_{j})\rangle,$

which also involves $O({N}Rd)$ computational time (Zhao and Meng,, 2015, Appendix A.1). In this work, we consider both $\text{r}\widehat{\mathrm{MMD}}_{b}^{2}$ and $\text{r}\widehat{\mathrm{MMD}}_{u}^{2}$ to demonstrate statistical and computational trade-offs in RFF-based two-sample testing.

There have been several methods proposed for determining the threshold for MMD tests, which ensures (asymptotic or non-asymptotic) type I error control. These methods include those using limiting distributions or concentration inequalities, Gamma approximations, and bootstrap/permutation methods (e.g., Gretton et al., 2012a, ; Schrab et al.,, 2023). Among these, permutation tests stand out for their unique strength: they maintain level $\alpha$ for any finite sample size and often achieve optimal power (e.g., Kim et al.,, 2022). This advantage has made permutation tests a popular choice in real-world applications despite extra computational costs. Given their practical relevance, this work focuses on permutation-based MMD tests and establishes their theoretical guarantees.

To explain the procedure, let us write the pooled sample as $\mathcal{Z}_{N}:=\left\{Z_{1},\ldots,Z_{N}\right\}=\left\{\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}}\right\}$, and denote the collection of all possible permutations of $(1,2,\ldots,{N})$ as $\Pi_{N}$. Given a permutation $\pi:=(\pi(1),\ldots,\pi({N}))\in\Pi_{N},$ we denote the permuted pooled samples as $\mathcal{Z}^{\pi}_{N}:=\left\{Z_{\pi(1)},\ldots,Z_{\pi({N})}\right\}.$ Then, for a generic test statistic $T_{{n_{1}},{n_{2}}},$ the permutation distribution of $T_{{n_{1}},{n_{2}}}$ is defined as

$$F^{\pi}_{T_{{n_{1}},{n_{2}}}}(t):=\frac{1}{{N}!}\sum_{\pi\in\Pi_{{N}}}\mathds{1}\{T_{{n_{1}},{n_{2}}}(\mathcal{Z}^{\pi}_{N})\leq t\}.$$

The permutation test rejects the null hypothesis when $T_{{n_{1}},{n_{2}}}(\mathcal{Z}_{N})\geq q_{{n_{1}},{n_{2}},1-\alpha}$ where $q_{{n_{1}},{n_{2}},1-\alpha}$ denotes the $1-\alpha$ quantile of $F^{\pi}_{T_{{n_{1}},{n_{2}}}}$ given as

$$q_{{n_{1}},{n_{2}},1-\alpha}:=\inf\big{\{}t:F_{T_{{n_{1}},{n_{2}}}}^{\pi}(t)\geq 1-\alpha\big{\}}.$$

It is well-known that the resulting permutation test maintains non-asymptotic type I error control under the exchangeability of random vectors (e.g., Hemerik and Goeman,, 2018, Theorem 1). This exchangeability condition is satisfied under the null hypothesis of two-sample testing where $\mathcal{Z}_{N}$ are assumed to be i.i.d. random vectors.

A more computationally efficient permutation test is defined through Monte Carlo simulations. Let $\pi_{1},\ldots,\pi_{B}$ be permutation vectors randomly drawn from $\Pi_{N}$ with replacement. We let $T_{{n_{1}},{n_{2}}}^{(1)},\ldots,T_{{n_{1}},{n_{2}}}^{(B)}$ denote the test statistics computed based on $\mathcal{Z}_{N}^{\pi_{1}},\dots,\mathcal{Z}_{N}^{\pi_{B}}$. Let $\hat{q}_{{n_{1}},{n_{2}},1-\alpha}$ be the $1-\alpha$ quantile of the empirical distribution of $\{T_{{n_{1}},{n_{2}}},T_{{n_{1}},{n_{2}}}^{(1)},\ldots,T_{{n_{1}},{n_{2}}}^{(B)}\}$, and reject the null when $T_{{n_{1}},{n_{2}}}>\hat{q}_{{n_{1}},{n_{2}},1-\alpha}$. The resulting Monte Carlo-based test is also valid in finite samples (Hemerik and Goeman,, 2018, Theorem 2) and has almost equivalent power behavior as the full permutation test for sufficiently large $B$.

An alternative formulation of $\text{r}\widehat{\mathrm{MMD}}_{b}^{2}$ in Equation (4) is in terms of the characteristic functions of $P_{X}$ and $P_{Y}$. This reformulation provides a key insight into our negative result in Section 3.2. To fix ideas, the squared MMD with a translation-invariant kernel $k$ can be represented as $\mathrm{MMD}^{2}(P_{X},P_{Y};\mathcal{H}_{k})=\int_{\mathbb{R}^{d}}|\phi_{X}(\omega)-\phi_{Y}(\omega)|^{2}d\Lambda(\omega)$ where $\phi_{X}$ and $\phi_{Y}$ are the characteristic functions of $P_{X}$ and $P_{Y}$, respectively (e.g., Sriperumbudur et al.,, 2010, Corollary 4). Letting $\hat{\phi}_{X}(\omega):=\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}e^{\sqrt{-1}\omega^{\top}X_{i}}$ and $\hat{\phi}_{Y}(\omega):=\frac{1}{n_{2}}\sum_{j=1}^{n_{2}}e^{\sqrt{-1}\omega^{\top}Y_{j}}$, we may represent the plug-in estimator in Equation (2) as

$\displaystyle\widehat{\mathrm{MMD}}_{b}^{2}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\mathcal{H}_{k})=\int_{\mathbb{R}^{d}}|\hat{\phi}_{X}(\omega)-\hat{\phi}_{Y}(\omega)|^{2}d\Lambda(\omega).$

With this identity in place, the RFF-MMD statistic $\text{r}\widehat{\mathrm{MMD}}_{b}^{2}$ can be regarded as an approximation of the above plug-in estimator via Monte Carlo simulations with $R$ random frequencies $\{\omega_{r}\}_{r=1}^{R}\overset{\mathrm{i.i.d.}}{\sim}\Lambda$. Specifically, the RFF-MMD statistic can be written in terms of the empirical characteristic functions as:

$\displaystyle\text{r}\widehat{\mathrm{MMD}}_{b}^{2}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})=\frac{1}{R}\sum^{R}_{r=1}|\hat{\phi}_{X}(\omega_{r})-\hat{\phi}_{Y}(\omega_{r})|^{2}.$

As it is well-known, the characteristic function uniquely determines the distribution of a random vector. Therefore, when the support of $\Lambda$ is the entire Euclidean space, the population MMD becomes zero if and only if $P_{X}$ and $P_{Y}$ coincide. However, the empirical MMD evaluated on a finite number of random points is unable to capture an arbitrary difference between $P_{X}$ and $P_{Y}$, even asymptotically. At a high-level, this happens due to a combination of two factors. First of all, it is possible that two distinct characteristic functions can be equal in an interval (e.g., Romano and Siegel,, 1986, page 74). Moreover, if random evaluation points $\{\omega_{r}\}_{r=1}^{R}$ fall within such interval with high probability, then the RFF-MMD statistic would behave similarly to the null case, resulting in a test that is inconsistent with a fixed number of random features. This observation was partly made in Chwialkowski et al., (2015, Proposition 1), which we generalize to $\mathbb{R}^{d}$ as below.

The above lemma implies that there exists a certain pair of $(P_{X},P_{Y})$ under the alternative such that the expectation of the RFF-MMD statistic (4) is approximately zero with high probability. Given that the same test statistic has an expectation approximately equal to zero under the null, one may argue that the power of an RFF-MMD test would be strictly less than one against that specific alternative. However, this argument is insufficient to correctly claim the lack of consistency. An instructive example would be the case where a test statistic $W$ is either $0$ or $1/n$ with probability $1-\alpha$ and $\alpha$, respectively, under the null, whereas it takes the value $\alpha/n$ with probability one. In this case, it is clear to see that the expectation of $W$ remains the same under $H_{0}$ and $H_{1}$, converging to zero as $n\rightarrow\infty$. Nevertheless, if we reject the null when $W>0$, the resulting test has size $\alpha$ and power one for any value of $n\geq 1$. This toy example suggests that Lemma 1 is insufficient to formally prove the inconsistency result and we indeed need a distribution-level understanding of the RFF-MMD statistic. Moreover, when the critical value is determined via the permutation procedure (Section 2.4), we further need to take care of random sources arising from permutations, which adds an additional layer of technical challenges. With this context in place, we next develop inconsistency results by carefully studying the limiting distribution of the RFF-MMD statistic and its permuted counterpart.

Consider a permutation test based on the test statistic in Equation (4) defined as follows:

$$\Delta_{{n_{1},n_{2}},R}^{\alpha}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R}):=\Delta_{{n_{1},n_{2}},R}^{\alpha}:=\mathds{1}\big{\{}V(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})>q_{{n_{1},n_{2}},1-\alpha}\big{\}},$$

(6)

where $q_{{n_{1}},{n_{2}},1-\alpha}:=\inf\{t:F_{V}^{\pi}(t)\geq 1-\alpha\}$ and $F_{V}^{\pi}(t)=\frac{1}{{N}!}\sum_{\pi\in\Pi_{{N}}}\mathds{1}\{V(Z_{\pi(1)},\dots,Z_{\pi({N})};\boldsymbol{\omega}_{R})\leq t\}$. Building on the intuition laid out in Section 3.1, we aim to prove that the asymptotic power of the test $\Delta_{{n_{1}},{n_{2}},R}^{\alpha}$ is strictly less than one with a fixed number of $R$ against certain fixed alternatives. To formally achieve this, let $\boldsymbol{\psi}_{\boldsymbol{x}}$ be defined similarly as $\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}$ by replacing $\boldsymbol{\omega}_{R}$ with $\boldsymbol{x}\in\mathbb{R}^{d\times R}$. Based on Euler’s formula, the event $\frac{1}{R}\sum^{R}_{r=1}|\phi_{X}(\omega_{r})-\phi_{Y}(\omega_{r})|^{2}=0$ is equivalent to $\boldsymbol{\omega}_{R}\in\mathcal{E}:=\mathcal{E}(X,Y)$ where

$$\mathcal{E}(X,Y):=\bigl{\{}\boldsymbol{x}\in\mathbb{R}^{d\times R}:\mathbb{E}_{X}[\boldsymbol{\psi}_{\boldsymbol{x}}(X)]=\mathbb{E}_{Y}[\boldsymbol{\psi}_{\boldsymbol{x}}(Y)]\bigr{\}}.$$

(7)

We call $\boldsymbol{\omega}_{R}\in\mathcal{E}$ as the first moment equivalence (1-ME) condition, which holds with high probability, say $1-\epsilon$, for some fixed $(P_{X},P_{Y})$ according to Lemma 1. As mentioned earlier, the 1-ME condition alone is insufficient to formally prove the inconsistency result, which prompts an extension of the 1-ME condition to include higher-order moments. Specifically, consider a subset $\mathcal{E}_{k}\subseteq\mathcal{E}$ where $\boldsymbol{\omega}_{R}\in\mathcal{E}_{k}$ implies equivalence up to the $k$-th moment, i.e., with denoting $\boldsymbol{i}=(i_{1},\ldots,i_{2R})\in\mathbb{R}^{2R},$

$$\mathcal{E}_{k}:=\biggl{\{}\boldsymbol{x}\in\mathbb{R}^{d\times R}:\mathbb{E}_{X}[\boldsymbol{\psi}_{\boldsymbol{x}}(X)^{\boldsymbol{i}}]=\mathbb{E}_{Y}[\boldsymbol{\psi}_{\boldsymbol{x}}(Y)^{\boldsymbol{i}}],~{}\forall\boldsymbol{i}~{}\text{such that}~{}\{i_{r}\}^{2R}_{r=1}\in\mathbb{N}\cup\{0\},~{}\sum^{2R}_{r=1}i_{r}\leq k\biggr{\}},$$

for $\boldsymbol{\psi}_{\boldsymbol{x}}(X)^{\boldsymbol{i}}:=(\boldsymbol{\psi}_{\boldsymbol{x}}(X)_{1}^{i_{1}},\ldots,\boldsymbol{\psi}_{\boldsymbol{x}}(X)_{2R}^{i_{2R}})$ and $\boldsymbol{\psi}_{\boldsymbol{x}}(Y)^{\boldsymbol{i}}:=(\boldsymbol{\psi}_{\boldsymbol{x}}(Y)_{1}^{i_{1}},\ldots,\boldsymbol{\psi}_{\boldsymbol{x}}(Y)_{2R}^{i_{2R}})$. We refer $\boldsymbol{\omega}_{R}\in\mathcal{E}_{k}$ as the first $k$ moments equivalence (k-ME) condition. In the following proposition, we prove a generalized version of Lemma 1 demonstrating that the $k$-ME condition holds with high probability for some fixed $(P_{X},P_{Y})$. The proof of this result can be found in Section B.1.

Suppose that $P_{X},P_{Y}\in\mathcal{A}_{k,\epsilon}$, specified in Proposition 2. The power of the considered test against this specific alternative is then upper bounded as

	$\displaystyle\mathbb{P}(\Delta_{{n_{1}},{n_{2}},R}^{\alpha}=1)$	$\displaystyle=\int_{\mathcal{E}_{k}}\mathbb{P}(\Delta_{{n_{1}},{n_{2}},R}^{\alpha}=1\,|\,\boldsymbol{\omega}_{R}=\boldsymbol{\omega})f_{\boldsymbol{\omega}_{R}}(\boldsymbol{\omega})d\boldsymbol{\omega}+\int_{\mathcal{E}^{c}_{k}}\mathbb{P}(\Delta_{{n_{1}},{n_{2}},R}^{\alpha}=1\,|\,\boldsymbol{\omega}_{R}=\boldsymbol{\omega})f_{\boldsymbol{\omega}_{R}}(\boldsymbol{\omega})d\boldsymbol{\omega}$		(8)
		$\displaystyle\leq\int_{\mathcal{E}_{k}}\mathbb{P}(\Delta_{{n_{1}},{n_{2}},R}^{\alpha}=1\,|\,\boldsymbol{\omega}_{R}=\boldsymbol{\omega})f_{\boldsymbol{\omega}_{R}}(\boldsymbol{\omega})d\boldsymbol{\omega}+\epsilon.$

Given this bound, our proof for inconsistency revolves around showing that $\int_{\mathcal{E}_{k}}\mathbb{P}(\Delta_{{n_{1}},{n_{2}},R}^{\alpha}=1\,|\,\boldsymbol{\omega}_{R}=\boldsymbol{\omega})f_{\boldsymbol{\omega}_{R}}(\boldsymbol{\omega})d\boldsymbol{\omega}$ is sufficiently small. This in turn requires understanding the limiting behavior of the test statistic $V(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})$ and the permutation critical value $q_{{n_{1}},{n_{2}},1-\alpha}$ under the $k$-ME condition. On the one hand, the limiting distribution of the test statistic can be derived using the standard asymptotic tools such as the central limit theorem. On the other hand, we leverage asymptotic results for permutation distributions in Chung and Romano, (2016) to show that the critical value $q_{{n_{1}},{n_{2}},1-\alpha}$ converges to the $1-\alpha$ quantile of a continuous distribution. We point out that both limiting and permutation distributions of $V(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})$ are determined by the first two moments of $\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(X)$ and $\boldsymbol{\psi}_{\boldsymbol{\omega}_{R}}(Y)$. Furthermore, both distributions become asymptotically identical when those moments are the same, implying the coincidence of both distributions under the 2-ME condition. Consequently, the power of the test $\Delta_{{n_{1}},{n_{2}},R}^{\alpha}$ under the 2-ME condition remains small even asymptotically, which together with inequality (8), leads to the inconsistency result. This negative result is formally stated in the following theorem and the proof can be found in Section B.2.

The underlying idea of the proof for Theorem 3 can be applied to the unbiased RFF-MMD statistic in Equation (5) as well. In particular, consider a permutation test

$$\Delta_{{n_{1}},{n_{2}},R}^{\alpha,u}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R}):=\mathds{1}\big{\{}{U(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})}>q_{{n_{1}},{n_{2}},1-\alpha}^{u}\big{\}},$$

(9)

where $U(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R}):=\text{r}\widehat{\mathrm{MMD}}_{u}^{2}(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})$ and $q_{{n_{1}},{n_{2}},1-\alpha}^{u}:=\inf\{t:F_{U}^{\pi}(t)\geq 1-\alpha\}$ is the corresponding critical value. Building on the observation that the difference between $U(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})$ and $V(\mathcal{X}_{n_{1}},\mathcal{Y}_{n_{2}};\boldsymbol{\omega}_{R})$ is asymptotically negligible, we derive a result analogous to Theorem 3, demonstrating that $\Delta_{{n_{1}},{n_{2}},R}^{\alpha,u}$ fails to be pointwise consistent.

The proof of Corollary 4 can be found in Section B.3. Our findings so far indicate that RFF-MMD tests with a fixed number of random features fail to be pointwise consistent. To address this issue, we naturally consider increasing $R$ with the sample size and show that the tests then become pointwise consistent. Moreover, in some cases, the RFF-MMD test can attain comparable power to the quadratic-time MMD test but in strictly less than quadratic time. These are the topics of the next section.

We start by examining the uniform separation rate of the RFF-MMD test over the Sobolev ball with respect to the $L_{2}$ distance. Let us denote $\mathcal{S}_{d}^{s}(M_{1})$ as the $s$th order Sobolev ball in $\mathbb{R}^{d}$ with radius $M_{1}>0,$ that is

$$\mathcal{S}_{d}^{s}(M_{1}):=\Big{\{}f\in L^{1}\left(\mathbb{R}^{d}\right)\cap L^{2}\left(\mathbb{R}^{d}\right):\int_{\mathbb{R}^{d}}\|\omega\|_{2}^{2s}|\widehat{f}(\omega)|^{2}\mathrm{~{}d}\omega\leq(2\pi)^{d}M_{1}^{2}\Big{\}},$$

where $s>0$ is the smoothness parameter, and $\hat{f}(\omega)=\frac{1}{(2\pi)^{d/2}}\int_{\mathbb{R}^{2}}f(x)e^{-i\langle x,\omega\rangle}d\omega$ is the Fourier transform of $f$. Furthermore, each of $L^{1}(\mathbb{R}^{d})$ and $L^{2}(\mathbb{R}^{d})$ denotes a set of functions that are integrable in absolute value and square integrable, respectively. Let $\mathcal{P}_{\mathrm{conti}}$ be the collection of distribution pairs on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ where each pair of distributions $(P_{X},P_{Y})\in\mathcal{P}_{\mathrm{conti}}$ admits the probability density functions $(p_{X},p_{Y})$ with respect to the Lebesgue measure. Defining the class of distribution pairs with some constant $M_{2}>0$ given as

$\displaystyle\widetilde{\mathcal{C}}_{L_{2}}:=\big{\{}(P_{X},P_{Y})\in\mathcal{P}_{\mathrm{conti}}\,\big{|}\,p_{X}-p_{Y}\in\mathcal{S}_{d}^{s}(M_{1}),~{}\max\{\|p_{X}\|_{\infty},\|p_{Y}\|_{\infty}\}\leq M_{2}\big{\}},$

Schrab et al., (2023) demonstrated that the minimax rate in terms of the $L_{2}$ distance is $\rho^{\star}(\alpha,\beta,\widetilde{\mathcal{C}}_{L_{2}},\delta_{L_{2}})\asymp{n}^{-2s/(4s+d)}$, where $a_{n}\asymp b_{n}$ indicates $c\leq|a_{n}/b_{n}|\leq C$ for some positive constants $c,C$. They further showed that the MMD test using a translation-invariant kernel is minimax rate optimal in a non-asymptotic sense. A similar but asymptotic result was obtained by Li and Yuan, (2019), focusing specifically on the Gaussian kernel.¹¹1Both Schrab et al., (2023) and Li and Yuan, (2019) assume that $n_{1}\asymp n_{2}$ under which the minimax rate against the $L_{2}$ alternative is given as $(n_{1}+n_{2})^{-2s/(4s+d)}$. Without this balanced sample size assumption, however, the minimax rate is dominated by the minimum sample size i.e., $n^{-2s/(4s+d)}$. It is intuitively clear that when the number of random features $R$ is sufficiently large, the RFF-MMD test will also attain the same minimax optimality as the law of large numbers guarantees that the approximated kernel $\hat{k}$ converges to the underlying kernel $k$ almost surely. Our next question is then to ask how rapidly $R$ should be increased to ensure the same minimax guarantee and whether it is possible to attain the same optimality in sub-quadratic time. We answer these questions in the affirmative.

Similarly to Schrab et al., (2023), our analysis assumes that the kernel $k$ can be represented as a product of $d$ one-dimensional translation-invariant characteristic kernels with a given bandwidth $\lambda$. More specifically, we assume that the kernel $k$ can be decomposed as

$$k(x,y)=k_{\lambda}(x,y):=\prod^{d}_{i=1}\frac{1}{\lambda_{i}}\kappa_{i}\left(\frac{x_{i}-y_{i}}{\lambda_{i}}\right)$$

for $\lambda=(\lambda_{1},\ldots,\lambda_{d})^{\top}\in(0,\infty)^{d},$ where $\kappa_{i}:\mathbb{R}\rightarrow\mathbb{R}$ are some non-negative functions in $L^{1}\left(\mathbb{R}\right)\cap L^{2}\left(\mathbb{R}\right)$ satisfying $\int_{\mathbb{R}}\kappa_{i}(x)dx=1$ for $i=1,\ldots,d$. We note that $k_{\lambda}$ is indeed a characteristic kernel on $\mathbb{R}^{d}\times\mathbb{R}^{d}$, and we treat the bandwidth $\lambda$ as a tuning parameter that varies with the sample size. In order to highlight the dependence on $\lambda$, we let $\Delta_{{n_{1}},{n_{2}},R}^{\alpha,\lambda}$ (resp. $\Delta_{{n_{1}},{n_{2}},R}^{\alpha,u,\lambda}$) denote the test $\Delta_{{n_{1}},{n_{2}},R}^{\alpha}$ (resp. $\Delta_{{n_{1}},{n_{2}},R}^{\alpha,u}$) equipped with the kernel $k_{\lambda}$. Given a constant $M_{3}>0$, let us now consider a subset of $\widetilde{\mathcal{C}}_{L_{2}}$ where the support of individual distributions lies within the $d$-dimensional hypercube $[-M_{3},M_{3}]^{d}$. In other words, we define

$$\mathcal{C}_{L_{2}}:=\bigl{\{}(P_{X},P_{Y})\in\widetilde{\mathcal{C}}_{L_{2}}\,\big{|}\,\mathrm{support}(P_{X}),\,\mathrm{support}(P_{Y})\subset[-M_{3},M_{3}]^{d}\bigr{\}}.$$

Recalling $N=n_{1}+n_{2}$ and $n=\min\{n_{1},n_{2}\}$, the following theorem discusses the choice of $R$ and $\lambda$ that allows $\Delta_{{n_{1}},{n_{2}},R}^{\alpha,\lambda}$ and $\Delta_{{n_{1}},{n_{2}},R}^{\alpha,u,\lambda}$ to achieve the minimax separation rate against the class of alternatives defined on $\mathcal{C}_{L_{2}}$.