The story of the algorithm drawing a ridiculous conclusion when trying to solve the traveling salesman problem properly

This article is the 13th day article of Furukawa Lab Advent_calendar. This article was written by a student at Furukawa Lab as part of his studies. The content may be ambiguous or the expression may be slightly different.

Introduction

Hello everyone Do you know the traveling salesman problem? It's a kind of common combination optimization problem, but in this article I'd like to write about a sad case when I tried to solve the traveling salesman problem properly. To conclude first, I'm afraid of rounding errors. Since it was a story when I was a beginner in the program, I hope you can read it with a low hurdle in everything.

Traveling salesman problem

First, I would like to explain the traveling salesman problem. The traveling salesman problem is "the problem of finding the shortest distance to go through all cities once and return to the starting point." For example, given the four cities A, B, C, and D as shown in the following figure, the shortest distance is as follows. jyunkai.png In this case, the number of cities is small and the arrangement is easy, so the shortest distance can be found immediately, but as the number of cities increases a little, the shortest distance cannot be seen at a glance. jyunkai2.png

Actually, the number of combinations of visit order of cities is (number of cities-1) !, so if it is 30 cities or 40 cities, there will be a great number of visit orders, so it is very difficult to check all visit orders. Is difficult, isn't it? Therefore, I tried to solve this combination optimization problem appropriately so that I could get a good solution.

SA method

The most orthodox general-purpose method for solving combinatorial optimization problems is the SA method (Simulated-anearing). Translated into Japanese, it is Simulated Annealing. Simulated Annealing sets the energy function $ E $. The energy function in this case is the total distance that goes around all the cities and returns to the starting point. We will change the order of visits to the cities so that this total distance will decrease. jyunkai3.png

In this case, the total distance will be longer than before the replacement. Therefore, I would like to reject this replacement, but in the SA method, the replacement is adopted even if the distance becomes long with the probability shown in the following formula.

p=\exp \left\{-\frac{1}{T}\left(E_{t+1}-E_{t}\right)\right\}

Here, $ T $ in the formula is a parameter called temperature, which acts as noise. The larger $ T $, the easier it is to adopt replacement. Also, if the distance becomes smaller as a result of replacement, it will be adopted with a probability of 1. This method is called the metropolis method.

To summarize the flow so far

  1. Change the order of visits
  2. Decide whether to adopt replacement by the Metropolis method After that, repeat this for the number of cities and make it one step. It will be.

Furthermore, in the SA method, the temperature $ T $ is gradually reduced according to the number of steps by the following formula.

T_{t+1}=\alpha T_{t}

Also, $ \ alpha $ takes the value of $ 0 <\ alpha <1 $ as an indicator of how quickly the temperature drops.

By lowering the temperature in this way, replacement is more likely to occur when the number of steps is small, and replacement is less likely to occur when the number of steps is large, that is, at the end of the game. Then, the first person can get out even if he / she falls into a local solution, and the last person will update so that he / she can decide exactly. This comes from the fact that the strength of steel is increased by starting from a high temperature and gradually cooling it during blacksmithing, which is why it is called Simulated Annealing. This method is very versatile in combinatorial optimization problems.

Algorithm rebellion

Well, this is where I wanted to talk. I was trying to solve the traveling salesman problem using Simulated Annealing above. The algorithm itself is simple, so I was able to build it immediately. That's why I tried to test it with a very simple city layout. The easiest city layout is probably circular. jyunkai4.png

This time, I would like to talk on the premise of this urban layout. By the way, as a review, the procedure of Simulated Annealing is

  1. Change the order of visits
  2. Decide whether to adopt replacement by the Metropolis method After that, repeat this for the number of cities and make it one step. It was. First, let's change the order of visiting the cities. The initial value is given as a random number. jyunkai5.png

You can see at a glance that it is not the shortest route. Next, let's decide the city to replace. This is also decided by random numbers, but at that time I was using C language, so by multiplying the uniform distribution that outputs the value of $ 0 <x <1 $ by the number of cities $ N $ (int type), $ 0 I was generating an int type random number of <x <N-1 $. This was the cause of the rebellion.

A simulation was performed using a random number that outputs this self-made integer. At first, the number of steps was about 100, but it works well enough. It's been replaced 600 times. But one day an incident happened. I decided to try 10000 steps (replacement 60,000 times) by chance, and turned the simulation. jyunkai5.png I got this result. I thought it was strange. I thought something was wrong when the total distance was output, which was less than the theoretical value obtained by manually calculating the coordinates. That's natural, because ** the city has disappeared **. Moreover, this result is not reproducible, and one city disappears, two cities disappear, and nothing disappears. If you are accustomed to programming, you may immediately think that "the result is not reproducible = the random number is bad", but at that time I started out and got stuck in my head. At the end of the story, I even had a delusion like a certain Skynet, saying, "Isn't the algorithm judging that the city should disappear in order to get the shortest distance?"

ending

