Neural Computation. vol.25, no.3, pp.725–758, 2013.

Suﬃcient Dimension Reduction
via Squared-loss Mutual Information Estimation
Taiji Suzuki
Department of Mathematical Informatics, The University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japant-suzuki@mist.i.u-tokyo.ac.jp
Masashi Sugiyama
Department of Computer Science, Tokyo Institute of Technology
2-12-1 O-okayama, Meguro-ku, Tokyo 152-8552, Japansugi@cs.titech.ac.jp
Abstract
The goal of suﬃcient dimension reduction in supervised learning is to ﬁnd the lowdimensional subspace of input features that contains the whole information about
the output values that the input features possess. In this paper, we propose anovel suﬃcient dimension reduction method using a squared-loss variant of mutual information as a dependency measure. We apply a density-ratio estimator for
approximating squared-loss mutual information that is formulated as a minimum
contrast estimator on parametric or non-parametric models. Since cross-validation
is available for choosing an appropriate model, our method does not require any
pre-speciﬁed structure on the underlying distributions. We elucidate the asymptotic bias of our estimator on parametric models and the asymptotic convergence
rate on non-parametric models. The convergence analysis utilizes the uniform tailbound of a U -process, and the convergence rate is characterized by the bracketing
entropy of the model. We then develop a natural gradient algorithm on the Grassmann manifold for suﬃcient subspace search. The analytic formula of our estimatorallows us to compute the gradient eﬃciently. Numerical experiments show that the
proposed method compares favorably with existing dimension reduction approaches
on artiﬁcial and benchmark datasets.
Keywords: dimension reduction, density ratio, convergence rate, mutual information.

1

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

1

2

Introduction

The purpose of dimension reduction in supervised learning is to ﬁnd the low-dimensional
subspace of input features which has ‘suﬃcient’ information for predicting output values.
Suﬃcient dimension reduction (SDR) initiated by Li (1991) is aimed at ﬁnding a lowrank projection matrix such that, given the relevant subspace of input features, the rest
becomes conditionally independent of output values (Cook, 1998b; Chiaromonte & Cook,
2002). Such a low-dimensional subspace contains all the information of the output thatthe covariate contains. Finding such a subspace not only allows us to use the dimensionreduced features for estimating input-output relations, but also gives insights about whichfeatures are important.
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 isuseful for Gaussian data. However, the Gaussian assumption may be rarely fulﬁlled inpractice.
Recently, kernel-based dimension reduction has been studied in order to overcome theweakness 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 by making use of the reproducing property of RKHSs (Aronszajn, 1950). Its
usefulness in feature selection scenarios has been shown in Song et al. (2007). However,
HSIC has several weaknesses both theoretically and practically. Theoretically, HSIC evaluates independence between random variables, not conditional independence. Thus HSIC
does not perform SDR in a 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 conditional covariance using the kernel trick, and thus KDR
directly performs SDR. Through experiments, KDR was demonstrated to outperform
other dimension reduction schemes such as canonical correlation analysis (Hotelling, 1936;
Breiman & Friedman, 1985), partial least-squares (Wold, 1966; Goutis & Fearn, 1996; Durand & Sabatier, 1997; Reiss & Ogden, 2007), sliced inverse regression (Li, 1991; Bura &Cook, 2001; Cook & Ni, 2005; Zhu et al., 2006), and the principal Hessian direction (Li,
1992; Cook, 1998a; Li et al., 2000). Theoretical properties of KDR such as consistency
have been studied thoroughly (Fukumizu et al., 2009). However, KDR still has a weakness
in practice—the performance of KDR (and also 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 .
1

In principle, it is possible to choose the Gaussian width and the regularization parameter by crossvalidation (CV) over a successive predictor. However, this is not preferable due to the following two

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

3

Table 1: Summary of existing and proposed dependency measures.
Methods
PCC
HSIC
KDR
HIST
KDE
NN
EDGE
MLMI
LSMI

Non-linear
dependence
Not detectable
Detectable
Detectable
Detectable
Detectable
Detectable
Detectable
Detectable
DetectableModelselection
Not necessary
Not available
Not available
Available
Available
Not available
Not necessary
Available
Available

Distribution
Gaussian
Free
Free
Free
Free
Free
Near Gaussian
Free
Free

Densityestimation
Not involved
Not involved
Not involved
Involved
Involved
Not involved
Not involved
Not involved
Not involved

Featureextraction
Possible
Possible
Possible
Not available
Possible
Not available
Possible
Not available
Possible

Another possible criterion for SDR is mutual information (MI) (Cover & Thomas,
2006). MI could be directly employed for SDR since maximizing MI between output andprojected input leads to conditional independence between output and input given the
projected input. So far, 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). Among them, MLMI has been shown to possess
various practical advantages.
As summarized in Table 1, MLMI aﬀords model selection by cross-variation, while
there is no systematic method to choose tuning parameters for HSIC, KDR, and NN.
MLMI does not require speciﬁc structures on the underlying distributions, while EDGE
requires that the distributions are near Gaussian. MLMI does not involve density estimation of the underlying distributions so that it shows a good performance in practice.
Based on the above comparison, we want to employ the MLMI method for dimension
reduction. However, this is not straightforward 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 (or feature extraction) scenarios,
the projection matrix needs to be optimized over an MI approximator. To cope with thisproblem, we adopt a squared-loss variant of MI called the squared-loss MI (SMI) as
our independence measure, and use an estimator of SMI called least-squares MI (LSMI)(Suzuki et al., 2009) for dimension reduction. LSMI inherits good properties from MLMI,
and moreover it provides an analytic SMI estimator that gives an analytic formula for its
derivative (see Table 1 again).
The goal of this paper is to develop a dimension reduction algorithm based on LSMI.
reasons. The ﬁrst is signiﬁcant increase of the computational cost. When CV is used, the tuning
parameters in KDR (or HSIC) and hyper-parameters in the target predictor (such as basis parameters
and the regularization parameter) should be optimized at the same time. This results in a deeply nested
CV procedure, and therefore this could be computationally very expensive. Another reason is thatfeatures extracted based on CV are no longer independent of predictors, which is not preferable from theviewpoint of interpretability.

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

4

Our ﬁrst contribution in this paper is to theoretically elucidate the rate of convergenceof the LSMI estimator in parametric and non-parametric settings. Then we developa practical dimension reduction algorithm based on LSMI, which we call Least-Squares
Dimension Reduction (LSDR). LSDR optimizes the projection matrix using a natural
gradient algorithm (Amari, 1998) on the Grassmann manifold. Finally, through numericalexperiments, we show the usefulness of the LSDR method.

2

Dimension Reduction via SMI Estimation

In this section, we ﬁrst formulate the problem of suﬃcient dimension reduction (SDR)(Cook, 1998b; Chiaromonte & Cook, 2002), and show how squared-loss mutual information (SMI) can be employed in the context of SDR. Then we introduce a method of
approximating SMI without going through density estimation, and we elucidate convergence properties of the SMI estimator. Finally, we develop a dimension reduction method
based on the SMI estimator.

2.1

Suﬃcient Dimension Reduction

Let DX (⊂ Rd ) be the domain of input feature x, and DY be the domain of output data2 y.
We suppose that DY is equipped with a σ-algebra BY and there is a base measure denotedby dy. As for DX , we consider the standard σ-algebra BX of the Lebesgue measurable sets
and denote by dx as the Lebesgue measure. We assume there is a joint density pxy (x, y)deﬁned on the product space (DX × DY , BX × BY ) with respect to dx × dy.
To search a subspace of input space containing suﬃcient information about the output,
we utilize the Grassmann manifold Grdm (R) that is the set of all m-dimensional subspaces
in Rd . The Grassmann manifold Grdm (R) is obtained by identifying those matrices ind × m orthonormal matrices whose columns span the same subspace:
Grdm (R) := {W ∈ Rm×d | W W ⊤ = Im }/ ∼,(1)

where ⊤ denotes the transpose, Im is the m-dimensional identity matrix, and ∼ is theequivalence relation such that W ∼ W ′ if the raws of both W and W ′ span the same
space.
Let W ∗ be any projection matrix corresponding to a member of the Grassmann manifold Grdm (R). Let z ∗ (∈ Rm ) be the orthogonal projection of input x given by W ∗ :z ∗ = W ∗ x.
Suppose that z ∗ satisﬁesy⊥⊥ x | z∗.

(2)

DY can be multi-dimensional and either continuous (i.e., regression) or categorical (i.e., classiﬁcation);
structured outputs such as strings, trees, and graphs can also be handled in our framework, as explainedlater.
2

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

5

That is, given the projected feature z ∗ , the (remaining) feature x is conditionally independent of output y and thus can be discarded without sacriﬁcing the predictabilityof y. Note that the conditionally independence is invariant against the choice of therepresentative W ∗ .
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},(3)drawn from a joint distribution with density pxy (x, y). The goal of SDR is, from data
Dn , to ﬁnd a projection matrix whose range agrees with that of W ∗ . For a projection
matrix W , we write
zi = W xi .
We assume that m is known throughout this paper.

