<jupyter_start><jupyter_text>Chapter 8 - Dimensionality Reduction  Unsupervised LearningIn supervised learning, we have an access to $p$ features measured on $n$ observations, and a response $y$ is given. The goal is then to predict $y$ using the $p$ features.In unsupervised learning, we only have a set of features $X_1, \cdots, X_p$ measured on $n$ observations. We are not interested in prediction because we do not have an associated response variable $y$. Rather, the goal is to discover interesting things about the measurements $X_1, \cdots, X_p$. Can we visualise the data? Can we discover subgroups among the variables or the observations?Unsupervised learning is much more challenging. The analysis tends to be more subjective / biased and there is no simple goal of the analysis. Unsupervised learning is part of exploratory data analysis. Furthermore, in unsupervised learning there is no way to check our work - we don't have tools like cross-validation to measure the performance of our technique.  The Curse of DimensionalityMany ML problems involve training on many features, for each training instance - $p$ can be very large. This process is slow and makes it harder to find a good solution. This is called the curse of dimensionality.Consider the MNIST example. The pixels on the image borders are almost always white (feature has low variation) so they can be removed. Neighbouring pixels usually have the same colour so they can be averaged to form one feature (features have high correlation). Such steps do not result in much information loss.In theory, one solution to overcome the curse of dimsensionality is to increase the size of the training set. However, in reality, the number of training instances required to reach a given density ($\frac np$) grows exponentially with the number of dimensions.<jupyter_code>import pandas as pd
import numpy as np

from sklearn.model_selection import KFold
from sklearn.decomposition import PCA

import matplotlib.pyplot as plt

def load(fname):
    import pickle
    mnist = None
    try:
        with open(fname, 'rb') as f:
            mnist = pickle.load(f)
            return mnist
    except FileNotFoundError:
        from sklearn.datasets import fetch_openml
        mnist = fetch_openml('mnist_784', version=1, cache=True)
        with open(fname, 'wb') as f:
            mnist = pickle.dump(mnist, f)
        return mnist<jupyter_output><empty_output><jupyter_text>Principal Component Analysis (PCA)Consider an ML problem that has a large set of correlated variables (e.g. the neighbouring pixels example in the MNIST dataset). We can summarize this large set of correlated variables with a smaller number of representative variables using principal components.Say we want to visualise $n$ observations with $p$ features $X_1, \cdots, X_p$. We can visualise the data using $n \choose 2$ scatterplots. If $p$ is large then we cannot possibly look at all of them. Also, most of them will likely be uninformative as they contain only a small fraction of the total information / variance in the dataset. A better method is to visualise the $n$ observations when $p$ is large. Particularly, we want to find a low-dimensional representation of the data / reduce the dimensions of the data, capturing as much of the information as possible.PCA allows us to do so. The approach is to pick the hyperplane that preserves the most amount of variance, as it will likely lose the least amount of information compared to other projections. Each of these hyperplanes is a linear combination of the $p$ features.The first principal component of a set of features $X_1, \cdots, X_p$ is the normalised linear combination of the features:$$Z_1 = \phi_{11}X_1 + \phi_{21}X_2 + \cdots + \phi_{p1}X_p$$that has the largest variance. The elements $\phi_{j1} \forall j \in \{1\cdots p\}$ are the loadings of the first principal component and together, they make the principal component loading vector $\phi_1$. Mathematically, the first principal component loading vector has the loadings:$$\phi_1 = \begin{pmatrix}\phi_{11}&\phi_{21}&\cdots&\phi_{p1}\end{pmatrix}^T$$Normalised means that the sum of the loadings $\sum_{j=1}^p\phi_{j1}^2=1$. This constraint is needed as setting these elements to be arbitrarily large would results in an arbitrary large variance. To find the first principal components of a $n\times p$ training set $\mathbf X$, we first center the data to have mean zero. Then, we find the linear combination of the feature values: $$z_{i1} = \phi_{11}x_{i1} +  \phi_{21}x_{i2} + \cdots +  \phi_{p1}x_{ip}\,\, \forall i \in \{1,\cdots,n\}$$that has the largest sample variance subject to the constraint $\sum_{j=1}^p\phi_{j1}^2=1$. In other words, the first principal component loading vector solves the optimisation problem:$$\underset{\phi_{11}, \cdots, \phi_{p1}}{\text{Maximise }} \left\{\frac 1n \sum_{i=1}^n\begin{pmatrix}\sum_{j=1}^p\phi_{j1}x_{ij}\end{pmatrix}^2\right\} \text{ s. t. }$$$$\sum_{j=1}^p \phi_{j1}^2=1$$Since $z_{i1} = \phi_{11}x_{i1} +  \phi_{21}x_{i2} + \cdots +  \phi_{p1}x_{ip}$ we can simplify the optimisation problem to:$$\underset{\phi_{11}, \cdots, \phi_{p1}}{\text{Maximise }} \left\{\frac 1n \sum_{i=1}^nz_{i1}^2\right\} \text{ s. t. }$$ $$\sum_{j=1}^p \phi_{j1}^2=1$$Furthermore, since we have a zero-ed mean, that means $\frac 1n \sum_{i=1}^nx_{ij}=0$, the mean of $z_{11}, \cdots, z_{n1}$ is zero as well. Hence, the objective we are maximising is just the sample variance of the $n$ values of $z_{i1}$. We refer $z_{11}, \cdots, z_{n1}$ as the scores of the first principal component.Solving the optimisation problem involves eigenvalue decomposition. In particular, there is a standard matrix factorization technique called Singlular Value Decomposition (SVD) that decomposes the training set matrix $\mathbf X$ to the dot product of three matrices:$$\mathbf X = \mathbf U \cdot \Sigma \cdot \mathbf V^T$$ where $\mathbf V^T$ contains all the principal components that we are looking for. Interpretation: The loadings of the first principal component, $\phi_1$ is the direction in feature space along which the data varies the most. If we project the $n$ training samples onto this direction, the projected values are the principal component scores $z_{11}, \cdots, z_{n1}$ themselves and they will lose the least amount of information compared to other projections. PCA identifies the axis that accounts for the largest amount of variance in the training set.In this example, the observations are in 2D. The first principal component loading vector is the green line. $\phi_1 = (\phi_{11}, \phi_{21}) = (0.839, 0.544)$<jupyter_code># Ingest
