Two-Sample Quadratic Tests with Unequal Covariances

Why Unequal Covariances Change the Test

A two-sample mean test asks whether two populations have the same mean vector. In high dimension, the squared distance between sample means is a natural starting point, but its null variability depends on both populations’ covariance matrices. If those matrices differ, replacing them with a pooled covariance or using a one-sample variance formula can give the wrong calibration.

Let $X_1,\ldots,X_{n_1}$ be independent observations from a population with mean $\mu_1$ and covariance matrix $\Sigma_1$, and let $Y_1,\ldots,Y_{n_2}$ be independent observations from a second population with mean $\mu_2$ and covariance matrix $\Sigma_2$. Assume the two samples are independent, and that both covariance matrices are finite. The target is $H_0:\mu_1=\mu_2$ against a difference in any direction. The quadratic statistic measures the full squared mean difference, so it is naturally suited to alternatives whose signal is spread across multiple coordinates.

Definition: Write $\bar X=n_1^{-1}\sum_i X_i$, $\bar Y=n_2^{-1}\sum_j Y_j$, and let $S_1,S_2$ be the usual sample covariance matrices with denominators $n_1-1,n_2-1$. The unbiased two-sample quadratic statistic is $$ T=\|\bar X-\bar Y\|^2-\frac{\operatorname{tr}(S_1)}{n_1}-\frac{\operatorname{tr}(S_2)}{n_2}. $$ It estimates $\|\mu_1-\mu_2\|^2$. Equivalently, it is the sum of the two within-sample cross-products and the between-sample cross-product, with their unbiased coefficients.

T= \frac{1}{n_1(n_1-1)}\sum_{i\ne k}X_i^\top X_k+ \frac{1}{n_2(n_2-1)}\sum_{j\ne l}Y_j^\top Y_l- \frac{2}{n_1n_2}\sum_{i,j}X_i^\top Y_j .

This form shows why the correction terms matter: the squared norm of a sample mean includes each observation’s squared length, which contributes sampling noise even when the population mean is zero. Removing the sample covariance trace terms makes the statistic unbiased for the squared population difference, not merely a noisy estimate of it.

Mean and Null Variance

Set $\delta=\mu_1-\mu_2$. Independence gives $E[\|\bar X-\bar Y\|^2]=\|\delta\|^2+\operatorname{tr}(\Sigma_1)/n_1+\operatorname{tr}(\Sigma_2)/n_2$, while $E[\operatorname{tr}(S_g)]=\operatorname{tr}(\Sigma_g)$ in each group. Therefore,

E[T]=\|\delta\|^2.

In particular, $T$ has mean zero under the null. Its variance under the null depends on three separate trace quantities: the squared covariance size in each group and the cross-trace between the groups. Define $V_0=\operatorname{Var}(T\mid H_0)$.

Theorem: If the observations are independent within and between groups, have finite covariance matrices, and $H_0$ holds, then $$ V_0=\frac{2\operatorname{tr}(\Sigma_1^2)}{n_1(n_1-1)} +\frac{2\operatorname{tr}(\Sigma_2^2)}{n_2(n_2-1)} +\frac{4\operatorname{tr}(\Sigma_1\Sigma_2)}{n_1n_2}. $$

Proof. Under the null, write $e_i=X_i-\mu_1$ and $f_j=Y_j-\mu_2$. The cross-product representation becomes

T=U_1+U_2-2\bar e^\top\bar f,\qquad U_1=\frac{\sum_{i\ne k}e_i^\top e_k}{n_1(n_1-1)},\quad U_2=\frac{\sum_{j\ne l}f_j^\top f_l}{n_2(n_2-1)}.

In $\operatorname{Var}(U_1)$, terms $e_i^\top e_k$ and $e_a^\top e_b$ have zero covariance unless their unordered index pairs agree: if an index appears only once, conditioning on all other observations leaves a factor with mean zero. For any distinct $i,k$,

E[(e_i^\top e_k)^2] =E[e_i^\top\Sigma_1 e_i] =\operatorname{tr}(\Sigma_1^2).

There are $n_1(n_1-1)$ ordered pairs and each unordered pair occurs twice. Thus the sum’s variance is $2n_1(n_1-1)\operatorname{tr}(\Sigma_1^2)$; dividing by the squared denominator gives $\operatorname{Var}(U_1)=2\operatorname{tr}(\Sigma_1^2)/[n_1(n_1-1)]$. The same argument gives $\operatorname{Var}(U_2)=2\operatorname{tr}(\Sigma_2^2)/[n_2(n_2-1)]$.

The two sample means are independent, centered, and have covariances $\Sigma_1/n_1$ and $\Sigma_2/n_2$. Consequently,

\operatorname{Var}(\bar e^\top\bar f) =E[\bar e^\top(\Sigma_2/n_2)\bar e] =\frac{\operatorname{tr}(\Sigma_1\Sigma_2)}{n_1n_2}.

