mhjensen
1245 lines
45 KiB
Plaintext
1245 lines
45 KiB
Plaintext
TITLE: Data Analysis and Machine Learning: Preprocessing and Dimensionality Reduction
|
||
AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
|
||
DATE: today
|
||
|
||
|
||
!split
|
||
===== Reducing the number of degrees of freedom, overarching view =====
|
||
!bblock
|
||
|
||
Many Machine Learning problems involve thousands or even millions of
|
||
features for each training instance. Not only does this make training
|
||
extremely slow, it can also make it much harder to find a good
|
||
solution, as we will see. This problem is often referred to as the
|
||
curse of dimensionality. Fortunately, in real-world problems, it is
|
||
often possible to reduce the number of features considerably, turning
|
||
an intractable problem into a tractable one.
|
||
|
||
Here we will discuss some of the most popular dimensionality reduction
|
||
techniques: the principal component analysis (PCA), Kernel PCA, and
|
||
Locally Linear Embedding (LLE). Furthermore, we will start by looking
|
||
at some simple preprocessing of the data which allow us to rescale the
|
||
data.
|
||
|
||
Principal component analysis and its various variants deal with the
|
||
problem of fitting a low-dimensional "affine
|
||
subspace":"https://en.wikipedia.org/wiki/Affine_space" to a set of of
|
||
data points in a high-dimensional space. With its family of methods it
|
||
is one of the most used tools in data modeling, compression and
|
||
visualization.
|
||
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Preprocessing our data =====
|
||
!bblock
|
||
|
||
Before we proceed however, we will discuss how to preprocess our
|
||
data. Till now and in connection with our previous examples we have
|
||
not met so many cases where we are too sensitive to the scaling of our
|
||
data. Normally the data may need a rescaling and/or may be sensitive
|
||
to extreme values. Scaling the data renders our inputs much more
|
||
suitable for the algorithms we want to employ.
|
||
|
||
_Scikit-Learn_ has several functions which allow us to rescale the
|
||
data, normally resulting in much better results in terms of various
|
||
accuracy scores. The _StandardScaler_ function in _Scikit-Learn_
|
||
ensures that for each feature/predictor we study the mean value is
|
||
zero and the variance is one (every column in the design/feature
|
||
matrix). This scaling has the drawback that it does not ensure that
|
||
we have a particular maximum or minimum in our data set. Another
|
||
function included in _Scikit-Learn_ is the _MinMaxScaler_ which
|
||
ensures that all features are exactly between $0$ and $1$. The
|
||
|
||
!split
|
||
===== More preprocessing =====
|
||
|
||
|
||
The _Normalizer_ scales each data
|
||
point such that the feature vector has a euclidean length of one. In other words, it
|
||
projects a data point on the circle (or sphere in the case of higher dimensions) with a
|
||
radius of 1. This means every data point is scaled by a different number (by the
|
||
inverse of it’s length).
|
||
This normalization is often used when only the direction (or angle) of the data matters,
|
||
not the length of the feature vector.
|
||
|
||
The _RobustScaler_ works similarly to the StandardScaler in that it
|
||
ensures statistical properties for each feature that guarantee that
|
||
they are on the same scale. However, the RobustScaler uses the median
|
||
and quartiles, instead of mean and variance. This makes the
|
||
RobustScaler ignore data points that are very different from the rest
|
||
(like measurement errors). These odd data points are also called
|
||
outliers, and might often lead to trouble for other scaling
|
||
techniques.
|
||
|
||
!eblock
|
||
|
||
!split
|
||
===== Simple preprocessing examples, Franke function and regression =====
|
||
|
||
!bc pycod
|
||
# Common imports
|
||
import os
|
||
import numpy as np
|
||
import pandas as pd
|
||
import matplotlib.pyplot as plt
|
||
import sklearn.linear_model as skl
|
||
from sklearn.metrics import mean_squared_error
|
||
from sklearn.model_selection import train_test_split
|
||
from sklearn.preprocessing import MinMaxScaler, StandardScaler, Normalizer
|
||
from sklearn.svm import SVR
|
||
|
||
# Where to save the figures and data files
|
||
PROJECT_ROOT_DIR = "Results"
|
||
FIGURE_ID = "Results/FigureFiles"
|
||
DATA_ID = "DataFiles/"
|
||
|
||
if not os.path.exists(PROJECT_ROOT_DIR):
|
||
os.mkdir(PROJECT_ROOT_DIR)
|
||
|
||
if not os.path.exists(FIGURE_ID):
|
||
os.makedirs(FIGURE_ID)
|
||
|
||
if not os.path.exists(DATA_ID):
|
||
os.makedirs(DATA_ID)
|
||
|
||
def image_path(fig_id):
|
||
return os.path.join(FIGURE_ID, fig_id)
|
||
|
||
def data_path(dat_id):
|
||
return os.path.join(DATA_ID, dat_id)
|
||
|
||
def save_fig(fig_id):
|
||
plt.savefig(image_path(fig_id) + ".png", format='png')
|
||
|
||
|
||
def FrankeFunction(x,y):
|
||
term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))
|
||
term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))
|
||
term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))
|
||
term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)
|
||
return term1 + term2 + term3 + term4
|
||
|
||
|
||
def create_X(x, y, n ):
|
||
if len(x.shape) > 1:
|
||
x = np.ravel(x)
|
||
y = np.ravel(y)
|
||
|
||
N = len(x)
|
||
l = int((n+1)*(n+2)/2) # Number of elements in beta
|
||
X = np.ones((N,l))
|
||
|
||
for i in range(1,n+1):
|
||
q = int((i)*(i+1)/2)
|
||
for k in range(i+1):
|
||
X[:,q+k] = (x**(i-k))*(y**k)
|
||
|
||
return X
|
||
|
||
|
||
# Making meshgrid of datapoints and compute Franke's function
|
||
n = 5
|
||
N = 1000
|
||
x = np.sort(np.random.uniform(0, 1, N))
|
||
y = np.sort(np.random.uniform(0, 1, N))
|
||
z = FrankeFunction(x, y)
|
||
X = create_X(x, y, n=n)
|
||
# split in training and test data
|
||
X_train, X_test, y_train, y_test = train_test_split(X,z,test_size=0.2)
|
||
|
||
|
||
svm = SVR(gamma='auto',C=10.0)
|
||
svm.fit(X_train, y_train)
|
||
|
||
# The mean squared error and R2 score
|
||
print("MSE before scaling: {:.2f}".format(mean_squared_error(svm.predict(X_test), y_test)))
|
||
print("R2 score before scaling {:.2f}".format(svm.score(X_test,y_test)))
|
||
|
||
scaler = StandardScaler()
|
||
scaler.fit(X_train)
|
||
X_train_scaled = scaler.transform(X_train)
|
||
X_test_scaled = scaler.transform(X_test)
|
||
|
||
print("Feature min values before scaling:\n {}".format(X_train.min(axis=0)))
|
||
print("Feature max values before scaling:\n {}".format(X_train.max(axis=0)))
|
||
|
||
print("Feature min values after scaling:\n {}".format(X_train_scaled.min(axis=0)))
|
||
print("Feature max values after scaling:\n {}".format(X_train_scaled.max(axis=0)))
|
||
|
||
svm = SVR(gamma='auto',C=10.0)
|
||
svm.fit(X_train_scaled, y_train)
|
||
|
||
print("MSE after scaling: {:.2f}".format(mean_squared_error(svm.predict(X_test_scaled), y_test)))
|
||
print("R2 score for scaled data: {:.2f}".format(svm.score(X_test_scaled,y_test)))
|
||
|
||
!ec
|
||
|
||
|
||
|
||
!split
|
||
===== Simple preprocessing examples, breast cancer data and classification, Support Vector Machines =====
|
||
|
||
We show here how we can use a simple regression case on the breast
|
||
cancer data using support vector machines (SVM) as algorithm for
|
||
classification.
|
||
|
||
|
||
!bc pycod
|
||
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.svm import SVC
|
||
cancer = load_breast_cancer()
|
||
|
||
X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
|
||
print(X_train.shape)
|
||
print(X_test.shape)
|
||
|
||
svm = SVC(C=100)
|
||
svm.fit(X_train, y_train)
|
||
print("Test set accuracy: {:.2f}".format(svm.score(X_test,y_test)))
|
||
|
||
from sklearn.preprocessing import MinMaxScaler, StandardScaler
|
||
scaler = MinMaxScaler()
|
||
scaler.fit(X_train)
|
||
X_train_scaled = scaler.transform(X_train)
|
||
X_test_scaled = scaler.transform(X_test)
|
||
|
||
print("Feature min values before scaling:\n {}".format(X_train.min(axis=0)))
|
||
print("Feature max values before scaling:\n {}".format(X_train.max(axis=0)))
|
||
|
||
print("Feature min values before scaling:\n {}".format(X_train_scaled.min(axis=0)))
|
||
print("Feature max values before scaling:\n {}".format(X_train_scaled.max(axis=0)))
|
||
|
||
|
||
svm.fit(X_train_scaled, y_train)
|
||
print("Test set accuracy scaled data with Min-Max scaling: {:.2f}".format(svm.score(X_test_scaled,y_test)))
|
||
|
||
scaler = StandardScaler()
|
||
scaler.fit(X_train)
|
||
X_train_scaled = scaler.transform(X_train)
|
||
X_test_scaled = scaler.transform(X_test)
|
||
|
||
svm.fit(X_train_scaled, y_train)
|
||
print("Test set accuracy scaled data with Standar Scaler: {:.2f}".format(svm.score(X_test_scaled,y_test)))
|
||
|
||
!ec
|
||
|
||
!split
|
||
===== More on Cancer Data, now with Logistic Regression =====
|
||
|
||
|
||
!bc pycod
|
||
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()
|
||
|
||
# Set up training data
|
||
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("Test set accuracy: {:.2f}".format(logreg.score(X_test,y_test)))
|
||
|
||
# Scale 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)
|
||
logreg.fit(X_train_scaled, y_train)
|
||
print("Test set accuracy scaled data: {:.2f}".format(logreg.score(X_test_scaled,y_test)))
|
||
|
||
!ec
|
||
|
||
|
||
|
||
|
||
!split
|
||
===== Why should we think of reducing the dimensionality =====
|
||
|
||
In addition to the plot of the features, we study now also the covariance (or rather the correlation matrix).
|
||
We use also _Pandas_ to compute the correlation matrix.
|
||
!bc pycod
|
||
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()
|
||
import pandas as pd
|
||
# Making a data frame
|
||
cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names)
|
||
|
||
fig, axes = plt.subplots(15,2,figsize=(10,20))
|
||
malignant = cancer.data[cancer.target == 0]
|
||
benign = cancer.data[cancer.target == 1]
|
||
ax = axes.ravel()
|
||
|
||
for i in range(30):
|
||
_, bins = np.histogram(cancer.data[:,i], bins =50)
|
||
ax[i].hist(malignant[:,i], bins = bins, alpha = 0.5)
|
||
ax[i].hist(benign[:,i], bins = bins, alpha = 0.5)
|
||
ax[i].set_title(cancer.feature_names[i])
|
||
ax[i].set_yticks(())
|
||
ax[0].set_xlabel("Feature magnitude")
|
||
ax[0].set_ylabel("Frequency")
|
||
ax[0].legend(["Malignant", "Benign"], loc ="best")
|
||
fig.tight_layout()
|
||
plt.show()
|
||
|
||
import seaborn as sns
|
||
correlation_matrix = cancerpd.corr().round(1)
|
||
# use the heatmap function from seaborn to plot the correlation matrix
|
||
# annot = True to print the values inside the square
|
||
sns.heatmap(data=correlation_matrix, annot=True)
|
||
plt.show()
|
||
|
||
#print eigvalues of correlation matrix
|
||
EigValues, EigVectors = np.linalg.eig(correlation_matrix)
|
||
print(EigValues)
|
||
!ec
|
||
|
||
In the above example we note two things. In the first plot we display
|
||
the overlap of benign and malignant tumors as functions of the various
|
||
features in the Wisconsing breast cancer data set. We see that for
|
||
some of the features we can distinguish clearly the benign and
|
||
malignant cases while for other features we cannot. This can point to
|
||
us which features may be of greater interest when we wish to classify
|
||
a benign or not benign tumour.
|
||
|
||
In the second figure we have computed the so-called correlation
|
||
matrix, which in our case with thirty features becomes a $30\times 30$
|
||
matrix.
|
||
|
||
We constructed this matrix using _pandas_ via the statements
|
||
!bc pycod
|
||
cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names)
|
||
!ec
|
||
and then
|
||
!bc pycod
|
||
correlation_matrix = cancerpd.corr().round(1)
|
||
!ec
|
||
|
||
Diagonalizing this matrix we can in turn say something about which
|
||
features are of relevance and which are not. But before we proceed we
|
||
need to define covariance and correlation matrices. This leads us to
|
||
the classical Principal Component Analysis (PCA) theorem with
|
||
applications.
|
||
|
||
|
||
|
||
!split
|
||
===== Basic ideas of the Principal Component Analysis (PCA) =====
|
||
|
||
The principal component analysis deals with the problem of fitting a
|
||
low-dimensional affine subspace $S$ of dimension $d$ much smaller than
|
||
the totaldimension $D$ of the problem at hand (our data
|
||
set). Mathematically it can be formulated as a statistical problem or
|
||
a geometric problem. In our discussion of the theorem for the
|
||
classical PCA, we will stay with a statistical approach. This is also
|
||
what set the scene historically which for the PCA.
|
||
|
||
We have a data set defined by a design/feature matrix $\bm{X}$ (see below for its definition)
|
||
* Each data point is determined by $p$ extrinsic (measurement) variables
|
||
* We may want to ask the following question: Are there fewer intrinsic variables (say $d << p$) that still approximately describe the data?
|
||
* If so, these intrinsic variables may tell us something important and finding these intrinsic variables is what dimension reduction methods do.
|
||
|
||
|
||
!split
|
||
===== Introducing the Covariance and Correlation functions =====
|
||
|
||
Before we discuss the PCA theorem, we need to remind ourselves about
|
||
the definition of the covariance and the correlation function. These are quantities
|
||
|
||
Suppose we have defined two vectors
|
||
$\hat{x}$ and $\hat{y}$ with $n$ elements each. The covariance matrix $\bm{C}$ is defined as
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x},\bm{y}] = \begin{bmatrix} \mathrm{cov}[\bm{x},\bm{x}] & \mathrm{cov}[\bm{x},\bm{y}] \\
|
||
\mathrm{cov}[\bm{y},\bm{x}] & \mathrm{cov}[\bm{y},\bm{y}] \\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
where for example
|
||
!bt
|
||
\[
|
||
\mathrm{cov}[\bm{x},\bm{y}] =\frac{1}{n} \sum_{i=0}^{n-1}(x_i- \overline{x})(y_i- \overline{y}).
|
||
\]
|
||
!et
|
||
With this definition and recalling that the variance is defined as
|
||
!bt
|
||
\[
|
||
\mathrm{var}[\bm{x}]=\frac{1}{n} \sum_{i=0}^{n-1}(x_i- \overline{x})^2,
|
||
\]
|
||
!et
|
||
we can rewrite the covariance matrix as
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x},\bm{y}] = \begin{bmatrix} \mathrm{var}[\bm{x}] & \mathrm{cov}[\bm{x},\bm{y}] \\
|
||
\mathrm{cov}[\bm{x},\bm{y}] & \mathrm{var}[\bm{y}] \\
|
||
\end{bmatrix}.
|
||
\]
|
||
!et
|
||
|
||
The covariance takes values between zero and infinity and may thus
|
||
lead to problems with loss of numerical precision for particularly
|
||
large values. It is common to scale the covariance matrix by
|
||
introducing instead the correlation matrix defined via the so-called
|
||
correlation function
|
||
|
||
!bt
|
||
\[
|
||
\mathrm{corr}[\bm{x},\bm{y}]=\frac{\mathrm{cov}[\bm{x},\bm{y}]}{\sqrt{\mathrm{var}[\bm{x}] \mathrm{var}[\bm{y}]}}.
|
||
\]
|
||
!et
|
||
|
||
The correlation function is then given by values $\mathrm{corr}[\bm{x},\bm{y}]
|
||
\in [-1,1]$. This avoids eventual problems with too large values. We
|
||
can then define the correlation matrix for the two vectors $\bm{x}$
|
||
and $\bm{y}$ as
|
||
|
||
!bt
|
||
\[
|
||
\bm{K}[\bm{x},\bm{y}] = \begin{bmatrix} 1 & \mathrm{corr}[\bm{x},\bm{y}] \\
|
||
\mathrm{corr}[\bm{y},\bm{x}] & 1 \\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
|
||
In the above example this is the function we constructed using _pandas_.
|
||
|
||
!split
|
||
===== Correlation Function and Design/Feature Matrix =====
|
||
|
||
In our derivation of the various regression algorithms like _Ordinary Least Squares_ or _Ridge regression_
|
||
we defined the design/feature matrix $\bm{X}$ as
|
||
|
||
!bt
|
||
\[
|
||
\bm{X}=\begin{bmatrix}
|
||
x_{0,0} & x_{0,1} & x_{0,2}& \dots & \dots x_{0,p-1}\\
|
||
x_{1,0} & x_{1,1} & x_{1,2}& \dots & \dots x_{1,p-1}\\
|
||
x_{2,0} & x_{2,1} & x_{2,2}& \dots & \dots x_{2,p-1}\\
|
||
\dots & \dots & \dots & \dots \dots & \dots \\
|
||
x_{n-2,0} & x_{n-2,1} & x_{n-2,2}& \dots & \dots x_{n-2,p-1}\\
|
||
x_{n-1,0} & x_{n-1,1} & x_{n-1,2}& \dots & \dots x_{n-1,p-1}\\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
with $\bm{X}\in {\mathbb{R}}^{n\times p}$, with the predictors/features $p$ refering to the column numbers and the
|
||
entries $n$ being the row elements.
|
||
We can rewrite the design/feature matrix in terms of its column vectors as
|
||
!bt
|
||
\[
|
||
\bm{X}=\begin{bmatrix} \bm{x}_0 & \bm{x}_1 & \bm{x}_2 & \dots & \dots & \bm{x}_{p-1}\end{bmatrix},
|
||
\]
|
||
!et
|
||
with a given vector
|
||
!bt
|
||
\[
|
||
\bm{x}_i^T = \begin{bmatrix}x_{0,i} & x_{1,i} & x_{2,i}& \dots & \dots x_{n-1,i}\end{bmatrix}.
|
||
\]
|
||
!et
|
||
|
||
With these definitions, we can now rewrite our $2\times 2$
|
||
correaltion/covariance matrix in terms of a moe general design/feature
|
||
matrix $\bm{X}\in {\mathbb{R}}^{n\times p}$. This leads to a $p\times p$
|
||
covariance matrix for the vectors $\bm{x}_i$ with $i=0,1,\dots,p-1$
|
||
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x}] = \begin{bmatrix}
|
||
\mathrm{var}[\bm{x}_0] & \mathrm{cov}[\bm{x}_0,\bm{x}_1] & \mathrm{cov}[\bm{x}_0,\bm{x}_2] & \dots & \dots & \mathrm{cov}[\bm{x}_0,\bm{x}_{p-1}]\\
|
||
\mathrm{cov}[\bm{x}_1,\bm{x}_0] & \mathrm{var}[\bm{x}_1] & \mathrm{cov}[\bm{x}_1,\bm{x}_2] & \dots & \dots & \mathrm{cov}[\bm{x}_1,\bm{x}_{p-1}]\\
|
||
\mathrm{cov}[\bm{x}_2,\bm{x}_0] & \mathrm{cov}[\bm{x}_2,\bm{x}_1] & \mathrm{var}[\bm{x}_2] & \dots & \dots & \mathrm{cov}[\bm{x}_2,\bm{x}_{p-1}]\\
|
||
\dots & \dots & \dots & \dots & \dots & \dots \\
|
||
\dots & \dots & \dots & \dots & \dots & \dots \\
|
||
\mathrm{cov}[\bm{x}_{p-1},\bm{x}_0] & \mathrm{cov}[\bm{x}_{p-1},\bm{x}_1] & \mathrm{cov}[\bm{x}_{p-1},\bm{x}_{2}] & \dots & \dots & \mathrm{var}[\bm{x}_{p-1}]\\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
and the correlation matrix
|
||
!bt
|
||
\[
|
||
\bm{K}[\bm{x}] = \begin{bmatrix}
|
||
1 & \mathrm{corr}[\bm{x}_0,\bm{x}_1] & \mathrm{corr}[\bm{x}_0,\bm{x}_2] & \dots & \dots & \mathrm{corr}[\bm{x}_0,\bm{x}_{p-1}]\\
|
||
\mathrm{corr}[\bm{x}_1,\bm{x}_0] & 1 & \mathrm{corr}[\bm{x}_1,\bm{x}_2] & \dots & \dots & \mathrm{corr}[\bm{x}_1,\bm{x}_{p-1}]\\
|
||
\mathrm{corr}[\bm{x}_2,\bm{x}_0] & \mathrm{corr}[\bm{x}_2,\bm{x}_1] & 1 & \dots & \dots & \mathrm{corr}[\bm{x}_2,\bm{x}_{p-1}]\\
|
||
\dots & \dots & \dots & \dots & \dots & \dots \\
|
||
\dots & \dots & \dots & \dots & \dots & \dots \\
|
||
\mathrm{corr}[\bm{x}_{p-1},\bm{x}_0] & \mathrm{corr}[\bm{x}_{p-1},\bm{x}_1] & \mathrm{corr}[\bm{x}_{p-1},\bm{x}_{2}] & \dots & \dots & 1\\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
|
||
|
||
!split
|
||
===== Covariance Matrix Examples =====
|
||
|
||
|
||
The Numpy function _np.cov_ calculates the covariance elements using
|
||
the factor $1/(n-1)$ instead of $1/n$ since it assumes we do not have
|
||
the exact mean values. The following simple function uses the
|
||
_np.vstack_ function which takes each vector of dimension $1\times n$
|
||
and produces a $2\times n$ matrix $\bm{W}$
|
||
|
||
|
||
!bt
|
||
\[
|
||
\bm{W} = \begin{bmatrix} x_0 & y_0 \\
|
||
x_1 & y_1 \\
|
||
x_2 & y_2\\
|
||
\dots & \dots \\
|
||
x_{n-2} & y_{n-2}\\
|
||
x_{n-1} & y_{n-1} &
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
|
||
which in turn is converted into into the $2\times 2$ covariance matrix
|
||
$\bm{C}$ via the Numpy function _np.cov()_. We note that we can also calculate
|
||
the mean value of each set of samples $\bm{x}$ etc using the Numpy
|
||
function _np.mean(x)_. We can also extract the eigenvalues of the
|
||
covariance matrix through the _np.linalg.eig()_ function.
|
||
|
||
!bc pycod
|
||
# Importing various packages
|
||
import numpy as np
|
||
n = 100
|
||
x = np.random.normal(size=n)
|
||
print(np.mean(x))
|
||
y = 4+3*x+np.random.normal(size=n)
|
||
print(np.mean(y))
|
||
W = np.vstack((x, y))
|
||
C = np.cov(W)
|
||
print(C)
|
||
!ec
|
||
|
||
!split
|
||
===== Correlation Matrix =====
|
||
|
||
The previous example can be converted into the correlation matrix by
|
||
simply scaling the matrix elements with the variances. We should also
|
||
subtract the mean values for each column. This leads to the following
|
||
code which sets up the correlations matrix for the previous example in
|
||
a more brute force way. Here we scale the mean values for each column of the design matrix, calculate the relevant mean values and variances and then finally set up the $2\times 2$ correlation matrix (since we have only two vectors).
|
||
|
||
!bc pycod
|
||
import numpy as np
|
||
n = 100
|
||
# define two vectors
|
||
x = np.random.random(size=n)
|
||
y = 4+3*x+np.random.normal(size=n)
|
||
#scaling the x and y vectors
|
||
x = x - np.mean(x)
|
||
y = y - np.mean(y)
|
||
variance_x = np.sum(x@x)/n
|
||
variance_y = np.sum(y@y)/n
|
||
print(variance_x)
|
||
print(variance_y)
|
||
cov_xy = np.sum(x@y)/n
|
||
cov_xx = np.sum(x@x)/n
|
||
cov_yy = np.sum(y@y)/n
|
||
C = np.zeros((2,2))
|
||
C[0,0]= cov_xx/variance_x
|
||
C[1,1]= cov_yy/variance_y
|
||
C[0,1]= cov_xy/np.sqrt(variance_y*variance_x)
|
||
C[1,0]= C[0,1]
|
||
print(C)
|
||
!ec
|
||
|
||
We see that the matrix elements along the diagonal are one as they
|
||
should be and that the matrix is symmetric. Furthermore, diagonalizing
|
||
this matrix we easily see that it is a positive definite matrix.
|
||
|
||
The above procedure with _numpy_ can be made more compact if we use _pandas_.
|
||
|
||
!split
|
||
===== Correlation Matrix with Pandas =====
|
||
|
||
We whow here how we can set up the correlation matrix using _pandas_, as done in this simple code
|
||
!bc pycod
|
||
import numpy as np
|
||
import pandas as pd
|
||
n = 10
|
||
x = np.random.normal(size=n)
|
||
x = x - np.mean(x)
|
||
y = 4+3*x+np.random.normal(size=n)
|
||
y = y - np.mean(y)
|
||
X = (np.vstack((x, y))).T
|
||
print(X)
|
||
Xpd = pd.DataFrame(X)
|
||
print(Xpd)
|
||
correlation_matrix = Xpd.corr()
|
||
print(correlation_matrix)
|
||
!ec
|
||
|
||
|
||
We expand this model to the Franke function discussed above.
|
||
|
||
!split
|
||
===== Correlation Matrix with Pandas and the Franke function =====
|
||
|
||
!bc pycod
|
||
# Common imports
|
||
import numpy as np
|
||
import pandas as pd
|
||
|
||
|
||
def FrankeFunction(x,y):
|
||
term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))
|
||
term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))
|
||
term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))
|
||
term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)
|
||
return term1 + term2 + term3 + term4
|
||
|
||
|
||
def create_X(x, y, n ):
|
||
if len(x.shape) > 1:
|
||
x = np.ravel(x)
|
||
y = np.ravel(y)
|
||
|
||
N = len(x)
|
||
l = int((n+1)*(n+2)/2) # Number of elements in beta
|
||
X = np.ones((N,l))
|
||
|
||
for i in range(1,n+1):
|
||
q = int((i)*(i+1)/2)
|
||
for k in range(i+1):
|
||
X[:,q+k] = (x**(i-k))*(y**k)
|
||
|
||
return X
|
||
|
||
|
||
# Making meshgrid of datapoints and compute Franke's function
|
||
n = 4
|
||
N = 100
|
||
x = np.sort(np.random.uniform(0, 1, N))
|
||
y = np.sort(np.random.uniform(0, 1, N))
|
||
z = FrankeFunction(x, y)
|
||
X = create_X(x, y, n=n)
|
||
|
||
Xpd = pd.DataFrame(X)
|
||
# subtract the mean values and set up the covariance matrix
|
||
Xpd = Xpd - Xpd.mean()
|
||
covariance_matrix = Xpd.cov()
|
||
print(covariance_matrix)
|
||
!ec
|
||
|
||
We note here that the covariance is zero for the first rows and
|
||
columns since all matrix elements in the design matrix were set to one
|
||
(we are fitting the function in terms of a polynomial of degree $n$).
|
||
|
||
This means that the variance for these elements will be zero and will
|
||
cause problems when we set up the correlation matrix. We can simply
|
||
drop these elements and construct a correlation
|
||
matrix without these elements.
|
||
|
||
|
||
!split
|
||
===== Rewriting the Covariance and/or Correlation Matrix =====
|
||
|
||
We can rewrite the covariance matrix in a more compact form in terms of the design/feature matrix $\bm{X}$ as
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x}] = \frac{1}{n}\bm{X}\bm{X}^T= \mathbb{E}[\bm{X}\bm{X}^T].
|
||
\]
|
||
!et
|
||
|
||
To see this let us simply look at a design matrix $\bm{X}\in {\mathbb{R}}^{2\times 2}$
|
||
!bt
|
||
\[
|
||
\bm{X}=\begin{bmatrix}
|
||
x_{00} & x_{01}\\
|
||
x_{10} & x_{11}\\
|
||
\end{bmatrix}=\begin{bmatrix}
|
||
\bm{x}_{0} & \bm{x}_{1}\\
|
||
\end{bmatrix}.
|
||
\]
|
||
!et
|
||
|
||
If we then compute the expectation value
|
||
!bt
|
||
\[
|
||
\mathbb{E}[\bm{X}\bm{X}^T] = \frac{1}{n}\bm{X}\bm{X}^T=\begin{bmatrix}
|
||
x_{00}^2+x_{01}^2 & x_{00}x_{10}+x_{01}x_{11}\\
|
||
x_{10}x_{00}+x_{11}x_{01} & x_{10}^2+x_{11}^2\\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
which is just
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x}_0,\bm{x}_1] = \bm{C}[\bm{x}]=\begin{bmatrix} \mathrm{var}[\bm{x}_0] & \mathrm{cov}[\bm{x}_0,\bm{x}_1] \\
|
||
\mathrm{cov}[\bm{x}_1,\bm{x}_0] & \mathrm{var}[\bm{x}_1] \\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
where we wrote $$\bm{C}[\bm{x}_0,\bm{x}_1] = \bm{C}[\bm{x}]$$ to indicate that this the covariance of the vectors $\bm{x}$ of the design/feature matrix $\bm{X}$.
|
||
|
||
It is easy to generalize this to a matrix $\bm{X}\in {\mathbb{R}}^{n\times p}$.
|
||
|
||
|
||
!split
|
||
===== Towards the PCA theorem =====
|
||
|
||
We have that the covariance matrix (the correlation matrix involves a simple rescaling) is given as
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x}] = \frac{1}{n}\bm{X}\bm{X}^T= \mathbb{E}[\bm{X}\bm{X}^T].
|
||
\]
|
||
!et
|
||
Let us now assume that we can perform a series of orthogonal transformations where we employ some orthogonal matrices $\bm{S}$.
|
||
These matrices are defined as $\bm{S}\in {\mathbb{R}}^{p\times p}$ and obey the orthogonality requirements $\bm{S}\bm{S}^T=\bm{S}^T\bm{S}=\bm{I}$. The matrix can be written out in terms of the column vectors $\bm{s}_i$ as $\bm{S}=[\bm{s}_0,\bm{s}_1,\dots,\bm{s}_{p-1}]$ and $\bm{s}_i \in {\mathbb{R}}^{p}$.
|
||
|
||
Assume also that there is a transformation $\bm{S}\bm{C}[\bm{x}]\bm{S}^T=\bm{C}[\bm{y}]$ such that the new matrix $\bm{C}[\bm{y}]$ is diagonal with elements $[\lambda_0,\lambda_1,\lambda_2,\dots,\lambda_{p-1}]$.
|
||
|
||
That is we have
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{y}] = \mathbb{E}[\bm{S}\bm{X}\bm{X}^T\bm{S}^T]=\bm{S}\bm{C}[\bm{x}]\bm{S}^T,
|
||
\]
|
||
!et
|
||
since the matrix $\bm{S}$ is not a data dependent matrix. Multiplying with $\bm{S}^T$ from the left we have
|
||
!bt
|
||
\[
|
||
\bm{S}^T\bm{C}[\bm{y}] = \bm{C}[\bm{x}]\bm{S}^T,
|
||
\]
|
||
!et
|
||
and since $\bm{C}[\bm{y}]$ is diagonal we have for a given eigenvalue $i$ of the covariance matrix that
|
||
|
||
!bt
|
||
\[
|
||
\bm{S}^T_i\lambda_i = \bm{C}[\bm{x}]\bm{S}^T_i.
|
||
\]
|
||
!et
|
||
|
||
In the derivation of the PCA theorem we will assume that the eigenvalues are ordered in descending order, that is
|
||
$\lambda_0 > \lambda_1 > \dots > \lambda_{p-1}$.
|
||
|
||
|
||
The eigenvalues tell us then how much we need to stretch the
|
||
corresponding eigenvectors. Dimensions with large eigenvalues have
|
||
thus large variations (large variance) and define therefore useful
|
||
dimensions. The data points are more spread out in the direction of
|
||
these eigenvectors. Smaller eigenvalues mean on the other hand that
|
||
the corresponding eigenvectors are shrunk accordingly and the data
|
||
points are tightly bunched together and there is not much variation in
|
||
these specific directions. Hopefully then we could leave it out
|
||
dimensions where the eigenvalues are very small. If $p$ is very large,
|
||
we could then aim at reducing $p$ to $l << p$ and handle only $l$
|
||
features/predictors.
|
||
|
||
!split
|
||
===== The Algorithm before the Theorem =====
|
||
|
||
Here's how we would proceed in setting up the algorithm for the PCA, see also discussion below here.
|
||
* Set up the datapoints for the design/feature matrix $\bm{X}$ with $\bm{X}\in {\mathbb{R}}^{n\times p}$, with the predictors/features $p$ referring to the column numbers and the entries $n$ being the row elements.
|
||
!bt
|
||
\[
|
||
\bm{X}=\begin{bmatrix}
|
||
x_{0,0} & x_{0,1} & x_{0,2}& \dots & \dots x_{0,p-1}\\
|
||
x_{1,0} & x_{1,1} & x_{1,2}& \dots & \dots x_{1,p-1}\\
|
||
x_{2,0} & x_{2,1} & x_{2,2}& \dots & \dots x_{2,p-1}\\
|
||
\dots & \dots & \dots & \dots \dots & \dots \\
|
||
x_{n-2,0} & x_{n-2,1} & x_{n-2,2}& \dots & \dots x_{n-2,p-1}\\
|
||
x_{n-1,0} & x_{n-1,1} & x_{n-1,2}& \dots & \dots x_{n-1,p-1}\\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
* Center the data by subtracting the mean value for each column. This leads to a new matrix $\bm{X}\rightarrow \overline{\bm{X}}$.
|
||
* Compute then the covariance/correlation matrix $\mathbb{E}[\overline{\bm{X}}\overline{\bm{X}}^T]$.
|
||
* Find the eigenpairs of $\bm{C}$ with eigenvalues $[\lambda_0,\lambda_1,\dots,\lambda_{p-1}]$ and eigenvectors $[\bm{s}_0,\bm{s}_1,\dots,\bm{s}_{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
|
||
Note that the mean refers to each column of data.
|
||
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 and writing out the design matrix.
|
||
Note that the function _multivariate_ returns also the covariance discussed above and that it is defined by dividing by $n-1$ instead of $n$.
|
||
!bc pycod
|
||
import numpy as np
|
||
import pandas as pd
|
||
from IPython.display import display
|
||
n = 100
|
||
mean = (-1, 2)
|
||
cov = [[4, 2], [2, 2]]
|
||
X = np.random.multivariate_normal(mean, cov, n)
|
||
# Print the X-matrix
|
||
print(X)
|
||
!ec
|
||
|
||
|
||
Try to add to this code your own calculation of the covariance matrix.
|
||
|
||
Now we are going to implement the PCA algorithm. We will break it down into various substeps.
|
||
|
||
=== 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 that 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.
|
||
The following code elements perform these operations using _pandas_ or your own functionality for doing so.
|
||
!bc pycod
|
||
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)
|
||
# test that we get the same as Pandas
|
||
print(X_centered-df)
|
||
!ec
|
||
Alternatively, you could also have used the functions we discussed earlier for scaling the data set.
|
||
That is, we could have used the _StandardScaler_ function in _Scikit-Learn_, a function which
|
||
ensures that for each feature/predictor we study the mean value is
|
||
zero and the variance is one (every column in the design/feature
|
||
matrix).
|
||
|
||
=== 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$.
|
||
We can write our own code or simply use either the functionaly of _numpy_ or that of _pandas_, as follows
|
||
!bc pycod
|
||
print(np.cov(X.T))
|
||
print(df.cov())
|
||
print(np.cov(X_centered.T))
|
||
!ec
|
||
Depending on the number of points $n$, we will get results that are close to the covariance values defined 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 we have these we can perform the
|
||
following tasks:
|
||
|
||
* We compute the percentage of the total variance captured by the first principal component
|
||
* We plot the mean centered data and lines along the first and second principal components
|
||
* Then we project the mean centered data onto the first and second principal components, and plot the projected data.
|
||
* Finally, we 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?
|
||
|
||
Collecting all these steps we can write our own PCA function and
|
||
compare this with the functionality included in _Scikit-Learn_.
|
||
|
||
The code here outlines some of the elements we could include in the analysis. Feel free to extend upon this.
|
||
!bc pycod
|
||
#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)
|
||
#thereafter we do a PCA with Scikit-learn
|
||
from sklearn.decomposition import PCA
|
||
pca = PCA(n_components = 2)
|
||
X2Dsl = pca.fit_transform(X)
|
||
print("Check that we get the same")
|
||
print(X2D-X2Dsl)
|
||
|
||
print(pca.components_.T[:, 0])
|
||
|
||
|
||
!ec
|
||
|
||
|
||
!split
|
||
===== Classical PCA Theorem =====
|
||
|
||
We assume now that we have a design matrix $\bm{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 $[\bm{x}_0,\bm{x}_1,\dots, \bm{x}_{p-1}]$ each with dimension
|
||
$\bm{x}\in {\mathbb{R}}^{n}$.
|
||
|
||
We assume also that we have an orthogonal transformation $\bm{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
|
||
!bt
|
||
\[
|
||
J(\bm{W},\bm{Z}) = \frac{1}{n}\sum_i (\bm{x}_i - \overline{\bm{x}}_i)^2,
|
||
\]
|
||
!et
|
||
with $\overline{\bm{x}}_i = \bm{W}\bm{z}_i$, where $\bm{z}_i$ is a row vector with dimension ${\mathbb{R}}^{n}$ of the matrix
|
||
$\bm{Z}\in{\mathbb{R}}^{p\times n}$. When doing PCA we want to reduce this dimensionality.
|
||
|
||
The PCA theorem states that minimizing the above reconstruction error
|
||
corresponds to setting $\bm{W}=\bm{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
|
||
$\bm{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 proof which follows will be updated by mid January 2020.
|
||
|
||
!split
|
||
===== Proof of the PCA Theorem =====
|
||
|
||
To show the PCA theorem let us start with the assumption that there is one vector $\bm{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 $\bm{w}_0$ and $\bm{z}_0$ as
|
||
!bt
|
||
\[
|
||
J(\bm{w}_0,\bm{z}_0)= \frac{1}{n}\sum_i (\bm{x}_i - z_{i0}\bm{w}_0)^2=\frac{1}{n}\sum_i (\bm{x}_i^T\bm{x}_i - 2z_{i0}\bm{w}_0^T\bm{x}_i+z_{i0}^2\bm{w}_0^T\bm{w}_0),
|
||
\]
|
||
!et
|
||
which we can rewrite due to the orthogonality of $\bm{w}_i$ as
|
||
!bt
|
||
\[
|
||
J(\bm{w}_0,\bm{z}_0)=\frac{1}{n}\sum_i (\bm{x}_i^T\bm{x}_i - 2z_{i0}\bm{w}_0^T\bm{x}_i+z_{i0}^2).
|
||
\]
|
||
!et
|
||
Minimizing $J$ with respect to the unknown parameters $z_{0i}$ we obtain that
|
||
!bt
|
||
\[
|
||
z_{i0}=\bm{w}_0^T\bm{x}_i,
|
||
\]
|
||
!et
|
||
where the vectors on the rhs are known.
|
||
|
||
|
||
!split
|
||
===== PCA Proof continued =====
|
||
|
||
We have now found the unknown parameters $z_{i0}$. These correspond to the projected coordinates and we can write
|
||
!bt
|
||
\[
|
||
J(\bm{w}_0)= \frac{1}{p}\sum_i (\bm{x}_i^T\bm{x}_i - z_{i0}^2)=\mathrm{const}-\frac{1}{n}\sum_i z_{i0}^2.
|
||
\]
|
||
!et
|
||
|
||
We can show that the variance of the projected coordinates defined by $\bm{w}_0^T\bm{x}_i$ are given by
|
||
!bt
|
||
\[
|
||
\mathrm{var}[\bm{w}_0^T\bm{x}_i] = \frac{1}{n}\sum_i z_{i0}^2,
|
||
\]
|
||
!et
|
||
since the expectation value of
|
||
!bt
|
||
\[
|
||
\mathbb{E}[\bm{w}_0^T\bm{x}_i] = \mathbb{E}[z_{i0}]= \bm{w}_0^T\mathbb{E}[\bm{x}_i]=0,
|
||
\]
|
||
!et
|
||
where we have used the fact that our data are centered.
|
||
|
||
Recalling our definition of the covariance as
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x}] = \frac{1}{n}\bm{X}\bm{X}^T=\mathbb{E}[\bm{X}\bm{X}^T],
|
||
\]
|
||
!et
|
||
we have thus that
|
||
!bt
|
||
\[
|
||
\mathrm{var}[\bm{w}_0^T\bm{x}_i] = \frac{1}{n}\sum_i z_{i0}^2=\bm{w}_0^T\bm{C}[\bm{x}]\bm{w}_0.
|
||
\]
|
||
!et
|
||
|
||
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.
|
||
|
||
!split
|
||
===== The final step =====
|
||
|
||
We could trivially maximize the variance of the projection (and
|
||
thereby minimize the error in the reconstruction function) by letting
|
||
the norm-2 of $\bm{w}_0$ go to infinity. However, this norm since we
|
||
want the matrix $\bm{W}$ to be an orthogonal matrix, is constrained by
|
||
$\vert\vert \bm{w}_0 \vert\vert_2^2=1$. Imposing this condition via a
|
||
Lagrange multiplier we can then in turn maximize
|
||
|
||
!bt
|
||
\[
|
||
J(\bm{w}_0)= \bm{w}_0^T\bm{C}[\bm{x}]\bm{w}_0+\lambda_0(1-\bm{w}_0^T\bm{w}_0).
|
||
\]
|
||
!et
|
||
Taking the derivative with respect to $\bm{w}_0$ we obtain
|
||
|
||
!bt
|
||
\[
|
||
\frac{\partial J(\bm{w}_0)}{\partial \bm{w}_0}= 2\bm{C}[\bm{x}]\bm{w}_0-2\lambda_0\bm{w}_0=0,
|
||
\]
|
||
!et
|
||
meaning that
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x}]\bm{w}_0=\lambda_0\bm{w}_0.
|
||
\]
|
||
!et
|
||
_The direction that maximizes the variance (or minimizes the construction error) is an eigenvector of the covariance matrix_! If we left multiply with $\bm{w}_0^T$ we have the variance of the projected data is
|
||
!bt
|
||
\[
|
||
\bm{w}_0^T\bm{C}[\bm{x}]\bm{w}_0=\lambda_0.
|
||
\]
|
||
!et
|
||
|
||
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 $\bm{X}$.
|
||
|
||
The proof
|
||
for the other eigenvectors $\bm{w}_1,\bm{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":"https://mitpress.mit.edu/books/machine-learning-1". 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.
|
||
|
||
Additional part of the proof for the other eigenvectors will be added by mid January 2020.
|
||
|
||
!split
|
||
===== Geometric Interpretation and link with Singular Value Decomposition =====
|
||
|
||
This material will be added by mid January 2020.
|
||
|
||
|
||
!split
|
||
===== Principal Component Analysis =====
|
||
|
||
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
|
||
!bc pycod
|
||
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)
|
||
!ec
|
||
|
||
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.
|
||
!bc pycod
|
||
W2 = V.T[:, :2]
|
||
X2D = X_centered.dot(W2)
|
||
!ec
|
||
|
||
!split
|
||
===== PCA and scikit-learn =====
|
||
|
||
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):
|
||
!bc pycod
|
||
#thereafter we do a PCA with Scikit-learn
|
||
from sklearn.decomposition import PCA
|
||
pca = PCA(n_components = 2)
|
||
X2D = pca.fit_transform(X)
|
||
print(X2D)
|
||
!ec
|
||
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
|
||
!bc pycod
|
||
pca.components_.T[:, 0].
|
||
!ec
|
||
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.
|
||
|
||
!split
|
||
===== 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.
|
||
!bc pycod
|
||
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
|
||
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)))
|
||
|
||
!ec
|
||
|
||
We see that our training data after the PCA decomposition has a performance similar to the non-scaled data.
|
||
|
||
!split
|
||
===== More on the PCA =====
|
||
|
||
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:
|
||
!bc pycod
|
||
pca = PCA()
|
||
pca.fit(X)
|
||
cumsum = np.cumsum(pca.explained_variance_ratio_)
|
||
d = np.argmax(cumsum >= 0.95) + 1
|
||
!ec
|
||
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:
|
||
!bc pycod
|
||
pca = PCA(n_components=0.95)
|
||
X_reduced = pca.fit_transform(X)
|
||
!ec
|
||
|
||
!split
|
||
===== Incremental 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).
|
||
|
||
!split
|
||
===== Randomized PCA =====
|
||
|
||
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$.
|
||
|
||
|
||
|
||
|
||
!split
|
||
===== Kernel PCA =====
|
||
!bblock
|
||
|
||
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
|
||
!bc pycod
|
||
from sklearn.decomposition import KernelPCA
|
||
rbf_pca = KernelPCA(n_components = 2, kernel="rbf", gamma=0.04)
|
||
X_reduced = rbf_pca.fit_transform(X)
|
||
!ec
|
||
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== LLE =====
|
||
|
||
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).
|
||
|
||
|
||
|
||
!split
|
||
===== Other techniques =====
|
||
|
||
|
||
There are many other dimensionality reduction techniques, several of which are available in Scikit-Learn.
|
||
|
||
Here are some of the most popular:
|
||
* _Multidimensional Scaling (MDS)_ reduces dimensionality while trying to preserve the distances between the instances.
|
||
* _Isomap_ creates a graph by connecting each instance to its nearest neighbors, then reduces dimensionality while trying to preserve the geodesic distances between the instances.
|
||
* _t-Distributed Stochastic Neighbor Embedding_ (t-SNE) reduces dimensionality while trying to keep similar instances close and dissimilar instances apart. It is mostly used for visualization, in particular to visualize clusters of instances in high-dimensional space (e.g., to visualize the MNIST images in 2D).
|
||
* Linear Discriminant Analysis (LDA) is actually a classification algorithm, but during training it learns the most discriminative axes between the classes, and these axes can then be used to define a hyperplane onto which to project the data. The benefit is that the projection will keep classes as far apart as possible, so LDA is a good technique to reduce dimensionality before running another classification algorithm such as a Support Vector Machine (SVM) classifier discussed in the SVM lectures.
|
||
|
||
|