# Kernel Sliced Inverse Regression: Regularization and Consistency

Authors: Qiang Wu; Feng Liang; Sayan Mukherjee
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

Kernel sliced inverse regression (KSIR) is a natural framework for nonlinear dimension reduction using the mapping induced by kernels. However, there are numeric, algorithmic, and conceptual subtleties in making the method robust and consistent. We apply two types of regularization in this framework to address computational stability and generalization performance. We also provide an interpretation of the algorithm and prove consistency. The utility of this approach is illustrated on simulated and real data.


---

## Body

## 1. Introduction

The goal of dimension reduction in the standard regression/classification setting is to summarize the information in thep-dimensional predictor variable X relevant to predicting the univariate response variable Y. The summary S(X) should have d≪p variates and ideally should satisfy the following conditional independence property:
(1)Y⫫X∣S(X).
Thus, any inference of  Y involves only the summary statistic S(X) which is of much lower dimension than the original data X.Linear methods for dimension reduction focus on linear summaries of the data,S(X)=(β1TX,…,βdTX). The d-dimensional subspace,  𝒮=span(β1,…,βd), is defined as the effective dimension reduction (e.d.r.) space in [1], since 𝒮 summarizes all the predictive information on Y. A key result in [1] is that under some mild conditions, the e.d.r. directions {βj}j=1d correspond to the eigenvectors of the matrix:
(2)T=[cov(X)]-1cov[𝔼(X∣Y)].
Thus, the e.d.r. directions or subspace can be estimated via an eigenanalysis of the matrix T, which is the foundation of the sliced inverse regression (SIR) algorithm proposed in [1, 2]. Further developments include sliced average variance estimation (SAVE)  [3] and Principal Hessian directions (PHDs)  [4]. The aforementioned algorithms cannot be applied in high-dimensional settings, where the number of covariates p is greater than the number of observations n, since the sample covariance matrix is singular. Recently, an extension of SIR has been proposed in [5], which can handle the case for p>n based on the idea of partial least squares.A common premise held in high-dimensional data analysis is that the intrinsic structure of data is in fact low dimensional, for example, the data is concentrated on a manifold. Linear methods such as SIR often fail to capture this nonlinear low-dimensional structure. However, there may exist a nonlinear embedding of the data into a Hilbert space, where a linear method can capture the low-dimensional structure. The basic idea in applying kernel methods is the application of a linear algorithm to the data mapped into a feature space induced by a kernel function. If projections onto this low-dimensional structure can be computed by inner products in this Hilbert space, the so-called kernel trick  [6, 7] can be used to obtain simple and efficient algorithms. Since the embedding is nonlinear, linear directions in the feature space correspond to nonlinear directions in the original data space. Nonlinear extensions of some classical linear dimensional reduction methods using this approach include kernel principle component analysis  [6], kernel Fisher discriminant analysis  [8], and kernel independent correlation analysis  [9]. This idea was applied to SIR in [10, 11] resulting in the kernel sliced inverse regression (KSIR) method which allows for the estimation of nonlinear e.d.r. directions.There are numeric, algorithmic, and conceptual subtleties to a direct application of this kernel idea to SIR, although it looks quite natural at first glance. In KSIR, thep-dimensional data are projected into a Hilbert space ℋ through a feature map: ϕ:X→ℋ and the nonlinear features are supposed to be recovered by the eigenfunctions of the following operator
(3)T=[cov(ϕ(X))]-1cov[𝔼(ϕ(X)∣Y)].
However, this operator T is actually not well defined in general, especially when ℋ is infinite dimensional and the covariance operator cov(ϕ(X)) is not invertible. In addition, the key utility of representation theorems in kernel methods is that optimization in a possibly infinite dimensional Hilbert space ℋ reduces to solving a finite dimensional optimization problem. In the KSIR algorithm developed in [10, 11], the representer theorem has been implicitly used and seems to work well in empirical studies. However, it is not theoretically justifiable since T is estimated empirically via observed data. Moreover, the computation of eigenfunctions of an empirical estimate of T from observations is often ill-conditioned and results in computational instability. In [11], a low rank approximation of the kernel matrix is used to overcome this instability and to reduce computational complexity. In this paper, our aim is to clarify the theoretical subtleties in KSIR and to motivate two types of regularization schemes to overcome the computational difficulties arising in KSIR. The consistency is proven and practical advantages of regularization are demonstrated via empirical experiments.


## 2. Mercer Kernels and Nonlinear e.d.r. Directions

