Year: 2018

Abstract:
Online learning has attracted great attention due to the increasing demand for systems that have the ability of learning and evolving. When the data to be processed is also high dimensional and dimension reduction is necessary for visualization or prediction enhancement, online dimension reduction will play an essential role. The purpose of this paper is to propose new online learning approaches for supervised dimension reduction. Our first algorithm is motivated by adapting the sliced inverse regression (SIR), a pioneer and effective algorithm for supervised dimension reduction, and making it implementable in an incremental manner. The new algorithm, called incremental sliced inverse regression (ISIR), is able to update the subspace of significant factors with intrinsic lower dimensionality fast and efficiently when new observations come in. We also refine the algorithm by using an overlapping technique and develop an incremental overlapping sliced inverse regression (IOSIR) algorithm. We verify the effectiveness and efficiency of both algorithms by simulations and real data applications.

heading: 1 INTRODUCTION
text: Dimension reduction aims to explore low dimensional representation for high dimensional data. It helps to promote our understanding of the data structure through visualization and enhance the predictive performance of machine learning algorithms by preventing the “curse of dimensionality”. Therefore, as high dimensional data become ubiquitous in modern sciences, dimension reduction methods are playing more and more important roles in data analysis. Dimension reduction algorithms can be either unsupervised or supervised. Principle component analysis (PCA) might be the most popular unsupervised dimension reduction method. Other unsupervised dimension reduction methods include the kernel PCA, multidimensional scaling, and manifold learning based methods such as isometric mapping and local linear embedding. Unlike unsupervised dimension reduction, supervised dimension reduction involves a response variable. It finds the intrinsic lower-dimensional representations that are relevant to the prediction of the response values. Supervised dimension reduction methods can date back to the well known linear discriminant analysis (LDA) while its blossom occurred in the last twenty years. Many approaches have been proposed and successfully applied in various scientific domains; see Li (1991); Cook & Weisberg (1991); Li (1992); Xia et al. (2002); Setodji & Cook (2004); Fukumizu et al. (2004); Li & Yin (2008); Wu (2008); Wu et al. (2009); Cook & Zhang (2014) and the references therein.
We are in a big data era and facing the challenges of big data processing, thanks to the fast development of modern information technology. Among others, two primary challenges are the big volume and fast velocity of the data. When a data set is too big to be stored in a single machine or when the data arrives in real time and information update is needed frequently, analysis of the data in an online manner is necessary and efficient. If the data is simultaneously big and high dimensional, it becomes necessary to develop online learning approaches for dimension reduction. As PCA and LDA are the most wildly used dimension reduction techniques, a bunch of PCA-based and LDA-based online dimension reduction algorithms has been proposed. Incremental PCA have been described in Hall et al. (1998; 2000); Weng et al. (2003); Zhao et al. (2006); Rodriguez & Wohlberg (2014). Incremental LDA have been developed in Pang et al. (2005); Zhao & Yuen (2008); Kim et al. (2007); Chu et al. (2015). Other strategies like QR decomposition or SVD have also been used in Chandrasekaran et al. (1997); Ye et al. (2005); Ren & Dai (2010).
In this paper, our purpose is to propose a new online learning approach for supervised dimension reduction. Our motivation is to implement the sliced inverse regression (SIR) in an incremental manner. SIR was proposed in Li (1991) and has become one of the most efficient supervised dimension reduction method. SIR and its refined versions have been found successful in many scientific areas such as bioinformatics, hyperspectral image analysis, and physics; see Cook (1994); Belhumeur et al. (1997); Becker & Fried (2003); Elnitski et al. (2003); Gannoun et al. (2004); He et al. (2003); Antoniadis et al. (2003); Li & Li (2004); Dai et al. (2006); Zhang et al. (2016). SIR can be implemented by solving an generalized eigen-decomposition problem, Γβ = λΣβ, where Γ is a matrix depending on the response variable (whose definition is described in the next section) and Σ is the covariance matrix. To make it implementable in an online manner we rewrite it as standard eigendecomposition problem Σ− 1 2 ΓΣ− 1 2η = λη where η = Σ 1 2β and adopt the ideas from incremental PCA. We need to overcome two main challenges in this process. First, how do we transform the data so that they are appropriate for the transformed PCA problem? Note that simply normalizing the data does not work. Second, online update of Σ− 1 2 , if not impossible, seems very difficult. The first contribution of this paper is to overcome these difficulties and design a workable incremental SIR method. Our second contribution will be to refine the method by an overlapping technique and design an incremental overlapping SIR algorithm.
The rest of this paper is arranged as follows. We review SIR algorithm in Section 2 and the incremental PCA algorithm in Section 3. We propose the incremental SIR algorithm in Section 4 and refine it in Section 5. Simulations are done in Section 6. We close with discussions in Section 7.

