User Tools

Site Tools


mds_class

Visual Analytics Course

MDS (Multi-dimensional scaling)

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

dissimilarityMatrix = matrix(...) # some matrix you built or read from disk
 
# k is the dimension of the space you project the data set onto
mds = MDS(dissimilarityMatrix, k) 
mds.embed() # embed data onto a kD space

MDS.py
from numpy import *
import random as rndm
 
class MDS():
    '''
    A simple class implementing MDS projection using the
    SMACOF algorithm from Borg and Groenen
    '''
 
    def __init__(self, dissimilarityMatrix, k):
        '''
        Dissimilarity matrix has to be given - make sure it is a numpy matrix.
        Of course, no need to say the matrix must be square and symmetric.
 
        N is the number of elements we are given
        k is the dimension of the Euclidean space we wish to embed the data in
        '''
        self.dissimilarityMatrix = dissimilarityMatrix
        self.N = self.dissimilarityMatrix.shape[0]
        self.k = k
        self.embeddingMatrix = matrix([[0.0] * self.k] * self.N)
        self.distanceMatrix = matrix([[0.0] * self.N] * self.N)
        self.stress = 0.0
        self.Vmatrix = matrix([[0.0] * self.N] * self.N)
 
    def __distance__(self, vector1, vector2):
        '''
        computes Euclidean distance between any two points
        assumes dimension of given vectors is self.k
 
        assumes vectors are described as arrays
        '''
        d = 0.0
        for i in range(self.k):
            d += (vector1[i] - vector2[i])**2
        return sqrt(d)
 
    def __computeDistanceMatrix__(self):
        '''
        Computes all d_ij distances used when computing stress, for instance.
        '''
        for i in range(self.N):
            self.distanceMatrix[i, i] = 0.0
            for j in range(i+1, self.N):
                xi = self.embeddingMatrix.A[i]
                xj = self.embeddingMatrix.A[j]
                self.distanceMatrix[i, j] = self.__distance__(xi, xj)
                self.distanceMatrix[j, i] = self.distanceMatrix[i, j]
        return True
 
    def __rawStress__(self):
        stress = 0.0
        for i in range(self.N):
            for j in range(i+1, self.N):
                 stress += (self.dissimilarityMatrix[i, j] - self.distanceMatrix[i, j])**2
        return sqrt(stress)
 
    def __eVector__(self, i):
        '''
        Integers i should belong to the set {1, 2, ..., N}
        to stick with notation sused inthe book
 
        No precaution is taken to make sire i is in the proper range
        '''
        arrayVector = [0] * self.N
        arrayVector[i-1] = 1
        return transpose(matrix(arrayVector))
 
    def __Aij__(self, i, j):
        '''
        Integers i should belong to the set {1, 2, ..., N}
        to stick with notation sused inthe book
 
        No precaution is taken to make sire i is in the proper range
        '''
        return dot((self.__eVector__(i)-self.__eVector__(j)), transpose((self.__eVector__(i)-self.__eVector__(j))))
 
 
    def __computeVmatrix__(self):
        '''
        Integers i should belong to the set {1, 2, ..., N}
        to stick with notation sused inthe book
 
        No precaution is taken to make sire i is in the proper range
        '''
        for i in range(self.N):
            for j in range(i+1, self.N):
                 self.Vmatrix += self.__Aij__(i, j)
        return True
 
    def __Guttman__(self, Z):
        '''
        computes the Guttman transform of a matrix
        encoding an embedding Z also given as a matrix
        '''
        I = identity(self.N)
        centeringMatrix = 1.0 / self.N * matrix([[1.0] * self.N] * self.N)
        MoorePenrose = 1.0 / self.N * (I - centeringMatrix)
 
        return dot(dot(MoorePenrose, self.__Bmatrix__(Z)), Z)
 
    def __Bmatrix__(self, Z):
        '''
        computes the B matrix from a given embedding Z given as a amtrix
        '''
        B = matrix([[0.0] * self.N] * self.N)
        for i in range(self.N):
            for j in range(i+1, self.N):
                zi = Z.A[i]
                zj = Z.A[j]
                dij =  self.__distance__(zi, zj)
                if dij != 0.0:
                    deltaij = self.dissimilarityMatrix[i, j]
                    B[i, j] = - deltaij / dij
                    B[j, i] = B[i, j]
        for i in range(self.N):
            for j in range(0,i):
                B[i, i] -= B[i, j]
            for j in range(i+1, self.N): 
                B[i, i] -= B[i, j]
        return B
 
    def __randomPositions__(self):
        '''
        Randomly positions data elements -- used to bootstrap the algorithm
        '''
        diameter = self.dissimilarityMatrix.max()
        side = diameter / sqrt(self.k)
        for i in range(self.N):
            for j in range(self.k):
                self.embeddingMatrix[i, j] = rndm.random() * side
 
    def embed(self, epsilon, maxIter):
        '''
        This is the SMACOF algorithm relying on all previous methods.
        '''
        self.__randomPositions__()
        self.__computeDistanceMatrix__()
        Z = self.embeddingMatrix
        prevStress = self.__rawStress__()
        self.embeddingMatrix = self.__Guttman__(Z)
        self.__computeDistanceMatrix__()
        stress = self.__rawStress__()
        nbIter = 0
        while abs(stress - prevStress) >= epsilon or nbIter == maxIter:
            print 'iteration ' + str(nbIter) + ' : ' + str(stress)
            Z = self.embeddingMatrix
            prevStress = stress
            self.embeddingMatrix = self.__Guttman__(Z)
            self.__computeDistanceMatrix__()
            stress = self.__rawStress__()
            nbIter += 1
        return True
/net/html/perso/melancon/Visual_Analytics_Course/data/pages/mds_class.txt · Last modified: 2012/10/25 14:52 by melancon