The extension of SIR to use kernels is based on properties of reproducing kernel Hilbert spaces (RKHSs) and in particular Mercer kernels  [12].Given predictor variablesX∈𝒳⊆ℝp, a Mercer kernel is a continuous, positive, and semidefinite function k(·,·):𝒳×𝒳→ℝ with the following spectral decomposition:
(4)k(x,z)=∑jλjϕj(x)ϕj(z),
where {ϕj} are the eigenfunctions and {λj} are the corresponding nonnegative, nonincreasing eigenvalues. An important property of Mercer kernels is that each kernel k uniquely corresponds to an RKHS as follows:
(5)ℋ={f∣f(x)=∑j∈Λajϕj(x)with∑j∈Λaj2λj<∞},
where the cardinality of Λ:={j:λj>0} is the dimension of the RKHS which may be infinite  [12, 13].Given a Mercer kernel, there exists a unique map or embeddingϕ from 𝒳 to a Hilbert space defined by the eigenvalues and eigenfunctions of the kernel. The map takes the following form:
(6)ϕ(x)=(λ1ϕ1(x),λ2ϕ2(x),…,λ|Λ|ϕ|Λ|(x)).
The Hilbert space induced by this map with the standard inner product k(x,z)=〈ϕ(x),ϕ(z)〉 is isomorphic to the RKHS (5), and we will denote both Hilbert spaces as ℋ  [13]. In the case where k is infinite dimensional, ϕ:𝒳→ℓ2.The random variableX∈𝒳 induces a random element ϕ(X) in the RKHS. Throughout this paper, we will use Hilbert space valued random variables; so we now recall some basic facts. Let Z be a random element in ℋ with 𝔼∥Z∥<∞, where ∥·∥ denotes the norm in ℋ induced by its inner product 〈·,·〉. The expectation 𝔼(Z) is defined to be an element in ℋ, satisfying 〈a,𝔼(Z)〉=𝔼〈a,Z〉, for all a∈ℋ. If 𝔼∥Z∥2≤∞, then the covariance operator of Z is defined as 𝔼[(Z-𝔼Z)⊗(Z-𝔼Z)], where
(7)(a⊗b)f=〈b,f〉a,foranyf∈ℋ.Let𝒫 denote the measure for random variable X. Throughout, we assume the following conditions.Assumption 1.
(1)
For allx∈𝒳,   k(x,·) is 𝒫-measurable.
(2)
There existsM>0 such that x∈𝒳,   k(X,X)≤M (a.s.) with respect to 𝒫.Under Assumption1, the random element ϕ(X) has a well-defined mean and a covariance operator because ∥ϕ(x)∥2=k(x,x) is bounded (a.s.). Without loss of generality, we assume 𝔼ϕ(X)=0, where 0 is the zero element in ℋ. The boundedness also implies that the covariance operator Σ=𝔼[ϕ(X)⊗ϕ(X)] is compact and has the following spectral decomposition:
(8)Σ=∑i=1∞wiei⊗ei,
where wi and ei∈ℋ are the eigenvalues and eigenfunctions, respectively.We assume the following model for the relationship betweenY and X:
(9)Y=F(〈β1,ϕ(X)〉,…,〈βd,ϕ(X)〉,ε),
with βj∈ℋ and the distribution of ε is independent of X. This model implies that the response variable Y depends on X only through a d-dimensional summary statistic
(10)S(X)=(〈β1,ϕ(X)〉,…,〈βd,ϕ(X)〉).
Although S(X) is a linear summary statistic in ℋ, it extracts nonlinear features in the space of the original predictor variables X. We call  {βj}j=1d  the nonlinear e.d.r. directions, and 𝒮=span(β1,…,βd) the nonlinear e.d.r. space. The following proposition  [10] extends the theoretical foundation of SIR to this nonlinear setting.Proposition 2.
Assume the following linear design condition forℋ that for any f∈ℋ, there exists a vector b∈ℝd such that
(11)𝔼[〈f,ϕ(X)〉∣S(X)]=bTS(X),withS(X)=(〈β1,ϕ(X)〉,…,〈βd,ϕ(X)〉)T.
Then, for the model specified in (9), the inverse regression curve 𝔼[ϕ(X)∣Y] is contained in the span of (Σβ1,…,Σβd), where Σ is the covariance operator of ϕ(X).Proposition2 is a straightforward extension of the multivariate case in [1] to a Hilbert space or a direct application of the functional SIR setting in [14]. Although the linear design condition (11) may be difficult to check in practice, it has been shown that such a condition usually holds approximately in a high-dimensional space  [15]. This conforms to the argument in [11] that the linearity in a reproducing kernel Hilbert space is less strict than the linearity in the Euclidean space. Moreover, it is pointed out in [1, 10] that even if the linear design condition is violated, the bias of SIR variate is usually not large.An immediate consequence of Proposition2 is that nonlinear e.d.r. directions are the eigenvectors corresponding to the largest eigenvalues of the following generalized eigendecomposition problem:
(12)Γβ=λΣβ,whereΣ=cov[ϕ(X)],Γ=cov[𝔼(ϕ(X)∣Y)],
or equivalently from an eigenanalysis of the operator T=Σ-1Γ. In the infinite dimensional case, a technical difficulty arises since the operator
(13)Σ-1=∑i=1∞‍wi-1ei⊗ei
is not defined on the entire Hilbert space ℋ. So for the operator T to be well defined, we need to show that the range of Γ is indeed in the range of Σ. A similar issue also arose in the analysis of dimension reduction and canonical analysis for functional data  [16, 17]. In these analyses, extra conditions are needed for operators like T to be well defined. In KSIR, this issue is resolved automatically by the linear design condition and extra conditions are not required as stated by the following Theorem; see Appendix A for the proof.Theorem 3.
Under Assumption1 and the linear design condition (11), the following hold:  (i)
the operatorΓ is of finite rank  dΓ≤d. Consequently, it is compact and has the following spectral decomposition:
(14)Γ=∑i=1dΓ‍τiui⊗ui,whereτi and ui are the eigenvalues and eigenvectors, respectively. Moreover,  ui∈range(Σ) for all  i=1,…,dΓ;(ii)
the eigendecomposition problem (12) is equivalent to the eigenanalysis of the operator T, which takes the following form
(15)T=∑i=1dΓτiui⊗Σ-1(ui).


## 3. Regularized Kernel Sliced Inverse Regression

The discussion in Section2 implies that nonlinear e.d.r. directions can be retrieved by applying the original SIR algorithm in the feature space induced by the Mercer kernel. There are some computational challenges to this idea such as estimating an infinite dimensional covariance operator and the fact that the feature map is often difficult or impossible to compute for many kernels. We address these issues by working with inner products of the feature map and adding a regularization term to kernel SIR.


### 3.1. Estimating the Nonlinear e.d.r. Directions