heading: 2 SLICED INVERSE REGRESSION
text: The goal of supervised dimension reduction is to find an intrinsic lower-dimensional subspace that contains all the information to predict the response variable. Assume a multivariate predictor x = (x1, x2, . . . , xp) > ∈ Rp and a scalar response y are linked by a semi-parametric regression model
y = f(β>1x,β > 2x, . . . ,β > Kx, ), (1)
where βk ∈ Rp is a p× 1 vector and is the error term independent of x. It implies
y |= x|B>x (2)
where |= denotes “statistical independence” and B = (β1,β2, . . . ,βK) is a p × K matrix. The column space of B is called the effective dimension reduction (EDR) space and each βi is an EDR direction. Note B>x contains all information for the prediction of y. The purpose of supervised dimension reduction is to learn the EDR directions from data.
Unlike the classical regression problem which regresses y against x, sliced inverse regression considers regressing x against y. With the semi-parametric model (1) and the assumption that x follows an elliptical contour distribution (e.g., normal distribution), it was proved in Li (1991) that the centered regression curve E[x|y] − E[x] falls into the K-dimensional subspace spanned by EDR directions Σβk, k = 1, . . . ,K. Consequently, all or part of the EDR directions can be recovered by solving a generalized eigenvalue decomposition problem:
Γβ = λΣβ (3)
where
Γ = E [( E[x|y]−E[x] )( E[x|y]−E[x] )>] is the covariance matrix of inverse regression curve E[x|y], Σ is the covariance matrix of x. Each eigenvector associated with a non-zero eigenvalue is an EDR direction.
Given i.i.d observations {(xi, yi)}ni=1, SIR algorithm can be implemented as follows:
1) Compute sample mean x̄ = 1n ∑n i=1 xi and the sample covariance matrix
Σ̂ = 1
n n∑ i=1 (xi − x̄)(xi − x̄)>.
2) Bin the observations into H slices according to y values. For each slice sh, h = 1, . . . ,H, compute the sample probability p̂h = nhn and the sample slice mean m̂h = 1 nh ∑ yi∈sh xi.
The matrix Γ is estimated by
Γ̂ = H∑ h=1 p̂h (m̂h − x̄) (m̂h − x̄)>.
3) Solve the generalized eigen-decomposition problem
Γ̂B̂ = Σ̂B̂Λ.
The EDR directions are estimated by the top K eigenvectors β̂k, k = 1, 2, . . . ,K.
This algorithm is not very sensitive to the choice of parameter H provided it is sufficiently larger than K while not greater than n2 . Root-n consistency is usually promised. It is suggested samples are evenly distributed into the H slices for the best performance.

