Nonparametric Independence Testing for Small Sample Sizes

1 Introduction

The problem of nonparametric independence testing deals with ascertaining if two random variables are independent or not, making no parametric assumptions about their underlying distributions. Formally, given n samples (xi, yi) for i ∈ {1, ..., n} where xi ∈ Rp , yi ∈ Rq , that are drawn from a joint distribution PXY supported on X × Y ⊆ Rp+q , we want to decide between the null and alternate hypotheses

where PX, PY are the marginals of PXY w.r.t. X, Y . A test is a function from the data to {0, 1}. Tests aim to have high

∗Both authors have equally contributed to this paper.

Leila Wehbe∗

Machine Learning Dept. Carnegie Mellon University lwehbe@cs.cmu.edu

power (probability of detecting dependence, when it exists) at a prespecified allowable type-1 error rate α (probability of detecting dependence when there isn't any).

Independence testing is often a precursor to further analysis. Consider for instance conditional independence testing for inferring causality, say by the PC algorithm [Spirtes et al., 2000], whose first step is (unconditional) independence testing. It is also useful for scientific discovery like in neuroscience, to see if a stimulus X (say an image) is independent of the brain activity Y (say fMRI) in a relevant part of the brain. Since detecting nonlinear correlations is much easier than estimating a nonparametric regression function (of Y onto X), it can be done at smaller sample sizes, with further samples collected for estimation only if an effect is detected by the hypothesis test. For such situations, correlation only tests for univariate linear independence, while other statistics like mutual information that do characterize multivariate independence are hard to estimate from data, suffering from a serious curse of dimensionality. A recent popular approach for this problem (and a related two-sample testing problem) involve the use of quantities defined in reproducing kernel Hilbert spaces (RKHSs) - see [Gretton et al., 2006; Harchaoui et al., 2007; Gretton et al., 2005b; 2005a].

This paper will concern itself with increasing the statistical power at small samples of a popular kernel statistic called HSIC, by using shrunk empirical estimators of the unknown population quantity (introduced below).

1.1 Hilbert Schmidt Independence Criterion

Due to limited space, familiarity with RKHS terminology is assumed - see [Scholkopf and Smola, 2002] for an introduction. Let k : X ×X → R and l : Y×Y → R be two positivedefinite reproducing kernels that correspond to RKHSs Hk and Hl respectively with inner-products 〈·, ·〉k and 〈·, ·〉l. Let k, l arise from (implicit) feature maps φ : X → Hk and ψ : Y → Hl. In other words, φ, ψ are not functions, but mappings to the Hilbert space. i.e. φ(x) ∈ Hk, ψ(y) ∈ Hl respectively. These functions, when evaluated at points in the original spaces, must satisfy φ(x)(x′ ) = 〈φ(x), φ(x′ )〉k = k(x, x′ ) and ψ(y)(y′ ) = 〈ψ(y), ψ(y′ )〉l = l(y, y′ ).

The mean embedding of PX and PY are defined as µX := Ex∼PX φ(x) ∈ Hk and µY := Ey∼PY ψ(y) ∈ Hl whose empirical estimates are ̂ µX := 1n ∑ni=1 φ(xi) and ̂ µY := 1n ∑ni=1 ψ(yi). Finally, the cross-covariance operator of

X, Y is defined as

where ⊗ is an outer-product. For unfamiliar readers, if we used the linear kernel k(x, x′ ) = xT x′ and l(y, y′ ) = yT y′ , then the cross-covariance operator is just the cross-covariance matrix. The plug-in empirical estimator of ΣXY is

For conciseness, define ˜ φ(xi) = φ(xi) - ̂ µX, ˜ ψ(yi) = ψ(yi) - ̂ µY , ˜k(x, x′ ) = 〈˜ φ(x), ˜ φ(x′ )〉k and ˜l(y, y′ ) = 〈 ˜ ψ(y), ˜ ψ(y′ )〉l. The test statistic Hilbert-Schmidt Independence Criterion (HSIC) defined in [Gretton et al., 2005a] is the squared Hilbert-Schmidt norm of SXY , and can be calculated using centered kernel matrices ˜K, ˜L, where ˜Kij =

For unfamiliar readers, if we used the linear kernel, this just corresponds to the Frobenius norm of the cross-covariance matrix. The most important property is: when the kernels k, l are "characteristic", then the corresponding population statistic ‖ΣXY ‖2HS is zero iff X, Y are independent [Gretton et al., 2005a]. This gives rise to a natural test - calculate ‖SXY ‖2HS and reject the null if it is large.

Examples of characteristic kernels include Gaussian k(x, x′ ) = exp (- ‖x-x′ ‖22 γ 2 ) and Laplace k(x, x′ ) = exp (- ‖x-x′ ‖1 γ ), for any bandwidth γ, while the aforementioned linear kernel is not characteristic - the corresponding HSIC tests only linear relationships, and a zero crosscovariance matrix characterizes independence only for multivariate Gaussian distributions. Working with the infinite dimensional operator with characteristic kernels, allows us to identify any general nonlinear dependence (in the limit) between any pair of distributions, not just Gaussians.

1.2 Independence Testing using HSIC

A permutation-based test is described in [Gretton et al., 2005a], and proceeds in the following manner. From the given data, calculate the test statistic T := ‖SXY ‖2HS . Keeping the order of x1, ..., xn fixed, randomly permute y1, ..., yn a large number of times, and recompute the permuted HSIC each time. This destroyed any dependence between x, y simulating a draw from the product of marginals, making the empirical distribution of the permuted HSICs behave like the null distribution of the test statistic (distribution of HSIC when H0 is true). For a pre-specified type-1 error α, calculate threshold tα in the right tail of the null distribution. Reject H0 if T > tα. This test was proved to be consistent against any fixed alternative, meaning for any fixed type-1 error α, the power goes to 1 as n → ∞. Empirically, the power can be calculated using simulations by repeating the above permutation test many times for a fixed PXY (for which dependence holds), and reporting the empirical probability of rejecting the null (detecting the dependence). Note that the power depends on PXY (unknown to the user of the test).

1.3 Shrunk Estimators of SXY

Even though SXY is an unbiased estimator of ΣXY , it typically has high variance at low sample sizes. The idea of Stein shrinkage [Stein, 1956] is to trade-off bias and variance, first introduced in the context of Gaussian mean estimation. This strategy of introducing some bias and decreasing the variance to get different estimators of ΣXY was followed by [Muandet et al., 2014] who define a linear shrinkage estimator of SXY called SCOSE (Simple Covariance Shrinkage Estimator) and a nonlinear shrinkage estimator called FCOSE (Flexible Covariance Shrinkage Estimator). When we refer to shrunk estimators, we implicitly mean SCOSE and FCOSE. We will describe these briefly in Section 2.

1.4 Contributions

Our first contribution is the following :

1. We provide evidence that employing shrunk estimators of ΣXY , instead of SXY , to calculate the aforementioned test statistic, can increase the power of the associated independence test at low false positive rates, when the sample size is small (there is higher variance in estimating infinitedimensional operators).

Our second contribution is to analyze the effect of shrinkage on the test statistic, to provide some practical insight.

2. The effect of shrinkage on the test-statistic is very similar to soft-thresholding (see Section 4), shrinking very small statistics to zero, and shrinking other values nearly (but not) linearly, and nearly (but not) monotonically.

Our last contribution is an insight on the two estimators considered in this paper, SCOSE and FCOSE.

3. We prove that SCOSE is (essentially, up to lower order terms) the optimal/oracle linear shrinkage estimator with respect to quadratic risk (see Section 5). However, we observe that FCOSE typically achieves higher power than SCOSE. This indicates that it may be useful to search for the optimal estimator in a larger class than linearly shrunk estimators, and also that quadratic loss may not be the right loss function for the purposes of test power.

The rest of this paper is organized as follows. Section 2 introduces SCOSE, FCOSE and their corresponding shrunk test statistics. Section 3 presents illuminating experiments that bring out the statistically significant improvement in power over HSIC. Section 4 conducts a deeper investigation into the effect of shrinkage and proves the oracle optimality of SCOSE under quadratic risk.

2 Shrunk Estimators and Test Statistics

Let HS(Hk,Hl) represent the set of Hilbert-Schmidt operators fromHk to Hl. We first note that SXY can be written as the solution to the following optimization problem.

Using this idea [Muandet et al., 2014] suggest the following two shrunk/regularized estimators.

From SCOSE to HSICS

This is derived in [Muandet et al., 2014] by solving

and the optimal solution (called SCOSE) is

where λ (and hence the shrinkage intensity) is estimated by leave-one-out cross-validation (LOOCV), in closed form as

Observing the expression for λCV in [Muandet et al., 2014], the denominator can be negative (for example, with the Gaussian kernel for small bandwidths, resulting in a kernel matrix close to the identity). This can cause λCV to be negative, and ρS to be (unintentionally) outside the range [0, 1]. Though not discussed in [Muandet et al., 2014], we shall follow the convention that when ρS < 0, we shall use ρS = 0 and if ρS > 1, we use ρS = 1. Indeed, one can show that (1- λ 1+λ )+ SXY dominates (1- λ 1+λ ) SXY where (x)+ = max{x, 0}. In Section 4, we prove that SSXY is (essentially) the optimal/oracle linear shrinkage estimator with respect to quadratic risk.

We can now calculate the corresponding shrunk statistic HSICS = ‖SSXY ‖2HS =

While the above expression looks daunting, one thing to note is that the amount that HSIC is shrunk (i.e. the multiplicative factor) depends on the value of HSIC. As we shall see in section 4, small HSIC values get shrunk to zero, but as can be seen above, the shrinkage of HSIC is non-monotonic.

From FCOSE to HSICF

The Flexible Covariance Shrinkage Estimator is derived by relying on the Representer theorem, see [Scholkopf and Smola, 2002], to instead minimize

over all β ∈ Rn , and the optimal solution (called FCOSE) is

where ◦ denotes elementwise (Hadamard) product, 1 is the vector [1, 1, ..., 1]T , and as before the best λ is determined by LOOCV. The procedure to evaluate the optimal λ efficiently is described by [Muandet et al., 2014] - a single eigenvalue decomposition of ˜K◦˜L costing O(n3 ) can be done, following which evaluating LOOCV is only O(n2 ) per λ, see [Muandet et al., 2014], section 3.1 for more details. As before, after picking the λ by LOOCV, we can derive the corresponding shrunk test statistic as

where M = ˜K ◦ ˜L. Note here that the shrinkage is not linear, and the effect on HSIC cannot be seen immediately. Similar to SCOSE, we shall see in section 4, small HSIC values get shrunk to zero (LOOCV chooses a large λ).

3 Linear Shrinkage and Quadratic Risk

In this section, we prove that SCOSE is (essentially) optimal within a particular class of estimators. Such "oracle" arguments also exist elsewhere in the literature, like [Ledoit and Wolf, 2004], so we provide only a brief proof outline.

Proposition 1. The oracle (with respect to quadratic risk) linear shrinkage estimator and intensity is defined as

and is given by S∗ := (1 - ρ∗ )SXY where

Proof. Define α2 = ‖ΣXY ‖2HS , β2 = E‖SXY -ΣXY ‖2HS , δ2 = E‖SXY ‖2 . Since E[SXY ] = ΣXY , it is easy to verify that α2 + β2 = δ2 . Substituting and expanding the objective, we get:

Differentiating and equating to zero, gives ρ∗ = β2 δ2 .

This ρ∗ appears in terms of quantities that depend on the unknown underlying distribution (hence the term oracle estimator). We use plugin estimates b, d for β, δ. Let d2 = ‖SXY ‖2HS = 1 n2 ∑ni,j=1 ˜Kij ˜Lij = HSIC. Since β2 is the variance of SXY , let b2 be the sample variance of SXY , i.e. b2 = 1n 1n ∑nk=1 ||˜ φ(xi) ⊗ ˜ ψxi - SXY ||2 = 1n [ 1n ∑ni=1 ˜Kii ˜Lii - 1 n2 ∑ni,j=1 ˜Kij ˜Lij ]. Plugging these into S∗ and simplifying, we see that HSIC∗ := ‖S∗ ‖2HS is

Comparing Eq.(3) with Eq.(2) shows that SCOSE is essentially S∗ , up to a factor in the denominator which is of the same order as the bias of the HSIC empirical estimator 1 (see Theorem 1 in [Gretton et al., 2005a]). In other words, SCOSE just corresponds to using a slightly different estimator for δ2 than the simple plugin d2 , which varies on the same order as the bias δ2 -Ed2 . Hence SCOSE, as estimated via regularization and LOOCV, is (essentially) the optimal linear shrinkage estimator under quadratic risk.

To the best of our knowledge, this is the first such characterization of optimality of an estimator achieved through leave-one-out cross-validation. We are only able to prove this because one can explicitly calculate both the oracle linear shrinkage intensity ρ∗ as well as the optimal λCV (as mentioned in Section 2). This raises a natural open question - can we find other situations where the LOOCV estimator is optimal with respect to some risk measure? (perhaps when explicit calculations are not possible, like ridge regression).

4 Experiments

In this section, we run three kinds of experiments: a) to verify that SCOSE has better quadratic risk than FCOSE and original sample estimator, b) detailed synthetic experiments to verify that shrinkage does improve power, across interesting regimes of α = {0.01, 0.05, 0.1}, and c) real data obtained from MNIST, to show that we shrinkage detect dependence at much lower samples than the original data size.

