-
@@ -227,7 +242,7 @@ MathJax.Hub.Config({
-We assume now that we have a design matrix \( \boldsymbol{X} \) which has been centered as discussed above. For the sake of simplicity we skip the overline symbol. The matrix is defined in terms of the various column vectors \( [\boldsymbol{x}_0,\boldsymbol{x}_1,\dots, \boldsymbol{x}_{p-1}] \) -each with dimension \( \boldsymbol{x}\in {\mathbb{R}}^{n} \). - -
-We assume also that we have an orthogonal transformation \( \boldsymbol{W}\in {\mathbb{R}}^{p\times p} \). We define the reconstruction error (which is similar to the mean squared error we have seen before) as +We will use a simple example first with two-dimensional data +drawn from a multivariate normal distribution with the following mean and covariance matrix: $$ -J(\boldsymbol{W},\boldsymbol{Z}) = \frac{1}{n}\sum_i (\boldsymbol{x}_i - \overline{\boldsymbol{x}}_i)^2, +\mu = (-1,2) \qquad \Sigma = \begin{bmatrix} 4 & 2 \\ +2 & 2 +\end{bmatrix} $$ -with \( \overline{\boldsymbol{x}}_i = \boldsymbol{W}\boldsymbol{z}_i \), where \( \boldsymbol{z}_i \) is a row vector with dimension \( {\mathbb{R}}^{n} \) of the matrix -\( \boldsymbol{Z}\in{\mathbb{R}}^{p\times n} \). When doing PCA we want to reduce this dimensionality. +
+We will generate \( N = 1000 \) points \( X = \{ x_1, \ldots, x_N \} \) from +this distribution, and store them in the \( 1000 \times 2 \) matrix \( \boldsymbol{X} \).
-The PCA theorem states that minimizing the above reconstruction error corresponds to setting \( \boldsymbol{W}=\boldsymbol{S} \), the orthogonal matrix which diagonalizes the empirical covariance(correlation) matrix. The optimal low-dimensional encoding of the data is then given by a set of vectors \( \boldsymbol{z}_i \) with at most \( l \) vectors, with \( l < < p \), defined by the orthogonal projection of the data onto the columns spanned by the eigenvectors of the covariance(correlations matrix). +The following Python code aids in setting up the data + +
+ + +
N = 1000
+mean = (-1, 2)
+cov = [[4, 2], [2, 2]]
+X = np.random.multivariate_normal(mean, cov, N)
++Make a small Python code which plots the data. + +
+Now we are going to implement the PCA algorithm. We will break it down into sub-steps and across multiple cells. + +
+The first step of PCA is to compute the sample mean of the data and use it to center the data. Recall the sample mean is +$$ +\mu_N = \frac{1}{N} \sum_{i=1}^N x_i +$$ + +and the mean-centered data \( \bar{X} = \{ \bar{x}_1, \ldots, \bar{x}_N \} \) takes the form +$$ +\bar{x}_i = x_i - \mu_N +$$ + +When you are done with these steps, print out \( \mu_N \) to verify it is +close to \( \mu \) and plot your mean centered data to verify it is +centered at the origin! Compare your code with the functionality from Scikit-Learn discussed above. + +
+Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by: +$$ +\begin{equation*} +\Sigma_N = \frac{1}{N-1} \sum_{i=1}^N \bar{x}_i^T \bar{x}_i = \frac{1}{N-1} \sum_{i=1}^N (x_i - \mu_N)^T (x_i - \mu_N) +\end{equation*} +$$ + +where the data points \( x_i \in \mathbb{R}^p \) (here in this example \( p = 2 \)) are column vectors and \( x^T \) is the transpose of \( x \). +Compare the computed covariance with the answer given above. + +
+Now we are ready to solve for the principal components! To do so we +diagonalize the sample covariance matrix \( \Sigma_N \). We can use the +function np.linalg.eig to do so. It will return the eigenvalues and +eigenvectors of \( \Sigma_N \). Once you have these, carry out the +following tasks: + +
+Finally, collect all these steps and write your own PCA function and +compare this with the functionality included in Scikit-Learn. +Have the input be the data and have the output be the principal components and their associated eigenvalues, sorted in descending order. Can you think of a way to make it more efficient than the algorithm outlined above? + +
+Finally, try out your own PCA function with other data sets. + +
+After this we ask ourselves how do we prove the link between the maximum variance and the feature reduction.
@@ -228,7 +322,7 @@ The PCA theorem states that minimizing the above reconstruction error correspond
-To show the PCA theorem let us start with the assumption that there is one vector \( \boldsymbol{w}_0 \) which corresponds to a solution which minimized the reconstruction error \( J \). This is an orthogonal vector. It means that we now approximate the reconstruction error in terms of \( \boldsymbol{w}_0 \) and \( \boldsymbol{z}_0 \) as +We assume now that we have a design matrix \( \boldsymbol{X} \) which has been centered as discussed above. For the sake of simplicity we skip the overline symbol. The matrix is defined in terms of the various column vectors \( [\boldsymbol{x}_0,\boldsymbol{x}_1,\dots, \boldsymbol{x}_{p-1}] \) +each with dimension \( \boldsymbol{x}\in {\mathbb{R}}^{n} \). + +
+We assume also that we have an orthogonal transformation \( \boldsymbol{W}\in {\mathbb{R}}^{p\times p} \). We define the reconstruction error (which is similar to the mean squared error we have seen before) as $$ -J(\boldsymbol{w}_0,\boldsymbol{z}_0)= \frac{1}{n}\sum_i (\boldsymbol{x}_i - z_{i0}\boldsymbol{w}_0)^2=\frac{1}{n}\sum_i (\boldsymbol{x}_i^T\boldsymbol{x}_i - 2z_{i0}\boldsymbol{w}_0^T\boldsymbol{x}_i+z_{i0}^2\boldsymbol{w}_0^T\boldsymbol{w}_0), +J(\boldsymbol{W},\boldsymbol{Z}) = \frac{1}{n}\sum_i (\boldsymbol{x}_i - \overline{\boldsymbol{x}}_i)^2, $$ -which we can rewrite due to the orthogonality of \( \boldsymbol{w}_i \) as -$$ -J(\boldsymbol{w}_0,\boldsymbol{z}_0)=\frac{1}{n}\sum_i (\boldsymbol{x}_i^T\boldsymbol{x}_i - 2z_{i0}\boldsymbol{w}_0^T\boldsymbol{x}_i+z_{i0}^2). -$$ +with \( \overline{\boldsymbol{x}}_i = \boldsymbol{W}\boldsymbol{z}_i \), where \( \boldsymbol{z}_i \) is a row vector with dimension \( {\mathbb{R}}^{n} \) of the matrix +\( \boldsymbol{Z}\in{\mathbb{R}}^{p\times n} \). When doing PCA we want to reduce this dimensionality. -Minimizing \( J \) with respect to the unknown parameters \( z_{0i} \) we obtain that -$$ -z_{i0}=\boldsymbol{w}_0^T\boldsymbol{x}_i, -$$ - -where the vectors on the rhs are known. +
+The PCA theorem states that minimizing the above reconstruction error corresponds to setting \( \boldsymbol{W}=\boldsymbol{S} \), the orthogonal matrix which diagonalizes the empirical covariance(correlation) matrix. The optimal low-dimensional encoding of the data is then given by a set of vectors \( \boldsymbol{z}_i \) with at most \( l \) vectors, with \( l < < p \), defined by the orthogonal projection of the data onto the columns spanned by the eigenvectors of the covariance(correlations matrix).
@@ -230,7 +243,7 @@ where the vectors on the rhs are known.
-We have now found the unknown parameters \( z_{i0} \). These correspond to the projected coordinates and we can write +To show the PCA theorem let us start with the assumption that there is one vector \( \boldsymbol{w}_0 \) which corresponds to a solution which minimized the reconstruction error \( J \). This is an orthogonal vector. It means that we now approximate the reconstruction error in terms of \( \boldsymbol{w}_0 \) and \( \boldsymbol{z}_0 \) as $$ -J(\boldsymbol{w}_0)= \frac{1}{p}\sum_i (\boldsymbol{x}_i^T\boldsymbol{x}_i - z_{i0}^2)=\mathrm{const}-\frac{1}{n}\sum_i z_{i0}^2. +J(\boldsymbol{w}_0,\boldsymbol{z}_0)= \frac{1}{n}\sum_i (\boldsymbol{x}_i - z_{i0}\boldsymbol{w}_0)^2=\frac{1}{n}\sum_i (\boldsymbol{x}_i^T\boldsymbol{x}_i - 2z_{i0}\boldsymbol{w}_0^T\boldsymbol{x}_i+z_{i0}^2\boldsymbol{w}_0^T\boldsymbol{w}_0), $$ -
-We can show that the variance of the projected coordinates defined by \( \boldsymbol{w}_0^T\boldsymbol{x}_i \) are given by +which we can rewrite due to the orthogonality of \( \boldsymbol{w}_i \) as $$ -\mathrm{var}[\boldsymbol{w}_0^T\boldsymbol{x}_i] = \frac{1}{n}\sum_i z_{i0}^2, +J(\boldsymbol{w}_0,\boldsymbol{z}_0)=\frac{1}{n}\sum_i (\boldsymbol{x}_i^T\boldsymbol{x}_i - 2z_{i0}\boldsymbol{w}_0^T\boldsymbol{x}_i+z_{i0}^2). $$ -since the expectation value of +Minimizing \( J \) with respect to the unknown parameters \( z_{0i} \) we obtain that $$ -\mathbb{E}[\boldsymbol{w}_0^T\boldsymbol{x}_i] = \mathbb{E}[z_{i0}]= \boldsymbol{w}_0^T\mathbb{E}[\boldsymbol{x}_i]=0, +z_{i0}=\boldsymbol{w}_0^T\boldsymbol{x}_i, $$ -where we have used the fact that our data are centered. - -
-Recalling our definition of the covariance as -$$ -\boldsymbol{C}[\boldsymbol{x}] = \frac{1}{n}\boldsymbol{X}\boldsymbol{X}^T=\mathbb{E}[\boldsymbol{X}\boldsymbol{X}^T], -$$ - -we have thus that -$$ -\mathrm{var}[\boldsymbol{w}_0^T\boldsymbol{x}_i] = \frac{1}{n}\sum_i z_{i0}^2=\boldsymbol{w}_0^T\boldsymbol{C}[\boldsymbol{x}]\boldsymbol{w}_0. -$$ - -
-We are almost there, we have obtained a relation between minimizing the reconstruction error and the variance and the covariance matrix. Minimizing the error is equivalent to maximizing the variance of the projected data. +where the vectors on the rhs are known.
@@ -245,7 +245,7 @@ We are almost there, we have obtained a relation between minimizing the reconstr
-We could trivially maximize the variance of the projection (and -thereby minimize the error in the reconstruction function) by letting -the norm-2 of \( \boldsymbol{w}_0 \) go to infinity. However, this norm since we -want the matrix \( \boldsymbol{W} \) to be an orthogonal matrix, is constrained by -\( \vert\vert \boldsymbol{w}_0 \vert\vert_2^2=1 \). Imposing this condition via a -Lagrange multiplier we can then in turn maximize - +We have now found the unknown parameters \( z_{i0} \). These correspond to the projected coordinates and we can write $$ -J(\boldsymbol{w}_0)= \boldsymbol{w}_0^T\boldsymbol{C}[\boldsymbol{x}]\boldsymbol{w}_0+\lambda_0(1-\boldsymbol{w}_0^T\boldsymbol{w}_0). -$$ - -Taking the derivative with respect to \( \boldsymbol{w}_0 \) we obtain - -$$ -\frac{\partial J(\boldsymbol{w}_0)}{\partial \boldsymbol{w}_0}= 2\boldsymbol{C}[\boldsymbol{x}]\boldsymbol{w}_0-2\lambda_0\boldsymbol{w}_0=0, -$$ - -meaning that -$$ -\boldsymbol{C}[\boldsymbol{x}]\boldsymbol{w}_0=\lambda_0\boldsymbol{w}_0. -$$ - -The direction that maximizes the variance (or minimizes the construction error) is an eigenvector of the covariance matrix! If we left multiply with \( \boldsymbol{w}_0^T \) we have the variance of the projected data is -$$ -\boldsymbol{w}_0^T\boldsymbol{C}[\boldsymbol{x}]\boldsymbol{w}_0=\lambda_0. +J(\boldsymbol{w}_0)= \frac{1}{p}\sum_i (\boldsymbol{x}_i^T\boldsymbol{x}_i - z_{i0}^2)=\mathrm{const}-\frac{1}{n}\sum_i z_{i0}^2. $$
-If we want to maximize the variance (minimize the construction error) -we simply pick the eigenvector of the covariance matrix with the -largest eigenvalue. This establishes the link between the minimization -of the reconstruction function \( J \) in terms of an orthogonal matrix -and the maximization of the variance and thereby the covariance of our -observations encoded in the design/feature matrix \( \boldsymbol{X} \). +We can show that the variance of the projected coordinates defined by \( \boldsymbol{w}_0^T\boldsymbol{x}_i \) are given by +$$ +\mathrm{var}[\boldsymbol{w}_0^T\boldsymbol{x}_i] = \frac{1}{n}\sum_i z_{i0}^2, +$$ + +since the expectation value of +$$ +\mathbb{E}[\boldsymbol{w}_0^T\boldsymbol{x}_i] = \mathbb{E}[z_{i0}]= \boldsymbol{w}_0^T\mathbb{E}[\boldsymbol{x}_i]=0, +$$ + +where we have used the fact that our data are centered.
-The proof -for the other eigenvectors \( \boldsymbol{w}_1,\boldsymbol{w}_2,\dots \) can be -established by applying the above arguments and using the fact that -our basis of eigenvectors is orthogonal, see Murphy chapter -12.2. The -discussion in chapter 12.2 of Murphy's text has also a nice link with -the Singular Value Decomposition theorem. For categorical data, see -chapter 12.4 and discussion therein. +Recalling our definition of the covariance as +$$ +\boldsymbol{C}[\boldsymbol{x}] = \frac{1}{n}\boldsymbol{X}\boldsymbol{X}^T=\mathbb{E}[\boldsymbol{X}\boldsymbol{X}^T], +$$ + +we have thus that +$$ +\mathrm{var}[\boldsymbol{w}_0^T\boldsymbol{x}_i] = \frac{1}{n}\sum_i z_{i0}^2=\boldsymbol{w}_0^T\boldsymbol{C}[\boldsymbol{x}]\boldsymbol{w}_0. +$$ + +
+We are almost there, we have obtained a relation between minimizing the reconstruction error and the variance and the covariance matrix. Minimizing the error is equivalent to maximizing the variance of the projected data.
@@ -257,6 +259,8 @@ chapter 12.4 and discussion therein.
-Principal Component Analysis (PCA) is by far the most popular dimensionality reduction algorithm. -First it identifies the hyperplane that lies closest to the data, and then it projects the data onto it. +
-The following Python code uses NumPy’s svd() function to obtain all the principal components of the -training set, then extracts the first two principal components. First we center the data using either pandas or our own code -
+We could trivially maximize the variance of the projection (and +thereby minimize the error in the reconstruction function) by letting +the norm-2 of \( \boldsymbol{w}_0 \) go to infinity. However, this norm since we +want the matrix \( \boldsymbol{W} \) to be an orthogonal matrix, is constrained by +\( \vert\vert \boldsymbol{w}_0 \vert\vert_2^2=1 \). Imposing this condition via a +Lagrange multiplier we can then in turn maximize - -
import numpy as np
-import pandas as pd
-from IPython.display import display
-np.random.seed(100)
-# setting up a 10 x 5 vanilla matrix
-rows = 10
-cols = 5
-X = np.random.randn(rows,cols)
-df = pd.DataFrame(X)
-# Pandas does the centering for us
-df = df -df.mean()
-display(df)
+$$
+J(\boldsymbol{w}_0)= \boldsymbol{w}_0^T\boldsymbol{C}[\boldsymbol{x}]\boldsymbol{w}_0+\lambda_0(1-\boldsymbol{w}_0^T\boldsymbol{w}_0).
+$$
-# we center it ourselves
-X_centered = X - X.mean(axis=0)
-# Then check the difference between pandas and our own set up
-print(X_centered-df)
-#Now we do an SVD
-U, s, V = np.linalg.svd(X_centered)
-c1 = V.T[:, 0]
-c2 = V.T[:, 1]
-W2 = V.T[:, :2]
-X2D = X_centered.dot(W2)
-print(X2D)
--PCA assumes that the dataset is centered around the origin. Scikit-Learn’s PCA classes take care of centering -the data for you. However, if you implement PCA yourself (as in the preceding example), or if you use other libraries, don’t -forget to center the data first. +Taking the derivative with respect to \( \boldsymbol{w}_0 \) we obtain + +$$ +\frac{\partial J(\boldsymbol{w}_0)}{\partial \boldsymbol{w}_0}= 2\boldsymbol{C}[\boldsymbol{x}]\boldsymbol{w}_0-2\lambda_0\boldsymbol{w}_0=0, +$$ + +meaning that +$$ +\boldsymbol{C}[\boldsymbol{x}]\boldsymbol{w}_0=\lambda_0\boldsymbol{w}_0. +$$ + +The direction that maximizes the variance (or minimizes the construction error) is an eigenvector of the covariance matrix! If we left multiply with \( \boldsymbol{w}_0^T \) we have the variance of the projected data is +$$ +\boldsymbol{w}_0^T\boldsymbol{C}[\boldsymbol{x}]\boldsymbol{w}_0=\lambda_0. +$$
-Once you have identified all the principal components, you can reduce the dimensionality of the dataset -down to \( d \) dimensions by projecting it onto the hyperplane defined by the first \( d \) principal components. -Selecting this hyperplane ensures that the projection will preserve as much variance as possible. -
+If we want to maximize the variance (minimize the construction error) +we simply pick the eigenvector of the covariance matrix with the +largest eigenvalue. This establishes the link between the minimization +of the reconstruction function \( J \) in terms of an orthogonal matrix +and the maximization of the variance and thereby the covariance of our +observations encoded in the design/feature matrix \( \boldsymbol{X} \). + +
+The proof +for the other eigenvectors \( \boldsymbol{w}_1,\boldsymbol{w}_2,\dots \) can be +established by applying the above arguments and using the fact that +our basis of eigenvectors is orthogonal, see Murphy chapter +12.2. The +discussion in chapter 12.2 of Murphy's text has also a nice link with +the Singular Value Decomposition theorem. For categorical data, see +chapter 12.4 and discussion therein. - -
W2 = V.T[:, :2]
-X2D = X_centered.dot(W2)
-
@@ -261,6 +271,7 @@ X2D = X_centered29
+Principal Component Analysis (PCA) is by far the most popular dimensionality reduction algorithm.
+First it identifies the hyperplane that lies closest to the data, and then it projects the data onto it.
-Scikit-Learn’s PCA class implements PCA using SVD decomposition just like we did before. The
-following code applies PCA to reduce the dimensionality of the dataset down to two dimensions (note
-that it automatically takes care of centering the data):
+The following Python code uses NumPy’s svd() function to obtain all the principal components of the
+training set, then extracts the first two principal components. First we center the data using either pandas or our own code
-
-After fitting the PCA transformer to the dataset, you can access the principal components using the
-components variable (note that it contains the PCs as horizontal vectors, so, for example, the first
-principal component is equal to
+PCA assumes that the dataset is centered around the origin. Scikit-Learn’s PCA classes take care of centering
+the data for you. However, if you implement PCA yourself (as in the preceding example), or if you use other libraries, don’t
+forget to center the data first.
+
+
+Once you have identified all the principal components, you can reduce the dimensionality of the dataset
+down to \( d \) dimensions by projecting it onto the hyperplane defined by the first \( d \) principal components.
+Selecting this hyperplane ensures that the projection will preserve as much variance as possible.
-
-Another very useful piece of information is the explained variance ratio of each principal component,
-available via the \( explained\_variance\_ratio \) variable. It indicates the proportion of the dataset’s
-variance that lies along the axis of each principal component.
-
@@ -236,6 +275,7 @@ variance that lies along the axis of each principal component.
+Scikit-Learn’s PCA class implements PCA using SVD decomposition just like we did before. The
+following code applies PCA to reduce the dimensionality of the dataset down to two dimensions (note
+that it automatically takes care of centering the data):
-
-We see that our training data after the PCA decomposition has a performance similar to the non-scaled data.
+After fitting the PCA transformer to the dataset, you can access the principal components using the
+components variable (note that it contains the PCs as horizontal vectors, so, for example, the first
+principal component is equal to
+
+
+
+
+Another very useful piece of information is the explained variance ratio of each principal component,
+available via the \( explained\_variance\_ratio \) variable. It indicates the proportion of the dataset’s
+variance that lies along the axis of each principal component.
@@ -244,6 +250,7 @@ We see that our training data after the PCA decomposition has a performance simi
-Instead of arbitrarily choosing the number of dimensions to reduce down to, it is generally preferable to
-choose the number of dimensions that add up to a sufficiently large portion of the variance (e.g., 95%).
-Unless, of course, you are reducing dimensionality for data visualization — in that case you will
-generally want to reduce the dimensionality down to 2 or 3.
-The following code computes PCA without reducing dimensionality, then computes the minimum number
-of dimensions required to preserve 95% of the training set’s variance:
+
-
-You could then set \( n\_components=d \) and run PCA again. However, there is a much better option: instead
-of specifying the number of principal components you want to preserve, you can set \( n\_components \) to be
-a float between 0.0 and 1.0, indicating the ratio of variance you wish to preserve:
-
+
+We see that our training data after the PCA decomposition has a performance similar to the non-scaled data.
+
@@ -232,6 +258,7 @@ X_reduced = pca
-One problem with the preceding implementation of PCA is that it requires the whole training set to fit in
-memory in order for the SVD algorithm to run. Fortunately, Incremental PCA (IPCA) algorithms have
-been developed: you can split the training set into mini-batches and feed an IPCA algorithm one minibatch
-at a time. This is useful for large training sets, and also to apply PCA online (i.e., on the fly, as new
-instances arrive).
+Instead of arbitrarily choosing the number of dimensions to reduce down to, it is generally preferable to
+choose the number of dimensions that add up to a sufficiently large portion of the variance (e.g., 95%).
+Unless, of course, you are reducing dimensionality for data visualization — in that case you will
+generally want to reduce the dimensionality down to 2 or 3.
+The following code computes PCA without reducing dimensionality, then computes the minimum number
+of dimensions required to preserve 95% of the training set’s variance:
+
+
+
+You could then set \( n\_components=d \) and run PCA again. However, there is a much better option: instead
+of specifying the number of principal components you want to preserve, you can set \( n\_components \) to be
+a float between 0.0 and 1.0, indicating the ratio of variance you wish to preserve:
+
+
+
+
@@ -213,6 +246,7 @@ instances arrive).
-Scikit-Learn offers yet another option to perform PCA, called Randomized PCA. This is a stochastic
-algorithm that quickly finds an approximation of the first d principal components. Its computational
-complexity is \( O(m \times d^2)+O(d^3) \), instead of \( O(m \times n^2) + O(n^3) \), so it is dramatically faster than the
-previous algorithms when \( d \) is much smaller than \( n \).
-
-
-
@@ -216,6 +227,7 @@ previous algorithms when \( d \) is much smaller than \( n \).
+
-The kernel trick is a mathematical technique that implicitly maps instances into a
-very high-dimensional space (called the feature space), enabling nonlinear classification and regression
-with Support Vector Machines. Recall that a linear decision boundary in the high-dimensional feature
-space corresponds to a complex nonlinear decision boundary in the original space.
-It turns out that the same trick can be applied to PCA, making it possible to perform complex nonlinear
-projections for dimensionality reduction. This is called Kernel PCA (kPCA). It is often good at
-preserving clusters of instances after projection, or sometimes even unrolling datasets that lie close to a
-twisted manifold.
-For example, the following code uses Scikit-Learn’s KernelPCA class to perform kPCA with an
-
+Scikit-Learn offers yet another option to perform PCA, called Randomized PCA. This is a stochastic
+algorithm that quickly finds an approximation of the first d principal components. Its computational
+complexity is \( O(m \times d^2)+O(d^3) \), instead of \( O(m \times n^2) + O(n^3) \), so it is dramatically faster than the
+previous algorithms when \( d \) is much smaller than \( n \).
-
-
-Locally Linear Embedding (LLE) is another very powerful nonlinear dimensionality reduction
-(NLDR) technique. It is a Manifold Learning technique that does not rely on projections like the previous
-algorithms. In a nutshell, LLE works by first measuring how each training instance linearly relates to its
-closest neighbors (c.n.), and then looking for a low-dimensional representation of the training set where
-these local relationships are best preserved (more details shortly).
+The kernel trick is a mathematical technique that implicitly maps instances into a
+very high-dimensional space (called the feature space), enabling nonlinear classification and regression
+with Support Vector Machines. Recall that a linear decision boundary in the high-dimensional feature
+space corresponds to a complex nonlinear decision boundary in the original space.
+It turns out that the same trick can be applied to PCA, making it possible to perform complex nonlinear
+projections for dimensionality reduction. This is called Kernel PCA (kPCA). It is often good at
+preserving clusters of instances after projection, or sometimes even unrolling datasets that lie close to a
+twisted manifold.
+For example, the following code uses Scikit-Learn’s KernelPCA class to perform kPCA with an
+
+
+
+
+
@@ -210,6 +243,7 @@ these local relationships are best preserved (more details shortly).
-There are many other dimensionality reduction techniques, several of which are available in Scikit-Learn.
+Locally Linear Embedding (LLE) is another very powerful nonlinear dimensionality reduction
+(NLDR) technique. It is a Manifold Learning technique that does not rely on projections like the previous
+algorithms. In a nutshell, LLE works by first measuring how each training instance linearly relates to its
+closest neighbors (c.n.), and then looking for a low-dimensional representation of the training set where
+these local relationships are best preserved (more details shortly).
-Here are some of the most popular:
-
-
-
@@ -227,7 +242,7 @@ MathJax.Hub.Config({
@@ -982,14 +982,120 @@ $$
+We will use a simple example first with two-dimensional data
+drawn from a multivariate normal distribution with the following mean and covariance matrix:
+
+We will generate \( N = 1000 \) points \( X = \{ x_1, \ldots, x_N \} \) from
+this distribution, and store them in the \( 1000 \times 2 \) matrix \( \boldsymbol{X} \).
+
+
+The following Python code aids in setting up the data
+
+
+
+Make a small Python code which plots the data.
+
+
+Now we are going to implement the PCA algorithm. We will break it down into sub-steps and across multiple cells.
+
+
+The first step of PCA is to compute the sample mean of the data and use it to center the data. Recall the sample mean is
+
+Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by:
+
+Now we are ready to solve for the principal components! To do so we
+diagonalize the sample covariance matrix \( \Sigma_N \). We can use the
+function np.linalg.eig to do so. It will return the eigenvalues and
+eigenvectors of \( \Sigma_N \). Once you have these, carry out the
+following tasks:
+
+
+Finally, collect all these steps and write your own PCA function and
+compare this with the functionality included in Scikit-Learn.
+Have the input be the data and have the output be the principal components and their associated eigenvalues, sorted in descending order. Can you think of a way to make it more efficient than the algorithm outlined above?
+
+
+Finally, try out your own PCA function with other data sets.
+
+
After this we ask ourselves how do we prove the link between the maximum variance and the feature reduction.
We assume now that we have a design matrix \( \boldsymbol{X} \) which has been centered as discussed above. For the sake of simplicity we skip the overline symbol. The matrix is defined in terms of the various column vectors \( [\boldsymbol{x}_0,\boldsymbol{x}_1,\dots, \boldsymbol{x}_{p-1}] \)
@@ -1012,7 +1118,7 @@ The PCA theorem states that minimizing the above reconstruction error correspond
To show the PCA theorem let us start with the assumption that there is one vector \( \boldsymbol{w}_0 \) which corresponds to a solution which minimized the reconstruction error \( J \). This is an orthogonal vector. It means that we now approximate the reconstruction error in terms of \( \boldsymbol{w}_0 \) and \( \boldsymbol{z}_0 \) as
@@ -1041,7 +1147,7 @@ where the vectors on the rhs are known.
We have now found the unknown parameters \( z_{i0} \). These correspond to the projected coordinates and we can write
@@ -1089,7 +1195,7 @@ We are almost there, we have obtained a relation between minimizing the reconstr
We could trivially maximize the variance of the projection (and
@@ -1148,7 +1254,7 @@ chapter 12.4 and discussion therein.
@@ -1205,7 +1311,7 @@ X2D = X_centered.dot(W2)
Scikit-Learn’s PCA class implements PCA using SVD decomposition just like we did before. The
@@ -1237,7 +1343,7 @@ variance that lies along the axis of each principal component.
@@ -1278,7 +1384,7 @@ We see that our training data after the PCA decomposition has a performance simi
Instead of arbitrarily choosing the number of dimensions to reduce down to, it is generally preferable to
@@ -1309,7 +1415,7 @@ X_reduced = pca.fit_transform(X)
One problem with the preceding implementation of PCA is that it requires the whole training set to fit in
@@ -1321,7 +1427,7 @@ instances arrive).
Scikit-Learn offers yet another option to perform PCA, called Randomized PCA. This is a stochastic
@@ -1335,7 +1441,7 @@ previous algorithms when \( d \) is much smaller than \( n \).
@@ -1361,7 +1467,7 @@ X_reduced = rbf_pca.fit_transform(X)
Locally Linear Embedding (LLE) is another very powerful nonlinear dimensionality reduction
@@ -1373,7 +1479,7 @@ these local relationships are best preserved (more details shortly).
There are many other dimensionality reduction techniques, several of which are available in Scikit-Learn.
diff --git a/doc/pub/DimRed/html/DimRed-solarized.html b/doc/pub/DimRed/html/DimRed-solarized.html
index 1708ad00f..686a4c32d 100644
--- a/doc/pub/DimRed/html/DimRed-solarized.html
+++ b/doc/pub/DimRed/html/DimRed-solarized.html
@@ -109,19 +109,30 @@ div { text-align: justify; text-justify: inter-word; }
'___sec14'),
('Towards the PCA theorem', 2, None, '___sec15'),
('The Algorithm before the Theorem', 2, None, '___sec16'),
- ('Classical PCA Theorem', 2, None, '___sec17'),
- ('Proof of the PCA Theorem', 2, None, '___sec18'),
- ('PCA Proof continued', 2, None, '___sec19'),
- ('The final step', 2, None, '___sec20'),
- ('Principal Component Analysis', 2, None, '___sec21'),
- ('PCA and scikit-learn', 2, None, '___sec22'),
- ('Back to the Cancer Data', 2, None, '___sec23'),
- ('More on the PCA', 2, None, '___sec24'),
- ('Incremental PCA', 2, None, '___sec25'),
- ('Randomized PCA', 2, None, '___sec26'),
- ('Kernel PCA', 2, None, '___sec27'),
- ('LLE', 2, None, '___sec28'),
- ('Other techniques', 2, None, '___sec29')]}
+ ('Writing our own PCA code', 2, None, '___sec17'),
+ ('Compute the sample mean and center the data',
+ 3,
+ None,
+ '___sec18'),
+ ('Compute the sample covariance', 3, None, '___sec19'),
+ ('Diagonalize the sample covariance matrix to obtain the '
+ 'principal components',
+ 3,
+ None,
+ '___sec20'),
+ ('Classical PCA Theorem', 2, None, '___sec21'),
+ ('Proof of the PCA Theorem', 2, None, '___sec22'),
+ ('PCA Proof continued', 2, None, '___sec23'),
+ ('The final step', 2, None, '___sec24'),
+ ('Principal Component Analysis', 2, None, '___sec25'),
+ ('PCA and scikit-learn', 2, None, '___sec26'),
+ ('Back to the Cancer Data', 2, None, '___sec27'),
+ ('More on the PCA', 2, None, '___sec28'),
+ ('Incremental PCA', 2, None, '___sec29'),
+ ('Randomized PCA', 2, None, '___sec30'),
+ ('Kernel PCA', 2, None, '___sec31'),
+ ('LLE', 2, None, '___sec32'),
+ ('Other techniques', 2, None, '___sec33')]}
end of tocinfo -->
-
+We will use a simple example first with two-dimensional data
+drawn from a multivariate normal distribution with the following mean and covariance matrix:
+$$
+\mu = (-1,2) \qquad \Sigma = \begin{bmatrix} 4 & 2 \\
+2 & 2
+\end{bmatrix}
+$$
+
+
+We will generate \( N = 1000 \) points \( X = \{ x_1, \ldots, x_N \} \) from
+this distribution, and store them in the \( 1000 \times 2 \) matrix \( \boldsymbol{X} \).
+
+
+The following Python code aids in setting up the data
+
+
+
+
+
+Make a small Python code which plots the data.
+
+
+Now we are going to implement the PCA algorithm. We will break it down into sub-steps and across multiple cells.
+
+
+The first step of PCA is to compute the sample mean of the data and use it to center the data. Recall the sample mean is
+$$
+\mu_N = \frac{1}{N} \sum_{i=1}^N x_i
+$$
+
+and the mean-centered data \( \bar{X} = \{ \bar{x}_1, \ldots, \bar{x}_N \} \) takes the form
+$$
+\bar{x}_i = x_i - \mu_N
+$$
+
+When you are done with these steps, print out \( \mu_N \) to verify it is
+close to \( \mu \) and plot your mean centered data to verify it is
+centered at the origin! Compare your code with the functionality from Scikit-Learn discussed above.
+
+
+Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by:
+$$
+\begin{equation*}
+\Sigma_N = \frac{1}{N-1} \sum_{i=1}^N \bar{x}_i^T \bar{x}_i = \frac{1}{N-1} \sum_{i=1}^N (x_i - \mu_N)^T (x_i - \mu_N)
+\end{equation*}
+$$
+
+where the data points \( x_i \in \mathbb{R}^p \) (here in this example \( p = 2 \)) are column vectors and \( x^T \) is the transpose of \( x \).
+Compare the computed covariance with the answer given above.
+
+
+Now we are ready to solve for the principal components! To do so we
+diagonalize the sample covariance matrix \( \Sigma_N \). We can use the
+function np.linalg.eig to do so. It will return the eigenvalues and
+eigenvectors of \( \Sigma_N \). Once you have these, carry out the
+following tasks:
+
+
+Finally, collect all these steps and write your own PCA function and
+compare this with the functionality included in Scikit-Learn.
+Have the input be the data and have the output be the principal components and their associated eigenvalues, sorted in descending order. Can you think of a way to make it more efficient than the algorithm outlined above?
+
+
+Finally, try out your own PCA function with other data sets.
+
+
After this we ask ourselves how do we prove the link between the maximum variance and the feature reduction.
We assume now that we have a design matrix \( \boldsymbol{X} \) which has been centered as discussed above. For the sake of simplicity we skip the overline symbol. The matrix is defined in terms of the various column vectors \( [\boldsymbol{x}_0,\boldsymbol{x}_1,\dots, \boldsymbol{x}_{p-1}] \)
@@ -976,7 +1084,7 @@ The PCA theorem states that minimizing the above reconstruction error correspond
To show the PCA theorem let us start with the assumption that there is one vector \( \boldsymbol{w}_0 \) which corresponds to a solution which minimized the reconstruction error \( J \). This is an orthogonal vector. It means that we now approximate the reconstruction error in terms of \( \boldsymbol{w}_0 \) and \( \boldsymbol{z}_0 \) as
@@ -999,7 +1107,7 @@ where the vectors on the rhs are known.
We have now found the unknown parameters \( z_{i0} \). These correspond to the projected coordinates and we can write
@@ -1037,7 +1145,7 @@ We are almost there, we have obtained a relation between minimizing the reconstr
We could trivially maximize the variance of the projection (and
@@ -1088,7 +1196,7 @@ chapter 12.4 and discussion therein.
@@ -1144,7 +1252,7 @@ X2D = X_centered.dot(W2)
-
Scikit-Learn’s PCA class implements PCA using SVD decomposition just like we did before. The
@@ -1176,7 +1284,7 @@ variance that lies along the axis of each principal component.
@@ -1217,7 +1325,7 @@ We see that our training data after the PCA decomposition has a performance simi
Instead of arbitrarily choosing the number of dimensions to reduce down to, it is generally preferable to
@@ -1247,7 +1355,7 @@ X_reduced = pca.fit_transform(X)
One problem with the preceding implementation of PCA is that it requires the whole training set to fit in
@@ -1259,7 +1367,7 @@ instances arrive).
Scikit-Learn offers yet another option to perform PCA, called Randomized PCA. This is a stochastic
@@ -1274,7 +1382,7 @@ previous algorithms when \( d \) is much smaller than \( n \).
@@ -1303,7 +1411,7 @@ X_reduced = rbf_pca.fit_transform(X)
Locally Linear Embedding (LLE) is another very powerful nonlinear dimensionality reduction
@@ -1315,7 +1423,7 @@ these local relationships are best preserved (more details shortly).
There are many other dimensionality reduction techniques, several of which are available in Scikit-Learn.
diff --git a/doc/pub/DimRed/html/DimRed.html b/doc/pub/DimRed/html/DimRed.html
index 08bc92a30..7b1ba18c7 100644
--- a/doc/pub/DimRed/html/DimRed.html
+++ b/doc/pub/DimRed/html/DimRed.html
@@ -114,19 +114,30 @@ div { text-align: justify; text-justify: inter-word; }
'___sec14'),
('Towards the PCA theorem', 2, None, '___sec15'),
('The Algorithm before the Theorem', 2, None, '___sec16'),
- ('Classical PCA Theorem', 2, None, '___sec17'),
- ('Proof of the PCA Theorem', 2, None, '___sec18'),
- ('PCA Proof continued', 2, None, '___sec19'),
- ('The final step', 2, None, '___sec20'),
- ('Principal Component Analysis', 2, None, '___sec21'),
- ('PCA and scikit-learn', 2, None, '___sec22'),
- ('Back to the Cancer Data', 2, None, '___sec23'),
- ('More on the PCA', 2, None, '___sec24'),
- ('Incremental PCA', 2, None, '___sec25'),
- ('Randomized PCA', 2, None, '___sec26'),
- ('Kernel PCA', 2, None, '___sec27'),
- ('LLE', 2, None, '___sec28'),
- ('Other techniques', 2, None, '___sec29')]}
+ ('Writing our own PCA code', 2, None, '___sec17'),
+ ('Compute the sample mean and center the data',
+ 3,
+ None,
+ '___sec18'),
+ ('Compute the sample covariance', 3, None, '___sec19'),
+ ('Diagonalize the sample covariance matrix to obtain the '
+ 'principal components',
+ 3,
+ None,
+ '___sec20'),
+ ('Classical PCA Theorem', 2, None, '___sec21'),
+ ('Proof of the PCA Theorem', 2, None, '___sec22'),
+ ('PCA Proof continued', 2, None, '___sec23'),
+ ('The final step', 2, None, '___sec24'),
+ ('Principal Component Analysis', 2, None, '___sec25'),
+ ('PCA and scikit-learn', 2, None, '___sec26'),
+ ('Back to the Cancer Data', 2, None, '___sec27'),
+ ('More on the PCA', 2, None, '___sec28'),
+ ('Incremental PCA', 2, None, '___sec29'),
+ ('Randomized PCA', 2, None, '___sec30'),
+ ('Kernel PCA', 2, None, '___sec31'),
+ ('LLE', 2, None, '___sec32'),
+ ('Other techniques', 2, None, '___sec33')]}
end of tocinfo -->
-
+We will use a simple example first with two-dimensional data
+drawn from a multivariate normal distribution with the following mean and covariance matrix:
+$$
+\mu = (-1,2) \qquad \Sigma = \begin{bmatrix} 4 & 2 \\
+2 & 2
+\end{bmatrix}
+$$
+
+
+We will generate \( N = 1000 \) points \( X = \{ x_1, \ldots, x_N \} \) from
+this distribution, and store them in the \( 1000 \times 2 \) matrix \( \boldsymbol{X} \).
+
+
+The following Python code aids in setting up the data
+
+
+
+
+
+Make a small Python code which plots the data.
+
+
+Now we are going to implement the PCA algorithm. We will break it down into sub-steps and across multiple cells.
+
+
+The first step of PCA is to compute the sample mean of the data and use it to center the data. Recall the sample mean is
+$$
+\mu_N = \frac{1}{N} \sum_{i=1}^N x_i
+$$
+
+and the mean-centered data \( \bar{X} = \{ \bar{x}_1, \ldots, \bar{x}_N \} \) takes the form
+$$
+\bar{x}_i = x_i - \mu_N
+$$
+
+When you are done with these steps, print out \( \mu_N \) to verify it is
+close to \( \mu \) and plot your mean centered data to verify it is
+centered at the origin! Compare your code with the functionality from Scikit-Learn discussed above.
+
+
+Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by:
+$$
+\begin{equation*}
+\Sigma_N = \frac{1}{N-1} \sum_{i=1}^N \bar{x}_i^T \bar{x}_i = \frac{1}{N-1} \sum_{i=1}^N (x_i - \mu_N)^T (x_i - \mu_N)
+\end{equation*}
+$$
+
+where the data points \( x_i \in \mathbb{R}^p \) (here in this example \( p = 2 \)) are column vectors and \( x^T \) is the transpose of \( x \).
+Compare the computed covariance with the answer given above.
+
+
+Now we are ready to solve for the principal components! To do so we
+diagonalize the sample covariance matrix \( \Sigma_N \). We can use the
+function np.linalg.eig to do so. It will return the eigenvalues and
+eigenvectors of \( \Sigma_N \). Once you have these, carry out the
+following tasks:
+
+
+Finally, collect all these steps and write your own PCA function and
+compare this with the functionality included in Scikit-Learn.
+Have the input be the data and have the output be the principal components and their associated eigenvalues, sorted in descending order. Can you think of a way to make it more efficient than the algorithm outlined above?
+
+
+Finally, try out your own PCA function with other data sets.
+
+
After this we ask ourselves how do we prove the link between the maximum variance and the feature reduction.
We assume now that we have a design matrix \( \boldsymbol{X} \) which has been centered as discussed above. For the sake of simplicity we skip the overline symbol. The matrix is defined in terms of the various column vectors \( [\boldsymbol{x}_0,\boldsymbol{x}_1,\dots, \boldsymbol{x}_{p-1}] \)
@@ -981,7 +1089,7 @@ The PCA theorem states that minimizing the above reconstruction error correspond
To show the PCA theorem let us start with the assumption that there is one vector \( \boldsymbol{w}_0 \) which corresponds to a solution which minimized the reconstruction error \( J \). This is an orthogonal vector. It means that we now approximate the reconstruction error in terms of \( \boldsymbol{w}_0 \) and \( \boldsymbol{z}_0 \) as
@@ -1004,7 +1112,7 @@ where the vectors on the rhs are known.
We have now found the unknown parameters \( z_{i0} \). These correspond to the projected coordinates and we can write
@@ -1042,7 +1150,7 @@ We are almost there, we have obtained a relation between minimizing the reconstr
We could trivially maximize the variance of the projection (and
@@ -1093,7 +1201,7 @@ chapter 12.4 and discussion therein.
@@ -1149,7 +1257,7 @@ X2D = X_centered
-
Scikit-Learn’s PCA class implements PCA using SVD decomposition just like we did before. The
@@ -1181,7 +1289,7 @@ variance that lies along the axis of each principal component.
@@ -1222,7 +1330,7 @@ We see that our training data after the PCA decomposition has a performance simi
Instead of arbitrarily choosing the number of dimensions to reduce down to, it is generally preferable to
@@ -1252,7 +1360,7 @@ X_reduced = pca
One problem with the preceding implementation of PCA is that it requires the whole training set to fit in
@@ -1264,7 +1372,7 @@ instances arrive).
Scikit-Learn offers yet another option to perform PCA, called Randomized PCA. This is a stochastic
@@ -1279,7 +1387,7 @@ previous algorithms when \( d \) is much smaller than \( n \).
@@ -1308,7 +1416,7 @@ X_reduced = rbf_pcaLLE
+
Locally Linear Embedding (LLE) is another very powerful nonlinear dimensionality reduction
@@ -1320,7 +1428,7 @@ these local relationships are best preserved (more details shortly).
There are many other dimensionality reduction techniques, several of which are available in Scikit-Learn.
diff --git a/doc/pub/DimRed/ipynb/DimRed.ipynb b/doc/pub/DimRed/ipynb/DimRed.ipynb
index 700af1f72..82b60b04e 100644
--- a/doc/pub/DimRed/ipynb/DimRed.ipynb
+++ b/doc/pub/DimRed/ipynb/DimRed.ipynb
@@ -10,7 +10,7 @@
" \n",
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
"\n",
- "Date: **Oct 25, 2019**\n",
+ "Date: **Dec 26, 2019**\n",
"\n",
"Copyright 1999-2019, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
"\n",
@@ -1049,6 +1049,156 @@
"\n",
"* Keep only those $l$ eigenvalues larger than a selected threshold value, discarding thus $p-l$ features since we expect small variations in the data here.\n",
"\n",
+ "## Writing our own PCA code\n",
+ "\n",
+ "We will use a simple example first with two-dimensional data\n",
+ "drawn from a multivariate normal distribution with the following mean and covariance matrix:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\mu = (-1,2) \\qquad \\Sigma = \\begin{bmatrix} 4 & 2 \\\\\n",
+ "2 & 2\n",
+ "\\end{bmatrix}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We will generate $N = 1000$ points $X = \\{ x_1, \\ldots, x_N \\}$ from\n",
+ "this distribution, and store them in the $1000 \\times 2$ matrix $\\boldsymbol{X}$.\n",
+ "\n",
+ "The following Python code aids in setting up the data"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "N = 1000\n",
+ "mean = (-1, 2)\n",
+ "cov = [[4, 2], [2, 2]]\n",
+ "X = np.random.multivariate_normal(mean, cov, N)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Make a small Python code which plots the data.\n",
+ "\n",
+ "Now we are going to implement the PCA algorithm. We will break it down into sub-steps and across multiple cells.\n",
+ "\n",
+ "### Compute the sample mean and center the data\n",
+ "\n",
+ "The first step of PCA is to compute the sample mean of the data and use it to center the data. Recall the sample mean is"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\mu_N = \\frac{1}{N} \\sum_{i=1}^N x_i\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and the mean-centered data $\\bar{X} = \\{ \\bar{x}_1, \\ldots, \\bar{x}_N \\}$ takes the form"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\bar{x}_i = x_i - \\mu_N\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "When you are done with these steps, print out $\\mu_N$ to verify it is\n",
+ "close to $\\mu$ and plot your mean centered data to verify it is\n",
+ "centered at the origin! Compare your code with the functionality from **Scikit-Learn** discussed above.\n",
+ "\n",
+ "\n",
+ "### Compute the sample covariance\n",
+ "\n",
+ "Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "\\Sigma_N = \\frac{1}{N-1} \\sum_{i=1}^N \\bar{x}_i^T \\bar{x}_i = \\frac{1}{N-1} \\sum_{i=1}^N (x_i - \\mu_N)^T (x_i - \\mu_N)\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where the data points $x_i \\in \\mathbb{R}^p$ (here in this example $p = 2$) are column vectors and $x^T$ is the transpose of $x$.\n",
+ "Compare the computed covariance with the answer given above.\n",
+ "\n",
+ "\n",
+ "### Diagonalize the sample covariance matrix to obtain the principal components\n",
+ "\n",
+ "Now we are ready to solve for the principal components! To do so we\n",
+ "diagonalize the sample covariance matrix $\\Sigma_N$. We can use the\n",
+ "function **np.linalg.eig** to do so. It will return the eigenvalues and\n",
+ "eigenvectors of $\\Sigma_N$. Once you have these, carry out the\n",
+ "following tasks:\n",
+ "\n",
+ "* Compute the percentage of the total variance captured by the first principal component\n",
+ "\n",
+ "* Plot the mean centered data and lines along the first and second principal components\n",
+ "\n",
+ "* Project the mean centered data onto the first and second principal components, and plot the projected data. What do you observe?\n",
+ "\n",
+ "* Approximate the data as"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "$$\n",
+ "x_i \\approx \\tilde{x}_i := \\mu_N + \\langle x_i, v_0 \\rangle v_0\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "where $v_0$ is the first principal component. What do you observe?\n",
+ "\n",
+ "Finally, collect all these steps and write your own PCA function and\n",
+ "compare this with the functionality included in **Scikit-Learn**.\n",
+ "Have the input be the data and have the output be the principal components and their associated eigenvalues, sorted in descending order. Can you think of a way to make it more efficient than the algorithm outlined above?\n",
+ "\n",
+ "Finally, try out your own PCA function with other data sets.\n",
+ "\n",
+ "\n",
+ "\n",
"After this we ask ourselves how do we prove the link between the maximum variance and the feature reduction.\n",
"\n",
"## Classical PCA Theorem\n",
@@ -1320,7 +1470,7 @@
},
{
"cell_type": "code",
- "execution_count": 11,
+ "execution_count": 12,
"metadata": {
"collapsed": false
},
@@ -1367,7 +1517,7 @@
},
{
"cell_type": "code",
- "execution_count": 12,
+ "execution_count": 13,
"metadata": {
"collapsed": false
},
@@ -1391,7 +1541,7 @@
},
{
"cell_type": "code",
- "execution_count": 13,
+ "execution_count": 14,
"metadata": {
"collapsed": false
},
@@ -1415,7 +1565,7 @@
},
{
"cell_type": "code",
- "execution_count": 14,
+ "execution_count": 15,
"metadata": {
"collapsed": false
},
@@ -1439,7 +1589,7 @@
},
{
"cell_type": "code",
- "execution_count": 15,
+ "execution_count": 16,
"metadata": {
"collapsed": false
},
@@ -1493,7 +1643,7 @@
},
{
"cell_type": "code",
- "execution_count": 16,
+ "execution_count": 17,
"metadata": {
"collapsed": false
},
@@ -1516,7 +1666,7 @@
},
{
"cell_type": "code",
- "execution_count": 17,
+ "execution_count": 18,
"metadata": {
"collapsed": false
},
@@ -1563,7 +1713,7 @@
},
{
"cell_type": "code",
- "execution_count": 18,
+ "execution_count": 19,
"metadata": {
"collapsed": false
},
diff --git a/doc/pub/DimRed/ipynb/ipynb-DimRed-src.tar.gz b/doc/pub/DimRed/ipynb/ipynb-DimRed-src.tar.gz
index cb7d0f3ef..55bb01bf7 100644
Binary files a/doc/pub/DimRed/ipynb/ipynb-DimRed-src.tar.gz and b/doc/pub/DimRed/ipynb/ipynb-DimRed-src.tar.gz differ
diff --git a/doc/pub/DimRed/pdf/DimRed-minted.pdf b/doc/pub/DimRed/pdf/DimRed-minted.pdf
index b934669d3..8712e83bf 100644
Binary files a/doc/pub/DimRed/pdf/DimRed-minted.pdf and b/doc/pub/DimRed/pdf/DimRed-minted.pdf differ
diff --git a/doc/src/DimRed/DimRed.do.txt b/doc/src/DimRed/DimRed.do.txt
index a8de47b75..0fdb494ba 100644
--- a/doc/src/DimRed/DimRed.do.txt
+++ b/doc/src/DimRed/DimRed.do.txt
@@ -739,6 +739,94 @@ x_{n-1,0} & x_{n-1,1} & x_{n-1,2}& \dots & \dots x_{n-1,p-1}\\
* Order the eigenvalue (and the eigenvectors accordingly) in order of decreasing eigenvalues.
* Keep only those $l$ eigenvalues larger than a selected threshold value, discarding thus $p-l$ features since we expect small variations in the data here.
+
+!split
+===== Writing our own PCA code =====
+
+We will use a simple example first with two-dimensional data
+drawn from a multivariate normal distribution with the following mean and covariance matrix:
+!bt
+\[
+\mu = (-1,2) \qquad \Sigma = \begin{bmatrix} 4 & 2 \\
+2 & 2
+\end{bmatrix}
+\]
+!et
+
+We will generate $N = 1000$ points $X = \{ x_1, \ldots, x_N \}$ from
+this distribution, and store them in the $1000 \times 2$ matrix $\bm{X}$.
+
+The following Python code aids in setting up the data
+
+!bc pycod
+N = 1000
+mean = (-1, 2)
+cov = [[4, 2], [2, 2]]
+X = np.random.multivariate_normal(mean, cov, N)
+!ec
+
+Make a small Python code which plots the data.
+
+Now we are going to implement the PCA algorithm. We will break it down into sub-steps and across multiple cells.
+
+=== Compute the sample mean and center the data ===
+
+The first step of PCA is to compute the sample mean of the data and use it to center the data. Recall the sample mean is
+!bt
+\[
+\mu_N = \frac{1}{N} \sum_{i=1}^N x_i
+\]
+!et
+and the mean-centered data $\bar{X} = \{ \bar{x}_1, \ldots, \bar{x}_N \}$ takes the form
+!bt
+\[
+\bar{x}_i = x_i - \mu_N
+\]
+!et
+When you are done with these steps, print out $\mu_N$ to verify it is
+close to $\mu$ and plot your mean centered data to verify it is
+centered at the origin! Compare your code with the functionality from _Scikit-Learn_ discussed above.
+
+
+=== Compute the sample covariance ===
+
+Now we are going to use the mean centered data to compute the sample covariance of the data. Recall it is given by:
+!bt
+\begin{equation*}
+\Sigma_N = \frac{1}{N-1} \sum_{i=1}^N \bar{x}_i^T \bar{x}_i = \frac{1}{N-1} \sum_{i=1}^N (x_i - \mu_N)^T (x_i - \mu_N)
+\end{equation*}
+!et
+where the data points $x_i \in \mathbb{R}^p$ (here in this example $p = 2$) are column vectors and $x^T$ is the transpose of $x$.
+Compare the computed covariance with the answer given above.
+
+
+=== Diagonalize the sample covariance matrix to obtain the principal components ===
+
+Now we are ready to solve for the principal components! To do so we
+diagonalize the sample covariance matrix $\Sigma_N$. We can use the
+function _np.linalg.eig_ to do so. It will return the eigenvalues and
+eigenvectors of $\Sigma_N$. Once you have these, carry out the
+following tasks:
+
+* Compute the percentage of the total variance captured by the first principal component
+* Plot the mean centered data and lines along the first and second principal components
+* Project the mean centered data onto the first and second principal components, and plot the projected data. What do you observe?
+* Approximate the data as
+!bt
+\begin{equation*}
+x_i \approx \tilde{x}_i := \mu_N + \langle x_i, v_0 \rangle v_0
+\end{equation*}
+!et
+where $v_0$ is the first principal component. What do you observe?
+
+Finally, collect all these steps and write your own PCA function and
+compare this with the functionality included in _Scikit-Learn_.
+Have the input be the data and have the output be the principal components and their associated eigenvalues, sorted in descending order. Can you think of a way to make it more efficient than the algorithm outlined above?
+
+Finally, try out your own PCA function with other data sets.
+
+
+
After this we ask ourselves how do we prove the link between the maximum variance and the feature reduction.
!split
PCA and scikit-learn
+Principal Component Analysis
+#thereafter we do a PCA with Scikit-learn
-from sklearn.decomposition import PCA
-pca = PCA(n_components = 2)
-X2D = pca.fit_transform(X)
+
import numpy as np
+import pandas as pd
+from IPython.display import display
+np.random.seed(100)
+# setting up a 10 x 5 vanilla matrix
+rows = 10
+cols = 5
+X = np.random.randn(rows,cols)
+df = pd.DataFrame(X)
+# Pandas does the centering for us
+df = df -df.mean()
+display(df)
+
+# we center it ourselves
+X_centered = X - X.mean(axis=0)
+# Then check the difference between pandas and our own set up
+print(X_centered-df)
+#Now we do an SVD
+U, s, V = np.linalg.svd(X_centered)
+c1 = V.T[:, 0]
+c2 = V.T[:, 1]
+W2 = V.T[:, :2]
+X2D = X_centered.dot(W2)
print(X2D)
pca.components_.T[:, 0].
+
W2 = V.T[:, :2]
+X2D = X_centered.dot(W2)
Back to the Cancer Data
-We can now repeat the above but applied to real data, in this case our breast cancer data.
-Here we compute performance scores on the training data using logistic regression.
+PCA and scikit-learn
+
+import matplotlib.pyplot as plt
-import numpy as np
-from sklearn.model_selection import train_test_split
-from sklearn.datasets import load_breast_cancer
-from sklearn.linear_model import LogisticRegression
-cancer = load_breast_cancer()
-
-X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
-
-logreg = LogisticRegression()
-logreg.fit(X_train, y_train)
-print("Train set accuracy from Logistic Regression: {:.2f}".format(logreg.score(X_train,y_train)))
-# We scale the data
-from sklearn.preprocessing import StandardScaler
-scaler = StandardScaler()
-scaler.fit(X_train)
-X_train_scaled = scaler.transform(X_train)
-X_test_scaled = scaler.transform(X_test)
-# Then perform again a log reg fit
-logreg.fit(X_train_scaled, y_train)
-print("Train set accuracy scaled data: {:.2f}".format(logreg.score(X_train_scaled,y_train)))
-#thereafter we do a PCA with Scikit-learn
+
#thereafter we do a PCA with Scikit-learn
from sklearn.decomposition import PCA
pca = PCA(n_components = 2)
-X2D_train = pca.fit_transform(X_train_scaled)
-# and finally compute the log reg fit and the score on the training data
-logreg.fit(X2D_train,y_train)
-print("Train set accuracy scaled and PCA data: {:.2f}".format(logreg.score(X2D_train,y_train)))
+X2D = pca.fit_transform(X)
+print(X2D)
pca.components_.T[:, 0].
+
More on the PCA
-
-Back to the Cancer Data
+We can now repeat the above but applied to real data, in this case our breast cancer data.
+Here we compute performance scores on the training data using logistic regression.
pca = PCA()
-pca.fit(X)
-cumsum = np.cumsum(pca.explained_variance_ratio_)
-d = np.argmax(cumsum >= 0.95) + 1
-
import matplotlib.pyplot as plt
+import numpy as np
+from sklearn.model_selection import train_test_split
+from sklearn.datasets import load_breast_cancer
+from sklearn.linear_model import LogisticRegression
+cancer = load_breast_cancer()
-
-
pca = PCA(n_components=0.95)
-X_reduced = pca.fit_transform(X)
+X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
+
+logreg = LogisticRegression()
+logreg.fit(X_train, y_train)
+print("Train set accuracy from Logistic Regression: {:.2f}".format(logreg.score(X_train,y_train)))
+# We scale the data
+from sklearn.preprocessing import StandardScaler
+scaler = StandardScaler()
+scaler.fit(X_train)
+X_train_scaled = scaler.transform(X_train)
+X_test_scaled = scaler.transform(X_test)
+# Then perform again a log reg fit
+logreg.fit(X_train_scaled, y_train)
+print("Train set accuracy scaled data: {:.2f}".format(logreg.score(X_train_scaled,y_train)))
+#thereafter we do a PCA with Scikit-learn
+from sklearn.decomposition import PCA
+pca = PCA(n_components = 2)
+X2D_train = pca.fit_transform(X_train_scaled)
+# and finally compute the log reg fit and the score on the training data
+logreg.fit(X2D_train,y_train)
+print("Train set accuracy scaled and PCA data: {:.2f}".format(logreg.score(X2D_train,y_train)))
Incremental PCA
+More on the PCA
pca = PCA()
+pca.fit(X)
+cumsum = np.cumsum(pca.explained_variance_ratio_)
+d = np.argmax(cumsum >= 0.95) + 1
+
pca = PCA(n_components=0.95)
+X_reduced = pca.fit_transform(X)
+
Randomized PCA
+Incremental PCA
Kernel PCA
-Randomized PCA
from sklearn.decomposition import KernelPCA
-rbf_pca = KernelPCA(n_components = 2, kernel="rbf", gamma=0.04)
-X_reduced = rbf_pca.fit_transform(X)
-
LLE
+Kernel PCA
+from sklearn.decomposition import KernelPCA
+rbf_pca = KernelPCA(n_components = 2, kernel="rbf", gamma=0.04)
+X_reduced = rbf_pca.fit_transform(X)
+
Other techniques
+LLE
-
-
-
Oct 25, 2019
Dec 26, 2019
-Oct 25, 2019
Dec 26, 2019
Writing our own PCA code
+
+
+$$
+\mu = (-1,2) \qquad \Sigma = \begin{bmatrix} 4 & 2 \\
+2 & 2
+\end{bmatrix}
+$$
+
+
+N = 1000
+mean = (-1, 2)
+cov = [[4, 2], [2, 2]]
+X = np.random.multivariate_normal(mean, cov, N)
+
Compute the sample mean and center the data
+
+
+$$
+\mu_N = \frac{1}{N} \sum_{i=1}^N x_i
+$$
+
+
+and the mean-centered data \( \bar{X} = \{ \bar{x}_1, \ldots, \bar{x}_N \} \) takes the form
+
+$$
+\bar{x}_i = x_i - \mu_N
+$$
+
+
+When you are done with these steps, print out \( \mu_N \) to verify it is
+close to \( \mu \) and plot your mean centered data to verify it is
+centered at the origin! Compare your code with the functionality from Scikit-Learn discussed above.
+
+Compute the sample covariance
+
+
+$$
+\begin{equation*}
+\Sigma_N = \frac{1}{N-1} \sum_{i=1}^N \bar{x}_i^T \bar{x}_i = \frac{1}{N-1} \sum_{i=1}^N (x_i - \mu_N)^T (x_i - \mu_N)
+\end{equation*}
+$$
+
+
+where the data points \( x_i \in \mathbb{R}^p \) (here in this example \( p = 2 \)) are column vectors and \( x^T \) is the transpose of \( x \).
+Compare the computed covariance with the answer given above.
+
+Diagonalize the sample covariance matrix to obtain the principal components
+
+
+
+
+$$
+\begin{equation*}
+x_i \approx \tilde{x}_i := \mu_N + \langle x_i, v_0 \rangle v_0
+\end{equation*}
+$$
+
+
+where \( v_0 \) is the first principal component. What do you observe?
+
+Classical PCA Theorem
+Classical PCA Theorem
Proof of the PCA Theorem
+Proof of the PCA Theorem
PCA Proof continued
+PCA Proof continued
The final step
+The final step
Principal Component Analysis
+Principal Component Analysis
PCA and scikit-learn
+PCA and scikit-learn
Back to the Cancer Data
+Back to the Cancer Data
We can now repeat the above but applied to real data, in this case our breast cancer data.
Here we compute performance scores on the training data using logistic regression.
More on the PCA
+More on the PCA
Incremental PCA
+Incremental PCA
Randomized PCA
+Randomized PCA
Kernel PCA
+Kernel PCA
LLE
+LLE
Other techniques
+Other techniques
Oct 25, 2019
Dec 26, 2019
@@ -950,12 +961,109 @@ $$
+
+Writing our own PCA code
+
+N = 1000
+mean = (-1, 2)
+cov = [[4, 2], [2, 2]]
+X = np.random.multivariate_normal(mean, cov, N)
+
Compute the sample mean and center the data
+
+Compute the sample covariance
+
+Diagonalize the sample covariance matrix to obtain the principal components
+
+
+
+
+$$
+\begin{equation*}
+x_i \approx \tilde{x}_i := \mu_N + \langle x_i, v_0 \rangle v_0
+\end{equation*}
+$$
+
+where \( v_0 \) is the first principal component. What do you observe?
+
+
-Classical PCA Theorem
+Classical PCA Theorem
-Proof of the PCA Theorem
+Proof of the PCA Theorem
-PCA Proof continued
+PCA Proof continued
-The final step
+The final step
-Principal Component Analysis
+Principal Component Analysis
PCA and scikit-learn
+PCA and scikit-learn
-Back to the Cancer Data
+Back to the Cancer Data
We can now repeat the above but applied to real data, in this case our breast cancer data.
Here we compute performance scores on the training data using logistic regression.
-More on the PCA
+More on the PCA
-Incremental PCA
+Incremental PCA
-Randomized PCA
+Randomized PCA
-Kernel PCA
+Kernel PCA
-LLE
+LLE
-Other techniques
+Other techniques
Oct 25, 2019
Dec 26, 2019
@@ -955,12 +966,109 @@ $$
+
+Writing our own PCA code
+
+N = 1000
+mean = (-1, 2)
+cov = [[4, 2], [2, 2]]
+X = np.random.multivariate_normal(mean, cov, N)
+
Compute the sample mean and center the data
+
+Compute the sample covariance
+
+Diagonalize the sample covariance matrix to obtain the principal components
+
+
+
+
+$$
+\begin{equation*}
+x_i \approx \tilde{x}_i := \mu_N + \langle x_i, v_0 \rangle v_0
+\end{equation*}
+$$
+
+where \( v_0 \) is the first principal component. What do you observe?
+
+
-Classical PCA Theorem
+Classical PCA Theorem
-Proof of the PCA Theorem
+Proof of the PCA Theorem
-PCA Proof continued
+PCA Proof continued
-The final step
+The final step
-Principal Component Analysis
+Principal Component Analysis
PCA and scikit-learn
+PCA and scikit-learn
-Back to the Cancer Data
+Back to the Cancer Data
We can now repeat the above but applied to real data, in this case our breast cancer data.
Here we compute performance scores on the training data using logistic regression.
-More on the PCA
+More on the PCA
-Incremental PCA
+Incremental PCA
-Randomized PCA
+Randomized PCA
-Kernel PCA
+Kernel PCA
LLE
-Other techniques
+Other techniques