As some of you may have noticed, the punch line was that 1 itself was output due to rounding error in the part that outputs a uniform random number of $ 0 <x <1 $. If 1 is output, it will be replaced with the Nth city. Since the array of c language is counted from 0, the array is prepared from 0 to N-1, so the Nth array does not exist. As a result of acquiring nothing from there and replacing it, the city disappeared.

SOM and Traveling Salesman Problem

From this point onward, it will be a completely different story, but recently I tried the traveling salesman problem again, so I would like to post the results there. This time I tried to solve the traveling salesman problem using SOM. SOM is very simply observation data $ \ mathbf {Y} = (\ mathbf {y} _1, \ mathbf {y} _2, ... \ mathbf {y} _N) \ in \ mathbb {R} ^ {D × N} $ and latent variables

$ \ mathbf {X} = (\ mathbf {x} _1, \ mathbf {x} _2, ... \ mathbf {x} _K) \ in \ mathbb {R} ^ {M × K} $ Mapping $ \ mathbf {f}: \ mathbb {R} ^ {M} → \ mathbb {R} ^ {D} $

For detailed algorithms, there is Document written by a professor in our laboratory, so please refer to that. I'm happy. For example, if you give 3D data and 2D latent space, you will get this result. Som.gif

The figure on the left shows the observation data and mapping, and the figure on the right shows the latent space.

You can also get this result by giving 2D data and 1D latent space. Som_TSP_Node48.gif Click here for an image showing the final result ああああああ.gif

By the way, the data given is the coordinate data of 48 cities in the United States called att48. It's kind of regrettable. I feel like I can do it. The bottleneck here is that the start and end points are not connected. If you can do that, you're done. Therefore, I will give a latent variable that connects the start point and the end point. That's right, we place latent variables in a circle on the latent space. Here is the source code and execution result of the SOM with the circular latent variables placed.

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from numpy.random import*
import matplotlib.animation as anm
import math
import random

%matplotlib nbagg
fig = plt.figure(1,figsize=(6, 4))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
class som():
    def __init__(self,N,K,sigmin,sigmax,r,w,zeta):
        self.N = N
        self.K = K
        self.sigmin = sigmin
        self.sigmax = sigmax
        self.r = r
        self.w = w
        self.zeta = zeta
        self.hist =[]
        self.hist2 = []
        self.hist3 = []
        self.sigma=0
    def fit(self,T):
        for i in range(T):
            #print(i)
            if i==0 :
                self.kn=np.random.randint(self.K,size=(self.N))
            #if i%10 == 0 :
            #    self.sigma = self.sigma+0.5
            self.sigma=max(self.sigmax-(i/self.r)*(self.sigmax-self.sigmin),self.sigmin)
            #print(self.kn.shape)
            self.test = (self.zeta[self.kn][None,:,:]-self.zeta[:,None,:])**2
            #print("self,test=",self.test.shape)
            #print("self,zeta=",self.zeta[self.kn][None,:,:].shape)
            #print("self,zeta=",self.zeta[:,None,:].shape)
            self.h_kn= np.exp(-1/(2*self.sigma**2)*np.sum((self.zeta[self.kn][None,:,:]-self.zeta[:,None,:])**2,axis=2))
            #print(self.h_kn)
            self.g_k = np.sum(self.h_kn,axis=1)
            self.y_test =(1/(self.g_k[:,None]))*self.h_kn @ self.w 
            #print("y_test:{}",self.y_test.shape)
            self.kn= np.argmin(np.sum((self.w[:,None,:]-self.y_test[None,:,:])**2,axis=2),axis=1)
            self.X1=np.reshape(self.y_test,(K,2))
            self.hist.append(self.X1)
            self.hist2.append(self.kn)
            self.hist3.append(self.sigma)
            print(self.sigma)
        print(np.array(self.hist).shape)
        return self.hist,self.hist2,self.hist3
            #self.hist.setdefault('y',[]).append(self.y_test) 
            #self.hist.setdefault('k',[]).append(self.kn) 
    #def history(self):
        #return self.hist