4.1 Quadratic Risk

Figure 1 shows that SCOSE is indeed much better than both SXY and FCOSE with respect to quadratic risk. Here, we calculate E‖Z - ΣXY ‖2HS for the distribution given in dataset (A) for Z ∈ {SXY , SSXY , SFXY }. The expectation is calculated by repeating the experiment 1000 times. Each time Z is calculated according to N ∈ {20, 50, 100} samples and ΣXY is approximated by the empirical cross-covariance matrix on 5,000 samples. The four panels use four different kernels which are linear, polynomial, Laplace and Gaussian from top to bottom. The shrunk estimators are always better than the unshrunk, with a larger difference between SCOSE and FCOSE for finite-dimensional feature spaces (top two). In infinite-dimensional feature spaces (bottom two), SCOSE and FCOSE are much better than the unshrunk estimator but very similar to each other. The differences between all estimators decreases with increasing n, since the sample crosscovariance operator itself becomes very accurate.

4.2 Synthetic Data

We perform synthetic experiments in a wide variety of settings to demonstrate that the shrunk test statistics achieve higher power than HSIC in a variety of settings. We follow the schema provided in the introduction for independence testing and calculating power. We only consider difficult distributions with nonlinear dependence between X, Y , on which linear methods like correlation are shown to fail to

