Suﬃcient Dimension Reduction
via Squared-loss Mutual Information Estimation

Taiji Suzuki
The University of Tokyos-taiji@stat.t.u-tokyo.ac.jp

Abstract
The goal of suﬃcient dimension reduction
in supervised learning is to ﬁnd the lowdimensional subspace of input features thatis ‘suﬃcient’ for predicting output values. In
this paper, we propose a novel suﬃcient dimension reduction method using a squaredloss variant of mutual information as a dependency measure. We utilize an analytic
approximator of squared-loss mutual information based on density ratio estimation,
which is shown to possess suitable convergence properties. We then develop a natural gradient algorithm for suﬃcient subspace
search. Numerical experiments show that the
proposed method compares favorably withexisting dimension reduction approaches.

1

Introduction

The purpose of dimension reduction in supervised
learning is to construct a map from input features totheir low-dimensional representation which has ‘suﬃcient’ information for predicting output values. Supervised dimension reduction methods can be divided
broadly into two types—wrappers and ﬁlters (Guyon& Elisseeﬀ, 2003). The wrapper approach performs
dimension reduction speciﬁcally for a particular predictor (such as support vector classiﬁcation or Gaussian process regression), while the ﬁlter approach is
independent of the choice of successive predictors.
If one wants to enhance the prediction accuracy, thewrapper approach would be a suitable choice since predictors’ characteristics could be taken into account in
the dimension reduction phase. On the other hand,
Appearing in Proceedings of the 13th International Conference on Artiﬁcial Intelligence and Statistics (AISTATS)
2010, Chia Laguna Resort, Sardinia, Italy. Volume 9 of
JMLR: W&CP 9. Copyright 2010 by the authors.

Masashi Sugiyama
Tokyo Institute of Technology
sugi@cs.titech.ac.jp

if one wants to interpret dimension-reduced features(e.g., in bioinformatics, computational chemistry, or
brain analysis), the ﬁlter approach would be more appropriate since the extracted features are independentof the choice of successive predictors and therefore reliable in terms of interpretability. In this paper, we
focus on the ﬁlter approach.
A standard formulation of ﬁlter-type dimension reduction is suﬃcient dimension reduction (SDR), which is
aimed at ﬁnding a low-rank projection matrix suchthat, given the relevant subspace of input features, therest becomes conditionally independent of output values (Cook, 1998; Chiaromonte & Cook, 2002; Fukumizu et al., 2009). A traditional dependency measure
between random variables would be the Pearson correlation coeﬃcient (PCC). PCC can be used for detecting linear dependency, so it is useful for Gaussian data.
However, the Gaussian assumption may be rarely fulﬁlled in practice.
Recently, kernel-based dimension reduction has beenstudied in order to overcome the weakness of PCC. The
Hilbert-Schmidt independence criterion (HSIC) (Gretton et al., 2005) utilizes cross-covariance operators on
universal reproducing kernel Hilbert spaces (RKHSs)(Steinwart, 2001). Cross-covariance operators are an
inﬁnite-dimensional generalization of covariance matrices. HSIC allows one to eﬃciently detect non-linear
dependency thanks to the kernel trick (Schölkopf &Smola, 2002). Its usefulness in feature selection scenarios has been shown in Song et al. (2007). However, HSIC has several weaknesses both theoreticallyand practically. Theoretically, HSIC evaluates independence between random variables, not conditional
independence. Thus HSIC does not perform SDR ina strict sense. From the practical point of view, HSIC
evaluates the covariance between random variables,
not the correlation. This means that the change of
input feature scaling aﬀects the dimension reduction
solution, which is not preferable in practice.
Kernel dimension reduction (KDR) (Fukumizu et al.,2004) can overcome these weaknesses. KDR evaluatesSuﬃcient Dimension Reduction via Squared-loss Mutual Information Estimation

