An approximate randomization test for high-dimensional two-sample Behrens-Fisher problem under arbitrary covariances

Two-sample mean testing is a fundamental problem in statistics with an enormous range of applications. In modern statistical applications, high-dimensional data, where the data dimension may be much larger than the sample size, is ubiquitous. However, most classical two-sample mean tests are designed for low-dimensional data, and may not be feasible, or may have suboptimal power, for high-dimensional data; see, e.g., Bai and Saranadasa, (1996). In recent years, the study of high-dimensional two-sample mean tests has attracted increasing attention.

Suppose that $X_{k,i}$, $i=1,\ldots,n_{k}$, $k=1,2$, are independent $p$-dimensional random vectors with $\operatorname{E}(X_{k,i})=\mu_{k}$, $\operatorname{var}(X_{k,i})=\bm{\Sigma}_{k,i}$. The hypothesis of interest is

$\displaystyle\mathcal{H}_{0}:\mu_{1}=\mu_{2}\quad\text{v.s.}\quad\mathcal{H}_{1}:\mu_{1}\neq\mu_{2}.$

(1)

Denote by $\bar{}\bm{\Sigma}_{k}=n_{k}^{-1}\sum_{i=1}^{n_{k}}\bm{\Sigma}_{k,i}$ the average covariance matrix within group $k$, $k=1,2$. Most existing methods on two-sample mean tests assumed that the observations within groups are identically distributed. In this case, $\bm{\Sigma}_{k,i}=\bar{}\bm{\Sigma}_{k}$, $i=1,\ldots,n_{k}$, $k=1,2$, and the considered testing problem reduces to the well-known Behrens-Fisher problem. In this paper, we consider the more general setting where the observations are allowed to have different distributions even within groups. Also, we consider the case where the data is possibly high-dimensional, that is, the data dimension $p$ may be much larger than the sample size $n_{k}$, $k=1,2$.

Let $\bar{X}_{k}=n_{k}^{-1}\sum_{i=1}^{n_{k}}X_{k,i}$ and $\mathbf{S}_{k}=(n_{k}-1)^{-1}\sum_{i=1}^{n_{k}}(X_{k,i}-\bar{X}_{k})(X_{k,i}-\bar{X}_{k})^{\intercal}$ denote the sample mean vector and the sample covariance matrix of group $k$, respectively, $k=1,2$. Denote $n=n_{1}+n_{2}$. A classical test statistic for hypothesis (1) is Hotelling’s $T^{2}$ statistic, defined as

$\displaystyle\frac{n_{1}n_{2}}{n}(\bar{X}_{1}-\bar{X}_{2})^{\intercal}\mathbf{S}^{-1}(\bar{X}_{1}-\bar{X}_{2}),$

where $\mathbf{S}=(n-2)^{-1}\{(n_{1}-1)\mathbf{S}_{1}+(n_{2}-1)\mathbf{S}_{2}\}$ is the pooled sample covariance matrix. When $p>n-2$, the matrix $\mathbf{S}$ is not invertible, and consequently Hotelling’s $T^{2}$ statistic is not well-defined. In a seminal work, Bai and Saranadasa, (1996) proposed the test statistic

$\displaystyle T_{\mathrm{BS}}(\mathbf{X}_{1},\mathbf{X}_{2})=\|\bar{X}_{1}-\bar{X}_{2}\|^{2}{-\frac{n}{n_{1}n_{2}}\operatorname{tr}(\mathbf{S})},$

where $\mathbf{X}_{k}=(X_{k,1},\ldots,X_{k,n_{k}})^{\intercal}$ is the data matrix of group $k$, $k=1,2$. The test statistic $T_{\mathrm{BS}}(\mathbf{X}_{1},\mathbf{X}_{2})$ is well-defined for arbitrary $p$. Bai and Saranadasa, (1996) assumed that the observations within groups are identically distributed and $\bar{}\bm{\Sigma}_{1}=\bar{}\bm{\Sigma}_{2}$. The main term of $T_{\mathrm{BS}}(\mathbf{X}_{1},\mathbf{X}_{2})$ is $\|\bar{X}_{1}-\bar{X}_{2}\|^{2}$. Chen and Qin, (2010) removed terms $X_{k,i}^{\intercal}X_{k,i}$ from $\|\bar{X}_{1}-\bar{X}_{2}\|^{2}$ and proposed the test statistic

$\displaystyle T_{\text{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})=\sum_{k=1}^{2}\frac{2\sum_{i=1}^{n_{k}}\sum_{j=i+1}^{n_{k}}X_{k,i}^{\intercal}X_{k,j}}{n_{k}(n_{k}-1)}-\frac{2\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}X_{1,i}^{\intercal}X_{2,j}}{n_{1}n_{2}}.$

The test method of Chen and Qin, (2010) can be applied when $\bar{}\bm{\Sigma}_{1}\neq\bar{}\bm{\Sigma}_{2}$. Both Bai and Saranadasa, (1996) and Chen and Qin, (2010) used the asymptotical normality of the statistics to determine the critical value of the test. However, the asymptotical normality only holds for a restricted class of covariance structures. For example, Chen and Qin, (2010) assumed that

