The implemented algorithm is described in more detail in [1,3].
This algorithm is based on the quantum bit (qubit) model. A qubit can have any value between 0 and 1 (superposition property) until it is observed, which is when the system collapses to either state. However, the probability with which the system collapses to either state may be different. The superposition property or linear combination of states can be expressed as
$$ [\psi] = \alpha[0] + \beta[1] $$where $\psi$ is an arbitrary state vector and $\alpha$, $\beta$ are the the probability amplitude coefficients of basis states $[0]$ and $[1]$, respectevely. The basis states correspond to the spin of the modeled particle (in this case, a ferminion, e.g. electron). The coefficients are subjected to the following normalization:
$$|\alpha|^2 + |\beta|^2 = 1$$where $|\alpha|^2$, $|\beta|^2$ are the probabilities of observing states $[0]$ and $[1]$, respectevely. $\alpha$ and $\beta$ are complex quantities and represent a qubit:
$$\begin{bmatrix} \alpha \\ \beta \end{bmatrix}$$Moreover, a qubit string may be represented by: $$ \begin{bmatrix} \left.\begin{matrix} \alpha_1\\ \beta_1 \end{matrix}\right| & \left.\begin{matrix} \alpha_2\\ \beta_2 \end{matrix}\right| & \begin{matrix} \alpha_3\\ \beta_3 \end{matrix} \end{bmatrix} $$
The probability of observing the state $[000]$ will be $|\alpha_1|^2 \times |\alpha_2|^2 \times |\alpha_3|^2$
To use this model for computing purposes, black-box objects called oracles are used. Oracles contain strings of qubits and generate their own input by observing the state of the qubits. After collapsing, the qubit value becomes analog to a classical bit. Each string of qubits represents a number, so the number of qubits in each string will define its precision. The number of strings chosen for the oracles depends on the number of clusters and dimensionality of the problem (e.g. for 3 clusters of 2 dimensions, 6 strings will be used since 6 numbers are required). Each oracle will represent a possible solution.
The algorithm has the following steps:
The oracles are created in this step and all qubit coefficients are initialized with $\frac{1}{\sqrt{2}}$, so that the system will observe either state with equal probability.
Collapsing the oracles implies making an observation of each qubit of each qubit string in each oracle. This is done by first choosing a coefficient to use (which is irrelevant), e.g. $\alpha$. Then, we generate a random value $r$ between 0 and 1. If $\alpha \ge r$ then the system collapses to $[0]$, otherwise to $[1]$.
In this step we convert the binary representation of the qubit strings to base 10 and use them those values as initial centroids for K-Means. For each oracle, classical K-Means is then executed until it stabilizes or reaches the iteration limit. The solution centroids are returned to the oracles in binary representation.
Cluster fitness is computed using the Davies-Bouldin index for each oracle. The score of each oracle is stored in the oracle itself.
The best scoring oracle is stored.
So far, we've had classical K-Means with a complex random number generation for the centroids and complicated datastructures. This is the step that fundamentally differs from the classical version. In this step a quantum gate (in this case a rotation gate) is applied to all oracles except the best one. The basic idea is to shift the qubit coefficients of the least scoring oracles so they'll have a higher probability of collapsing into initial centroid values closer to the best solution so far. This way, in future generations, we'll not initiate with the best centroids so far (which will not converge further into a better solution) but we'll be closer while still ensuring diversity (which is also a desired property of the genetic computing paradigm). In conclusion, we want to look for better solutions than the one we got before in each oracle while moving in the direction of the best we found so far.
In the original formulation of this algorithm two extra step existed to further increase diversity: quantum crossover and quantum mutation inversion. Both are part of the genetic algorithms toolbox, but were not implemented due to the suggestion from [1] that they are unnecessary steps with the careful choice of the rotatin angle.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import numpy as np
import sys
import os
import pickle
from datetime import datetime
from sklearn.cluster import KMeans
This should be rerun every time the code changes.
In [245]:
# add path for custom modules
mypath=os.path.abspath('/home/chiroptera/workspace/thesis/quantum k-means/implementation')
if not mypath in sys.path:
sys.path.insert(1, mypath)
del mypath
import oracle
import qubitLib
import DaviesBouldin
reload(DaviesBouldin)
reload(oracle)
reload(qubitLib)
Out[245]:
To test the algorithm, a mixture of 4 Gaussian distribution is generated.
In [20]:
class gaussMix:
def __init__(self,**kwargs):
if kwargs is None:
raise Exception('Some arguments must be supplied. dim=[dimension] \\
or pack=[dictionary with a previously created package]')
else:
if 'dim' in kwargs:
self.dim=kwargs['dim']
self.idNum=0
self.gaussList=list()
self.cardinality=0
elif 'pack' in kwargs:
self.unpack(kwargs['pack'])
else:
raise Exception('dim=[dimension] or pack=[dictionary with a previously created package]')
def unpack(self,pack):
self.dim=pack['dim']
self.gaussList=pack['all']
self.idNum=len(self.gaussList)
self.cardinality=pack['cardinality']
def add(self,mean,variance,numPoints):
if len(mean)!=self.dim:
raise Exception("wrong dimension")
self.cardinality += numPoints
g={'id':self.idNum,'mean':mean,'var':variance,
'card':numPoints,'dist':np.random.normal(mean,variance,(numPoints,self.dim))}
"""
if self.mixture== None:
self.mixture=g['dist']
else:
self.mixture=np.concatenate((self.mixture,g['dist']))
"""
self.gaussList.append(g)
self.updateWeights()
self.idNum+=1
def getMixture(self):
returnMix=[]
for i,g in enumerate(self.gaussList):
returnMix.append(g['dist'])
returnMix=np.concatenate(tuple(returnMix))
return returnMix
def updateWeights(self):
for g in self.gaussList:
g['weight']=float(g['card'])/self.cardinality
def sample(self,numPoints):
returnMix=[]
if numPoints > self.cardinality:
return self.getMixture()
for i,g in enumerate(self.gaussList):
num=int(numPoints*g['weight'])
returnMix.append(g['dist'][0:num,:])
returnMix=np.concatenate(tuple(returnMix))
return returnMix
def getWeights(self):
w=[0]*len(self.gaussList)
for i,g in enumerate(self.gaussList):
w[i]=g['weight']
return w
def package(self):
pack=dict()
pack['all']=self.gaussList
pack['data']=self.getMixture()
pack['dim']=self.dim
pack['cardinality']=self.cardinality
return pack
def gaussAsList(self):
m=list()
for i,g in enumerate(self.gaussList):
m.append(g['dist'])
return m
"""
# better function without redundacy
def package(self):
pack=dict()
pack['all']=[None]*len(self.gaussList)
pack['all']
for i,g in enumerate(self.gaussList):
pack['all']=self.gaussList
pack['data']=self.getMixture()
pack['dim']=self.dim
pack['cardinality']=self.cardinality
return pack
"""
In [24]:
name='bi36'
m=gaussMix(dim=2)
m.add((-5,3),[0.25,2],200)
m.add((2,-2),(0.25,0.25),400)
m.add((2,10),[0.75,0.25],300)
m.add((-3,2),[0.3,0.3],300)
m.add((3,7),[0.3,0.3],300)
m.add((-2,7),[0.3,1],300)
for i,g in enumerate(m.gaussAsList()):
plt.title(name)
plt.plot(g[:,0],g[:,1],'.')
plt.savefig("figure.png")
In [253]:
output = open(name+'.pkl', 'wb')
pickle.dump(m.package(), output)
output.close()
!ls *pkl
In [14]:
name='bi46'
m=gaussMix(2)
m.add((-5,3),[0.25,2],200)
m.add((2,-2),(0.25,0.25),400)
m.add((2,10),[0.75,0.25],300)
m.add((-3,2),[0.3,0.3],300)
m.add((3,7),[0.3,0.3],300)
m.add((-2,7),[0.3,1],300)
s=m.getMixture()
plt.title(name)
plt.plot(s[:,0],s[:,1],'.')
Out[14]:
In [5]:
cd workspace/QCThesis/QuantumGeneticKMeans/datasets/
In [6]:
output = open(name+'.pkl', 'wb')
pickle.dump(m.package(), output)
output.close()
!ls *pkl
In [271]:
name='tri36'
m=gaussMix(3)
m.add((-5,3,3),[0.25,2,1],200)
m.add((2,-2,4),(0.25,0.25,0.25),400)
m.add((2,10,5),[0.75,0.25,0.5],300)
m.add((-3,2,3),[0.3,0.3,0.3],300)
m.add((3,7,4),[0.3,0.3,0.3],300)
m.add((-2,7,5),[0.3,1,0.6],300)
fig = plt.figure()
ax = fig.gca(projection='3d')
# plot points in 3D
for i,g in enumerate(m.gaussList):
s=g['dist']
ax.plot(s[:,0],s[:,1],s[:,2],'.')
plt.title(name)
#ax.view_init(elev=10., azim=45)
Out[271]:
In [272]:
output = open(name+'.pkl', 'wb')
pickle.dump(m.package(), output)
output.close()
!du -sh *pkl
In [269]:
name='tri56'
m=gaussMix(3)
m.add((-5,3,3),[0.25,2,1],20000)
m.add((2,-2,4),(0.25,0.25,0.25),40000)
m.add((2,10,5),[0.75,0.25,0.5],30000)
m.add((-3,2,3),[0.3,0.3,0.3],30000)
m.add((3,7,4),[0.3,0.3,0.3],30000)
m.add((-2,7,5),[0.3,1,0.6],30000)
fig = plt.figure()
ax = fig.gca(projection='3d')
# plot points in 3D
for i,g in enumerate(m.gaussList):
s=g['dist']
ax.plot(s[:,0],s[:,1],s[:,2],'.')
plt.title(name)
#ax.view_init(elev=10., azim=45)
Out[269]:
In [270]:
output = open(name+'.pkl', 'wb')
pickle.dump(m.package(), output)
output.close()
!du -sh *pkl
In [7]:
name='sex46'
m=gaussMix(6)
m.add((-5,3,3,3,3,3),[0.25,2,1,0.25,0.25,1],2000)
m.add((2,-2,4,4,4,4),(0.25,0.25,0.25,0.25,0.25,0.25),4000)
m.add((2,10,5,5,5,5),[0.75,0.25,0.5,0.75,0.25,0.5],3000)
m.add((-3,2,3,3,3,3),[0.3,0.3,0.3,0.3,0.3,0.3],3000)
m.add((3,7,4,4,4,4),[0.3,0.3,0.3,0.3,0.3,0.3],3000)
m.add((-2,7,5,5,5,5),[0.3,1,0.6,0.3,1,0.6],3000)
"""
fig = plt.figure()
ax = fig.gca(projection='3d')
# plot points in 3D
for i,g in enumerate(m.gaussList):
s=g['dist']
ax.plot(s[:,0],s[:,1],s[:,2],'.')
plt.title(name)"""
#ax.view_init(elev=10., azim=45)
Out[7]:
In [8]:
output = open(name+'.pkl', 'wb')
pickle.dump(m.package(), output)
output.close()
!du -sh *pkl
In [276]:
name='sex56'
m=gaussMix(6)
m.add((-5,3,3,3,3,3),[0.25,2,1,0.25,0.25,1],20000)
m.add((2,-2,4,4,4,4),(0.25,0.25,0.25,0.25,0.25,0.25),40000)
m.add((2,10,5,5,5,5),[0.75,0.25,0.5,0.75,0.25,0.5],30000)
m.add((-3,2,3,3,3,3),[0.3,0.3,0.3,0.3,0.3,0.3],30000)
m.add((3,7,4,4,4,4),[0.3,0.3,0.3,0.3,0.3,0.3],30000)
m.add((-2,7,5,5,5,5),[0.3,1,0.6,0.3,1,0.6],30000)
"""
fig = plt.figure()
ax = fig.gca(projection='3d')
# plot points in 3D
for i,g in enumerate(m.gaussList):
s=g['dist']
ax.plot(s[:,0],s[:,1],s[:,2],'.')
plt.title(name)"""
#ax.view_init(elev=10., azim=45)
Out[276]:
In [277]:
output = open(name+'.pkl', 'wb')
pickle.dump(m.package(), output)
output.close()
!du -sh *pkl
In [22]:
200+400+300+300+300+300
Out[22]:
In [12]:
name='sex310'
m=gaussMix(10)
m.add((-5,3,3,3,3,3,3,3,3,3),0.6,200)
m.add((2,-2,4,4,4,4,4,4,4,4),0.25,400)
m.add((2,10,5,5,5,5,5,5,5,5),0.4,300)
m.add((-3,2,3,3,3,3,3,3,3,3),0.3,300)
m.add((3,7,4,4,4,4,4,4,4,4),0.3,300)
m.add((-2,7,5,5,5,5,5,5,5,5),0.6,300)
"""
fig = plt.figure()
ax = fig.gca(projection='3d')
# plot points in 3D
for i,g in enumerate(m.gaussList):
s=g['dist']
ax.plot(s[:,0],s[:,1],s[:,2],'.')
plt.title(name)"""
#ax.view_init(elev=10., azim=45)
Out[12]:
In [13]:
output = open(name+'.pkl', 'wb')
pickle.dump(m.package(), output)
output.close()
!du -sh *pkl
In [100]:
# generate gaussians
numPoints=1000
numGaussians=6
pointsPerGaussian=int(numPoints/numGaussians)
gMix=list()
gMix.append(np.random.normal((-5,3),[0.25,2],(pointsPerGaussian,2)))
gMix.append(np.random.normal((2,-2),[0.25,0.25],(pointsPerGaussian,2)))
gMix.append(np.random.normal((2,10),[0.75,0.25],(pointsPerGaussian,2)))
gMix.append(np.random.normal((-3,2),[0.3,0.3],(pointsPerGaussian,2)))
gMix.append(np.random.normal((3,7),[0.3,0.3],(pointsPerGaussian,2)))
gMix.append(np.random.normal((-2,7),[0.3,1],(pointsPerGaussian,2)))
mixture=np.concatenate(tuple(gMix))
# add outliers
# plot data points
for i in range(0,len(gMix)):
plt.plot(gMix[i][:,0],gMix[i][:,1],'.')
In [77]:
output = open('data.pkl', 'wb')
pickle.dump(mixture, output)
output.close()
In [ ]:
In [7]:
## Initialization step
numClusters=15
numOracles=5
qubitStringLen=5
qGenerations=100
In [241]:
%load /home/chiroptera/workspace/QCThesis/QuantumGeneticKMeans/QK-Means.py
In [8]:
# needs:
# - mixture
# - numOracles
# - numClusters
# - qubitStringLen
# - qGenerations
# -
# returns:
# - qk_timings_cg
# - qk_centroids (centroids of oracles from last generation)
# - qk_assignment(assignment of oracles from last generation)
# - fitnessEvolution
# - qk_total
# - qk_genTimes
# - qk_assignedData
def qk_means(mixture,numOracles,numClusters,qubitStringLen,qGenerations):
fitnessEvolution = np.zeros((qGenerations,numOracles+1))
qk_timings_cg = list()
start = datetime.now()
best = 0 #index of best oracle (starts at 0)
oras = list()
qk_centroids = [0]*numOracles
qk_estimator = [0]*numOracles
qk_assignment = [0]*numOracles
for i in range(0,numOracles):
oras.append(oracle.Oracle())
oras[i].initialization(numClusters*2,qubitStringLen)
oras[i].collapse()
qk_timings_cg.append((datetime.now() - start).total_seconds())
start = datetime.now()
for qGen_ in range(0,qGenerations):
## Clustering step
for i,ora in enumerate(oras):
if qGen_ != 0 and i == best: # current best shouldn't be modified
continue
qk_centroids[i] = np.vstack(np.hsplit(ora.getIntArrays(),numClusters))
qk_estimator[i] = KMeans(n_clusters=numClusters,init=qk_centroids[i],n_init=1)
qk_assignment[i] = qk_estimator[i].fit_predict(mixture)
qk_centroids[i] = qk_estimator[i].cluster_centers_
ora.setIntArrays(np.concatenate(qk_centroids[i]))
## Compute fitness
score = DaviesBouldin.DaviesBouldin(mixture,qk_centroids[i],qk_assignment[i])
ora.score = score.eval()
## Store best from this generation
for i in range(1,numOracles):
if oras[i].score < oras[best].score:
best = i
## Quantum Rotation Gate
for i in range(0,numOracles):
if i == best:
continue
oras[i].QuantumGateStep(oras[best])
## Collapse qubits
oras[i].collapse()
qk_timings_cg.append((datetime.now() - start).total_seconds())
for i in range(0,numOracles):
fitnessEvolution[qGen_,i]=oras[i].score
fitnessEvolution[qGen_,-1]=best
print '.', # simple "progress bar"
start = datetime.now()
print 'Done!'
return qk_centroids,qk_assignment,fitnessEvolution,qk_timings_cg
# - qk_timings_cg
# - qk_centroids (centroids of oracles from last generation)
# - qk_assignment(assignment of oracles from last generation)
# - fitnessEvolution
# - qk_total
# - qk_genTimes
# - qk_assignedData
In [9]:
numRounds=5
qk_results=list()
for i in range(numRounds):
print 'round ',i
qk_centroids,qk_assignment,fitnessEvolution,qk_timings_cg=qk_means(mixture,numOracles,numClusters,qubitStringLen,qGenerations)
qk_results.append([qk_centroids,qk_assignment,fitnessEvolution,qk_timings_cg])
Organize data in dictionary.
In [10]:
qk_rounds=dict()
qk_rounds['centroids']=list()
qk_rounds['assignment']=list()
qk_rounds['fitness']=list()
qk_rounds['times']=list()
for i in range(numRounds):
qk_rounds['centroids'].append(qk_results[i][0])
qk_rounds['assignment'].append(qk_results[i][1])
qk_rounds['fitness'].append(qk_results[i][2])
qk_rounds['times'].append(qk_results[i][3])
print 'Done!'
In [71]:
len(qk_rounds['centroids'][0]),type(qk_rounds['centroids'][0])
Out[71]:
Assign points to clusters.
In [12]:
qk_rounds['assignedData']=list() #assigned data for the best solution in each round
for i in range(numRounds):
# assign clusters to data
best=int(qk_rounds['fitness'][i][-1,-1])
qk_assignment=qk_rounds['assignment'][i]
qk_assignedData = [None]*numClusters
for i,j in enumerate(qk_assignment[best]):
if qk_assignedData[j] != None:
qk_assignedData[j] = np.vstack((qk_assignedData[j],mixture[i]))
else:
qk_assignedData[j] = mixture[i]
qk_rounds['assignedData'].append(qk_assignedData)
print 'Done!'
In [67]:
qk_rounds['assignedData'][0]
Out[67]:
In [13]:
qk_rounds['total time'] = list()
for i in range(numRounds):
qk_times=qk_rounds['times'][i]
qk_total = np.sum(np.array(qk_times))
qk_rounds['total time'].append(qk_total)
qk_time_mean = np.mean(np.array(qk_rounds['total time']))
qk_time_var = np.var(np.array(qk_rounds['total time']))
qk_time_best = np.min(np.array(qk_rounds['total time']))
qk_time_worst = np.max(np.array(qk_rounds['total time']))
print 'Total execution times'
print 'mean\t\tvariance\tbest\t\tworst'
print qk_time_mean,'\t',qk_time_var,'\t',qk_time_best,'\t',qk_time_worst
In [14]:
qk_rounds['best evolution'] = list()
qk_rounds['pop variance'] = list()
qk_rounds['pop mean'] = list()
for i in range(numRounds):
bestSeq=[]
bestSeqIndex=[]
fitnessEvolution=qk_rounds['fitness'][i]
for i in range(0,fitnessEvolution.shape[0]):
bestSeq.append(fitnessEvolution[i,fitnessEvolution[i,-1]])
bestSeq=np.array(bestSeq)
genVar=np.zeros(qGenerations)
genMean=np.zeros(qGenerations)
for i,ar in enumerate(fitnessEvolution):
genVar[i]=np.var(ar[:-1])
genMean[i]=np.mean(ar[:-1])
qk_rounds['best evolution'].append(bestSeq)
qk_rounds['pop variance'].append(genVar)
qk_rounds['pop mean'].append(genMean)
In [15]:
for i in range(numRounds):
plt.plot(range(qGenerations),qk_rounds['best evolution'][i],label="round "+str(i))
plt.legend(loc='best', framealpha=0.3,prop={'size':'small'})
Out[15]:
In [14]:
for i in range(numRounds):
plt.plot(range(qGenerations),qk_rounds['pop mean'][i],label="round "+str(i))
plt.legend(loc='best', framealpha=0.3,prop={'size':'small'})
Out[14]:
In [15]:
for i in range(numRounds):
plt.plot(range(qGenerations),qk_rounds['pop variance'][i],label="round "+str(i))
plt.legend(loc='best', framealpha=0.3,prop={'size':'small'})
Out[15]:
Save QK centroids. Each row has the format:
[number of round, number of oracle, centroids in line]
In [16]:
qk_saveCentroids=list()
for i in range(numRounds):
for j in range(numOracles):
line=np.concatenate(qk_rounds['centroids'][i][j]).tolist()
line.insert(0,j)
line.insert(0,i)
qk_saveCentroids.append(line)
qk_saveCentroids=np.array(qk_saveCentroids)
qk_saveCentroids.shape
Out[16]:
To make the run of classic K-Means approximate to that of QK-Means, we'll set the number of initializations to $qGenerations \times numOracles$. This is chosen since the classical K-Means is executed once for each oracle in each generation.
In [16]:
def k_means(mixture,numClusters,numInits):
k_timings_cg=list()
start=datetime.now()
k_assignment=list()
k_centroids=list()
for i in range(numInits):
estimator = KMeans(n_clusters=numClusters,init='k-means++',n_init=1)
assignment = estimator.fit_predict(mixture)
centroids = estimator.cluster_centers_
k_centroids.append(centroids)
k_assignment.append(assignment)
k_timings_cg.append((datetime.now() - start).total_seconds())
start=datetime.now()
return k_centroids,k_assignment,k_timings_cg
In [17]:
initsPercentage=1.15
numInits=np.int(qGenerations*numOracles*initsPercentage)
k_results=list()
for i in range(numRounds):
print 'round ',i
k_centroids,k_assignment,k_timings=k_means(mixture,numClusters,numInits)
k_results.append([k_centroids,k_assignment,k_timings])
print 'done!'
In [18]:
k_rounds=dict()
k_rounds['centroids']=list()
k_rounds['assignment']=list()
k_rounds['times']=list()
for i in range(numRounds):
k_rounds['centroids'].append(k_results[i][0])
k_rounds['assignment'].append(k_results[i][1])
k_rounds['times'].append(k_results[i][2])
print 'Done!'
In [62]:
k_rounds['assignedData']=list() #assigned data for the best solution in each round
for i in range(numRounds):
# assign clusters to data
best=int(qk_rounds['fitness'][i][-1,-1])
k_assignment=k_rounds['assignment'][i]
k_assignedData = [None]*numClusters
for i,j in enumerate(k_assignment[best]):
if k_assignedData[j] != None:
k_assignedData[j] = np.vstack((k_assignedData[j],mixture[i]))
else:
k_assignedData[j] = mixture[i]
k_rounds['assignedData'].append(k_assignedData)
print 'Done!'
In [63]:
k_rounds['assignedData'][0]
Out[63]:
In [20]:
k_rounds['total time'] = list()
for i in range(numRounds):
k_times=k_rounds['times'][i]
k_total = np.sum(np.array(k_times))
k_rounds['total time'].append(k_total)
k_time_mean = np.mean(np.array(k_rounds['total time']))
k_time_var = np.var(np.array(k_rounds['total time']))
k_rounds['fastest round']=np.argmin(np.array(k_rounds['total time']))
k_rounds['slowest round']=np.argmax(np.array(k_rounds['total time']))
k_time_best = k_rounds['total time'][k_rounds['fastest round']]
k_time_worst = k_rounds['total time'][k_rounds['slowest round']]
print 'Total execution times'
print 'mean\t\tvariance\tbest\t\tworst'
print k_time_mean,'\t',k_time_var,'\t',k_time_best,'\t',k_time_worst
In [60]:
# assign data to clusters
k_rounds['assignedData']=[None]*numRounds #assigned data for the best solution in each round
for i in range(numRounds):
# assign clusters to data
k_assignment=k_rounds['assignment'][i][k_rounds['fitness best index'][i]]
k_assignedData = [None]*numClusters
for point,j in enumerate(k_assignment):
if k_assignedData[j] != None:
k_assignedData[j] = np.vstack((k_assignedData[j],mixture[point]))
else:
k_assignedData[j] = mixture[i]
k_rounds['assignedData'][i]=k_assignedData
In [64]:
k_rounds['assignedData'][0]
Out[64]:
In [38]:
k_rounds['total time'] = list()
for i in range(numRounds):
k_times=k_rounds['times'][i]
k_total = np.sum(np.array(k_times))
k_rounds['total time'].append(k_total)
k_rounds['time mean'] = np.mean(np.array(k_rounds['total time']))
k_rounds['time variance'] = np.var(np.array(k_rounds['total time']))
k_rounds['time best']=np.min(np.array(k_rounds['total time']))
k_rounds['time worst']=np.max(np.array(k_rounds['total time']))
print 'Total execution times'
print 'mean\t\tvariance\t\tbest\t\tworst'
print k_rounds['time mean'],'\t',k_rounds['time variance'],'\t',k_rounds['time best'],'\t',k_rounds['time worst']
In [30]:
k_score=DaviesBouldin.DaviesBouldin(mixture,k_rounds['centroids'][i][j],k_rounds['assignment'][i][j])
print k_score.eval()
print numInits
print type(k_score.eval())
print numInits*64/1024
In [52]:
#compute fitnesses
k_rounds['fitness']=[np.inf]*numRounds
k_rounds['fitness best index']=[None]*numRounds
k_bestScore=np.inf
for i in range(numRounds):
print 'round ',i
init_scores=list()
k_best=np.inf
for j in range(numInits):
k_score=DaviesBouldin.DaviesBouldin(mixture,k_rounds['centroids'][i][j],k_rounds['assignment'][i][j])
s=k_score.eval()
if s<k_rounds['fitness'][i]:
k_rounds['fitness'][i]=s
k_rounds['fitness best index'][i]=j
In [55]:
print k_rounds['fitness']
print k_rounds['fitness best index']
In [21]:
k_rounds['fitness']=list()
start=datetime.now()
k_bestScore=np.inf
for i in range(numRounds):
print 'round ',i
init_scores=list()
for j in range(numInits):
# compute cluster fitness
k_score=DaviesBouldin.DaviesBouldin(mixture,k_rounds['centroids'][i][j],k_rounds['assignment'][i][j])
init_scores.append(k_score.eval())
k_rounds['fitness'].append(np.array(init_scores))
m=np.min(init_scores)
if m < k_bestScore:
k_bestScore=m
k_rounds['best fit round']=i
print 'done!',(datetime.now() - start).total_seconds()
In [48]:
k_rounds['fitness best']=list()
k_rounds['fitness worst']=list()
k_rounds['fitness mean']=list()
k_rounds['fitness variance']=list()
for fit in k_rounds['fitness']:
k_rounds['fitness best'].append(np.min(fit))
k_rounds['fitness worst'].append(np.max(fit))
k_rounds['fitness mean'].append(np.mean(fit))
k_rounds['fitness variance'].append(np.var(fit))
for i in range(len(k_rounds['fitness'])):
print k_rounds['fitness best'][i],k_rounds['fitness worst'][i],k_rounds['fitness mean'][i],k_rounds['fitness variance'][i]
print np.mean(k_rounds['fitness best']),np.mean(k_rounds['fitness worst']),np.mean(k_rounds['fitness mean']),np.mean(k_rounds['fitness variance'])
In [25]:
k_rounds.keys()
#choose only best/worst assignment
#choose only best/worst centroid
Out[25]:
In [30]:
# save
output = open('k_ass.pkl', 'wb')
pickle.dump(k_rounds['centroids'], output)
output.close()
In [54]:
# print DB scores
print "Davies-Bouldin cluster fitness (less is better)"
qk_best=np.inf
for i in range(numRounds):
m=np.min(qk_rounds['best evolution'][i])
if m < qk_best:
qk_best = m
print "QK-Means cluster fitness (DB score):\t",qk_best
print "K-Means cluster fitness (DB score):\t",np.min(k_rounds['fitness'])
In [182]:
print 'Total execution times'
print '\t\tmean\t\tvariance\t\tbest\t\tworst'
print 'K-Means:\t',k_time_mean,'\t',k_time_var,'\t',k_time_best,'\t',k_time_worst
print 'QK-Means:\t',qk_time_mean,'\t',qk_time_var,'\t\t',qk_time_best,'\t',qk_time_worst
print "\nQK-Means / K-Means ratio:\t",qk_time_mean/k_time_mean
In [ ]:
# display results
fig1=plt.figure()
ax1 = fig1.add_subplot(3,2,1, adjustable='box', aspect=0.5, position=[0.1,0.1, 1, 1])
ax2 = fig1.add_subplot(3,2,2,adjustable='box', aspect=0.5,position=[1.2,0.1, 1, 1])
ax3 = fig1.add_subplot(3,2,3,adjustable='box', aspect=0.5, position=[0.1,-1.2, 1, 1])
ax1.set_title("Gaussian mixture")
for i in range(0,len(gMix)):
ax1.plot(gMix[i][:,0],gMix[i][:,1],'.')
ax2.set_title("QK-Means")
for i in range(0,numClusters):
ax2.plot(qk_assignedData[i][:,0],qk_assignedData[i][:,1],'.')
ax2.plot(qk_centroids[best][i][0],qk_centroids[best][i][1],'ko')
ax3.set_title("K-Means")
for i in range(0,numClusters):
ax3.plot(k_assignedData[i][:,0],k_assignedData[i][:,1],'.')
ax3.plot(k_centroids[i][0],k_centroids[i][1],'ko')
fig1.savefig('out.png', dpi=600)
The results in [1] suggest QK-Means presents a moderate improvement on the accuracy of the results relative to K-Means, while the same algorithm in [3] doesn't show significant improvement. The results above suggest that QK-Means presents a marginal accuracy improvement but a big handicap concerning computation time, taking significantly more time than classical K-Means to compute the same number of K-Means executions (each with a 300 iteration limit and $1\times10^{-4}$ convergence tolerance). It should be noted, however, that the algorithm in [1,3] uses the two quantum genetic operations not implemented here. Plus, their application was in image segmentation, while here we're using data from a mixture of Gaussians.
In conclusion, QK-Means doesn't offer a good tradeoff between accuracy improvement and computation time.
[1] E. Casper, C. Hung, E. Jung, and M. Yang, “A Quantum-Modeled K-Means Clustering Algorithm for Multi-band Image Segmentation”
[2] W. Liu, H. Chen, Q. Yan, Z. Liu, J. Xu, and Y. Zheng, “A novel quantum-inspired evolutionary algorithm based on variable angle-distance rotation”
[3] E. Casper and C. Hung, “Quantum Modeled Clustering Algorithms for Image Segmentation”