In [1]:

import numpy as np
from bokeh.plotting import figure, show, output_notebook

output_notebook()

Loading BokehJS ...

algorithm¶

In [2]:

def PCA(X, n_components):
    # normalize to zero mean
    mu = X.mean(axis=0)
    X = X - mu
    
    # eigenvectors of covariance matrix
    sigma = X.T @ X
    eigvals, eigvecs = np.linalg.eig(sigma)
    
    # principal components
    order = np.argsort(eigvals)[::-1]
    components = eigvecs[:, order[:n_components]]
    
    # projection
    Z = X @ components
    
    # result
    return Z, components

2D data & principal components¶

In [3]:

# generate points
x = np.linspace(0, 13, num=100)
y = x + np.sin(x) - np.cos(x)

# 2D data
X = np.c_[x, y]

In [4]:

# PCA
projection, components = PCA(X, n_components=2)

# principal components
components

Out[4]:

array([[ 0.72299657, -0.69085162],
       [ 0.69085162,  0.72299657]])

In [5]:

# convariance matrix of projected data
(projection.T @ projection).round(3)

Out[5]:

array([[ 2705.081,     0.   ],
       [    0.   ,    47.716]])

In [6]:

# prepare plot data
mean = np.mean(X, axis=0)
extent = projection.min(), projection.max()
angle = np.arctan(components[1] / components[0]) + np.pi * (components[0] < 0)

In [7]:

# plot original data & principal components
plot = figure()

plot.scatter(x, y)
plot.ray(*mean, length=0, angle=angle[0], line_width=2, line_color='red')
plot.ray(*mean, length=0, angle=angle[1], line_width=2, line_color='green')

show(plot)

In [8]:

# plot projected data
plot = figure(x_range=extent, y_range=extent)
plot.scatter(projection[:, 0], projection[:, 1])
show(plot)

binary vectors & dimensionality reduction¶

In [9]:

# generate binary vectors
X = np.random.rand(90, 30)
X[:30, :] = X[:30, :] < ([.4] * 10 + [.1] * 10 + [.1] * 10)
X[30:60, :] = X[30:60, :] < ([.1] * 10 + [.4] * 10 + [.1] * 10)
X[60:, :] = X[60:, :] < ([.1] * 10 + [.1] * 10 + [.4] * 10)

# define 3 classes
Y = ['red'] * 30 + ['green'] * 30 + ['blue'] * 30

In [10]:

# PCA
projection, _ = PCA(X, n_components=2)

In [11]:

# plot projected data: 30D -> 2D
plot = figure()
plot.scatter(projection[:, 0], projection[:, 1], color=Y)
show(plot)

In [ ]: