<jupyter_start><jupyter_text>Principal Component Analysis - Key Ideas PCA finds principal components as linear combinations of the existing features and uses these components to represent the original data. The number of components is a hyperparameter that determines the target dimensionality and needs to be equal to or smaller than the number of observations or columns, whichever is smaller. PCA aims to capture most of the variance in the data, to make it easy to recover the original features, and that each component adds information. It reduces dimensionality by projecting the original data into the principal component space.   Imports & Settings<jupyter_code>%matplotlib inline
import pandas as pd
import numpy as np
from numpy.linalg import lstsq
from numpy.random import randn, seed
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import seaborn as sns
from sklearn.decomposition import PCA
# plt.style.use('fivethirtyeight')
sns.set_style('whitegrid')
pd.options.display.float_format = '{:,.2f}'.format
seed(42)<jupyter_output><empty_output><jupyter_text>Linear Projection: Principal Component Analysis The PCA algorithm works by identifying a sequence of principal components, each of which aligns with the direction of maximum variance in the data after accounting for variation captured by previously computed components. The sequential optimization also ensures that new components are not correlated with existing components so that the resulting set constitutes an orthogonal basis for a vector space.This new basis corresponds to a rotated version of the original basis so that the new axes point in the direction of successively decreasing variance. The decline in the amount of variance of the original data explained by each principal component reflects the extent of correlation among the original features. The number of components that captures, e.g., 95% of the original variation relative to the total number of features provides insights into the linearly independent information in the original data. The below figures illustrates several aspects of PCA for a two-dimensional random data set.  Create Noisy, Correlated Data from Signal<jupyter_code>n_signals = 250
x1 = np.linspace(-10, 10, n_signals) + .1 * randn(n_signals)
x2 = 1.5 * x1 + 2 * randn(n_signals)
data = pd.DataFrame({'$x_1$': x1, '$x_2$': x2})
ax = data.plot.scatter(x=0, y=1, s=10, title='2D Noisy Data', figsize=(7,8))
ax.set_aspect('equal')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Compute Principal Components<jupyter_code>pca = PCA()
pca.fit(data)
pca.components_
mean = pca.mean_
mean
pc1, pc2 = np.split(pca.components_.T, 2, axis=1)
pc1<jupyter_output><empty_output><jupyter_text>Check share of explained variance<jupyter_code>pca.explained_variance_ratio_<jupyter_output><empty_output><jupyter_text>Components are orthogonal to each other<jupyter_code>np.dot(pc1.T, pc2)<jupyter_output><empty_output><jupyter_text>Plot Principal Components as new Basis Vectors The below figure shows how the first and second principal components align with the directions of maximum variance while being orthogonal.<jupyter_code>l1, l2 = pca.singular_values_ / 10
ax = data.plot.scatter(x=0, y=1, s=15, title='Principal Component Vectors', figsize=(6,8))
ax.set_aspect('equal')
origin_x, origin_y = pca.mean_
dx1, dy1 = np.squeeze(pc1.T) * l1
dx2, dy2 = np.squeeze(pc2.T) * l2
pc1_arrow = ax.arrow(origin_x, origin_y, dx1, dy1,  width=.2, color='k')
pc2_arrow = ax.arrow(origin_x, origin_y, dx2, dy2,  width=.2, color='r')
plt.legend([pc1_arrow, pc2_arrow], 
           ['Principal Component 1', 'Principal Component 2'], 
           fontsize='small')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Project 2D data onto the first Principal Component The following figure demonstrates how the first principal component minimizes the reconstruction error, measured as the sum of the distances between the data points and the new axis.<jupyter_code># de-mean data, convert to numpy array
data_ = data.sub(data.mean())
X_ = data_.values
x_, y_ = X_.T
ax = pd.DataFrame({'$x_1$': x_, '$x_2$': y_}).plot.scatter(x='$x_1$', 
                                                           y='$x_2$', 
                                                           s=15, 
                                                           title='1D Projection', 
                                                           figsize=(6, 12))
ax.set_aspect('equal')

# plot first component
t = np.linspace(-25, 25, n_signals)
pc_x, pc_y = t * pc1
ax.plot(pc_x, pc_y, c='k', lw=1)

# project original data on first component
proj_x, proj_y = (X_.dot(pc1) * pc1.T).T
ax.scatter(proj_x, proj_y, s=15, c='k')

# plot link from data to projected points
lines_x, lines_y = np.c_[x_, proj_x], np.c_[y_, proj_y]
ax.plot(lines_x.T, lines_y.T, lw=1, c='darkgrey')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Plot 1D Representation<jupyter_code>projection1D = data_.dot(pc1)
ax = projection1D.rename(columns={0: '$z_1$'})\
    .assign(x2=0).plot.scatter(x='$z_1$', y='x2', s=10, title='1D Signal')
ax.get_yaxis().set_visible(False)
plt.tight_layout();<jupyter_output><empty_output><jupyter_text>Compare to Linear Regression Finally, the below figure contrasts PCA with supervised OLS, which approximates the outcome variable (here we choose x2) by a (one-dimensional) hyperplane computed from the (single) feature. The vertical lines highlight how OLS minimizes the distance along the outcome axis, in contrast with PCA that minimizes the distances orthogonal to the hyperplane.<jupyter_code>ax = pd.DataFrame({'$x_1$': x_, '$x_2$': y_}).plot.scatter(x='$x_1$', 
                                                           y='$x_2$', 
                                                           s=15, 
                                                           title='PCA vs OLS', 
                                                           figsize=(6, 12))
