Ceci est l'article du 23ème jour du Calendrier de l'Avent Nextremer.

Actuellement, ** Nextremer Co., Ltd. ** de l'Université Waseda ** [Mune Tanaka](http://www.slideshare.net/shu-t/ss- 57178026? Ref = http: //www.shutanaka.com/) ** et ** Recherche conjointe sur le recuit quantique * * Nous faisons.

introduction

On dit que le recuit quantique est un algorithme efficace pour résoudre les problèmes d'optimisation de combinaison.

Dans cet article, nous utilisons le recuit quantique par la méthode quantique de Monte Carlo, qui est un problème d'optimisation de combinaison typique 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% Je voudrais trouver une solution pour E9% A1% 8C).

Bases de la théorie quantique

Avant d'entrer dans la discussion TSP, laissez-moi aborder les bases de la théorie quantique.

Les grandeurs physiques (comme l'énergie et la quantité de mouvement) correspondent aux opérateurs du monde quantique. Par exemple, la quantité physique d'énergie correspondante est l'opérateur hamiltonien (représenté par $ \ hat {H} $).

Cet opérateur hamiltonien $ \ hat {H} $ agit sur un état quantique $ \ psi $.

\hat{H}\psi = E\psi

Il a une structure dans laquelle l'énergie $ E $ est visible (observable) sous la forme d'une valeur propre.

De cette manière, l'opérateur peut être considéré comme un "sortilège" pour connaître le monde quantique. Dans le même temps, les opérateurs sont un outil indispensable pour décrire la théorie quantique.

Formulation de TSP

TSP est le problème de trouver l'itinéraire qui minimise la distance totale lors du retour au point de départ via toutes les villes une fois, étant donné l'ensemble de n villes et la distance entre chacune des deux villes.

Supposons maintenant que vous ayez $ n $ villes $ 1,2, \ cdots, n $.

La distance entre les villes $ i et j $ est $ d_ {ij} $, la variable $ qui représente $ 1 $ si la ville $ i $ est dans le pas $ a $ de la ville $ 1 $, et $ 0 $ sinon. En utilisant n_ {i, a} $, la distance totale $ L $ est

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

Il peut être exprimé par. (L'étape au point de départ est la première étape.)

Cependant, en raison de la condition qu'une ville ne peut s'arrêter qu'une seule fois et de la condition qu'elle patrouille,

\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)

Rencontrer. Si $ n_ {i, a} $ est représenté graphiquement,

Il est représenté par une combinaison de 0 et 1.

** Le but de ce temps est de trouver une paire de $ n_ {i, a} $ qui minimise ce $ L $. ** **

Modèle en hausse

Les électrons ont une quantité physique appelée spin. Comme je l'ai mentionné précédemment, puisque les quantités physiques correspondent aux opérateurs de la théorie quantique, il existe des opérateurs correspondant aux spins, et [Pauli matrix](https://ja.wikipedia.org/wiki/ Il est appelé% E3% 83% 91% E3% 82% A6% E3% 83% AA% E8% A1% 8C% E5% 88% 97).

Il existe trois matrices de Pauli, et vous pouvez utiliser n'importe laquelle d'entre elles, mais en général, vous pouvez utiliser $ \ hat {\ sigma} ^ {z} $ (un opérateur qui exprime la dynamique de l'angle de rotation dans la direction z). Nous obtenons deux valeurs propres de $ 1, -1 $ comme observations.

Puisqu'une variable binaire peut être exprimée à l'aide de spin, elle est utilisée comme variable binaire dans l'optimisation de combinaison.

De plus, pour résoudre réellement l'optimisation des combinaisons, nous devons préparer un grand nombre de variables binaires.

Par conséquent, un modèle dans lequel les spins sont disposés sur des points de grille (** [Ising model](https://ja.wikipedia.org/wiki/%E3%82%A4%E3%82%B8%E3%83%B3%] E3% 82% B0% E6% A8% A1% E5% 9E% 8B) **) est utilisé.

L'hamiltonien du système dans ce modèle est

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

Est décrit. ($ J_ {i, j} $ représente la force de l'interaction entre les spins, $ \ hat {\ sigma} _ {i} $ représente la variable de spin au point de grille $ i $)

Par conséquent, l'énergie correspondant à cet hamiltonien est

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

Peut être exprimé comme. ($ \ Sigma_ {i} ^ {z} $ est une variable de spin, représentant $ 1, -1 $)

Si vous regardez de plus près, cela ressemble beaucoup à l'équation (1) ci-dessus.

Par conséquent, si TSP est converti avec succès dans le modèle ascendant, la combinaison optimale $ n_ {i, a} $ peut être trouvée dans TSP, et l'ensemble de variables de spin qui minimise l'énergie dans le modèle ascendant $ \ sigma_ {i, a Vous pouvez voir que trouver $ est équivalent.

Convertir TSP en modèle Rising

Ce que nous voulons trouver dans TSP est un ensemble de $ n_ {i, a} $, qui doit prendre $ 0,1 $ et le convertir en la variable spin $ 1, -1 $ dans le modèle d'Ising.

Donc, $ n_ {i, a} $

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

Ensuite, l'expression $ \ (1 ) $ est

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

Peut être exprimé comme. (Cependant, le multiple constant du premier terme ($ \ frac {1} {4} $) et la constante du deuxième terme sont des termes qui n'affectent pas l'optimisation, ils sont donc ignorés.)

Il pourrait donc être réécrit par $ \ sigma_ {i, a} $ où L prend deux valeurs de $ 1, -1 $.

De plus, l'hamiltonien pour obtenir ce L en tant que quantité physique est

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

On voit que c'est un hamiltonien équivalent au modèle 2D Rising.

Champ magnétique transversal

J'ai pu poser le problème en trouvant l'ensemble de $ \ sigma_ {i, a} $ qui minimise l'énergie obtenue comme valeur propre de (4).

Maintenant, utilisons enfin la méthode de recuit quantique.

S'il reste (3), il peut être obtenu par la méthode dite de recuit simulé. Ici, comme effet quantique, le terme de champ magnétique transverse $ \ color {red} {- \ Gamma \ sum_ {i, a = 1} ^ {n} \ hat {\ sigma} _ {i, a} ^ {x}} $ Ajouter.

\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} $ est le composant x de la matrice de Pauli, $ \ Gamma $ est appelé le coefficient de recuit, et reviendra plus tard.

Fonction de distribution

Utilisons (5) pour trouver la fonction de distribution $ Z $ requise pour la méthode quantique de Monte Carlo. La définition est

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

Donné dans. ($ \ Beta $ est l'inverse de la température, $ tr $ est la somme des états quantiques de tous les spins.)

Maintenant, parce qu'il existe un terme pour le champ magnétique transversal, un terme hors diagonale apparaît, et il ne peut pas être calculé directement.

Par conséquent, nous allons approximer la fonction de distribution en utilisant une technique appelée décomposition Suzuki Trotter. Avec un assez gros $ m $,

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

ça ira. ($ \ Hat {H} _c, \ hat {H} _q $ représentent les premier et deuxième termes de l'équation (5))

Parallèlement à cela, le système de rotation|\sigma>Une nouvelle dimension dans la dimension de(Dimension trotta)Avec l'ajout den\times n\times mSystème de rotation 3D|\sigma_{i,a,k}> (kはDimension trottaにおけるスピンの番号)Penser à.

Si vous faites un petit calcul long à partir d'ici, la fonction de distribution sera

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}

Peut être obtenu. (Veuillez noter que nous ne pouvons garantir que ce résultat de calcul sera disponible.)

Simulation quantique de Monte Carlo

Si vous regardez de près l'équation (6), vous pouvez voir qu'elle équivaut à la fonction de distribution du modèle Ising tridimensionnel, qui a la fonction d'énergie suivante.

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}

Cependant, $ \ sigma_ {i, a, k} ^ {z} \ in \ (1, -1 ) $. Un algorithme de recuit quantique utilisant la méthode quantique de Monte Carlo est appliqué à ce modèle 3D ing.

L'algorithme est

Initialisez aléatoirement 1, -1 pour que chaque ligne et colonne n'ait qu'un 1 et que toutes les autres soient -1 dans chaque dimension de trotteur.
Sélectionnez au hasard la dimension du trotteur c et les numéros de pas $ a, b (a \ neq b, a, b \ geq2) $, et soit p et q les numéros de colonne où l'élément est 1 pour chaque ligne de a et b, respectivement. .. Si la différence d'énergie avant et après le retournement des quatre spins à l'intersection des lignes a et b et des colonnes p et q est ∆E, la probabilité est $ min (1, e ^ {- \ beta \ Delta E}) $. Acceptez le flip.
Après avoir répété l'étape 2 pour les étapes Monte Carlo, réduisez le paramètre de recuit $ \ Gamma $.
Répétez les étapes 2 et 3 pour rapprocher $ \ Gamma $ de 0.

Où $ \ Delta {E} $ est

\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})\}

En utilisant,

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

Peut être exprimé comme. Selon cet algorithme, la coordination de la dimension Trotter qui minimise la valeur de coût par la fonction de coût d'origine (1) est la combinaison que vous voulez trouver.

Après cela, si vous remplacez la variable -1,1 par la variable 0,1, vous obtiendrez la combinaison optimale que vous vouliez à l'origine et le but initial sera atteint.

Courir

Comme référence, nous avons utilisé les données de la ville de "Djibouti" (République de Dijibouti) sur le site ici.

Les paramètres utilisés pour la simulation sont les suivants.

--Température inverse initiale $ \ beta $: 37 --Nombre de trottas $ m $: 10 --Valeur initiale du paramètre de recuit $ \ Gamma_ {init} $: 1.0 --Quantum Monte Carlo étapes: 13320

Taux de décroissance du paramètre d'annelage: 0,99

Le code (stupide) ressemble à ceci.

`QuantumAnnealing.py`


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

#Paramètres d'entrée
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

"""
Données sur la ville(Numéro de ville et x,coordonnée y)Avoir
"""
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(" "))

#Données de la ville
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])