Givenn observations {(x1,y1),…,(xn,yn)}, our objective is to obtain an estimate of the e.d.r. directions (β^1,…,β^d). We first formulate a procedure almost identical to the standard SIR procedure except that it operates in the feature space ℋ. This highlights the immediate relation between the SIR and KSIR procedures.(1)
Without loss of generality, we assume that the mapped predictor variables are mean zero,  that is,∑i=1nϕ(xi)=0,  for otherwise we can subtract ϕ-=(1/n)∑i=1nϕ(xi) from ϕ(xi). The sample covariance is estimated by
(16)Σ^=1n∑i=1nϕ(xi)⊗ϕ(xi).(2)
Bin theY variables into H slices G1,…,GH and compute mean vectors of the corresponding mapped predictor variables for each slice
(17)ψh=1nh∑i∈Ghϕ(xi),h=1,…,H.Compute the sample between-group covariance matrix(18)Γ^=∑h=1Hnhnψh⊗ψh.(3)
Estimate the SIR directionsβ^j by solving the generalized eigendecomposition problem:
(19)Γ^β=λΣ^β.This procedure is computationally impossible if the RKHS is infinite dimensional or the feature map cannot be computed (which is the usual case). However, the model given in (9) requires not the e.d.r. directions but the projection onto these directions, that is, the d summary statistics
(20)v1=〈β1,ϕ(x)〉,…,vd=〈βd,ϕ(x)〉,
which we call the KSIR variates. Like other kernel algorithms, the kernel trick enables KSIR variates to be efficiently computed from only the kernel k, not the map ϕ.The key quantity in this alternative formulation is the centred Gram matrixK defined by the kernel function  k(·,·), where
(21)Kij=〈ϕ(xi),ϕ(xj)〉=k(xi,xj).
Note that the rank of K is less than n; so K is always singular.Given the centered Gram matrixK, the following generalized eigendecomposition problem can be used to compute the KSIR variates:
(22)KJKc=λK2c,
where c denotes the n-dimensional generalized eigenvector and J denotes an n×n matrix with Jij=1/nm if i,j are in the mth group consisting of nm observations and zero otherwise. The following proposition establishes the equivalence between two eigendecomposition problems, (22) and (19), in the recovery of KSIR variates (v1,…,vd).Proposition 4.
Given the observations{(x1,y1),…,(xn,yn)}, let (β^1,…,β^d) and (c^1,…,c^d) denote the generalized eigenvectors of (22) and (19), respectively. Then, for any x∈𝒳 and j=1,…,d, the following holds:
(23)vj=〈β^j,ϕ(x)〉=c^jTKx,Kx=(k(x,x1),…,k(x,xn))T,
provided that Σ^  is invertible. When  Σ^  is not invertible, the conclusion holds modulo the null space of  Σ^.This result was proven in [10], in which the algorithm was further reduced to solving
(24)JKc=λKc
by canceling K from both sides of (22).Remark 5.
It is important to remark that whenΣ^ is not invertible, Proposition 4 states that the equivalence between (19) and (22) holds modulo the null space of Σ^ that is a requirement of the representer theorem, that is, β^j are linear combinations of ϕ(xi). Without this mandated condition, we will see that each eigenvector of (22) produces an eigenvector of (19), while the eigenvector of (19) is not necessarily recovered by (22).Remark 6.
It is necessary to clarify the difference between (12) and (19). In (12), it is natural to assume that β is orthogonal to the null space of  Σ. To see this, let  β=β0+β⋆ with β0 belonging to the null space and  β⋆ orthogonal to the null space. Then,
(25)𝔼[〈β,ϕ(x)〉2]=〈β,Σβ〉=〈β⋆,Σβ⋆〉=𝔼[〈β⋆,ϕ(x)〉2],
that is, 〈β,ϕ(x)〉=〈β⋆,ϕ(x)〉 (a.s.). It means that β0 does not contribute to the KSIR variates and thus could be set as 0. However, in (19), if β^ is an eigenvector and β^⋆ is its orthogonal component relative to the null space of Σ^, the identity 〈β^,ϕ(x)〉=〈β^⋆,ϕ(x)〉 is in general not true for a new point x which is different from the observations. Thus, from a theoretical perspective, it is not as natural to assume the representer theorem, although it works well in practice. In this sense, the KSIR algorithm based on (22) or (24) does not have a thorough mathematical foundation.


### 3.2. Regularization and Stability

Except for the theoretical subtleties, in applications with relatively small samples, the eigendecomposition in (22) is often ill-conditioned resulting in overfitting as well as numerically unstable estimates of the e.d.r. space. This can be addressed by either thresholding eigenvalues of the estimated covariance matrix Σ^ or by adding a regularization term to (22) or (24).We motivate two types of regularization schemes. The first one is the traditional ridge regularization. It is used in both linear SIR and functional SIR  [18–20], which solves the eigendecomposition problem
(26)Γ^β=λ(Σ^+sI)β.
Here, and in the sequel, I denotes the identity operator and s is a tuning parameter. Assuming the representer theorem, its kernel form is given as
(27)JKc=λ(K+nsI)c.Another type of regularization is to regularize (22) directly:
(28)KJKc=λ(K2+n2sI)c.
The following proposition, whose proof is in Appendix B, states that solving the generalized eigendecomposition problem (28) is equivalent to finding the eigenvectors of
(29)(Σ^2+sI)-1Σ^Γ^.Proposition 7.
Letc^j be the eigenvectors of (28), and let β^j be the eigenvectors of (29). Then, the following holds for the regularized KSIR variates:
(30)v^j=c^jT[k(x,x1),…,k(x,xn)]=〈β^j,ϕ(x)〉.This algorithm is termed as the Tikhonov regularization. For linear SIR, it is shown in [21] that the Tikhonov regularization is more efficient than the ridge regularization.Except for the computational stability, regularization also makes the matrix forms of KSIR, (27) and (28), interpretable by justifiable representer theorem.Proposition 8.
For both ridge and the Tikhonov regularization scheme of KSIR, the eigenfunctionsβ^j are linear combinations of  ϕ(xi),i=1,…,n.The conclusion follows from the observation thatβ^j=(1/λj)(Γ^β^j-Σ^β^j) for the ridge regularization and β^j=(1/λj)(Σ^Γ^β^j-Σ^2β^j) for the Tikhonov regularization.To close, we remark that KSIR is computationally advantageous even for the case of linear models whenp≫n due to the fact that the eigendecomposition problem is for n×n matrices rather than the p×p matrices in the standard SIR formulation.