heading: 3 INCREMENTAL PCA
text: PCA looks for directions along which the data have the largest variances. It is implemented by solving an eigen-decomposition problem
Σ̂u = λu. (4)
The principal components are the eigenvectors corresponding to largest eigenvalues. Throughout this paper, we assume all eigenvalues are arranged in a descending order, i.e., λ1 ≥ λ2 ≥ . . . ≥ λp. Suppose that we need to retain the top K principal components. Denote UK = [u1, . . . ,uK ], ΛK = diag(λ1, . . . , λK). We have a reduced system Σ̂UK = UKΛK .
In incremental PCA, after receiving a new coming observation x0, we need to update the reduced eigen-system to a new one
Σ̂′U′K ≈ U′KΛ′K . (5) A “=” is generally impossible unless λK+1 = . . . = λp = 0 in (4).
The idea of updating the system in Hall et al. (1998) is as follows. Compute a residual vector
v = (x0 − x̄′)−UKU>K(x0 − x̄′)
where x̄′ is the mean of all observations (including x0). It defines the component of x0 that is perpendicular with the subspace defined by UK . If x0 lies exactly within the current eigenspace, then the residual vector is zero and there is no need to update the system. Otherwise, we normalize v to obtain v̄ = v‖v‖ .We may reasonably assume each column vector of U ′ K is a linear combination of column vectors of UK and v̄. (Note this is exactly true if λK+1 = . . . = λp = 0.) This allows us to write
[U′K ,u ′ K+1] = [UK , v̄]R.
where R is a (K + 1) × (K + 1) rotation matrix and u′K+1 is an approximation of the (K + 1)th eigenvector of Σ̂′. So we have
Σ̂′[UK , v̄]R = [UK , v̄]RΛ ′ K+1,
which is equivalent to
[UK , v̄] >Σ̂′[UK , v̂]R = RΛ ′ K+1.
This is an eigen-decomposition problem of dimensionality K + 1 p. It solves the rotation matrix R and allows us to update principal components to U′K , given by the first K column of [UK , v̄]R. If we need to increase the number of principal components, we can just update the system to K ′ = K+ 1 and U′K′ = [U ′ K ,u ′ K+1]. This incremental PCA algorithm was shown convergent to a stable solution when the sample size increases (Hall et al., 1998).