$\displaystyle\operatorname{tr}(\bar{}\bm{\Sigma}_{k}^{4})=o\left[[\operatorname{tr}\{(\bar{}\bm{\Sigma}_{1}+\bar{}\bm{\Sigma}_{2})^{2}\}]^{2}\right],\quad k=1,2.$

(2)

However, (2) may not hold when $\bar{}\bm{\Sigma}_{k}$ has a few eigenvalues significantly larger than the others, excluding some important senarios in practice. For example, (2) is often violated when the variables are affected by some common factors; see, e.g., Fan et al., (2021) and the references therein. If the condition (2) is violated, the test procedures of Bai and Saranadasa, (1996) and Chen and Qin, (2010) may have incorrect test level; see, e.g., Wang and Xu, (2019). For years, how to construct test procedures that are valid under general covariances has been an important open problem in the field of high-dimensional hypothesis testing; see the review paper of Hu and Bai, (2016).

Intuitively, one may resort to the bootstrap method to control the test level under general covariances. Surprisingly, as shown by our numerical experiments, the empirical bootstrap method does not work well for $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$. Also, the wild bootstrap method, which is popular in high-dimensional statistics, does not work well either. We will give heuristic arguments to understand this phenomenon in Section 3. Recently, several methods were proposed to give a better control of the test level of $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ and related test statistics. In the setting of $\bar{}\bm{\Sigma}_{1}=\bar{}\bm{\Sigma}_{2}$, Zhang et al., (2020) considered the statistic $\|\bar{X}_{1}-\bar{X}_{2}\|^{2}$ and proposed to use the Welch-Satterthwaite $\chi^{2}$-approximation to determine the critical value. Subsequently, Zhang et al., (2021) extended the test procedure of Zhang et al., (2020) to the two-sample Behrens-Fisher problem, and Zhang and Zhu, (2021) considered a $\chi^{2}$-approximation of $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$. However, the test methods of Zhang et al., (2020), Zhang et al., (2021) and Zhang and Zhu, (2021) still require certain conditions on the eigenstructure of $\bar{}\bm{\Sigma}_{k}$, $k=1,2$. In another work, Wu et al., (2018) investigated the general distributional behavior of $\|\bar{X}_{1}\|^{2}$ for one-sample mean testing problem. They proposed a half-sampling procedure to determine the critical value of the test statistic. A generalization of the method of Wu et al., (2018) to the statistc $\|\bar{X}_{1}-\bar{X}_{2}\|^{2}$ for the two-sample Behrens-Fisher problem when both $n_{1}$ and $n_{2}$ are even numbers was presented in Chapter 2 of Lou, (2020). Unfortunately, Wu et al., (2018) and Lou, (2020) did not include detailed proofs of their theoretical results. Recently, Wang and Xu, (2019) considered a randomization test based on the test statistic of Chen and Qin, (2010) for the one-sample mean testing problem. Their randomization test has exact test level if the distributions of the observations are symmetric and has asymptotically correct level in the general setting. Although many efforts have been made, no existing test procedure based on $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ is valid for arbitrary $\bar{}\bm{\Sigma}_{k}$ with rigorous theoretical guarantees.

The statistics $T_{\mathrm{BS}}(\mathbf{X}_{1},\mathbf{X}_{2})$ and $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ are based on sum-of-squares. Some researchers investigated variants of sum-of-squares statistics that are scalar-invariant; see, e.g., Srivastava and Du, (2008), Srivastava et al., (2013) and Feng et al., (2015). These test methods also impose some nontrivial assumptions on the eigenstructure of the covariance matrices. There is another line of research initiated by Cai et al., (2013) which utilizes extreme value type statistics to test hypothesis (1). Chang et al., (2017) considered a parametric bootstrap method to determine the critical value of the extreme value type statistics. The method of Chang et al., (2017) allows for general covariance structures. Recently, Xue and Yao, (2020) investigated a wild bootstrap method. The test procedure in Xue and Yao, (2020) does not require that the observations within groups are identically distributed. The theoretical analyses in Chang et al., (2017) and Xue and Yao, (2020) are based on recent results about high-dimensional bootstrap methods developed by Chernozhukov et al., (2013) and Chernozhukov et al., (2017). While their theoretical framework can be applied to extreme value type statistics, it can not be applied to sum-of-squares type statistics like $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$. Indeed, we shall see that the wild bootstrap method does not work well for $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$.