### 3.3. Consistency of Regularized KSIR

In this subsection, we prove the asymptotic consistency of the e.d.r. directions estimated by regularized KSIR and provide conditions under which the rate of convergence isOp(n-1/4). An important observation from the proof is that the rate of convergence of the e.d.r. directions depends on the contribution of the small principal components. The rate can be arbitrarily slow if the e.d.r. space depends heavily on eigenvectors corresponding to small eigenvalues of the covariance operator.Note that various consistency results are available for linear SIR  [22–24]. These results hold only for the finite dimensional setting and cannot be adapted to KSIR where the RKHS is often infinite dimensional. Consistency of functional SIR has also been studied before. In [14], a thresholding method is considered, which selects a finite number of eigenvectors and uses results from finite rank operators. Their proof of consistency requires stronger and more complicated conditions than ours. The consistency for functional SIR with ridge regularization is proven in [19], but it is of a weaker form than our result. We remark that the consistency results for functional SIR can be improved using a similar argument in this paper.In the following, we state the consistency results for the Tikhonov regularization. A similar result can be proved for the ridge regularization while the details are omitted.Theorem 9.
Assume𝔼k(X,X)2<∞, limn→∞s(n)=0 and  limn→∞sn=∞; then
(31)|〈β^j,ϕ(·)〉-〈βj,ϕ(·)〉|=op(1),j=1,…,dΓ,
where dΓ is the rank of Γ, 〈βj,ϕ(·)〉 is the projection onto the jth e.d.r., and 〈β^j,ϕ(·)〉 is the projection onto the jth e.d.r. as estimated by regularized KSIR.
If the e.d.r. directions{βj}j=1dΓ depend only on a finite number of eigenvectors of the covariance operator Σ, the rate of convergence is O(n-1/4).This theorem is a direct corollary of the following theorem which is proven in AppendixC.Theorem 10.
Define the projection operator and its complement for eachN≥1(32)ΠN=∑i=1N‍ei⊗ei,ΠN⊥=I-ΠN=∑i=N+1∞‍ei⊗ei,
where {ei}i=1∞ are the eigenvectors of the covariance operator Σ as defined in (8), with the corresponding eigenvalues denoted by wi.
Assume𝔼k(X,X)2<∞. For each N≥1, the following holds:
(33)∥(Σ^2+sI)-1Σ^Γ^-T∥HS=Op(1sn)+∑j=1dΓ‍(swN2∥ΠN(u~j)∥+∥ΠN⊥(u~j)∥),
where u~j=Σ-1uj and {uj}i=1dΓ are the eigenvectors of Γ as defined in (14) and ∥·∥HS denotes the Hilbert-Schmidt norm of a linear operator.
Ifs=s(n) satisfy s→0 and sn→∞ as n→∞, then
(34)∥(Σ^2+sI)-1Σ^Γ^-T∥HS=op(1).


## 3.1. Estimating the Nonlinear e.d.r. Directions

Givenn observations {(x1,y1),…,(xn,yn)}, our objective is to obtain an estimate of the e.d.r. directions (β^1,…,β^d). We first formulate a procedure almost identical to the standard SIR procedure except that it operates in the feature space ℋ. This highlights the immediate relation between the SIR and KSIR procedures.(1)
Without loss of generality, we assume that the mapped predictor variables are mean zero,  that is,∑i=1nϕ(xi)=0,  for otherwise we can subtract ϕ-=(1/n)∑i=1nϕ(xi) from ϕ(xi). The sample covariance is estimated by
(16)Σ^=1n∑i=1nϕ(xi)⊗ϕ(xi).(2)
Bin theY variables into H slices G1,…,GH and compute mean vectors of the corresponding mapped predictor variables for each slice
(17)ψh=1nh∑i∈Ghϕ(xi),h=1,…,H.Compute the sample between-group covariance matrix(18)Γ^=∑h=1Hnhnψh⊗ψh.(3)
Estimate the SIR directionsβ^j by solving the generalized eigendecomposition problem:
(19)Γ^β=λΣ^β.This procedure is computationally impossible if the RKHS is infinite dimensional or the feature map cannot be computed (which is the usual case). However, the model given in (9) requires not the e.d.r. directions but the projection onto these directions, that is, the d summary statistics
(20)v1=〈β1,ϕ(x)〉,…,vd=〈βd,ϕ(x)〉,
which we call the KSIR variates. Like other kernel algorithms, the kernel trick enables KSIR variates to be efficiently computed from only the kernel k, not the map ϕ.The key quantity in this alternative formulation is the centred Gram matrixK defined by the kernel function  k(·,·), where
(21)Kij=〈ϕ(xi),ϕ(xj)〉=k(xi,xj).
Note that the rank of K is less than n; so K is always singular.Given the centered Gram matrixK, the following generalized eigendecomposition problem can be used to compute the KSIR variates:
(22)KJKc=λK2c,
where c denotes the n-dimensional generalized eigenvector and J denotes an n×n matrix with Jij=1/nm if i,j are in the mth group consisting of nm observations and zero otherwise. The following proposition establishes the equivalence between two eigendecomposition problems, (22) and (19), in the recovery of KSIR variates (v1,…,vd).Proposition 4.
Given the observations{(x1,y1),…,(xn,yn)}, let (β^1,…,β^d) and (c^1,…,c^d) denote the generalized eigenvectors of (22) and (19), respectively. Then, for any x∈𝒳 and j=1,…,d, the following holds:
(23)vj=〈β^j,ϕ(x)〉=c^jTKx,Kx=(k(x,x1),…,k(x,xn))T,
provided that Σ^  is invertible. When  Σ^  is not invertible, the conclusion holds modulo the null space of  Σ^.This result was proven in [10], in which the algorithm was further reduced to solving
(24)JKc=λKc
by canceling K from both sides of (22).Remark 5.
It is important to remark that whenΣ^ is not invertible, Proposition 4 states that the equivalence between (19) and (22) holds modulo the null space of Σ^ that is a requirement of the representer theorem, that is, β^j are linear combinations of ϕ(xi). Without this mandated condition, we will see that each eigenvector of (22) produces an eigenvector of (19), while the eigenvector of (19) is not necessarily recovered by (22).Remark 6.
It is necessary to clarify the difference between (12) and (19). In (12), it is natural to assume that β is orthogonal to the null space of  Σ. To see this, let  β=β0+β⋆ with β0 belonging to the null space and  β⋆ orthogonal to the null space. Then,
(25)𝔼[〈β,ϕ(x)〉2]=〈β,Σβ〉=〈β⋆,Σβ⋆〉=𝔼[〈β⋆,ϕ(x)〉2],
that is, 〈β,ϕ(x)〉=〈β⋆,ϕ(x)〉 (a.s.). It means that β0 does not contribute to the KSIR variates and thus could be set as 0. However, in (19), if β^ is an eigenvector and β^⋆ is its orthogonal component relative to the null space of Σ^, the identity 〈β^,ϕ(x)〉=〈β^⋆,ϕ(x)〉 is in general not true for a new point x which is different from the observations. Thus, from a theoretical perspective, it is not as natural to assume the representer theorem, although it works well in practice. In this sense, the KSIR algorithm based on (22) or (24) does not have a thorough mathematical foundation.


