IEICE Transactions on Information and Systems,
vol.E94-D, no.6, pp.1333-1336, 2011.

1

Least-Squares Independence Test
Masashi Sugiyama
Tokyo Institute of Technology, JapanPRESTO, Japan Science and Technology Agency, Tokyo, Japansugi@cs.titech.ac.jp
http://sugiyama-www.cs.titech.ac.jp/~sugi/
Taiji Suzuki
The University of Tokyo, Tokyo, Japans-taiji@stat.t.u-tokyo.ac.jp

Abstract
Identifying the statistical independence of random variables is one of the importanttasks in statistical data analysis. In this paper, we propose a novel non-parametricindependence test based on a least-squares density ratio estimator. Our method,
called least-squares independence test (LSIT), is distribution-free, and thus it ismore flexible than parametric approaches. Furthermore, it is equipped with a model
selection procedure based on cross-validation. This is a significant advantage overexisting non-parametric approaches which often require manual parameter tuning.
The usefulness of the proposed method is shown through numerical experiments.
Keywords
independence test, density ratio estimation, unconstrained least-squares importancefitting, squared-loss mutual information.

1IntroductionIdentifying the statistical independence of random variables is one of the fundamental
tasks in statistical data analysis. Independence tests can be used for various purposessuch as feature selection [9] and causal inference [10].
A traditional independence measure is the Pearson correlation coeﬃcient, which canbe used for detecting linear dependency. Thus, it is useful for Gaussian data, although
the Gaussian assumption is rarely fulﬁlled in practice. Recently, kernel-based independence measures have been studied in order to overcome the weakness of the Pearson
correlation coeﬃcient. The Hilbert-Schmidt independence criterion (HSIC) [4] utilizes

Least-Squares Independence Test

2cross-covariance operators on universal reproducing kernel Hilbert spaces (RKHSs) [7],which is an inﬁnite-dimensional generalization of covariance matrices. HSIC allows efﬁcient detection of non-linear dependency by making use of the reproducing propertyof RKHSs [1]. However, HSIC has a critical weakness that its performance depends onthe choice of RKHSs and there is no theoretically justiﬁed way to determine the RKHS
properly. In practice, using the Gaussian RKHS with width set to the median distance
between samples is a popular heuristic [4].
Another popular independence criterion would be mutual information [2]. In thispaper, we consider a squared-loss variant of mutual information and use its estimator leastsquares mutual information (LSMI) [9] for independence test. LSMI is also distributionfree as HSIC, but it is equipped with a natural model selection procedure based on crossvalidation, which is an advantage over HSIC. Through experiments, we show the usefulnessof the LSMI-based independence test called least-squares independence test (LSIT).

2

Least-Squares Independence Test

In this section, we propose a novel non-parametric independence test.

2.1

Formulation

Let x (∈ X ⊂ Rdx ) be an input feature and y (∈ Y ⊂ Rdy ) be an output feature, which
follow a joint probability distribution with density p(x, y). Suppose we are given a set ofi.i.d. paired samples {(xi , yi )}ni=1 . Our goal is to test whether x and y are statisticallyindependent or not, based on the samples {(xi , yi )}ni=1 .
For independence test, we employ the squared-loss mutual information (SMI) deﬁned
as follows:()2∫∫p(x, y)1p(x)p(y)− 1 dxdy,(1)SMI :=
2p(x)p(y)where p(x) and p(y) are marginal densities of x and y, respectively. SMI is the Pearson
divergence [6] from the joint density p(x, y) to the product of marginals p(x)p(y), andSMI is zero if and only if x and y are statistically independent. Hence, SMI can be used
for detecting the statistical independence of random variables.
SMI includes unknown probability densities p(x, y), p(x), and p(y), and thus it cannot
be directly computed. A naive approach is to estimate the densities p(x, y), p(x), andp(y), and plug the estimated densities in Eq.(1). However, since density estimation is
known to be a hard task and division by estimated densities can magnify the estimation
error, we consider estimating the following density ratio function directly:r(x, y) :=p(x, y).
p(x)p(y)(2)

Least-Squares Independence Test

3

Once an estimator of the density ratio rb(x, y) is obtained, SMI can be approximated using
samples as follows:nn1∑
1 ∑1dSMI :=rb(xi , yi ) − 2rb(xi , yj )2 − .
n i=1
2n i,j=1
2

2.2

(3)

Least-Squares Mutual Information

Here we explain the least-squares mutual information (LSMI) method [9], which directly
learns r(x, y) from data samples without going through density estimation of p(x, y),p(x), and p(y).
Let us approximate the density ratio (2) using the following model:
rb(x, y) =b∑αℓ ψℓ (x, y) = α⊤ ψ(x, y),ℓ=1