In the present paper, we propose a new test procedure based on $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$. In the proposed test procedure, we use a new randomization method to determine the critical value of the test statistic. The proposed randomization method is motivated by the randomization method proposed in Fisher, (1935), Section 21. While this randomization method is widely used in one-sample mean testing problem, it can not be directly applied to two-sample Behrens-Fisher problem. In comparison, the proposed randomization method has satisfactory performance for two-sample Behrens-Fisher problem. We investigate the asymptotic properties of the proposed test procedure and rigorously prove that it has correct test level asymptotically under fairly weak conditions. Compared with existing test procedures based on sum-of-squares, the present work has the following favorable features. First, the proposed test procedure has correct test level asymptotically without any restriction on $\bar{}\bm{\Sigma}_{k}$, $k=1,2$. As a consequence, the proposed test procedure serves as a rigorous solution to the open problem that how to construct a valid test based on $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ under general covariances. Second, the proposed test procedure is valid even if the observations within groups are not identically distributed. To the best of our knowledge, the only existing method that can work in such a general setting is the test procedure in Xue and Yao, (2020). However, the test procedure in Xue and Yao, (2020) is based on extreme value statistic, which is different from the present work. Third, our theoretical results are valid for arbitrary $n_{1}$ and $n_{2}$ as long as $\min(n_{1},n_{2})\to\infty$. In comparison, all the above mentioned methods only considered the balanced data, that is, the sample sizes $n_{1}$ and $n_{2}$ are of the same order. We also derive the asymptotic local power of the proposed method. It shows that the asymptotic local power of the proposed test procedure is the same as that of the oracle test procedure where the distribution of the test statistic is known. A key tool for the proofs of our theoretical results is a universality property of the generalized quadratic forms. This result may be of independent interest. We also conduct numerical experiments to compare the proposed test procedure with some existing methods. Our theoretical and numerical results show that the proposed test procedure has good performance in both level and power.

In this section, we investigate the asymptotic behavior of $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ under general conditions. We begin with some notations. For two random variables $\xi$ and $\zeta$, let $\mathcal{L}(\xi)$ denote the distribution of $\xi$, and $\mathcal{L}(\xi\mid\zeta)$ denote the conditional distribution of $\xi$ given $\zeta$. For two probability measures $\nu_{1}$ and $\nu_{2}$ on $\mathbb{R}$, let $\nu_{1}-\nu_{2}$ be the signed measure such that for any Borel set $\mathcal{A}$, $(\nu_{1}-\nu_{2})(\mathcal{A})=\nu_{1}(\mathcal{A})-\nu_{2}(\mathcal{A})$. Let $\mathscr{C}_{b}^{3}(\mathbb{R})$ denote the class of bounded functions with bounded and continuous derivatives up to order $3$. It is known that a sequence of random variables $\{\xi_{n}\}_{n=1}^{\infty}$ converges weakly to a random variable $\xi$ if and only if for every $f\in\mathscr{C}_{b}^{3}(\mathbb{R})$, $\operatorname{E}f(\xi_{n})\to\operatorname{E}f(\xi)$; see, e.g., Pollard, (1984), Chapter III, Theorem 12. We use this property to give a metrization of the weak convergence in $\mathbb{R}$. For a function $f\in\mathscr{C}_{b}^{3}(\mathbb{R})$, let $f^{(i)}$ denote the $i$th derivative of $f$, $i=1,2,3$. For a finite signed measure $\nu$ on $\mathbb{R}$, we define the norm $\|\nu\|_{3}$ as $\sup_{f}\left|\int_{\mathbb{R}^{d}}f(\mathbf{x})\,\nu(\mathrm{d}\mathbf{x})\right|$, where the supremum is taken over all $f\in\mathscr{C}_{b}^{3}(\mathbb{R})$ such that $\sup_{\mathbf{x}\in\mathbb{R}}|f(x)|\leq 1$ and $\sup_{\mathbf{x}\in\mathbb{R}}|f^{(\ell)}(x)|\leq 1$, $\ell=1,2,3$. It is straightforward to verify that $\|\cdot\|_{3}$ is indeed a norm for finite signed measures. Also, a sequence of probability measures $\{\nu_{n}\}_{n=1}^{\infty}$ converges weakly to a probability measure $\nu$ if and only if $\|\nu_{n}-\nu\|_{3}\to 0$; see, e.g., Dudley, (2002), Corollary 11.3.4.

The test procedure of Chen and Qin, (2010) determines the critical value based on the asymptotic normality of $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ under the null hypothesis. Define $Y_{k,i}=X_{k,i}-\mu_{k}$, $i=1,\ldots,n_{k}$, $k=1,2$. Then $\operatorname{E}(Y_{k,i})=\mathbf{0}_{p}$. Under the null hypothesis, we have $T_{\text{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})=T_{\text{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})$ where $\mathbf{Y}_{k}=(Y_{k,1},\ldots,Y_{k,n_{k}})^{\intercal}$, $k=1,2$. Denote by $\sigma_{T,n}^{2}$ the variance of $T_{\text{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})$. Under certain conditions, Chen and Qin, (2010) proved that $T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})/\sigma_{T,n}$ converges weakly to $\mathcal{N}(0,1)$. They reject the null hypothesis if $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})/\hat{\sigma}_{T,n}$ is greater than the $1-\alpha$ quantile of $\mathcal{N}(0,1)$ where $\hat{\sigma}_{T,n}$ is an estimate of $\sigma_{T,n}$ and $\alpha\in(0,1)$ is the test level. However, we shall see that normal distribution is just one of the possible asymptotic distributions of $T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})$ in the general setting.

