<jupyter_start><jupyter_text>Principal Components AnalysisPrincipal Components Analysis (PCA) is a widely used technique for dimensionality reduction. Learning Objectives* Gain an understanding of PCA and how it relates to covariance matrices.  Variance and Covariance MatricesThe cell below plots a small two-dimensional dataset.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

X = np.array([[-4.40394652,  5.7203612 ],
       [-5.4420319 ,  1.93587109],
       [-4.56570001,  4.2530108 ],
       [-5.07827013,  1.51210658],
       [-5.39695317,  1.25160109],
       [-4.28688083,  4.07841531],
       [-3.58299231,  6.30384398],
       [-5.60152663,  2.24470856],
       [-4.47414671,  3.54114625],
       [-5.00867538,  3.37502121],
       [-3.61698704,  6.57347194],
       [-7.14400941, -1.24682603],
       [-8.34978399, -3.85127918],
       [-7.35710211, -1.33622093],
       [-2.82445297,  7.90898657],
       [-4.65304533,  3.10063785],
       [-3.41097162,  4.96990604],
       [-5.91298693,  1.85557319],
       [-6.65790912, -0.91626863],
       [-3.67245142,  5.35009044]])


plt.plot(X[:, 0], X[:, 1], '.')
plt.axis('equal')
plt.show()<jupyter_output><empty_output><jupyter_text>Run the cell above, then answer the following questions. Questions:*   Approximately where is the mean of this data set?  In other words, what is $[\overline{x}_0, \overline{x}_1]^T$ (where the subscripts corrspond to features)?*   If $\mathbf{S}$ is the sample covariance matrix of this data set,     $$\mathbf{S} = \begin{bmatrix} a & b \\ c & d \end{bmatrix}$$    * Which entry in $\mathbf{S}$ corresponds to the variance of feature 0?     * Which entry in $\mathbf{S}$ corresponds to the variance of feature 1?     * What is the value of $a$ relative to $d$?  Larger?  Smaller? Equal?    * What is the sign of $b$?    * What is the value of $b$ relative to $c$?  Larger? Smaller? Equal?*   If you could only keep one of the two features in this data set, which one would you pick.  Why? Answers:*     *        *     *     *     *     *   Centering the dataThe first step in the PCA algorithm is to center the data by subtracting out the mean value of each feature:<jupyter_code>X_mean = np.mean(X, axis=0)
print("The mean is: ", X_mean)

X = X - X_mean # Center the data by subtracting out the mean.

plt.plot(X[:, 0], X[:, 1], '.')
plt.axis('equal')
plt.show()<jupyter_output>The mean is:  [-5.07204118  2.83120787]<jupyter_text>Finding the Sample Covariance MatrixThe next step is to calculate the sample covariance matrix.  The sample covariance matrix can be expressed as: $$\mathbf{S} = \frac{1}{n-1} \sum_1^n (X_i - \overline{X})(X_i - \overline{X})^T $$Where $X_i$ represents the $i$th row of $X$ and $\overline{X}$ is the sample mean. Of course... In this case, we've already centered the data, so $\overline{X}$ is $\mathbf{0}$ and the formula reduces to:$$\mathbf{S} = \frac{1}{n-1} \sum_1^n X_i X_i^T $$This is equivalent to:$$\mathbf{S} = \frac{1}{n-1} X^T X $$Where $X$ is the full, centered, data matrix.  Lets see if it works!Once you've executed the cell below, go back and check the answers you entered above.<jupyter_code>n = X.shape[0] # grab the number of samples

S = np.dot(X.T, X) / (n-1)

print("X.T times X:\n", S)

print("\nNumpy's version:\n", np.cov(X.T))<jupyter_output>X.T times X:
 [[2.1033442  4.26897913]
 [4.26897913 9.18908776]]

Numpy's version:
 [[2.1033442  4.26897913]
 [4.26897913 9.18908776]]<jupyter_text>EigenvectorsThe *eigenvectors* of a matrix are special vectors that don't change direction when multiplied by that matrix. In other words, multiplying a matrix by one of its eigenvectors gives the same answer as multiplying by some scalar:$$\mathbf{A} \mathbf{v} = \lambda \mathbf{v}$$Here the scalar $\lambda$ is the known as an *eigenvalue*.Consider the matrix $$\mathbf{A} = \begin{bmatrix} 1 & .5 \\ .5 & 1 \end{bmatrix}$$Its eigenvectors are:$$\mathbf{v}_1 = \begin{bmatrix} \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix}, ~~\mathbf{v}_2 = \begin{bmatrix} -\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix}$$  Question* Calculate by hand $\mathbf{A} \mathbf{v}_1$ and $\mathbf{A} \mathbf{v}_1$.  What are $\lambda_1$ and $\lambda_2$? Answer* (Note that these two eigenvectors aren't unique.  We can rescale them, and the results will still be eigenvectors.  It is common to use eigenvectors that are scaled to unit length.)  Principal ComponentsIt turns out that we can find the directions of maximum variance in a dataset by finding the eigenvectors of the covariance matrix.  The eigenvector with the largest eigenvalue points in the direction of maximum variance. Before executing the next cell, take a look back at the dataset above to predict the direction of the first principal component.<jupyter_code># Numpy is able to calculate eigenvalues and eigenvectors.
#   `lambdas` contains the eigenvalues (in sorted order)
#    The columms of `W` contain the corresponding eigenvectors.
lambdas, W = np.linalg.eigh(S)

# Lets reverse the order to get the largest first...

lambdas = lambdas[::-1]
W = W[:, ::-1]

print("Eigenvalues:", lambdas)
print("Eigenvectors:\n", W)

# Replot the data points
plt.plot(X[:, 0], X[:, 1], '.')
plt.axis('equal')

# Show the first principal component.
plt.arrow(0, 0, W[0,0] * 3, W[1,0] * 3, head_width=0.5, head_length=0.5)
plt.title('First Principal Component')
plt.show()

# Show the second principal component.
plt.figure()
plt.plot(X[:, 0], X[:, 1], '.')
plt.axis('equal')
plt.arrow(0, 0, W[0,1] * 3, W[1,1] * 3, head_width=0.5, head_length=0.5)
plt.title('Second Principal Component')
plt.show()<jupyter_output>Eigenvalues: [11.19383918  0.09859277]
Eigenvectors:
 [[ 0.42507131 -0.90515987]
 [ 0.90515987  0.42507131]]<jupyter_text>Percentage VarianceThe magnitude of the eigenvectors are exactly equal to the amount of variance along the corresponding eigenvector.  If we normalize the eigenvalues, they give us the fraction of the total variance explained by each principal component:<jupyter_code>print(lambdas / np.sum(lambdas))<jupyter_output>[0.99126913 0.00873087]<jupyter_text>In this case, around 99% of the variance is captured by the first principal component.  Projecting Data Onto Principal ComponentsOnce we have our principal components, we would like to know were each data point "lands" on each component.  In other words, we would like to project our data onto the principal components:Here **A** corresponds to a data point and **B** corresponds to a principal component.  In the case where **B** is a unit vector, $|\mathbf{A}| \cos(\theta) = \mathbf{A} \cdot \mathbf{B}$.  In other words, the projection of **A** onto **B** is just the dot product between the two.  Using this information, we can project any one of our data points onto our first principal component. For example:<jupyter_code># Plot all the original data and the first PC
plt.axis('equal')
plt.plot(X[:, 0], X[:, 1], '.', alpha=.2)
plt.arrow(0, 0, W[0,0]*5, W[1,0]*5, head_width=0.5, head_length=0.5)

# Pull out a single poing and project it onto the first PC:
x0 = X[0,:]
x0_projected = np.dot(x0, W[:,0])
print("Projection: ", x0_projected)

# Plot the original point
plt.plot(x0[0], x0[1], 'og')

# Plot where the original point is projected onto the PC:
plt.plot(W[0,0] * x0_projected, W[1,0] * x0_projected, 'or')
plt.legend(['', 'Original Point', 'Projected Point'])

plt.show()

print(x0_projected)<jupyter_output>Projection:  2.899133515105061<jupyter_text>Projecting ALL the dataWe can get a rotated version of our original, centered, dataset by projecting every data point onto our two principal components.  This can be accomplished at one go by multiplying our data matrix by our matrix of principal components:<jupyter_code>X_new = np.dot(X, W)
# Replot the data points
plt.plot(X_new[:, 0], X_new[:, 1], '.')
plt.axis('equal')
plt.show()<jupyter_output><empty_output><jupyter_text>In this figure, the x-axis shows the projection of the data onto the first principal component, and the y-axis shows the projection of the data onto the second principal component. Questions:You should be able to predict all of the values in the covariance matrix of this rotated version of the data without doing any additional calculations.  What are they?$$\mathbf{S} = \begin{bmatrix} a & b \\ c & d \end{bmatrix}$$ * $a$:  * $b$:  * $c$:  * $d$:       Answers* * * *   Reducing DimensionalityIn the example above, we multiplied our data matrix by a matrix containing *all* of the principal components.  We can reduce the dimensionality of the data by leaving out the principal components with smaller eigenvalues:<jupyter_code>Z = np.dot(X, W[:,0]) # Only the first column!

# Replot the data points
plt.plot(Z, np.zeros(Z.shape), '.')
plt.axis('equal')
plt.show()<jupyter_output><empty_output><jupyter_text>Moving Back to the Original SpaceOne way of viewing PCA is as a lossy compression algorithm: we compress the data by removing dimensions that don't carry much information.  We can then decompress by projecting our data back into the original space.  This can be done by multiplying our compressed data by the traspose of the vector of principal components that were used for compression:<jupyter_code>X_hat = np.outer(Z, W[:, 0].T) # need to use `outer `because stupid numpy 
                               # treats 1d arrays differently from column vectors
X_hat += X_mean # restore the mean

plt.plot(X_hat[:,0], X_hat[:,1], '.')
plt.axis('equal')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12785196088739353, 'fraction_numerical': 0.06760592997581748, 'mean_word_length': 4.1695652173913045, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '332242', 'license': 'MIT', 'n_tokens_mistral': 3526, 'n_tokens_neox': 3096, 'n_words': 1147, 'path': 'pca.ipynb', 'repo_name': 'spragunr/jmu_ml_pca', 'size': '9511'}