Table 1: Summary of existing and proposed dependency measures.
Methods Non-linear dependency Model selection Distribution Density estimation Feature extraction
PCC
Not detectable
Not necessary
Gaussian
Not involved
Possible
HSIC
Detectable
Not availableFree
Not involved
Possible
KDR
Detectable
Not availableFree
Not involved
Possible
HIST
Detectable
Available
Free
Involved
Not available
KDE
Detectable
Available
Free
Involved
PossibleNN
Detectable
Not availableFree
Not involved
Not available
EDGE
Detectable
Not necessary Near Gaussian
Not involved
Possible
MLMI
Detectable
AvailableFree
Not involved
Not available
LSMI
Detectable
AvailableFree
Not involved
Possibleconditional covariance using the kernel trick. Therefore, KDR directly performs SDR; furthermore, its
theoretical properties such as consistency have beenstudied thoroughly (Fukumizu et al., 2009). However,
KDR still has a weakness in practice—the performanceof KDR (and HSIC) depends on the choice of kernelparameters (e.g., the Gaussian width) and the regularization parameter. So far, there seems no model
selection method for KDR and HSIC (as discussed in
Fukumizu et al. (2009))1 .
Another popular criterion for SDR is mutual information (MI) (Cover & Thomas, 1991). MI could
be directly employed in the context of SDR sincemaximizing MI between output and projected input leads to conditional independence between output and input given the projected input. A great
deal of eﬀort has been made to estimate MI accurately, e.g., based on an adaptive histogram (HIST)(Darbellay & Vajda, 1999), kernel density estimation(KDE) (Torkkola, 2003), the nearest neighbor distance(NN) (Kraskov et al., 2004), the Edgeworth expansion (EDGE) (Hulle, 2005), and maximum likelihood
MI estimation (MLMI) (Suzuki et al., 2008). Amongthem, MLMI has been shown to possess various advantages as summarized in Table 1.
So we want to employ the MLMI method for dimension
reduction. However, this may not be possible since
the MLMI estimator is not explicit (i.e., the MLMI
estimator is implicitly deﬁned as the solution of an
optimization problem and is computed numerically)—in the dimension reduction scenarios, the projection
1

In principle, it is possible to choose the Gaussian width
and the regularization parameter by cross validation overa successive predictor. However, this is not preferable dueto the following two reasons. The ﬁrst is signiﬁcant increase of the computational cost. When cross validation
is used, the tuning parameters in KDR (or HSIC) andhyper-parameters in the target predictor (such as the kernel parameters and the regularization parameter) shouldbe optimized at the same time. This results in a deeply
nested cross validation procedure and therefore this could
be computationally very expensive. Another reason is thatfeatures extracted based on cross validation are no longer
independent of predictors. Thus a merit of the ﬁlter approach (i.e., the obtained features are ‘reliable’) is lost.

matrix needs to be optimized over an MI approximator. To cope with this problem, we adopt a squaredloss variant of MI called the squared-loss MI (SMI) as
our independence measure, and apply an SMI approximator called least-squares MI (LSMI) (Suzuki et al.,2009). LSMI inherits the good properties of MLMI,
and moreover it provides an analytic MI estimator (seeTable 1 again).
Based on LSMI, we develop a dimension reduction
algorithm called Least-squares dimension reduction(LDR). LDR optimizes the projection matrix using anatural gradient algorithm (Amari, 1998) on the Stiefel
manifold. Through numerical experiments, we show
the usefulness of the LDR method.

2

Dimension Reduction via SMI
Estimation

In this section, we ﬁrst formulate the problem ofsuﬃcient dimension reduction (SDR) (Cook, 1998;
Chiaromonte & Cook, 2002; Fukumizu et al., 2009)and show how squared-loss mutual information (SMI)
could be employed in the context of SDR. Then we introduce a method of approximating SMI without going
through density estimation and develop a dimension
reduction method. Finally, several theoretical issues
such as convergence properties of the proposed SMI
estimator are investigated.
2.1

Suﬃcient Dimension Reduction

Let DX (⊂ Rm ) be the domain of input features and
DY be the domain of output data. In the following, DY
could be multi-dimensional and either continuous (i.e.,
regression) or categorical (i.e., classiﬁcation); structured outputs can also be handled in our framework
as shown later.
The purpose of dimension reduction is to ﬁnd a good
low-dimensional representation of x which ‘describes’
output y. Here we focus on linear dimension reduction, i.e.,z = W x,where W is a projection matrix onto a d-dimensional

Suzuki, Sugiyama

