<jupyter_start><jupyter_text>Diffusion maps 1 Author: <NAME>, 2019  AbstractThis is my first sight to diffusion maps [1],  I study [Ann Lee's Matlab code](https://reurl.cc/E3Ykv) and simply convert it into Python. In order to understand the meaning of the code, I also read papers that describe this method [1, 2]. In the following article, I will show how Ann Lee constructs a diffusion map step by step, then demonstrate my Python code. Also, there's a paper provided a better way to build one map [3], it leads me to the [next project](), in which I will implement it's algorithm.  FrameworkGiven a data set $X=\{x_1,...,x_n\} \in \mathbb{R}^p$, Ann Lee construct the diffusion maps by using following steps:  1. Define a kernel matrixLet $K$ be a kernel matrix with entries $ K_{ij}=k(x_{i},x_{j}) $, where $k$ is the Guassian kernel$$ k(x,y) = \exp(-\frac{\lVert x-y \rVert^2}{4\sigma}), $$and $\sigma$ is a manually picked constant.  2. Construct a matrix that related to diffusion matrixDiffusion matrix $P$ is a matrix that describe the connectivity between two data points, which can be consider as the probability of jumping from one point to another in one step of random walk [2].  Consider a matrix $Q$, with it's entries $$ Q_{ij}=\frac{K_{ij}}{\sqrt{\sum_{j=1}^n K_{ij}} \sqrt{\sum_{i=1}^n K_{ij}}}. $$We can find the relation between $Q$ and $P$ in the next step.  3. Find eigenvectors and eigenvalues of the diffusion matrixThe eigenvalues of $P$ is equal to the eigenvalues of $Q$, which is arrange in descending power.  The right and left eigenvectors of $P$ are $$ \psi_k=e_k * e_1 $$ $$\phi_k=e_k / e_1 $$respectively, where $e_k$ is the corresponding eigenvector of the k-th eigenvalue of $Q$, $e_1$ is the eigenvector of the first eigenvalue of $Q$, and the multiplication and division is applied to each element, not the matrix operations.  4. Define a diffusion map and treat each row as a point coordinateLet $Y$ be a diffusion map, which map corrdinates between data and diffusion space [2], with its column$$ Y_i=[\frac{\lambda_i}{1-\lambda_i} \psi_i] $$Dimensional reduction is done by neglecting certain dimensions in the diffusion space, since the orthogonal left eigenvectors of $P$ form a basis for the diffusion space, and the associated coefficient indicate the importance of each dimension [2].  5. Cluster via k-meansUsing k-means method to get the labeling of the datas, here k is a manually select number of groups.  ImplementationPython code that implement diffusion maps using the steps above.  1. Import some libraries for programming<jupyter_code>%pylab inline<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>2. Define functions The functioning of get_K is to get the Guassian kernel matrix.<jupyter_code>def get_K(D, sigma):
    list1=[]

    for i in range(n):
        for j in range(n):
            list1.append( exp(-D[i][j]*D[i][j]/(4*sigma)) )
    K = array(list1).reshape(n, n)
    
    return K    #Guassian kernel matrix<jupyter_output><empty_output><jupyter_text>The functioning of get_Q is to get the matrix $Q$ we mention above.<jupyter_code>from scipy.sparse import csr_matrix

def get_Q(K):
    list1=[]
    list2=[]
    
    for i in range(n): list1.append( sqrt(sum(K[i])) )
    v = array(list1).reshape(n,1)    #v_i=(sum_j(K[i][j]))^(1/2)
    
    Q = K/(v*v.T)    #v*v.T is a n*n symmetric matrix.
                     #the multiplication and division is applied to each element,
                     #not the matrix operations.
    
    #turn Q into a sparse matrix to speed up calculation
    threshold = 5e-6
    for i in range(n):
        for j in range(n):
            if Q[i][j]>threshold: list2.append(1.)
            else: list2.append(0.)
    a = array(list2).reshape(n, n)
    Q = csr_matrix(Q*a, dtype=float)
    
    return Q    #Q is the matrix we define earlier, i.e., Q=D^(-1/2)KD^(-1/2)<jupyter_output><empty_output><jupyter_text>The functioning of get_eig is to get eigenvalues and eigenvectors of the diffusion matrix.<jupyter_code>from scipy.sparse.linalg import svds

def get_eig(Q, c):
    list1=[]
    list2=[]
    
    [U, S, V] = svds(Q, c+1)    #compute the largest c singular values/vectors for a sparse matrix
                               #since Q is a real symmetric matrix, singular values and non-negative
                               #eigenvalues coincide

    eigenvals = flip(S).reshape(c+1,1)
    
    for i in range(n): U[i]=flip(U[i])
    for i in range(n):
        for j in range(c+1):
            list1.append(U[i][j]/U[i][0])
            list2.append(U[i][j]*U[i][0])
            
    psi = array(list1).reshape(n, c+1)
    phi = array(list2).reshape(n, c+1)

    
    return [psi, phi, eigenvals]    #eigenvectors and eigenvalues of diff matrix<jupyter_output><empty_output><jupyter_text>The functioning of get_Y is to get the diffusion map.<jupyter_code>def get_Y(eigenvals, c, psi, t):
    list1=[]
    
    if t!=0:
        for i in range(n):
            for j in range(c):
                list1.append(psi[i][j+1]*pow(eigenvals[j+1], t))
        Y = array(list1).reshape(n, c)
    else:
        for i in range(n):
            for j in range(c):
                list1.append( psi[i][j+1]*(eigenvals[j+1]/(1-eigenvals[j+1])) )
        Y = array(list1).reshape(n, c)
        
    return Y    #Y is the diffusion map<jupyter_output><empty_output><jupyter_text>The functioning of diffusion_kmeans and distortionMinimization is to cluster via k-means.<jupyter_code>from math import ceil
import random as rnd

def diffusion_kmeans(Y, clusters, phi0, Niter, epsilon):
    
    d=len(Y[0])
    aD=Inf

    for i in range(Niter):
        tmp_ind = rnd.sample(range(n), clusters)
        
        list1=[]
        for p in tmp_ind:
            for q in range(d):
                list1.append(Y[p-1][q])
        c_0 = array(list1).reshape(clusters, d)
    
        idx, c, cindex, D, DX = distortionMinimization(Y, phi0, clusters, c_0, 0, epsilon)

        if D<aD:
            aD=D
            aDX=DX
            a_idx=idx
            ac=c
    D=aD
    DX=aDX
    idx=a_idx
    C=ac
    
    return idx, C, D, DX
def distortionMinimization(Y, phi0, clusters, c_0, DspFlag, epsilon):
    list1=[]
    list2=[]
    list3=[]
    
    d=len(Y[0])
    c=c_0
    oldD=Inf
    MaxIter=1000
    
    for i in range(MaxIter):
        for j in range(clusters):
            for p in range(n):
                for q in range(d):
                    list1.append(Y[p][q]-c[j][q])
            dX = (array(list1).reshape(n,d)).T    #dX=dX.T
            list1.clear()
            
            list2.append(sum(dX.conj()*dX,axis=0))
        DX = (array(list2).reshape(clusters,n)).T
        list2.clear()
        
        for j in range(n):
            list1.append(min(DX[j]))
            list2.append(argmin(DX[j])+1)
        Dtmp = array(list1).reshape(1,n)
        S = array(list2).reshape(n,1)
        
        list1.clear()
        list2.clear()
        
        ind=[]
        for j in range(clusters):
            for p in range(n):
                if S[p]==j+1: ind.append(p)
            if len(ind)==0:
                mx = max(Dtmp[0])
                m = argmax(Dtmp[0])
                S[m][0] = j+1
                Dtmp[0][m] = 0
            ind.clear()

        for j in range(clusters):
            for p in range(n):
                if S[p]==j+1: ind.append(p)
            for p in range(d):
                tempsum=0
                tempsum2=0
                for q in ind: tempsum += phi0[q][0]
                for q in ind: tempsum2 += phi0[q][0]*Y[q][p]/tempsum
                list3.append(tempsum2)
            ind.clear()
        c = array(list3).reshape(clusters,d)
        
        list3.clear()
        
        D = dot(Dtmp, phi0)
        
        
        if (oldD-D)/D < epsilon: break
        oldD=D
        
        
    for j in range(clusters):
        for p in range(n):
            for q in range(d):
                list1.append(Y[p][q]-c[j][q])
        dX = (array(list1).reshape(n,d)).T    #dX=dX.T
        list1.clear()
        list2.append(sum(dX.conj()*dX,axis=0))
        DX = array(list2).reshape(n,1)
        list2.clear()
            
        dummy = min(DX.T[0])
        list3.append( argmin(DX.T[0]) )
    cindex = array(list3).reshape(clusters,1)
    
    
    return S, c, cindex, D, DX<jupyter_output><empty_output><jupyter_text>3. Load data<jupyter_code>from scipy.io import loadmat

ori = loadmat('Data.mat')

n = len(ori['data'])       #number of points
p = len(ori['data'][0])    #dimension of points

for i in range(n):
    plt.scatter(ori['data'][i][0],ori['data'][i][1], s=2, c="blue")
plt.title("Original data")
plt.gcf().set_size_inches(6,6)
plt.show()<jupyter_output><empty_output><jupyter_text>4. Set parameters and do the calculation  4.1. Set parameters for constructing a diffusion mapD : a n*n matrix with pairwise Euclidean distance as its entries  sigma : scaling parameter $\sigma$ in the Guassian kernel function  c : the largest $c$ eigenvalues, that is, the largest dimension we consider  t : 'time'By taking powers of $P$, we increase the number of steps taken, that is, $P_{ij}^{t}$ sum all paths that take $t$ steps from $x_i$ to $x_j$, this is called the diffusion process.  When increasing value $t$, the probability of following a path along the underlying geometric structure of the data set increases, since along the geometric structure, points are dense and therefore highly connected, pathways form along short and high probability jumps [2].<jupyter_code>from scipy.spatial.distance import squareform, pdist
D = squareform(pdist(ori['data']))
c = 10
t = 0
flag_t = 0<jupyter_output><empty_output><jupyter_text>Ann Lee guess the correct sigma by considering the distance between one point to its k-nearest neighbor point.<jupyter_code>#guess the correct sigma
D_sort = sort(D,1)
k = 30       #k-nearest point of one point

dist_k=[]    #distance between one point and its k-nearest neighbor point
for i in range(len(D_sort)): dist_k.append( D_sort[i][k+1] )
    
sigma = median(dist_k)*median(dist_k)/2    #guessed sigma
print("suggest sigma = ", sigma)<jupyter_output>suggest sigma =  0.00436767844472376<jupyter_text>Manually pick one sigma.<jupyter_code>sigma = 0.00025<jupyter_output><empty_output><jupyter_text>4.2. Compute the eigenvalues and eigenvectors of diffusion matrix and get the diffusion map<jupyter_code>def diffuse(D, sigma, c, t):
    
    K = get_K(D, sigma)
    Q = get_Q(K)
    [psi, phi, eigenvals] = get_eig(Q, c)
    Y = get_Y(eigenvals, c, psi, t)

    return Y, eigenvals, psi, phi
Y, eigenvals, psi, phi = diffuse(D, sigma, c, t)<jupyter_output><empty_output><jupyter_text>4.3. Set parameters for clustering via k-meansclusters : number of clusters  Niter : number of iterations to repeat the clustering with new IC  epsilon : relative distortion  phi0 : first left eigenvector of P<jupyter_code>clusters=3
Niter=100
epsilon=1e-3

list1=[]
for i in range(n): list1.append(phi[i][0])
phi0=array(list1).reshape(n,1)<jupyter_output><empty_output><jupyter_text>4.4. Diffusion k-means clustering<jupyter_code>idx, C, ERR, DX = diffusion_kmeans(Y, clusters, phi0, Niter, epsilon)<jupyter_output><empty_output><jupyter_text>5. Display the result  5.1. Show the fall-off of eigenvalues of diffusion matrix<jupyter_code>plt.plot(range(1,c+1),eigenvals[1:c+1])
plt.title('Fall-off of non-trivial eigenvalues')
plt.ylabel('$\lambda$')
plt.gcf().set_size_inches(6,3)
plt.show()

if flag_t:
    plt.plot(range(1,c+1), pow(eigenvals[1:c+1], t))
    plt.ylabel('$\lambda$^t')
else:
    plt.plot(range(1,c+1), eigenvals[1:c+1]/(1-eigenvals[1:c+1]) )
    ylabel('$\lambda$ / (1-$\lambda$)')
plt.gcf().set_size_inches(6,3)
plt.show()<jupyter_output><empty_output><jupyter_text>5.2. Show the diffusion space<jupyter_code>from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = Axes3D(fig)
for i in range(n): ax.scatter(Y[i][0],Y[i][1],Y[i][2],s=2, c='b')
ax.set_title('Embedding with first 3 diffusion coordinates')
plt.show()<jupyter_output><empty_output><jupyter_text>5.3. Show the labeling<jupyter_code>cmap = cm.brg
norm = Normalize(1, clusters)

for i in range(n):
    plt.scatter(ori['data'][i][0], ori['data'][i][1], s=2, c=cmap(norm(idx[i])) )
plt.title('sigma=0.00025')
plt.gcf().set_size_inches(6,6)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12581818181818183, 'fraction_numerical': 0.02101010101010101, 'mean_word_length': 3.006474587245063, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 65, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '317920', 'license': 'MIT', 'n_tokens_mistral': 4446, 'n_tokens_neox': 4243, 'n_words': 1260, 'path': 'python/Diff_Map_1.ipynb', 'repo_name': 'teshenglin/diffusion_maps', 'size': '12378'}