a=0.1
c=-0.1
N=48#The number of data
D=2#Number of dimensions of data
L=1#Number of dimensions of latent space
K=300#Number of nodes
T=100
zeta=np.zeros((K,2))
distance = 0
sigmax=5.0#Initial value of sigma
sigmin=0.01#Minimum value of sigma
latent_space=np.random.normal
r=T-5#Tilt
rad = np.linspace(-np.pi, np.pi, K)
zetax = np.sin(rad)+1
zetay = np.cos(rad)+1
zeta[:,0] = zetax
zeta[:,1] =zetay
train_data = np.loadtxt('att48.txt')
w = np.reshape(train_data,(48,D))
i=0
SOM = som(N,K,sigmin,sigmax,r,w,zeta)
kekka,k,sigma = SOM.fit(T)
kekka = np.array(kekka)
k = np.array(k)
sigma = np.array(sigma)
#print(k.shape)
for i in range(K):
    #print(("x1",kekka[T-1,i,0]))
    #print(("x2",kekka[T-1,i+1,0]))
    #print(("y1",kekka[T-1,i,1]))
    #print(("y2",kekka[T-1,i+1,1]))
    #print(("x_delta",kekka[T-1,i,0]-kekka[T-1,i+1,0]))
    #print(("y_delta",kekka[T-1,i,1]-kekka[T-1,i+1,1]))
    if(i==K-1):
        distance += np.sqrt(abs(kekka[T-1,i,0]-kekka[T-1,0,0])**2+abs(kekka[T-1,i,1]-kekka[T-1,0,1])**2)
    else:
        distance += np.sqrt(abs(kekka[T-1,i,0]-kekka[T-1,i+1,0])**2+abs(kekka[T-1,i,1]-kekka[T-1,i+1,1])**2)
    #print("delta",np.sqrt(abs(kekka[T-1,i,0]-kekka[T-1,i+1,0])**2+abs(kekka[T-1,i,1]-kekka[T-1,i+1,1])**2))
    print(distance)
#print(kekka.shape)
def update(i, zeta,w):
#for i in range(T):
    print(i)
    ax1.cla()
    ax2.cla()
    ax1.scatter(w[:,0],w[:,1],s=100,color='r')
    ax1.scatter(kekka[i,:,0],kekka[i,:,1],color = 'k',marker = '*')
    ax1.plot(kekka[i,:,0],kekka[i,:,1],color='k')
    ax1.plot(kekka[i,:,0].T,kekka[i,:,1].T,color='k')
    ax2.scatter(zeta[k[i,:]][:,0],zeta[k[i,:]][:,1],c=w[:,0])
    ax1.set_xlabel("X_axis")
    ax1.set_ylabel("Y_axis")
    ax2.set_title('latentspace')
    ax2.set_xlabel("X_axis")
    ax2.set_ylabel("Y_axis")

ani = anm.FuncAnimation(fig, update, fargs = (zeta,w), interval = 100, frames = T-1,repeat = False)


ani.save("RSom_TSP_Node300.gif", writer = 'imagemagick')

RSom_TSP_Node300.gif

If it is a gif, the final result will disappear in an instant, so I also prepare an image. aaa.jpg it's here

It feels really good. The calculation time is fast, so I think it's a pretty good line. After that, I think it will be even better if we can get out of the local solution by adding noise.

Recommended Posts

The story of the algorithm drawing a ridiculous conclusion when trying to solve the traveling salesman problem properly
Try to solve the traveling salesman problem with a genetic algorithm (Theory)
Try to solve the traveling salesman problem with a genetic algorithm (Python code)
Try to solve the traveling salesman problem with a genetic algorithm (execution result)
I tried "Implementing a genetic algorithm (GA) in python to solve the traveling salesman problem (TSP)"
The story of trying to reconnect the client
I want to solve the problem of memory leak when outputting a large number of images with Matplotlib
I tried to implement the traveling salesman problem
Solve the Python asymmetric traveling salesman problem with a branch and bound method
A story of a high school graduate technician trying to predict the survival of the Titanic
A story that failed when trying to remove the suffix from the string with rstrip
A story that got stuck when trying to upgrade the Python version on GCE
I can't find the clocksource tsc! ?? The story of trying to write a kernel patch
Zip 4 Gbyte problem is a story of the past
A memo of misunderstanding when trying to load the entire self-made module with Python3
A story about trying to introduce Linter in the middle of a Python (Flask) project
Finding a solution to the N-Queen problem with a genetic algorithm (2)
A story about how to deal with the CORS problem
Finding a solution to the N-Queen problem with a genetic algorithm (1)
About the traveling salesman problem
Let's write a program to solve the 4x4x4 Rubik's Cube! 2. Algorithm
A story about trying to automate a chot when cooking for yourself
A story that struggled to handle the Python package of PocketSphinx
I wanted to solve the ABC164 A ~ D problem with Python
About the Ordered Traveling Salesman Problem
A story about trying to improve the testing process of a system written in C language for 20 years
The story of writing a program
Use Raspberry Pi to solve the problem of insufficient mobile Wi-Fi connection
Solve a simple traveling salesman problem using a Boltzmann machine with simulated annealing
A story of a deep learning beginner trying to classify guitars on CNN
Try to solve a set problem of high school math with Python
The story of trying to push SSH_AUTH_SOCK obsolete on screen with LD_PRELOAD
A story that suffered from OS differences when trying to implement a dissertation
Want to solve a simple classification problem?
How to solve the bin packing problem
The story of adding MeCab to ubuntu 16.04
Python: I tried the traveling salesman problem
The story of trying deep3d and losing
The story of blackjack A processing (python)
The story of pep8 changing to pycodestyle
The story of IPv6 address that I want to keep at a minimum
Solving the traveling salesman problem with the genetic algorithm (GA) and its library (vcopt)
[Python] Solution to the problem that elements are linked when copying a list