classical mds python

1from __future__ import division
2 
3import numpy as np
4 
5def cmdscale(D):
6    """                                                                                       
7    Classical multidimensional scaling (MDS)                                                  
8                                                                                               
9    Parameters                                                                                
10    ----------                                                                                
11    D : (n, n) array                                                                          
12        Symmetric distance matrix.                                                            
13                                                                                               
14    Returns                                                                                   
15    -------                                                                                   
16    Y : (n, p) array                                                                          
17        Configuration matrix. Each column represents a dimension. Only the                    
18        p dimensions corresponding to positive eigenvalues of B are returned.                 
19        Note that each dimension is only determined up to an overall sign,                    
20        corresponding to a reflection.                                                        
21                                                                                               
22    e : (n,) array                                                                            
23        Eigenvalues of B.                                                                     
24                                                                                               
25    """
26    # Number of points                                                                        
27    n = len(D)
28 
29    # Centering matrix                                                                        
30    H = np.eye(n) - np.ones((n, n))/n
31 
32    # YY^T                                                                                    
33    B = -H.dot(D**2).dot(H)/2
34 
35    # Diagonalize                                                                             
36    evals, evecs = np.linalg.eigh(B)
37 
38    # Sort by eigenvalue in descending order                                                  
39    idx   = np.argsort(evals)[::-1]
40    evals = evals[idx]
41    evecs = evecs[:,idx]
42 
43    # Compute the coordinates using positive-eigenvalued components only                      
44    w, = np.where(evals > 0)
45    L  = np.diag(np.sqrt(evals[w]))
46    V  = evecs[:,w]
47    Y  = V.dot(L)
48 
49    return Y, evals
50