Homotopic Convex Transformation: A New Landscape Smoothing Method for the Traveling Salesman Problem

TL;DR

This paper introduces a homotopic convex transformation method that smooths the TSP landscape, reducing local optima and enhancing heuristic search effectiveness, leading to significant performance improvements on large instances.

Contribution

It proposes a novel landscape smoothing technique for TSP using homotopic convex transformation and integrates it into heuristics to improve global search capability.

Findings

Smoothed TSP landscapes have fewer local optima.

Increased fitness distance correlation improves search.

Significant performance gains on large TSP instances.

Abstract

This paper proposes a novel landscape smoothing method for the symmetric Traveling Salesman Problem (TSP). We first define the Homotopic Convex (HC) transformation of a TSP as a convex combination of a well-constructed simple TSP and the original TSP. The simple TSP, called the convex-hull TSP, is constructed by transforming a known local or global optimum. We observe that controlled by the coefficient of the convex combination, with local or global optimum, (i) the landscape of the HC transformed TSP is smoothed in terms that its number of local optima is reduced compared to the original TSP; (ii) the fitness distance correlation of the HC transformed TSP is increased. Further, we observe that the smoothing effect of the HC transformation depends highly on the quality of the used optimum. A high-quality optimum leads to a better smoothing effect than a low-quality optimum. We then propose an iterative algorithmic framework in which the proposed HC transformation is combined within a heuristic TSP solver. It works as an escaping scheme from local optima aiming to improve the global search ability of the combined heuristic. Case studies using the 3-Opt and the Lin-Kernighan local search as the heuristic solver show that the resultant algorithms significantly outperform their counterparts and two other smoothing-based TSP heuristic solvers on most of the test instances with up to 20,000 cities.

Tables

Table 1. TABLE I: Test Instances

Instance	eil51	berlin52	st70	pr76	rat99	rd100	ch130	kroA150	gr96	brazil58	gr120	si175
Size	51	52	70	76	99	100	130	150	96	58	120	175
Edge type	EUC_2D	EUC_2D	EUC_2D	EUC_2D	EUC_2D	EUC_2D	EUC_2D	EUC_2D	GEO	EXPLICIT	EXPLICIT	EXPLICIT

Table 2. TABLE II: Costs of the selected optima, where x o ​ p ​ t subscript 𝑥 𝑜 𝑝 𝑡 x_{opt} and x L ​ O subscript 𝑥 𝐿 𝑂 x_{LO} denotes the global and local optimum, respectively.

Instance	$f\left(x_{opt}\right)$	$f\left(x_{LO}\right)$	excess (%)
eil51	426	428	0.4695
berlin52	7542	7902	4.7733
st70	675	684	1.3333
pr76	108159	109190	0.9532
rat99	1211	1217	0.4955
rd100	7910	8054	1.8205
ch130	6110	6363	4.1408
kroA150	26524	27040	1.9454
gr96	55209	55507	0.5398
brazil58	25395	25643	0.9766
gr120	6942	7066	1.7862
si175	21407	21485	0.3644

Table 3. TABLE III: Settings of λ 𝜆 \lambda in LSILS Implementations

Function Evaluation Period ($\times {10}^7$):	$0\to 1$	$1\to 2$	$2\to 3$	$3\to 4$	$4\to 5$	$5\to 6$	$6\to 7$	$7\to 8$	$8\to 9$	$9\to 10$
Setting1	0.02	0.02	0.02	0.02	0.02	0.02	0.02	0.02	0.02	0.02
Setting2	0.04	0.04	0.04	0.04	0.04	0.04	0.04	0.04	0.04	0.04
Setting3	0.06	0.06	0.06	0.06	0.06	0.06	0.06	0.06	0.06	0.06
Setting4	0.00	0.00	0.01	0.01	0.02	0.02	0.03	0.03	0.04	0.04
Setting5	0.00	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09

Table 4. TABLE IV: Areas Under the Excess Curves

Instance	ILS	GH	SSA	LSILS
				Setting1	Setting2	Setting3	Setting4	Setting5
rd400	6.930e+05	6.692e+05	7.005e+05	2.368e+05	1.684e+05	2.266e+05	3.900e+05	2.706e+05
p654	8.422e+05	6.072e+05	2.227e+05	5.856e+05	7.223e+05	8.285e+05	6.388e+05	6.084e+05
u724	1.307e+06	1.321e+06	1.283e+06	1.004e+06	6.724e+05	6.701e+05	1.085e+06	7.906e+05
pcb1173	2.149e+06	2.012e+06	1.910e+06	1.801e+06	1.832e+06	2.038e+06	1.862e+06	1.690e+06
rl1304	1.109e+06	1.093e+06	1.039e+06	7.521e+05	8.690e+05	1.229e+06	8.329e+05	6.904e+05
vm1748	1.922e+06	1.894e+06	1.767e+06	1.856e+06	1.734e+06	1.751e+06	1.820e+06	1.634e+06
u1817	2.753e+06	2.263e+06	2.801e+06	1.922e+06	2.121e+06	2.487e+06	2.106e+06	1.920e+06

Table 5. TABLE V: Mean CPU time

Instance	ILS	GH	SSA	LSILS
				Setting1	Setting2	Setting3	Setting4	Setting5
rd400	48.04s	81.37s	68.88s	42.75s	34.58s	28.55s	79.98s	65.87s
p654	7.96s	13.57s	10.85s	8.11s	8.14s	8.36s	16.38s	16.32s
u724	76.40s	138.50s	119.41s	60.86s	41.52s	28.71s	115.37s	85.00s
pcb1173	119.94s	229.56s	201.21s	56.18s	26.49s	18.59s	117.19s	80.67s
rl1304	73.89s	165.36s	127.96s	61.47s	43.28s	31.10s	109.97s	85.77s
vm1748	116.19s	229.42s	196.85s	78.50s	49.04s	31.64s	151.44s	103.60s
u1817	123.23s	230.96s	212.64s	53.15s	24.99s	17.79s	122.85s	79.61s

Equations

$$\begin{array}{rl}\text{mbox}minimize& f(x)=c_{x_{(n)},x_{(1)}}+\sum _{i=1}^{n-1}{c}_{x_{(i)},x_{(i+1)}},\\ \text{mbox}subjectto& x=(x_{(1)},x_{(2)},\dots ,x_{(n)})\in {\mathcal{P}}_n,\end{array}$$