2.2

Squared-Loss Mutual Information

A direct approach to SDR would be to determine W so that Eq.(2) is fulﬁlled. Let us
denote by z = W x for some projection matrix W . To this end, we adopt SMI as ourcriterion to be maximized with respect to W :)2∫ (1
pyz (y, z)
SMI(Y, Z) :=− 1 py (y)pz (z)dydz,(4)2
py (y)pz (z)where pyz (y, z) denotes the joint density of y and z, and py (y) and pz (z) denote themarginal densities of y and z, respectively. SMI(Y, Z) allows us to evaluate independence
between y and z since SMI(Y, Z) vanishes if and only if
pyz (y, z) = py (y)pz (z).
Note that Eq.(4) corresponds to the f -divergence (Ali & Silvey, 1966; Csiszár, 1967)
from pyz (y, z) to py (y)pz (z) with the squared loss, while ordinary MI corresponds to thef -divergence with the log loss, i.e., the Kullback-Leibler (KL) divergence (Kullback &Leibler, 1951):()∫
pyz (y, z)MI(Y, Z) := log
pyz (y, z)dydz.
py (y)pz (z)Thus SMI could be regarded as a natural alternative to ordinary MI.
The rationale behind the use of SMI in the context of SDR relies on the following
lemma:
Lemma 1. Decompose x into z and the component orthogonal to z as x = (z, z⊥ ), i.e.,
z is a member of the image of W and z⊥ is a member of the subspace perpendicular toSuﬃcient Dimension Reduction via Squared-loss Mutual Information Estimation

6

the image of W . Let pz⊥ y|z (z⊥ , y|z), pz⊥ |z (z⊥ |z), and py|z (y|z) be conditional densities.
Then we have
SMI(X, Y ) − SMI(Z, Y ))2∫ (pz⊥ y|z (z⊥ , y|z)1pyz (y, z)2 px (x)=1−dzdz⊥ dy2pz⊥ |z (z⊥ |z)py|z (y|z)
pz (z)2 py (y)≥ 0.
A proof of this lemma is given in Appendix A. Lemma 1 implies
SMI(X, Y ) ≥ SMI(Z, Y ),and the equality holds if and only ifpz⊥ y|z (z⊥ , y|z) = pz⊥ |z (z⊥ |z)py|z (y|z),which is equivalent to Eq.(2). Thus, Eq.(2) can be achieved by maximizing SMI(Z, Y )with respect to W ; then the ‘suﬃcient’ subspace can be identiﬁed.
Now we want to ﬁnd the projection matrix W that maximizes SMI(Z, Y ). However,
SMI is inaccessible in practice since densities pyz (y, z), py (y), and pz (z) are unknown.
Thus SMI needs to be estimated from data samples. Below, we introduce an SMI estimator.

2.3

SMI Approximation via Density-Ratio Estimation

Here, we consider a ﬁxed projection matrix W , and discuss the problem of approximating
SMI from samples. The convex duality (Boyd & Vandenberghe, 2004, p.91) gives thevariational representation (Keziou, 2003; Nguyen et al., 2010) of SMI as
1
SMI(Y, Z) = − inf J(g) − ,g2
where inf g is taken over all measurable functions on (DX × DY , BX × BY ), and∫∫12g(y, z) py (y)pz (z)dydz − g(y, z)pyz (y, z)dydz.
J(g) =
2This can be checked as follows: For f (u) = 21 (u2 − 1), we have)∫ (pyz (y, z)
SMI(Y, Z) = fpy (y)pz (z)dydz
py (y)pz (z)∫= sup pyz (y, z)g(y, z) − f ∗ (g(y, z)) py (y)pz (z)dydz
g∫1
1= sup pyz (y, z)g(y, z) − g 2 (y, z)py (y)pz (z)dydz −22g1= − inf J(g) − ,
g
2

(5)

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

7

where f ∗ is the convex conjugate of f that satisﬁes f (u) = supv∈R {uv − f ∗ (v)} (Boyd& Vandenberghe, 2004). Thus computing SMI is reduced to ﬁnding the minimizer g ∗ ofJ(g). We can show that g ∗ is given byg ∗ (y, z) :=pyz (y, z).
py (y)pz (z)(6)

Thus, estimating SMI(Y, Z) is reduced to estimating the above density ratio3 . We do notchoose a strategy to plug in density estimators of pyz , py , and pz into the formula (6).
This is because, in a region with small py (y)pz (z), small estimation error of pyz (y, z) isstrongly ampliﬁed. To avoid the unstable behavior around the tail, we directly modelthe density ratio g ∗ itself and impose regularization to control instability of density-ratio
estimators when needed.
Below, we consider parametric and non-parametric methods for estimating SMI.
2.3.1

Parametric Convergence Analysis

Let us consider the case where the function class G from which the function g is searched
is a parametric model:
G = {gθ (y, z) | θ ∈ Θ ⊂ Rb }.
Suppose that the true density-ratio g ∗ is contained in the model G, i.e., there exists θ ∗(∈ Θ) such that g ∗ = gθ∗ . Approximating the probability densities pyz (y, z), py (y), andpz (z) in Eq.(5) by their empirical distributions, we obtain the following optimization
problem.
[
]
n
n
∑
∑
1
1
θb := argmingθ (yi , zj )2 −gθ (yi , zi ) .
(7)2n2 i,j=1
n i=1
θ∈Θ
d
Then an SMI approximator SMI(Y,Z) can be constructed as
n
n
1∑
1 ∑
1
d
SMI(Y,Z) :=gθb(yi , zi ) − 2gθb(yi , zj )2 − .
n i=1
2n i,j=12(8)pb ∗∥ →
Suppose the standard regularity conditions for the consistency ∥θ−θ
0 is satisﬁed(see, for example, Section 3.2.1 of van der Vaart & Wellner, 1996). Let A and B be b × b
matrices deﬁned as

Aℓ,ℓ′ := Epx pz [∂ℓ gθ∗ (y, z)∂ℓ′ gθ∗ (y, z)],)[(Bℓ,ℓ′ := Epyz ∂ℓ gθ∗ (y, z) − Epz′ |y [∂ℓ gθ∗ (y, z ′ )] − 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 ′ )] ,3

This result can be generalized to a general f -divergence with small modiﬁcation (Keziou, 2003;
Nguyen et al., 2010).

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

8

where y ′ and z ′ are copies of y and z, and the partial derivative ∂ℓ is taken with respect
to the ℓ-th element θℓ of the parameter θ. Then we have the following theorem.
Theorem 1. Suppose that the matrix A is positive deﬁnite, then the SMI estimator (8)satisﬁesd
SMI(Y,Z) − SMI(Y, Z) = Op (n−1/2 ),(9)where Op denotes the asymptotic order in probability. Furthermore, we have1dtr(A−1 B) + o(n−1 ),EDn [SMI(Y,Z) − SMI(Y, Z)] =
2n

(10)

where EDn denotes the expectation over data samples Dn (see Eq.(3)).
A proof of Theorem 1 can be found in Appendix B. This theorem means that theabove SMI estimator retains the optimality in terms of the order of convergence in n,since Op (n−1/2 ) is the optimal convergence rate in the parametric setup (van der Vaart,2000).
2.3.2

Non-Parametric Convergence Analysis

Next, we consider non-parametric cases. Let the function class G be a general set offunctions on DY × DZ , where DZ = W DX . Let us consider a non-parametric version of
the empirical problem (cf. Eq.(7)):
[
]
n
n
∑
∑
1
1
λn
gb := argmin
g(yi , zj )2 −g(yi , zi ) + R(g)2 ,(11)
22n i,j=1
n i=1
2
g∈G
where R(g) is a non-negative regularization functional such thatsup[g(y, z)] ≤ R(g).

(12)

y,z

d
Then a non-parametric version of SMI approximator SMI(Y,Z) is given as
n
n
1
1∑
1 ∑
d
gb(yi , zj )2 − .
SMI(Y,Z) :=gb(yi , zi ) − 2n i=1
2n i,j=1
2

A useful example is to use a reproducing kernel Hilbert space (RKHS) (Aronszajn,1950) as G and the RKHS norm as R(g). Suppose G is an RKHS associated with boundedkernel k(·, ·):sup[k((y, z), (y, z))] ≤ C.
(13)
y,z

√
Let ∥ · ∥G denote the norm in the RKHS G. Then R(g) = C∥g∥G satisﬁes Eq.(12):√√g(y, z) = ⟨k((y, z), ·), g(·)⟩ ≤ k((y, z), (y, z))∥g∥G ≤ C∥g∥G ,Suﬃcient Dimension Reduction via Squared-loss Mutual Information Estimation

9

where we used the reproducing property of the kernel and Schwartz’s inequality. Notethat the Gaussian kernel satisﬁes this with C = 1. It is known that Gaussian kernel
RKHS spans a dense subset in the set of continuous functions. Another example of
RKHSs is Sobolev space. The canonical norm for this space is the integral of the squared
derivatives of functions. Thus the regularization term R(g) = ∥g∥G imposes the solution
to be smooth. The RKHS technique in Sobolev space has been well exploited in thecontext of spline models (Wahba, 1990). We intend that the regularization term R(g)is a generalization of the RKHS norm. Roughly speaking, R(g) is like a “norm” of thefunction space G.
We assume that the true density-ratio function g ∗ (y, z) is contained in the model G
and is bounded from above:g ∗ (y, z) ≤ M0 for all (y, z) ∈ DY × DZ .
Let GM be a ball of G with radius M > 0:
GM := {g ∈ G | R(g) ≤ M }.
To derive the convergence rate of our estimator, we utilize the bracketing entropy that is
a complexity measure of a function class (van der Vaart & Wellner, 1996, p. 83).
Deﬁnition 1. Given two functions l and u, the bracket [l, u] is the set of all functions
f with l(y, z) ≤ f (y, z) ≤ u(y, z) for all y and z. An ϵ-bracket is a bracket [l, u] with∥l − u∥L2 (py pz ) < ϵ. The bracketing entropy H[] (F, ϵ, L2 (py pz )) is the logarithm of theminimum number of ϵ-brackets needed to cover a function set F.
We assume that there exists γ (0 < γ < 2) such that, for all M > 0,(( )γ )M
H[] (GM , ϵ, L2 (py pz )) = O
.
ϵ

(14)

