Compressing Data via Dimensionality Reduction

In [1]:

    

import pandas as pd
import matplotlib.pyplot as plt


    

In [2]:

    

df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data', header=None)


    

In [3]:

    

from sklearn.cross_validation import train_test_split
from sklearn.preprocessing import StandardScaler


    

    




/data/bergeric/miniconda3/envs/s2rnai/lib/python3.6/site-packages/sklearn/cross_validation.py:41: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.
  "This module will be removed in 0.20.", DeprecationWarning)

In [4]:

    

X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values


    

In [5]:

    

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)


    

In [6]:

    

#first standardize data
sc = StandardScaler()


    

In [7]:

    

X_train_std = sc.fit_transform(X_train)


    

In [8]:

    

X_test_std = sc.transform(X_test)


    

In [9]:

    

import numpy as np


    

In [10]:

    

#Obtain covariance matrix 
cov_mat = np.cov(X_train_std.T)
#find eigenvalues/vectors
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
print('Eigenvalues \n%s' % eigen_vals)


    

    




Eigenvalues 
[ 4.8923083   2.46635032  1.42809973  1.01233462  0.84906459  0.60181514
  0.52251546  0.08414846  0.33051429  0.29595018  0.16831254  0.21432212
  0.2399553 ]

In [11]:

    

tot = sum(eigen_vals)


    

In [12]:

    

#Plot explained variance
var_exp = [(i/tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)

plt.bar(range(1,14), var_exp, alpha=0.5, align='center', label='individual explained variance')
plt.step(range(1,14), cum_var_exp, where='mid', label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.show()


    

In [13]:

    

#Sort eigenpairs by descending order of the eigenvalues
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]
eigen_pairs.sort(reverse=True)


    

In [14]:

    

#use two eigenvectors that corresepond to two largest e.values to create projection matrix W
w = np.hstack((eigen_pairs[0][1][:, np.newaxis], eigen_pairs[1][1][:, np.newaxis]))
print('Matrix W:\n', w)


    

    




Matrix W:
 [[ 0.14669811  0.50417079]
 [-0.24224554  0.24216889]
 [-0.02993442  0.28698484]
 [-0.25519002 -0.06468718]
 [ 0.12079772  0.22995385]
 [ 0.38934455  0.09363991]
 [ 0.42326486  0.01088622]
 [-0.30634956  0.01870216]
 [ 0.30572219  0.03040352]
 [-0.09869191  0.54527081]
 [ 0.30032535 -0.27924322]
 [ 0.36821154 -0.174365  ]
 [ 0.29259713  0.36315461]]

In [15]:

    

#transform a sample x onto PCA subspace (x' = xW)
X_train_std[0].dot(w)


    

    
Out[15]:





array([ 2.59891628,  0.00484089])

In [16]:

    

#transform entire training dataset (X' = XW)
X_train_pca = X_train_std.dot(w)


    

In [17]:

    

colors = ['r','b','g']
markers = ['s','x','o']
for l,c,m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train==l, 0], X_train_pca[y_train==l, 1], c=c, label=l, marker=m)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.show()


    

In [18]:

    

from matplotlib.colors import ListedColormap


    

In [19]:

    

def plot_decision_regions(X,y, classifier, resolution=0.02):
    #setup marker generator and color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])
    
    #plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution), np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())
    
    #plot class samples
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1], alpha=0.8, c=cmap(idx), marker=markers[idx], label=cl)


    

In [20]:

    

from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA


    

In [21]:

    

pca = PCA(n_components=2)


    

In [22]:

    

lr = LogisticRegression()


    

In [23]:

    

X_train_pca = pca.fit_transform(X_train_std)


    

In [24]:

    

X_test_pca = pca.transform(X_test_std)


    

In [25]:

    

lr.fit(X_train_pca, y_train)


    

    
Out[25]:





LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,
          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,
          verbose=0, warm_start=False)

In [26]:

    

plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc='lower left')
plt.show()


    

In [27]:

    

plot_decision_regions(X_test_pca, y_test, classifier=lr)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc='lower left')
plt.show()


    

In [28]:

    

#initialize PCA class with n_components parameter set to None in order to access explained_variance_ratio: 

pca = PCA(n_components=None)
X_train_pca = pca.fit_transform(X_train_std)
pca.explained_variance_ratio_


    

    
Out[28]:





