Tutorial on PCA and approximate PCA and approximate kernel PCA

Principal Component Analysis (PCA) is one of the most widely used data analysis methods in machine learning and AI. This manuscript focuses on the mathematical foundation of classical PCA and its application to a small-sample-size scenario and a large dataset in a high-dimensional space scenario. In particular, we discuss a simple method that can be used to approximate PCA in the latter case. This method can also help approximate kernel PCA or kernel PCA (KPCA) for a large-scale dataset. We hope this manuscript will give readers a solid foundation on PCA, approximate PCA, and approximate KPCA.


Introduction
What would be the description if you were asked to describe a person sitting next to you? Male or female, wear glasses or not, big mouth, small nose, long hair, tall, overweight, BMI score, etc. These descriptions can be put together in the form of a "vector". In this manuscript, we consider "column vector", a matrix of one column and several rows, each for a specific description. With ten descriptions, we have a vector of 10 rows. We say that this vector lives in a 10-dimensional space. With 100 descriptions, we have a vector of 100 rows that lives in a 100-dimensional space.
A vector is a mathematical object composed of "length" or "magnitude" and "direction". Nonetheless, in this paper, we will use the term "vector" and "point" interchangeably since, given a point, we can always assign the vector from the origin point to that given point. Note that it is common to work in a very high-dimensional space in AI, machine learning, and data science. Some intuitions we have gained from our daily life experience in two or three dimensions may not be correct in this very high-dimensional space. Readers interested in properties of high-dimensional data may refer to Blum et al. (2017) for more detail.
Suppose we have prepared 100 descriptions or features to describe each person. Suppose that these features were all numerical, e.g., the height, the weight, the BMI score, etc. Suppose 1 3 that we sat at a coffee shop and tried to describe each client using our 100 features. For each client, we obtain a feature vector that contains all descriptions assigned to him/her. Given a collection of feature vectors, what shall we see ? We will see that the observed height values do not go from 0 to infinity. They will be bounded within some range, for example, between 150 to 180 cm. The observed weights and BMI scores are bounded as well. The observed feature vectors will not be scattered all over the whole space. They are clustered around some area, or subspace, within the big multi-dimensional space. With a deeper inspection, you may find correlations between features such as height, weight, and BMI score. Therefore, it is tempting to combine these features into a new one. The question is then how to combine these features and measure if the combination is good.
Principal Component Analysis (PCA) provides such a combination method. Indeed, PCA relies on linear combination of these features to construct principal subspace that is the main subspace on which most of the feature vectors lie. PCA uses variance as a measure of the information content of the subspace. The quality of the selected subspace can be measured by comparing the variance of data within it to the total variance of the whole dataset. These are two keys idea of PCA. In practice, however, one may run into practical problems such as having too few observations compared to the data dimension or having too much computational cost due to the high dimensionality of data and a large number of observations, etc. This paper describes some practical algorithms for handling these different scenarios. Some material presented hereafter is extended from our previously published paper (Marukatat 2016). It should be noted that this tutorial focuses mainly on the algorithmic aspect of PCA. Readers interested in its applications are invited to consult the paper (Ringnér 2008;Gewers et al. 2018;Atchade-Adelomou et al. 2022) and references therein to see different use-cases of PCA on different domains including bioinformatics and quantum computing.
In the following, Sect. 2 first describes PCA in the usual setting where the number of available data points is larger than the dimension of feature vectors. Then Sect. 3 describes PCA that is designed for small-sample-size case. This case could happen; for example, when the number of available subjects with a particular disease is limited, the researchers may extract all possible features from their blood or tissue sample. The total amount of features could be easily more significant than the total number of samples. In this case, classical PCA will not be suitable. This section describes a variation of PCA designed to handle this case. Section 4 describes another variation of PCA designed for large-scale datasets in a very high-dimension space. Indeed, current AI and data science datasets are often large with very high-dimension feature vectors. This section explains how to mitigate the computational cost of PCA in this case. The technique introduced in this section, namely dot product preserving transformation, will be used in the subsequent section on approximate kernel PCA. Section 5 discusses kernel PCA or KPCA, a non-linear extension of PCA using kernel trick. The classical KPCA is presented in Sect. 5.1. In this section, one will see that the main crux of KPCA is its computational cost. Two variations of KPCA designed to reduce the cost of KPCA are presented in Sects. 5.2 and 5.3 with some examples shown in Sect. 5.5.

Classical case: n > d
This section describes a classical scenario for applying PCA, namely when the number of available data points is larger than the dimension of feature vectors. In this section, we assume that the dimension of feature vectors is not very large, so all calculations can be done on a modest computer.