$$d(x_1,x_2)=|\{e\in \mathcal{E}|e\in x_1\wedge e\notin x_2\}|$$

$$c_{i,j}^{\prime }(\alpha )=\{\begin{array}{rl}\bar{c}+(c_{i,j}-\bar{c}{)}^{\alpha }& \text{mbox}if{c}_{i,j}\ge \bar{c},\\ \bar{c}-(\bar{c}-c_{i,j}{)}^{\alpha }& \text{mbox}if{c}_{i,j}<\bar{c},\end{array}$$

$$\text{mbox}cityinterval=\frac{1}{n}\sum _{i=1}^n{c}_{i,m(i)},$$

$$m(i)=\underset{j\in \{1,\dots ,n\},j\ne i}{\mathrm{argmin}}{c}_{i,j}.$$

$$c_{i,j}^{\prime }(\lambda )=(1-\lambda )c_{i,j}+\lambda {\hat{c}}_{i,j}.$$

$$g(x|\lambda )=(1-\lambda )f_o(x)+\lambda f_c(x).$$

$$\text{mbox}FDC(f(x_{LO}),d_{opt})=\frac{\text{mbox}cov(f(x_{LO}),d_{opt})}{\sigma (f(x_{LO}))\sigma (d_{opt})},$$

$$\text{mbox}excess=\frac{f(x_{LO})-f(x_{opt})}{f(x_{opt})}\times 100\%,$$

Full text

I Introduction

As a well-known ${\cal NP}$-hard combinatorial optimization problem, the Traveling Salesman Problem (TSP) [1] has been extensively studied for years. Algorithms, including both exact and metaheuristics, have been developed in a number of literatures. It’s been well acknowledged that exact algorithms struggle to find the global optima of large size TSPs within a reasonable time. Metaheuristics, while not guaranteeing to find the global optimum, can provide near-optimal solutions within tolerable time for medium and large size TSPs. The hardness of the TSP is due to the ruggedness and irregularity of its fitness landscape. It is extremely hard for a metaheuristic to find the global optimum of a TSP when there are many local optima in its search landscape.

In this paper, we propose a novel method to smooth the TSP landscape. For a TSP with $n$ cities, we first construct a 2D Euclidean TSP in which the $n$ cities lie on the convex hull of the TSP graph and the order of the cities follows the order they appear in a known local optimum of the original TSP. The generated TSP is named as the convex-hull TSP. It can be proved that the convex-hull TSP is unimodal, which means that a search process starting from any initial solution will always end in the global optimum of the convex-hull TSP. The original TSP can then be smoothed by a convex combination of the convex-hull TSP and the original TSP with coefficient $\lambda\in[0,1]$. When $\lambda=0$, the transformed TSP is just the original TSP, while $\lambda=1$ means the convex-hull TSP. This actually defines a homotopic transformation from the original TSP to the transformed TSP. We thus call the proposed combination as the Homotopic Convex (HC) transformation.

To investigate the smoothing effects of the HC transformation, we conduct landscape analysis experiments on 12 TSP instances from the well-known TSPLIB [2]. In our experiment, the ruggedness of the TSP landscape is measured by four metrics, including local optimum density, escaping rate, Fitness Distance Correlation (FDC) and runtime. These metrics are computed from the results of 1000 runs of Iterated Local Search (ILS) [3] with the 3-Opt local search and the double bridge perturbation. In our experiments, different values of $\lambda$ are tested and their influence to the smoothing effect is analyzed.

We observe from the experimental results that the HC transformation exhibits the following three features. First, the landscape of the transformed TSP has less number of local optima than the original TSP, which implies a smoother TSP landscape. Second, in the transformed TSP, the FDC is higher than that of the original TSP. As stated in [4, 5], a high FDC implies that a search heuristic can find a path to the global optimum by observing the fitness change of the encountered local optima. Third, the strength of smoothing can be controlled by the coefficient $\lambda$. In addition, we observe that the runtime for ILS to find the global optimum is significantly reduced after the transformation.

We thereby make use of these features to propose a general iterative algorithmic framework. In the framework, the HC transformation is combined within a TSP heuristic. During the search, the heuristic works on the transformed TSP, while the transformed TSP is iteratively updated by the newly-found best solution with respect to the original TSP. As case studies, we embed the 3-Opt local search and the Lin-Kernighan (LK) local search [14] into the proposed framework. The instantiated algorithms are called Landscape Smoothing Iterated Local Search (LSILS) and LSILS-LK, respectively. In the experimental study, LSILS and LSILS-LK are compared against the original ILS and two existing landscape smoothing based algorithms on 17 TSP instances with up to 20,000 cities. Comparison results show that LSILS and LSILS-LK significantly outperform their counterparts on most of the test instances. Note that we do not claim the proposed LSILS is better than the state-of-the-art TSP solvers. Our goal is to show the potential of the HC transformation on improving the existing TSP heuristics.

The rest of this paper is organized as follows. Section II presents the related work. Section III presents the definition of the HC transformation and the experimental studies on the investigation of the characteristics of the HC transformation. In Section IV the framework and the instantiated algorithms are presented. Experimental results against its counterparts are also provided in this section. Section V concludes the paper and discusses future work.

II Related Work

Given a set of cities and the cost between every pair of cities, the TSP is to find the most cost-effective tour that visits every city exactly once and returns to the starting city. It is ${\cal NP}$-hard [6] and one of the most widely used testbeds in combinatorial optimization. The formal definition of the TSP can be stated as follows. Let ${\cal G}=({\cal V},{\cal E})$ be a fully connected graph with cities as vertexes, where $\cal V$ is the vertex set and $\cal E$ the edge set. Denote $c_{i,j}>0$ the cost of the edge between vertex $i$ and vertex $j$, the objective function of a TSP is defined as

[TABLE]

where $f:\mathcal{P}_{n}\to\mathbb{R}$ is the objective function and $\mathcal{P}_{n}$ is the permutation space of $\{1,2,\cdots,n\}$. In this paper we focus on symmetric TSPs, i.e., $c_{i,j}=c_{j,i}$ for all $i,j\in\{1,2,\dots,n\}$.

The hardness of a TSP is usually characterized by its fitness landscape. Formally, the fitness landscape of a combinatorial optimization problem is defined by $(\mathcal{S},f,d)$, where $\mathcal{S}$ is the solution space, $f:\mathcal{S}\to\mathbb{R}$ is the objective function (fitness) and $d$ is a distance measure on $\mathcal{S}$. For the TSP, the distance between any two solutions $x_{1}$ and $x_{2}$ is defined as the number of edges contained in one solution but not in the other:

[TABLE]

where $e\in x$ means edge $e$ is in the solution $x$.

A great amount of efforts have been made on the TSP fitness landscape analysis. Among these studies, Stadler and Schnabl [7] first investigated the landscape of symmetric and asymmetric TSPs by random walk. Boese [8] first observed that there is a strong correlation between the distance of a solution to the global optimum and its cost when he investigated the landscape of a specific TSP instance att532. He called this phenomenon the big valley structure. Hains et al. [9] confirmed the existence of the big valley structure on some TSP instances. However, they also argued that the big valley breaks down into several “funnels” around some local optima with function values very close to that of the global optimum. Ochoa and Veerapen [10] analyzed the big valley structure assumption by their local optima network model. They concluded that in the four studied instances the big valley decomposes into a number of sub-valleys. However, Ochoa and Veerapen’s conclusion depends strongly on the solution perturbation method used in their study. Shi et al. [11] gave a more detailed definition of the big valley structure. In their definition, a TSP instance exhibits a big valley structure if (i) all the global optima of a TSP are located in a relatively small region in the solution space (this defines the “bottom” of the big valley), and (ii) there is a strong correlation between the fitness of a solution and its distance to the nearest global optimum (this defines the “slope” of the big valley). Empirical analysis on the landscape of 10 TSPLIB instances showed that 9 of the 10 test instances meet the detailed definition of the big valley structure.

In the work of Fonlupt et al. [12], it is found that the TSP landscape defined by 2-Opt-move has a higher FDC value than the TSP landscape defined by city-swap. Merz and Freisleben [13] studied the ruggedness and the FDC of the TSP landscape using the 3-Opt local search and the Lin-Kernighan (LK) local search [14]. They observed that local optima are frequently close to each other and also close to the global optimum. Tayarani-N and Prügel-Bennett [15, 16] also used the 3-Opt local search as a tool to analyze the characteristics of the fitness landscape on different classes of TSPs. They found that high-quality local optima are more likely to be found by a local search process than low-quality local optima and the probability of finding a global optimum decreases exponentially with increasing problem size. They also found that the Euclidean TSPs have relatively high FDC values and the value is relatively low for random TSP. For more detailed survey, interested readers please see [17, 18, 19].

Existing TSP metaheuristics fall basically into four categories. One kind is constructive based where a solution to a TSP is constructed gradually from partial to whole. Examples of such including Ant Colony Optimization (ACO) [20], Greedy Randomized Adaptive Search Procedure (GRASP) [21] and others. The second kind is local search based. Local search first defines a neighborhood structure on the solution space of the TSP. Starting from an initial solution, local search iteratively explores the neighborhood of the current solution for a better solution. Local search can only guarantee to find local optima. The most widely-used local search for the TSP is the $k$-Opt local search ($k\geq 2$), in which the neighborhood structure is defined by edge exchange. In each descent move of a $k$-Opt local search, $k$ edges are replaced by another $k$ edges. The well-known 2-Opt local search, 3-Opt local search [1] and the LK local search all belong to this class. The third category is the population based metaheuristics, e.g. Genetic Algorithm (GA) [22], Particle Swarm Optimization (PSO) [23] and others. In this category, solutions are created through recombination operators working in the permutation space. The fourth category is the Memetic Algorithm (MA) [24], in which the local search based heuristics are hybridized with heuristics from other categories.

To escape from local optima, existing local search based metaheuristics for the TSP are usually embedded with global optimization strategies such as criterion, solution or problem perturbation and landscape smoothing. In criterion perturbation, rather than only accepting better solution in the neighborhood, a worse solution could be accepted in probability (e.g. Simulated Annealing [25, 26]) or if it prevents from returning to previously visited areas (e.g. Tabu Search [27]). In the solution perturbation methods, the current local optimum is perturbed to obtain a new solution and the local search process restarts from that new solution. The well-known ILS [3, 28] belongs to this category.

The problem perturbation methods modify the problem itself so that the current local optimum is no longer locally optimal to the modified problem. In the Guided Local Search (GLS) proposed by Voudouris and Tsang [29], the original TSP is changed by adding penalties on some selected edges. Whenever the GLS is trapped in a local optimum, edges from the local optimum are evaluated and penalized so that the local optimum is no longer optimal in the next round of search. Shi et al. [11] improved the GLS based on the big valley assumption. Their GLS maintains an elite solution and the edges in the elite solution will be protected from the penalization. Walshaw [30] introduced a multilevel refinement approach to simplify the TSP. In each level, several selected cities are matched and the edges between them is fixed so that the problem is coarsened to an easier problem. In the backbone-guided local search proposed by Zhang and Looks [31], the edge cost is reduced proportionally to the frequency that the edge appears in the historical local optima set.

Very few global search approaches have been proposed based on landscape smoothing. Gu and Huang [32] proposed to smooth the landscape of a TSP by edge cost manipulation. They normalize and transform the edge cost to

[TABLE]

where $\bar{c}$ is the average cost of all edges and $\alpha\geq 1$ is the smoothing factor. When $\alpha=1$, it is just the original TSP. Note that in Eq. (3) the edge costs are normalized, hence $\lim_{\alpha\rightarrow+\infty}c^{\prime}_{i,j}=\bar{c}$. That is, with increasing $\alpha$, all edge costs approach to a fixed value $\bar{c}$. This implies that all solutions will have the same fitness, i.e., the TSP landscape is smoothed to a plane, as sketched in Fig. 1.

Schneider et al. [33] proposed several new smoothing functions (exponential, hyperbolic, sigmoidal, and logarithmic) for the TSP. Based on the work of Gu and Huang and Schneider et al., Coy et al. [34] carried out systematic experiments to investigate the smoothing effect of the smoothing functions to the fitness landscape. They claimed that the smoothing method can help local search escape from poor local optima by moving uphill occasionally and by not taking a downhill move occasionally. Further study by Coy et al. [35] shows that a sequential smoothing algorithm with alternated convex and concave smoothing function can achieve satisfactory performance. Dong et al. [36] proposed a stochastic local search metaheuristic based on the exponential smoothing method proposed by Schneider et al.. Hasegawa and Hiramatsu [37] combined Gu and Huang’s smoothing method within a MA and found that it is beneficial to apply the smoothing method in the MA.

III Homotopic Convex Transformation

In this section, we present a new landscape smoothing method called the Homotopic Convex (HC) Transformation. It transforms the original TSP by combining it with a well-designed unimodal TSP, which is named as the convex-hull TSP. In the following, we first define the convex-hull TSP and the HC transformation. The effects of the HC transformation to landscape smoothing are illustrated afterwards.

III-A Convex-Hull TSP

Although the TSP is ${\cal NP}$-hard, some TSP instances are easy to solve. For example, a 2D Euclidean TSP with all the cities lying on the convex hull of the TSP graph, which we name it as the convex-hull TSP, is very easy to solve [38]. Fig. 2 shows an example of the convex-hull TSP and its global optimum. The only global optimum of a convex-hull TSP is the convex hull itself.

Here we prove that the convex-hull TSP is unimodal for any $k$-Opt local search, i.e., the TSP tour formed by the edges on the convex hull is the only $k$-optimal tour [39] of such TSP.

Definition**.**

A tour is said to be $k$-optimal if it is impossible to obtain a tour with smaller cost by replacing any $k$ of its edges by any other set of $k$ edges.

The property of the $k$-optimal tour is summarized in Theorem 1 reproduced from [39].

Theorem 1**.**

Denote $C_{k}$ the set of all $k$-optimal tours, we have $C_{1}\supset C_{2}\supset\dots\supset C_{n}$. In other words, a $k$-optimal tour is also a $k^{\prime}$-optimal tour for $k^{\prime}<k$.

Based on Theorem 1, we prove that the convex-hull TSP is unimodal to any $k$-Opt local search. Theorem 2 summarizes the result.

Theorem 2**.**

For a convex-hull TSP, let $x_{c}$ denote the convex hull tour, then $x_{c}$ is the only $k$-optimal tour in the solution space of the convex-hull TSP for any $k\in\{2,3,\dots,n\}$.

Proof.

It is sure that $x_{c}$ is $k$-optimal for any $k\in\{2,3,\dots,n\}$. We here prove the uniqueness of $x_{c}$ by contradiction. Assume there is another $k$-optimal tour $x_{\ast}$ for $f_{c}$. By Theorem 1, we known that $x_{\ast}$ is also 2-optimal. Since $x_{\ast}\neq x_{c}$, there are at least two edges in $x_{\ast}$ that cross each other. According to the triangle inequality, it is always possible to replace the crossed edges with non-interacting edges which are less cost. This indicates that $x_{\ast}$ is not 2-optimal, which contradict the assumption. ∎

III-B Definition of the HC Transformation

Given a TSP $f_{o}$ 111In the sequel, without causing confusion, we do not differentiate the name of a TSP and its objective function. with known local or global optimum $x^{*}$, a convex-hull TSP, denoted by $f_{c}$, is firstly constructed based on $x^{*}$. The construction makes sure that $f_{c}$ has the same optimum $x^{*}$ as $f_{o}$. Assume that the original TSP $f_{o}$ has $n$ cities, to construct the corresponding convex-hull TSP $f_{c}$, the HC transformation simply let $n$ cities uniformly distribute on a circle following the order they appear in $x^{*}$ of $f_{o}$. To level the scale of $f_{c}$ and $f_{o}$, the city interval on the circle of $f_{c}$ is set to be

[TABLE]

where $c_{i,m(i)}$ denotes the edge cost in the original TSP $f_{o}$ and $m(i)$ denotes the index of the nearest city to city $i$ in all the cities of $f_{o}$, i.e.,

[TABLE]

The city interval is used to determine the distance between adjoint pair of cities in the circle. Once all the cities are arranged on the circle, the distance, denoted as $\hat{c}_{i,j}$, between each city pair $(i,j)$ is computed and considered as the edge cost in the convex-hull TSP. Fig. 3(b) shows an example convex-hull TSP (and its global optimum) constructed for the TSP instance eil51 (Fig. 3(a)).

Once the convex-hull TSP is constructed, the HC transformation is defined by combining the convex-hull TSP with the original TSP by a coefficient $\lambda\in[0,1]$. In the transformed TSP, the edge cost $c^{\prime}_{i,j}$ between city $i$ and city $j$ is calculated by

[TABLE]

The objective function of the transformed TSP, denoted as $g$, can then be expressed as

[TABLE]

It is seen that when $\lambda=0$ the transformed TSP $g$ degenerates to the original TSP $f_{o}$ and when $\lambda=1$ it is smoothed to the convex-hull TSP $f_{c}$. With increasing $\lambda$, the original TSP $f_{o}$ can continuously deform into the convex-hull TSP $f_{c}$. This actually defines a homotopy transformation from $f_{o}$ to $f_{c}$ (cf. homotopy in topology [40]). We thus call it the homotopic convex transformation. The sketch of the HC transformation is illustrated in Fig. 4.

Comparing Fig. 1 with Fig. 4, it is seen that in the TSP smoothing approach proposed in [32], with $\alpha\rightarrow\infty$, the landscape becomes flat. While in the HC transformation, with $\lambda\rightarrow 1$, the landscape becomes unimodal. The key difference is that the optimum can be preserved after smoothing by the HC transformation. Note that the transformed TSP is still $\cal{NP}$-hard. For a detailed description, an example of the HC transformation can be found in the supplementary material of this paper.

III-C Effects of the HC transformation

To illustrate the effects of the HC transformation on landscape smoothing, in the following we conduct landscape analysis experiments to study the changes of the TSP landscape after the HC transformation. In the experiments, we execute ILS on the transformed TSP and use the following four metrics to characterize the landscape:

• •

Local Optimum Density: It is the number of local optima encountered by an ILS process per 100 moves. Here one move means the local search algorithm moves from the current solution to a new solution.

• •

Escaping Rate: It is the success rate that a new local optimum is reached by ILS starting from a new solution obtained by perturbing the current local optimum.

• •

Fitness Distance Correlation (FDC): To calculate the FDC, 1000 local optima (denoted as $x_{LO}$) obtained by ILS are randomly selected. Their function values $f(x_{LO})$ and their distances to the nearest global optimum $d_{opt}$ are counted. Then the FDC is define as

[TABLE]

where $\mbox{cov}(\cdot)$ denotes the covariance and $\sigma(\cdot)$ denotes the standard deviation.

• •

Runtime: It is the average runtime (CPU time) for ILS to find the global optimum.

The first two metrics, local optimum density and escaping rate, measure the ruggedness of the landscape of a TSP. A lower local optimum density means a smoother TSP landscape. It is intuitive that the escaping rate is negatively correlated to the average diameter of the attractive basins of local optima. A large size of attractive basin implies a few local optima on the TSP landscape. Hence a lower escaping rate means a smoother TSP landscape. The third metric, FDC, measures the overall trend of the TSP landscape. A higher FDC means a more regular TSP landscape which implies that a search algorithm can find a path to the global optimum by observing the fitness change of the encountered local optima [4, 5]. The last metric, runtime, directly reflects the hardness of a TSP for ILS. A larger runtime means a harder TSP.

In the experiments, ILS is used to sample local optima on the TSP landscape. As shown in Alg. 1, ILS iteratively executes a local search procedure and a perturbation procedure till the stopping criterion is met. In our experiments, the 3-Opt local search and double bridge perturbation (please see Fig. 5 for its illustration) are selected as the local search method and the perturbation method of ILS, respectively. In Alg. 1, LocalSearch($x_{j}^{\prime}$, $x_{best}$) takes $x_{j}^{\prime}$ and the current best $x_{best}$ as input and returns a local optimum $x_{j+1}$ and a new $x_{best}=\arg\min\{f(x_{j+1}),f(x_{best})\}$.

A total of 12 instances are selected from the TSPLIB as the test instances. The characteristics of the test instances are shown in Table I. The third column of Table I shows the edge cost type of each test instance, in which EUC_2D means the edge costs are Euclidean distances in a 2D space, GEO means the edge costs are geographical distances and EXPLICIT means the edge costs are listed explicitly in the TSP description file. The reasons we select these 12 instances are: (1) their sizes are not too small to lose the complexity and ruggedness of the landscape; (2) their sizes are not too large to find their global optima; (3) they belong to different TSP types. To conduct the analysis, $\lambda$ is set to be $0,0.01,\dots,0.09,0.1$ respectively, for each TSP instance. Here we set $\lambda\leq 0.1$ since our preliminary results show that when $\lambda>0.1$ the property of the transformed TSPs are nearly the same to the property of the convex-hull TSP, which means ILS can find the global optimum immediately.

In the following, we investigate the effects of the HC transformation by using either global optimum or local optimum. The code is implemented in GNU C++ with O2 optimizing compilation. The computing platform is two 6-core 2.00GHz Intel Xeon E5-2620 CPUs (24 Logical Processors) under Ubuntu OS.

III-D Effects of the HC transformation with global optimum

Since the test instances are from the TSPLIB, their global optima are known. First, we use the global optima of the test instances to construct the convex-hull TSPs and conduct HC transformation with different $\lambda$ values. On each transformed TSP, 1000 runs of ILS are executed from random solutions and terminate only when the global optima are reached. We record all the local optima ever encountered by ILS and other relating information in the entire search history. Then we compute the four metrics for each transformed TSP, respectively.

Fig. 6(a) shows how the local optimum density of the transformed TSP landscape changes against $\lambda$. From Fig. 6(a) we can see that, on the instances with Euclidean edge cost type (denoted by “EUC”), in general the local optimum density decrease as $\lambda$ increases, which means that the HC transformation can smooth the TSP landscape and its smoothing effect is controlled by $\lambda$. However, on the instances with non-Euclidean edge cost type (“GEO” and “EXP”), the local optimum density is not always negatively related to $\lambda$. For example, it is seen that the local optimum density of the TSP landscape for brazil58 increases when $0.01\leq\lambda\leq 0.08$.

The escaping rates of ILS on different transformed TSPs are shown in Fig. 6(b). As stated before, a lower escaping rate means a smoother TSP landscape. From Fig. 6(b) we can see that, on all of the Euclidean TSP instances, in general the escaping rate decrease as $\lambda$ increases, which exhibits similar phenomena as the local optimum density. On the non-Euclidean TSP instances, it seems that the escaping rate also decrease as $\lambda$ increases. However, the decline of the escaping rate is not very obvious on brazil58 and si175 when $\lambda\geq 0.01$. Based on the above observations, we may conclude that the HC transformation can indeed smooth the landscapes of all the Euclidean TSP and some non-Euclidean TSP instances.

Fig. 6(c) shows how the FDC of the transformed TSP changes against $\lambda$. From Fig. 6(c) we can see that, although the FDC curve has fluctuations, in general the FDC increases as $\lambda$ increases on most test instances. An exception is found on the non-Euclidean instance si175, on which the lowest FDC value appears when $\lambda=0.02$ and when $\lambda\geq 0.04$ the FDC is always 1. Based on the above observations, we may conclude that the HC transformation can increase the FDC value on most of the test instances.

On different transformed TSPs, the runtime for ILS to find the global optimum is shown in Fig. 6(d). From Fig. 6(d) we can see that, in general, the runtime of ILS decreases as $\lambda$ increases on most of the test instances. An exception is found on the non-Euclidean instance brazil58, on which the runtime first increases then decreases. Based on the above observation, we conclude that the HC transformation can reduce the running time of ILS to find the global optimum simply because the fitness landscape of the transformed TSP is smoothed with fewer number of local optima and higher FDC than the original TSP. Of course, this claim is based on the premise that the instance size is relatively small.

In summary, we may conclude that in most cases the proposed HC transformation can indeed smooth the TSP landscape and the strength on the smoothing can be controlled by $\lambda$. This effect is more significant in the instances with the Euclidean edge cost type, which is the most common edge cost type in the TSP, than with other TSP types.

III-E Effects of the HC Transformation with local optimum

In practice, the global optimum is not known. We therefore investigate the effects of the HC transformation when a local optimum is used to construct the convex-hull TSP. The same test instances in Table I are used. For each instance, we run a round of 3-Opt local search to obtain a local optimum $x_{LO}$. Then we use this local optimum to construct the convex-hull TSP with different $\lambda$ values. The function value of the used local optima are listed in Table II. Table II also lists the excess of the used local optima, which is defined in Eq. (9):

[TABLE]

where $x_{opt}$ is the global optimum. The other experimental settings are the same as those in Section III-D.

Since the HC transformation is based on local optima, the global optima of the transformed TSPs may be different from the original global optima. To conduct landscape analysis experiments, we first use the LKH software [41] to obtain the global optima of the transformed TSPs. On each transformed TSP, we run the LKH 100 times and record the best solution as the global optimum. Since the LKH is one of the state-of-the-art algorithms for the TSP and the size of the transformed TSPs is relatively small, we believe the TSP solutions found by the LKH are indeed the global optima of the transformed TSPs.

The local optimum density values of different transformed TSPs are shown in Fig. 7(a). From Fig. 7(a) we can see that in 11 of the 12 test instances, the local optimum density is approximately negatively related to $\lambda$. Specifically, on seven instances (eil51, st70, pr76, rat99, rd100, gr96 and si175) the local optimum density decreases as $\lambda$ increases. On four instances (berlin52, ch130, kroA150 and gr120), the local optimum density first increases a little then starts to decrease.

Experimental results for the escaping rate are shown in Fig. 7(b). From Fig. 7(b) we can see that in 9 of the 12 test instances, the escaping rate is approximately negatively related to $\lambda$. Specifically, on five instances (eil51, rat99, ch130, gr96 and gr120), in general the escaping rate decreases as $\lambda$ increases in spite of the fact that on some instances the curve of the escaping rate is not very smooth. On four instances (st170, pr76, rd100 and kroA150), the escaping rate first increases a little then starts to decrease. Through comparing Fig. 7(a) and Fig. 7(b) with Fig. 6(a) and Fig. 6(b), we conclude that, on most of the test instances, the HC transformation based on local optima can achieve approximately the same smoothing effect as the HC transformation based on global optima.

Fig. 7(c) shows the FDC versus $\lambda$ when the HC transformation is based on local optima. From Fig. 7(c) we can see that in 9 of the 12 test instances, the FDC is approximately positively related to $\lambda$. Specifically, on four instances (eil51, kroA150, gr96 and brazil58), there is a general trend that the FDC increases as $\lambda$ increases, in spite of the fact that on some instances the FDC curve is not very smooth. On four instances (st70, rd100, ch130 and gr120), the FDC first decreases a little then starts to increase. On rat99, although the lowest FDC appears when $\lambda=0.05$, but in general the FDC is increase as $\lambda$ increases. Through comparing Fig. 7(c) to Fig. 6(c) we can conclude that, on most of the test instances, the HC transformation based on local optima can achieve approximately the same FDC-increasing effect as the HC transformation based on global optima.

Fig. 7(d) shows the runtime for ILS to find the global optimum when the HC transformation is based on local optima. From Fig. 7(d) we can see that in 10 of the 12 test instances, the runtime of ILS is approximately positively related to $\lambda$. Specifically, on five instances (eil51, rat99, ch130, gr120 and si175), the runtime decreases as $\lambda$ increases, while on another five instances (st70, pr76, rd100, kroA150 and gr96), the runtime first increases a little then starts to decrease. Through comparing Fig. 7(d) to Fig. 6(d) we can conclude that, in most of the test instances, the HC transformation based on local optima can achieve approximately the same runtime reducing effect as the HC transformation based on global optima.

Based on the above observations, we may conclude that, in most test instances, the HC transformation based on local optima can achieve approximately the same effects as the the HC transformation based on global optima.

III-F Influences of different local optima

In previous experiments, on some test instances, especially on berlin52, the HC transformation based on local optima does not achieve as good smoothing and FDC-increasing effect as the HC transformation based on global optima. Recall that in Table II, the local optimum used to construct the convex-hull TSP for berlin52 has a relatively low quality compared to the local optima of other instances. We thus conjecture that the effects depend highly on the quality of the local optimum used in the HC transformation. To confirm, we conduct the following experiments.

For berlin52, using the 3-Opt local search, we find totally 89 local optima with different function values ranging from 7658 to 8385. Based on these 89 local optima and the global optimum (with function value 7542), the experimental configurations used in previous sections are applied on berlin52. The obtained results on the four metrics against different local optima (and the global optimum) and $\lambda$ values are shown in Fig. 8 while Fig. 8(e)-8(h) show the cubic spline interpolation of the results.

From the results we can see that, when the local optima have relatively high quality (the blue curves in Fig. 8(a)-8(d)) the general tendency of the curves achieved by the local optima are similar to that achieved by the global optimum (the curves with the lowest function value in Fig. 8(a)-8(d)), while when the local optima have relatively low quality (the red curves in Fig. 8(a)-8(d)), the difference between the curves achieved by local optima and global optimum are quite clear. The cubic spline interpolation results (Fig. 8(e)-8(h)) also show the same tendency. This means that the HC transformation based on a high-quality local optimum can achieve approximately the same transformation performance as the HC transformation based on the global optimum. Hence, we conclude that the HC transformation using a high-quality local optimum can smooth the TSP landscape similar to as using the global optimum.

IV Algorithmic Framework

We claim that the proposed HC transformation can be used to improve the global search ability of existing TSP heuristics. To justify, in this section, we first propose a general framework in which the HC transformation is combined with a local search for the TSP. Then two instantiated algorithms of the framework with the 3-Opt local search and the LK local search, respectively, are proposed and compared against some known TSP smoothing algorithms. Note here that our objective is not to propose a state-of-the-art algorithm for the TSP, but to show that the proposed HC transformation can be used to improve existing TSP heuristics.

IV-A The general framework

The general framework is summarized in Alg. 2. Similar to solution perturbation based algorithms, the framework iteratively executes a local search procedure and a perturbation procedure. The key is that the local search is performed on the transformed TSP $g=(1-\lambda)f_{o}+\lambda f_{c}$. If we set $\lambda=0$, then $g=f_{o}$ and Alg. 2 degenerates to ILS in Alg. 1. A local optimum, denoted as $x_{0}$, is first obtained by any local search algorithm (line 2) from a random initial solution $x_{ini}$ (line 2). The algorithm repeats the following procedure until stop. At each iteration, the current best solution w.r.t. $f_{o}$, denoted as $x_{f_{o},best}$, is used to construct a convex-hull TSP (line 2) for a transformed TSP $g$ (line 2). The local search is then applied on the transformed TSP from a perturbed solution (line 2) for a new best $x_{f_{o},best}$. The coefficient $\lambda$ is updated (line 2). Note that LocalSearch($x^{\prime}_{j},x_{f_{o},best}|g$) means executing a local search from $x^{\prime}_{j}$ on the transformed TSP $g$, meanwhile keeping updating $x_{f_{o},best}$ by tracking the value change of $f_{o}$. It returns a local optimum $x_{j+1}$ w.r.t. the transformed TSP $g$ and the updated $x_{f_{o},best}$ w.r.t. the original TSP $f_{o}$. In the next iteration, $x_{j+1}$ will be perturbed to generate an initial solution for the next local search. $\lambda$ is updated after local search (line 2).

IV-B Performance Comparison based on 3-Opt Local Search

We first use the 3-Opt local search as the local search operator and the double bridge perturbation as the perturbation operator to realize the proposed framework. The resultant algorithm is named as Landscape Smoothing Iterated Local Search (LSILS). The developed LSILS is compared against ILS, the smoothing algorithm proposed by Gu and Huang [32] (denoted as GH) and the sequential smoothing algorithm proposed by Coy et al. [35] (denoted as SSA).

Seven instances are chosen from the TSPLIB as the test instances, including rd400, p654, u724, pcb1173, rl1304, vm1748 and u1817. Each algorithm is executed 100 runs on each test instance from different random initial solutions. The stopping criterion of each run is $10^{8}$ function evaluations. For LSILS, we test five different settings of $\lambda$. In the first three settings constant values of $\lambda$ are applied during the search and in the last two settings, $\lambda=0$ at the beginning and increases during the search, as shown in Table III. The parameter settings are the same as in [35]. That is, for GH the smoothing factor $\alpha=6\to 5\to 4\to 3\to 2\to 1$ in the first six local search rounds and for SSA $\alpha=7\to 5\to 3\to 1$ in the first four local search rounds. In both GH and SSA, if $\alpha$ has dropped to $1$ and the function evaluation budget has not been run out, the algorithm will execute ILS on the original TSP landscape (i.e. $\alpha=1$) until the function evaluation budget is run out. The 3-Opt local search and the double bridge perturbation are also used in the implementations of ILS, LSILS, GH and SSA. In addition, the 3-Opt local search implementation is speeded up by the techniques of near neighbor search and don’t look bits [42]. The other experimental settings are the same to the settings in the previous section.

Fig. 9 shows the mean excess values to the best solutions found so far by different smoothing algorithms against time. To properly compare among TSP solvers, a reasonable method is to compare the areas under curves (AUC) [43, 44]. The AUC values of the excess curves are listed in Table IV. From Fig. 9 and Table IV we can see that, in six out of seven test instances, the minimum AUC is obtained by LSILS, which means that it performs better than ILS, GH and SSA on these six instances. From Fig. 9 and Table IV we can also see that, in four instances, LSILS with increasing $\lambda$ values (Setting5 in Table III) achieves the best performance. On p654 (Fig. 9(b)), SSA achieves the best performance.

However, we notice that on p654 the performance of LSILS is close to that of SSA at termination. We conjecture LSILS can perform better if there are more computational budget. To justify, we run LSILS on p654 with an increased function evaluation budget ($5\times 10^{8}$). The results are shown in Fig. 10, from which we can see that at termination LSILS outperforms SSA in all settings. This confirms our conjecture.

On average, we see that LSILS performs significantly better than ILS, GH and SSA in terms of excess on most of the test instances. This clearly indicates that the proposed HC transformation can truly improve the performance of ILS.

However, we can also see that, the performance of LSILS is relatively poor at the beginning of the search. In fact, in most cases, the descent rate of GH and SSA is much higher than that of LSILS at the beginning of the search. This is because that the smoothing effect of the HC transformation depends highly on the quality of the local optimum used to construct the convex-hull TSP. As stated in Section III-F, with a low-quality local optimum, the smoothing effect of the HC transformation is not as good as a high-quality local optimum. As search goes by, better solutions are found, the performance of LSILS becomes better and exceeds the performance of GH and SSA eventually in most cases.

Table V lists the mean CPU time used by each algorithm in the experiments. Table V shows that, LSILS with Setting3 (i.e. $\lambda$ is constant and equals to $0.06$) has the shortest CPU time in six out of seven test instance and on most instances, the runtime of LSILS with increasing $\lambda$ value (i.e. Setting4 and Setting5) is lower than that of GH and SSA. This means that LSILS with increasing $\lambda$ value can find better solutions and requires less time than those of GH and SSA in most cases.

IV-C Performance Comparison based on the LK Local Search

In previous sections, we only test the performance of the proposed framework based on the 3-Opt local search and the test instances is relatively small with city number $n<2000$. We conduct another experiment to test the performance of the proposed framework based on the most powerful TSP local search method, the Lin-Kernighan (LK) local search [14], on ten middle-size and large-size TSP instances. The resultant algorithm is named as LSILS-LK.

In our experiment, the implementation of the LK local search is from the Concorde software package 222http://www.math.uwaterloo.ca/tsp/concorde/. In the LK local search, the edge exchange is restricted in a sub-graph of the original TSP graph. In the sub-graph, each vertex (city) only connects with its 20 nearest vertexes (cities). Based on the LK local search and the double bridge perturbation, we compare LSILS-LK against ILS with LK local search (ILS-LK), GH with LK local search (GH-LK) and SSA with LK local search (SSA-LK). In ILS-LK, GH-LK and SSA-LK, the LK local search is used as the local search and double bridge perturbation is used as the perturbation method.

Further, to test the perturbation strategy to the performance of LSILS-LK, we include another two algorithm instantiations, i.e. LSILS-LK with three double bridge perturbation (LSILS-LK-3DBP) and ILS-LK with three double bridge perturbation (ILS-LK-3DBP). The difference between LSILS-LK-3DBP and LSILS-LK is that in LSILS-LK-3DBP the double bridge perturbation is replaced by three successive double bridge perturbations. Here we do not apply the 3DBP strategy to GH-LK and SSA-LK, because except the first few iterations, these two algorithms are exactly the same to ILS-LK.

Six middle-size and large-size instances from the TSPLIB: {pcb3038, fnl4461, pla7397, rl11849, usa13509, d18512} and four randomly generate instances: {rand5000, rand10000, rand15000, rand20000} are used as benchmarks. Among the selected TSPLIB instances, pcb3038, pla7397 and rl11849 come from printed circuit board, while fnl4461, usa13509 and d18512 come from real world maps. Since it is very hard to count the number of function evaluations in the LK local search, we use the CPU runtime as the stopping criterion for all the algorithms. The maximal runtimes in seconds allowed for optimizing the test instances are {pcb3038: 600s, fnl4461: 900s, pla7397: 1500s, rl11849: 2400s, usa13509: 2700s, d18512: 3700s, rand5000:1000s, rand10000:2000s, rand15000:3000s, rand20000:4000s}, respectively. For LSILS-LK and LSILS-LK-3DBP, we use Setting4 and Setting5 since with these settings the proposed framework has shown its superiority in previous study.

The time point to increase $\lambda$ is based on the CPU time. For example, for LSILS-LK running on pcb3038 with Setting5, $\lambda$ will change from 0 to 0.01 at the 60 second, from 0.01 to 0.02 at the 120 second and so on. Further, to save computational effort, in the implementations of LSILS-LK and LSILS-LK-3DBP the transformed TSP $g$ (which guides the LK local search) is only updated at the time points when $\lambda$ increases (i.e., each time $\lambda$ increases, $g$ will be updated based on the new $\lambda$ value and the current best solution $x_{f_{o},best}$). Each algorithm is executed 50 runs on each test instance from different random initial solutions. The other experimental settings are the same to the previous experiments.

Fig. 11 shows the mean excess values of the best-so-far solutions found by the compared algorithms against time. Since the globally optimal function values of the randomly generated instances are unknown, we use the best found function value $\times 0.99$ as the globally optimal function value to calculate the excess on the random instances. From Fig. 11 we can see that, on all the test instances, the best performance is achieved by LSILS-LK-3DBP with Setting5.

Fig. 11 also shows that the curves of LSILS-LK and LSILS-LK-3DBP are regularly jagged. We observed that a jag happens when the $\lambda$ increases. This phenomenon confirms again that the HC transformation with better local optima can result in better smoothing effect and hence a higher $\lambda$ value is preferred. We can also observe that, although three double bridge perturbation reduces the performance of ILS-LK, in LSILS-LK, when the perturbation strength becomes larger, the algorithm performance significantly improves. It means that LSILS-LK can work better with three double bridge perturbation strategy. A possible reason is that after smoothing, the number of local optima in the TSP landscape reduces and the size of the attractive basin increases. Hence a larger jump perturbation strategy can be more helpful for LSILS-LK to escape from current local optimum.

We further conduct some experiments for LSILS-LK-3DBP on some medium- and large size TSPs with increased runtime budget. Please see supplementary material for details. The results show that LSILS-LK-3DBP can further improve its performance in case extra budget is available.

In summary, we may conclude that the proposed HC transformation based on local optima can indeed improve the global search ability of existing TSP algorithms. The good performance of LSILS-LK and LSILS-LK-3DBP is due to the HC transformation, which can reduce the number of local optima in the transformed landscape meanwhile preserve the useful information in the optimum of the original TSP.

V Conclusions

In this paper, a new TSP landscape smoothing method called the Homotopic Convex (HC) transformation was proposed. The HC transformation is defined as the convex combination of the original TSP and a convex-hull TSP. For a given TSP, the convex-hull TSP is constructed based on a known optimum. We proved that the convex-hull TSP is unimodal for any $k$-Opt local search ($k\geq 2$). Our empirical study showed that controlled by the convex coefficient, the proposed HC transformation can smooth the landscape in terms that in the transformed TSP, the number of local optima is reduced and the fitness distance correlation of the TSP landscape is increased. Systematic experiments were conducted to verify our statements, in which the HC transformation was used to transform 12 TSP instances with different $\lambda$ values and the characteristics of the landscape of the transformed TSPs were analyzed. In addition, to investigate the relationship between the smoothing effect of the HC transformation and the fitness of the local optimum it based on, we empirically analyzed the effect of the HC transformation on the TSP instance berlin52 based on 89 local optima and the global optimum. The experimental results showed that the HC transformation do have the claimed advantages and using a high-quality local optimum can help the HC transformation achieve a smoothing effect similar to the smoothing effect using the global optimum

Based on the observations, we proposed a landscape smoothing based iterative algorithmic framework, in which the HC transformation is combined with a local search algorithm. An instantiation of the algorithmic framework with 3-Opt local search and double bridge perturbation, named as Landscape Smoothing Iterated Local Search (LSILS), was studied on some widely-used TSP test instances. The experimental results showed that LSILS significantly outperforms ILS and two existing TSP smoothing algorithms on most test instances. In addition, we embedded the LK local search and three double bridge perturbation strategy within the framework, named as LSILS-LK and LSILS-LK-3DBP respectively and tested the resultant algorithm on some middle-size and large-size TSP instances (with up to 20,000 cities). The experimental results showed that LSILS-LK performs the best on all the middle-size and large-size TSP instances against the ILS-LK and other algorithms and the three double bridge perturbation strategy can further improve the performance of LSILS-LK. In addition, we found that LSILS-LK-3DBP performs significantly better than LSILS-LK, which implies that a larger jump strategy is helpful to further improve LSILS. In conclusion, we may conclude that the the proposed HC transformation is very promising to improve the global search ability of existing TSP algorithms.

In the future, we will apply the proposed HC transformation to the asymmetric TSP and other combinatorial optimization problems. It is not easy since the HC transformation is based on a unimodal toy problem with known optimum. For other combinatorial optimization problems, finding a way to build this unimodal toy problem is very challenging. We will try some problems that are similar to the TSP, for example, the Vehicle Routing Problem (VRP). Another possible research avenue is to use the HC transformation to improve other kinds of TSP heuristics, for example, ACO and GA.