# -*- coding: utf-8 -*- from __future__ import division import numpy as np def cmdscale(D): """ Classical multidimensional scaling (MDS) Code by Francis Song (song.francis@gmail.com) http://www.nervouscomputer.com/hfs/cmdscale-in-python/ Parameters ---------- D : (n, n) array Symmetric distance matrix. Returns ------- Y : (n, p) array Configuration matrix. Each column represents a dimension. Only the p dimensions corresponding to positive eigenvalues of B are returned. Note that each dimension is only determined up to an overall sign, corresponding to a reflection. e : (n,) array Eigenvalues of B. """ # Number of points n = len(D) # Centering matrix H = np.eye(n) - np.ones((n, n))/n # YY^T B = -H.dot(D**2).dot(H)/2 # Diagonalize evals, evecs = np.linalg.eigh(B) # Sort by eigenvalue in descending order idx = np.argsort(evals)[::-1] evals = evals[idx] evecs = evecs[:,idx] # Compute the coordinates using positive-eigenvalued components only w, = np.where(evals > 0) L = np.diag(np.sqrt(evals[w])) V = evecs[:,w] Y = V.dot(L) return Y, evals