1 HSIC and HSIC - 2HSIC/n - C/n2 both converge to population HSIC at same rate determined by the dominant term (HSIC).

N

N

Figure 1: All panels show quadratic risk E‖X - ΣXY ‖2HS for X ∈ {SXY , SSXY , SFXY }. Dataset (A) was used in all four panels, but the kernels were varied - from top to bottom is the linear, quadratic, Gaussian and Laplace kernel.

detect dependence (some of them were used in previous papers on independence testing like [Gretton et al., 2007] and [Chwialkowski and Gretton, 2014]).

For all experiments, α ∈ {0.01, 0.05, 0.1} is chosen as the type-1 error (for choosing the threshold level of the null distribution's right tail). For every setting of parameters of each experiment, power is calculated as the percentage of rejection over 200 repetitions (independent trials), with 2000 permutations per repetition (permutation testing to find the null distribution threshold at level α). We use the Gaussian kernel where the bandwidth is chosen by the common median heuristic [Scholkopf and Smola, 2002].

Table 1 is a representative sample from what we saw on other examples - either large, small or no improvement in power was seen but almost never a worsening of power. The improvements in power may not always be huge, but they are statistically significant - it is difficult to detect such non-linear dependencies at low sample sizes, so any increase in power can be important in scientific applications.

Remark. A more appropriate way than using error bars to assess significance is by the Wilcoxon rank sum test, omitted for lack of space, though it yields more favorable results.

4.3 Real Data

We use two real datasets - the first is a good example where shrinkage helps a lot, but in the second it does not help (we show it on purpose). Like the synthetic datasets, for most real datasets it either helps or does not hurt (being very rarely worse; see remark in the discussion section).

The first is the Eckerle dataset [Eckerle, 1979] from the NIST Statistical Reference Datasets (NIST StRD) for Nonlinear Regression, data from a NIST study of circular interference transmittance (n=35, Y is transmittance, X is wavelength). A plot of the data in Figure 2 reveals a nonlinear relationship between X, Y (though the correlation is 0.035 with p-value 0.84). We subsample the data to see how often we can detect a relationship at 10%, 20%, 30% of the original data size, when the false positive level is always controlled at 0.05. The second is the Aircraft dataset [Bowman and Azzalini, 2014] (n=709, X is log(speed), Y is log(span)). Once again, correlation is low, with a p-value of over 0.8, and we subsample the data to 5%, 10%, 20% of the original data size.

5 Discussion

Why might shrinkage improve power? Let us examine the net effect of using shrunk estimators on the value of HSIC, i.e. let us compare HSICS and HSICF to HSIC by computing these over all the repetitions of the permutation testing procedure described in the introduction. In Fig. 3, both estimators are visually similar in transforming the actual test statistic. Perhaps the more interesting phenomenon is that Fig. 3 is reminiscent of the graph of a soft-thresholding operator STt(x) = max{0, x - t}. Intuitively, if the unshrunk HSIC value is small, the shrinkage methods deem it to be "noise" and it is shrunk to zero. Looking at the X-axis scaling of the top and bottom row, the size of the region that gets shrunk to zero decreases with n - as expected, shrinkage has less effect when SXY has low variance). The shrinkage being non-monotone (more so for n = 20 than n = 50 in Figure 3) is key to achieving an improvement in power.