Now we list the main conditions imposed by Chen and Qin, (2010) to prove the asymptotic normality of $T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})$. First, Chen and Qin, (2010) imposed the condition (2) on the eigenstructure of $\bar{}\bm{\Sigma}_{k}$, $k=1,2$. Second, Chen and Qin, (2010) assumed that as $n\to\infty$,

$\displaystyle{n_{1}}/{n}\to b\in(0,1).$

(3)

That is, Chen and Qin, (2010) assumed the sample sizes in two groups are balanced. This condition is commonly adopted by existing test procedures for hypothesis (1). Third, Chen and Qin, (2010) assumed the general multivariate model

$\displaystyle Y_{k,i}=\bm{\Gamma}_{k}Z_{k,i},\quad i=1,\ldots,n_{k},\quad k=1,2,$

(4)

where $\bm{\Gamma}_{k}$ is a $p\times m$ matrix for some $m\geq p$, $k=1,2$, and $\{Z_{k,i}\}_{i=1}^{n_{k}}$ are $m$-dimensional independent and identically distributed random vectors such that $\operatorname{E}(Z_{k,i})=\mathbf{0}_{m}$, $\operatorname{var}(Z_{k,i})=\mathbf{I}_{m}$, and the elements of $Z_{k,i}=(z_{k,i,1},\ldots,z_{k,i,m})^{\intercal}$ satisfy $\operatorname{E}(z_{k,i,j}^{4})=3+\Delta<\infty$, and $\operatorname{E}(z_{k,i,\ell_{1}}^{\alpha_{1}}z_{k,i,\ell_{2}}^{\alpha_{2}}\cdots z_{k,i,\ell_{q}}^{\alpha_{q}})=\operatorname{E}(z_{k,i,\ell_{1}}^{\alpha_{1}})\operatorname{E}(z_{k,i,\ell_{2}}^{\alpha_{2}})\cdots\operatorname{E}(z_{k,i,\ell_{q}}^{\alpha_{q}})$ for a positive integer $q$ such that $\sum_{\ell=1}^{q}\alpha_{\ell}\leq 8$ and distinct $\ell_{1},\ldots,\ell_{q}\in\{1,\ldots,m\}$. We note that the general multivariate model (4) assumes that the observations within group $k$ are identically distributed, $k=1,2$. Also, the $m$ elements $z_{k,i,1},\ldots,z_{k,i,m}$ of $Z_{k,i}$ have finite moments up to order $8$ and behave as if they are independent.

As we have noted, the condition (2) can be violated when the variables are affected by some common factors. The condition (3) is assumed by most existing test methods. However, this condition is not reasonable if the sample sizes are unbalanced. The general multivariate model (4) is commonly adopted by many high-dimensional test methods; see, e.g., Chen and Qin, (2010), Zhang et al., (2020) and Zhang et al., (2021). However, the conditions on $Z_{k,i}$ may be difficult to verify and are not generally satisfied by the elliptical symmetric distributions. Also, the model (4) is not valid if the observations within groups are not identically distributed. Thus, we would like to investigate the asymptotic behavior of $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ beyond the conditions (2), (3) and (4). We consider the asymptotic setting where $n\to\infty$ and all quantities except for absolute constants are indexed by $n$, a subscript we often suppress. We make the following assumption on $n_{1}$ and $n_{2}$.

Assumption 1 only requires that both $n_{1}$ and $n_{2}$ tend to infinity, which allows for the unbalanced sample sizes. This relaxes the condition (3). We make the following assumption on the distributions of $Y_{k,i}$, $i=1,\ldots,n_{k}$, $k=1,2$.

Intuitively, Assumption 2 requires that the fourth moments of $Y_{k,i}$ are of the same order as the squared second moments of $Y_{k,i}$. We shall see that this assumption is fairly weak. However, the above inequality is required to be satisfied for all positive semi-definite matrix $\mathbf{B}$, which may not be straightforward to check in some cases. The following lemma gives a sufficient condition of Assumption 2.