heading: 4 INCREMENTAL SIR
text: Our idea to develop the incremental sliced inverse regression (ISIR) is motivated by reformulating SIR problem to a PCA problem. To this end, we define η = Σ 1 2β, called the standardized EDR direction, and rewrite the generalized eigen-decomposition problem (3) as an eigen-decomposition problem
Σ− 1 2 ΓΣ− 1 2η = λη. (6)
Note that Σ− 1 2 ΓΣ− 1 2 is the covariance matrix of Σ− 1 2 E[x|y]. So (6) can be regarded as a PCA problem with data collected for Σ− 1 2 E[x|y]. To apply the ideas from IPCA to this transformed PCA problem, however, is not as direct as it looks like. We face two main challenges. First, when a new observation (x0, y0) is received, we need to transform it to an observation for the standardized inverse regression curve. This is different from simply standardizing the data. Second, conceptually, we need to update Σ− 1 2 in an online manner in order to standardize the data. This does not seem feasible. In the following, we will describe in detail how we address these challenges and make the ISIR implementable.
Suppose we have n observations in hand with well defined sample slice probabilities p̂h and means (m̂h, ȳh) for h = 1, . . . ,H, and the eigenvectors B̂ = [β̂1, . . . , β̂K ] of the generalized eigendecomposition problem Γ̂ = λΣ̂β. Then with ΛK = diag(λ1, . . . , λK), we have
Γ̂B̂ = Σ̂B̂ΛK .
Denote Ξ = Σ̂ 1 2 B̂. We have
Σ̂ − 12 Γ̂Σ̂ − 12 Ξ = ΞΛK . (7)
When we have a new observation (x0, y0), we first locate which slice it belongs to according to the distances from y0 to sample slice mean values ȳh of the response variable. Let us suppose the distance from y0 to ȳk is the smallest. So we place the new observation into the slice k and update sample slice probabilities by p̂′h = np̂h n+1 for h 6= k and p ′ k = np′k+1 n+1 . Let nk = npk be the number of observations in slice k before receiving the new observation. For slice mean values we update
m̂′k = nk
nk + 1 m̂k +
1
nk + 1 x0
for slice k only. We can regard z0 = Σ̂ − 12 m′k as a new observation for the standardized inverse regression curve Σ− 1 2 E[x|y]. Following the idea of IPCA, we define a residual vector
v = ( z0 − Σ̂ − 12 x̄′ ) −ΞΞ> ( z0 − Σ̂ − 12 x̄′ )
and normalize it to v̄ = v‖v‖ when v is not zero. To update the eigen-decomposition system to
Σ̂′ − 12 Γ̂′ Σ̂′ − 12
Ξ′ = Ξ′Λ′K , (8)
we assume [Ξ′,η′K+1] = [Ξ, v̄]R with R being a (K+ 1)× (K+ 1) rotation matrix and η′K+1 the (K + 1)th eigenvector of Σ̂′ − 12 Γ̂′ Σ̂′ − 12 . So we have
Σ̂′ − 12 Γ̂′ Σ̂′ − 12
[Ξ, v̄]R = [Ξ, v̄]RΛ′K+1
where Λ′K+1 = diag(Λ ′ K , λ ′ K+1) and λ ′ K+1 is the (K + 1)th eigenvalue. Multiplying both sides by [Ξ, v̄]>, we obtain ( Σ̂′ − 12 [Ξ, v̄] )> Γ̂′ ( Σ̂′ − 12 [Ξ, v̄] ) R = RΛ′K+1. (9)
Note, however, since Σ̂′ − 12
cannot be easily updated, we have to avoid using it. To overcome this challenge, we notice that
Σ̂′ = n
n+ 1 Σ̂′ +
n
(n+ 1)2 (x0 − x̄)(x0 − x̄)> (10)
and the well known Sherman-Morisson formula allows us to update the inverse covariance matrix
Σ̂′ −1 = n+ 1 n Σ̂ −1 −
1 nΣ̂ −1 (x0 − x̄)(x0 − x̄)>Σ̂ −1
1 + 1n+1 (x0 − x̄)>Σ̂ −1 (x0 − x̄) . (11)
If we store Σ̂ −1 and update it incrementally, we can approximate the quantities in (9) as follows:
Σ̂′ − 12 Ξ ≈ Σ̂ − 12 Ξ = B̂,
Σ̂′ − 12 v ≈ Σ̂ −1
(m̂′k − x̄′)− B̂B̂> (m̂′k − x̄′) ,
‖v‖2 = (m̂′k − x̄′) > ( Σ̂ −1 − B̂B̂> ) (m̂′k − x̄′) ,
ṽ = Σ̂′ − 12 v̄ = Σ̂′ − 12 v
‖v‖
≈
( Σ̂ −1 − B̂B̂> ) (m̂′k − x̄′)√
(m̂′k − x̄′) > ( Σ̂ −1 − B̂B̂> ) (m̂′k − x̄′) .
So the problem (9) is approximated by[ B̂, ṽ ]> Γ̂′ [ B̂, ṽ ] R = RΛ′K+1. (12)
Finally notice that the new EDR space B̂′ = Σ̂′ − 12 Ξ′ is the first K columns of Σ̂′ − 12
[Ξ′,η′K+1] =
Σ̂′ − 12 [Ξ, v̄]R and can be approximated by the first K columns of [B̂, ṽ]R.
Note that we avoided updating the inverse square root of the covariance matrix by using the approx-
imation Σ̂ − 12 ≈ Σ̂′ − 12 . This approximation can be very accurate when n is large enough because both converge to Σ− 1 2 . Therefore, we may expect the convergence of ISIR as a corollary of the convergence of IPCA. However, when n is small, the approximation may be less accurate and result in larger difference between EDR spaces estimated by ISIR and SIR. So we recommend that ISIR be used with a warm start, that is, using SIR first on a small amount of data before using ISIR.
In terms of memory, the primary requirement is the storage of Σ̂ −1
, the slice mean matrix M̂ = [m̂1, . . . , m̂H ], and the EDR space B̂. So the memory requirement is O(p2 + pH + pK). As for the computational complexity, notice that the update of Σ̂′ −1
in (11) requires the calculation of Σ̂ −1
(x0 − x̄) and matrix addition and thus has a complexity of O(p2). Since we need to store M̂ and update it sequentially, it is not efficient to store and update Γ̂′ for either memory or computation consideration. Instead, we use the fact Γ̂′ = M̂′P̂′M̂′> where P̂′ = diag(p̂′1, . . . , p̂ ′ H) and write
[B̂, ṽ]>Γ̂′ [B̂, ṽ] = [ B̂>M̂′, ṽ>M̂′ ] P̂′ [ B̂>M̂′, ṽ>M̂′ ]> .
Notice that
ṽ>M̂′ =
( Σ̂ −1 (m̂′k − x̄′) )> M̂′ − ( B̂>(m̂′k − x̄′) )>( B̂>M̂′ ) √
(m̂′k − x̄′)> ( Σ̂ −1 (m̂′k − x̄′) ) − ( B̂>(m̂′k − x̄′) )>( B̂>(m̂′k − x̄′) ) and B̂>(m̂′k− x̄′) is just the kth column of B̂>M̂′. The most time consuming computations for the matrix [B̂, ṽ]>Γ̂′ [B̂, ṽ] are B̂>M̂′ and Σ̂ −1 (m̂′k− x̄′) which have a complexity of O(pKH+ p2).
The complexity of the eigen-decomposition in (12) is O((K+ 1)3) and it requires O(p(K+ 1)2) to update B̂ to B̂′. So the computational complexity for each ISIR update isO(p2+pKH+pK2+K3). For a high dimensional problem, this is much smaller than the complexity O(p3 + p2n) of SIR.