Notation
In the following, let ∈ ℝ d denotes a feature vector in d-dimensional space with A feature vector contains values of features that are observed on actual data. We assume that there are no missing values. We further assume that these features are all numerical.
In theory, the features should be selected adequately by domain experts. Nonetheless, in practice, we often rely on low-level features, such as the intensity of each pixel in an image, that can be computed easily. High-level features, like loops and cusps in handwritten strokes, contain more information. However, their detection depends on image processing functions that must be put in place first. These functions may contain some heuristic or programming bias. As a result, for some poorly scanned images, the detection of high-level features may not be done accurately. Low-level features, such as pixel intensities, are more robust because they can always be extracted from all input images. Even if each low-level feature may not convey a meaningful description by itself, combining lots of them with lots of data makes it possible to derive meaningful meaning using machine learning models. PCA is one method to achieve this goal.

Linear projection and dot product
The main idea of PCA is to use variance as a measure of information content and to identify linear subspace maximizing the variance. The linear subspace can be constructed from multiple unit vectors that are orthogonal to each other. Each unit vector represents a projection axis. Consider a projection axis ∈ ℝ d , the projection of ∈ ℝ d onto the axis is the point on this axis that is closest to . This projection can be computed by ‖ ‖ cos , . Recall that the "dot product" or the "scalar product" or the "inner product" between any two vectors and denotes as ⟨ , ⟩ satisfies In the following we will use ⟨ , ⟩ or T interchangeably when talking about dot product. The projection of onto the unit vector axis is then given by a simple dot product T ∈ ℝ. Note that the "projection" refers to the location on the projection axis. The same location can be described relative to the input space as well by multiplying it with the 1 3 projection axis, i.e., ( T ) . This quantity is often referred to as the "reconstruction" of from its projection on the axis .

Variance on projection axis and covariance matrix
Given a set of centered real values z 1 , z 2 , … , z n z t ∈ ℝ , the variance of this dataset is Given a set of feature vectors 1 , … , n , t ∈ ℝ d , and a projection axis (that is a unit vector, i.e. ‖ ‖ = 1 ). The projections on are T 1 , … , T n . Suppose that the feature vectors are centered around zero vectors. Thus, their linear projections will also be centered appropriately. Note that each dot product yields a simple real value, i.e., T t ∈ ℝ . Hence, the Eq. (2) can be used to compute the variance: The last line gives us a shortcut to compute variance on any axis . Indeed, instead of projecting all data onto the axis and then computing the variance, we can pre-compute the covariance matrix first, then given any projection axis the variance of the projection on is given by T .

PCA problem and Lagrangian
From the above calculation, we can formulate the main idea of PCA as the following constrained optimization problem: For standard optimization problems, we may set the partial derivatives to zero and solve the resulting equations or follow the gradient. However, the above problem has an additional constraint that must be considered. To this end, the standard approach is to use the Lagrangian method (see (Dimitri 1996) for more detail on the Lagrangian method).
In short, to solve the following constrained optimization problem: (3) max T subject to T = 1.
the following Lagrangian function is considered: with the Lagrangian coefficient describing the important of the constraint. Then, given the solution ( * , * ) optimizing the function L, * will be the solution of the initial problem.
For the PCA problem, its Lagrangian is given by Then, if we compute the partial derivative of L with respect to and set it to zero, we obtain: The above derivation uses matrix calculus. Readers unfamiliar with this calculation could consult reference (Petersen and Pedersen 2012) for more detail. The Eq. (10) means that the axis maximizing variance must be an eigenvector of , and is its corresponding eigenvalue. The eigenvectors and eigenvalues are, in fact, mathematical objects that have been studied for a very long time. PCA represents another application of this well-founded subject. An example of a commonly used method to extract eigenvector and eigenvalue pair, or eigen-decomposition algorithm called the Jacobi method is shown in Section A.

Eigenvalue, variance and axis selection
From Sect. 2.4, we have seen that the projection axis maximizing variance must be an eigenvector of the covariance matrix . However, as is a d × d matrix, it has d eigenvectors; which one should be selected? To answer this question, we can observed that In other words, the eigenvalue corresponds to the variance of the data projected onto its corresponding eigenvector.
Given eigenvalues 1 ≥ 2 ≥ ... ≥ d with 1 , 2 , … , d the corresponding eigenvectors. From the discussion above, we know that if we select only one projection axis, then 1 should be the best choice since it corresponds to the largest eigenvalue, thus the largest variance of projected data.
Note that from these eigenvalues we can easily compute the total variance, i.e. ∑ d i=1 i . This total variance represents the entire information content of the dataset. By comparing 1 to the total variance, we can judge if a single projection axis 1 is enough to retain the most information or not. In practice, a single projection axis is not enough. Fortunately, as the eigenvectors are already orthogonal, it is possible to select the next "largest" eigenvectors to form the basis of the principal subspace. The next question is how many axes should be selected. There are three approaches to this question: 1. User gives the desired number of axis m. 2. If we assume that p portion of the information content was, in fact, noise, then m should be the smallest number such that We can also assume that all axis with variance smaller than a pre-defined threshold are noise and should be discarded. Thus, we select all axis that i ≥ 1 .
The PCA algorithm is summarized in Algorithm 1. Given principal axis 1 , … , m , each new data point can be represented within the principal subspace as: and the reconstruction of from its projection on the principal subspace is given by: It is pretty straightforward to prove that the selected principal subspace also minimizes the reconstruction error ∑ t ‖ t −̂ t ‖ 2 and that the reconstruction error is indeed the sum of all eigenvalues outside the principal subspace.