It can be seen that the conditions of Lemma 1 is strictly weaker than the multivariate model (4). In fact, Lemma 1 does not require $m_{k,i}\geq p$, nor does it require the finiteness of $8$th moments of $z_{k,i,j}$. Also, the moment conditions in (5) are much weaker than that required by the multivariate model (4). In addition to the multivariate model (4), the conditions of Lemma 1 also allow for an important class of elliptical symmetric distributions. In fact, if $Y_{k,i}$ has an elliptical symmetric distribution, then $Y_{k,i}$ can be written as $Y_{k,i}=\eta_{k,i}\bm{\Gamma}_{k,i}U_{k,i}$, where $\bm{\Gamma}_{k,i}$ is a $p\times m_{k,i}$ matrix, $U_{k,i}$ is a random vector distributed uniformly on the unit sphere in $\mathbb{R}^{m_{k,i}}$, and $\eta_{k,i}$ is a nonnegative random variable which is independent of $U_{k,i}$; see, e.g., Fang et al., (1990), Theorem 2.5. Suppose there is an absolute constant $C$ such that $\operatorname{E}(\eta_{k,i}^{4})\leq C\{\operatorname{E}(\eta_{k,i}^{2})\}^{2}<\infty$. Then from the symmetric property of $U_{k,i}$ and the independence of $\eta_{k,i}$ and $U_{k,i}$, the conditions of Lemma 1 hold. In comparison, the multivariate model (4) does not allow for elliptical symmetric distributions in general.

Most existing test methods for the hypothesis (1) assume that the observations within groups are identically distributed. In this case, Assumptions 1 and 2 are all we need, and we completely avoid an assumption on the eigenstucture of $\bar{}\bm{\Sigma}_{k}$ like (2). In the general setting, $\bm{\Sigma}_{k,i}$ may be different within groups, and we make the following assumption to avoid the case in which there exist certain observations with significantly larger variance than the others.

If the covariance matrices within groups are equal, i.e., $\bm{\Sigma}_{k,i}=\bar{}\bm{\Sigma}_{k}$, $i=1,\ldots,n_{k}$, $k=1,2$, then Assumption 3 holds for arbitrary $\bar{}\bm{\Sigma}_{k}$, $k=1,2$, provided $\min(n_{1},n_{2})\to\infty$ as $n\to\infty$. In this view, Assumption 3 imposes no restriction on $\bar{}\bm{\Sigma}_{k}$, $k=1,2$.

Define $\bm{\Psi}_{n}=n_{1}^{-1}\bar{}\bm{\Sigma}_{1}+n_{2}^{-1}\bar{}\bm{\Sigma}_{2}$. Let $\bm{\xi}_{p}$ be a $p$-dimensional standard normal random vector. We have the following theorem.

Theorem 1 characterizes the general distributional behavior of $T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})$. It implies that the distributions of $T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})/\sigma_{T,n}$ and ${\{\bm{\xi}_{p}^{\intercal}\bm{\Psi}_{n}\bm{\xi}_{p}-\operatorname{tr}(\bm{\Psi}_{n})\}}/{\{2\operatorname{tr}(\bm{\Psi}_{n}^{2})\}^{1/2}}$ are equivalent asymptotically. To gain further insights on the distributional behavior of $T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})$, we would like to derive the asymptotic distributions of $T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})/\sigma_{T,n}$. However, Theorem 1 implies that $T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})/\sigma_{T,n}$ may not converge weakly in general. Nevertheless, $\mathcal{L}\{T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})/\sigma_{T,n}\}$ is uniformly tight and we can use Theorem 1 to derive all possible asympotitc distributions of $T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})/\sigma_{T,n}$.

Corollary 1 gives a full characteristic of the possible asymptotic distributions of $T_{\mathrm{CQ}}(\mathbf{Y}_{1},\mathbf{Y}_{2})/\sigma_{T,n}$. In general, these possible asymptotic distributions are weighted sums of an independent normal random variable and centered $\chi^{2}$ random variables. From the proof of Corollary 1, the parameters $\{\kappa_{i}\}_{i=1}^{\infty}$ are in fact the limit of the eigenvalues of the matrix $\bm{\Psi}_{n}/\{\operatorname{tr}(\bm{\Psi}_{n}^{2})\}^{1/2}$ along a subsequence of $\{n\}$. If $\sum_{i=1}^{\infty}\kappa_{i}^{2}=0$, then (6) becomes the standard normal distribution, and the test procedure of Chen and Qin, (2010) has correct level asymptotically. On the other hand, if $\sum_{i=1}^{\infty}\kappa_{i}^{2}=1$, then (6) becomes the distribution of a weighted sum of independent centered $\chi^{2}$ random variables. This case was considered in Zhang et al., (2021). However, these two settings are just special cases among all possible asymptotic distributions where $\sum_{i=1}^{\infty}\kappa_{i}^{2}\in[0,1]$. In general, the distribution (6) relies on the nuisance paramters $\{\kappa_{i}\}_{i=1}^{\infty}$. To construct a test procedure based on the asymptotic distributions, one needs to estimate these nuisance parameters consistently. Unfortunately, the estimation of the eigenvalues of high-dimensional covariance matrices may be a highly nontrivial task; see, e.g., Kong and Valiant, (2017) and the references therein. Hence in general, it may not be a good choice to construct test procedures based on the asymptotic distributions.