The three terms have zero pairwise covariance. For example, conditioning on the first sample makes the conditional expectation of $\bar e^\top\bar f$ zero; conditioning on the second sample gives the corresponding result for its covariance with $U_1$. The other cases follow in the same way, using independence and zero means. Adding the three variances, including the factor $4$ from $-2\bar e^\top\bar f$, proves the formula. $\square$

Worked Example: Computing the Null Variance With Different Covariances

Suppose $p=2$, $n_1=20$, $n_2=25$, and the covariance matrices are known to be $\Sigma_1=\operatorname{diag}(2,1)$ and $\Sigma_2=\operatorname{diag}(1,3)$. Calculate the null variance and standard deviation of $T$.

First, $\operatorname{tr}(\Sigma_1^2)=2^2+1^2=5$, and $\operatorname{tr}(\Sigma_2^2)=1^2+3^2=10$. The cross-trace is $\operatorname{tr}(\Sigma_1\Sigma_2)=2(1)+1(3)=5$. Hence

V_0=\frac{2(5)}{20(19)}+\frac{2(10)}{25(24)} +\frac{4(5)}{20(25)} =\frac{1}{38}+\frac{1}{30}+\frac{1}{25} \approx 0.09965.

The null standard deviation is $\sqrt{V_0}\approx 0.3157$. The first two contributions reflect within-group variation; the third reflects the interaction of the two covariance structures in the difference of sample means. Equal sample sizes would not make the covariance matrices equal, and the formula does not require them to be equal.

Estimating the Variance Without Pooling

The variance formula requires $\operatorname{tr}(\Sigma_1^2)$, $\operatorname{tr}(\Sigma_2^2)$, and $\operatorname{tr}(\Sigma_1\Sigma_2)$. A useful construction uses differences of disjoint observations. In group 1, form $D_r=X_{2r-1}-X_{2r}$ for $r=1,\ldots,m_1$, where $m_1=\lfloor n_1/2\rfloor$. These differences are independent and have mean zero and covariance $2\Sigma_1$. Define $E_s=Y_{2s-1}-Y_{2s}$, $s=1,\ldots,m_2$, in group 2; these are independent of the $D_r$ and have covariance $2\Sigma_2$.

Definition: For $m_1,m_2\ge2$, unbiased estimators of the three trace terms are $$ \widehat A_1=\frac{1}{4m_1(m_1-1)}\sum_{r\ne s}(D_r^\top D_s)^2,\quad \widehat A_2=\frac{1}{4m_2(m_2-1)}\sum_{r\ne s}(E_r^\top E_s)^2, $$ $$ \widehat B=\frac{1}{4m_1m_2}\sum_{r=1}^{m_1}\sum_{s=1}^{m_2}(D_r^\top E_s)^2. $$ The resulting variance estimate is $$ \widehat V=\frac{2\widehat A_1}{n_1(n_1-1)} +\frac{2\widehat A_2}{n_2(n_2-1)} +\frac{4\widehat B}{n_1n_2}. $$

To verify the first estimator, take two independent differences $D_r,D_s$. Conditional on $D_r$, the covariance of $D_s$ is $2\Sigma_1$, so

E[(D_r^\top D_s)^2\mid D_r] =D_r^\top(2\Sigma_1)D_r.

Taking expectation again gives $2\operatorname{tr}(\Sigma_1 E[D_rD_r^\top])=2\operatorname{tr}(\Sigma_1(2\Sigma_1))=4\operatorname{tr}(\Sigma_1^2)$. This is the expectation of every summand in $\widehat A_1$, proving $E[\widehat A_1]=\operatorname{tr}(\Sigma_1^2)$. The same argument applies to $\widehat A_2$. For independent $D_r,E_s$, conditioning on $D_r$ gives $E[(D_r^\top E_s)^2\mid D_r]=D_r^\top(2\Sigma_2)D_r$; its expectation is $4\operatorname{tr}(\Sigma_1\Sigma_2)$, proving the claim for $\widehat B$.

Moment condition: For unbiasedness of these estimators at fixed finite dimension, finite second moments suffice. In the conditional calculation, $E\|D_r\|^2<\infty$ ensures the final expectation is finite; no fourth-moment assumption is needed for unbiasedness. This does not by itself guarantee that the estimators concentrate around their expectations as dimension and sample size grow. Consistency and normal calibration require additional conditions.

Worked Example: Disjoint-Difference Trace Estimation

Consider two disjoint differences in one group, $D_1=(1,1)^\top$ and $D_2=(2,-1)^\top$. With $m_1=2$, the ordered-pair sum has two identical terms. Their inner product is $D_1^\top D_2=2-1=1$, so the trace estimate from these observed differences is