array([ 0.37329648,  0.18818926,  0.10896791,  0.07724389,  0.06478595,
        0.04592014,  0.03986936,  0.02521914,  0.02258181,  0.01830924,
        0.01635336,  0.01284271,  0.00642076])

In [29]:

    

#compute mean vector for each class
np.set_printoptions(precision=4)
mean_vecs = []
for label in range(1,4):
    mean_vecs.append(np.mean(X_train_std[y_train==label], axis=0))
    print('MV %s: %s\n' %(label, mean_vecs[label-1]))


    

    




MV 1: [ 0.9259 -0.3091  0.2592 -0.7989  0.3039  0.9608  1.0515 -0.6306  0.5354
  0.2209  0.4855  0.798   1.2017]

MV 2: [-0.8727 -0.3854 -0.4437  0.2481 -0.2409 -0.1059  0.0187 -0.0164  0.1095
 -0.8796  0.4392  0.2776 -0.7016]

MV 3: [ 0.1637  0.8929  0.3249  0.5658 -0.01   -0.9499 -1.228   0.7436 -0.7652
  0.979  -1.1698 -1.3007 -0.3912]

In [30]:

    

d = 13  # number of features
S_W = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vecs):
    class_scatter = np.zeros((d, d))  # scatter matrix for each class
    for row in X_train_std[y_train == label]:
        row, mv = row.reshape(d, 1), mv.reshape(d, 1)  # make column vectors
        class_scatter += (row - mv).dot((row - mv).T)
    S_W += class_scatter                          # sum class scatter matrices

print('Within-class scatter matrix: %sx%s' % (S_W.shape[0], S_W.shape[1]))


    

    




Within-class scatter matrix: 13x13

In [31]:

    

print('Class label distribution: %s' % np.bincount(y_train)[1:])


    

    




Class label distribution: [40 49 35]

In [32]:

    

#scale scatter matrices first
d = 13  # number of features
S_W = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vecs):
    class_scatter = np.cov(X_train_std[y_train == label].T)
    S_W += class_scatter
print('Scaled within-class scatter matrix: %sx%s' % (S_W.shape[0],
                                                     S_W.shape[1]))


    

    




Scaled within-class scatter matrix: 13x13

In [33]:

    

#compute between-class scatter matrix:
mean_overall = np.mean(X_train_std, axis=0)
d = 13  # number of features
S_B = np.zeros((d, d))
for i, mean_vec in enumerate(mean_vecs):
    n = X_train[y_train == i + 1, :].shape[0]
    mean_vec = mean_vec.reshape(d, 1)  # make column vector
    mean_overall = mean_overall.reshape(d, 1)  # make column vector
    S_B += n * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)

print('Between-class scatter matrix: %sx%s' % (S_B.shape[0], S_B.shape[1]))


    

    




Between-class scatter matrix: 13x13

In [34]:

    

#find e.pairs of (Sw)-1(Sb)
eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))


    

In [35]:

    

eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:,i]) for i in range(len(eigen_vals))]


    

In [36]:

    

eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)
print('Eigenvalues in decreasing order:\n')
for eigen_val in eigen_pairs:
    print(eigen_val[0])


    

    




Eigenvalues in decreasing order:

452.721581245
156.43636122
4.32282196963e-14
3.34763921076e-14
2.84217094304e-14
2.11943360486e-14
2.11943360486e-14
2.05859329275e-14
1.20711099662e-14
1.20711099662e-14
1.08538563383e-14
1.08538563383e-14
2.96988251897e-15

In [37]:

    

#plot 'discriminability' of linear discriminants
tot = sum(eigen_vals.real)
discr = [(i / tot) for i in sorted(eigen_vals.real, reverse=True)]
cum_discr = np.cumsum(discr)
plt.bar(range(1, 14), discr, alpha=0.5, align='center', label='individual "discriminability"')
plt.step(range(1, 14), cum_discr, where='mid', label='cumulative "discriminability"')
plt.ylabel('"discriminability" ratio')
plt.xlabel('Linear Discriminants')
plt.ylim([-0.1, 1.1])
plt.legend(loc='best')
plt.show()


    

In [38]:

    

#stack two most discriminative eigenvector columns to create proj matrix W: 
w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real, eigen_pairs[1][1][:, np.newaxis].real))
print('Matrix W:\n', w)


    

    




Matrix W:
 [[-0.0662 -0.3797]
 [ 0.0386 -0.2206]
 [-0.0217 -0.3816]
 [ 0.184   0.3018]
 [-0.0034  0.0141]
 [ 0.2326  0.0234]
 [-0.7747  0.1869]
 [-0.0811  0.0696]
 [ 0.0875  0.1796]
 [ 0.185  -0.284 ]
 [-0.066   0.2349]
 [-0.3805  0.073 ]
 [-0.3285 -0.5971]]