heading: 5 REFINEMENT BY OVERLAPPING
text: In Zhang et al. (2018) an overlapping technique was introduced to SIR algorithm and shown effectively improving the accuracy of EDR space estimation. It is motivated by placing each observation in two or more adjacent slices to reduce the deviations of the sample slice means m̂h from the EDR subspace. This is equivalent to using each observation two or more times. In this section, we adopt the overlapping technique to ISIR algorithm above to develop an incremental overlapping sliced inverse regression (IOSIR) algorithm and wish it refines ISIR.
To apply the overlapping idea, we use each observation twice. So when we have n observations, we duplicate them and assume we have N = 2n observations. When a new observation (x0, y0) is received, we duplicate it and assume we receive two identical observations. Based on the y0 value we place the first copy into the slice sk if ȳk is closest to y0 and run ISIR update as described in Section 4. Note that if ȳ1 < y0 < ȳH , then y0 must fall into the interval [ȳk′ , ȳk] with k′ = k− 1 or it falls into [ȳk, ȳk′ ] with k′ = k + 1. So we place the second copy of the new observation into the slice sk′ , an adjacent slice to sk, and run ISIR algorithm again. If y0 ≤ ȳ1 or y0 > ȳH , the second copy will be still placed into sk to guarantee all observations are weighted equally.
As OSIR has superior performance over SIR, we expect IOSIR will perform better than ISIR by a price of double calculation time.
We remark that SIR and ISIR can be used for both regression problems and classification problems. But since the concept of “adjacent slice” cannot be defined for categorical values (as is the case in classification problems), IOSIR can only be used for regression problems where the response variable is numeric.

heading: 6 SIMULATIONS
text: In this section, we will verify the effectiveness of ISIR and IOSIR with simulations on artificial and real-world data. Comparisons will be made between them and SIR.

heading: 6.1 ARTIFICIAL DATA
text: In the simulations with artificial data, since we know the true model, we measure the performance by the accuracy of the estimated EDR space. We adopt the trace correlation r(K) = trace(PBPB̂)/K used in Ferré (1998) as the criterion, where PB and PB̂ are the projection operators onto the true EDR space B and the estimated EDR space B̂, respectively.
We consider the following model from Li (1991)
y = x1(x1 + x2 + 1) + , (13)
where x = [x1, x2, . . . , xp]> follows a multivariate normal distribution, follows standard normal distribution and is independent of x. It has K = 2 effective dimensions with β1 = (1, 0, 0, . . . , 0) > and β2 = (0, 1, 0 . . . , 0) >. We conduct the simulation in p = 10 dimensional space and select the number of slices as H = 10. We give the algorithm a warm start with the initial guess of the EDR space obtained by applying SIR algorithm to a small data set of 40 observations. Then a total of 400 new observations will be fed to update the EDR space one by one. SIR, ISIR, and IOSIR are applied when each observation was fed in and we calculate their trace correlation and cumulative computation time. We repeat this process 100 times. The mean trace correlations of all three methods are reported in Figure 1(a) and the mean cumulative time is in Figure 1(b). We see that ISIR performs quite similar to SIR. IOSIR slightly outperforms ISIR and SIR. ISIR is much faster than SIR. IOSIR gains higher accuracy by sacrificing on computation time. This verifies the convergence and efficiency of ISIR and IOSIR.