mnsit = load('mnist.data.pkl')
mnsit_X, mnsit_y = mnsit['data'], mnsit['target']

kf = KFold(n_splits=10, shuffle=True, random_state=0)
kf.get_n_splits()
t1 = []
for train_index, test_index in kf.split(mnsit_X, mnsit_y):
    t1 = test_index
    break
X = mnsit_X[t1]
y_test = mnsit_y[t1]
# For testing
# print(pd.Series(y_test).value_counts())<jupyter_output><empty_output><jupyter_text>The following is an implementation of PCA using the SVD method from `numpy`.<jupyter_code># PCA using SVD Algorithm

# Center the data
X_centered = X - X.mean(axis=0)

# SVD algorithm to obtain the loadings
u, s, Vt = np.linalg.svd(X_centered)

# Obtain the principal components
W = Vt.T[:,:196]
X_d1 = X_centered.dot(W)
# For testing
# print(X_test.shape)
# print(X_test_centered.shape)
# X_test_centered[:2]<jupyter_output><empty_output><jupyter_text>After the first principal component $Z_1$ of the features are determined, we can find the second principal component $Z_2$. The second principal component is the linear combination $X_1, \cdots, X_p$ that has maximal variance out of all linear combinations that are uncorrelated with $Z_1$. The second principal component scores $z_{12}, \cdots, z_{n2}$ take the form:$$z_{i2} = \phi_{12}x_{i1} +  \phi_{22}x_{i2} + \cdots +  \phi_{p2}x_{ip}$$where $\phi_{2}$ is the second principal component loading vector, with elements $\phi_{12}, \phi_{22}, \cdots, \phi_{p2}$. Note that this loading vector is constrained such that the direction must be orthogonal (perpendicular) to the direction of $\phi_1$. In 3D space, once we have found $\phi_1$, there is only one possibility for $\phi_2$, which is the blue dashed line.But in a larger dataset with $p>2$ variables, there are multiple candidates for principal components, and they are defined in a similar manner. To find $\phi_2$, we solve the same maximisation problem, but with the additional constraint that $\phi_2$ is orthogonal to $\phi_1$.Once all the principal components are identified, you can reduce the dimensionality of the dataset 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. To do so, simply compute the dot product of the training sest matrix $\mathbf X$ by the matrix $\mathbf W_d$.$$\mathbf X_{d\text{-proj}} = \mathbf X \cdot \mathbf W_d$$ The following is the `sklearn` implementation of PCA.<jupyter_code>pca = PCA(n_components=196, svd_solver='full')
X_d2 = pca.fit_transform(X)
# For testing
# Validate that both numpy and sklearn approaches yield the same results
# i_X = 19
# Xv_d1 = X_d1[i_X]
# Xv_d2 = X_d2[i_X]
# for i,j in zip(Xv_d1, Xv_d2):
#     print('{:2f}'.format(abs(float(abs(j) - abs(i)))))
# For testing
# Validate that the loadings for both numpy and sklearn approaches yield the same results
# print(pca.components_[100][20:50])
# print(Vt[100][20:50])<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup