User Tools

Site Tools


Visual Analytics Course

PCA (Principal Components Analysis)

Here is a piece of python code that implements the PCA algorithm using a data matrix . You would typically use it as:

data = matrix(...) # some matrix you built or read from disk
pca = PCA(data)
pca.computeDataProjection() # projects data onto a 2D space
from numpy import *
class PCA():
    a number of functions computing PCA
    and related objects
    vectors a line vectors
    def __init__(self, dataMatrix):
        The method expects a matrix containing all tabular data
        self.dataMatrix = array(dataMatrix)
        # uno is a 1d vector filled with 1's = array([1] * self.dataMatrix.shape[0])
            weights is a diagonal matrix with data weights for each data element
            initialized with equal weights
        self.weights = identity(self.dataMatrix.shape[0], float) * 1.0 / self.dataMatrix.shape[0]
        self.barycenter = None
        self.covarianceMatrix = None
        self.normalizedDataMatrix = None    
    def __computeBarycenter__(self):
        self.barycenter = dot(transpose(self.dataMatrix), dot(self.weights,
        return True
    def __dotVector__(self, v1, v2):
        v1 and v2 are 1d vectors, the array class dot function would compute
        their scalar product
        array class does not know how to handle the case
        where one wants to multiply a column vector with a row vector
        to output a matrix
        this is what this function does
        m = []
        for i in range(len(v1)):
            for j in range(len(v2)):
                m.append(v1[i] * v2[j])
        m = array(m)
        return m.reshape(len(v1), len(v2))
    def __computeCovarianceMatrix__(self):
        Data normalization
        if self.barycenter == None:
        Y = self.dataMatrix - self.__dotVector__(, self.barycenter)
        V = dot(dot(transpose(Y), self.weights), Y)
        self.centeredDataMatrix = Y
        d = zeros([self.dataMatrix.shape[1], self.dataMatrix.shape[1]])
        for i in range(self.dataMatrix.shape[1]):
            d[i, i] = 1.0 / std(self.dataMatrix[:, i])
        self.normalizedDataMatrix = dot(Y, d)
        self.covarianceMatrix = dot(dot(d, V), d)
        return True
    def computePrincipalComponents(self):
        Computes *all* eigenValues and eigenVectors of
        the D^T D matrix
        if self.covarianceMatrix == None:
        self.eigValues, self.eigVectors = linalg.eig(self.covarianceMatrix)
        return True
    def computeDataProjection(self):
        Computes a 2D projection of the data usgin the eigenvectors
        associated with the two largest eigenValues
        Assumes method computeCovarianceMatrix has been previously called
        m = zeros([self.dataMatrix.shape[1], self.dataMatrix.shape[1]])
        for i in range(self.dataMatrix.shape[1]):
            m[i, i] = 1.0 / var(self.dataMatrix[:, i])
        return [dot(dot(self.normalizedDataMatrix, m), self.eigVectors[:, 0]), dot(dot(self.normalizedDataMatrix, m), self.eigVectors[:, 1])]
    def cumulativeVariance(self):
        Computes all ratios
        (lambda_1 + ... + lambda_i)/(lambda_1 + ... + lambda_N)
        Looking at the first two values is relevant when projecting in 2D
        s = 0.0
        for x in self.eigValues:
            s += abs(x)
        v = []
        for x in self.eigValues:
            v.append(abs(x) / s)
        return v
    def __project__(self, point, v1 = None, v2 = None):
        Projects original data points onto the 2D space
        defined by the two given vectors
        When no vector are given, the method uses the
        eigenvectors associated with the two largest eigen values
        if v1 == None:
            v1 = self.eigVectors[:, 0]
        if v2 == None:
            v2 = self.eigVectors[:, 1]
        return [dot(p, v1), dot(p, v2)]
/net/html/perso/melancon/Visual_Analytics_Course/data/pages/pca_class.txt · Last modified: 2012/10/25 14:57 by melancon