heading: 6.2 REAL DATA APPLICATION
text: We validate the reliability of ISIR on two real data sets: the Concrete Compressive Strength and Cpusmall (available on https://www.csie.ntu.edu.tw/˜cjlin/libsvmtools/ datasets/regression.html). There have been many proposed algorithms to increase the prediction accuracy on these data sets (Yeh, 1998; 1999; 2007; Öztaş et al., 2006; Gu et al., 2015). We do not intend to outperform those methods. Our goal is to compare the performance of supervised dimension reduction algorithms and verify the effectiveness of our incremental methods.
The Concrete Compressive Strength data has p = 8 predictors and 1030 samples. We use H = 10 and K = 3 to run SIR, ISIR, and IOSIR. We select 50 observations to warm start ISIR and IOSIR algorithms, then 700 observations are fed sequentially. The left 280 observations are used as test data. After each new observation is received we estimate the EDR space, project the available training set to the estimated EDR space, build a regression model using the k-nearest neighbor method, and compute the MSE on the test data set. This process is repeated 100 times and the average MSE was reported in Figure 2(a). For the Cpusmall data, which has p = 12 predictors and 8192 samples, we do the experiment with H = 10, K = 3, 50 observations to warm start ISIR and IOSIR, 2000 observations for sequential training, and 6142 observations for testing. The average MSE was plotted in Figure 2(b). The results indicate both ISIR and IOSIR are as effective as SIR.

heading: 7 CONCLUSIONS AND DISCUSSIONS
text: We proposed two online learning approaches for supervised dimension reduction, namely, ISIR and IOSIR. They are motivated by standardizing the data and reformulate the SIR algorithm to a PCA problem. However, data standardization is only used to motivate the algorithm while not explicitly calculated in the algorithms. We proposed to use Sherman Morrison formula to online update Σ−1 and some approximated calculations to circumvent explicit data standardization. This novel idea played a key role in our algorithm design. Both algorithms are shown effective and efficient. While IOSIR does not apply to classification problems, it is usually superior over ISIR in regression problems.
We remark that the purpose of ISIR and IOSIR is to keep the dimension reduction accuracy in the situation that a batch learning is not suitable. This is especially the case for streaming data where information update and system involving is necessary whenever new data becomes available. When the whole data set is given and one only needs the EDR space from batch learning, ISIR or IOSIR is not necessarily more efficient than SIR because their complexity to run over the whole sample path is O(p2N), comparable to the complexity O(p3 + p2N) of SIR.
There are two open problems worth further investigation. First, the need to store and use Σ−1 during the updating process is the main bottleneck for ISIR and IOSIR when the dimensionality of
the data is ultrahigh. Second, for SIR and other batch dimension reduction methods, many methods have been proposed to determine the intrinsic dimension K; see e.g. Li (1991); Schott (1994); Bura & Cook (2001); Bai & He (2004); Barrios & Velilla (2007); Nkiet (2008). They depend on all p eigenvalues of the generalized eigen-decomposition problem and are impractical for incremental learning. We do not have obvious solutions to these problems at this moment and would like to leave them for future research.

Review summaries:
This paper proposes an online learning algorithm for supervised dimension reduction, called incremental sliced inverse regression (ISIR). The key idea is converting the SIR problem into PCA problem by using the inverse of covariance matrix. After the transformation, we can use incremental PCA to compute the top eigenvector and obtain the approximate solution of SIR in streaming way. The authors also extend ISIR to overlapping case.

Sliced Inverse Regression is a well-known technique for finding EDR space in supervised dimension reduction problems, under condition that input X is normally disributed. When the number of dimensions is large and an access to observations is online, finding eigenvalues of covariance matrices (online) could become computationally costly. So, the paper focuses the problem of updating few principal components of covariance matrices as new examples come.