Using the intuition from the above figure, we can finally piece together why shrinkage may yield benefits. A rejection of H0 occurs when the test statistic stands out in the right tail of its null distribution. Typically, when the alternative is true (this is when rejecting the null improves power) the unshrunk test statistics calculated from the permuted samples is smaller than the unshrunk HSIC calculated on the original sample. However, the effect of shrinking the small statistics towards zero, and setting the smallest ones to zero, is that the unpermuted test statistic under the alternative distribution stands out more in the right tail of the null.

In other words, relative to the unshrunk null distribution and the unshrunk test statistic, the tail of the null distribution is shrunk more towards zero than the unpermuted test statistic, causing the latter to have a higher quantile in the right tail of the former (relative to the quantile before shrinkage). Let us verify this experimentally. In Fig.4 we plot for each of the datasets in Table 1, the average ratio of unpermuted statistic T to the 95th percentile of the permuted statistics, for T ∈ {HSIC,HSICS ,HSICF }. Recall that for dataset (C), we didn't see much of an improvement in power, but for (A),(B),(D) it is clear from Fig. 4 that the unpermuted statistic is shrunk less than its null distribution's 95th quantile.

Remark. In our experiments, real and synthetic, shrinkage usually improves (and almost never worsens) power in falsepositive regimes that we usually care about. Will shrinkage always improve power? Possibly not. Even though shrunk the shrunk SXY dominates SXY for estimation error, it may not be the case that shrunk HSIC always dominates unshrunk HSIC for test power (i.e. the latter may not be inadmissible). However, just as no single classifier always outperforms another, it is still beneficial to add techniques like shrinkage, that seem to consistently yield benefits in practice, to the practitioner's array of tools.

