Dies ist der Artikel zum 23. Tag von Nextremer Adventskalender.

Derzeit ** Nextremer Co., Ltd. ** der Waseda University ** [Mune Tanaka](http://www.slideshare.net/shu-t/ss- 57178026? Ref = http: //www.shutanaka.com/) ** und ** Gemeinsame Forschung zum Quantenglühen * * Wir machen.

Einführung

Das Quantenglühen soll ein effektiver Algorithmus zur Lösung von Kombinationsoptimierungsproblemen sein.

In diesem Artikel verwenden wir das Quantenglühen nach der Quanten-Monte-Carlo-Methode, einem typischen Problem der Kombinationsoptimierung TSP (Circuit Salesman Problem). A1% E5% 9B% 9E% E3% 82% BB% E3% 83% BC% E3% 83% AB% E3% 82% B9% E3% 83% 9E% E3% 83% B3% E5% 95% 8F% Ich möchte eine Lösung für E9% A1% 8C finden.

Grundlagen der Quantentheorie

Bevor wir zur TSP-Diskussion kommen, möchte ich auf die Grundlagen der Quantentheorie eingehen.

Physikalische Größen (wie Energie und Impuls) entsprechen Operatoren in der Quantenwelt. Beispielsweise ist die entsprechende physikalische Energiemenge der Hamilton-Operator (dargestellt als $ \ hat {H} $).

Dieser Hamilton-Operator $ \ hat {H} $ wirkt auf einen Quantenzustand $ \ psi $.

\hat{H}\psi = E\psi

Es hat eine Struktur, in der die Energie $ E $ in Form eines Eigenwerts sichtbar (beobachtbar) ist.

Auf diese Weise kann der Operator als "Zauber" für die Kenntnis der Quantenwelt bezeichnet werden. Gleichzeitig sind Operatoren ein unverzichtbares Werkzeug zur Beschreibung der Quantentheorie.

Formulierung von TSP

TSP ist das Problem, die Route zu finden, die die Gesamtentfernung minimiert, wenn über alle Städte einmal zum Startpunkt zurückgekehrt wird, angesichts der Menge von n Städten und der Entfernung zwischen jeder der beiden Städte.

Angenommen, Sie haben $ n $ Städte $ 1,2, \ cdots, n $.

Der Abstand zwischen den Städten $ i und j $ beträgt $ d_ {ij} $, die Variable $, die $ 1 $ darstellt, wenn sich die Stadt $ i $ im Schritt $ a $ von der Stadt $ 1 $ befindet, und $ 0 $, wenn nicht. Mit n_ {i, a} $ beträgt der Gesamtabstand $ L $

L = \sum_{i,j,a=1}^{n}d_{i,j} n_{i,a} n_{j,a+1} \tag{1}

Es kann ausgedrückt werden durch. (Der Schritt am Abfahrtspunkt ist der erste Schritt.)

Aufgrund der Bedingung, dass eine Stadt nur einmal anhalten kann, und der Bedingung, dass sie patrouilliert,

\sum_{i=1}^{n}n_{i,a}=\sum_{a=1}^{n}n_{i,a}=1 \\
n_{i,n+1}=n_{i,1}  \ \ \ \ \ \ (１\leq i \leq n)

Treffen. Wenn $ n_ {i, a} $ grafisch dargestellt wird,

Es wird durch eine Kombination von 0 und 1 dargestellt.

** Der Zweck dieser Zeit ist es, ein Paar von $ n_ {i, a} $ zu finden, das dieses $ L $ minimiert. ** **.

Steigendes Modell

Elektronen haben eine physikalische Größe, die als Spin bezeichnet wird. Wie bereits erwähnt, gibt es Operatoren, die Spins und [Pauli-Matrix] entsprechen, da physikalische Größen Operatoren in der Quantentheorie entsprechen (https://ja.wikipedia.org/wiki/). Es heißt% E3% 83% 91% E3% 82% A6% E3% 83% AA% E8% A1% 8C% E5% 88% 97).

Es gibt drei Pauli-Matrizen, und Sie können jede von ihnen verwenden, aber im Allgemeinen können Sie $ \ hat {\ sigma} ^ {z} $ verwenden (ein Operator, der den Drehwinkelimpuls in z-Richtung ausdrückt). Wir erhalten zwei Eigenwerte von $ 1, -1 $ als Beobachtungen.

Da eine binäre Variable mit Spin ausgedrückt werden kann, wird sie bei der Kombinationsoptimierung als binäre Variable verwendet.

Um die Kombinationsoptimierung tatsächlich zu lösen, müssen wir viele binäre Variablen vorbereiten.

Daher ein Modell, in dem Drehungen auf Gitterpunkten angeordnet sind (** [Ising-Modell](https://ja.wikipedia.org/wiki/%E3%82%A4%E3%82%B8%E3%83%B3%] E3% 82% B0% E6% A8% A1% E5% 9E% 8B) **) wird verwendet.

Der Hamilton-Operator des Systems in diesem Modell ist

\hat{H}=-\sum_{i,j}J_{i,j}\hat{\sigma}^{z}_{i}\hat{\sigma}^{z}_{j} \tag{2}

Wird beschrieben. ($ J_ {i, j} $ repräsentiert die Stärke der Interaktion zwischen Spins, $ \ hat {\ sigma} _ {i} $ repräsentiert die Spinvariable am Gitterpunkt $ i $)

Daher ist die diesem Hamilton-Operator entsprechende Energie

H = -\sum_{i,j}J_{i,j}\sigma^{z}_{i}\sigma^{z}_{j}

Kann ausgedrückt werden als. ($ \ Sigma_ {i} ^ {z} $ ist eine Spinvariable, die $ 1, -1 $ darstellt.)

Wenn Sie genau hinschauen, sieht es der obigen Gleichung (1) sehr ähnlich.

Wenn TSP erfolgreich in das Ising-Modell konvertiert wurde, kann daher die optimale Kombination $ n_ {i, a} $ in TSP und der Menge von Spinvariablen gefunden werden, die die Energie im Ising-Modell $ \ sigma_ {i, a minimieren. Sie können sehen, dass das Finden von $ gleichwertig ist.

Konvertieren Sie TSP in Rising-Modell

Was wir in TSP finden wollen, ist eine Menge von $ n_ {i, a} $, die $ 0,1 $ nehmen und im Ising-Modell in die Spinvariable $ 1, -1 $ konvertieren muss.

Also $ n_ {i, a} $

n_{i,a}\equiv\frac{\sigma_{i,a}+1}{2}

Dann ist der Ausdruck $ \ (1 ) $

L=\frac{1}{4}\sum_{a,i,j=1}^{n} d_{i,j} \sigma_{i,a}\sigma_{j,a+1}+(const.) \tag{3}

Kann ausgedrückt werden als. (Das konstante Vielfache des ersten Terms ($ \ frac {1} {4} $) und die Konstante des zweiten Terms sind jedoch Terme, die die Optimierung nicht beeinflussen, sodass sie ignoriert werden.)

Es könnte also durch $ \ sigma_ {i, a} $ umgeschrieben werden, wobei L zwei Werte von $ 1, -1 $ annimmt.

Auch der Hamilton-Operator zum Erhalten dieses L als physikalische Größe ist

\hat{H}=\sum_{a,i,j=1}^{n} d_{i,j} \hat{\sigma}^{z}_{i,a}\hat{\sigma}^{z}_{j,a+1}　\tag{4}

Es ist ersichtlich, dass es sich um ein Hamilton-Äquivalent zum 2D-Rising-Modell handelt.

Transversales Magnetfeld

Die Problemeinstellung war, dass das Paar von $ \ sigma_ {i, a} $ gefunden werden sollte, das die als Eigenwert von (4) erhaltene Energie minimiert.

Verwenden wir nun endlich die Quantenglühmethode.

Bleibt es (3), kann es durch das sogenannte simulierte Glühverfahren erhalten werden. Hier wird als Quanteneffekt der transversale Magnetfeldterm $ \ color {red} {- \ Gamma \ sum_ {i, a = 1} ^ {n} \ hat {\ sigma} _ {i, a} ^ {x}} $ Hinzufügen.

\hat{H}=\sum_{a,i,j=1}^{n} d_{i,j} \hat{\sigma}^{z}_{i,a}\hat{\sigma}^{z}_{j,a+1}\color{red}{-\Gamma\sum_{i,a=1}^{n}\hat{\sigma}_{i,a}^{x}} \tag{5}

$ \ hat {\ sigma} _ {i, a} ^ {x} $ ist die x-Komponente der Pauli-Matrix, $ \ Gamma $ wird als Glühkoeffizient bezeichnet und wird später wieder herauskommen.

Verteilungsfunktion

Verwenden wir (5), um die Verteilungsfunktion $ Z $ zu finden, die für die Quanten-Monte-Carlo-Methode erforderlich ist. Die Definition ist

Z\equiv tr e^{-\beta\hat{H}}=\sum_{\sigma}<\sigma|e^{-\beta\hat{H}}|\sigma>

Gegeben in. ($ \ Beta $ ist die Umkehrung der Temperatur, $ tr $ ist die Summe für die Quantenzustände aller Spins.)

Da es nun einen Term für das transversale Magnetfeld gibt, erscheint ein nicht diagonaler Term, der nicht einfach berechnet werden kann.

Daher werden wir die Verteilungsfunktion unter Verwendung einer Technik approximieren, die als Suzuki Trotter-Zerlegung bezeichnet wird. Mit einem ausreichend großen $ m $,

Z=tre^{-\beta\hat{H}}\simeq tr(e^{-\frac{\beta}{m}\hat{H}_{c}}e^{-\frac{\beta}{m}\hat{H}_{q}})^m

Wird besorgt. ($ \ Hat {H} _c, \ hat {H} _q $ repräsentieren den ersten und zweiten Term von Gleichung (5))

Dazu das Spin-System|\sigma>Eine neue Dimension in der Dimension von(Trotta Dimension)Mit dem Zusatz von$n\times n\times m3D-Spin-System|\sigma_{i,a,k}>$ (kはTrotta Dimensionにおけるスピンの番号)Nachdenken über.

Wenn Sie von hier aus eine etwas lange Berechnung durchführen, ist die Verteilungsfunktion

Z \simeq \biggl[\frac{1}{2}sinh\biggl(\frac{2\beta\Gamma}{m}\biggr)\biggr]^{\frac{mn^2}{2}}\sum_{\sigma_{i,a,k}^z=\pm 1}exp\biggl[-\frac{\beta}{m} \sum_{a=1}^{n}\sum_{i,j=1}^{n}\sum_{k=1}^{m}d_{i,j}\sigma_{i,a,k}^{z}\sigma_{j,a+1,k}^{z} +\frac{1}{2}\log{coth\biggl(\frac{\beta \Gamma}{m}\biggr)}\sum_{i,a=1}^{n}\sum_{k=1}^{m}\sigma_{i,a,k}^{z}\sigma_{i,a,k+1}^{z}\biggr] \tag{6}

Kann erhalten werden. (Bitte beachten Sie, dass wir nicht garantieren können, dass dieses Berechnungsergebnis verfügbar ist.)

Quanten-Monte-Carlo-Simulation

Wenn Sie sich Gleichung (6) genau ansehen, können Sie sehen, dass sie der Verteilungsfunktion des dreidimensionalen Ising-Modells entspricht, das die folgende Energiefunktion hat.

E = \frac{1}{m} \sum_{a=1}^{n}\sum_{i,j=1}^{n}\sum_{k=1}^{m}d_{i,j}\sigma_{i,a,k}^{z}\sigma_{j,a+1,k}^{z} -\frac{1}{2\beta}\log{coth\biggl(\frac{\beta \Gamma}{m}\biggr)}\sum_{i,a=1}^{n}\sum_{k=1}^{m}\sigma_{i,a,k}^{z}\sigma_{i,a,k+1}^{z} \tag{7}

$ \ Sigma_ {i, a, k} ^ {z} \ in \ (1, -1 ) $. Ein Quantenglühalgorithmus unter Verwendung der Quanten-Monte-Carlo-Methode wird auf dieses 3D-Modell angewendet.

Der Algorithmus ist

Initialisieren Sie in jeder Traber-Dimension zufällig 1, -1, sodass jede Zeile und jede Spalte nur eine 1 hat und alle anderen -1 sind.
Wählen Sie zufällig die Traberdimension c und die Schrittnummern $ a, b (a \ neq b, a, b \ geq2) $ aus und lassen Sie p und q die Spaltennummern sein, wobei das Element für jede Zeile von a bzw. b 1 ist. .. Wenn die Energiedifferenz vor und nach dem Umdrehen der vier Drehungen am Schnittpunkt der Zeilen a und b und der Spalten p und q ∆E beträgt, beträgt die Wahrscheinlichkeit $ min (1, e ^ {- \ beta \ Delta E}) $. Akzeptiere Flip.
Reduzieren Sie nach dem Wiederholen von Schritt 2 für Monte-Carlo-Schritte den Glühparameter $ \ Gamma $.
Wiederholen Sie die Schritte 2 und 3, um $ \ Gamma $ näher an 0 zu bringen.

Wo $ \ Delta {E} $ ist

\Delta E_{1}=\frac{2}{m}\sum_{i=1}^{n}(\sigma_{i,a-1,c}^{z}+\sigma_{i,a+1,c}^{z}-\sigma_{i,b-1,c}^{z}-\sigma_{i,b+1,c}^{z})(d_{q,i}-d_{p,i}) \\
\Delta E_{2}=\frac{1}{\beta}\log{coth\biggl(\frac{\beta \Gamma}{m}\biggr)}\{\sigma_{p,a,c}^{z}(\sigma_{p,a,c-1}^{z}+\sigma_{p,a,c+1}^{z})+\sigma_{q,a,c}^{z}(\sigma_{q,a,c-1}^{z}+\sigma_{q,a,c+1}^{z})\} \\
\Delta E_{3}=\frac{1}{\beta}\log{coth\biggl(\frac{\beta \Gamma}{m}\biggr)}\{\sigma_{p,b,c}^{z}(\sigma_{p,b,c-1}^{z}+\sigma_{p,b,c+1}^{z})+\sigma_{q,b,c}^{z}(\sigma_{q,b,c-1}^{z}+\sigma_{q,b,c+1}^{z})\}

Verwenden von,

\Delta E=\Delta E_{1}+\Delta E_{2}+\Delta E_{3}

Kann ausgedrückt werden als. Nach diesem Algorithmus ist die Koordination der Trotta-Dimension, die den Kostenwert durch die ursprüngliche Kostenfunktion (1) minimiert, die Kombination, die Sie finden möchten.

Wenn Sie danach die Variable -1,1 in die Variable 0,1 ändern, erhalten Sie die optimale Kombination, die Sie ursprünglich wollten, und der ursprüngliche Zweck wird erreicht.

Lauf

Als Benchmark haben wir die Stadtdaten von "Dschibuti" (Republik Dischibuti) auf der Website hier verwendet.

Die für die Simulation verwendeten Parameter sind wie folgt.

Anfängliche Umkehrtemperatur $ \ beta $: 37
Anzahl der Trottas $ m $: 10
Anfangswert des Glühparameters $ \ Gamma_ {init} $: 1.0
Quantum Monte Carlo Schritte: 13320
Abklingrate der Glühparameter: 0,99

Der (dumme) Code sieht so aus.

`QuantumAnnealing.py`


#coding:utf-8
import time
import math
import numpy as np
import os
import random
import matplotlib.pyplot as plt

#Eingabeparameter
TROTTER_DIM = int(input("Trotter dimension: "))
ANN_PARA =  float(input("initial annealing parameter: "))
ANN_STEP = int(input("Annealing Step: "))
MC_STEP = int(input("MC step: "))
BETA = float(input("inverse Temperature: "))
REDUC_PARA = 0.99

"""
Daten über die Stadt(Stadtnummer und x,y-Koordinate)Erhalten
"""
FILE_NAME = 'FILE_NAME'

f = open(os.path.dirname(os.path.abspath(FILE_NAME))+FILE_NAME).read().split("\n")

POINT = []
for i in f:
	POINT.append(i.split(" "))

#Stadtdaten
NCITY = len(POINT)
TOTAL_TIME = NCITY
for i in range(NCITY):
	POINT[i].remove(POINT[i][0])
for i in range(NCITY):
	for j in range(2):
		POINT[i][j] = float(POINT[i][j])

##Entfernung zwischen zwei Städten
def distance(point1, point2):
	return math.sqrt((point1[1]-point2[1])**2 + (point1[0]-point2[0])**2)


##Koordination des gesamten Spins
##Jede Dimension Nummer des Spins[TROTTER_DIM, TOTAL_TIME, NCITY]Vertreten durch
def getSpinConfig():

	##Koordination des Spins zu einem bestimmten Zeitpunkt in einer bestimmten Trottadimension
	def spin_config_at_a_time_in_a_TROTTER_DIM(tag):
		config = list(-np.ones(NCITY, dtype = np.int))
		config[tag] = 1
		return config


	##Koordination des Spins in einer Trotta-Dimension
	def spin_config_in_a_TROTTER_DIM(tag):
		spin = []
		spin.append(config_at_init_time)
		for i in xrange(TOTAL_TIME-1):
			spin.append(list(spin_config_at_a_time_in_a_TROTTER_DIM(tag[i])))
		return spin

	spin = []
	for i in xrange(TROTTER_DIM):
		tag = np.arange(1,NCITY)
		np.random.shuffle(tag)
		spin.append(spin_config_in_a_TROTTER_DIM(tag)) 
	return spin


#Wählen Sie die Trotter-Dimension aus, die die kürzeste Entfernung darstellt, und geben Sie die Route zu diesem Zeitpunkt aus
def getBestRoute(config):	
	length = []
	for i in xrange(TROTTER_DIM):
		route = []
		for j in xrange(TOTAL_TIME):
			route.append(config[i][j].index(1))
		length.append(getTotaldistance(route))

	min_Tro_dim = np.argmin(length)
	Best_Route = []
	for i in xrange(TOTAL_TIME):
		Best_Route.append(config[min_Tro_dim][i].index(1))
	return Best_Route


##Gesamtentfernung geteilt durch den Maximalwert für eine Route
def getTotaldistance(route):
	Total_distance = 0
	for i in xrange(TOTAL_TIME):
		Total_distance += distance(POINT[route[i]],POINT[route[(i+1)%TOTAL_TIME]])/max_distance
	return Total_distance


##Gesamtentfernung zu einer Route
def getRealTotaldistance(route):
	Total_distance = 0
	for i in xrange(TOTAL_TIME):
		Total_distance += distance(POINT[route[i]], POINT[route[(i+1)%TOTAL_TIME]])
	return Total_distance

	
##Quanten-Monte-Carlo-Schritt
def QMC_move(config, ann_para):
	#Zwei verschiedene Zeiten a,wähle b
	c = np.random.randint(0,TROTTER_DIM)
	a_ = range(1,TOTAL_TIME)
	a = np.random.choice(a_)
	a_.remove(a)
	b = np.random.choice(a_)

	#Bei einer bestimmten Trottazahl c, Zeit a,Stadt in b p,q
	p = config[c][a].index(1)
	q = config[c][b].index(1)

	#Initialisieren Sie den Wert der Energiedifferenz
	delta_cost = delta_costc = delta_costq_1 = delta_costq_2 = delta_costq_3 = delta_costq_4 = 0

	# (7)Energiedifferenz vor und nach dem Umdrehen des Spins für den ersten Term der Gleichung
	for j in range(NCITY):
		l_p_j = distance(POINT[p], POINT[j])/max_distance
		l_q_j = distance(POINT[q], POINT[j])/max_distance
		delta_costc += 2*(-l_p_j*config[c][a][p] - l_q_j*config[c][a][q])*(config[c][a-1][j]+config[c][(a+1)%TOTAL_TIME][j])+2*(-l_p_j*config[c][b][p] - l_q_j*config[c][b][q])*(config[c][b-1][j]+config[c][(b+1)%TOTAL_TIME][j])

	# (7)Energiedifferenz vor und nach dem Umdrehen des Spins für den zweiten Term der Gleichung
	para = (1/BETA)*math.log(math.cosh(BETA*ann_para/TROTTER_DIM)/math.sinh(BETA*ann_para/TROTTER_DIM))
	delta_costq_1 = config[c][a][p]*(config[(c-1)%TROTTER_DIM][a][p]+config[(c+1)%TROTTER_DIM][a][p])
	delta_costq_2 = config[c][a][q]*(config[(c-1)%TROTTER_DIM][a][q]+config[(c+1)%TROTTER_DIM][a][q])
	delta_costq_3 = config[c][b][p]*(config[(c-1)%TROTTER_DIM][b][p]+config[(c+1)%TROTTER_DIM][b][p])
	delta_costq_4 = config[c][b][q]*(config[(c-1)%TROTTER_DIM][b][q]+config[(c+1)%TROTTER_DIM][b][q])

    # (7)Über die Formel Energiedifferenz vor und nach dem Umdrehen des Spins
	delta_cost = delta_costc/TROTTER_DIM+para*(delta_costq_1+delta_costq_2+delta_costq_3+delta_costq_4)

	#Wahrscheinlichkeit min(1,exp(-BETA*delta_cost))Mit Flip
	if delta_cost <= 0:
		config[c][a][p]*=-1
		config[c][a][q]*=-1
		config[c][b][p]*=-1
		config[c][b][q]*=-1
	elif np.random.random() < np.exp(-BETA*delta_cost):
		config[c][a][p]*=-1
		config[c][a][q]*=-1
		config[c][b][p]*=-1
 		config[c][b][q]*=-1

	return config


"""
Quantenglühensimulation
"""
if __name__ == '__main__':

	#Maximale Entfernung zwischen zwei Städten
	max_distance = 0
	for i in range(NCITY):
		for j in range(NCITY):
			if max_distance <= distance(POINT[i], POINT[j]):
				max_distance = distance(POINT[i], POINT[j])


	#Koordination des Spins zum Anfangszeitpunkt(Zum ersten Mal immer in Stadt 0)
	config_at_init_time = list(-np.ones(NCITY,dtype=np.int))
	config_at_init_time[0] = 1

	print "start..."
	t0 = time.clock()

	np.random.seed(100)
	spin = getSpinConfig()
	LengthList = []
	for t in range(ANN_STEP):
		for i in range(MC_STEP):
			con = QMC_move(spin, ANN_PARA)
			rou = getBestRoute(con)
			length = getRealTotaldistance(rou)
		LengthList.append(length)
		print "Step:{}, Annealing Parameter:{}, length:{}".format(t+1,ANN_PARA, length)
		ANN_PARA *= REDUC_PARA

	Route = getBestRoute(spin)
	Total_Length = getRealTotaldistance(Route)
	elapsed_time = time.clock()-t0

	print "Der kürzeste Weg ist:{}".format(Route)
	print "Die kürzeste Entfernung ist{}".format(Total_Length)
	print "Die Bearbeitungszeit beträgt{}s".format(elapsed_time)
	
	plt.plot(LengthList)
	plt.show()

Ergebnis

Die vertikale Achse ist der Gesamtabstand und die horizontale Achse ist der Glühschritt.

Die Lösung zum Zeitpunkt der Konvergenz war 6737. Da die optimale Lösung im Benchmark 6656 ist, wurde eine relativ gute ungefähre Lösung erhalten.

Schließlich

Dieses Mal haben wir mithilfe des TSP-Benchmarks das Quantenglühen mit der Quanten-Monte-Carlo-Methode simuliert.

Auf den Straßen wird gesagt, dass das Problem der Kombinationsoptimierung durch Quantenglühen gelöst werden kann, aber ich habe es erstellt, weil es nicht viele spezifische Informationen gab. Ich hoffe es hilft dir.

[PYTHON] Kombinationsoptimierung mit Quantenglühen