In [39]:

    

#project samples onto new feature space
X_train_lda = X_train_std.dot(w)


    

In [40]:

    

colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_lda[y_train==l, 0]*(-1), X_train_lda[y_train==l, 1]*(-1), c=c, label=l, marker=m)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower right')
plt.show()


    

In [41]:

    

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA


    

In [42]:

    

lda = LDA(n_components=2)


    

In [43]:

    

X_train_lda = lda.fit_transform(X_train_std, y_train)


    

In [44]:

    

#see how logistic regression handles lower-dimensional training set
lr = LogisticRegression()
lr.fit(X_train_lda, y_train)


    

    
Out[44]:





LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,
          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,
          verbose=0, warm_start=False)

In [45]:

    

plot_decision_regions(X_train_lda, y_train, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.show()


    

In [46]:

    

X_test_lda = lda.transform(X_test_std)


    

In [47]:

    

plot_decision_regions(X_test_lda,y_test, classifier=lr)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend(loc='lower left')
plt.show()


    

In [48]:

    

from scipy.spatial.distance import pdist, squareform
from scipy import exp
from scipy.linalg import eigh


    

In [49]:

    

def rbf_kernel_pca(X, gamma, n_components):
    """
    RBF kernel PCA implementation.
    
    Parameters
    -----------
    X: {NumPy ndarray}, shape = [n_samples, n_features]
    
    gamma: float, Tuning parameter of RBF kernel
    
    n_components: int, Number of principal components to return
    
    Returns
    -----------
    X_pc: {NumPy ndarray}, shape = [n_samples, k_features], Projected dataset
    """
    
    #Calculate pairwise squared Euclidean distances
    sq_dists = pdist(X, 'sqeuclidean')
    
    #convert pairwise distances into square matrix
    mat_sq_dists = squareform(sq_dists)
    
    #Compute symmetric kernel matrix
    K = exp(-gamma * mat_sq_dists)
    
    #Center kernel matrix
    N = K.shape[0]
    one_n = np.ones((N, N)) / N
    K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)
    
    #Obtain eigenpairs from centered kernel matrix
    eigvals, eigvecs = eigh(K)
    
    # Collect top k eigenvectors(projected samples)
    X_pc = np.column_stack((eigvecs[:, -i] for i in range(1, n_components + 1)))
    
    return X_pc


    

In [50]:

    

# EXAMPLE 1 : separate half-moons

from sklearn.datasets import make_moons
X, y = make_moons(n_samples=100, random_state=123)


    

In [51]:

    