where ψ(x, y) : Rdx × Rdy → Rb is a non-negative basis function vector, α (∈ Rb ) is aparameter vector, and ⊤ denotes the transpose.
We determine the parameter α so that the following squared-error J0 is minimized:∫∫1(br(x, y) − r(x, y))2 p(x)p(y)dxdy
J0 (α) :=2∫∫1=rb(x, y)2 p(x)p(y)dxdy2∫∫−rb(x, y)p(x, y)dxdy + Const.
Let us denote the ﬁrst two terms by J. Since J contains the expectations over unknown
densities p(x, y), p(x), and p(y), we approximate the expectations by empirical averages.
By including an ℓ2 -regularizer, the LSMI optimization problem is formulated as follows.
[]1 ⊤cλ ⊤⊤bb := argmin
αα Hα − α h + α α ,(4)22α∈Rbwhere λ (≥ 0) is the regularization parameter that controls the strength of regularization,andn∑c := 1ψ(xi , yj )ψ(xi , yj )⊤ ,Hn2 i,j=1n∑b := 1hψ(xi , yi ).
n i=1b can be analytically obtained as
The solution αbc + λIb )−1 h,b = (Hα(5)Least-Squares Independence Test4

where Ib is the b-dimensional identity matrix. Finally, the density ratio estimator rb(x) isgiven byb ⊤ ψ(x).
rb(x) := αOnce a density ratio estimator rb(x, y) is obtained, SMI can be approximated by Eq.(3).
Thanks to the analytic-form expression, LSMI is computationally very eﬃcient. Furthermore, the above least-squares density ratio estimator was shown to possess the optimal
non-parametric convergence rate and optimal numerical stability [8, 5].

2.3

Basis Function Choice and Model Selection

Given that both {xi }ni=1 and {yi }ni=1 are normalized so that every element of x and y has
unit variance, we use the following basis functions:()()
∥x − xℓ ∥2∥y − yℓ ∥2


exp −
exp −
for regression,


2σ 2
2σ 2

ψℓ (x, y) =()∥x − xℓ ∥2exp −
δ(y = yℓ )for classiﬁcation,
2σ 2where δ(c) = 1 if the condition c is true; otherwise δ(c) = 0. We may also include a
constant basis function ϕ0 (x, y) = 1 to the above kernel basis functions.
The practical performance of LSMI depends on the choice of the kernel width σ and theregularization parameter λ. Model selection of LSMI is possible based on cross-validation
with respect to the criterion J. More speciﬁcally, the sample set Z = {(xi , yi )}ni=1 isdivided into M disjoint sets {Zm }Mbm (x) is obtained using
m=1 . Then an LSMI solution rZ\Zm (i.e., all samples without Zm ), and its J-score for the hold-out samples Zm iscomputed asCVJbm:=

12|Zm |∑
(x,y)∈Zm

rbm (x, y)2 −1|Zm |∑rbm (x, y),(x,y)∈Zm

where |Z| denotes the number of elements in the set Z. This procedure is repeated for∑M bCVm = 1, . . . , M , and the average score JbCV := M1
m=1 Jm is computed. Finally, themodel (the kernel width σ and the regularization parameter λ in the current setup) thatminimizes JbCV is chosen as the most suitable one.

2.4

Permutation Test

Our independence test procedure is based on the permutation test [3].
We ﬁrst run LSMI using the original datasets Z = {(xi , yi )}ni=1 , and obtain an SMI
destimate SMI(Z).
Next, we randomly permute {yi }ni=1 and form a shuﬄed dataset Ze ={(xi , yτ (i) )}ni=1 , where τ (·) is a randomly chosen permutation function. Then we run LSMI

Least-Squares Independence Test

5e and obtain an SMI estimate SMI(d Z).
e Noteagain using the randomly shuﬄed dataset Z,
that the random permutation eliminates the dependency between x and y (if exists), sod Z)
e would take a value close to zero.
SMI(This random permutation procedure is repeated many times, and the distribution ofde under the null-hypothesis (i.e., x and y are independent) is constructed. FiSMI(Z)dnally, the p-value is approximated by evaluating the relative ranking of SMI(Z)in thede
distribution of SMI(Z). We refer to this procedure as the least-squares independence test(LSIT).
RA MATLAB⃝
implementation of LSIT is available from‘http://sugiyama-www.cs.titech.ac.jp/∼sugi/software/LSIT/’.

3Experiments

In this section, we report experimental results.

3.1

Numerical Illustration

First, we illustrate how the proposed LSIT method behaves using the following toy
datasets with one-dimensional input x and one-dimensional output y:(A) Regression: For x ∼ U (−20, 20) where U (a, b) denotes the uniform distribution on(a, b),U (−1, 1)(Independent),y∼N (sin(20x/π), 1) (Dependent),
where N (µ, σ 2 ) denotes the normal distribution with mean µ and variance σ 2 .
(B) Classification: For y ∼ B(0.5) where B(p) denotes the binomial distribution on{−1, +1} with probability of having +1 being p,
0.5N (−1, 1) + 0.5N (1, 1) (Independent),x∼N (y, 1)(Dependent).
Examples of realized samples are plotted in Figure 1 for n = 100, where {xi }ni=1 andd Z)e{yi }ni=1 are normalized to have unit variance. Figure 2 depicts the distributions of SMI(
dand the value of SMI(Z).
The graphs show that reasonable p-values were obtained for
all the four cases. Figure 3 depicts the p-values and the frequency of accepting the nullhypothesis (i.e., x and y are independent) as functions of the sample size n. The graphs
show that LSIT works reasonably well.

Least-Squares Independence Test

6232

10yy

1
0−1−1−2−2−2−10x

1−3−22−1

12(b) Regression/Dependent

1

10.5

0.500yy(a) Regression/Independent

0x

−0.5

−0.5

−1−3−2−10x12

3−1−2−10123x(c) Classification/Independent(d) Classification/Dependent

Figure 1: Toy datasets.

3.2Performance Comparison

Here we compare the performance of LSMI and HSIC under the common permutationtest framework. HSIC is the state-of-the-art measure of statistical independence utilizing
Gaussian kernels [4]. The performance of HSIC depends on the choice of the Gaussian
width, and to the best of our knowledge, there is no theoretically justiﬁed method todetermine the kernel width. Here we use a standard heuristic of setting the Gaussian
width to the median distance between samples, which was also adopted in the originalpaper [4].
We generate data samples by[ ] [][ ]xcos θ − sin θ x′,=ysin θcos θ y ′wherex′ ∼ 0.5N (−1, 1) + 0.5N (1, 1),y ′ ∼ 0.5N (−2, 1) + 0.5N (2, 1).

Least-Squares Independence Test

7p−value=0.00

p−value=0.74
60

60

40

40

20

20

000.05

00.1

^

00.05

0.1

^

SMI

(a) Regression/Independent(b) Regression/Dependentp−value=0.22

p−value=0.00

60

60

40

40

20

20

0000.02

0.04

0.15

SMI0.06

0.08

00.05

0.1

0.15

0.20.25^^SMISMI(c) Classification/Independent(d) Classification/Dependentd Z)
e (randomly shuﬄed samples) for the toy datasets. ‘×’
Figure 2: Distributions of SMI(ddenotes the value of SMI(Z)(original samples).

Reg/dep
Cla/dep

Reg/ind
Cla/ind

1
0.8

10.8
Reg/ind

0.6

0.6

0.4

0.4

0.2

0.2

Cla/ind
Reg/dep

0Cla/dep

0102030405060
n7080

90 10010203040

5060

708090 100

nFigure 3: Experimental results for the toy datasets. Left: Mean and standard deviationof p-values over 100 runs. Right: The frequency of accepting the null hypothesis over 100
runs under the signiﬁcance level 0.05.

Least-Squares Independence Test

8Thus, (x, y) are rotation of (x′ , y ′ ) by angle θ. We conduct experiments for θ = 0 (i.e., xand y are independent) and θ = π/8 (i.e., x and y are dependent). Data samples {xi }ni=1
and {yi }ni=1 are normalized to have unit variance (see Figure 4).
The results plotted in Figure 5 show that the proposed LSIT has comparable type-Ierror (rejecting correct null-hypotheses) to HSIC, with far smaller type-II error (acceptingincorrect null-hypotheses).
Next, we vary the Gaussian width of HSIC and evaluate how the performance ischanged. In Figure 5, ‘HSIC(c)’ denotes HSIC with Gaussian width multiplied by c.
The results show that, although the type-I error does not really change with respect
to c, the type-II error is heavily aﬀected by the choice of the Gaussian kernel width.
For this dataset, the median heuristic does not work well, and c = 1/2 works the best.
However, the optimally-tuned HSIC is still outperformed by the proposed LSIT, which
is automatically tuned based on cross-validation and does not involve manual parametertuning.
Finally, we use the MNIST handwritten digit dataset for further performance evaluation. Each digit image (representing an integer in {0, 1, 2, . . . , 9}) consists of 784(= 28 × 28) pixels, each of which takes an integer value between 0 to 255 representing its intensity level in gray-scale. Here, we label the data as ‘small’ for digits ‘0’, ‘1’,‘2’, ‘3’, and ‘4’, and ‘large’ for digits ‘5’, ‘6’, ‘7’, ‘8’, and ‘9’. We randomly choose 250
samples and randomly shuﬄe the label of 250(1 − η) samples. Thus, increasing η from 0to 1 corresponds to increasing the dependence between digit patterns and labels.
The results are plotted in Figure 6, showing that the proposed LSIT has slightly larger
type-I error (i.e., lower acceptance rate when η = 0) than HSIC, but the type-II errorof LSIT is slightly smaller than HSIC (i.e., lower acceptance rate when η > 0). For thisdataset, the optimally-tuned HSIC (c = 1/2) slightly outperforms automatically-tuned
LSIT.

4Discussions and Conclusions

We proposed a novel non-parametric method of independence test based on an estimatorof a squared-loss variant of mutual information called least-squares mutual information [9].
The proposed method, which we called the least-squares independence test (LSIT), canovercome the limitation of the state-of-the-art method, the Hilbert-Schmidt independence
criterion (HSIC) [4], which is not equipped with a model selection procedure of theGaussian kernel width. Through experiments, we conﬁrmed that the proposed LSIT
compares favorably with the HSIC-based independence test.

Acknowledgment
MS was supported by SCAT, AOARD, and the JST PRESTO program. TS was supportedby MEXT Grant-in-Aid for Young Scientists (B) 22700289.

Least-Squares Independence Test

92θ=0
θ=π/8

1.5
10.5

y0−0.5−1−1.5−2−2.5−3−2−10x

123Figure 4: Rotation datasets. x and y are independent if θ = 0, and they are dependentwhen θ = π/8.

10.90.80.70.60.50.4

θ=π/8

θ=0LSITHSIC
HSIC(1/5)
HSIC(1/2)
HSIC(2)

0.3
0.2
0.1
02030

405060
n7080

90100

Figure 5: Experimental results for the rotation datasets. Frequency of accepting the null
hypothesis over 100 runs under the signiﬁcance level 0.05 is depicted.

10.9
0.8
0.7
0.6
0.5
0.4
0.3

LSIT
HSIC
HSIC(1/5)
HSIC(1/2)
HSIC(2)

0.2
0.1
00

0.2

0.4

0.6

0.8

1η

Figure 6: Experimental results for the MNIST datasets. Frequency of accepting the nullhypothesis over 50 runs under the signiﬁcance level 0.05 is depicted.

Least-Squares Independence Test

10References[1] N. Aronszajn, “Theory of reproducing kernels,” Trans. the American MathematicalSociety, vol.68, pp.337–404, 1950.
[2] T.M. Cover and J.A. Thomas, Elements of Information Theory, 2nd ed., John Wiley& Sons, Inc., 2006.
[3] B. Efron and R.J. Tibshirani, An Introduction to the Bootstrap, Chapman & Hall,
1993.
[4] A. Gretton, K. Fukumizu, C.H. Teo, L. Song, B. Schölkopf, and A. Smola, “A kernel statistical test of independence,” in Advances in Neural Information Processing
Systems 20, pp.585–592, 2008.
[5] T. Kanamori, T. Suzuki, and M. Sugiyama, “Condition number analysis of kernelbased density ratio estimation,” tech. rep., arXiv, 2009.
[6] K. Pearson, “On the criterion that a given system of deviations from the probable inthe case of a correlated system of variables is such that it can be reasonably supposed
to have arisen from random sampling,” Philosophical Magazine, vol.50, pp.157–175,
1900.
[7] I. Steinwart, “On the inﬂuence of the kernel on the consistency of support vectormachines,” J. Machine Learning Research, vol.2, pp.67–93, 2001.
[8] T. Suzuki and M. Sugiyama, “Suﬃcient dimension reduction via squared-loss mutual information estimation,” International Conference on Artiﬁcial Intelligence andStatistics, pp.804–811, 2010.
[9] T. Suzuki, M. Sugiyama, T. Kanamori, and J. Sese, “Mutual information estimation
reveals global associations between stimuli and biological processes,” BMC Bioinformatics, vol.10, no.1, p.S52, 2009.
[10] M. Yamada and M. Sugiyama, “Dependence minimizing regression with model selection for non-linear causal inference under non-Gaussian noise,” AAAI Conferenceon Artiﬁcial Intelligence pp.643–648, 2010.