Examples
This section shows examples of the application of PCA to various datasets, namely Diabetes 1 , Wisconsin breast cancer 2 , Iris 3 , MNIST 4 , and CIFAR-10 5 . Table 1 summarizes the detail of these datasets Diabetes, Wisconsin breast cancer, and Iris contain multi-dimensional vectors, whereas MNIST and CIFAR contain grayscale images in 28 × 28 pixels and RGB images in 32 × 32 pixels, respectively. These images were reorganized into vectors in 28 × 28 = 784 dimensions and 32 × 32 × 3 = 3072 dimensions respectively. Z-score normalization was applied to all datasets. Thus all data were properly centered. Figure 1, from the top to the bottom row, shows results from five datasets: Diabetes, Wisconsin breast cancer, Iris, MNIST, and CIFAR, respectively. The left column shows the eigenvalues plot from the four datasets. The second column shows the cumulated proportion of variances or eigenvalues. The two right columns show the projection of the five datasets onto the two first eigenvectors, and onto the first and the last eigenvectors, respectively. In both cases, the projection onto the first eigenvector corresponds to the abscissa. One can observe that the variances on the horizontal axis are higher than that on the vertical axis since the eigenvectors are sorted in decreasing order of their eigenvalues, and each eigenvalue corresponds to variance on the corresponding eigenvector. For the MNIST dataset, the rightmost plot appears to be quite strange. This is due to an outlier in the dataset shown in Fig. 2a. Indeed, this image has a small white dot on the left side. This dot disappears in the reconstruction using 50 first principal components (Fig. 2b). However, this reconstruction is blurred. Using more principal components results in a sharper image but with some noise. Hence, it is possible to denoise the input data using the PCA method. Table 2 depicted the number of axes needed to retain 90%, 95%, and 99% of the total information compared with the dimension of these datasets. One can observe that for some  datasets, e.g., Iris and CIFAR, the number of the required axis to retain 90% of information is only a small fraction of the total dimension. Other datasets, like Diabetes, requires a larger number of axis.

Discussion
The principal subspace obtained from PCA describes the area in the whole space on which most of the data lies. This knowledge is helpful. For example, consider the problem of face verification. We are given a pair of face images and have to decide if both images come from the same person or not. If we compare both images directly, in pixel space, so to speak, the observed difference may be due to the background area in the images, which does not reflect the identity of the person. By working on the principal subspace of the face images, we can be sure that every comparison is meaningful. This was the idea behind Eigenface that was the foundation of early face recognition methods (Sirovich and Kirby 1987; Turk and Pentland 1991). Today's face verification methods are based on deep features extracted from deep learning models (Wang and Deng 2018). This approach is possible now because we have more data to train deep neural networks. In the early period of face recognition, the classical face dataset called ORL (orl database) contained only 40 subjects with ten images per subject. The size of each image was 92 × 112 pixels, hence 10,304 pixels in total. This dataset is tiny compared to today's datasets. Besides, the main problem, in this case, is that we have more dimensions than the number of observations ( d > n ). Consequently, the covariance matrix will most likely be ill conditioned and cannot be eigen-decomposed. The Eigenface method resolves this problem via dot matrix instead of the covariance matrix. This procedure will be described in the next section.

Small-sample-size case: n ≤ d
This section describes the PCA method designed for the case where the number of data points (n) is smaller than the dimension of feature vectors (d). Even if the method described in this section originated from image application, it could also be used with general feature vectors when n ≤ d . This scenario appears, for example, in bioinformatics (Ringnér 2008). In Ringnér (2008), the author studied 8,534 probes on the microarrays with expression measurements extracted from 105 samples. In other words, the authors consider a dataset of 105 feature vectors in 8,534 dimensions. This dataset can be processed using the PCA algorithm described in this section. In the following, Sects. 3.1 and 3.2 discuss the dot matrix, its relation to the covariance matrix, and the conversion of eigenvectors obtained from both matrices. The use of the dot matrix is the foundation of the small-sample-size PCA method described in Sect. 3.3.

Covariance matrix and Dot matrix
Given feature vectors 1 , … , n with t ∈ ℝ d . We assume that these vectors are centered around zero. We can construct a data matrix = [ 1 , … , n ] of size d × n . The column i of corresponds to the feature vector i . Using this data matrix, the covariance matrix can be rewritten as We often use T to denote the covariance matrix since maximizing T or T T will result in the same eigenvectors (if we constraint these eigenvectors to be unit vectors).
The dot matrix is defined as T . It is easy to see that the element i, j of this matrix is indeed the dot product between example i and j, i.e. T i j .

