[PYTHON] Combinatorial optimization with quantum annealing

This is the 23rd day article of Nextremer Advent Calendar.

Currently, ** Nextremer Co., Ltd. ** of Waseda University ** [Mune Tanaka](http://www.slideshare.net/shu-t/ss- 57178026? Ref = http: //www.shutanaka.com/) ** and ** Joint research on quantum annealing * * We are doing.

Introduction

Quantum annealing is said to be an effective algorithm for solving combinatorial optimization problems.

In this article, we use quantum annealing by the quantum Monte Carlo method, which is a typical combinatorial optimization problem TSP (Traveling 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% I would like to find a solution for E9% A1% 8C).

Basics of quantum theory

Before we get into the TSP discussion, let me touch on the basics of quantum theory.

Physical quantities (energy, momentum, etc.) correspond to operators in the quantum world. For example, the operator corresponding to the physical quantity of energy is the Hamiltonian (expressed as $ \ hat {H} $).

This Hamiltonian $ \ hat {H} $ operator acts on a quantum state $ \ psi $.

\hat{H}\psi = E\psi

It has a structure in which the energy $ E $ is visible (observable) in the form of an eigenvalue.

In this way, operators can be said to be "spells" for knowing the quantum world. At the same time, operators are an indispensable tool for describing quantum theory.

Formulation of TSP

TSP is the problem of finding the route that minimizes the total distance when returning to the starting point via all cities once, given the set of n cities and the distance between each of the two cities.

Now suppose you have $ n $ cities $ 1,2, \ cdots, n $.

The distance between the cities $ i, j $ is $ d_ {ij} $, the variable $ that represents $ 1 $ if the city $ i $ is in the $ a $ step from the city $ 1 $, and $ 0 $ if it is not. Using n_ {i, a} $, the total distance $ L $ is

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

It can be expressed by. (The step at the departure point is the first step.)

However, due to the condition that one city can only stop once and the condition that it patrols,

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

Meet. If $ n_ {i, a} $ is represented graphically,

It is represented by a combination of 0 and 1.

** The purpose of this time is to find a pair of $ n_ {i, a} $ that minimizes this $ L $. ** **

Ising model

Electrons have a physical quantity called spin. As I mentioned earlier, since physical quantities correspond to operators in quantum theory, there are operators corresponding to spins, and [Pauli matrices](https://ja.wikipedia.org/wiki/ It is called% E3% 83% 91% E3% 82% A6% E3% 83% AA% E8% A1% 8C% E5% 88% 97).

There are three Pauli matrices, and you can use any of them, but in general, you can use $ \ hat {\ sigma} ^ {z} $ (an operator that expresses the spin angular momentum in the z direction). We get two eigenvalues $ 1, -1 $ as observations.

Since a binary variable can be expressed by using spin, it is used as a binary variable in combinatorial optimization.

Also, in order to actually solve combinatorial optimization, we have to prepare a lot of binary variables.

Therefore, a model in which spins are arranged on a grid point (** [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) **) is used.

The Hamiltonian of the system in this model is

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

Is described. ($ J_ {i, j} $ represents the strength of the interaction between spins, $ \ hat {\ sigma} _ {i} $ represents the variable of the spin at the grid point $ i $)

Therefore, the energy corresponding to this Hamiltonian is

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

Can be expressed as. ($ \ Sigma_ {i} ^ {z} $ is a spin variable, representing $ 1, -1 $)

If you look closely, it looks very similar to equation (1) above.

Therefore, if TSP is successfully converted to the Ising model, the optimum combination $ n_ {i, a} $ can be found in the TSP, and the set of spin variables that minimizes energy in the Ising model $ \ sigma_ {i, a You can see that finding $ is equivalent.

Convert TSP to Ising model

What we want to find in TSP is a set of $ n_ {i, a} $, which must take $ 0,1 $ and convert it to the spin variable $ 1, -1 $ in the Ising model.

So, $ n_ {i, a} $

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

Then, the $ \ (1 ) $ expression is

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

Can be expressed as. (However, the constant multiple of the first term ($ \ frac {1} {4} $) and the constant of the second term are terms that do not affect optimization, so they are ignored.)

So it could be rewritten by $ \ sigma_ {i, a} $ where L takes two values of $ 1, -1 $.

Also, the Hamiltonian for obtaining this L as a physical quantity is

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

It turns out that it is a Hamiltonian equivalent to the 2D Ising model.

Transverse magnetic field

I was able to set the problem by finding the set of $ \ sigma_ {i, a} $ that minimizes the energy obtained as the eigenvalue of (4).

Now, let's finally use the quantum annealing method.

If it remains (3), it can be obtained by the so-called simulated annealing method. Here, as a quantum effect, the transverse magnetic field term $ \ color {red} {-\ Gamma \ sum_ {i, a = 1} ^ {n} \ hat {\ sigma} _ {i, a} ^ {x}} $ To add.

\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} $ is the x component of Pauli matrices, $ \ Gamma $ is called the annealing coefficient, and will come out again later.

Partition function

Let's use (5) to find the partition function $ Z $ required for the quantum Monte Carlo method. The definition is

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

Given in. ($ \ Beta $ is the reciprocal of temperature, $ tr $ is the sum of the quantum states of all spins.)

Now, because there is a term for the transverse magnetic field, an off-diagonalization term appears, and it cannot be calculated straightforwardly.

Therefore, we approximate the partition function using a technique called Suzuki Trotter decomposition. With a sufficiently large $ m $,

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

will do. ($ \ Hat {H} _c, \ hat {H} _q $ represent the first and second terms of equation (5))

Along with that, the spin system|\sigma>A new dimension in the dimension of(Trotta dimension)With the addition of$n\times n\times m3D spin system|\sigma_{i,a,k}>$ (kはTrotta dimensionにおけるスピンの番号)Think about.

If you do a little long calculation from here, the partition function will be

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}

Can be obtained. (Please note that we cannot guarantee that this calculation result will be available.)

Quantum Monte Carlo simulation

If you look closely at equation (6), you can see that it is equivalent to the partition function of the three-dimensional Ising model, which has the following energy function.

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}

However, $ \ sigma_ {i, a, k} ^ {z} \ in \ (1, -1 ) $. A quantum annealing algorithm using the quantum Monte Carlo method is applied to this 3D Ising model.

The algorithm is

1. Randomly initialize 1, -1 so that each row and column has only one 1 and all others are -1 in each trotter dimension.
2. Randomly select the trotter dimension c and step numbers $ a, b (a \ neq b, a, b \ geq2) $, and let p and q be the column numbers where the element is 1 for each row of a and b, respectively. .. If the energy difference before and after flipping the four spins at the intersection of rows a and b and columns p and q is ∆E, the probability is $ min (1, e ^ {-\ beta \ Delta E}) $. Accept flip.
3. After repeating step 2 for Monte Carlo steps, reduce the annealing parameter $ \ Gamma $.
4. Repeat steps 2 and 3 to bring $ \ Gamma $ closer to 0.

Where $ \ Delta {E} $ is

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

Using,

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

Can be expressed as. According to this algorithm, the coordination of the Trotter dimension that minimizes the cost value by the original cost function (1) is the combination you want to find.

After that, if you change the -1,1 variable to the 0,1 variable, you will get the optimum combination you want, and the original purpose will be achieved.

Run

As a benchmark, we used the city data of "Djibouti" (Republic of Djibouti) on the site here.

The parameters used in the simulation are as follows.

--Initial reverse temperature $ \ beta $: 37 --Number of trottas $ m $: 10 --Initial value of annealing parameter $ \ Gamma_ {init} $: 1.0 --Quantum Monte Carlo Steps: 13320 --Attenuation factor of annealing parameter: 0.99

The (stupid) code looks like this.

QuantumAnnealing.py

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

#Parameter input
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

"""
Data about the city(City number and x,y coordinate)Get
"""
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(" "))

#City data
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 between two cities
def distance(point1, point2):
	return math.sqrt((point1[1]-point2[1])**2 + (point1[0]-point2[0])**2)


##Coordination of the entire spin
##Each dimension number of spin is[TROTTER_DIM, TOTAL_TIME, NCITY]Represented by
def getSpinConfig():

	##Coordination of spin at a certain time in a certain trotta dimension
	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 of spin in a 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


#Select the Trotter dimension that is the shortest distance and output the route at that time
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


##Total distance divided by max value for a 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


##Total distance to a 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

	
##Quantum Monte Carlo step
def QMC_move(config, ann_para):
	#Two different times a,choose 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_)

	#At a certain trotta number c, time a,City in b p,q
	p = config[c][a].index(1)
	q = config[c][b].index(1)

	#Initialize the energy difference value
	delta_cost = delta_costc = delta_costq_1 = delta_costq_2 = delta_costq_3 = delta_costq_4 = 0

	# (7)Difference in energy before and after flipping spin for the first term of the equation
	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)Energy difference before and after flipping spin for the second term of the equation
	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)About the formula Energy difference before and after flipping the spin
	delta_cost = delta_costc/TROTTER_DIM+para*(delta_costq_1+delta_costq_2+delta_costq_3+delta_costq_4)

	#Probability min(1,exp(-BETA*delta_cost))With 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


"""
Quantum annealing simulation
"""
if __name__ == '__main__':

	#Maximum distance between two cities
	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 of spin at initial time(Always in city 0 at initial time)
	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 "The shortest route is:{}".format(Route)
	print "The shortest distance is{}".format(Total_Length)
	print "The processing time is{}s".format(elapsed_time)
	
	plt.plot(LengthList)
	plt.show()

result

The vertical axis is the total distance and the horizontal axis is the annealing step.

The solution at convergence was 6737. Since the optimal solution in the benchmark is 6656, a relatively good approximate solution was obtained.

Finally

This time, using the TSP benchmark, we simulated quantum annealing using the quantum Monte Carlo method.

It is said in the streets that quantum annealing can solve combinatorial optimization problems, but I created it because I didn't have much specific information. I hope it helps you.