6 Conclusion

We presented evidence for an important phenomenon - using biased but lower variance shrunk estimators of crosscovariance operators can often significantly improve test power of HSIC at small sample sizes. This observation (that shrinkage can improve power) has rarely been made in the statistics and machine learning testing literature. We think the reason is that most test statistics for independence testing cannot be immediately expressed as the norm of an empirical operator, making it less obvious how to apply shrinkage to improve their power at low sample sizes.

We also showed the optimality (among linear shrinkage estimators) of SCOSE, but observe that the nonlinear shrinkage of FCOSE usually yields higher power. To the best of our knowledge, there seems to be no current literature showing that the choice made by leave-one-out cross-validation (SCOSE) explicitly leads to an estimator that is "optimal" in some sense (among linear shrinkage estimators). This may be because it is often not possible to explicitly calculate the form of the LOOCV estimator, nor the explicit form of the best linear shrinkage estimator, as can both be done in this simple setting.

Since even the best possible linear shrinkage estimator (as represented by SCOSE) is usually worse than FCOSE, this result indicates that in order to improve upon FCOSE, it will be necessary to further study the class of non-linear shrinkage estimators for our infinite dimensional operators, as done for finite dimensional covariance matrices in [Ledoit and Wolf, 2011] and other papers by the same authors.

We ended with a brief investigation into the effect of shrinkage on HSIC and why shrinkage may intuitively improve power. We think that our work will be important for more powerful nonparametric detection of subtle nonlinear dependencies at low sample sizes, a common problem in scientific applications.

Acknowledgments

We would like to thank Arthur Gretton and Jessica Chemali for useful feedback on an earlier draft of the paper, Larry Wasserman for pointing out some useful references, and Krikamol Muandet for sharing his code. AR was supported in part by ONR MURI grant N000140911052.