# Two-Phase Incremental Kernel PCA for Learning Massive or Online Datasets

Authors: Feng Zhao; Islem Rekik; Seong-Whan Lee; Jing Liu; Junying Zhang; Dinggang Shen
Journal: Complexity
 (2019)

---

## Abstract

As a powerful nonlinear feature extractor, kernel principal component analysis (KPCA) has been widely adopted in many machine learning applications. However, KPCA is usually performed in a batch mode, leading to some potential problems when handling massive or online datasets. To overcome this drawback of KPCA, in this paper, we propose a two-phase incremental KPCA (TP-IKPCA) algorithm which can incorporate data into KPCA in an incremental fashion. In the first phase, an incremental algorithm is developed to explicitly express the data in the kernel space. In the second phase, we extend an incremental principal component analysis (IPCA) to estimate the kernel principal components. Extensive experimental results on both synthesized and real datasets showed that the proposed TP-IKPCA produces similar principal components as conventional batch-based KPCA but is computationally faster than KPCA and its several incremental variants. Therefore, our algorithm can be applied to massive or online datasets where the batch method is not available.


---

## Body

## 1. Introduction

As a conventional linear subspace analysis method, principal component analysis (PCA) can only produce linear subspace feature extractors [1], which are unsuitable for highly complex and nonlinear data distributions. In contrast, as a nonlinear extension of PCA, kernel principal component analysis (KPCA) [2] can capture the higher-order statistical information contained in data, thus producing nonlinear subspaces for better feature extraction performance. This has propelled the use of KPCA in a wide range of applications such as pattern recognition, statistical analysis, image processing, and so on [3–8]. Basically, KPCA firstly projects all samples from the input space into a kernel space using nonlinear mapping and then extracts the principal components (PCs) in the kernel space. In practice, such nonlinear mapping is performed implicitly via the “kernel trick”, where an appropriately chosen kernel function is used to evaluate the dot products of mapped samples without having to explicitly carry out the mapping. As a result, the extracted kernel principal component (KPC) of the mapped data is nonlinear with respect to the original input space.Standard KPCA has some drawbacks which limit its practical applications when handling big or online datasets.Firstly, in the training stage, KPCA needs to store and compute the eigenvectors of a N×N kernel matrix, where N is the total number of samples. This computation results in a space complexity of ON2 and a time complexity of ON3, thus rendering the evaluation of KPCA on large-scale datasets very time-consuming.Secondly, in the testing stage, the resulting kernel principal components have to be defined implicitly by the linear expression of the training data, and thus all the training data must be saved after training. For a massive dataset, this translates into high costs for storage resources and increases the computational burden during the utilization of kernel principal components (KPCs). Furthermore, KPCA is impractical for many real-world applications where online samples are progressively collected since it is used in a batch manner. This implies that each time new data arrive, KPCA has to be conducted from scratch.To overcome these limitations, many promising methods have been proposed in the past few years. These methods can be grouped into two classes. The first class is the batch-based modeling method, which requires that all training data is available for estimating KPCs. Rosipal and Girolami proposed an EM algorithm for reducing the computational cost of KPCA [9]. However, the convergence behavior of the EM algorithm to KPCA cannot be guaranteed in theory. In [6], the kernel Hebbian algorithm (KHA) was proposed as an iterative variant of KPCA algorithm. By kernelizing the generalized Hebbian algorithm (GHA), KHA computes KPCA without storing the kernel matrix, such that large-scale datasets with high dimensionality can be processed. Nonetheless, KHA has a scalar gain parameter which is either held constant or decreased according to a predetermined annealing schedule, leading to slow convergence during the training stage. To improve the convergence of KHA, gain adaptation methods were developed by providing a separate gain for each eigenvector estimate [10]. An improved version of KPCA was proposed based on the eigenvalue decomposition of a symmetric matrix [11], where datasets are divided into multiple subsets, each of which is processed separately. One of the major drawbacks of this approach is that it requires storing the kernel matrix, which means that the space complexity could be extremely large for large-scale dataset. Another variant of conventional KPCA is greedy KPCA [12, 13], which was employed to approximate the KPCs by a prior filtering of the training data. However, prior filtering of the training data could be computationally expensive. Overall, compared with standard KPCA, these batch-based modeling methods can potentially reduce the time or space complexity to some degree. Unfortunately, such methods cannot handle online data.The second class is incremental methods, which can compute KPCs incrementally to handle online data processing. Chin and Suter proposed an incremental version of KPCA [14, 15], which is called IKPCA-RS for the notational simplicity. In IKPCA-RS, singular value decomposition is used to update an eigenfeature space incrementally for incoming data. However, IKPCA-RS may lead to high time complexity especially when dealing with high-dimensional data. In [16, 17], an incremental KPCA was presented based on the empirical kernel map. It is more efficient in memory requirement than the standard KPCA. However, it is only an approximate method and only suitable for polynomial kernel function. Inspired by the incremental PCA algorithm proposed by Hall et al. [18], Kimura et al. presented an incremental KPCA algorithm [19] in which an incremental updating algorithm for eigenaxes is derived based on a set of linearly independent data. Subsequently, some modified versions are proposed by Ozawa and Takeuchi et al. [20, 21]. Furthermore, in order to incrementally deal with data streams which are given in a chunk of multiple samples at one time, other extensions of KPCA were also successively presented [22–24]. Hallgren and Northrop [25] proposed an incremental KPCA (INKPCA) by applying rank one updates to the eigendecomposition of kernel matrix. However, INKPCA needs to store the whole data when evaluated on a new sample. Notably, incremental methods have the capacity of integrating new data, initially unavailable, in some way that maintains nonincreasing memory. However, to the best of our knowledge, most of these methods operate in the kernel space where all the samples areimplicitly represented. This has two key limitations.First, a number of incremental methods may suffer from high computational cost.Second, the others can only capture the approximate KPCs rather than the accurate ones, which may affect the accuracy of its subsequent process.Before continuing, a note on mathematical notations is given as follows. We use lower case and upper case letters (e.g.,i,j,l,N) to denote scalars, lower case letters with the subscript (e.g., kij,αi) to denote an element from a matrix or a vector, lower case bold letters (e.g., x,y,α,β) to denote vectors, and upper case bold letters (e.g., A,C,M) to denote matrices. We use xT (CT) to denote the transpose of a vector (matrix) and · to denote the L2-norm of a vector. Furthermore, we adopt xii=1N to denote a set, x1,x2,…,xN to denote a matrix with N column vectors and [kij]1≤i≤M,1≤j≤N to denote a M×N matrix composed of the corresponding element kij. In this paper, xi always denotes a column vector and the inner product between xi and xj is expressed as xiTxj. The lower case bold letter φ denotes a nonlinear mapping. The mapped sample φ(x) of the input sample x is a column vector.To address these limitations, we propose a two-phase incremental KPCA (TP-IKPCA), where the mapped data is represented in anexplicit form and KPCs are updated in anexplicit space. The computational cost of the whole process is very low and the accuracy of KPCs can be theoretically guaranteed. An overview of TP-IKPCA is briefly illustrated in Figure 1. In this figure, xii=1N denotes the sample set in ad-dimensional input space and N denotes the total number of available samples. Let φ denote the nonlinear mapping which maps the sample set xii=1N into anh-dimensional implicit kernel space, resulting in the mapped sample set φxii=1N. Here,h may be very large or even infinite, depending on the specific mapping. The TP-IKPCA includes two phases.In the first phase, we develop an incremental algorithm to capture standard orthogonal basis βjj=1r of the subspace spanned by φxii=1N and thenexplicitly obtain the projection vectors yii=1N of φxii=1N by (1)yi=β1,β2,…,βrTφxi,where r denotes the number of a standard orthogonal basis βjj=1r.In the second phase, the existing incremental method of PCA is employed to capture KPCs based on the explicit data yii=1N in the projection space. In the following sections, we will detail how to incrementally express the implicit mapped data φxii=1N using an explicit form. We will also theoretically verify that performing PCA based on the implicitly mapped samples φxii=1N is equivalent to that of based on the explicit projection vectors yii=1N.Figure 1
Overview of TP-IKPCA.Here, we should clarify the relationship among some important quantities, includingd,h,r (see Figure 1). In the case of KPCA or TP-IKPCA, the sample set xii=1N in ad-dimensional input space is firstly mapped into anh-dimensional kernel space by a nonlinear mapping φ and then a linear PCA is performed based on the mapped set φxii=1N. Usually,h may be very large or even infinite, depending on the specific mapping φ. So, we usually have d≤h, which implies that the dimension of each mapped sample φxi is larger than its original dimension d. In the case of TP-IKPCA, βjj=1r denote an orthonormal basis of the subspace spanned by φxii=1N, and yii=1N is the corresponding projection vectors of the mapped set φxii=1N on the basis βjj=1r, which means that r is the dimension of the subspace spanned by φxii=1N and equals the dimension of each projected vector yii=1,2,…,N. Generally, since the mapped data φxii=1N have strong linear correlation in theh-dimensional kernel space, we have r≤h, which means that a few components generally suffice to capture the nonlinear distribution of the data. Furthermore, if the dimension r of the subspace spanned by φxii=1N is high,r may be larger than the dimension d of the input space. However, if the mapped data φxii=1N have strong linear correlation, which means r may be low and the dimension d of the input space is very high, we may get the contrast conclusion.The main contributions of our work are fourfold:(1) Presenting an algorithm to express the mapped data in an explicit form. This will help for better understanding the distribution of the mapped data in the implicit kernel space. (2) Endowing KPCA with the capacity of handling dynamic dataset. (3) Compared to the standard KPCA, the computational complexity of TP-IKPCA is reduced from ON3 to Or3 and the storage complexity from ON2 to Or2, where N denotes the number of training samples and r is the number of bases of the subspace spanned by nonlinear mapped samples. Usually the assumption that r≪N is valid, which makes TP-IKPCA very convenient for processing large-scale datasets [26]. (4) In the testing stage, the feature extraction from one sample is faster than that of the batch KPCA, since TP-IKPCA only needs to calculate the kernel functions between the new sample and r selected training samples which forms the orthonormal basis.The rest of the paper is organized as follows. Section2 briefly introduces KPCA. In Section 3, we provide a theoretical analysis of the proposed TP-IKPCA method and elucidate the concrete steps for incrementally capturing KPCs based on the projection vectors in an explicit space. The effectiveness of TP-IKPCA is demonstrated in Section 4. Finally, the conclusions of our study are given in Section 5.


## 2. Kernel Principal Component Analysis (KPCA)

In this section, we briefly outline the standard procedure of KPCA. As mentioned above, in KPCA, the input sample setxii=1N is mapped into a kernel space by a nonlinear mapping φ and then a linear PCA is performed based on φxii=1N in the kernel space.To obtain the eigenvectors in the kernel space, the covariance matrix is defined as(2)C=1N∑i=1Nφxi-c)φxi-c)T,where c=(1/N)∑i=1Nφxi. However, the eigendecomposition of C is hindered by the fact that the mapping function φ is implicit. To avoid the explicit calculation in the kernel space, a N×N kernel matrix K is defined, whose elements kij are determined by the virtue of the following kernel trick:(3)kij=φxiTφxi=kxi,xji,j=1,2,…,N.where k(·,·) is a kernel function that allows us to compute inner products in the kernel space [2].Combining (2) and (3), Schölkopf et al. [2] derived the equivalent eigenvalue problem as follows:(4)Kα=Nλα,where α=α1,α2,…,αNT denotes the column vector such that the orthogonal eigenvector ν of the covariance matrix C satisfies(5)v=∑i=1Nαiφxi.Letλ1≥λ2≥⋯≥λk>0k≤N denote the first k nonzero eigenvalues of K and α1,α2,…,αk the corresponding complete set of eigenvectors (see (4)). We can obtain the corresponding eigenvectors v1,v2,…,vk of C using (5).Consideringv1,v2,…,vk need to be normalized, α1,α2,…,αk need to satisfy(6)λlαlTαl=1l=1,2,…,k.For a test samplex, the projection of φx on thel-th nonlinear principal component can be computed by(7)dlx=vlTφx=∑i=1NαilφxiTφx=∑i=1Nαilkxi,x.where αil is theith element of αl; in other words, αl=α1l,α2l,…,αNlT.For the sake of simplicity, we assume that the mapped dataφxii=1,2,…,N is zero-centered (see (2)). The detailed description of the centering processing is given in [2].Of note, for KPCA, the kernel matrixK needs to be predefined before performing eigendecompositions. Since the size of K scales with N2, a large memory space is required for a massive dataset. Additionally, the eigendecomposition of K involves a time complexity of ON3. This can severely handicap the computation of KPCA on large datasets. In online processing applications, the arrival of a new sample requires adding a new row and a new column in K, and eigendecomposition has to be constantly reevaluated for an ever-growing kernel matrix to update the kernel subspaces. Hence, the batch KPCA is not convenient for such applications.


## 3. Explicit Representation of the Mapped Data

At present, there have been many incremental algorithms for PCA [27–34]. However, it is difficult to directly extend them to KPCA because all mapped samples φxii=1N are expressed implicitly in the kernel space. Obviously, once φxii=1N can be expressed using an explicit form, it will be straightforward to extend incremental PCA to KPCA. In fact, the geometrical structure of φxii=1N can be captured by using a standard orthogonal basis of the subspace spanned by all the samples φxii=1N [26]. Hence, we aim to explicitly express φxii=1Nusing an indirect way. This motivation comes from the following property shown in Theorem 1.Letβjj=1r denote an orthonormal basis of the subspace spanned by φxii=1N, and yii=1N is the corresponding projection vectors of φxii=1N under βjj=1r (see (1)). Theorem 1 is established as follows.Theorem 1.
If linear PCA is performed based onφxii=1N and yii=1N, respectively, then their covariance matrices have the same nonzero eigenvalues, and those corresponding eigenvectors satisfy the following relationship:(8)vl=β1,β2,…,βrul,where vl is thel-th eigenvector of the covariance matrix of φxii=1N and ul is thel-th eigenvector of the covariance matrix of yii=1N. The proof is given in Appendix A.Based on Theorem1, linear PCA based on φxii=1N can be converted into a linear PCA based on yii=1N. So, if we can write yii=1N in an explicit form, then it will be easy to further extend KPCA using existing linear incremental algorithms. To incrementally obtain the orthonormal basis βjj=1r and the projection vectors yii=1N, we firstly introduce two correlativelemmas.Lemma 2.
Letφxbjj=1r denote a basis of the subspace spanned by φxii=1N, then the orthonormal basis βjj=1r can be determined using(9)β1,β2,…,βr=φxb1,φxb2,…,φxbrD,where D = γ1/ε1,γ2/ε2,…,γr/εr. γj(j=1,2,…,r) is the eigenvector of the kernel matrix Krr = φxb1,φxb2,…,φxbrTφxb1,φxb2,…,φxbr, scilicet Krr = [k(xbs,xbt)]1≤s,t≤r. εj(j=1,2,…,r) is the corresponding eigenvalue of γj. The proof is given in Appendix B.Based on Lemma2, for any mapped sample φx∈φxii=1N, we can explicitly define its projection vector y under the orthonormal basis βjj=1r as(10)y=β1,β2,…,βrTφx=DTφxb1,φxb2,…,φxbrTφx=DTkbx,where kbx=kxb1,x,kxb2,x,…,kxbr,xT. Obviously, using the kernel function k·,· can complete the computation of y.However, the orthogonalization process using Lemma2 is a batch-based method. Subsequently, when samples are added one by one, its computational cost is still very expensive. So, inspired by the Gram-Schmidt orthogonalization process [35], we designed an online algorithm for incrementally finding the orthonormal basis and the projection vectors.Lemma 3.
Letφxbjj=1r denote a basis of the subspace spanned by φxii=1N and βjj=1r is the orthonormal basis obtained by (9). Suppose xN+1 denotes a new sample we have just included into our dataset. We derive the following properties.
(1) Ifδ = kN+1-kbN+1TDDTkbN+1 = 0, then βjj=1r is the orthonormal basis of the subspace spanned by φxii=1N+1 and the projection vector yN+1 of φxN+1 can be computed using (10). Here, kN+1 = k(xN+1,xN+1) and kbN+1 = kxb1,xN+1,kxb2,xN+1,…,kxbr,xN+1T.
(2) Ifδ = kN+1-kbN+1TDDTkbN+1≠0, then the orthonormal basis βjj=1r+1 of the subspace spanned by φxii=1N+1 can be obtained by(11)β1,β2,…,βr,βr+1=φxb1,φxb2,…,φxbr,φxbr+1D-DDTkbN+1δ01δ,where xb(r+1)=xN+1. The projection vector yN+1 of φxN+1 can be computed by(12)yN+1=kbN+1TD,δT.Obviously, based on Lemma3, it is straightforward to incrementally estimate the projection vector. Notably, the dimensionality of yii=1,2,…,N is smaller than that of yN+1 in the case of δ≠0. In fact, based on the Gram-Schmidt orthogonalization process, let yi′i=1,2,…,N denote the projection vector of φxi on the orthonormal basis βjj=1r+1, then yi′=yiT,0T. The proof of Lemma 3 is provided in Appendix C.Combining both Lemmas2 and 3, we summarize the online algorithm, which incrementally finds the orthonormal basis and the projection vectors as follows.Algorithm 4.
An online algorithm for incrementally finding the orthonormal basis and the projection vectors.
Step 0 (initialization). For the timeN=1, we found a sample x1. We suppose kx1,x1≠0. Let the set S=x1, D=1/kx1,x1, and y1=kx1,x1.
Step 1. Calculate δ for a new sample xN+1 according to (13)δ=kN+1-kbN+1TDDTkbN+1,where kN+1,D and kbN+1 have the same definition as in Lemma 3.
Step 2. If δ=0, then yN+1=DTkbN+1 and return to Step 1.
Step 3. If δ≠0, then S=S∪xN+1 and update yN+1 using (12). Finally, update D using (14) and return to Step 1. (14)D=D-DDTkbN+1δ01δ.
Obviously, if we map all the samples ofS into the kernel space and get the dataset φxbj∣xbj∈S, then the mapped samples φxbj∣xbj∈S are linearly independent. Furthermore, we can get an orthonormal basis based on φxbj∣xbj∈S and D (see (9) or (11)). In fact, taking into account the actual calculation error, we usually use a very small threshold value θ to decide performing Step 2 or Step 3 in Algorithm 4. In other words, if δ<θ, we perform step 2; otherwise, we perform Step 3.


## 4. Incremental Learning of KPCA

In this section, we will outline the incremental learning method of KPCA based on the incremental version of PCA (IPCA) proposed by Hall et al. [18]. The key difference between our method and Hall et al.’s method is that the dimensionality of the projection vector yii=1,2,…,N in our case is not constant; hence we further adapt IPCA to our aim and address this limitation.


### 4.1. Description of IKPCA

Given a sample setφxii=1N and its corresponding projection vector set yii=1N (see the Section 3), we assume we have already built a set of eigenvectors ull=1p and its corresponding matrix U=u1,u2,…,up with the yii=1N set as an input. Note that we have p≤r where r denotes the dimension of yii=1,2,…,N. Let Λ=diag⁡λ1,λ2,…,λp denote the corresponding matrix of eigenvalues and y- the mean vector. Incremental building requires updating these eigenvectors when a new input sample yN+1 is obtained, which is the projection vector of φxN+1. Obviously, the dimensionality of yN+1 may be larger than that of yii=1,2,…,N (see (12)). When their dimensionalities are identical, we denote rN+1=r, otherwise, rN+1≠r. Firstly, we update the mean:(15)y-=1N+1Ny-+yN+1rN+1=r1N+1Ny-0+yN+1rN+1≠r,where y-,0T means adding one zero to the original vector y-. Then we update the set of eigenvectors ull=1pby adding a new vector yN+1 and applying a rotational transformation. In order to do this, we first compute the orthogonal residual vector:(16)hN+1=UηN+1+y--yN+1rN+1=rU0ηN+1+y--yN+1rN+1≠r,where ηN+1 is computed by(17)ηN+1=UTyN+1-y-rN+1=rU0TyN+1-y-rN+1≠r.Subsequently, we normalizehN+1 to obtain h^N+1=hN+1/hN+1 for hN+1>0 and h^N+1=0 otherwise. The new matrix of eigenvectors U′ is computed by(18)U′=U,h^N+1TrN+1=rNU0,h^N+1TrN+1≠rN,where T∈Rp+1×p+1 is a rotation matrix with dimensionp+1. T is the solution of the eigenproblem of the following form [18, 36]:(19)HT=TΛ′.We composeH∈Rp+1×p+1 as(20)H=NN+1Λ00T0+NN+12ηηTτητηTτ2.where τ=h^N+1TyN+1-y- and η=ηN+1.Broadly, the procedure of our incremental method is similar to IPCA presented by Hall et al. [18]. Only under the condition rN+1≠rN, we add one zero (zero vector) to the corresponding variant (matrix).Once we determined the principal direction setull=1p, for a test sample x, the projection of φx onto thel-th nonlinear principal direction vl can be obtained using the formulas in (8) and (10):(21)dlx=vlTφx=ulTβ1,β2,…,βrTφx=ulTDTkbx=ulTy.


### 4.2. Framework of TP-IKPCA

Based on the analysis in Section3 and in Section 4.1, we present the flowchart of TP-IKPCA in Figure 2. Here, xii=1Ndenotes the input sample set and φ(xi)i=1N represents its corresponding mapped set by an implicit nonlinear mapping φ. Firstly, an orthonormal basis βjj=1r of the subspace spanned by φ(xi)i=1N (see (9) or (11)) and the corresponding projection vectors set yii=1N (see (10)) are obtained. Subsequently, a set of eigenvectors ull=1p of yii=1N are built. Then, for a new coming sample xN+1, its mapped sample φ(xN+1) is used to update the orthonormal basis βjj=1r and its projection vector yN+1 is computed using Algorithm 4. Finally, based on the IKPCA described in Section 4.1, the eigenvector set ull=1p is updated using yN+1. It can be seen from Figure 2 that TP-IKPCA includes two main steps. The first step is based on incremental learning algorithm that represents the mapped data using an explicit form (see Algorithm 4). The second step incrementally computes the principal components of yii=1N+1 using IKPCA (see Section 4.1).Figure 2
Flowchart of TP-IKPCA.The dimensionsd,h, and r of the input, kernel, and projection spaces, respectively, generally satisfy the following inequalities: d≤h and r≤h (see Section 1). Here, we need to focus on the meaning of p. In this paper,pdenotes the number of the eigenvectors derived from the covariance matrix of the projection vector set yii=1N. In other words,p is the number of the principal components of yii=1N+1. We note that the number of the principal components, that isp,should be also smaller than the dimension r of the projection space. In summary, we have p≤r≤h.


### 4.3. Complexity Analysis of TP-IKPCA

Suppose that given the current sample setxii=1N, the dimension of the subspace spanned by φxii=1N is r. When a new sample xN+1 arrives, TP-IKPCA first conducts Algorithm 4 where the most time-consuming step is to compute DTkbN+1 appearing in (11)-(14). The time complexity for this step is O(r2) since D is an r×r orthogonal transformation matrix and kbN+1 is anN-dimensional vector. Then, TP-IKPCA incrementally computes the principal components using IKPCA (see Section 4.1). In this step, the computation of the eigendecomposition of matrix H∈Rp+1×p+1 in (20) is most time-consuming. We can define the upper bound for its time complexity as O(p3)≤O(r3) since we have p≤r. Considering the above two steps, the overall time complexity of TP-IKPCA in worst case can be estimated as O(r3). In fact, r≪N is usually tenable [26, 37], which makes TP-IKPCA convenient for processing large-scale or online datasets. Meanwhile, TP-IKPCA requires storing an r×r orthogonal transformation matrix D (see (9)), which only involves a space complexity of Or2. Especially in the testing stage, obtaining the principal component of a new sample only needs to calculate the kernel functions between the new sample and the r old training samples that compose the orthonormal basis (see (21)), which results in improving the computing speed.


## 4.1. Description of IKPCA

Given a sample setφxii=1N and its corresponding projection vector set yii=1N (see the Section 3), we assume we have already built a set of eigenvectors ull=1p and its corresponding matrix U=u1,u2,…,up with the yii=1N set as an input. Note that we have p≤r where r denotes the dimension of yii=1,2,…,N. Let Λ=diag⁡λ1,λ2,…,λp denote the corresponding matrix of eigenvalues and y- the mean vector. Incremental building requires updating these eigenvectors when a new input sample yN+1 is obtained, which is the projection vector of φxN+1. Obviously, the dimensionality of yN+1 may be larger than that of yii=1,2,…,N (see (12)). When their dimensionalities are identical, we denote rN+1=r, otherwise, rN+1≠r. Firstly, we update the mean:(15)y-=1N+1Ny-+yN+1rN+1=r1N+1Ny-0+yN+1rN+1≠r,where y-,0T means adding one zero to the original vector y-. Then we update the set of eigenvectors ull=1pby adding a new vector yN+1 and applying a rotational transformation. In order to do this, we first compute the orthogonal residual vector:(16)hN+1=UηN+1+y--yN+1rN+1=rU0ηN+1+y--yN+1rN+1≠r,where ηN+1 is computed by(17)ηN+1=UTyN+1-y-rN+1=rU0TyN+1-y-rN+1≠r.Subsequently, we normalizehN+1 to obtain h^N+1=hN+1/hN+1 for hN+1>0 and h^N+1=0 otherwise. The new matrix of eigenvectors U′ is computed by(18)U′=U,h^N+1TrN+1=rNU0,h^N+1TrN+1≠rN,where T∈Rp+1×p+1 is a rotation matrix with dimensionp+1. T is the solution of the eigenproblem of the following form [18, 36]:(19)HT=TΛ′.We composeH∈Rp+1×p+1 as(20)H=NN+1Λ00T0+NN+12ηηTτητηTτ2.where τ=h^N+1TyN+1-y- and η=ηN+1.Broadly, the procedure of our incremental method is similar to IPCA presented by Hall et al. [18]. Only under the condition rN+1≠rN, we add one zero (zero vector) to the corresponding variant (matrix).Once we determined the principal direction setull=1p, for a test sample x, the projection of φx onto thel-th nonlinear principal direction vl can be obtained using the formulas in (8) and (10):(21)dlx=vlTφx=ulTβ1,β2,…,βrTφx=ulTDTkbx=ulTy.


## 4.2. Framework of TP-IKPCA

Based on the analysis in Section3 and in Section 4.1, we present the flowchart of TP-IKPCA in Figure 2. Here, xii=1Ndenotes the input sample set and φ(xi)i=1N represents its corresponding mapped set by an implicit nonlinear mapping φ. Firstly, an orthonormal basis βjj=1r of the subspace spanned by φ(xi)i=1N (see (9) or (11)) and the corresponding projection vectors set yii=1N (see (10)) are obtained. Subsequently, a set of eigenvectors ull=1p of yii=1N are built. Then, for a new coming sample xN+1, its mapped sample φ(xN+1) is used to update the orthonormal basis βjj=1r and its projection vector yN+1 is computed using Algorithm 4. Finally, based on the IKPCA described in Section 4.1, the eigenvector set ull=1p is updated using yN+1. It can be seen from Figure 2 that TP-IKPCA includes two main steps. The first step is based on incremental learning algorithm that represents the mapped data using an explicit form (see Algorithm 4). The second step incrementally computes the principal components of yii=1N+1 using IKPCA (see Section 4.1).Figure 2
Flowchart of TP-IKPCA.The dimensionsd,h, and r of the input, kernel, and projection spaces, respectively, generally satisfy the following inequalities: d≤h and r≤h (see Section 1). Here, we need to focus on the meaning of p. In this paper,pdenotes the number of the eigenvectors derived from the covariance matrix of the projection vector set yii=1N. In other words,p is the number of the principal components of yii=1N+1. We note that the number of the principal components, that isp,should be also smaller than the dimension r of the projection space. In summary, we have p≤r≤h.


## 4.3. Complexity Analysis of TP-IKPCA

Suppose that given the current sample setxii=1N, the dimension of the subspace spanned by φxii=1N is r. When a new sample xN+1 arrives, TP-IKPCA first conducts Algorithm 4 where the most time-consuming step is to compute DTkbN+1 appearing in (11)-(14). The time complexity for this step is O(r2) since D is an r×r orthogonal transformation matrix and kbN+1 is anN-dimensional vector. Then, TP-IKPCA incrementally computes the principal components using IKPCA (see Section 4.1). In this step, the computation of the eigendecomposition of matrix H∈Rp+1×p+1 in (20) is most time-consuming. We can define the upper bound for its time complexity as O(p3)≤O(r3) since we have p≤r. Considering the above two steps, the overall time complexity of TP-IKPCA in worst case can be estimated as O(r3). In fact, r≪N is usually tenable [26, 37], which makes TP-IKPCA convenient for processing large-scale or online datasets. Meanwhile, TP-IKPCA requires storing an r×r orthogonal transformation matrix D (see (9)), which only involves a space complexity of Or2. Especially in the testing stage, obtaining the principal component of a new sample only needs to calculate the kernel functions between the new sample and the r old training samples that compose the orthonormal basis (see (21)), which results in improving the computing speed.


## 5. Experiments

We evaluated and compared the performance of TP-IKPCA on synthetic and real datasets with several typical KPCA-based approaches in terms of accuracy and time complexity. The comparison methods include(1) conventional batch mode KPCA, (2) incremental KPCA [14, 15] with reduced-set (IKPCA-RS), and (3) the recently developed incremental KPCA [25] based on rank one updates to the eigendecomposition of kernel matrix (INKPCA). The time complexity can be captured by two aspects: (1) time required for learning the training data; (2)r, i.e., the number of orthonormal basis elements of the subspace spanned by the mapping training data. Usually, a smaller r indicates a reduced time complexity of TP-IKPCA (see Section 4.3). At the same time, we also used two different measures to evaluate the accuracy of TP-IKPCA. The first one is the correlation coefficient between two corresponding principal components (PCs) of TP-IKPCA and KPCA. Since KPCA is performed using all training data in a batch learning model, the PCs of KPCA are the target that TP-IKPCA needs to capture. Ideally, the PCs of TP-IKPCA should be identical with those of KPCA. Therefore, the correlation coefficient between two corresponding PCs is evaluated to show how accurate is TP-IKPCA in comparison to batch KPCA. Specifically, let vl denote thel-th PC of KPCA after learning all of the N samples and vlTPi is thel-th PC of TP-IKPCA after learning ii≤N samples, we define the correlation coefficient (corr) by (22)corrvl,vlTPi=vlTvlTPivlvlTPi,The specific computation in (22) can be deduced from (5) and (8). The second measure is to compare the effectiveness of TP-IKPCA and KPCA. Here, for two-dimensional synthesized datasets, we adopt the contour lines of PCs as an evaluation measure. For real datasets, we adopt the practical denoising effect.


### 5.1. Synthesized Data

In this experiment, we use two-dimensional nonlinear synthetized data to evaluate the accuracy and memory space efficiency of KPCA, IKPCA-RS, INKPCA, and our proposed TP-IKPCA. The data is generated by:(23)y=x22+1+0.2ξ,ξ~N0,1,x~U0,1,where N0,1 denotes the standard Gaussian distribution and U0,1 is the uniform distribution in [0,1]. In this experiment, the kernel function is the polynomial kernel form, namely, kxi,xj=xiTxjm. Here, the parameterm=2.The contour lines of the first three PCs obtained by each method forN=500 training samples are illustrated in Figure 3 (Top to bottom: KPCA, TP-IKPCA, IKPCA-RS, and INKPCA), where red dots represent samples and green lines represent the contour lines of the corresponding PCs. λii=1,2,3 denotes the eigenvalue of the corresponding PC. Figure 3 shows no visually discernible differences between our method and the 3 comparison methods. In other words, their contour lines are very similar. Furthermore, the differences in eigenvalues λi are also very small. Therefore, it can be concluded that all of the results derived from the three incremental algorithms: IKPCA-RS, INKPCA, and our proposed TP-IKPCA, follow the ground truth results closely.Figure 3
Synthesized data including 500 samples and the contours of the first three principal components drawn using a polynomial kernel. The first row is from KPCA, the second row is from TP-IKPCA, the third row is from IKPCA-RS, and the forth row is from INKPCA. Data points are represented by red dots “∗” and the green lines are the contour lines of constant value of the first three principal components.Figure4 shows the evolution curves of the correlation coefficients between the first three PCs obtained by TP-IKPCA, IKPCA-RS, INKPCA, and their corresponding PCs obtained by KPCA when increasing the number of training samples. Figure 4 shows for different incremental algorithms that the resulting correlation coefficient gradually converges to 1 as the number of training samples increases. It indicates that the PCs obtained by TP-IKPCA, IKPCA-RS, and INKPCA gradually approximate those obtained by batch mode KPCA with high accuracy.Figure 4
Evolution of the correlation coefficients between the first three PCs of three incremental algorithms and KPCA when increasing the number of training samples.
(a)
TP-IKPCA (b)
IKPCA-RS (c)
INKPCATable1 shows the average training time for learning from 500 training samples and testing time for extracting features from 100 testing samples (in seconds). We repeated this procedure 20 times and reported the averaged training and testing times as well as the corresponding standard deviations.Table 1
Learning time (sec) for proposed and comparison methods.
Training stage Testing stage KPCA 1.067 ± 0.03 0.203 ± 0.006 IKPCA-RS 0.201±0.0032 0.007±0.0003 INKPCA 0.145±0.0041 0.201 ± 0.004 TP-IKPCA 0.087 ± 0.0026 0.004 ± 0.0002 Notice The training numberN=500 and the basis number r=9 in TP-IKPCA.From Table1, we can clearly see that when using the synthesized data, all of the three incremental KPCA algorithms have faster training speed than KPCA due to low-dimensional features as well as simple distribution of this data. Among these incremental variants, our proposed TP-IKPCA leads to fastest training speed. From the perspective of testing speed, our method can perform prediction in the least time. This can be explained by the fact that the learning time of TP-IKPCA depends on the size of orthonormal basisr (=9) while KPCA and INKPCA both depend on the total number of training samples (=500 in our case). IKPCA-RS also leads to fast testing speed since it uses a reduced set. However, it performs slower than our method in training speed since seeking a set of approximate preimages when new samples arrive using certain optimization techniques or fixed-point iteration is time-consuming.In what follows, we design two experiments to investigate the behavior of our algorithm when the number of training samples increases. Firstly, in Figure5, we plot the variation curve of the number of basis, i.e.,r, against the number of training samples (N). From this result, we notice that in the beginning,r gradually increases with N. However, after N≥21,r stops increasing. This shows that the mapped data have strong linear correlation in the kernel space. Hence, although the number of training samples continues to increase, the number of basis remains stable. More importantly, the computational complexity of TP-IKPCA becomes a constant O93 when N≥21, which is a significant improvement over standard KPCA (in the order of N3) —particularly when the number of training samples becomes very large.Figure 5
Variation of the number of basisr with respect to the number of training samples N.Then, we compute the acceleration ratio which represents the ratio of the time consumed by KPCA to extract features from 100 test samples to the time consumed by TP-IKPCA. The resulting variation of acceleration ratio for testing speed with respect to the number of training samples is shown in Figure6. Obviously, a larger ratio indicates a faster test speed of TP-IKPCA compared with KPCA. Figure 6 also shows that the larger the number of training samples, the larger the ratio, implying that TP-IKPCA can significantly improve test speed compared with KPCA.Figure 6
Changes of the ratio of the test speed with respect to the training sample numbersN.


### 5.2. MNIST Data