## 3.2. Regularization and Stability

Except for the theoretical subtleties, in applications with relatively small samples, the eigendecomposition in (22) is often ill-conditioned resulting in overfitting as well as numerically unstable estimates of the e.d.r. space. This can be addressed by either thresholding eigenvalues of the estimated covariance matrix Σ^ or by adding a regularization term to (22) or (24).We motivate two types of regularization schemes. The first one is the traditional ridge regularization. It is used in both linear SIR and functional SIR  [18–20], which solves the eigendecomposition problem
(26)Γ^β=λ(Σ^+sI)β.
Here, and in the sequel, I denotes the identity operator and s is a tuning parameter. Assuming the representer theorem, its kernel form is given as
(27)JKc=λ(K+nsI)c.Another type of regularization is to regularize (22) directly:
(28)KJKc=λ(K2+n2sI)c.
The following proposition, whose proof is in Appendix B, states that solving the generalized eigendecomposition problem (28) is equivalent to finding the eigenvectors of
(29)(Σ^2+sI)-1Σ^Γ^.Proposition 7.
Letc^j be the eigenvectors of (28), and let β^j be the eigenvectors of (29). Then, the following holds for the regularized KSIR variates:
(30)v^j=c^jT[k(x,x1),…,k(x,xn)]=〈β^j,ϕ(x)〉.This algorithm is termed as the Tikhonov regularization. For linear SIR, it is shown in [21] that the Tikhonov regularization is more efficient than the ridge regularization.Except for the computational stability, regularization also makes the matrix forms of KSIR, (27) and (28), interpretable by justifiable representer theorem.Proposition 8.
For both ridge and the Tikhonov regularization scheme of KSIR, the eigenfunctionsβ^j are linear combinations of  ϕ(xi),i=1,…,n.The conclusion follows from the observation thatβ^j=(1/λj)(Γ^β^j-Σ^β^j) for the ridge regularization and β^j=(1/λj)(Σ^Γ^β^j-Σ^2β^j) for the Tikhonov regularization.To close, we remark that KSIR is computationally advantageous even for the case of linear models whenp≫n due to the fact that the eigendecomposition problem is for n×n matrices rather than the p×p matrices in the standard SIR formulation.


## 3.3. Consistency of Regularized KSIR

In this subsection, we prove the asymptotic consistency of the e.d.r. directions estimated by regularized KSIR and provide conditions under which the rate of convergence isOp(n-1/4). An important observation from the proof is that the rate of convergence of the e.d.r. directions depends on the contribution of the small principal components. The rate can be arbitrarily slow if the e.d.r. space depends heavily on eigenvectors corresponding to small eigenvalues of the covariance operator.Note that various consistency results are available for linear SIR  [22–24]. These results hold only for the finite dimensional setting and cannot be adapted to KSIR where the RKHS is often infinite dimensional. Consistency of functional SIR has also been studied before. In [14], a thresholding method is considered, which selects a finite number of eigenvectors and uses results from finite rank operators. Their proof of consistency requires stronger and more complicated conditions than ours. The consistency for functional SIR with ridge regularization is proven in [19], but it is of a weaker form than our result. We remark that the consistency results for functional SIR can be improved using a similar argument in this paper.In the following, we state the consistency results for the Tikhonov regularization. A similar result can be proved for the ridge regularization while the details are omitted.Theorem 9.
Assume𝔼k(X,X)2<∞, limn→∞s(n)=0 and  limn→∞sn=∞; then
(31)|〈β^j,ϕ(·)〉-〈βj,ϕ(·)〉|=op(1),j=1,…,dΓ,
where dΓ is the rank of Γ, 〈βj,ϕ(·)〉 is the projection onto the jth e.d.r., and 〈β^j,ϕ(·)〉 is the projection onto the jth e.d.r. as estimated by regularized KSIR.
If the e.d.r. directions{βj}j=1dΓ depend only on a finite number of eigenvectors of the covariance operator Σ, the rate of convergence is O(n-1/4).This theorem is a direct corollary of the following theorem which is proven in AppendixC.Theorem 10.
Define the projection operator and its complement for eachN≥1(32)ΠN=∑i=1N‍ei⊗ei,ΠN⊥=I-ΠN=∑i=N+1∞‍ei⊗ei,
where {ei}i=1∞ are the eigenvectors of the covariance operator Σ as defined in (8), with the corresponding eigenvalues denoted by wi.
Assume𝔼k(X,X)2<∞. For each N≥1, the following holds:
(33)∥(Σ^2+sI)-1Σ^Γ^-T∥HS=Op(1sn)+∑j=1dΓ‍(swN2∥ΠN(u~j)∥+∥ΠN⊥(u~j)∥),
where u~j=Σ-1uj and {uj}i=1dΓ are the eigenvectors of Γ as defined in (14) and ∥·∥HS denotes the Hilbert-Schmidt norm of a linear operator.
Ifs=s(n) satisfy s→0 and sn→∞ as n→∞, then
(34)∥(Σ^2+sI)-1Σ^Γ^-T∥HS=op(1).


