Diffusion maps, Gaussian kernel

src/soap/soapy/dimred.py

@@ -33,6 +33,10 @@ def dimred_matrix(method, kmat=None, distmat=None, outfile=None, ix=None, symmet
             max_iter=None,
             remove_zero_eig=True)
         positions = kernelpca.fit_transform(kmat)
+    elif method == 'diffmap':
+        positions = diffusion_map(
+            kmat=kmat,
+            n_components=prj_dimension)
     else: raise NotImplementedError(method)
     if outfile: np.savetxt(outfile, positions)
     return positions
@@ -44,6 +48,41 @@ def dimred(kernelpot, method, outfile, symmetrize=False):
     positions = dimred_matrix(method, kmat, distmat, outfile, ix, symmetrize=symmetrize)
     return positions
 
+def diffusion_map(kmat, n_components, alpha=0.5):
+    # Normalize kernel matrix 
+    # (correct for non-unity diagonal components)
+    kmat_diag = kmat.diagonal()
+    kmat_norm = kmat/kmat_diag
+    kmat_norm = kmat_norm.T/kmat_diag
+    np.fill_diagonal(kmat_norm, 1.)
+    #print kmat_norm
+
+    # Rescale kernel
+    D0 = np.sum(kmat_norm, axis=1)
+    D0_diag = np.zeros(kmat_norm.shape)
+    np.fill_diagonal(D0_diag, D0)
+    kmat_norm = kmat_norm/D0**alpha
+    kmat_norm = kmat_norm.T/D0**alpha
+    #print kmat_norm
+
+    # Form Markov matrix
+    D1 = np.sum(kmat_norm, axis=1)
+    kmat_norm = (kmat_norm/D1).T
+    #print kmat_norm
+    #print np.sum(kmat_norm, axis=1)
+
+    # Decompose
+    import scipy.sparse
+    evals, evecs = scipy.sparse.linalg.eigs(kmat_norm, n_components+1, which="LM")
+    order = evals.argsort()[::-1]
+    evals = evals[order]
+    evecs = evecs.T[order].T
+    #print evals
+    #print evecs
+
+    # Project
+    return kmat_norm.dot(evecs[:,1:])
+
 
 
 

src/soap/soapy/kernel.py

@@ -342,6 +342,18 @@ class KernelAdaptorGlobalGeneric(object):
         return dX_dx, dX_dy, dX_dz
 
 
+class KernelFunctionGaussianDiff(object):
+    def __init__(self, options):
+        self.delta = float(options['kernel.gaussian_delta'])
+        self.sigma = float(options['kernel.gaussian_sigma'])
+    def evaluate(self, IX, IY):
+        #res = np.zeros((IX.shape[0],IY.shape[0]))
+        #for i in range(IX.shape[0]):
+        #    for j in range(IY.shape[0]):
+        #        res[i][j] = self.delta**2 * np.exp(-0.5/(self.sigma)**2 * (IX[i]-IY[j]).dot(IX[i]-IY[j]))
+        #return res
+        return self.delta**2 * np.exp(-(1. - IX.dot(IY.T))/self.sigma**2)
+
 class KernelFunctionDotShifted(object):
     """
     For descriptors X,Y where X.dot(Y) e [-1,+1]
@@ -505,6 +517,7 @@ KernelAdaptorFactory = {
 }
 
 KernelFunctionFactory = { 
+'gaussian-diff': KernelFunctionGaussianDiff,
 'dot': KernelFunctionDot, 
 'dot-shifted': KernelFunctionDotShifted,
 'dot-harmonic': KernelFunctionDotHarmonic,