pca.py

# -----------------------------------------------------------------------------
# Copyright 2019 (C) Jeremy Fix
# Released under a BSD two-clauses license
#
# Reference: Pearson, K. (1901). "On Lines and Planes of Closest Fit to Systems
#            of Points in Space". Philosophical Magazine. 2 (11): 559–572.
#            doi:10.1080/14786440109462720
# -----------------------------------------------------------------------------
import numpy as np

import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import matplotlib.cm as cmx
import matplotlib.colors as colors

def get_cmap(N):
    '''Returns a function that maps each index in 0, 1, ... N-1 to a distinct
        RGB color.'''
    color_norm  = colors.Normalize(vmin=0, vmax=N)
    scalar_map = cmx.ScalarMappable(norm=color_norm, cmap='hsv')
    def map_index_to_rgb_color(index):
        return scalar_map.to_rgba(index)
    return map_index_to_rgb_color


class PCA:
    ''' PCA class '''

    def __init__(self, n_components):
        self.n_components = n_components

    def fit(self, X):
        '''
        Performs the PCA of X and stores the principal components.
        The datapoints are supposed to be stored in the row vectors of X.
        It keeps only the n_components projection vectors, associated
        with the n_components largest eigenvalues of the covariance matrix X.T X
        '''

        # Center the datapoints
        self.centroid = np.mean(X, axis=0)

        # Computes the covariance matrix
        sigma = np.dot((X - self.centroid).T, X - self.centroid)

        # Compute the eigenvalue/eigenvector decomposition of sigma
        eigvals, eigvecs = np.linalg.eigh(sigma)

        # Note :The eigenvalues returned by eigh are ordered in ascending order.

        # Stores the n_components eigvenvectors/eigenvalues associated
        # with the largest eigen values
        #self.eigvals = eigvals[-self.n_components:]
        #self.eigvecs = eigvecs[:, -self.n_components:]

        # Stores all the eigenvectors/eigenvalues
        # Usefull for latter computing some statistics of the variance we keep
        self.eigvals = eigvals
        self.eigvecs = eigvecs

    def transform(self, Z):
        '''
        Uses a fitted PCA to transform the row vectors of Z
        Remember the eigen vectors are ordred by ascending order
        Denoting Z_trans = transform(Z),
        The first  component is Z_trans[:, -1]
        The second component is Z_trans[:, -2]
        ...
        '''
        return np.dot(Z - self.centroid, self.eigvecs[:, -self.n_components:])


if __name__ == '__main__':


    # Example 1 : Random 2D data
    # -------------------------------------------------------------------------
    print("Random data example")
    print('-' * 70)

    N = 500

    ## Generate some random normal data rotated and translated
    X = np.column_stack([np.random.normal(0.0, 1.0, (N, )),
                         np.random.normal(0.0, 0.2, (N, ))])
    X = np.dot(X, np.array([[np.cos(1.0), np.sin(1.0)],[-np.sin(1.0), np.cos(1.0)]]))
    X = np.array([0.5, 1.0]) + X

    ## Extract the 2 first principal component vectors
    ## With a 2D dataset, these are the only interesting components
    n_components = 2
    pca = PCA(n_components=n_components)
    pca.fit(X)

    ## Project the original data
    X_trans = pca.transform(X)

    print("{:.2f}% of the variance along the first axis (green)".format(
          100 * pca.eigvals[-1] / pca.eigvals.sum()))

    print("{:.2f}% of the variance along the second axis (red)".format(
          100 * pca.eigvals[-2] / pca.eigvals.sum()))

    ## Plot
    plt.figure()
    plt.scatter(X[:,0], X[:,1], alpha = .25)
    ## Plot the first projecton axis in green
    plt.plot([pca.centroid[0]-2*pca.eigvecs[0, -1], pca.centroid[0]+2*pca.eigvecs[0,-1]],
             [pca.centroid[1]-2*pca.eigvecs[1, -1], pca.centroid[1]+2*pca.eigvecs[1,-1]],
             'g', linewidth=3)
    ## Plot the second projection axis in red
    plt.plot([pca.centroid[0]-2*pca.eigvecs[0, 0], pca.centroid[0]+2*pca.eigvecs[0,0]],
             [pca.centroid[1]-2*pca.eigvecs[1, 0], pca.centroid[1]+2*pca.eigvecs[1,0]],
             'r', linewidth=3)
    plt.xlim([-3, 3])
    plt.ylim([-3, 3])
    plt.gca().set_aspect('equal')
    plt.savefig("random_points.png")
    print()

    # Example 2 : 28 x 28 images of digits
    # -------------------------------------------------------------------------
    print("Digits example")
    print('-' * 70)

    samples = np.load('dig_app_text.cb.npy')
    X, y = samples['input'], samples['label']

    ## Extract the 10 first principal component vectors
    n_components = 10
    pca = PCA(n_components = n_components)
    pca.fit(X)

    ## Project the original data
    X_trans = pca.transform(X)

    print("{:.2f}% of the variance is kept with {} components".format(
          100 * pca.eigvals[-n_components:].sum()/pca.eigvals.sum(),
          n_components))

    ## Plot
    fig = plt.figure(figsize=(10,4), tight_layout=True)
    gs  = GridSpec(3, 10)

    ### Coordinates in the projected space
    ### the color shows how the digits get separated by the principal vectors
    ax = fig.add_subplot(gs[0:2, :-2])
    cmap = get_cmap(10)
    ax.scatter(X_trans[:, -1], X_trans[:,-2] , c=[cmap(l) for l in y], alpha = .25)
    ax = fig.add_subplot(gs[0:2, -2:])
    ax.set_axis_off()
    for i in range(10):
        col = cmap(i)
        rect = plt.Rectangle((0.1, i/10.), 0.4, 1.0/14., facecolor=col)
        ax.add_artist(rect)
        ax.text(0.6, i/10.+1./28., str(i), fontsize=12, color="k", **{'ha':'center', 'va':'center'})


    ### Projection vectors
    for i in range(n_components):
        ax = fig.add_subplot(gs[2, i])
        ax.imshow(pca.eigvecs[:, -(i+1)].reshape((28, 28)), cmap='gray_r')
        ax.set_xticks([])
        ax.set_yticks([])

    plt.savefig("digits.png")
    plt.show()