## 4. Application to Simulated and Real Data

In this section, we compare regularized kernel sliced inverse regression (RKSIR) with several other SIR-related dimension reduction methods. The comparisons are used to address two questions:(1) does regularization improve the performance of kernel sliced inverse regression, and (2) does the nonlinearity of kernel sliced inverse regression improve the prediction accuracy?We would like to remark that the assessment of nonlinear dimension reduction methods could be more difficult than that of linear ones. When the feature mappingϕ for an RKHS is not available, we do not know the true e.d.r. directions or subspace. So in that case, we will use the prediction accuracy to evaluate the goodness of RKSIR.


### 4.1. Importance of Nonlinearity and Regularization

Our first example illustrates that both the nonlinearity and regularization of RKSIR can significantly improve prediction accuracy.The regression model has ten predictor variablesX=(X1,…,X10) with each one following a normal distribution Xi~N(0,1). A univariate response is given as
(35)Y=(sin(X1)+sin(X2))(1+sin(X3))+ε
with noise ε~N(0,0.12). We compare the effectiveness of the linear dimension reduction methods SIR, SAVE, and PHD with RKSIR, by examining the predictive accuracy of a nonlinear kernel regression model on the reduced space. We generate 100 training samples and apply the above methods to compute the e.d.r. directions. In RKSIR, we used an additive Gaussian kernel as follows:
(36)k(x,z)=∑j=1d‍exp(-(xj-zj)22σ2),
where σ=2. After projecting the training samples on the estimated e.d.r. directions, we train a Gaussian kernel regression model based on the new variates. Then, the mean square error is computed on 1000 independent test samples. This experiment is repeated 100 times with all parameters set by cross-validation. The results are summarized in Figure 1.Dimension reduction for model (35): (a) mean square error for RKSIR with different regularization parameters; (b) mean square error for various dimension reduction methods. The blue line represents the mean square error without using dimension reduction. (c)-(d) RKSIR variates versus nonlinear factors.
(a)
(b)
(c)
(d)Figure1(a) displays the accuracy for RKSIR as a function of the regularization parameter, illustrating the importance of selecting regularization parameters. KSIR without regularization performs much worse than the RKSIR. Figure 1(b) displays the prediction accuracy of various dimension reduction methods.RKSIR outperforms all the linear dimension reduction methods, which illustrates the power of nonlinearity introduced in RKSIR. It also suggests that there are essentially two nonlinear e.d.r. directions. This observation seems to agree with the model in (35). Indeed, Figures 1(c) and 1(d) show that the first two e.d.r. directions from RKSIR estimate the two nonlinear factors well.


### 4.2. Effect of Regularization

This example illustrates the effect of regularization on the performance of KSIR as a function of the anisotropy of the predictors.The regression model has ten predictor variablesX=(X1,…,X10) and a univariate response specified by
(37)Y=X1+X22+ε,ε~N(0,0.12),
where X~N(0,ΣX) and ΣX=QΔQ with Q a randomly chosen orthogonal matrix and Δ=diag(1θ,2θ,…,10θ). We will see that increasing the parameter θ∈[0,∞) increases the anisotropy of the data that increases the difficulty of identifying the correct e.d.r. directions.For this model, it is known that SIR will miss the direction along the second variableX2. So we focus on the comparison of KSIR and RKSIR in this example.If we use a second-order polynomial kernelk(x,z)=(1+xTz)2 that corresponds to the feature space
(38)Φ(X)={1,Xi,(XiXj)i≤j},i,j=1,…,10,
then X1+X22 can be captured in one e.d.r. direction. Ideally the first KSIR variate v=〈β1,ϕ(X)〉 should be equivalent to X1+X22 modulo shift and scale
(39)v-𝔼v∝X1+X22-𝔼(X1+X22).
So for this example given estimates of KSIR variates at the n data points {v^i}i=1n={〈β^1,ϕ(xi)〉}i=1n, the error of the first e.d.r. direction can be measured by the least squares fitting of v with respect to (X1+X22)(40)error=mina,b∈ℝ1n∑i=1n(v^i-(a(xi,1+xi,22)+b))2.We drew200 observations from the model specified in (37), and then we applied the two dimension reduction methods KSIR and RKSIR. The mean and standard errors of 100 repetitions of this procedure are reported in Figure 2. The result shows that KSIR becomes more and more unstable as θ increases and the regularization helps to reduce this instability.Figure 2
Error in e.d.r. as a function ofθ.


### 4.3. Importance of Nonlinearity and Regularization in Real Data