In this section, we consider an image processing application where we process the MNIST database of handwritten digits (http://yann.lecun.com/exdb/mnist/). The database consists of handwritten digits from 0 to 9. Each digit set includes training set and testing set. Each digit is represented by a 28-by-28 image. In order to evaluate the performance of different approaches, we carried out image denoising experiment to the even number. Firstly, for each digit, we randomly select 500 training samples and 100 testing samples and then add the Gaussian noise and the salt-and-pepper noise, respectively, to the testing images. The mean and variance for the Gaussian noise are 0 and 0.2, respectively, while the level of the salt-and-pepper noise is 0.3. Next, we estimate the first sixteen principal components using KPCA, IKPCA-RS, INKPCA, and our proposed TP-IKPCA, respectively. Finally, we perform denoising experiments on all corrupted testing samples and reconstruct them using the reconstructive scheme presented in [38]. In the experiment, we use the Gaussian function kxi,xj=exp-xi-xj2/σ2 as the kernel function in each method and set the bandwidth σ=0.2.Figure7 shows some restoration results from the corrupted images by conducting KPCA, TP-IKPCA, IKPCA-RS and INKPCA, respectively. It can be seen from these figures that all of these methods can eliminate noise and reconstruct the images well. Therefore, we can draw the conclusion that different incremental KPCA algorithms can closely approximate batch mode KPCA in reconstruction performance with high accuracy.Figure 7
Restoration results by TP-IKPCA, IKPCA-RS, INKPCA, and KPCA.Figure8 shows the evolution curves of the correlation coefficients between the first three PCs by TP-IKPCA, IKPCA-RS, INKPCA, and their corresponding ones by KPCA when the number of training samples increases, where the horizontal axis denotes the number of training samples and the vertical axis is the correlation coefficient computed by (22). Figure 8 shows that all incremental KPCA algorithms lead to good approximation accuracies for the first three PCs since the correlation coefficients gradually converge to 1 as the number of training samples increases. The results indicate that incremental KPCA algorithms can gradually capture PCs in real high-dimensional datasets with good approximation accuracy.Figure 8
Evolution curves of the correlation coefficients between the first three PCs of (a) TP-IKPCA, (b) IKPCA-RS, (c) INKPCA, and KPCA as the number of training samples increases (x-axis).
(a)
TP-IKPCA (b)
IKPCA-RS (c)
INKPCAWe also display in Table2 the average training and testing times (in seconds) for training 500 samples and extracting features from 100 testing samples, respectively. We repeated this process 20 times and reported the average values and the corresponding standard deviations.Table 2
Experimental results of learning time on different digit (sec).
Training stage Testing stage “0” “2” “4” “6” “8” “0” “2” “4” “6” “8”  KPCA 1.75 ± 0.055 1.74 ± 0.018 1.74 ± 0.021 1.76 ± 0.036 1.74 ± 0.030 0.358 ± 0.022 0.351 ± 0.013 0.357 ± 0.012 0.352 ± 0.011 0.354 ± 0.011  IKPCA-RS 43.68 ± 0.049 41.32 ± 0.038 43.21 ± 0.068 44.75 ± 0.074 42.87 ± 0.059 0.221 ± 0.009 0.243± 0.011 0.208 ± 0.014 0.169 ± 0.017 0.214 ± 0.008  INKPCA 23.67 ± 0.032 22.13 ± 0.038 23.91 ± 0.041 23.10 ± 0.047 22.13 ± 0.035 0.360 ± 0.023 0.358 ± 0.014 0.353 ± 0.012 0.355 ± 0.012 0.352 ± 0.011  TP-IKPCA 1.89± 0.019 3.06± 0.017 1.67± 0.043 1.40± 0.024 2.27± 0.031 0.181± 0.012 0.251 ± 0.008 0.162± 0.016 0.144± 0.013 0.199± 0.010  r 242 351 219 190 282 242 351 219 190 282Firstly, we analyze the training results shown in Table2. We can derive the following observations: (1) IKPCA-RS consumes much more time in training than other approaches, including batch model KPCA. These results differ from those obtained when using the synthesized data where IKPCA-RS ran faster than batch mode KPCA. Through inspecting the experiments, we found that IKPCA-RS iterates reduced-set expansions many times when computing the preimages for compression. As a result, when handling data with a large number of features (=784 in our case), the calculation of preimages increases the computational load. This observation is in line with the conclusion drawn in [15]. (2) INKPCA needs to apply rank one update to iteratively calculate eigenvalues and eigenvectors associated with kernel matrix when new data arrive. Although INKPCA consumes more time than batch mode KPCA, it does not need to store the whole kernel matrix in memory and has the advantage of better handling massive data where KPCA rapidly becomes infeasible. (3) As for our proposed TP-IKPCA, it runs much faster than other incremental algorithms, including IKPCA-RS and INKPCA, excluding a few cases where our algorithm may be slower than batch mode KPCA. This can be explained by the computational complexity of TP-IKPCA of the order of r3. In this experiment, the linear correlation between the mapped digital images is very weak, which results in a very large r and even approximates the number of training samples.From the testing results shown in Table2, we can conclude the following: (1) The testing speeds of KPCA and INKPCA are similar since INKPCA cannot select a few yet important samples from the whole data but still makes use of all available samples when calculating the projections of new data. Therefore, their speed is proportional to the total number of the training samples. (2) Regardless of specific digit, the testing speed of IKPCA-RS and TP-IKPCA are much faster than that of KPCA and INKPCA since both methods are able to reduce the number of samples used for kernel evaluation although they adopt different strategies. The testing time of TP-IKPCA is proportional to the size of basis r. In this experiment, the size of basis r is smaller than the number of training samples, thus leading to an improvement in the test speed compared with KPCA. Considering the training time, our proposed TP-IKPCA is obviously preferred over IKPCA-RS.In what follows, we gradually increase the number of training samples and summarize the training and testing time required by TP-IKPCA and standard KPCA. We find from extensive experiments that the computational superiority of TP-IKPCA over KPCA increases with the number of training samples. Taking the experiments on digit “0” as an example, we repeated this evaluation 20 times and recorded the resulting training and testing time (in seconds) required by TP-IKPCA and KPCA under different training sample sizeN. Finally, the averaged time and the standard deviation are shown in Table 3 where the training ratio in Table 3 denotes the ratio of training time of KPCA to that of TP-IKPCA, given a total number of N samples. In a similar way, the testing ratio represents the ratio of testing time of KPCA to that of TP-IKPCA when extracting features from 100 test samples.Table 3
Training and testing time (sec) on digit “0” when increasing the number of training samples from 500 to 5000.
N 500 1000 1500 2000 3000 4000 5000 r 242 366 453 514 605 670 735  N/r 2.07 2.73 3.31 3.89 4.96 5.97 6.80  training stage KPCA 1.75 ± 0.055 6.98 ± 0.066 18.72 ± 0.140 35.46 ± 0.185 79.97 ± 0.570 154.65 ± 0.687 255.01 ± 2.083 TP-IKPCA 1.89 ± 0.019 5.34 ± 0.064 9.86 ± 0.085 14.92 ± 0.078 25.40 ± 0.015 37.46 ± 0.031 50.43 ± 0.411 ratio 0.93 1.31 1.90 2.38 3.15 4.13 5.06  testing stage KPCA 0.358 ± 0.022 0.650 ± 0.013 1.131 ± 0.018 1.600 ± 0.026 2.406 ± 0.070 3.657 ± 0.170 5.234 ± 0.293 TP-IKPCA 0.181 ± 0.012 0.266 ± 0.014 0.323 ± 0.017 0.371 ± 0.016 0.429 ± 0.429 0.465 ± 0.011 0.511 ± 0.010 ratio 1.98 2.44 4.06 4.31 5.61 8.02 10.24Based on Table3, we can make the following observations: (1) As the number of training samples continues to increase, the number of basis r tends to increase as well. However, the increasing speed of r gradually decreases, which indicates that the larger the size of the training set, the stronger the correlation between the samples. (2) With the increasing of the training set sizeN, the training time of both KPCA and TP-IKPCA also increases gradually. However, the increasing speed of TP-IKPCA is much slower than that of KPCA. The reason lies in that the time complexity of TP-IKPCA has a close relationship with the number of basis r while that of KPCA depends on the total number of training samples N. AsN increases, the increasing speed of r progressively decreases since most of the correlation structure among data has been revealed. We also derive a similar conclusion from the ratio’s evolution in the training stage with respect to the changes of the number of training samples, where the ratio gradually increases with N in the training stage. (3) As N increases, the testing time of KPCA and TP-IKPCA both increase gradually. However, the increasing speed of TP-IKPCA is much slower than that of KPCA, which is also reflected by the ratio’s evolution in the testing stage. The reason is that the test speed of TP-IKPCA is closely related to the number of basis r and the increasing speed of r gradually becomes slower with the increasing of N. Based on the above analysis, we conclude that TP-IKPCA does significantly improve the computational complexity of KPCA. Moreover, TP-IKPCA can deal with dynamic dataset due to its “incremental” nature.


## 5.1. Synthesized Data

In this experiment, we use two-dimensional nonlinear synthetized data to evaluate the accuracy and memory space efficiency of KPCA, IKPCA-RS, INKPCA, and our proposed TP-IKPCA. The data is generated by:(23)y=x22+1+0.2ξ,ξ~N0,1,x~U0,1,where N0,1 denotes the standard Gaussian distribution and U0,1 is the uniform distribution in [0,1]. In this experiment, the kernel function is the polynomial kernel form, namely, kxi,xj=xiTxjm. Here, the parameterm=2.The contour lines of the first three PCs obtained by each method forN=500 training samples are illustrated in Figure 3 (Top to bottom: KPCA, TP-IKPCA, IKPCA-RS, and INKPCA), where red dots represent samples and green lines represent the contour lines of the corresponding PCs. λii=1,2,3 denotes the eigenvalue of the corresponding PC. Figure 3 shows no visually discernible differences between our method and the 3 comparison methods. In other words, their contour lines are very similar. Furthermore, the differences in eigenvalues λi are also very small. Therefore, it can be concluded that all of the results derived from the three incremental algorithms: IKPCA-RS, INKPCA, and our proposed TP-IKPCA, follow the ground truth results closely.Figure 3
Synthesized data including 500 samples and the contours of the first three principal components drawn using a polynomial kernel. The first row is from KPCA, the second row is from TP-IKPCA, the third row is from IKPCA-RS, and the forth row is from INKPCA. Data points are represented by red dots “∗” and the green lines are the contour lines of constant value of the first three principal components.Figure4 shows the evolution curves of the correlation coefficients between the first three PCs obtained by TP-IKPCA, IKPCA-RS, INKPCA, and their corresponding PCs obtained by KPCA when increasing the number of training samples. Figure 4 shows for different incremental algorithms that the resulting correlation coefficient gradually converges to 1 as the number of training samples increases. It indicates that the PCs obtained by TP-IKPCA, IKPCA-RS, and INKPCA gradually approximate those obtained by batch mode KPCA with high accuracy.Figure 4
Evolution of the correlation coefficients between the first three PCs of three incremental algorithms and KPCA when increasing the number of training samples.
(a)
TP-IKPCA (b)
IKPCA-RS (c)
INKPCATable1 shows the average training time for learning from 500 training samples and testing time for extracting features from 100 testing samples (in seconds). We repeated this procedure 20 times and reported the averaged training and testing times as well as the corresponding standard deviations.Table 1
Learning time (sec) for proposed and comparison methods.
Training stage Testing stage KPCA 1.067 ± 0.03 0.203 ± 0.006 IKPCA-RS 0.201±0.0032 0.007±0.0003 INKPCA 0.145±0.0041 0.201 ± 0.004 TP-IKPCA 0.087 ± 0.0026 0.004 ± 0.0002 Notice The training numberN=500 and the basis number r=9 in TP-IKPCA.From Table1, we can clearly see that when using the synthesized data, all of the three incremental KPCA algorithms have faster training speed than KPCA due to low-dimensional features as well as simple distribution of this data. Among these incremental variants, our proposed TP-IKPCA leads to fastest training speed. From the perspective of testing speed, our method can perform prediction in the least time. This can be explained by the fact that the learning time of TP-IKPCA depends on the size of orthonormal basisr (=9) while KPCA and INKPCA both depend on the total number of training samples (=500 in our case). IKPCA-RS also leads to fast testing speed since it uses a reduced set. However, it performs slower than our method in training speed since seeking a set of approximate preimages when new samples arrive using certain optimization techniques or fixed-point iteration is time-consuming.In what follows, we design two experiments to investigate the behavior of our algorithm when the number of training samples increases. Firstly, in Figure5, we plot the variation curve of the number of basis, i.e.,r, against the number of training samples (N). From this result, we notice that in the beginning,r gradually increases with N. However, after N≥21,r stops increasing. This shows that the mapped data have strong linear correlation in the kernel space. Hence, although the number of training samples continues to increase, the number of basis remains stable. More importantly, the computational complexity of TP-IKPCA becomes a constant O93 when N≥21, which is a significant improvement over standard KPCA (in the order of N3) —particularly when the number of training samples becomes very large.Figure 5
Variation of the number of basisr with respect to the number of training samples N.Then, we compute the acceleration ratio which represents the ratio of the time consumed by KPCA to extract features from 100 test samples to the time consumed by TP-IKPCA. The resulting variation of acceleration ratio for testing speed with respect to the number of training samples is shown in Figure6. Obviously, a larger ratio indicates a faster test speed of TP-IKPCA compared with KPCA. Figure 6 also shows that the larger the number of training samples, the larger the ratio, implying that TP-IKPCA can significantly improve test speed compared with KPCA.Figure 6
Changes of the ratio of the test speed with respect to the training sample numbersN.


## 5.2. MNIST Data

In this section, we consider an image processing application where we process the MNIST database of handwritten digits (http://yann.lecun.com/exdb/mnist/). The database consists of handwritten digits from 0 to 9. Each digit set includes training set and testing set. Each digit is represented by a 28-by-28 image. In order to evaluate the performance of different approaches, we carried out image denoising experiment to the even number. Firstly, for each digit, we randomly select 500 training samples and 100 testing samples and then add the Gaussian noise and the salt-and-pepper noise, respectively, to the testing images. The mean and variance for the Gaussian noise are 0 and 0.2, respectively, while the level of the salt-and-pepper noise is 0.3. Next, we estimate the first sixteen principal components using KPCA, IKPCA-RS, INKPCA, and our proposed TP-IKPCA, respectively. Finally, we perform denoising experiments on all corrupted testing samples and reconstruct them using the reconstructive scheme presented in [38]. In the experiment, we use the Gaussian function kxi,xj=exp-xi-xj2/σ2 as the kernel function in each method and set the bandwidth σ=0.2.Figure7 shows some restoration results from the corrupted images by conducting KPCA, TP-IKPCA, IKPCA-RS and INKPCA, respectively. It can be seen from these figures that all of these methods can eliminate noise and reconstruct the images well. Therefore, we can draw the conclusion that different incremental KPCA algorithms can closely approximate batch mode KPCA in reconstruction performance with high accuracy.Figure 7
Restoration results by TP-IKPCA, IKPCA-RS, INKPCA, and KPCA.Figure8 shows the evolution curves of the correlation coefficients between the first three PCs by TP-IKPCA, IKPCA-RS, INKPCA, and their corresponding ones by KPCA when the number of training samples increases, where the horizontal axis denotes the number of training samples and the vertical axis is the correlation coefficient computed by (22). Figure 8 shows that all incremental KPCA algorithms lead to good approximation accuracies for the first three PCs since the correlation coefficients gradually converge to 1 as the number of training samples increases. The results indicate that incremental KPCA algorithms can gradually capture PCs in real high-dimensional datasets with good approximation accuracy.Figure 8
Evolution curves of the correlation coefficients between the first three PCs of (a) TP-IKPCA, (b) IKPCA-RS, (c) INKPCA, and KPCA as the number of training samples increases (x-axis).
(a)
TP-IKPCA (b)
IKPCA-RS (c)
INKPCAWe also display in Table2 the average training and testing times (in seconds) for training 500 samples and extracting features from 100 testing samples, respectively. We repeated this process 20 times and reported the average values and the corresponding standard deviations.Table 2
Experimental results of learning time on different digit (sec).
Training stage Testing stage “0” “2” “4” “6” “8” “0” “2” “4” “6” “8”  KPCA 1.75 ± 0.055 1.74 ± 0.018 1.74 ± 0.021 1.76 ± 0.036 1.74 ± 0.030 0.358 ± 0.022 0.351 ± 0.013 0.357 ± 0.012 0.352 ± 0.011 0.354 ± 0.011  IKPCA-RS 43.68 ± 0.049 41.32 ± 0.038 43.21 ± 0.068 44.75 ± 0.074 42.87 ± 0.059 0.221 ± 0.009 0.243± 0.011 0.208 ± 0.014 0.169 ± 0.017 0.214 ± 0.008  INKPCA 23.67 ± 0.032 22.13 ± 0.038 23.91 ± 0.041 23.10 ± 0.047 22.13 ± 0.035 0.360 ± 0.023 0.358 ± 0.014 0.353 ± 0.012 0.355 ± 0.012 0.352 ± 0.011  TP-IKPCA 1.89± 0.019 3.06± 0.017 1.67± 0.043 1.40± 0.024 2.27± 0.031 0.181± 0.012 0.251 ± 0.008 0.162± 0.016 0.144± 0.013 0.199± 0.010  r 242 351 219 190 282 242 351 219 190 282Firstly, we analyze the training results shown in Table2. We can derive the following observations: (1) IKPCA-RS consumes much more time in training than other approaches, including batch model KPCA. These results differ from those obtained when using the synthesized data where IKPCA-RS ran faster than batch mode KPCA. Through inspecting the experiments, we found that IKPCA-RS iterates reduced-set expansions many times when computing the preimages for compression. As a result, when handling data with a large number of features (=784 in our case), the calculation of preimages increases the computational load. This observation is in line with the conclusion drawn in [15]. (2) INKPCA needs to apply rank one update to iteratively calculate eigenvalues and eigenvectors associated with kernel matrix when new data arrive. Although INKPCA consumes more time than batch mode KPCA, it does not need to store the whole kernel matrix in memory and has the advantage of better handling massive data where KPCA rapidly becomes infeasible. (3) As for our proposed TP-IKPCA, it runs much faster than other incremental algorithms, including IKPCA-RS and INKPCA, excluding a few cases where our algorithm may be slower than batch mode KPCA. This can be explained by the computational complexity of TP-IKPCA of the order of r3. In this experiment, the linear correlation between the mapped digital images is very weak, which results in a very large r and even approximates the number of training samples.From the testing results shown in Table2, we can conclude the following: (1) The testing speeds of KPCA and INKPCA are similar since INKPCA cannot select a few yet important samples from the whole data but still makes use of all available samples when calculating the projections of new data. Therefore, their speed is proportional to the total number of the training samples. (2) Regardless of specific digit, the testing speed of IKPCA-RS and TP-IKPCA are much faster than that of KPCA and INKPCA since both methods are able to reduce the number of samples used for kernel evaluation although they adopt different strategies. The testing time of TP-IKPCA is proportional to the size of basis r. In this experiment, the size of basis r is smaller than the number of training samples, thus leading to an improvement in the test speed compared with KPCA. Considering the training time, our proposed TP-IKPCA is obviously preferred over IKPCA-RS.In what follows, we gradually increase the number of training samples and summarize the training and testing time required by TP-IKPCA and standard KPCA. We find from extensive experiments that the computational superiority of TP-IKPCA over KPCA increases with the number of training samples. Taking the experiments on digit “0” as an example, we repeated this evaluation 20 times and recorded the resulting training and testing time (in seconds) required by TP-IKPCA and KPCA under different training sample sizeN. Finally, the averaged time and the standard deviation are shown in Table 3 where the training ratio in Table 3 denotes the ratio of training time of KPCA to that of TP-IKPCA, given a total number of N samples. In a similar way, the testing ratio represents the ratio of testing time of KPCA to that of TP-IKPCA when extracting features from 100 test samples.Table 3
Training and testing time (sec) on digit “0” when increasing the number of training samples from 500 to 5000.
N 500 1000 1500 2000 3000 4000 5000 r 242 366 453 514 605 670 735  N/r 2.07 2.73 3.31 3.89 4.96 5.97 6.80  training stage KPCA 1.75 ± 0.055 6.98 ± 0.066 18.72 ± 0.140 35.46 ± 0.185 79.97 ± 0.570 154.65 ± 0.687 255.01 ± 2.083 TP-IKPCA 1.89 ± 0.019 5.34 ± 0.064 9.86 ± 0.085 14.92 ± 0.078 25.40 ± 0.015 37.46 ± 0.031 50.43 ± 0.411 ratio 0.93 1.31 1.90 2.38 3.15 4.13 5.06  testing stage KPCA 0.358 ± 0.022 0.650 ± 0.013 1.131 ± 0.018 1.600 ± 0.026 2.406 ± 0.070 3.657 ± 0.170 5.234 ± 0.293 TP-IKPCA 0.181 ± 0.012 0.266 ± 0.014 0.323 ± 0.017 0.371 ± 0.016 0.429 ± 0.429 0.465 ± 0.011 0.511 ± 0.010 ratio 1.98 2.44 4.06 4.31 5.61 8.02 10.24Based on Table3, we can make the following observations: (1) As the number of training samples continues to increase, the number of basis r tends to increase as well. However, the increasing speed of r gradually decreases, which indicates that the larger the size of the training set, the stronger the correlation between the samples. (2) With the increasing of the training set sizeN, the training time of both KPCA and TP-IKPCA also increases gradually. However, the increasing speed of TP-IKPCA is much slower than that of KPCA. The reason lies in that the time complexity of TP-IKPCA has a close relationship with the number of basis r while that of KPCA depends on the total number of training samples N. AsN increases, the increasing speed of r progressively decreases since most of the correlation structure among data has been revealed. We also derive a similar conclusion from the ratio’s evolution in the training stage with respect to the changes of the number of training samples, where the ratio gradually increases with N in the training stage. (3) As N increases, the testing time of KPCA and TP-IKPCA both increase gradually. However, the increasing speed of TP-IKPCA is much slower than that of KPCA, which is also reflected by the ratio’s evolution in the testing stage. The reason is that the test speed of TP-IKPCA is closely related to the number of basis r and the increasing speed of r gradually becomes slower with the increasing of N. Based on the above analysis, we conclude that TP-IKPCA does significantly improve the computational complexity of KPCA. Moreover, TP-IKPCA can deal with dynamic dataset due to its “incremental” nature.


## 6. Conclusion

In this paper, we proposed a novel incremental feature extraction method termed as TP-IKPCA which endowed KPCA with the capability of handling dynamic or large-scale datasets. The proposed TP-IKPCA differs from the existing incremental approaches in providing an explicit form of the mapped data and the updating process of KPCs is also performed in an explicit space. Specifically, TP-IKPCA is implemented in two phases. First, an incremental algorithm is given to explicitly project the mapped samples in the kernel space. Second, we employed the existing incremental method of PCA to capture KPCs based on the explicit data in the projection space. The computational complexity of TP-IKPCA has a close relationship with the size of basisr of the subspace spanned by the mapped training samples. Usually,r is much smaller than the number of training samplesN, and thus TP-IKPCA can greatly improve the computational complexity of KPCA. In the case of large-scale or online dataset, the computational superiority of our approach is remarkable. Experimental results on synthetic and real datasets demonstrate that TP-IKPCA can significantly improve the time complexity of KPCA while preserving a high accuracy as standard KPCA. In comparison with two incremental KPCA algorithms, TP-IKPCA also illustrates superiority in terms of training and testing speed.TP-IKPCA can be utilized in any application where KPCA needs to be conducted, especially when training data is of large scale, or can only be collected one by one, where the conventional batch-based KPCA cannot be applied. The idea of this study can be extended to other kernel-based methods, such as Kernel Fisher discriminant analysis (KFDA), Kernel independent component analysis (KICA), and so on.



---

*Source: 5937274-2019-02-11.xml*