##Distance entre deux villes
def distance(point1, point2):
	return math.sqrt((point1[1]-point2[1])**2 + (point1[0]-point2[0])**2)


##Coordination de l'ensemble du spin
##Chaque numéro de dimension de spin[TROTTER_DIM, TOTAL_TIME, NCITY]Représenté par
def getSpinConfig():

	##Coordination de spin à un certain moment dans une certaine dimension trotta
	def spin_config_at_a_time_in_a_TROTTER_DIM(tag):
		config = list(-np.ones(NCITY, dtype = np.int))
		config[tag] = 1
		return config


	##Coordination du spin dans une dimension trotta
	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


#Sélectionnez la dimension Trotter qui est la distance la plus courte et sortez l'itinéraire à ce moment
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


##Distance totale divisée par la valeur maximale d'un itinéraire
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


##Distance totale jusqu'à un itinéraire
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

	
##Étape Quantum Monte Carlo
def QMC_move(config, ann_para):
	#Deux fois un,choisissez 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_)

	#À un certain nombre trotta c, temps a,Ville en b p,q
	p = config[c][a].index(1)
	q = config[c][b].index(1)

	#Initialiser la valeur de la différence d'énergie
	delta_cost = delta_costc = delta_costq_1 = delta_costq_2 = delta_costq_3 = delta_costq_4 = 0

	# (7)Différence d'énergie avant et après le retournement du spin pour le premier terme de l'équation
	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)Différence d'énergie avant et après le retournement du spin pour le deuxième terme de l'équation
	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)À propos de la formule Différence d'énergie avant et après avoir inversé la rotation
	delta_cost = delta_costc/TROTTER_DIM+para*(delta_costq_1+delta_costq_2+delta_costq_3+delta_costq_4)

	#Probabilité min(1,exp(-BETA*delta_cost))Avec 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


"""
Simulation de recuit quantique
"""
if __name__ == '__main__':

	#Distance maximale entre deux villes
	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])


	#Coordination de la rotation au moment initial(Toujours dans la ville 0 à l'heure initiale)
	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 "Le trajet le plus court est:{}".format(Route)
	print "La distance la plus courte est{}".format(Total_Length)
	print "Le temps de traitement est{}s".format(elapsed_time)
	
	plt.plot(LengthList)
	plt.show()

résultat

L'axe vertical est la distance totale et l'axe horizontal est l'étape de recuit.

La solution au moment de la convergence était 6737. Puisque la solution optimale dans la référence est 6656, une solution approximative relativement bonne a été obtenue.

finalement

Cette fois, en utilisant le benchmark TSP, nous avons simulé le recuit quantique en utilisant la méthode quantique de Monte Carlo.

On dit dans la rue que le problème d'optimisation des combinaisons peut être résolu par recuit quantique, mais je l'ai créé parce qu'il n'y avait pas beaucoup d'informations spécifiques. J'espère que ça t'aide.

[PYTHON] Optimisation combinée avec recuit quantique