subspace. That is, W is a member of the Stiefel manifold Sm
d (R):d×m
Sm| W W ⊤ = Id },d (R) := {W ∈ R

where Id is the d-dimensional identity matrix. We
assume that d is known when developing a theory andan algorithm; practically d may be chosen by crossvalidation.
SDR is the problem of ﬁnding a projection matrix W
such thaty ⊥⊥ x | z.
(1)
This means that, given the projected feature z, the(remaining) feature x is conditionally independent of
output y and therefore can be discarded without sacriﬁcing the prediction ability.
Suppose that we are given n independent and identically distributed (i.i.d.) paired samples
Dn = {(xi , yi ) | xi ∈ DX , yi ∈ DY , i = 1, . . . , n}drawn from a joint distribution with density pxy (x, y).
Our goal is to estimate the subspace (or to ﬁnd the
projection matrix on it) such that Eq.(1) is fulﬁlled.
We write zi = W xi .
A direct approach to SDR would be to determine W
so that Eq.(1) is fulﬁlled. To this end, we adopt SMI
as our criterion to be maximized with respect to W :´2
R ³ pyz (y,z)
Is (Y, Z) := 12py (y)pz (z)dydz,
py (y)pz (z) − 1(2)
where pyz (y, z) denotes the joint density of y and z,
and py (y) and pz (z) denote the marginal densities ofy and z, respectively. Is (Y, Z) allows us to evaluate
independence between y and z since Is (Y, Z) vanishes
if and only if
pyz (y, z) = py (y)pz (z).
Note that Eq.(2) corresponds to the f -divergence(Ali & Silvey, 1966; Csiszár, 1967) from pyz (y, z) topy (y)pz (z) with the squared loss, while ordinary MI
corresponds to the f -divergence with the log loss (i.e.,
the Kullback-Leibler divergence). Thus SMI could be
regarded as a natural alternative to ordinary MI.
The rationale behind SMI in the context of SDR relieson the following lemma:
Lemma 1 Let pxy|z (x, y|z), px|z (x|z), and py|z (y|z)
be conditional densities. Then we haveIs (X, Y ) − Is (Z, Y )¶2Zµpxy|z (x, y|z)1pyz (y, z)2 px (x)
=
1−
dxdy
2
px|z (x|z)py|z (y|z)
pz (z)2 py (y)≥ 0.

A proof of this lemma is given in Appendix. Lemma 1
implies that Is (X, Y ) ≥ Is (Z, Y ) and the equality
holds if and only if
pxy|z (x, y|z) = px|z (x|z)py|z (y|z),which is equivalent to Eq.(1). Thus, Eq.(1) can beachieved by maximizing Is (Z, Y ) with respect to W ;then the ‘suﬃcient’ subspace can be identiﬁed.
Now we want to ﬁnd the projection matrix W thatmaximizes Is (Z, Y ). However, SMI is inaccessible inpractice since densities pyz (y, z), py (y), and pz (z) are
unknown. Thus SMI needs to be estimated from data
samples. Our key constraint when estimating SMI isthat we want to avoid density estimation since thiswould be harder than dimension reduction itself. To
accomplish this requirement, we employ an estimatorof SMI proposed recently in Suzuki et al. (2009). The
estimator utilizes density ratio estimation instead ofdensity estimation itself based on the density ratio estimator proposed by Kanamori et al. (2009). Below,
we explain the details.
2.2

SMI Approximation via Density Ratio
Estimation

For the moment, we consider a ﬁxed projection matrix
W and let DZ = W DX . Using convex duality (Boyd& Vandenberghe, 2004), we can express SMI as
Is (Y, Z) = − inf g J(g) − 21 ,Z
1J(g) =g(y, z)2 py (y)pz (z)dydz
2Z− g(y, z)pyz (y, z)dydz,
where inf g is taken over all measurable functions. Thus
computing Is is reduced to ﬁnding the minimizer g ∗ ofJ(g)—it was shown that g ∗ is given as follows (Nguyen
et al., 2008)2 :yzg ∗ (y, z) := py (y)p
.
z (z)

p

(y,z)

(3)

Thus, estimating Is (Y, Z) amounts to estimating thedensity ratio (3).
However, directly minimizing J(g) is not possible dueto the following two reasons. The ﬁrst reason is thatﬁnding the minimizer over all measurable functions is
not tractable in practice since the search space is toovast. To overcome this problem, we restrict the search
space to some linear subspace G:
G := {α⊤ φ(y, z) | α = (α1 , . . . , αb )⊤ ∈ Rb },2(4)

A more general result—the solution is given by Eq.(3)for any f -divergence—was obtained in Nguyen et al.
(2008).

Suﬃcient Dimension Reduction via Squared-loss Mutual Information Estimation

where α is a parameter to be learned from samples,⊤
denotes the transpose of a matrix or a vector, andφ(y, z) is a basis function such that, for 0b being theb-dimensional vector with all zeros,φ(y, z) := (φ1 (y, z), . . . , φb (y, z))⊤ ≥ 0b for all y, z.
Note that φ(y, z) could be dependent on the samples {(yi , zi )}ni=1 , i.e., kernel models are also allowed.
Later, we explain how the basis functions φ(y, z) aredesigned in practice.
The second reason why directly minimizing J(g) is notpossible is that the true probability densities pyz (y, z),py (y), and pz (z) contained in the density ratio (3) areunavailable. To cope with this problem, we approximate them by their empirical distributions—then we
have
¸
·
1 ⊤c
b ⊤ α + λ α⊤ Rα , (5)b := argminαα Hα − h22α∈Rbwhere we included λα⊤ Rα (λ > 0) for regularization
purposes; R is some positive deﬁnite matrix,

andc := 12 Pn ′ φ(yi , zi′ )φ(yi , zi′ )⊤ ,Hi,i =1nPn
1bh := n i=1 φ(yi , zi ).

Diﬀerentiating the objective function (5) with respect
to α and equating it to zero, we obtainbc + λR)−1 h.
b = (Hαb ⊤αcαb − 12 αb ⊤Hb − 12 .
Ibs (Y, Z) := h(6)

Least-squares Dimension Reduction

Next, we show how the LSMI approximator could be
employed in dimension reduction scenarios. Since W
is a projection matrix, dimension reduction involves anoptimization problem over the Stiefel manifold Sm
d (R).
Here we employ a gradient ascent algorithm to ﬁnd
the maximizer of the LSMI approximator with respect
to W . After a few lines of calculations, we can showthat the gradient is given by∂ Ibs∂Wℓ,ℓ′

∂h
∂H
b
b −α
b − β)
b − β)b ⊤ ∂W
= ∂W
(2α
(3α
′
′ 2
b⊤

c

ℓ,ℓ

Then the geodesic from W to the direction of the natural gradient ∇Ibs (W ) over Sm
d (R) can be expressedusing t (∈ R) as
´
³
⊤
∂ Ibs
∂ Ibs
Wt := W exp t(W ⊤ ∂W
− ∂W
W) ,
where ‘exp’ for a matrix denotes the matrix exponential. Thus line search along the geodesic in the natural
gradient direction is equivalent to ﬁnding the maximizer from {Wt | t ≥ 0}. More details of the geometric structure of the Stiefel manifold can be found in
Nishimori and Akaho (2005).
For choosing the step size of each gradient update, we
may use Armijo’s rule (Patriksson, 1999). We call the
proposed dimension reduction algorithm Least-squares
Dimension Reduction (LDR). The entire algorithm issummarized in Figure 1.
2.4

Thus, the solution can be computed analytically bysolving a system of linear equations.
Then wecan analytically approximate SMI as follows, whichis called least-squares mutual information (LSMI)(Suzuki et al., 2009):

2.3

b

∂ Is
In the Euclidean space, the ordinary gradient ∂Wgives the steepest direction. However, on a manifold,the natural gradient (Amari, 1998) gives the steepest
direction. The natural gradient ∇Ibs (W ) at W is the∂ Ibs
projection of the ordinary gradient ∂W
to the tangentmspace of Sd (R) at W . If the tangent space is equippedwith the canonical metric ⟨G1 , G2 ⟩ = 21 tr(G⊤1 G2 ), thenatural gradient is given by
´
³
⊤
∂ Ibs
∂ Ibs
W .
∇Ibs (W ) = 12 ∂W− W ∂W

Convergence Analysis

Here, we analyze convergence properties of LSMI forparametric and non-parametric setups.
Let us begin with the case where the function class G
is a parametric model:
G = {gθ (y, z) | θ ∈ Θ ⊂ Rb }.
Suppose that the true density ratio g ∗ is containedin the model G, i.e., there exists θ ∗ (∈ Θ) such thatg ∗ = gθ∗ . Then we have the following theorem (itsproof is omitted due to lack of space).
Theorem 1 We have√Ibs (Y, Z) − Is (Y, Z) = Op (1/ n),where Op denotes the asymptotic order in probability.
Furthermore, we have1
EDn [Ibs (Y, Z) − Is (Y, Z)] = 2ntr(A−1 B) + o(1/n),ℓ,ℓ∂Rb ⊤ ∂Wb+ λα(βb − α),′ℓ,ℓc + λR)−1 Hcα.
bwhere βb := (Hwhere EDn denotes the expectation over data samples
Dn . A and B are b × b matrices deﬁned as
Aℓ,ℓ′ :=Epx pz [∂ℓ gθ∗ (y, z)∂ℓ′ gθ∗ (y, z)],Suzuki, Sugiyama

1.
2.
3.
4.

Initialize projection matrix W .
Optimize Gaussian width σ and regularization parameter λ by CV (explained later).
Update W by W ← Wε , where step-size ε may be chosen using Armijo’s rule.
Repeat 2. and 3. until W converges.
Figure 1: The LDR algorithm.

Bℓ,ℓ′ :=Epyz [(∂ℓ gθ∗ (y, z) − Epz′ |y [∂ℓ gθ∗ (y, z ′ )]GM := {g ∈ G | R(g) ≤ M },− Epy′ |z [∂ℓ gθ∗ (y ′ , z)] + Epy′ z′ [∂ℓ gθ∗ (y ′ , z ′ )])× (∂ℓ′ gθ∗ (y, z) − Epz′ |y [∂ℓ′ gθ∗ (y, z ′ )]− Epy′ |z [∂ℓ′ gθ∗ (y ′ , z)] + Epy′ z′ [∂ℓ′ gθ∗ (y ′ , z ′ )])],
where y ′ and z ′ are copies of y and z. The partial
derivative ∂ℓ is taken with respect to the ℓ-th element
θℓ of the parameter θ.
This theorem means that LSMI retains optimality in1terms of the order of convergence in n since Op (n− 2 ) is
the optimal convergence rate in the parametric setup.
Next, we consider non-parametric cases. Let G be a
general set of functions on DY × DZ . For a functiong (∈ G), let us consider a non-negative regularization
functional R(g) such thatsupy,z [g(y, z)] ≤ R(g).

(7)

If G is an RKHS with kernel k(·, ·) and there exists C such that supy,z k((y, z), (y, z)) ≤ C, R(g) :=√C∥g∥G satisﬁes Eq.(7):g(y, z) = ⟨k((y, z), ·), g(·)⟩
p√≤ k((y, z), (y, z))∥g∥G ≤ C∥g∥G ,where we used the reproducing property of the kernel (Aronszajn, 1950) and Schwartz’s inequality. Notethat the Gaussian RKHS satisﬁes this with C = 1.
Let us consider a non-parametric version of the problem (5):
w
b :=
n
n
h 1 X
i
1X
λn
2
2
argmin
g(y
,
z
)
−
R(g)
.
g(y
,
z
)+
i
j
i
i
2n2i,j=1
n i=1
2
g∈G

Note that if G is an RKHS, the above optimization
problem is reduced to a form of the ﬁnite dimensional
optimization problem (5). We assume that the true
density ratio function g ∗ (y, z) is contained in model G
and is bounded from above:g ∗ (y, z) < M0 for all (y, z) ∈ DY × DZ .
We also assume that there exists γ (0 < γ < 2) suchthat
H[] (GM , ϵ, L2 (py pz )) = O((M/ϵ)γ ),where H[] is the bracketing entropy of GM with respect
to the L2 (py pz )-norm (van der Vaart & Wellner, 1996).
This quantity represents a complexity of function class
G—the larger γ is, the more complex the function class
G is. The Gaussian RKHS satisﬁes this condition for
arbitrarily small γ (Steinwart & Scovel, 2007). Then
we have the following theorem (its proof is omitted due
to lack of space; we used Theorem 5.11 in van de Geer
(2000)).
Theorem 2 Under the above setting, if λn → 0 and2/(2+γ)λ−1), then we haven = o(nIbs (Y, Z) − Is (Y, Z) = Op (max(λn , n−1/2 )).

(8)

2/(2+γ)The conditions λn → 0 and λ−1) roughlyn = o(nmeans that the regularization parameter λn should
be suﬃciently small but not too small. This theorem shows that the convergence rate of non-parametric
LSMI is slightly slower than the parametric counterpart (Op (n−1/2 )), but the non-parametric method
would require a milder model assumption for eliminating the modeling error. According to Nguyen et al.
(2008) where a log-loss version of the above theorem
has been proven in the context of KL-divergence estimation, the above convergence rate achieves the optimal minimax rate under some setup. Thus the convergence property of non-parametric LSMI would also
be optimal in the same sense.

2.5

Model Selection by Cross Validation

As shown above, LSMI has preferable convergenceproperties. Nevertheless, its practical performance depends on the choice of basis functions and the regularization parameter. In order to determine basis
functions φ(y, z) and the regularization parameter λ,cross validation (CV) is available for the LSMI estimator: First, the samples {(yi , zi )}ni=1 are divided into K
disjoint subsets {Sk }K
k=1 of (approximately) the sameb Sk is obtained using {Sj }j̸=ksize. Then an estimator α(i.e., without Sk ) and the approximation error for the
hold-out samples Sk is computed; this procedure isrepeated for k = 1, . . . , K and its mean Jb(K-CV) is
outputted:´
PK ³ 1 ⊤ c
1b ⊤αb Sk − hb Sk . (9)b Sk H α
Jb(K-CV) := K
k=1 2 αSuﬃcient Dimension Reduction via Squared-loss Mutual Information Estimation

For model selection, we compute Jb(K-CV) for all modelcandidates (the basis function φ(y, z) and the regularization parameter λ) and choose the best model thatminimizes Jb(K-CV) . We can show that Jb(K-CV) is analmost unbiased estimator of the objective function J,
where ‘almost’-ness comes from the fact that the sample size is reduced in the CV procedure due to data
splitting (Schölkopf & Smola, 2002).
For the parametric setup, we can derive an asymptotic unbiased estimator of J (a.k.a. an information
criterion (Akaike, 1974)) based on Theorem 1, which
could be employed for model selection. However, weomit the detail due to lack of space.
2.6

Design of Basis Functions

The above CV procedure would be useful when goodcandidates of basis functions are prepared. Here we
propose to use the product kernel of the following formas basis functions:φℓ (y, z) = ϕyℓ (y)ϕzℓ (z)since the number of kernel evaluation when computing
b ℓ,ℓ′ is reduced from n2 to 2n:H¡¢¡¢
b ℓ,ℓ′ = 12 Pn ϕy (yi )ϕy′ (yi ) Pn ϕz (zj )ϕz′ (zj ) .
H
ℓ
i=1 ℓ
j=1 ℓ
ℓ
n
In the regression scenarios where y is continuous, we use the Gaussian
kernel as ¢the ‘base’ ker¡22z
nels:¡ ϕyℓ (y) := exp −∥y−uℓ ∥ /(2σ ) and ϕℓ (z) :=¢exp −∥z − vℓ ∥2 /(2σ 2 ) , where {(uℓ , vℓ )}bℓ=1 are Gaussian centers randomly chosen from {(yi , zi )}ni=1 —more
precisely, we set uℓ := yc(ℓ) and vℓ := zc(ℓ) , where{c(ℓ)}bℓ=1 are randomly chosen from {1, . . . , n} without replacement.
The rationale behind this basis function choice is asfollows: The density ratio tends to take large values
if py (y)pz (z) is small and pyz (y, z) is large (and viseversa). When a non-negative function is approximatedby a Gaussian kernel model, many kernels may beneeded in the region where the output of the targetfunction is large; on the other hand, only a small number of kernels would be enough in the region where theoutput of the target function is close to zero. Following this heuristic, we decided to allocate many
kernels in the regions where pyz (y, z) is large; thiscan be achieved by setting the Gaussian centers at3{(yi , zi )}ni=1 .
In the classiﬁcation scenarios where y is categorical,
we use the delta kernel for y: ϕyℓ (y) := δ(y = uℓ ),3
Alternatively, we may locate n2 Gaussian kernels at{(yi , zj )}ni,j=1 . However, in our preliminary experiments,
this did not further improve the performance, but signiﬁcantly increased the computational cost.

where δ(y = uℓ ) is 1 if y = uℓ and 0 otherwise. Moregenerally, when y is structured (e.g., strings, trees, andgraphs), we may employ kernels for structured data(Gärtner, 2003) as ϕyℓ (y).

3

Numerical Experiments

In this section, we experimentally investigate the performance of the proposed and existing dimension reduction methods using artiﬁcial and real datasets. In
the proposed method, we use the Gaussian kernel
as basis functions and employ the regularized kernel
Gram matrix as the regularization matrix R: R =f + ϵIb , where K
f is the kernel Gram matrix for theKe ℓ,ℓ′ := ϕy (uℓ′ )ϕz (vℓ′ ). ϵIb is added tochosen centers: Kℓℓf for avoiding non-degeneracy; we set ϵ = 0.01. WeK
ﬁx the number of basis functions at b = min(100, n),and choose the Gaussian width σ and the regularization parameter λ based on 5-fold CV with grid search.
we restart the natural gradient search 10 times withrandom initial points and choose the one having theminimum CV score (9).
3.1

Dimension Reduction for Artiﬁcial
Datasets

We use 6 artiﬁcial datasets (3 datasets designed by us
and 3 datasets borrowed from Fukumizu et al. (2009);
see Figure 2), and compare LDR with kernel dimension reduction (KDR) (Fukumizu et al., 2009), theHilbert-Schmidt independence criterion (HSIC) (Gretton et al., 2005), sliced inverse regression (SIR) (Li,1991), and sliced average variance estimation (SAVE)(Cook, 2000). In KDR and HSIC, the Gaussian widthis set to the median sample distance, following the suggestions in the original papers (Gretton et al., 2005;
Fukumizu et al., 2009). We evaluate the performanceof each method byc ⊤W
c − W ∗⊤ W ∗ ∥Frobenius ,
√1 ∥W
2d
c is
where ∥ · ∥Frobenius denotes the Frobenius norm, W∗
an estimated projection matrix, and W is the optimal
projection matrix. Note that the above error measuretakes its value in [0, 1].
The performance of each method is summarized in Table 2, which depicts the mean and standard deviationof the above Frobenius-norm error over 50 trials when
the number of samples is n = 100. LDR overall shows
good performance; in particular, it performs the bestfor datasets (b), (c), and (e). KDR also tends to workreasonably well, but it sometimes performs poorly; thisseems to be caused by the inappropriate choice of theGaussian kernel width, implying that the heuristic ofusing the median sample distance as the kernel width isSuzuki, Sugiyama

8

0

6
y

10

1
y

2

−1

4
2

−2

0

−3

−2
−2

−1

x1

0

1

2

−2

(a) Linear

0
x1

2

(b) Non-linear dependence 1

1

20

0

10

y

y

Table 2: Mean and standard deviation of Frobeniusnorm error for toy datasets. The best method in termsof the mean error and comparable ones based on the
one-sided t-test at the signiﬁcance level 1% are indicated by boldface.
LDR
KDR
HSIC
SIR
SAVE
.13(.04) .13(.05) .17(.07) .11(.04) .19(.10).15(.06) .25(.21) .44(.36) .83(.20) .24(.08).10(.05) .44(.32) .68(.32) .85(.20) .31(.11).20(.14) .16(.06) .18(.08) .29(.16) .40(.17).09(.06) .13(.06) .16(.07) .20(.10) .21(.14).35(.12) .40(.12) .49(.17) .64(.22) .73(.20)

0

−1

2

−2

m d
(a) 5 1(b) 5 1(c) 5 1(d) 4 2(e) 4 1(f) 10 1

−0.5

0
x1

0

0.5

(c) Non-linear
dence 2x2

depen-

2

0

−2

−2

x1

(d)
Fukumizual. (2009)

Table 3: Mean and standard deviation of misclassiﬁcation rates for benchmark datasets.

et

image
4
1

y

y

2
0.5
0

waveform
0
−2

−0.5

pima

−4
0

0.2

0.4

x1

0.6

0.8

(e)
Fukumizual. (2009)

1

et

−2

−1

0
x1

(f)
Fukumizual. (2009)

1

2

et

letter(a, b & c)

d
5
9
14
6
11
16
2
4
6
4
8
12

LDR
KDR
HSIC
.078(.017) .115(.034) .154(.040).091(.015) .100(.031) .107(.025).090(.019) .086(.017) .088(.019).133(.015) .132(.013) .157(.017).123(.012) .138(.013) .162(.013).117(.009) .137(.010) .159(.016).249(.022) .246(.018) .252(.021).251(.019) .255(.020) .261(.027).249(.021) .246(.021) .253(.021).029(.007) .025(.009) .034(.010).026(.008) .018(.009) .019(.007).015(.007) .015(.007) .015(.007)

Figure 2: Artiﬁcial datasets.
not always appropriate. On the other hand, LDR with
CV performs stably well for various types of datasets.

density estimation is not involved, it is distributionfree, and model selection by cross validation is available. The eﬀectiveness of the proposed method overexisting methods was shown through experiments.

3.2

Acknowledgements

Classiﬁcation for Benchmark Datasets

Finally, we evaluate the classiﬁcation performanceafter dimension reduction for several benchmark
datasets. We use ‘image’, ‘waveform’, ‘pima-indiansdiabetes’, and ‘letter recognition’ in the UCI repository. We randomly choose 200 samples from thedataset and apply LDR, KDR, and HSIC to obtainprojections onto low-dimension subspaces with d =⌈m/4⌉, ⌈m/2⌉, and ⌈3m/4⌉. Then we train the support vector machine on the projected 200 training samples.
The generalization error is computed for the samples
not used for training. Table 3 summarizes the meanand standard deviation of the classiﬁcation error over
20 iterations. This shows that the proposed method
overall compares favorably with the other methods.

4

Conclusions

In this paper, we proposed a new dimension reduction
method utilizing a squared-loss variant of mutual information (SMI). The proposed method inherits several preferable properties of the SMI estimator, e.g.,

T.S. was supported in part by Global COE Program“The research and training center for new development
in mathematics”, MEXT, Japan. M.S. was supported
by AOARD, SCAT, and the JST PRESTO program.
Proof of Lemma 1
Let x = (z, z⊥ ). By dx = dzdz⊥ , we haveIs (X, Y ) − Is (Z, Y )R ³ pxy (x,y) ´2= 12px (x)py (y)dxdy
px (x)py (y)³´2Rpyz (y,z)− 12
pz (z)py (y)dzdy
pz (z)py (y)³´2Rpxy (x,y)pyz (y,z)px (x)py (x)dxdy,−= 21
px (x)py (y)pz (z)py (y)R pxy (x,y)pyz (y,z)
becausepx (x)py (y)py (y)pz (z) px (x)py (y)dxdy
R ³ pyz (y,z) ´2
Noticingpy (y)pz (z) py (y)pz (z)dydz.
pxy
px py

pxy|z pyz
px|z py|z py pz

=
that

=
where we used the relationpx (x) = px|z (x|z)pz (z), Then we have
´2 2
R³
p
pyz px
1 − px|zxy|z
Is (X, Y ) − Is (Z, Y ) = 12py|zpz 2 py dxdy.

Suﬃcient Dimension Reduction via Squared-loss Mutual Information Estimation

References
Akaike, H. (1974). A new look at the statistical model
identiﬁcation. IEEE Transactions on Automatic Control, AC-19, 716–723.
Ali, S. M., & Silvey, S. D. (1966). A general class of coeﬃcients of divergence of one distribution from another.
Journal of the Royal Statistical Society, Series B, 28,
131–142.
Amari, S. (1998). Natural gradient works eﬃciently inlearning. Neural Computation, 10, 251–276.
Aronszajn, N. (1950). Theory of reproducing kernels.
Transactions of the American Mathematical Society, 68,
337–404.
Boyd, S., & Vandenberghe, L. (2004). Convex optimization.
Cambridge: Cambridge University Press.
Chiaromonte, F., & Cook, R. D. (2002). Suﬃcient dimension reduction and graphics in regression. Annals of theInstitute of Statistical Mathematics, 54, 768–795.
Cook, R. D. (1998). Regression graphics: Ideas for studying
regressions through graphics. New York: Wiley.
Cook, R. D. (2000). SAVE: A method for dimension reduction and graphics in regression. Communications in
Statistics - Theory and Methods, 29, 2109–2121.
Cover, T. M., & Thomas, J. A. (1991). Elements of information theory. N. Y.: John Wiley & Sons, Inc.
Csiszár, I. (1967). Information-type measures of diﬀerence of probability distributions and indirect observation. Studia Scientiarum Mathematicarum Hungarica,
2, 229–318.
Darbellay, G. A., & Vajda, I. (1999). Estimation of theinformation by an adaptive partitioning of the observation space. IEEE Transactions on Information Theory,45, 1315–1321.
Fukumizu, K., Bach, F. R., & Jordan, M. I. (2004). Dimensionality reduction for supervised learning with reproducing kernel Hilbert spaces. Journal of Machine
Learning Research, 5, 73–99.
Fukumizu, K., Bach, F. R., & Jordan, M. I. (2009). Kernel dimension reduction in regression. The Annals of
Statistics, 37, 1871–1905.
Gärtner, T. (2003). A survey of kernels for structured data.
ACM SIGKDD Explorations Newsletter, 5, 49–58.
Gretton, A., Bousquet, O., Smola, A., & Schölkopf, B.
(2005). Measuring statistical dependence with HilbertSchmidt norms. Algorithmic Learning Theory (pp. 63–
77). Berlin: Springer-Verlag.
Guyon, I., & Elisseeﬀ, A. (2003). An introduction to variable feature selection. Journal of Machine Learning Research, 3, 1157–1182.
Hulle, M. M. V. (2005). Edgeworth approximation of multivariate diﬀerential entropy. Neural Computation, 17,
1903–1910.

Kanamori, T., Hido, S., & Sugiyama, M. (2009). A leastsquares approach to direct importance estimation. Journal of Machine Learning Research, 10, 1391–1445.
Kraskov, A., Stögbauer, H., & Grassberger, P. (2004). Estimating mutual information. Physical Review E, 69,066138.
Li, K.-C. (1991). Sliced inverse regression for dimension
reduction. Journal of American Statistical Association,
86, 316–342.
Nguyen, X., Wainwright, M. J., & Jordan, M. I. (2008). Estimating divergence functions and the likelihood ratio by
penalized convex risk minimization. Advances in Neural Information Processing Systems 20 (pp. 1089–1096).
Cambridge, MA: MIT Press.
Nishimori, Y., & Akaho, S. (2005). Learning algorithms
utilizing quasi-geodesic ﬂows on the Stiefel manifold.
Neurocomputing, 67, 106–135.
Patriksson, M. (1999). Nonlinear programming and variational inequality problems. Dredrecht: Kluwer Academic.
Schölkopf, B., & Smola, A. J. (2002). Learning with kernels. Cambridge, MA: MIT Press.
Song, L., Smola, A., Gretton, A., Borgwardt, K. M., &Bedo, J. (2007). Supervised feature selection via dependence estimation. Proceedings of the 24th International
Conference on Machine learning (pp. 823–830). New
York, NY, USA: ACM.
Steinwart, I. (2001). On the inﬂuence of the kernel on
the consistency of support vector machines. Journal of
Machine Learning Research, 2, 67–93.
Steinwart, I., & Scovel, C. (2007). Fast rates for support
vector machines using Gaussian kernels. The Annals of
Statistics, 35, 575–607.
Suzuki, T., Sugiyama, M., Kanamori, T., & Sese, J. (2009).
Mutual information estimation reveals global associations between stimuli and biological processes. BMC
Bioinformatics, 10, S52.
Suzuki, T., Sugiyama, M., Sese, J., & Kanamori, T. (2008).
Approximating mutual information by maximum likelihood density ratio estimation. JMLR Workshop andConference Proceedings, 4, 5–20.
Torkkola, K. (2003). Feature extraction by non-parametricmutual information maximization. Journal of Machine
Learning Research, 3, 1415–1438.
van de Geer, S. (2000).
Empirical processes in Mestimation. Cambridge University Press.
van der Vaart, A. W., & Wellner, J. A. (1996). Weakconvergence and empirical processes: With applications
to statistics. Springer, New York.