Eigenvector conversions
Suppose that and are eigenvalue and corresponding eigenvector of dot matrix, then we have Hence will be an eigenvector of the covariance matrix under the same eigenvalue . This relation allows us to convert the dot matrix's eigenvectors into the covariance matrix's eigenvectors. After the conversion, the newly obtained eigenvector should be normalized as usual. Note that we have Hence the normalized version of the eigenvector is simply (1∕ √ ) On the other hand, it is also possible to convert eigenvectors of the covariance matrix into eigenvectors of the dot matrix. Indeed, if and are eigenvalue and corresponding eigenvector of dot matrix, then we have Hence T will be an eigenvector of the dot matrix under the same eigenvalue . To normalize this eigenvector, we have Hence the normalized version of the eigenvector is (1∕

PCA from dot matrix
From the above discussion, the small-sample-size PCA can be done by eigen-decompose the dot matrix first, then converting the result into eigenvectors of the covariance matrix. The whole procedure is summarized in Algorithm 2. Note that the size of the covariance matrix is d × d , and the size of the dot matrix is n × n . From the computational point of view, the rule of thumb is to eigen-decompose a smaller matrix. If d ≤ n , then we should eigen-decompose the dot matrix, if d > n we eigen-decompose the covariance matrix.

Examples
To illustrate the application of PCA in this case ( n ≤ d ), ORL dataset 6 is considered along with CIFAR-10 dataset. The ORL dataset contains 400 images of 92 × 112 pixels. These images were vectorized into 10,304-dimensional vectors. For CIFAR-10, a dataset of 1,000 examples is constructed for classes 0, 1, and 2. In both cases, the dimension of feature vectors is much higher than the number of examples. Thus, the size of the covariance matrix will be larger than that of the dot matrix. Figure 3 shows results from four datasets namely 'ORL' and 'CIFAR' with class 0, 1, and 2 respectively. The left column shows the eigenvalues plot from the four datasets. The second column shows the cumulative proportion of variance in the principal subspace. The two right columns show a projection of the four datasets onto the two first eigenvectors, and onto the first and the last eigenvectors, respectively. In both cases, the projection onto the first eigenvector corresponds to the abscissa. These 2D plots show the same trend as in the standard PCA case. Indeed, the variance on the horizontal axis is higher than that on the vertical axis since the eigenvectors are sorted in decreasing order of their eigenvalues corresponding to variance on the axis. Moreover, it is interesting to see that most variance is due to the first few principal axes.

Large-scale dataset in high-dimensional space: large n and large d
As time goes by, the size of the face dataset increases from a few hundred images to more than ten thousand images. Consequently, both covariance matrix and dot matrix will be huge and become difficult to eigen-decompose. This section describes how to deal with this problem, starting from a simple heuristic designed for face images (Sect. 4.1), then we discuss why it works (Sect. 4.2) before giving a general method to handle this case in Sect. 4.4.

Motivating idea from face images
Consider the three images in Fig. 4. They were principal axis or Eigenfaces obtained from ORL dataset at 3 resolutions namely 96 × 112 in Fig. 4a, 77 × 89 in Fig. 4b, and 46 × 56 in Fig. 4c. Each face image was resized into the desired resolution; then, their pixels were rearranged to form a feature vector. We can obtain the face image corresponding to each Eigenface by inverting this pixel arrangement.
One can easily see the similarity between the three sets of Eigenfaces. It is worth noting that the negative of an eigenvector is also an eigenvector with the same variance. Thus, given two Eigenfaces, with one that appears to be the inverted color of another, they are, in fact, colinear and can be converted into the same Eigenface. An interesting question is why eigenvectors obtained from different resolutions look so similar.
Recall that an eigenvector i of the dot matrix can be converted into an eigenvector of the covariance matrix by i . Note that Hence, the coordinates of an eigenvector of the dot matrix can be interpreted as relative importance of different examples. This is why Google's Pagerank uses the coordinate of the leading eigenvector of the Internet to describe the importance of each web page (Brin and Page 1998).
In general, resizing does not change the face image's global aspect. Thus, the relative importance between images should be affected slightly. This could explain why resizing produces similar eigenvectors.
Indeed, from the above discussion, the PCA can be done by approximating the eigenstructure, i.e., eigenvalues and eigenvectors, of the original dot matrix by eigenstructure from reduced-size images. The whole PCA process is summarized in Algorithm 3.

Matrix perturbation
From the previous Section, we may wonder what kind of aspect of face images impacts the relative importance, thus the eigenstructure. If we can somehow identify this aspect, maybe a similar approach can also be used for non-image data.
To answer the above question, we first investigate the effect of resizing images. Suppose we scaled up the face images from 46 × 56 back to 92 × 112 . We can see that these images look similar to the original image but with some blocking artifacts. The same effect can be produced by replacing neighbor pixels with the same value, such as the value of the top-left pixel. We call this pixel-grouping version of the original image. Figure 5c shows an image of Fig. 5a with pixel-grouping strategy using 2 × 2 neighborhood. Figure 5b is obtained using downscale + upsample strategy. From Fig. 5, we can see that the downscale + upsample method and the pixel-grouping method produce similar images that are also similar to the original image. Hence, in this subsection, we will consider the feature grouping strategy for the general feature vector.
Indeed, we consider the simple case of replacing one feature with another one; for example, replace feature a by b. Given a vector t , let ̃ tk = tk if k ≠ a and equals to tb if k = a . Following the discussion in the previous subsection, we investigate the dot product of these modified feature vectors: Hence, the dot matrix can be rewritten as with ij = x ia x ja − x ib x jb From matrix perturbation theory, we know that the difference between the eigenvalue of T and that of ̃ T̃ is bounded by the size of . That is, if 1 ≥ 2 ≥ ... ≥ d are the eigenvalues of T and ̃1 ≥ ... ≥̃d the eigenvalues of ̃ T̃ , then we have The first inequality comes from corollary 8.1.6 in Gene (1996) The second inequality is because norm-2 of the matrix is the square root of the maximum eigenvalue, whereas the Frobenius norm is the square root of the trace that is the sum over all eigenvalues. In the case of matrix we have: It is worth noting that ( T ) ab ∕ √ ( T ) aa ( T ) bb is indeed the correlation between features a and b.
From the Eqs. (30) and (35) we can conclude that if a and b are highly correlated, then replacing one by another will only slightly alter the value of obtained eigenvalues. As neighbor pixels are often correlated with each other, the pixels replacing strategy or resizing strategy allows for maintaining dot product, thus the eigenstructure.

Dot product preserving transformation
From the previous analysis, we know that resizing works because it groups features that are highly correlated with each other. Nonetheless, during the analysis, we have seen another critical component: the dot product preserving transformation. Indeed, if we can somehow transform t ∈ ℝ d into ̃ t ∈ ℝ q with q ≪ d such that T i j ≈̃ T ĩ j , then we may follow the same procedure to approximate eigenvalues of the original dataset. The key idea is to preserve dot product and fast calculation. This section reviews some transformations that can be used to this end.

Gaussian random projection
An example of such transformation is the Gaussian random projection. Indeed, given Gaussian random vectors 1 , We know that if q is large enough, then its sample covariance matrix will converge to d , i.e. The dot product between two hash vectors is then Observe that if we take the expectation over hash function s, then the second term on the right-hand side will vanish. As a result h,s [⟨f (h,s)  , but with different conditions on the hash functions. This example underlines that dot product preserving transformations could also be found in other domains that may appear to be distant subjects at first, such as stream processing.

Approximate PCA
In summary, the PCA for large-scale data in high dimensional space can be done by adapting the Algorithm 3 using dot product preserving transformation instead of simple downsizing. This approximate PCA algorithm is summarized in Algorithm 4.

Examples
To illustrate the application of PCA in this case (large n and large d), a larger dataset is considered namely Dogs VS Cats dataset 7 . The training split of this dataset contains 25,000 images of cats and dogs in various sizes. These images were resized into 224 × 224 pixels, then fed into several deep neural networks to extract visual feature vectors. Table 3 summarizes the dimensions of feature vector from different deep neural structures. Figure 6 shows Mean Absolute Proportion Error (MAPE) between the dot product in the input space and that in the space defined by Gaussian random projection (top row) and feature hashing (bottom row). These plots are obtained using q equal to 500, 1000, 2000, 3000, and 5000 dimensions. Figure 6a-d are obtained from VGG16, DenseNet121, InceptionV3, and ResNet50 respectively. Similarly, Fig. 6e-h are obtained from VGG16, DenseNet121, InceptionV3, and ResNet50. This Figure shows that the MAPE reduces as q increases. Moreover, comparing the two dot product preserving transformations, we can observe that Gaussian random projection outperforms feature hashing in several cases. Thus, we shall consider only Gaussian random projection in the following examples. Figure 7 shows the plot of eigenvalues obtained by approximate PCA using Gaussian random projection with 500 dimensions in the top row. This Figure also shows the plot of the cumulative proportion of total variance obtained from different neural structures in the second row. The third row of this Figure shows the projection onto two first principal axis with red color for "Dog" and blue for "Cat". One can observe that even if InceptionV3 has the fastest increase in cumulative variance, its scatter plot shows a mixed area between the two classes. For VGG16 and ResNet50, the two classes seem to be better separated. Nonetheless, it should be emphasized that PCA is not designed for classification tasks. Hence, one should not expect the principal subspace to be suited for classification in all cases. Another subspace projection, such as linear discriminant analysis or LDA (Fukunaga 1990), that is designed for classification, could yield better classification accuracy.

Kernel PCA and its approximation methods
So far, we have considered a simple linear projection of numerical feature vectors. In practice, most data could be cluttered around an area that may not form a linear subspace but a non-linear one. A convenient way of performing non-linear PCA is to use kernel trick (Schölkopf et al. , 2018. The kernel trick idea is to reformulate any linear analysis method to involve only the dot product, then replace the dot product with a kernel function. Et Voilà! You obtain a non-linear version of the original linear method. The kernel PCA relies on this kernel trick to perform non-linear projection. The main crux of kernel PCA lies in the eigen-decomposition of the kernel matrix whose size is n × n where n is the size of the dataset. This section briefly reviews the kernel method, and the Kernel PCA  in Sect. 5.1. Then the approximate KPCA is discussed in Sects. 5.2 and 5.3.

Kernel PCA
The kernel considered in this manuscript is the Mercer's kernel Shawe-Taylor and Cristianini 2004). In short, if is a valid Mercer's kernel, then there exists a mapping function Φ such that Recall that ⟨., .⟩ denotes the dot product. This bracket notation will be used in the following to render formulas easier to read.
The range of the mapping function Φ is called the reproducing kernel Hilbert space or RKHS. To simplify the discussion, we shall refer to it as the kernel space instead of the input space.
The key of the kernel method is that Φ can be unknown meaning that we know it does exist, but we do not know its analytical form. Fortunately, several properties can be derived implicitly from the dot product without accessing the vector's coordinate. For example, the Euclidean distance in kernel space is given by: This implicit calculation allows transforming the classical linear method into a non-linear one at the expense of a higher computational cost.

Implicit centering in kernel space
Indeed as Φ is unknown, we cannot explicitly compute the mean vector nor subtract it from each data point in the kernel space. However, the mean vector can be defined implicitly via Φ as follows: Even if we cannot compute the mean vector explicitly in kernel space, we can compute the dot product between two data points after mean subtraction as follows: In summary, data centering in kernel space can be done implicitly via the modified kernel ̂ . This modified kernel will be used instead of the dot product in dot matrix calculation.
To reduce confusion, we shall refer to the dot matrix in this case as the kernel matrix. The kernel PCA relies on the eigen-decomposition of this kernel matrix : For n data points, the size of this kernel matrix will be n × n.

Projection axis in kernel space
Before writing down the kernel PCA algorithm, it should be noted that the projection axis is also defined in the kernel space, similar to the mean vector above. Hence, it is impossible to write it down explicitly since Φ is unknown. We have to work with its definition. Indeed, by adopting the Eq. (25) to the kernel space with normalization, an eigenvector i of the covariance matrix in kernel space must be of the form: T is the i th eigenvector of the kernel matrix and i its corresponding eigenvalue.
Even if the numerical values of the vector i cannot be computed explicitly, we may still compute the projection of Φ( ) onto i as follows: Given a new data point , the projection of Φ( ) onto i th eigenvector is given by: Note that we may replace in the Eq. (56) by ̃ for projection of centered data in kernel space.

KPCA
The whole kernel PCA can be summarized as shown in Algorithm 5.

Approximate Kernel PCA
The main crux of kernel PCA lies in the eigen-decomposition of the kernel matrix. Indeed, the size of the kernel matrix is n × n , where n is the size of the dataset. As the mapping function, Φ is unknown. It is not possible applying the method described earlier in Sect. 4 directly. In the following, we first present empirical kernel map and its relation to the kernel matrix. Then we combine it with the approximation method described in Sect. 4.

Empirical kernel map and kernel matrix
Given a dataset of n points 1 , ..., n and a kernel , the empirical kernel map (EKM) transforms any data point into n-dimensional space as follows: It is straightforward, seeing that the dot product between two EKM vectors is not equal to the kernel function. Thus EKM does not preserve the dot product in kernel space. Some .
authors have proposed additional methods to restore the kernel value . Our method will not follow this line. Indeed, we consider the outer product used to construct the covariance matrix instead.
Let ̂ t = Φ( t ) be the mapped value of t in the kernel space via the kernel with its mapping function Φ . In fact ̂ t is not computable since Φ is unknown. Let ̂ = [̂ 1 , ...,̂ n ] be the obtained data matrix. The kernel matrix can be expressed as ̂ T̂ , and the EKM can be rewritten by Hence where is the kernel matrix.

Large-scale approximate KPCA algorithm
The last equation in the previous section implies that the covariance matrix of EKM vectors is exactly squared of the kernel matrix. Hence, we may decompose the former instead of the latter. However, one may wonder about this advantage since both matrices are n × n with large n.
The advantage is that we now have access to these EKM vectors' coordinates. Therefore, it is possible to reduce the computational cost using dot product preserving transformation described in Sect. 4. The whole approximate kernel PCA can be summarized in Algorithm 6

Nyström method
It should be noted that there exists another approximate kernel matrix eigen-decomposition called the Nyström method (Williams and Seeger 2000). The basic idea of the Nyström algorithm is to eigen-decompose a small selected working set first. Then to expand the resulting eigenvectors/eigenvalues to cover the whole dataset.
Let 1 , ..., n be given data. Suppose, without loss of generality, that q first points were selected. Let q i and q i , i = 1, … , m be the eigenvalues and the corresponding eigenvectors obtained from q data points. Nyström method approximate the eigenvalues/eigenvectors using kernel matrix as follows: where n,q is the appropriate n × q sub-matrix of the whole kernel matrix .
As the eigen-decomposition is performed on the working set, its selection is crucial for the correctness of the approximated eigenstructure. To this end, several selection methods have been proposed, including uniform sampling, diagonal sampling, and column-norm sampling (Kumar et al. April 2012). For the uniform sampling, each example has an equal probability of being selected. The diagonal sampling selects working point i according to probability p d i such that: The column-norm sampling selects working set with probability p c i such that:

Complexity analysis
For the approximate KPCA, the main computational cost is in EKM calculation ( O(n 2 ) ), the dot product preserving transformation ( O(qnd) ), the eigen-decomposition ( O(q 3 ) ), and eigen-conversion ( O(qn 2 ) ). The highest cost is in eigen-decomposition and eigen-conversion. The total complexity of this method is O((q + 1)n 2 + qnd + q 3 ). For Nyström algorithm, the main computational cost is in working set selection, eigendecomposition ( O(q 3 ) ), and eigen-conversion ( O(q 2 n) ). For working set selection, the complexity varies depending on the selection method. In fact, the complexity of uniform, diagonal, and column-norm sampling are O(1) , O(n) , and O(n 2 ) respectively. The complexity of the Nyström method is smaller than that of Approx-KPCA, especially when uniform sampling is used. It is worth noting that in some cases, this simple uniform sampling could outperform diagonal and column-norm sampling. This is due to the randomized nature of the sampling procedure. Thus, we would suggest using uniform or diagonal sampling first since they have lower complexity costs.

Examples
To illustrate the application of KPCA and its approximation methods, the ORL dataset was considered again. The dimension of the lower subspace in Approx-KPCA is set to the same value as the size of the working set in the Nyström method. Thus Approx-KPCA and Nyström algorithm will need to eigen-decompose matrices of the same size, i.e., q × q . To compare the two approximation methods, it is worth emphasizing that an eigenvector represents the projection axis but not the direction. Indeed, an eigenvector and its negation − will have the same eigenvalue, i.e., the same variance of projected data. Hence, to compare eigenvectors obtained from the different methods, we consider the absolute correlation between them. Let and ′ be two eigenvectors obtained from two different methods, and the similarity measure is defined as follows: A suitable approximation method should have a high similarity, close to 1. Furthermore, we also consider the difference in eigenvalues. Indeed, let and ′ be two eigenvalues obtained from vanilla KPCA and one of the two approximation methods, respectively, The difference between and ′ is measured as: A suitable approximation method should have a smaller difference.
We perform KPCA, Approx-KPCA and KPCA with Nyström using various kernels, namely Tangent hyperbolic kernel with a = 0.0001, r = 0 , RBF kernel with = 0.0001 , Polynomial kernel degree 5, and the sum of these three kernels. The resulting similarities and difference are shown in Figs. 8 and 9. Approx-KPCA (g) and (f) correspond to approximate KPCA using Gaussian random projection and feature hashing, respectively. Nyström (r) refers to Nyström method with random working set selection. Nyström (d) and (c) refer to Nyström method with diagonal sampling and column-norm sampling, respectively The experiment was repeated ten times. Figures 8 and 9 show average similarity between 10 first eigenvectors and average difference between 10 first eigenvalues respectively. From the two Figures, one can observe that both Approx-KPCA and Nyström method yield high similarity, especially for the first eigenvectors. This similarity increases as q increases. Therefore, one can expect that the principal space obtained from these approximation methods lies close to the real one. For eigenvalues, one can see that Approx-KPCA with Gaussian random projection is the best in several cases. Approx-KPCA with feature hashing yields a large difference for tangent hyperbolic and RBF kernels. Nyström method with diagonal and column sampling yield a large difference for polynomial and sum kernels. From these results, we would suggest using Approx-KPCA with Gaussian random projection when precise eigenvalues are needed.
Tanh RBF Poly Sum Fig. 8 Average similarity between first 10 eigenvectors obtained from KPCA and that obtained from difference approximation methods q = 50 q = 200 Tanh RBF Poly Sum Fig. 9 Average difference between first 10 eigenvalues obtained from KPCA and that obtained from difference approximation methods

Summary
In this tutorial, we reviewed the classical PCA method and the problem that may arise when applying it to very small or very large high-dimensional datasets. We have also discussed different methods that may be used to handle these cases. We have shown that a good PCA approximation could be achieved with dot product preserving transformation. This paper considers two particular transformations: Gaussian random projection and the feature hashing method. Other transformations could also be used if they provide dot product preservation and fast calculation. We also showed how this approximation method could be applied to alleviate the computational problem that occurred in kernel PCA. We have made our code in Python that is used in the experiment available at https:// github. com/ peune/ pca. Readers with a background in Python programming can study this code to better understand this subject. All data used in this paper are from public datasets. They are also available upon request.

A Jacobi method for eigen-decomposition
Given a square matrix of size d × d , its eigen-decomposition is given by: The eigen-decomposition of , gives a powerful method to compute power of , e.g. p = p T where p = diag ( p 1 , ..., p d ) . In particular, when p = −1 , this eigen-decomposition tells us that is not invertible if some of its eigenvalues are zeros. It is also possible to compute fancy things like 1∕2 .
There are several eigen-decomposition algorithms. One of the most widely used is the Jacobi method (Jacobi 1846), which relies on rotation matrices. The obtained eigenvectors form an orthonormal basis of the d-dimensional space. Observed that = T can be rearranged as T = . Therefore, eigen-decomposition can be found by searching for a matrix such that that the product T is a diagonal matrix. To this end, the Jacobi method will successively rotate the default basis of the d-dimensional space until the product becomes diagonal. In fact, the d-dimensional space, the default basis is composed of The rotation is performed on a couple of axes. For example given two axis i, j and a rotation angle we define the rotation matrix P(i, j, ): where c = cos and s = sin It can be proved that and T will have the same eigenstructure. Therefore, Jacobi method tries to find P 1 , P 2 , … , P n such that P n ...P 1 AP T 1 ...P T n is diagonal. The main question is how to select (i, j) and . The Jacobi method answers this question by considering � = T As we want ′ to become closer and closer to the diagonal matrix, we should select the largest off-diagonal element ′ ij and choose proper to shrink its value to 0. For this task, we select such that If ii = jj we select = ∕4 . Hence, the Jacobi method iteratively selects (i, j) from largest off-diagonal element, then construct rotation matrix using as defined in 72. The whole process is repeated until all off-diagonal elements become small enough. This requires computational cost around O(d 3 ).

B Kernel function
This section introduces the basic idea behind kernel and kernel methods (Schölkopf et al. , 2018Shawe-Taylor and Cristianini 2004;Hofmann et al. 2008). The kernel that we are interested in is the Mercer's kernel. A real-valued function satisfies the Mercer's condition if for all square-integrable function g we have ∫ ∫ g( ) ( , )g( ) ≥ 0 . In practice, given a dataset 1 , … , and its kernel matrix , this condition implies 0 =(c 2 − s 2 ) ij + sc( ii − jj ) = cos(2 ) ij + 1 2 sin(2 )( ii − jj ) that the kernel matrix is positive-semidefinite matrix, i.e. for all vector ∈ ℝ n we have T ≥ 0 The Mercer condition guarantees that there exists a mapping function Φ such that Recall that ⟨., .⟩ denotes the dot product. This bracket notation will be used in the following to render formulas easier to read.
The range of the mapping function Φ is called the reproducing kernel Hilbert space or RKHS. To simplify the discussion, we shall refer to it as the kernel space as opposed to the input space.
The key of the kernel method is that Φ can be unknown meaning that we know it does exist, but we do not know how it is. Fortunately, several properties can be derived directly from the dot product without accessing the vector's coordinate. Thus, these properties can also be computed from kernel function, even if Φ remains unknown. An example of such properties is the Euclidean distance. Indeed, as we have so by replacing the dot product with a kernel function we obtain: where Φ is the mapping function of the kernel .
-String kernel (Lodhi et al. 2002) between two strings , defined over finite alphabet Σ : where n is the given length of sub-string considered. For a string = 1 , … , | | with i ∈ Σ , let [i ∶ j] = i , ..., j . The mapping function ( ) is defined as the weighted sum of sub-string found within , i.e.
with < 1 weighting coefficient. -Fisher kernel (Jaakkola and Haussler 1999): for a given probabilistic model p: where the explicit mapping function Φ is the gradient of log likelihood: and is the Fisher's information matrix, i.e.
In practice, we often use an approximate version of this kernel, i.e. ( , ) = Φ( ) T Φ( ). This kernel allows nicely integrating a generative model into a discriminative framework such as SVM. It is also possible constructing a new kernel from existing kernels. Let 1 , 2 ∶ X × X → ℝ be two valid kernels. The following kernels are also valid -( , ) = 1 ( , ) + 2 ( , ) -( , ) = c 1 ( , ), c ∈ ℝ + -( , ) = 1 ( , ) + c, c ∈ ℝ + -( , ) = 1 ( , ) 2 ( , ) -( , ) = f ( )f ( ) for all f ∶ X → ℝ -( , ) = exp − ( 1 ( , ) + 1 ( , ) − 2 1 ( , )) (note: this is an RBF kernel in another kernel space) -( , ) = Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.