When SIR is applied to classification problems, it is equivalent to a Fisher discriminant analysis. For the case of multiclass classification, it is natural to use SIR and consider each class as a slice. Kernel forms of Fisher discriminant analysis (KFDA)  [8] have been used to construct nonlinear discriminant surfaces and the regularization has improved performance of KFDA  [25]. In this example, we show that this idea of adding a nonlinearity and a regularization term improves predictive accuracy in a real multiclass classification data set, the classification of handwritten digits.The MNIST data set (Y. LeCun,http://yann.lecun.com/exdb/mnist/) contains 60,000 images of handwritten digits {0,1,2,…,9} as training data and 10,000 images as test data. Each image consists of p=28×28=784 gray-scale pixel intensities. It is commonly believed that there is clear nonlinear structure in this 784-dimensional space.We compared regularized SIR (RSIR) as in (26), KSIR, and RKSIR, on this data to examine the effect of regularization and nonlinearity. Each draw of the training set consisted of 100 observations of each digit. We then computed the top 10 e.d.r. directions using these 1000 observations and 10 slices, one for each digit. We projected the 10,000 test observations onto the e.d.r. directions and used a k-nearest neighbor (kNN) classifier with k=5 to classify the test data. The accuracy of the kNN classifier without dimension reduction was used as a baseline. For KSIR and RKSIR we used a Gaussian kernel with the bandwith parameter set as the median pairwise distance between observations. The regularization parameter was set by cross-validation.The mean and standard deviation of the classification accuracy over100 iterations of this procedure are reported in Table 1. The first interesting observation is that the linear dimension reduction does not capture discriminative information, as the classification accuracy without dimension reduction is better. Nonlinearity does increase classification accuracy and coupling regularization with nonlinearity increases accuracy more. This improvement is dramatic for 2, 3,5, and 8.Table 1
Mean and standard deviations for error rates in classification of digits.
Digit
RKSIR
KSIR
RSIR
kNN
0
0.0273 (0.0089)
0.0472 (0.0191)
0.0487 (0.0128)
0.0291 (0.0071)
1
0.0150 (0.0049)
0.0177 (0.0051)
0.0292 (0.0113)
0.0052 (0.0012)
2
0.1039 (0.0207)
0.1475 (0.0497)
0.1921 (0.0238)
0.2008 (0.0186)
3
0.0845 (0.0208)
0.1279 (0.0494)
0.1723 (0.0283)
0.1092 (0.0130)
4
0.0784 (0.0240)
0.1044 (0.0461)
0.1327 (0.0327)
0.1617 (0.0213)
5
0.0877 (0.0209)
0.1327 (0.0540)
0.2146 (0.0294)
0.1419 (0.0193)
6
0.0472 (0.0108)
0.0804 (0.0383)
0.0816 (0.0172)
0.0446 (0.0081)
7
0.0887 (0.0169)
0.1119 (0.0357)
0.1354 (0.0172)
0.1140 (0.0125)
8
0.0981 (0.0259)
0.1490 (0.0699)
0.1981 (0.0286)
0.1140 (0.0156)
9
0.0774 (0.0251)
0.1095 (0.0398)
0.1533 (0.0212)
0.2006 (0.0153)
Average
0.0708 (0.0105)
0.1016 (0.0190)
0.1358 (0.0093)
0.1177 (0.0039)


## 4.1. Importance of Nonlinearity and Regularization

Our first example illustrates that both the nonlinearity and regularization of RKSIR can significantly improve prediction accuracy.The regression model has ten predictor variablesX=(X1,…,X10) with each one following a normal distribution Xi~N(0,1). A univariate response is given as
(35)Y=(sin(X1)+sin(X2))(1+sin(X3))+ε
with noise ε~N(0,0.12). We compare the effectiveness of the linear dimension reduction methods SIR, SAVE, and PHD with RKSIR, by examining the predictive accuracy of a nonlinear kernel regression model on the reduced space. We generate 100 training samples and apply the above methods to compute the e.d.r. directions. In RKSIR, we used an additive Gaussian kernel as follows:
(36)k(x,z)=∑j=1d‍exp(-(xj-zj)22σ2),
where σ=2. After projecting the training samples on the estimated e.d.r. directions, we train a Gaussian kernel regression model based on the new variates. Then, the mean square error is computed on 1000 independent test samples. This experiment is repeated 100 times with all parameters set by cross-validation. The results are summarized in Figure 1.Dimension reduction for model (35): (a) mean square error for RKSIR with different regularization parameters; (b) mean square error for various dimension reduction methods. The blue line represents the mean square error without using dimension reduction. (c)-(d) RKSIR variates versus nonlinear factors.
(a)
(b)
(c)
(d)Figure1(a) displays the accuracy for RKSIR as a function of the regularization parameter, illustrating the importance of selecting regularization parameters. KSIR without regularization performs much worse than the RKSIR. Figure 1(b) displays the prediction accuracy of various dimension reduction methods.RKSIR outperforms all the linear dimension reduction methods, which illustrates the power of nonlinearity introduced in RKSIR. It also suggests that there are essentially two nonlinear e.d.r. directions. This observation seems to agree with the model in (35). Indeed, Figures 1(c) and 1(d) show that the first two e.d.r. directions from RKSIR estimate the two nonlinear factors well.


## 4.2. Effect of Regularization

This example illustrates the effect of regularization on the performance of KSIR as a function of the anisotropy of the predictors.The regression model has ten predictor variablesX=(X1,…,X10) and a univariate response specified by
(37)Y=X1+X22+ε,ε~N(0,0.12),
where X~N(0,ΣX) and ΣX=QΔQ with Q a randomly chosen orthogonal matrix and Δ=diag(1θ,2θ,…,10θ). We will see that increasing the parameter θ∈[0,∞) increases the anisotropy of the data that increases the difficulty of identifying the correct e.d.r. directions.For this model, it is known that SIR will miss the direction along the second variableX2. So we focus on the comparison of KSIR and RKSIR in this example.If we use a second-order polynomial kernelk(x,z)=(1+xTz)2 that corresponds to the feature space
(38)Φ(X)={1,Xi,(XiXj)i≤j},i,j=1,…,10,
then X1+X22 can be captured in one e.d.r. direction. Ideally the first KSIR variate v=〈β1,ϕ(X)〉 should be equivalent to X1+X22 modulo shift and scale
(39)v-𝔼v∝X1+X22-𝔼(X1+X22).
So for this example given estimates of KSIR variates at the n data points {v^i}i=1n={〈β^1,ϕ(xi)〉}i=1n, the error of the first e.d.r. direction can be measured by the least squares fitting of v with respect to (X1+X22)(40)error=mina,b∈ℝ1n∑i=1n(v^i-(a(xi,1+xi,22)+b))2.We drew200 observations from the model specified in (37), and then we applied the two dimension reduction methods KSIR and RKSIR. The mean and standard errors of 100 repetitions of this procedure are reported in Figure 2. The result shows that KSIR becomes more and more unstable as θ increases and the regularization helps to reduce this instability.Figure 2
Error in e.d.r. as a function ofθ.


## 4.3. Importance of Nonlinearity and Regularization in Real Data

When SIR is applied to classification problems, it is equivalent to a Fisher discriminant analysis. For the case of multiclass classification, it is natural to use SIR and consider each class as a slice. Kernel forms of Fisher discriminant analysis (KFDA)  [8] have been used to construct nonlinear discriminant surfaces and the regularization has improved performance of KFDA  [25]. In this example, we show that this idea of adding a nonlinearity and a regularization term improves predictive accuracy in a real multiclass classification data set, the classification of handwritten digits.The MNIST data set (Y. LeCun,http://yann.lecun.com/exdb/mnist/) contains 60,000 images of handwritten digits {0,1,2,…,9} as training data and 10,000 images as test data. Each image consists of p=28×28=784 gray-scale pixel intensities. It is commonly believed that there is clear nonlinear structure in this 784-dimensional space.We compared regularized SIR (RSIR) as in (26), KSIR, and RKSIR, on this data to examine the effect of regularization and nonlinearity. Each draw of the training set consisted of 100 observations of each digit. We then computed the top 10 e.d.r. directions using these 1000 observations and 10 slices, one for each digit. We projected the 10,000 test observations onto the e.d.r. directions and used a k-nearest neighbor (kNN) classifier with k=5 to classify the test data. The accuracy of the kNN classifier without dimension reduction was used as a baseline. For KSIR and RKSIR we used a Gaussian kernel with the bandwith parameter set as the median pairwise distance between observations. The regularization parameter was set by cross-validation.The mean and standard deviation of the classification accuracy over100 iterations of this procedure are reported in Table 1. The first interesting observation is that the linear dimension reduction does not capture discriminative information, as the classification accuracy without dimension reduction is better. Nonlinearity does increase classification accuracy and coupling regularization with nonlinearity increases accuracy more. This improvement is dramatic for 2, 3,5, and 8.Table 1
Mean and standard deviations for error rates in classification of digits.
Digit
RKSIR
KSIR
RSIR
kNN
0
0.0273 (0.0089)
0.0472 (0.0191)
0.0487 (0.0128)
0.0291 (0.0071)
1
0.0150 (0.0049)
0.0177 (0.0051)
0.0292 (0.0113)
0.0052 (0.0012)
2
0.1039 (0.0207)
0.1475 (0.0497)
0.1921 (0.0238)
0.2008 (0.0186)
3
0.0845 (0.0208)
0.1279 (0.0494)
0.1723 (0.0283)
0.1092 (0.0130)
4
0.0784 (0.0240)
0.1044 (0.0461)
0.1327 (0.0327)
0.1617 (0.0213)
5
0.0877 (0.0209)
0.1327 (0.0540)
0.2146 (0.0294)
0.1419 (0.0193)
6
0.0472 (0.0108)
0.0804 (0.0383)
0.0816 (0.0172)
0.0446 (0.0081)
7
0.0887 (0.0169)
0.1119 (0.0357)
0.1354 (0.0172)
0.1140 (0.0125)
8
0.0981 (0.0259)
0.1490 (0.0699)
0.1981 (0.0286)
0.1140 (0.0156)
9
0.0774 (0.0251)
0.1095 (0.0398)
0.1533 (0.0212)
0.2006 (0.0153)
Average
0.0708 (0.0105)
0.1016 (0.0190)
0.1358 (0.0093)
0.1177 (0.0039)


## 5. Discussion

The interest in manifold learning and nonlinear dimension reduction in both statistics and machine learning has led to a variety of statistical models and algorithms. However, most of these methods are developed in the unsupervised learning framework. Therefore, the estimated dimensions may not be optimal for the regression models. Our work incorporates nonlinearity and regularization to inverse regression approaches and results in a robust response driven nonlinear dimension reduction method.RKHS has also been introduced into supervised dimension reduction in [26], where the conditional covariance operators on such kernel spaces are used to characterize the conditional independence between linear projections and the response variable. Therefore, their method estimates linear e.d.r. subspaces, while in [10, 11] and in this paper, RKHS is used to model the e.d.r. subspace, which leads to nonlinear dimension reduction.There are several open issues in regularized kernel SIR method, such as the selection of kernels, regularization parameters, and number of dimensions. A direct assessment of the nonlinear e.d.r directions is expected to reduce the computational burden in procedures based on cross validation. While these are well established in linear dimension reduction, however, little is known for nonlinear dimension reduction. We would like to leave them for future research.There are some interesting connections between KSIR and functional SIR, which are developed by Ferré and his coauthors in a series of papers  [14, 17, 19]. In functional SIR, the observable data are functions and the goal is to find linear e.d.r. directions for functional data analysis. In KSIR, the observable data are typically not functions but mapped into a function space in order to characterize nonlinear structures. This suggests that computations involved in functional SIR can be simplified by a parametrization with respect to an RKHS or using a linear kernel in the parametrized function space. On the other hand, from a theoretical point of view, KSIR can be viewed as a special case of functional SIR, although our current theoretical results on KSIR are different from the ones for functional SIR.



---

*Source: 540725-2013-07-17.xml*