# Code from Chapter 6 of Machine Learning: An Algorithmic Perspective (2nd Edition) # by Stephen Marsland (http://stephenmonika.net) # You are free to use, change, or redistribute the code in any way you wish for # non-commercial purposes, but please maintain the name of the original author. # This code comes with no warranty of any kind. # Stephen Marsland, 2008, 2014 # The Locally Linear Embedding algorithm, and the swissroll example import pylab as pl import numpy as np def swissroll(): # Make the swiss roll dataset N = 1000 noise = 0.05 t = 3*np.pi/2 * (1 + 2*np.random.rand(1,N)) h = 21 * np.random.rand(1,N) data = np.concatenate((t*np.cos(t),h,t*np.sin(t))) + noise*np.random.randn(3,N) return np.transpose(data), np.squeeze(t) def lle(data,nRedDim=2,K=12): ndata = np.shape(data)[0] ndim = np.shape(data)[1] d = np.zeros((ndata,ndata),dtype=float) # Inefficient -- not matrices for i in range(ndata): for j in range(i+1,ndata): for k in range(ndim): d[i,j] += (data[i,k] - data[j,k])**2 d[i,j] = np.sqrt(d[i,j]) d[j,i] = d[i,j] indices = d.argsort(axis=1) neighbours = indices[:,1:K+1] W = np.zeros((K,ndata),dtype=float) for i in range(ndata): Z = data[neighbours[i,:],:] - np.kron(np.ones((K,1)),data[i,:]) C = np.dot(Z,np.transpose(Z)) C = C+np.identity(K)*1e-3*np.trace(C) W[:,i] = np.transpose(np.linalg.solve(C,np.ones((K,1)))) W[:,i] = W[:,i]/np.sum(W[:,i]) M = np.eye(ndata,dtype=float) for i in range(ndata): w = np.transpose(np.ones((1,np.shape(W)[0]))*np.transpose(W[:,i])) j = neighbours[i,:] #print shape(w), np.shape(np.dot(w,np.transpose(w))), np.shape(M[i,j]) ww = np.dot(w,np.transpose(w)) for k in range(K): M[i,j[k]] -= w[k] M[j[k],i] -= w[k] for l in range(K): M[j[k],j[l]] += ww[k,l] evals,evecs = np.linalg.eig(M) ind = np.argsort(evals) y = evecs[:,ind[1:nRedDim+1]]*np.sqrt(ndata) return evals,evecs,y data,t = swissroll() evals,evecs,y = lle(data) t -= t.min() t /= t.max() pl.scatter(y[:,0],y[:,1],s=50,c=t,cmap=pl.cm.gray) pl.axis('off') pl.show()