An intuitive idea to control the test level of $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ in the general setting is to use the bootstrap method. Surprisingly, the empirical bootstrap method and wild bootstrap method may not work for $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$. This phenomenon will be shown by our numerical experiments. For now, we give a heuristic argument to understand this phenomenon. First we consider the empirical bootstrap method. Suppose the resampled observations $\{X_{k,i}^{*}\}_{i=1}^{n_{k}}$ are uniformly sampled from $\{X_{k,i}-\bar{X}_{k}\}_{i=1}^{n_{k}}$ with replacement, $k=1,2$. Denote $\mathbf{X}_{k}^{*}=(X_{k,1}^{*},\ldots,X_{k,n_{k}}^{*})^{\intercal}$, $k=1,2$. The empirical bootstrap method uses the conditional distribution $\mathcal{L}\{T_{\mathrm{CQ}}(\mathbf{X}_{1}^{*},\mathbf{X}_{2}^{*})\mid\mathbf{X}_{1},\mathbf{X}_{2}\}$ to approximate the null distribution of $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$. If this bootstrap method can work, one may expect that the first two moments of $\mathcal{L}\{T_{\mathrm{CQ}}(\mathbf{X}_{1}^{*},\mathbf{X}_{2}^{*})\mid\mathbf{X}_{1},\mathbf{X}_{2}\}$ can approximately match the first two moments of $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ under the null hypothesis. Under the null hypothesis, we have $\operatorname{E}\{T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})\}=0$. Lemma S.5 implies that under Assumptions 1 and 3 and under the null hypothesis,

$\displaystyle\operatorname{var}\{T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})\}=\sigma_{T,n}^{2}=\{1+o(1)\}2\operatorname{tr}(\bm{\Psi}_{n}^{2}).$

On the othen hand, it is straightforward to show that $\operatorname{E}\{T_{\mathrm{CQ}}(\mathbf{X}_{1}^{*},\mathbf{X}_{2}^{*})\mid\mathbf{X}_{1},\mathbf{X}_{2}\}=0$. Also, under Assumptions 1 and 3,

$\displaystyle\operatorname{var}\{T_{\mathrm{CQ}}(\mathbf{X}_{1}^{*},\mathbf{X}_{2}^{*})\mid\mathbf{X}_{1},\mathbf{X}_{2}\}=\{1+o_{p}(1)\}2\operatorname{tr}\{(n_{1}^{-1}\mathbf{S}_{1}+n_{2}^{-1}\mathbf{S}_{2})^{2}\}.$

Unfortunately, $\operatorname{tr}\{(n_{1}^{-1}\mathbf{S}_{1}+n_{2}^{-1}\mathbf{S}_{2})^{2}\}$ is not a ratio-consistent estimator of $\operatorname{tr}(\bm{\Psi}_{n}^{2})$ even in the settings where $X_{k,i}$ is normally distributed, the covariance matrices $\bm{\Sigma}_{k,i}$, $i=1,\ldots,n_{k}$, $k=1,2$, are all equal and $n_{1}=n_{2}$; see, e.g., Bai and Saranadasa, (1996) and Zhou and Guo, (2017). Consequently, the empirical bootstrap method may not be valid for $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$.

We turn to the wild bootstrap method. Recently, the wild bootstrap method has been widely used for extreme value type statistics in the high-dimensional setting; see, e.g., Chernozhukov et al., (2013), Chernozhukov et al., (2017), Xue and Yao, (2020) and Deng and Zhang, (2020). For the wild bootstrap method, the resampled observations are defined as $X_{k,i}^{*}=\varepsilon_{k,i}(X_{k,i}-\bar{X}_{k})$, $i=1,\ldots,n_{k}$, $k=1,2$, where $\{\varepsilon_{k,i}\}$ are independent and identically distributed random variables with $\operatorname{E}(\varepsilon_{k,i})=0$ and $\operatorname{var}(\varepsilon_{k,i})=1$, and are independent of the original data $\mathbf{X}_{1}$, $\mathbf{X}_{2}$. We have $\operatorname{E}\{T_{\mathrm{CQ}}(\mathbf{X}_{1}^{*},\mathbf{X}_{2}^{*})\mid\mathbf{X}_{1},\mathbf{X}_{2}\}=0$. With some tedious but straightforward derivations, it can be seen that under Assumptions 1 and 3,

$\displaystyle\operatorname{E}\{\operatorname{var}(T_{\mathrm{CQ}}(\mathbf{X}_{1}^{*},\mathbf{X}_{2}^{*})\mid\mathbf{X}_{1},\mathbf{X}_{2})\}=$

$\displaystyle\{1+o(1)\}2\operatorname{tr}(\bm{\Psi}_{n}^{2})+b,$

where $b$ is the bias term and satisfies

	$\displaystyle b=$	$\displaystyle\{1+o(1)\}\sum_{k=1}^{2}\sum_{i=1}^{n_{k}}\frac{4}{n_{k}^{5}}\operatorname{E}\{(Y_{k,i}^{\intercal}Y_{k,i})^{2}\}-\{1+o(1)\}\sum_{k=1}^{2}\frac{2}{n_{k}^{4}}\{\operatorname{tr}(\bar{}\bm{\Sigma}_{k})\}^{2}$
	$\displaystyle\geq$	$\displaystyle\{1+o(1)\}\sum_{k=1}^{2}\sum_{i=1}^{n_{k}}\frac{2}{n_{k}^{5}}\{\operatorname{tr}(\bm{\Sigma}_{k,i})\}^{2}.$