ax.set_aspect('equal')

# get OLS line
reg_X = np.column_stack((x_, np.ones_like(x_)))
(m, b), _, _, _ = lstsq(reg_X, y_)
reg_y = m * x_ + b
ax.plot(x_, reg_y, c='k')

# plot residuals
lines_x, lines_y = np.c_[x_, x_], np.c_[y_, reg_y]
ax.plot(lines_x.T, lines_y.T, lw=1)
plt.tight_layout()<jupyter_output>/home/stefan/.pyenv/versions/miniconda3-latest/envs/ml4t/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  # Remove the CWD from sys.path while we load stuff.<jupyter_text>Combined Figure The below figure is included in the book.<jupyter_code>fig, axes = plt.subplots(ncols=3, figsize=(14, 6), sharex=True, sharey=True)

title='Principal Component Vectors'
data.plot.scatter(x=0, y=1, s=15, ax=axes[0], title=title)
axes[0].set_aspect('equal')

origin_x, origin_y = pca.mean_
dx1, dy1 = np.squeeze(pc1.T) * l1
dx2, dy2 = np.squeeze(pc2.T) * l2
pc1_arrow = axes[0].arrow(origin_x, origin_y, dx1, dy1,  width=.2, color='k')
pc2_arrow = axes[0].arrow(origin_x, origin_y, dx2, dy2,  width=.2, color='r')
axes[0].legend([pc1_arrow, pc2_arrow],
               ['Principal Component 1', 'Principal Component 2'],
               fontsize='small')

# de-mean data, convert to numpy array
data_ = data.sub(data.mean())
X_ = data_.values
x_, y_ = X_.T
pd.DataFrame({'$x_1$': x_, '$x_2$': y_}).plot.scatter(x='$x_1$',
                                                      y='$x_2$',
                                                      s=15,
                                                      title='1D Projection',
                                                      ax=axes[1])
axes[1].set_aspect('equal')

# plot first component
t = np.linspace(-20, 20, n_signals)
pc_x, pc_y = t * pc1
axes[1].plot(pc_x, pc_y, c='k', lw=1)

# project original data on first component
proj_x, proj_y = (X_.dot(pc1) * pc1.T).T
axes[1].scatter(proj_x, proj_y, s=15, c='k')

# plot link from data to projected points
lines_x, lines_y = np.c_[x_, proj_x], np.c_[y_, proj_y]
axes[1].plot(lines_x.T, lines_y.T, lw=1, c='darkgrey')

pd.DataFrame({'$x_1$': x_, '$x_2$': y_}).plot.scatter(x='$x_1$',
                                                      y='$x_2$',
                                                      s=15,
                                                      title='PCA vs OLS',
                                                      ax=axes[2])
ax = axes[2].set_aspect('equal')

# get OLS line
reg_X = np.column_stack((x_, np.ones_like(x_)))
(m, b), _, _, _ = lstsq(reg_X, y_)
reg_y = m * x_ + b
ax = axes[2].plot(x_, reg_y, c='k')

# plot residuals
lines_x, lines_y = np.c_[x_, x_], np.c_[y_, reg_y]
ax = axes[2].plot(lines_x.T, lines_y.T, lw=1, c='darkgrey')

fig.tight_layout()
# fig.savefig('figures/pc_comparison', dpi=600);<jupyter_output>/home/stefan/.pyenv/versions/miniconda3-latest/envs/ml4t/lib/python3.7/site-packages/ipykernel_launcher.py:49: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.<jupyter_text>Recover Data using Inverse Transformation 1D => 2D<jupyter_code>recovered_data = projection1D.dot(pc1.T).rename(columns={0: '$x_1$', 1: '$x_2$'})

rms_reconstruction_error = np.sqrt(np.mean(np.sum(np.square(recovered_data-data_), axis=1)))

rss_data = np.sqrt(np.sum(data_.values**2))

relative_loss = rms_reconstruction_error / rss_data

title='Reconstructed Data | Error: {:.2%}'.format(relative_loss)
ax = recovered_data.plot.scatter(x=0, y=1, title=title, c='k')
ax.set_aspect('equal')
data_.plot.scatter(x=0, y=1, s=10, ax=ax, figsize=(4,8))
plt.legend(handles=[Patch(color='k', label='Recovered'),  
                    Patch(label='Original Data')])
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Projection and inverse transformation lead to the same result<jupyter_code>np.allclose(recovered_data, X_.dot(pc1) * pc1.T)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11827016520894072, 'fraction_numerical': 0.02118561710398445, 'mean_word_length': 3.466579861111111, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '81620', 'license': 'MIT', 'n_tokens_mistral': 3481, 'n_tokens_neox': 3316, 'n_words': 959, 'path': 'Chapter12/01_linear_dimensionality_reduction/01_pca_key_ideas.ipynb', 'repo_name': 'gabrielmahia/pesAI', 'size': '10290'}