Title: ONLINE LEARNING FOR SUPERVISED DIMENSION RE-DUCTION

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.

1. INTRODUCTION:

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.

2. SLICED INVERSE REGRESSION:

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 = (x 1 , x 2 , . . . , x p ) ∈ R p and a scalar response y are linked by a semi-parametric regression model

EQUATION

where β k ∈ R p is a p × 1 vector and is the error term independent of x. It implies

EQUATION

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:

EQUATION

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 {(x i , y i )} n i=1 , SIR algorithm can be implemented as follows:

1) Compute sample mean x = 1 n n i=1 x i and the sample covariance matrix

Σ = 1 n n i=1 (x i -x)(x i -x) .

2) Bin the observations into H slices according to y values. For each slice s h , h = 1, . . . , H, compute the sample probability p h = n h n and the sample slice mean

m h = 1 n h yi∈s h x i . 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 n 2 . Root-n consistency is usually promised. It is suggested samples are evenly distributed into the H slices for the best performance.

3. INCREMENTAL PCA:

PCA looks for directions along which the data have the largest variances. It is implemented by solving an eigen-decomposition problem

EQUATION

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

U K = [u 1 , . . . , u K ], Λ K = diag(λ 1 , . . . , λ K ). We have a reduced system ΣU K = U K Λ K .

In incremental PCA, after receiving a new coming observation x 0 , we need to update the reduced eigen-system to a new one

EQUATION

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 = (x 0 -x ) -U K U K (x 0 -x )

where x is the mean of all observations (including x 0 ). It defines the component of x 0 that is perpendicular with the subspace defined by U K . If x 0 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 U K and v. (Note this is exactly true if λ K+1 = . . . = λ p = 0.) This allows us to write

[U K , u K+1 ] = [U K , 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

Σ [U K , v]R = [U K , v]RΛ K+1 ,

which is equivalent to

[U K , v] Σ [U K , 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 [U K , 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) .

4. INCREMENTAL SIR:

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 (x 0 , y 0 ) 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 , ȳh ) for h = 1, . . . , H, and the eigenvectors B = [ β 1 , . . . , β K ] of the generalized eigendecomposition problem Γ = λ Σβ. Then with Λ K = diag(λ 1 , . . . , λ K ), we have 

EQUATION

m k = n k n k + 1 m k + 1 n k + 1 x 0

for slice k only. We can regard z 0 = Σ -1 2 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 = z 0 -Σ -1 2 x -ΞΞ z 0 -Σ -1 2 x

and normalize it to v = v v when v is not zero. To update the eigen-decomposition system to

EQUATION

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 Σ -1 2 Γ Σ -1 2 . So we have Σ -1 2 Γ Σ -1 2 [Ξ, 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

EQUATION

Note, however, since Σ -1 2 cannot be easily updated, we have to avoid using it. To overcome this challenge, we notice that

EQUATION

and the well known Sherman-Morisson formula allows us to update the inverse covariance matrix

EQUATION

)

If we store Σ -1

and update it incrementally, we can approximate the quantities in (9) as follows:

Σ -1 2 Ξ ≈ Σ -1 2 Ξ = B, Σ -1 2 v ≈ Σ -1 ( m k -x ) -B B ( m k -x ) , v 2 = ( m k -x ) Σ -1 -B B ( m k -x ) , v = Σ -1 2 v = Σ -1 2 v v ≈ Σ -1 -B B ( m k -x ) ( m k -x ) Σ -1 -B B ( m k -x )

.

So the problem ( 9) is approximated by

EQUATION

Finally notice that the new EDR space

B = Σ -1 2 Ξ is the first K columns of Σ -1 2 [Ξ , η K+1 ] = Σ -1 2 [Ξ, 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 Σ -1 2 ≈ Σ -1

2 . 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(p 2 + pH + pK). As for the computational complexity, notice that the update of Σ -1 in (11) requires the calculation of Σ -1 (x 0 -x) and matrix addition and thus has a complexity of O(p 2 ). 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, ṽ] Γ [ B, ṽ] = B M , ṽ M P B M , ṽ 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, ṽ] Γ [ B, ṽ] are B M and Σ -1

( m k -x ) which have a complexity of O(pKH + p 2 ).

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 is O(p 2 +pKH +pK 2 +K 3 ).

For a high dimensional problem, this is much smaller than the complexity O(p 3 + p 2 n) of SIR.

5. REFINEMENT BY OVERLAPPING:

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 (x 0 , y 0 ) is received, we duplicate it and assume we receive two identical observations. Based on the y 0 value we place the first copy into the slice s k if ȳk is closest to y 0 and run ISIR update as described in Section 4. Note that if ȳ1 < y 0 < ȳH , then y 0 must fall into the interval

[ȳ k , ȳk ] with k = k -1 or it falls into [ȳ k , ȳk ] with k = k + 1.

So we place the second copy of the new observation into the slice s k , an adjacent slice to s k , and run ISIR algorithm again. If y 0 ≤ ȳ1 or y 0 > ȳH , the second copy will be still placed into s k 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.

6. SIMULATIONS:

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.

6.1. ARTIFICIAL DATA:

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(P B P B )/K used in Ferré (1998) as the criterion, where P B and P B 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)

EQUATION

where x = [x 1 , x 2 , . . . , x p ] 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. 

6.2. REAL DATA APPLICATION:

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; Ö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 

7. CONCLUSIONS AND DISCUSSIONS:

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(p 2 N ), comparable to the complexity O(p 3 + p 2 N ) 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.

:



Figures and Tables:

When we have a new observation (x 0 , y 0 ), we first locate which slice it belongs to according to the distances from y 0 to sample slice mean values ȳh of the response variable. Let us suppose the distance from y 0 to ȳk is the smallest. So we place the new observation into the slice k and update sample slice probabilities by p h = n p h n+1 for h = k and p k = np k +1 n+1 . Let n k = np k be the number of observations in slice k before receiving the new observation. For slice mean values we update

Figure 1: Simulation result for the model (13). (a) Trace correlation and (b) cumulative calculation time by SIR, ISIR, and IOSIR.

(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.

Figure 2: Mean square errors (MSE) for two real data applications: (a) for Concrete Compressive Strength data and (b) for Cpusmall data.

Peer review (iclr2019):
 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. The key idea is that in classical SIR two problems are solved sequentially (estimation of Sigma = cov(X) and finding PCA of vectors Sigma^{-0.5}X_s where s corresponds to average X of s-th slice). The second part can be reduced to an online version of PCA (Peter M Hall et al.). The paper seems consistent and experiments convince that the proposed method works. The lack of theoretical analysis is a disadvantage. Also, some typos: 1. Γb = λΣbβ (p.4 par.3) --- beta omitted 2. We can regard z0 = Σ^{-1/2}m*_k (p.4) --- bar over m omitted 3. Formula 10 looks weird (Sigma instead of Sigma prime)