The above inequality implies that the bias term $b$ may not be negligible compared with $2\operatorname{tr}(\bm{\Psi}_{n}^{2})$. For example, if $\bm{\Sigma}_{k,i}=\mathbf{I}_{p}$, $i=1,\ldots,n_{k}$, $k=1,2$, then we have $2\operatorname{tr}(\bm{\Psi}_{n}^{2})=2(n_{1}^{-1}+n_{2}^{-1})^{2}p$ and $b\geq\{1+o(1)\}2(n_{1}^{-4}+n_{2}^{-4})p^{2}$. In this case, the bias term $b$ is not negligible provided $n^{2}=O(p)$. Thus, the wild bootstrap method may not be valid for $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ either.

We have seen that the methods based on asymptotic distributions and the bootstrap methods may not work well for the test statistic $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$. These phenomenons imply that it is highly nontrivial to construct a valid test procedure based on $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$. To construct a valid test procedure, we resort to the idea of the randomization test, a powerful tool in statistical hypothesis testing. The randomization test is an old idea and can at least date back to Fisher, (1935), Section 21; see Hoeffding, (1952), Lehmann and Romano, (2005), Section 15.2, Zhu, (2005) and Hemerik and Goeman, (2018) for general frameworks and extensions of the randomization tests. The original randomization method considered in Fisher, (1935), Section 21 can be abstracted into the following general form. Suppose $\xi_{1},\ldots,\xi_{n}$ are independent $p$-dimensional random vectors such that $\mathcal{L}(\xi_{i})=\mathcal{L}(-\xi_{i})$. Let $T(\xi_{1},\ldots,\xi_{n})$ be any statistics taking values in $\mathbb{R}$. Suppose $\epsilon_{1},\ldots,\epsilon_{n}$ are independent and identically distributed Rademacher random variables, i.e., $\mathrm{pr}\,(\epsilon_{i}=1)=\mathrm{pr}\,(\epsilon_{i}=-1)=1/2$, and are independent of $\xi_{1},\ldots,\xi_{n}$. Define the conditional cumulative distribution function $\hat{F}(\cdot)$ as $\hat{F}(x)=\mathrm{pr}\,\{T(\epsilon_{1}\xi_{1},\ldots,\epsilon_{n}\xi_{n})\leq x\mid\xi_{1},\ldots,\xi_{n}\}$. Then from the theory of randomization test (see, e.g., Lehmann and Romano, (2005), Section 15.2), for any $\alpha\in(0,1)$,

$\displaystyle\mathrm{pr}\,\left\{T(\xi_{1},\ldots,\xi_{n})>\hat{F}^{-1}(1-\alpha)\right\}\leq\alpha,$

where for any right continuous cumulative distribution function $F(\cdot)$ on $\mathbb{R}$, $F^{-1}(q)=\min\{{x\in\mathbb{R}}:F(x)\geq q\}$ for $q\in(0,1)$. Also, under mild conditions, the difference between the above probability and $\alpha$ is negligible.

To apply Fisher’s randomization method to specific problems, the key is to construct random variables $\xi_{i}$ such that $\mathcal{L}(\xi_{i})=\mathcal{L}(-\xi_{i})$ under the null hypothesis. This randomization method can be directly applied to one-sample mean testing problem, as did in Wang and Xu, (2019). However, it can not be readily applied to the testing problem (1). In fact, the mean vector $\mu_{k}$ of $X_{k,i}$ is unknown under the null hypothesis, and consequently, one can not expect that $\mathcal{L}(X_{k,i})=\mathcal{L}(-X_{k,i})$ holds under the null hypothesis. As a result, $X_{k,i}$ can not serve as $\xi_{i}$ in Fisher’s randomization method.

We observe that the difference $X_{k,i}-X_{k,i+1}$ has zero means and hence is free of the mean vector $\mu_{k}$. Also, if $\mathcal{L}(X_{k,i})=\mathcal{L}(X_{k,i+1})$, then $\mathcal{L}(X_{k,i}-X_{k,i+1})=\mathcal{L}(X_{k,i+1}-X_{k,i})$. These facts imply that Fisher’s randomization method may be applied to the random vectors $\tilde{X}_{k,i}=(X_{k,2i}-X_{k,2i-1})/2$, $i=1,\ldots,m_{k}$, $k=1,2$ where $m_{k}=\lfloor n_{k}/2\rfloor$, $k=1,2$. Define

$\displaystyle T_{\mathrm{CQ}}(\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})=$

$\displaystyle\sum_{k=1}^{2}\frac{2\sum_{i=1}^{m_{k}}\sum_{j=i+1}^{m_{k}}\tilde{X}_{k,i}^{\intercal}\tilde{X}_{k,j}}{m_{k}(m_{k}-1)}-\frac{2\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}\tilde{X}_{1,i}^{\intercal}\tilde{X}_{2,j}}{m_{1}m_{2}},$