\widehat A_1=\frac{(D_1^\top D_2)^2+(D_2^\top D_1)^2}{4(2)(1)} =\frac{1+1}{8}=\frac14.

The value $1/4$ is a realization of an estimator, not a claim that the population trace equals $1/4$. Its unbiasedness follows from the expectation calculation above. With just two differences it can be highly variable, which is why an unbiased estimate should not be mistaken for a precise one.

Normal Calibration and Its Conditions

A common calibration uses $Z=T/\sqrt{\widehat V}$ and compares it with a standard normal distribution. This is an asymptotic approximation, not an exact finite-sample result. One precise sufficient route to a normal limit is available under Gaussian sampling. Stack the centered observations into a vector $W$, whose covariance matrix is block diagonal with $n_1$ blocks $\Sigma_1$ and $n_2$ blocks $\Sigma_2$. Since $T$ is a quadratic cross-product statistic, there is a symmetric matrix $M$ such that $T=W^\top M W$ under the null. Let $C=W$'s covariance square root times $M$ on both sides, that is $C=\Omega^{1/2}M\Omega^{1/2}$, where $\Omega=\operatorname{Cov}(W)$.

Theorem: Under the null, suppose the observations are Gaussian, $V_0>0$, and along the sequence of problems $$ \frac{\max_k|\lambda_k(C)|}{\{\operatorname{tr}(C^2)\}^{1/2}}\longrightarrow 0. $$ If also $\widehat V/V_0\xrightarrow{p}1$, then $T/\sqrt{\widehat V}\xrightarrow{d}N(0,1)$.

Proof. Diagonalizing $C$ expresses $T$ in distribution as $\sum_k\lambda_k(C)(Z_k^2-1)$, where the $Z_k$ are independent $N(0,1)$ random variables. The centering is valid because $E[T]=0$ under the null. Its variance is $2\sum_k\lambda_k(C)^2=2\operatorname{tr}(C^2)=V_0$, which also equals the variance formula already derived. The displayed eigenvalue condition makes the largest summand’s scale negligible relative to the total standard deviation; the central limit theorem for these independent centered Gaussian-square terms then gives $T/\sqrt{V_0}\xrightarrow{d}N(0,1)$. By $\widehat V/V_0\xrightarrow{p}1$, Slutsky’s theorem gives the stated studentized limit. $\square$

The two conditions play different roles. The eigenvalue condition controls the distribution of the statistic; relative consistency of $\widehat V$ controls the estimated scale. Unbiasedness alone does not establish the latter. In applications, the required consistency depends on the sampling regime and moment assumptions, so it must be justified for the particular estimator and data model rather than inferred from its expectation.

Worked Example: Turning a Statistic Into a Test

Suppose a study reports $T=0.80$, and its variance estimate is $\widehat V=0.09965$. Assuming the conditions for normal calibration are reasonable, calculate the test statistic and a two-sided approximate $p$-value.

The estimated standard deviation is $\sqrt{0.09965}\approx0.3157$, giving $Z=0.80/0.3157\approx2.53$. The two-sided normal approximation is $2\{1-\Phi(2.53)\}\approx0.0114$. At a 5% level this rejects equality of the mean vectors. This conclusion is conditional on adequate calibration: the observed value and variance estimate alone cannot verify the asymptotic assumptions.

Common Pitfall: Pooling Away the Difference

A tempting shortcut is to assume equal covariances and use a pooled estimate in the null variance. But the two within-group terms involve $\operatorname{tr}(\Sigma_1^2)$ and $\operatorname{tr}(\Sigma_2^2)$, while the cross term involves $\operatorname{tr}(\Sigma_1\Sigma_2)$. These are not interchangeable. For instance, the covariance matrices in the earlier numerical example have squared traces $5$ and $10$, and cross-trace $5$; substituting one pooled quantity changes the variance calculation even though the null hypothesis about the means is true.

The correction is to estimate each trace component separately, using independent groups and the appropriate within- or between-group differences. Also distinguish three claims: $T$ is unbiased for the squared mean difference; the difference-based trace estimates are unbiased under finite second moments; and a normal-calibrated test is justified only when the distributional and variance-estimation conditions for its limiting approximation hold.

Check Your Understanding

Use the definitions and results in this tutorial to answer the following questions.

Why does subtracting the two sample covariance trace terms make $T$ unbiased for $\|\mu_1-\mu_2\|^2$?
Which three population trace quantities determine the null variance when $\Sigma_1$ and $\Sigma_2$ may differ?
For unbiasedness of the disjoint-differences estimator of $\operatorname{tr}(\Sigma_1^2)$, is any condition beyond finite second moments needed at fixed finite dimension? Explain.
Why does unbiasedness of $\widehat V$ not, by itself, justify replacing $V_0$ by $\widehat V$ in a normal calibration?
What role does the eigenvalue condition in the Gaussian normal-limit theorem play?