This quantity represents a complexity of function class G—the larger γ is, the more
complex the function class G is because, for larger γ, more brackets are needed to cover thefunction class. Gaussian RKHS satisﬁes this condition for arbitrarily small γ (Steinwart& Scovel, 2007). Note that when R(g) is the RKHS norm, the condition (14) holds for
all M > 0 if that holds for M = 1.
Then we have the following theorem.
2/(2+γ)), then we have
Theorem 2. Under the above setting, if λn → 0 and λ−1n = o(ndSMI(Y,Z) − SMI(Y, Z) = Op (max(λn , n−1/2 )).

(15)

A proof of Theorem 2 can be found in Appendix C. The conditions λn → 0 and λ−1n =o(n2/(2+γ) ) roughly mean that the regularization parameter λn should be suﬃciently small
but not too small. This theorem shows that the convergence rate of the non-parametricSuﬃcient Dimension Reduction via Squared-loss Mutual Information Estimation

10

version is also Op (n−1/2 ) if we take λn as n−2/(2+γ)+ϵ ≲ λn ≲ n−1/2 for suﬃciently small
ϵ. However, the non-parametric method requires a milder model assumption.
According to Nguyen et al. (2010) where a log-loss version of the above theorem
has been proven in the context of KL-divergence estimation, the above convergence rateachieves the optimal minimax rate under some setup. Thus the convergence property of
the above non-parametric method would also be optimal in the same sense.
As stated above, Gaussian RKHS satisﬁes the bracketing number condition (14) for
arbitrary γ > 0. Thus SMI with the Gaussian kernel achieves the convergence rate(15) with arbitrary γ > 0. However Gaussian RKHS is not suﬃciently rich to estimate afunction in Sobolev spaces with the optimal rate. To estimate a function in a Sobolev spacewith the optimal rate, we need to adjust the Gaussian width appropriately depending on
the sample size and the regularity of the Sobolev space. To analyze the convergence rate
for varying Gaussian widths, another technique than that used in this paper is required.
What we have shown in the theorem works only for a ﬁxed Gaussian width. To cope witha situation of varying Gaussian width, the techniques recently developed by Eberts andSteinwart (2011) are useful.
2.3.3

Practical Implementation

Here we describe a practical version of the above SMI approximation method, called
least-squares mutual information (LSMI) (Suzuki et al., 2009).
Let us restrict the search space of function g to some linear subspace G:
G := {α⊤ φ(y, z) | α = (α1 , . . . , αb )⊤ ∈ Rb },(16)

where α is a parameter to be learned from samples. φ(y, z) is a basis function such thatφ(y, z) := (φ1 (y, z), . . . , φb (y, z))⊤ ≥ 0b for all y, z,
where 0b is the b-dimensional vector with all zeros, and the inequality for vectors is appliedin the element-wise manner. φ(y, z) may be dependent on the samples {(yi , zi )}ni=1 , i.e.,kernel models are also allowed. Later in Section 2.3.5, we explain how the basis functions
φ(y, z) are designed.
Let us approximate the probability densities pyz (y, z), py (y), and pz (z) in Eq.(5) by
their empirical distributions. Then we have
[
]
1 ⊤c
λ ⊤
⊤
b
b := argminαα Hα − h α + α Rα ,
(17)
2
2
α∈Rb
where we included λα⊤ Rα (λ > 0) for regularization purposes, andn∑c := 1Hφ(yi , zj )φ(yi , zj )⊤ ,n2 i,j=1n∑b := 1hφ(yi , zi ).
n i=1

(18)

Note
∑n that when G is an RKHS corresponding to a kernel k, g(y, z) =ℓ=1 αℓ k((yℓ , zℓ ), (y, z)) is a member of the RKHS and the RKHS norm of g satisﬁes

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

11∥g∥2G = α⊤ Rα where R is the Gram matrix,√i.e., Rℓ,ℓ′ = k((yℓ , zℓ ), (yℓ′ , zℓ′ )). Thus the√√regularization term R(g) = C∥g∥G = C α⊤ Rα satisﬁes the condition (12) if thekernel function is bounded as Eq.(13).
Diﬀerentiating the objective function (17) with respect to α and equating it to zero,
we obtainbc + λR)α = h.
(H
(19)
Thus, the solution can be obtained just by solving the above system of linear equations.
b is given analytically as
The solution αbc + λR)−1 h.
b = (HαThen we can analytically approximate SMI as follows:
1 ⊤c
1
b ⊤α
d
b− α
b Hα
b− .
SMI(Y,Z) := h
2
2
2.3.4

(20)

Model Selection by Cross-Validation

As shown in Sections 2.3.1 and 2.3.2, our SMI estimator was shown to possess preferable
convergence properties. Nevertheless, its practical performance depends on the choice ofbasis 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 S = {(yi , zi )}ni=1 are divided into K disjoint subsets
b Sk is obtained using S\Sk{Sk }K
k=1 of (approximately) the same size. Then an estimator α(i.e., without Sk ) and the approximation error for the hold-out samples Sk is computed.
This procedure is repeated for k = 1, . . . , K, and its mean JbCV is outputted:
1 ∑JbCV :=
K k=1
K

(

)
1 ⊤c
⊤
b α
b Sk − hb HSk ααSk b Sk ,
2 Sk

(21)

b S denote H
b computed on the hold-out samples Sk in (18). For
cS and hc and hwhere Hkkmodel selection, we compute JbCV for all model candidates (the basis function φ(y, z) andthe regularization parameter λ), and choose the best model that minimizes JbCV . We canshow that JbCV is an almost unbiased estimator of the objective function J, where the‘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 may 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, we do not pursue this direction in this paper.
2.3.5

Design of Basis Functions

The above CV procedure would be useful when good candidates of basis functions areprepared. Here we propose to use the product kernel of the following form as basis

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

12

functions:φℓ (y, z) = ϕyℓ (y)ϕzℓ (z),b ℓ,ℓ′ is reduced from n2 to 2n:since the number of kernel evaluation when computing H
)( n
)
( n
∑ y
∑
1
b ℓ,ℓ′ =ϕ (yi )ϕyℓ′ (yi )ϕzℓ (zj )ϕzℓ′ (zj ) .
H
n2 i=1 ℓ
j=1
In the regression scenarios where y is continuous, we use the Gaussian kernel as the‘base’ kernels:()()∥y − uℓ ∥2∥z − vℓ ∥2yz
ϕℓ (y) := exp −, ϕℓ (z) := exp −
,
2σ 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 as follows: The density ratio (6) tends
to take large values if py (y)pz (z) is small and pyz (y, z) is large. When a non-negative
function is approximated by a Gaussian kernel model, many kernels may be needed in
the region where the output of the target function is large. On the other hand, only asmall number of kernels would be enough in the region where the output of the targetfunction is close to zero. Following this heuristic, we decided to allocate many kernels in
the regions where pyz (y, z) is large; this can be achieved by setting the Gaussian centers
at4 {(yi , zi )}ni=1 .
In the classiﬁcation scenarios where y is categorical, we use the delta kernel for y:ϕyℓ (y) := δ(y = uℓ ),cwhere δ(y = uℓ ) is 1 if y = uℓ and 0 otherwise. Note that, in this case, the matrix H
becomes block-diagonal, given that the samples are sorted according to the class labels.
Then the linear equation (19) can be solved eﬃciently.
More generally, when y is structured (e.g., strings, trees, and graphs), we may employ
kernels for structured data as ϕyℓ (y) (Lodhi et al., 2002; Collins & Duﬀy, 2002; Kashima& Koyanagi, 2002; Kondor & Laﬀerty, 2002; Kashima et al., 2003; Gärtner et al., 2003;
Gärtner, 2003).

2.4

Least-Squares Dimension Reduction

Finally, we show how the SMI approximator is employed for dimension reduction. To
ﬁnd a suﬃcient subspace, the dimension reduction problem is casted as an optimization
problem over the Grassmann manifold Grdm (R) (see Eq.(1)).
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 computationalcost.
4

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

13

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 show that thegradient is given by
b⊤
c 3
d
∂h
∂H
∂R b
∂ SMI
b −α
b + λαb − β)b⊤b − β)
b⊤
b
=
(2α
( α
(β − α),∂Wℓ,ℓ′ ∂Wℓ,ℓ′∂Wℓ,ℓ′ 2∂Wℓ,ℓ′
c +λR)−1 H
cαb and we used the formula ∂X∂t = −X −1 ∂X
where βb := (HX −1 for a positive∂tsymmetric matrix X each element of which is a function of t.
d
SMI
In the Euclidean space, the ordinary gradient ∂∂Wgives the steepest direction. However, on a manifold, the natural gradient (Amari, 1998) gives the steepest direction. ThedSMIdnatural gradient ∇SMI(W) at W is the projection of the ordinary gradient ∂∂W
to thetangent space of Grm
(R)
at
W
.
If
the
tangent
space
isequippedwiththecanonical
metricd⟨G1 , G2 ⟩ = 12 tr(G⊤1 G2 ), the natural gradient is given as follows (Edelman et al., 1998):−1d∇SMI(W
)=

d
d
d
∂ SMI
∂ SMI
∂ SMI
−
W ⊤W =W⊥⊤ W⊥ ,
∂W
∂W
∂W

(22)

where W⊥ is any (d − m) × d matrix such that [W ⊤ W⊥⊤ ] is orthogonal. Then the geodesicdfrom W to the direction of the natural gradient ∇SMI(W) over Grm
d (R) can be expressedusing t (∈ R) as
( [
]) [
]
d
∂ SMI
⊤
[
]
Om
W
W
⊥
∂W
Wt := Id Od−m exp t,d⊤W⊥−W⊥ ∂ SMI
Od−m
∂W

where ‘exp’ for a matrix denotes the matrix exponential, and Ob is the b × b zero matrix(note that the derivative ∂t Wt |t=0 coincides with the natural gradient (22), see Edelman
et al. (1998) for detailed derivation of the geodesic). Thus line search along the geodesicin the natural gradient direction is equivalent to ﬁnding the maximizer from {Wt | t ≥ 0}.
For choosing the step size of each gradient update, we may use some approximate line
search method such as Armijo’s rule (Patriksson, 1999, p.50) or backtracking line search(Boyd & Vandenberghe, 2004, p.464). In our setting, Armijo’s rule ﬁnds the step size as
the maximum tk that satisﬁes⊤d∂ SMI(W ) 
d
d
 ∂Wt
SMI(W
tk ) − SMI(W ) ≥ tk µ tr
∂t
∂W
t=0
[
]
⊤ d
d
∂
SMI(W
)
∂
SMI(W
)
= tk µ tr W⊥⊤ W⊥,∂W∂W
where t0 = 1, tk = αtk−1 (k = 1, 2, . . . ), α ∈ (0, 1), and µ ∈ (0, 1) are given parameters.
We call the proposed dimension reduction algorithm Least-Squares Dimension Reduction (LSDR). The entire algorithm is summarized in Figure 1. In practice, we performed
CV once in several steps because executing CV at every step is computationally expensive.

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

1.
2.
3.
4.

14

Initialize projection matrix W .
Optimize Gaussian width σ and regularization parameter λ by CV.
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 LSDR algorithm.
.

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 Gramf + ϵIb , where K
f is the kernel Grammatrix as the regularization matrix R: R = K
yf for avoiding
e ℓ,ℓ′ := ϕ (uℓ′ )ϕz (vℓ′ ). ϵIb is added to K
matrix for the chosen centers: Kℓℓ
non-degeneracy; we set ϵ = 0.01. We ﬁx the number of basis functions to b = min(100, n),and choose the Gaussian width σ and the regularization parameter λ based on 5-fold CVwith grid search. We restart the natural gradient search 10 times with random initialpoints, and choose the one having the minimum CV score (21).

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):(a) Linear dependence: d = 5, m = 1. y has a linear dependence on x as
y = x(1) + ϵ,
where x(k) denotes the k-th element of x, x ∼ N (x; 0, I5 ) and ϵ ∼ N (y; 0, 0.25).
Here N (x; µ, Σ) denotes the normal density with mean µ and variance-covariancematrix Σ. The optimal projection is given by
W ∗ = [1 0 0 0 0] .

(23)(b) Non-linear dependence 1: d = 5, m = 1. y has a quadratic dependence on x as
y = (x(1) )2 + ϵ,
where x ∼ N (x; 0, I5 ) and ϵ ∼ N (y; 0, 1). The optimal projection is given by
Eq.(23).

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

15(c) Non-linear dependence 2: d = 5, m = 1. y has a lattice-structured dependence
on x asx ∼ U (x; [−0.5, 0.5]5 ),{N (y; 0, 0.25)
if x(1) ≤ | 16 |,y|x ∼1
N (y; 1, 0.25) + 12 N (y; −1, 0.25) otherwise,2
where U(x; S) denotes the uniform density on a set S. The optimal projection is
given by Eq.(23).
(d) Fukumizu et al. (2009): d = 4, m = 2. y has a non-linear dependence on x asy=x(1)+ (1 + x(2) )2 + 0.4ϵ,
0.5 + (x(2) + 1.5)2

wherex ∼ )N (x; 0, I4 ) and ϵ ∼ N (1; 0, 1). The optimal projection is W ∗ =(1 0 0 0.
0 1 0 0(e) Fukumizu et al. (2009): d = 4, m = 1. y has a non-linear dependence on x as
y = sin2 (πx(1) + 1) + 0.4ϵ,
where x ∼ U(x; [0, 1]4 \ {x ∈ R4 | x(i) ≤ 0.7 (i = 1, . . . , 4)}) and ϵ ∼ N (1; 0, 1). The
optimal projection is W ∗ = [1 0 0 0].
(f ) Fukumizu et al. (2009): d = 10, m = 1. y has a non-linear dependence on x as
1
y = (x(1) − 1)2 ϵ,
2where x ∼ N (x; 0, I10 ) and ϵ ∼ N (1; 0, 1). The optimal projection is W ∗ =[1 0 · · · 0].
Let us compare the proposed LSDR with the following methods.
• Kernel dimension reduction (KDR) (Fukumizu et al., 2009),• The Hilbert-Schmidt independence criterion (HSIC) (Gretton et al., 2005),• Sliced inverse regression (SIR) (Li, 1991),• Sliced average variance estimation (SAVE) (Cook, 2000),• Principal Hessian direction (pHd) (Li, 1992).
• Minimum Average (conditional) Variance Estimation (dMAVE) (Xia, 2007).

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

0

6

−1

1

y

8

y

10

1
y

2

16

4

0

2

−2

−1

0

−3

−2
−2

−1

x1

0

1

2

−2

(a) Linear

0
x1

−2

2

(b) Non-linear dependence 1

−0.5

0
x1

0.5

(c) Non-linear dependence 2

4
20

1

10

0.5

0

y

y

y

2

0

2
0
x2

−2

0

−2
x1

(d) Fukumizu et al. (2009)

2

0
−2

−0.5

−4
0

0.2

0.4

x1

0.6

0.8

−2

1

(e) Fukumizu et al. (2009)

−1

0
x1

1

2

(f) Fukumizu et al. (2009)

Figure 2: Artiﬁcial datasets.
In KDR and HSIC, the Gaussian width is set to the median sample distance, following
the suggestions in the original papers (Gretton et al., 2005; Fukumizu et al., 2009). We
used the dimension reduction package dr included in R for SIR, SAVE, and pHd. The
parameters for these methods such as the number of slices were set to be the default
values. To execute dMAVE, we used the publicly available code5 . The principal directions
estimated by SIR, SAVE, pHd, and dMAVE do not necessarily form an orthogonal system,
i.e., if we let F be the matrix each row of which corresponds to each principal direction,then F is not necessarily a projection matrix. To recover a projection matrix W , weperformed singular decomposition of F as F = V SU , and set W = U .
We evaluate the performance of each method by
c − W ∗⊤ W ∗
c ⊤W
W

,

(24)

Frobenius

c is an estimated projection matrix, andwhere ∥ · ∥Frobenius denotes the Frobenius norm, W∗
W is an optimal projection matrix.
The performance of each method is summarized in Table 2, which depicts the meanand standard deviation of the Frobenius-norm error (24) over 50 trials when the number ofsamples is n = 100. LSDR overall shows good performance; in particular, it performs thebest for datasets (b) and (c). KDR also tends to work reasonably well, but it sometimes
performs poorly; this seems to be caused by an inappropriate choice of the Gaussian kernel
5

http://www.stat.nus.edu.sg/˜staxyc/dMAVE.m

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

17

Table 2: Mean (and standard deviation in the bracket) of the Frobenius-norm error (24)for toy datasets. The best method in terms of the mean error and comparable ones basedon the one-sided t-test at the signiﬁcance level 1% are indicated by boldface.
Data d m
(a)
5 1
(b)
5 1
(c)
5 1
(d)
4 2
(e)
4 1
(f) 10 1

LSDR.13(.04).15(.06).10(.05).20(.14).09(.06).35(.12)KDR.13(.05).25(.21).44(.32).16(.06).13(.06).40(.12)HSIC
SIR.17(.07) .11(.05).44(.36) .83(.19).68(.32) .89(.14).18(.08) .30(.15).16(.07) .21(.10).49(.17) .68(.22)

SAVE.37(.27).31(.11).48(.20).44(.18).34(.19).91(.13)pHd.89(.12).24(.07).86(.12).50(.18).36(.14).83(.12)dMAVE.13(.04).20(.09).25(.28).10(.05).08(.04).26(.06)

width, implying that the heuristic of using the median sample distance as the kernel widthis not always appropriate. On the other hand, LSDR with CV performs stably well forvarious types of datasets. dMAVE also works well, and is competitive to LSDR for these
artiﬁcial datasets.

3.2

Classiﬁcation for Benchmark Datasets

Finally, we evaluate the classiﬁcation performance after dimension reduction for several
benchmark datasets. We use ‘image’, ‘waveform’, ‘pima-indians-diabetes’, and ‘letter
recognition’ in the UCI repository6 . We randomly choose 200 samples from the dataset,and apply LSDR, KDR, HSIC, and dMAVE to obtain projections onto low-dimension
subspaces with m = ⌈d/4⌉, ⌈d/2⌉, ⌈3d/4⌉, where ⌈c⌉ denotes the smallest integer not
smaller than c. Then we train the support vector machine (Schölkopf & Smola, 2002) on
the 200 projected training samples.
The misclassiﬁcation rate is computed for samples not used for training. Table 3
summarizes the mean and 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). Our contributions were parametric and non-parametricanalyses of the rate of convergence of the SMI approximator, and the proposal of a dimension reduction algorithm based on the SMI approximator. The proposed method isadvantageous in several respects, e.g., density estimation is not involved, it is distributionfree, and model selection by cross-validation is available. The eﬀectiveness of the proposed
method over existing methods was shown through experiments.

6

http://www.ics.uci.edu/~mlearn/MLRepository.html.

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

18

Table 3: Mean (and standard deviation in the bracket) of misclassiﬁcation rates for benchmark datasets. The best method in terms of the mean error and comparable ones basedon the one-sided t-test at the signiﬁcance level 1% are indicated by boldface.
Data set

d
18
image
18
18
21
waveform
21
21
8
pima
8
8
16
letter (a,b, & c) 16
16

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

LSDR
.083(.019).088(.022).093(.018).130(.014).119(.013).116(.007).249(.022).260(.016).244(.020).031(.009).026(.008).016(.006)

KDR.125(.038).106(.026).091(.019).127(.008).135(.010).131(.008).247(.024).250(.021).243(.019).028(.012).017(.007).014(.006)

HSIC
dMAVE
.158(.044).501(.134).115(.035).468(.130).095(.023) .468(.130).160(.016).183(.016).163(.016).184(.015).159(.014).182(.021).252(.020) .257(.017).252(.017) .265(.027).251(.021) .252(.019).035(.014) .018(.008).020(.006) .016(.006).017(.008) .013(.008)

Acknowledgements
T.S. was supported in part by MEXT Kakenhi 22700289 and the Aihara Project, the
FIRST program from JSPS, initiated by CSTP. M.S. was supported by AOARD, SCAT,
and the JST PRESTO program.

A

Proof of Lemma 1

Let x = (z, z⊥ ). By the relations dx = dzdz⊥ and px (x) = pz⊥ |z (z⊥ |z)pz (z), we have∫pxy (x, y)pyz (y, z)
px (x)py (y)dxdy
px (x)py (y)py (y)pz (z)∫pz⊥ |zy (z⊥ |z, y)pyz (y, z)pyz (y, z)dzdz⊥ dy=py (y)pz (z))2∫∫ (pyz (y, z)2
pyz (y, z)=dzdy =py (y)pz (z)dzdy.
py (y)pz (z)py (y)pz (z)Thus we obtain
SMI(X, Y ) − SMI(Z, Y )
)2
)2
∫ (
∫ (
1
pxy (x, y)1pyz (y, z)=
px (x)py (y)dxdy −pz (z)py (y)dzdy
2px (x)py (y)2pz (z)py (y))2∫ (1
pxy (x, y)pyz (y, z)=−
px (x)py (x)dxdy.
2px (x)py (y) pz (z)py (y)Suﬃcient Dimension Reduction via Squared-loss Mutual Information Estimation

19

Noticing thatpxy (x, y)pxy (x, y)pyz (y, z)
pxy (x, y)pyz (y, z)==
px (x)py (y)px (x)py (y)pyz (y, z)
pz⊥ |z (z⊥ |z)pz (z)py (y)py|z (y|z)pz (z)pz⊥ y|z (z⊥ , y|z)pyz (y, z)=,
pz⊥ |z (z⊥ |z)py|z (y|z) py (y)pz (z)we have
1
SMI(X, Y ) − SMI(Z, Y ) =
2

∫ (
1−

pz⊥ y|z (z⊥ , y|z)pz⊥ |z (z⊥ |z)py|z (y|z))2pyz (y, z)2 px (x)dzdz⊥ dy,pz (z)2 py (y)which concludes the proof of Lemma 1.

B

Proof of Theorem 1

en , P, Pn as
For notational simplicity, we deﬁne linear operators Q, Qn , Q∑n∑i,j=1 f (yi , zj )1≤i̸=j≤n f (yi , zj )e
Qf := Epy pz f, Qn f :=, Qn :=,
2nn(n − 1)∑nf (yi , zi ).
P f := Epyz f, Pn f := i=1nLet ∇gθ′ denote ∇θ gθ |θ=θ′ = (∂ℓ gθ |θ=θ′ )ℓ for θ ′ ∈ Θ. Since θb is the optimizer of the
problem (7), we have
)
(
1
1
2
Qn gθb − Pn gθb = Qn ∇(gθ2b) − Pn ∇gθb.
(25)
0=∇
2
2
Therefore, as in the standard asymptotic expansion for maximum likelihood estimators(van der Vaart, 2000), we have
(
)
1
1
2
⊤ 2
⊤
∗ ⊤
b
Qn ∇∇ (gθ∗ ) − Pn ∇∇ gθ∗
0 = Qn ∇(gθ∗ ) − Pn ∇gθ∗ + (θ − θ )22)(+ Op ∥θb − θ ∗ ∥2 .
This implies(

1
Qn ∇∇⊤ (gθ2∗ ) − Pn ∇∇⊤ gθ∗2)(+ Op ∥θb − θ ∗ ∥2 .

θb − θ ∗ = −

)−1 (

1
Qn ∇(gθ2∗ ) − Pn ∇gθ∗2

Since Q (gθ∗ ∇gθ∗ ) − P ∇gθ∗ = 0, we have1
Qn ∇(gθ2∗ ) − Pn ∇gθ∗ = Op (n−1/2 ),
2

)

(26)

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

20

and
1
Qn ∇∇⊤ (gθ2∗ ) − Pn ∇∇⊤ gθ∗2
1= Q∇∇⊤ (gθ2∗ ) − P ∇∇⊤ gθ∗ + Op (n−1/2 )
2(
)
)
(
)
(
= Q ∇gθ∗ ∇⊤ gθ∗ + Q gθ∗ ∇∇⊤ gθ∗ − P ∇∇⊤ gθ∗ + Op (n−1/2 )
)
)
(
)
(
(
= Q ∇gθ∗ ∇⊤ gθ∗ + P ∇∇⊤ gθ∗ − P ∇∇⊤ gθ∗ + Op (n−1/2 )()= Q ∇gθ∗ ∇⊤ gθ∗ + Op (n−1/2 ).
Thus Eq.(26) implies(27)

θb − θ ∗ = Op (n−1/2 ) + Op (∥θb − θ ∗ ∥2 ),in particular, θb − θ ∗ = Op (n−1/2 ). Moreover Eq.(26) becomes
(
)
(
)−1 1
∗
⊤
2
Qn ∇(gθ∗ ) − Pn ∇gθ∗ + Op (n−1 ).
θb − θ = −Q ∇gθ∗ ∇ gθ∗
2(28)

Eqs.(5) and (8) give

(
)
1
1 2
2
d
SMI(Y, Z) − SMI(Y, Z) = − Qn gθb + Pn gθb − − Qgθ∗ + P gθ∗22
1
1= − Qn gθ2b + Pn gθb + Qn gθ2∗ − Pn gθ∗
2
2
(
)
1 2
1
2
− − Qgθ∗ + P gθ∗ + Qn gθ∗ − Pn gθ∗ .
2
2

(29)

The ﬁrst four terms of the RHS can be expanded as
1
1− Qn gθ2b + Pn gθb + Qn gθ2∗ − Pn gθ∗
2
( 2
)
1
2
∗
⊤
b
= (θ − θ) ∇Qn gθb − Pn gθb
2()
1
1 ∗ b⊤⊤
2b + Op (n−3/2 )+ (θ − θ) ∇∇Qn gθb − Pn gθb (θ ∗ − θ)
2
2
)
(
1
1 ∗ b⊤⊤
2b + Op (n−3/2 )Qn gθb − Pn gθb (θ ∗ − θ)
= (θ − θ) ∇∇2
2)(
1 ∗ b⊤
1⊤
2
b + Op (n−3/2 )= (θ − θ) ∇∇Qn gθ∗ − Pn gθ∗ (θ ∗ − θ)22)
1 ∗ b⊤ (b + Op (n−3/2 )= (θ − θ) Q ∇gθ∗ ∇⊤ gθ∗ (θ ∗ − θ)2
= Op (n−1 ).

(∵ Eq.(25))(∵ Eq.(27))(30)

On the other hand, by the central limit theorem,
11− Qgθ2∗ + P gθ∗ + Qn gθ2∗ − Pn gθ∗ = Op (n−1/2 ).
2
2

(31)

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

21

Substituting Eqs.(30) and (31) into Eq.(29), we have the ﬁrst assertion (9).
Next we prove the second assertion (10). Based on Eqs.(29) and (30), we can evaluatedthe expectation of SMI(Y,Z) − SMI(Y, Z) asdEDn [SMI(Y,Z) − SMI(Y, Z)]
[
]
1
1
2
2
= EDn − Qn gθb + Pn gθb + Qn gθ∗ − Pn gθ∗
2
2
{
[
]}
1
= tr AEDn (θb − θ ∗ )(θb − θ ∗ )⊤ + o(n−1 ),2where we used the fact that)(Q ∇gθ∗ ∇⊤ gθ∗ = A.

Below, we will show that[] 1
EDn (θb − θ ∗ )(θb − θ ∗ )⊤ = A−1 BA−1 + o(n−1 ).
n(32)Obviously this gives Eq.(10).
Eq.(28) implies[]
EDn (θb − θ ∗ )(θb − θ ∗ )⊤
)(
)]
[(
1
1
2
⊤ 2
⊤
−1
Qn ∇(gθ∗ ) − Pn ∇gθ∗Qn ∇ (gθ∗ ) − Pn ∇ gθ∗= A ED n
A−1 + o(n−1 )
2
2
[
(
)]
−1
= A EDn (Qn (gθ∗ ∇gθ∗ ) − Pn ∇gθ∗ ) Qn (gθ∗ ∇⊤ gθ∗ ) − Pn ∇⊤ gθ∗ A−1 + o(n−1 )
[(
)(
)]
−1
⊤
⊤
e
e
= A EDn Qn (gθ∗ ∇gθ∗ ) − Pn ∇gθ∗ Qn (gθ∗ ∇ gθ∗ ) − Pn ∇ gθ∗ A−1 + o(n−1 ).
(33)Leth(w, w′ ) =)1(gθ∗ (y, z ′ )∇gθ∗ (y, z ′ ) + gθ∗ (y ′ , z)∇gθ∗ (y ′ , z) − ∇gθ∗ (y, z) − ∇gθ∗ (y ′ , z ′ )2

for w = (y, z) and w′ = (y ′ , z ′ ). Then h is symmetric (i.e., h(w, w′ ) = h(w′ , w)), hasmean 0 (i.e., Epw pw′ [h(w, w′ )] = 0), and satisﬁes
en (gθ∗ ∇gθ∗ ) − Pn ∇gθ∗ =Q∑
1h(wi , wj ) =: Un .
n(n − 1) 1≤i̸=j≤n

en (gθ∗ ∇gθ∗ ) − Pn ∇gθ∗ is a U -statistic with the symmetric kernel h(wi , wj ).
Therefore, Q
It is known (Theorem 7.1 in Hoeﬀding, 1948) that the variance of U -statistic is given by[] 4
EDn Un Un⊤ = E[h(w1 , w2 )h(w1 , w2′ )⊤ ] + o(n−1 ),n(34)

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

22

where w1 = (y1 , z1 ), w2 = (y2 , z2 ), and w2′ = (y2′ , z2′ ) are i.i.d. variables with probability
density pyz . For notational simplicity, we write
f1 (y, z) := gθ∗ (y, z)∇gθ∗ (y, z),Then
4E[h(w1 , w2 )h(w1 , w2′ )⊤ ] = E

f2 (y, z) := ∇gθ∗ (y, z).

[

f1 (y1 , z2 )f1 (y1 , z2′ )⊤+ f1 (y1 , z2 )f1 (y2′ , z1 )⊤− f1 (y1 , z2 )f2 (y1 , z1 )⊤− f1 (y1 , z2 )f2 (y2′ , z2′ )⊤
|
{z
} |
{z
} |
{z
} |
{z
}
(A)

(B)

(C)

(D)

+ f1 (y2 , z1 )f1 (y1 , z2′ )⊤+ f1 (y2 , z1 )f1 (y2′ , z1 )⊤− f1 (y2 , z1 )f2 (y1 , z1 )⊤− f1 (y2 , z1 )f2 (y2′ , z2′ )⊤
|

{z

} |

(E)

{z

(F )

} |

{z

} |

(G)

{z

}

(H)

− f2 (y1 , z1 )f1 (y1 , z2′ )⊤− f2 (y1 , z1 )f1 (y2′ , z1 )⊤+ f2 (y1 , z1 )f2 (y1 , z1 )⊤+ f2 (y1 , z1 )f2 (y2′ , z2′ )⊤
|
{z
} |
{z
} |
{z
} |
{z
}
(I)

(J)

(K)

(L)

(M )

(N )

(O)

(P )

]
− f2 (y2 , z2 )f1 (y1 , z2′ )⊤− f2 (y2 , z2 )f1 (y2′ , z1 )⊤+ f2 (y2 , z2 )f2 (y1 , z1 )⊤+ f2 (y2 , z2 )f2 (y2′ , z2′ )⊤ .
|
{z
} |
{z
} |
{z
} |
{z
}
Going through simple calculations, each term in the RHS of the above equation can beevaluated as(A), (C), (I) = Epz|y pz′ |y py [∇gθ∗ (y, z)∇⊤ gθ∗ (y, z ′ )],(B) = Epz′ |y py′ |z pyz [∇gθ∗ (y, z ′ )∇⊤ gθ∗ (y ′ , z)],(D), (H), (L), (M), (N), (O), (P) = Epyz [∇gθ∗ (y, z)]Epyz [∇⊤ gθ∗ (y, z)],(E) = Epz′ |y py′ |z pyz [∇gθ∗ (y ′ , z)∇⊤ gθ∗ (y, z ′ )],(F), (G), (J) = Epy|z py′ |z pz [∇gθ∗ (y, z)∇⊤ gθ∗ (y ′ , z)],(K) = Epyz [∇gθ∗ (y, z)∇⊤ gθ∗ (y, z)].
Therefore,
4E[h(w1 , w2 )h(w1 , w2′ )⊤ ]= Epyz [∇gθ∗ (y, z)∇⊤ gθ∗ (y, z)] − Epz|y pz′ |y py [∇gθ∗ (y, z)∇⊤ gθ∗ (y, z ′ )]+ Epz′ |y py′ |z pyz [∇gθ∗ (y, z ′ )∇⊤ gθ∗ (y ′ , z)] + Epz′ |y py′ |z pyz [∇gθ∗ (y ′ , z)∇⊤ gθ∗ (y, z ′ )]− Epy|z py′ |z pz [∇gθ∗ (y, z)∇⊤ gθ∗ (y ′ , z)] − Epzy [∇gθ∗ (y, z)]Epzy [∇⊤ gθ∗ (y, z)]
)
[(
′
′
′
′
= Epyz ∇gθ∗ (y, z) − Epz′ |y [∇gθ∗ (y, z )] − 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 ′ )]= B.
This with Eqs.(33) and (34) yields Eq.(32), and thus we have the second assertion (10).

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

C

23

Proof of Theorem 2

Before proving Theorem 2, we show the following auxiliary lemma.
2/(2+γ)Lemma 2. Under the setting of Theorem 2, if λn → 0 and λ−1), then wen = o(nhaveR(bg ) = Op (1),(35)∥bg − g ∗ ∥2 = Op (λ1/2n ),where ∥ · ∥2 means the L2 (px py )-norm and Op denotes the asymptotic order in probability.
Proof of Lemma 2. Hoeﬀding (1963) derived Bernstein’s inequality for double sum, i.e.,
(
)
(
)
x
x2
e
√
P |(Qn − Q)f | > √≤ 2 exp −n
4(∥f ∥∞ x/(3 n) + Qf 2 )for x ≥ 0. By applying the above inequality to the proof of Theorem 2 in Birgé and
Massart (1993) or Theorem 5.11 in van de Geer (2000) instead of Bernstein’s inequality
for i.i.d. sum, we obtain a “double√ sum” version of those theorems. That is, exponential
decay of the tail probability of n supf ∈F |(Qn − Q)f |, where F is a class of uniformly
bounded measurable functions and satisﬁes a polynomial bracketing condition (14) as GM(see Theorem 3 in the supplementary material). Later, this is used to obtain Eqs.(37)
and (38).
From the deﬁnition, we obtain⇒

1
1
Qn gb2 − Pn gb + λn R(bg )2 ≤ Qn g ∗ 2 − Pn w + λn R(g ∗ )2
2
2
1
∗ 2
∗ ∗
Qn (b
g − g ) − Qn (g (g − gb)) − Pn (bg − g ∗ ) + λn (R(bg )2 − R(g ∗ )2 ) ≤ 0.
2

On the other hand, Q(g ∗ (g ∗ − gb)) = P (g ∗ − gb) indicates
1(Q − Qn )(bg − g ∗ )2 − (Q − Qn )(g ∗ (g ∗ − gb)) − (P − Pn )(bg − g ∗ ) − λn (R(bg )2 − R(g ∗ )2 )
2
1
≥ Q(b
g − g ∗ )2 .
2Therefore, to bound ∥b
g − g ∗ ∥2 , it suﬃces to bound the LHS of the above inequality.
2
Deﬁne FM and FM as
FM := {g − g ∗ | g ∈ GM },2:= {f 2 | f ∈ FM }.
FM

For g, g ′ ∈ GM such that ∥g − g ′ ∥2 ≤ δ, the L2 (px py )-norm of the diﬀerence between(g − g ∗ )2 and (g ′ − g ∗ )2 is bounded by∥(g − g ∗ )2 − (g ′ − g ∗ )2 ∥2 = ∥(g − g ′ )(g + g ′ − 2g ∗ )∥2 ≤ 2δ(M + M0 ).
2Thus, the bracketing entropy of FM
has the following order:(()γ )(M + M0 )22
H[] (FM , δ, L2 (px py )) = Oδas δ → 0.

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

Let f := gb − g ∗ . Then, as in Lemma 5.14 and Theorem 10.6 in van de Geer (2000), we
obtain
(
)
γ
2
1
− 2+γ
2
2 1− 222 γ2
22|(Qn − Q)(f )| = Op √ ∥f ∥2 (1 + R(bg ) + M0 ) ∨ n(1 + R(bg ) + M0 ) ,n(36)
where a ∨ b = max(a, b). Here, we have used the double sum version of Theorem 5.11 invan de Geer (2000), which is needed to obtain the same formula as Lemma 5.14 in van deGeer (2000). Since√∥f 2 ∥2 ≤ ∥f ∥2

2(1 + R(bg )2 + M02 ),the RHS of Eq.(36) is further bounded by
(
)
2
1
1− γ2
− 2+γ
2
2
2 12 + γ422|(Qn − Q)(f )| = Op √ ∥f ∥2 (1 + R(bg ) + M0 )∨n
(1 + R(bg ) + M0 ) .
n(37)Similarly, we can show that
(
)
2
γ
1
1− γ2
−
∗ ∗
|(Qn − Q)(g (g − gb))| = Op √ ∥f ∥2 (1 + R(bg ) + M0 ) 2 ∨ n 2+γ (1 + R(bg ) + M0 ) .
n(38)′Note that, for g, g ∈ G,∫∫
′ 2
(g − g ) dpxy = (g − g ′ )2 wdpx dpy ≤ M0 ∥g − g ′ ∥22 ,which implies((
H[] (FM , δ, L2 (pxy )) = O

M + M0
δ

)γ )
.

Thus,
(
∗

|(Pn − P )(g − gb)| = Op

)
2
γ
1
1− γ2
−
√ ∥f ∥2 (1 + R(bg ) + M0 ) 2 ∨ n 2+γ (1 + R(bg ) + M0 ) . (39)n

Combining Eqs.(37), (38), and (39), we can bound the L2 (px py )-norm of f as1g )2
∥f ∥22 + λn R(b2(∗ 2≤ λn R(g ) + Op

)
2
1
1− γ2
− 2+γ
2
2 21 + γ422√ ∥f ∥2 (1 + R(bg ) + M0 )∨n
(1 + R(bg ) + M0 ) .
n

The rest can be proved by following a similar line to Theorem 10.6 in van de Geer (2000).
We redeﬁne M0 ← max{R(g ∗ ), M0 }, and deﬁne J = 1 + R(bg )2 + M02 . There are two
possible situations, namely, (a) R(bg ) ≥ M0 + 1 or (b) R(bg ) < M0 + 1. We show thestochastic order of ∥f ∥2 and R(bg ) for the two situations separately.

Suﬃcient Dimension Reduction via Squared-loss Mutual Information Estimation(a) (i) If ∥f ∥2 ≥ n−1/(2+γ) J, then1
∥f ∥22 + λn R(bg )2 ≤ λn M02 + Op
2

(

)
1
1−γ/2 1/2+γ/4√ ∥f ∥2J.
n

In this case1g )2 ≤ 2λn M02∥f ∥22 + λn R(b
2

or

25

1
g )2 ≤ Op∥f ∥22 + λn R(b
2

(

)
1
1−γ/2 1/2+γ/4√ ∥f ∥2J.
n

For the ﬁrst event, we have∥f ∥22 ≤ 4λn M02 = O(λn ),

R(b
g) ≤

√
2M0 = Op (1).

On the other hand, for the second event, we have
)
)
(
(
1
1
1
1−γ/2 1/2+γ/4
1−γ/2
2
2
1+γ/2
∥f ∥2 + λn R(bg ) ≤ Op √ ∥f ∥2
J≤ Op √ ∥f ∥2
R(b
g)
,
2
n
n
which indicates1∥f ∥22 ≤ Op
2
and

(

1
√ ∥f ∥1−γ/2R(bg )1+γ/2
2
n

)

(
)
∥f ∥2 ≤ Op n−1/(2+γ) R(b
g) ,

⇒

)
1
1−γ/2
1+γ/2
λn R(bg ) ≤ Op √ ∥f ∥2
R(b
g)
n
)
(
g)
−2/(2−γ)2/(2+γ)−1/(2−γ) R(b⇒ R(bg ) ≤ λnOp n= op (R(bg )) (∵ λ−1)).
n = o(nn1/(2+γ)(

2

Thus, the probability of this event goes to 0.
(ii) If ∥f ∥2 ≤ n−1/(2+γ) J, then()
1
∥f ∥22 + λn R(bg )2 ≤ λn M02 + Op n−2/(2+γ) R(bg )2 .
2
In this case,)(g )2 .
g )2 ≤ Op n−2/(2+γ) R(bλn R(bg )2 ≤ 2λn M02 or λn R(bThen, similarly to the above case, we can show that the probability of the second event
goes to 0. On the other hand, for the ﬁrst event, ∥f ∥2 ≤ n−1/(2+γ) J ≤ n−1/(2+γ) (1+2M02 ) =√1/2
Op (λn ) and R(bg ) ≤ 2M0 = Op (1).
(b) In this situation, R(bg ) < M0 + 1 = Op (1).
−1/(2+γ)(i) If ∥f ∥2 ≥ n
J, then1g )2 ≤ λn M02 + Op
∥f ∥22 + λn R(b
2

(

)
1
1−γ/2
√ ∥f ∥2.
n

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

26

In this case,
(
)
(
)
1
1
1−γ/2
2
−1/(2+γ)∥f ∥2 ≤ Op λn ∨ √ ∥f ∥2⇒ ∥f ∥2 ≤ Op λ1/2
= Op (λ1/2
n ∨n
n ).
2
n
(ii) If ∥f ∥2 ≤ n−1/(2+γ) J, it is obvious that ∥f ∥2 = Op (λn ).
Consequently, the proof of Lemma 2 was completed.
1/2

Based on Lemma 2, we prove Theorem 2 below.
As in the proof of Lemma 2, let f := gb − g ∗ . Since Q(f g ∗ ) = P f , we have1
1
Qn gb2 − Pn gb − ( Qg ∗ 2 − P g ∗ )
2
2
(
)
1
1 ∗2
∗ 2
∗
∗
= Qn (f + g ) − Pn (f + g ) −Qg − P g2
2
1
1= Qf 2 + (Qn − Q)f 2 + (Qn − Q)(g ∗ f ) − (Pn − P )f2
2
1∗2
+ (Qn g − Qg ∗ 2 ) − (Pn g ∗ − P g ∗ ).
2

(40)

Below, we show that each term of the RHS of the above equation is Op (λn ). By thecentral limit theorem, we have1
(Qn g ∗ 2 − Qn g ∗ 2 ) − (Pn g ∗ − P g ∗ ) = Op (n−1/2 ).
2
Since Lemma 2 gives ∥f ∥22 = Op (λn ) and R(bg ) = Op (1), Eqs.(37), (38), and (39) in the
proof of Lemma 2 imply
)
(
2
1
1− γ2
− 2+γ
2
≤ Op (λn ),|(Qn − Q)f | = Op √ ∥f ∥2 ∨ nn|(Qn − Q)(g ∗ f )| = Op (λn ),|(Pn − P )g ∗ | = Op (λn ).
2/(2+γ)) was used. Lemma 2 also implies
In the above derivation, λ−1n = o(nQf 2 = ∥f ∥22 = Op (λn ).
Combining these inequalities with Eq.(40) implies
1
1d
SMI(Y,Z) − SMI(Y, Z) = − Qn gb2 + Pn gb − (− Qg ∗ 2 + P g ∗ )2
2= Op (λn + n−1/2 ) = Op (max(λn , n−1/2 )),which concludes the proof of Theorem 2.

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

27

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 in learning. Neural Computation, 10,
251–276.
Aronszajn, N. (1950). Theory of reproducing kernels. Transactions of the American
Mathematical Society, 68, 337–404.
Birgé, L., & Massart, P. (1993). Rates of convergence for minimum contrast estimators.
Probability Theory and Related Fields, 97, 113–150.
Boyd, S., & Vandenberghe, L. (2004). Convex optimization. Cambridge, UK: Cambridge
University Press.
Breiman, L., & Friedman, J. H. (1985). Estimating optimal transformations for multiple
regression and correlation. Journal of the American Statistical Association, 80, 580–598.
Bura, E., & Cook, R. D. (2001). Extending sliced inverse regression. Journal of theAmerican Statistical Association, 96, 996–1003.
Chiaromonte, F., & Cook, R. D. (2002). Suﬃcient dimension reduction and graphics inregression. Annals of the Institute of Statistical Mathematics, 54, 768–795.
Collins, M., & Duﬀy, N. (2002). Convolution kernels for natural language. Advances in
Neural Information Processing Systems 14. Cambridge, MA: MIT Press.
Cook, R. D. (1998a). Principal Hessian directions revisited. Journal of the American
Statistical Association, 93, 84–100.
Cook, R. D. (1998b). Regression graphics: Ideas for studying regressions through graphics.
New York, NY, USA: Wiley.
Cook, R. D. (2000). SAVE: A method for dimension reduction and graphics in regression.
Communications in Statistics - Theory and Methods, 29, 2109–2121.
Cook, R. D., & Ni, L. (2005). Suﬃcient dimension reduction via inverse regression.
Journal of the American Statistical Association, 100, 410–428.
Cover, T. M., & Thomas, J. A. (2006). Elements of information theory. Hoboken, NJ,
USA: John Wiley & Sons, Inc. 2nd edition.

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

28

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 the information by an adaptive
partitioning of the observation space. IEEE Transactions on Information Theory, 45,
1315–1321.
Durand, J., & Sabatier, R. (1997). Additive splines for partial least squares regression.
Journal of the American Statistical Association, 92, 1546–1554.
Eberts, M., & Steinwart, I. (2011). Optimal learning rates for least squares svms using
gaussian kernels. Advances in Neural Information Processing Systems 24 (pp. 1539–
1547).
Edelman, A., Arias, T. A., & Smith, S. T. (1998). The geometry of algorithms withorthogonality constraints. SIAM J. Matrix Anal. Appl., 20, 303–353.
Fukumizu, K., Bach, F. R., & Jordan, M. (2009). Kernel dimension reduction in regression. The Annals of Statistics, 37, 1871–1905.
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.
Gärtner, T. (2003). A survey of kernels for structured data. ACM SIGKDD Explorations
Newsletter, 5, 49–58.
Gärtner, T., Flach, P., & Wrobel, S. (2003). On graph kernels: Hardness results andeﬃcient alternatives. Proceedings of the Sixteenth Annual Conference on Computational
Learning Theory.
Goutis, C., & Fearn, T. (1996). Partial least squares regression on smooth factors. Journalof the American Statistical Association, 91, 627–632.
Gretton, A., Bousquet, O., Smola, A., & Schölkopf, B. (2005). Measuring statisticaldependence with Hilbert-Schmidt norms. Algorithmic Learning Theory (pp. 63–77).
Berlin: Springer-Verlag.
Hoeﬀding, W. (1948). A class of statistics with asymptotically normal distribution. The
Annals of Mathematical Statistics, 19, 293–325.
Hoeﬀding, W. (1963). Probability inequalities for sums of bounded random variables.
Journal of the American Statistical Association, 58, 13–30.
Hotelling, H. (1936). Relations between two sets of cariants. Biometrika, 28, 321–377.
Hulle, M. M. V. (2005). Edgeworth approximation of multivariate diﬀerential entropy.
Neural Computation, 17, 1903–1910.

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

29

Kashima, H., & Koyanagi, T. (2002). Kernels for semi-structured data. Proceedingsof the Nineteenth International Conference on Machine Learning (pp. 291–298). San
Francisco, CA: Morgan Kaufmann.
Kashima, H., Tsuda, K., & Inokuchi, A. (2003). Marginalized kernels between labeledgraphs. Proceedings of the 20th International Conference on Machine Learning. San
Francisco, CA.
Keziou, A. (2003). Dual representation of φ-divergences and applications. C. R. Acad.
Sci. Paris, Ser. I, 336, 857–862.
Kondor, R. I., & Laﬀerty, J. (2002). Diﬀusion kernels on graphs and other discrete input
spaces. Proceedings of the Nineteenth International Conference on Machine Learning(pp. 315–322).
Kraskov, A., Stögbauer, H., & Grassberger, P. (2004). Estimating mutual information.
Physical Review E, 69, 066138.
Kullback, S., & Leibler, R. A. (1951). On information and suﬃciency. Annals of Mathematical Statistics, 22, 79–86.
Li, K. C. (1991). Sliced inverse regression for dimension reduction. Journal of American
Statistical Association, 86, 316–342.
Li, K. C. (1992). On principal Hessian directions for data visualization and dimension
reduction: Another application of Stein’s lemma. Journal of the American Statistical
Association, 87, 1025–1039.
Li, K. C., Lue, H. H., & Chen, C. H. (2000). Interactive tree-structured regression via
principal Hessian directions. Journal of the American Statistical Association, 95, 547–
560.
Lodhi, H., Saunders, C., Shawe-Taylor, J., Cristianini, N., & Watkins, C. (2002). Textclassiﬁcation using string kernels. Journal of Machine Learning Research, 2, 419–444.
Nguyen, X., Wainwright, M. J., & Jordan, M. I. (2010). Estimating divergence functionals
and the likelihood ratio by convex risk minimization. IEEE Transactions on Information
Theory. to appear.
Patriksson, M. (1999). Nonlinear programming and variational inequality problems.
Dredrecht: Kluwer Academic.
Reiss, P. T., & Ogden, R. T. (2007). Functional principal component regression andfunctional partial least squares. Journal of the American Statistical Association, 102,
984–996.
Schölkopf, B., & Smola, A. J. (2002). Learning with kernels. Cambridge, MA: MIT Press.

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

30

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. New Challenges for Feature
Selection in Data Mining and Knowledge Discovery (pp. 5–20).
Torkkola, K. (2003). Feature extraction by non-parametric mutual information maximization. Journal of Machine Learning Research, 3, 1415–1438.
van de Geer, S. (2000). Empirical processes in M-estimation. Cambridge University Press.
van der Vaart, A. W. (2000). Asymptotic statistics. Cambridge University Press.
van der Vaart, A. W., & Wellner, J. A. (1996). Weak convergence and empirical processes:
With applications to statistics. New York: Springer.
Wahba, G. (1990). Spline model for observational data. Philadelphia and Pennsylvania:
Society for Industrial and Applied Mathematics.
Wold, H. (1966). Estimation of principal components and related models by iterative
least squares. In P. R. Krishnaiah (Ed.), Multivariate analysis, 391–420. New York,NY, USA: Academic Press.
Xia, Y. (2007). A constructive approach to the estimation of dimension reduction directions. The Annals of Statistics, 35, 2654–2690.
Zhu, L., Miao, B., & Peng, H. (2006). On sliced inverse regression with high-dimensional
covariates. Journal of the American Statistical Association, 101, 630–643.