where $\tilde{\mathbf{X}}_{k}=(\tilde{X}_{k,1},\ldots,\tilde{X}_{k,m_{k}})^{\intercal}$, $k=1,2$. Let $E=(\epsilon_{1,1},\ldots,\epsilon_{1,m_{1}},\epsilon_{2,1},\ldots,\epsilon_{2,m_{2}})^{\intercal}$ where $\{\epsilon_{k,i}\}$ are independent and identically distributed Rademacher random variables which are independent of $\tilde{\mathbf{X}}_{1}$, $\tilde{\mathbf{X}}_{2}$. Define randomized statistic

$\displaystyle T_{\mathrm{CQ}}(E;\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})=$

$\displaystyle\sum_{k=1}^{2}\frac{2\sum_{i=1}^{m_{k}}\sum_{j=i+1}^{m_{k}}\epsilon_{k,i}\epsilon_{k,j}\tilde{X}_{k,i}^{\intercal}\tilde{X}_{k,j}}{m_{k}(m_{k}-1)}-\frac{2\sum_{i=1}^{m_{1}}\sum_{j=1}^{m_{2}}\epsilon_{1,i}\epsilon_{2,j}\tilde{X}_{1,i}^{\intercal}\tilde{X}_{2,j}}{m_{1}m_{2}}.$

Define the conditional distribution function $\hat{F}_{\mathrm{CQ}}(x)=\mathrm{pr}\,\{T_{\mathrm{CQ}}(E;\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})\leq x\mid\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2}\}$. From Fisher’s randomization method, if $\mathcal{L}(\tilde{X}_{k,i})=\mathcal{L}(-\tilde{X}_{k,i})$, then for $\alpha\in(0,1)$, $\mathrm{pr}\,\{T_{\mathrm{CQ}}(\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})>\hat{F}_{\mathrm{CQ}}^{-1}(1-\alpha)\}\leq\alpha$. It can be seen that $T_{\mathrm{CQ}}(\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})$ and $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$ take a similar form. Also, under the null hypothesis, $\operatorname{E}\{T_{\mathrm{CQ}}(\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})\}=\operatorname{E}\{T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})\}=0$ and $\operatorname{var}\{T_{\mathrm{CQ}}(\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})\}=\{1+o(1)\}\operatorname{var}\{T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})\}$. Thus, it may be expected that $\mathcal{L}\{T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})\}\approx\mathcal{L}\{T_{\mathrm{CQ}}(\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})\}$ under the null hypothesis. On the other hand, the classical results of Hoeffding, (1952) on randomization tests give the insight that the randomness of $\hat{F}_{\mathrm{CQ}}^{-1}(1-\alpha)$ may be negligible for large samples. From the above insights, it can be expected that under the null hypothesis, for $\alpha\in(0,1)$,

$\displaystyle\mathrm{pr}\,\{T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})>\hat{F}_{\mathrm{CQ}}^{-1}(1-\alpha)\}\approx\mathrm{pr}\,\{T_{\mathrm{CQ}}(\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})>\hat{F}_{\mathrm{CQ}}^{-1}(1-\alpha)\}\approx\alpha.$

Motivated by the above heuristics, we propose a new test procedure which rejects the null hypothesis if

$\displaystyle T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})>\hat{F}_{\mathrm{CQ}}^{-1}(1-\alpha).$

While the assumption $\mathcal{L}(\tilde{X}_{k,i})=\mathcal{L}(-\tilde{X}_{k,i})$ is used in the above heuristic arguments, it will not be assumed in our theoretical analysis. This generality may not be surprising. Indeed, for low-dimensional testing problems, it is known that such symmetry conditions can often be relaxed for randomization tests; see, e.g., Romano, (1990), Chung and Romano, (2013) and Canay et al., (2017). In the proposed procedure, the conditional distribution $\mathcal{L}\{T_{\mathrm{CQ}}(E;\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})\mid\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2}\}$ is used to approximate the null distribution of $T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})$. We have $\operatorname{E}\{T_{\mathrm{CQ}}(E;\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})\mid\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2}\}=0$. From Lemma S.10, under Assumptions 1-3, we have

$\displaystyle\operatorname{var}\{T_{\mathrm{CQ}}(E;\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})\mid\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2}\}=\{1+o_{P}(1)\}\sigma_{T,n}^{2}.$

That is, the first two moments of $\mathcal{L}\{T_{\mathrm{CQ}}(E;\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2})\mid\tilde{\mathbf{X}}_{1},\tilde{\mathbf{X}}_{2}\}$ can match those of $\mathcal{L}\{T_{\mathrm{CQ}}(\mathbf{X}_{1},\mathbf{X}_{2})\}$. As we have seen, this favorable property is not shared by the empirical bootstrap method and wild bootstrap method.