plt.scatter(X[y==0,0], X[y==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X[y==1,0], X[y==1, 1], color='blue', marker='o', alpha=0.5)
plt.show()


    

In [52]:

    

#Clearly not linearly separable
#What does standard PCA look like? 

scikit_pca = PCA(n_components=2)


    

In [53]:

    

X_spca = scikit_pca.fit_transform(X)


    

In [54]:

    

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7,3))
ax[0].scatter(X_spca[y==0, 0], X_spca[y==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(X_spca[y==1, 0], X_spca[y==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_spca[y==0, 0], np.zeros((50,1))+0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(X_spca[y==1, 0], np.zeros((50, 1))-0.02, color='blue', marker='o', alpha=0.5)
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.show()


    

In [55]:

    

#Now try our kernel PCA

from matplotlib.ticker import FormatStrFormatter

X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)


    

In [56]:

    

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7,3))
ax[0].scatter(X_kpca[y==0, 0], X_kpca[y==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(X_kpca[y==1, 0], X_kpca[y==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_kpca[y==0, 0], np.zeros((50,1)) + 0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(X_kpca[y==1, 0], np.zeros((50,1)) - 0.02, color='blue', marker='o', alpha=0.5)
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
ax[0].xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))
ax[1].xaxis.set_major_formatter(FormatStrFormatter('%0.1f'))
plt.show()


    

In [57]:

    

## EXAMPLE 2: separate concentric circles

from sklearn.datasets import make_circles

X,y = make_circles(n_samples=1000, random_state=123, noise=0.1, factor=0.2)


    

In [58]:

    

plt.scatter(X[y==0, 0], X[y==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X[y==1, 0], X[y==1, 1], color='blue', marker='o', alpha=0.5)
plt.show()


    

In [59]:

    

#with std PCA
scikit_pca = PCA(n_components=2)
X_spca = scikit_pca.fit_transform(X)


    

In [60]:

    

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7,3))
ax[0].scatter(X_spca[y==0, 0], X_spca[y==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(X_spca[y==1, 0], X_spca[y==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_spca[y==0, 0], np.zeros((500,1))+0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(X_spca[y==1, 0], np.zeros((500, 1))-0.02, color='blue', marker='o', alpha=0.5)
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.show()


    

In [61]:

    

#with RBF kernel
X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)


    

In [62]:

    

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(7,3))
ax[0].scatter(X_kpca[y==0, 0], X_kpca[y==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(X_kpca[y==1, 0], X_kpca[y==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(X_kpca[y==0, 0], np.zeros((500,1)) + 0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(X_kpca[y==1, 0], np.zeros((500,1)) - 0.02, color='blue', marker='o', alpha=0.5)
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.show()


    

In [63]:

    

#Projecting data points not part of training dataset
#modify function to return eigenvalues of kernel matrix

def rbf_kernel_pca(X, gamma, n_components):
    """
    RBF kernel PCA implementation.
    
    Parameters
    -----------
    X: {NumPy ndarray}, shape = [n_samples, n_features]
    
    gamma: float, Tuning parameter of RBF kernel
    
    n_components: int, Number of principal components to return
    
    Returns
    -----------
    X_pc: {NumPy ndarray}, shape = [n_samples, k_features], Projected dataset
    
    lambdas: list, Eigenvalues
    """
    
    #Calculate pairwise squared Euclidean distances
    sq_dists = pdist(X, 'sqeuclidean')
    
    #convert pairwise distances into square matrix
    mat_sq_dists = squareform(sq_dists)
    
    #Compute symmetric kernel matrix
    K = exp(-gamma * mat_sq_dists)
    
    #Center kernel matrix
    N = K.shape[0]
    one_n = np.ones((N, N)) / N
    K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)
    
    #Obtain eigenpairs from centered kernel matrix
    eigvals, eigvecs = eigh(K)
    
    # Collect top k eigenvectors(projected samples)
    alphas = np.column_stack((eigvecs[:, -i] for i in range(1, n_components + 1)))
    
    #Collect corresponding eigenvalues
    lambdas = [eigvals[-i] for i in range(1, n_components + 1)]
    
    return alphas, lambdas


    

In [64]:

    

X, y = make_moons(n_samples=100, random_state=123)


    

In [65]:

    

alphas, lambdas= rbf_kernel_pca(X, gamma=15, n_components=1)


    

In [66]:

    

#test with projecting 26th datapoint
x_new = X[25]
x_new


    

    
Out[66]:





array([ 1.8713,  0.0093])

In [67]:

    

x_proj = alphas[25] #original projection
x_proj


    

    
Out[67]:





array([ 0.0788])

In [69]:

    

def project_x(x_new, X, gamma, alphas, lambdas):
    pair_dist = np.array([np.sum((x_new - row)**2) for row in X])
    k = np.exp(-gamma * pair_dist)
    return k.dot(alphas / lambdas)


    

In [70]:

    

#reproduce original projection
x_reproj = project_x(x_new, X, gamma=15, alphas=alphas, lambdas=lambdas)
x_reproj


    

    
Out[70]:





array([ 0.0788])

In [71]:

    

#visualize projection on first principal component

plt.scatter(alphas[y==0, 0], np.zeros((50)), color='red', marker='^', alpha=0.5)
plt.scatter(alphas[y==1, 0], np.zeros((50)), color='blue', marker='o', alpha=0.5)
plt.scatter(x_proj, 0, color='black', label='original projection of point X[25]', marker='^', s=100)
plt.scatter(x_reproj, 0, color='green', label='remapped point X[25]', marker='x', s=500)
plt.legend(scatterpoints=1)
plt.show()


    

In [72]:

    

from sklearn.decomposition import KernelPCA


    

In [73]:

    

X, y = make_moons(n_samples=100, random_state=123)
scikit_kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)


    

In [74]:

    

X_skernpca = scikit_kpca.fit_transform(X)


    

In [75]:

    

plt.scatter(X_skernpca[y==0, 0], X_skernpca[y==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(X_skernpca[y==1, 0], X_skernpca[y==1, 1], color='blue', marker='o', alpha=0.5)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.show()


    

In [ ]: