IEICE TRANS. INF. & SYST., VOL.E102–D, NO.12 DECEMBER 2019

2329

PAPER

Special Section on Parallel and Distributed Computing and Networking

Accelerating the Held-Karp Algorithm for the Symmetric Traveling
Salesman Problem
Kazuro KIMURA†a) , Shinya HIGA† , Nonmembers, Masao OKITA†† , and Fumihiko INO†† , Members

SUMMARY In this paper, we propose an acceleration method for the
Held-Karp algorithm that solves the symmetric traveling salesman problem by dynamic programming. The proposed method achieves acceleration
with two techniques. First, we locate data-independent subproblems so thatthe subproblems can be solved in parallel. Second, we reduce the number ofsubproblems by a meet in the middle (MITM) technique, which computes
the optimal path from both clockwise and counterclockwise directions. We
show theoretical analysis on the impact of MITM in terms of the time andspace complexities. In experiments, we compared the proposed methodwith a previous method running on a single-core CPU. Experimental results show that the proposed method on an 8-core CPU was 9.5–10.5 times
faster than the previous method on a single-core CPU. Moreover, the proposed method on a graphics processing unit (GPU) was 30–40 times fasterthan that on an 8-core CPU. As a side eﬀect, the proposed method reduced
the memory usage by 48%.
key words: symmetric traveling salesman problem, Held-Karp algorithm,
parallelization, meet in the middle, GPU

1.

Introduction

The traveling salesman problem (TSP) is a problem thatfinds the minimum Hamiltonian cycle for a weighted complete graph G = (V, E), where V and E are the sets of vertices and edges, respectively. A Hamiltonian cycle here is acycle that visits each vertex exactly once and the minimum
Hamiltonian cycle has the lowest-weight cycle, i.e., the lowest sum of edge weights, in the given graph G. Hereafter, wecall the minimum Hamiltonian cycle the optimal path. We
also call the weight of the optimal path the optimal value.
One special case of the TSP is the symmetric TSP(sTSP), which deals with an undirected graph. That is, thegiven graph has the same weight between two vertices ineach opposite direction: w(i, j) = w( j, i) holds for all vertices i, j ∈ V, where w(i, j) represents the weight of theedge (i, j) ∈ E from vertices i to j. Many practical applications include the sTSP as an underlying fundamental
problem. For example, welding robots can minimize theirpower consumption by solving the sTSP [1]. In more detail,an optimal motion plan can be obtained by solving the sTSP
for a graph in which vertices represent the welding points
Manuscript received January 7, 2019.
Manuscript revised March 22, 2019.
Manuscript publicized August 23, 2019.
†
The authors are with the School of Engineering Science,
Osaka University, Toyonaka-shi, 560–8531 Japan.
††
The authors are with the Graduate School of Information
Science and Technology, Osaka University, Suita-shi, 565–0871
Japan.
a) E-mail: u779795i@alumni.osaka-u.ac.jp
DOI: 10.1587/transinf.2019PAP0008

and edge weights represent the costs required for moving
between the welding points.
Algorithms for the sTSP can be classified into three
approaches: exact, approximate, and heuristic approaches.
Exact algorithms in [2], [3] find the optimal path but the
problem size that can be solved in a reasonable time islimited due to high demand for both compute and memory capacities. By contrast, approximate and heuristic algorithms [4], [5] rapidly find an approximate path but the pathis not optimal. Consequently, approximate and heuristic algorithms are typically used if the given graph has many vertices or the approximation error is acceptable for the overlying application. Exact algorithms are preferred otherwise.
The latter case is of our interest where applications require
an exact solution for the sTSP.
The Held-Karp algorithm [2] is an exact algorithm that
deploys dynamic programming for the sTSP. To the best of
our knowledge, the Held-Karp algorithm is the fastest exact
algorithm [6] and its worst-case time complexity is O(n2 2n ),where n is the number of vertices. Therefore, accelerating
the Held-Karp algorithm increases the problem size n thatcan be solved in a reasonable time. The original paper [2]is a theoretical paper that presented the Held-Karp algorithm with its recursive equations and complexity analysis.
By contrast, Kubo [7] presented a real implementation withspecific data and loop structures. The contribution of [7]over [2] are twofold: the paper (1) shows a data structurethat represents a set of vertices with a bit sequence; and (2)finds a loop structure that correctly solves data-dependent
subproblems according to sorted bit sequences.
In this paper, we propose an acceleration method for
the Held-Karp algorithm, which solves the sTSP. The proposed method achieves acceleration with two techniques,parallelization and meet in the middle (MITM). The former
identifies data-independent subproblems from those generated by the Held-Karp algorithm. On the other hand, thelatter approximately halves the number of subproblems by
exploiting the symmetric attribute embedded in the given
graph. Moreover, the proposed method reduces the memory usage as a side eﬀect. We also show theoretical analysis for understanding the impact of MITM in terms of thetime and space complexities. The source code of our implementation is available at http://www-ppl.ist.osaka-u.ac.
jp/research/code/.
The paper is structured as follows. Section 2 introduces
related work to clarify the contribution of the present paper.
Section 3 presents an overview of the Held-Karp algorithmc 2019 The Institute of Electronics, Information and Communication Engineers
Copyright ⃝IEICE TRANS. INF. & SYST., VOL.E102–D, NO.12 DECEMBER 20192330

with [7]. Section 4 describes the proposed method with theoretical analysis. Section 5 shows experiments that evaluate
the performance of the proposed method. Finally, Sect. 6
concludes the paper and discusses future work.
2.

Related Work

Concorde [3] is an exact sTSP solver that relies on a cuttingplane method. This solver is capable of finding the optimal
path for a large graph of n = 85, 900. However, its execution time depends not only on the number n of vertices, i.e.,
the problem size, but also on the weight of edges. In fact,
Ahammed et al. [8] found that the execution time diﬀeredby 30,000 times though n was fixed. On the other hand, theexecution time of the proposed method depends only on thenumber n of vertices. Therefore, the proposed method isuseful for strict real-time situations where the optimal pathmust be computed within a certain period of time.
Moﬀat [9] parallelized the Held-Karp algorithm for thesTSP using a JavaScript application programming interface(API). Similar to the proposed method, this method exploited the data parallelism inherent in the Held-Karp algorithm. However, the maximum problem size was limited
by n ≤ 16 because the method simply stored subproblemsin O(n2n ) space. By contrast, the proposed method saves
the memory consumption by using a bit representation thatis useful for identifying data-independent subproblems with
O(2n ) space. Furthermore, the proposed method realizes further acceleration with MITM.
As for approximation approaches, these approaches return a solution that is theoretically guaranteed to be close
to the optimal solution. Sahni and Gonzalez [10] provedthat there is no algorithm that solves the sTSP in polynomial time. However, a polynomial time algorithm exists forthe sTSP if edge weights in the graph satisfy the triangle
inequality. Christofides [4] presented a 1.5 approximation
algorithm for the sTSP. The worst time complexity of thisalgorithm is O(n3 ), which is better than that of the HeldKarp algorithm. Therefore, approximation approaches areuseful when the overlying application accepts solutions being within a certain bound.
In contrast to the approaches mentioned above, heuristic approaches are usually based on a rule of thumb, which
returns a practical solution but not guaranteed to be optimal
or within certain error bounds. The Lin-Kernighan algorithm [5] deploys an iterative approach that improves an initial solution by removing two or more edges and adding thesame number of edges. The error to the optimal solution wasusually 1–2%, making heuristic approaches more attractive
than approximation approaches [7]. However, Chandra√etal. [11] pointed out that the computed weight could be 4 ntimes larger than the optimal value. Consequently, approximation approaches must be selected if the overlying application requires guarantee on the maximum error.

3. Held-Karp Algorithm
Let V = {0, 1, · · · , n−1} be the set of vertices that compose agraph G. Each vertex then can be identified with an integer
value ranging from 0 to n − 1. The optimal path can beobtained by determining the best visiting order starting froma vertex and ending at the same vertex. In the following,
let vertex n − 1 be the starting/ending vertex. Accordingly,
we determine the visiting order for the remaining vertices,which compose a set V ′ = V \ {n − 1} = {0, 1, · · · , n − 2}.
To describe how the Held-Karp algorithm computes the
optimal value OPT, we define the function tsp as follows.
Definition 1. For S ⊆ V ′ and x ∈ S , let tsp(S , x) be the minimum weight of the path that starts from the starting vertex
n − 1, visits all vertices in set S , and reaches vertex x.
The optimal value OPT and tsp(S , x) then can be given
by Eqs. (1) and (2), respectively [2].
OPT = min′ (tsp(V ′ , x) + w(x, n − 1)),x∈V

tsp(S , x) =
if |S | = 1,w(n − 1, x),min(tsp(S
\
{x},
y)
+
w(y,
x)),
otherwise.
y∈S \{x}

(1)

(2)

According to Eqs. (1) and (2), each value of tsp(S , x), where
S ∈ P(V ′ ) \ ϕ and x ∈ S , is accessed max(n − 1 − |S |, 1)
times to obtain OPT. P(·) here is the power set and ϕ is theempty set. Figure 1 shows how these references occur during OPT computation, where the values of tsp(S , x) are accessed multiple times. This data dependence relation motivates the Held-Karp algorithm to deploy dynamic programming to avoid redundant computation; the algorithm storesthe values of tsp(S , x) in table Dtsp to reuse the computed
values without duplicated computation. This
table can be∑n−1 (n−1)
k k cells, whereindexed by pair (S , x), which have k=1(n )k is the number of k-combinations from n elements, because x can take k diﬀerent vertices when |S | = k, where1 ≤ k ≤ n − 1 (Fig. 1).
Next, the optimal path relies on the following definitions.
Definition 2. Let R = r1 → r2 → · · · → rn be the path that
visits n vertices in the following order: r1 , r2 , · · · , rn .

Fig. 1 Data dependence relation inherent in OPT computation for n = 4.
The value of tsp({0}, 0) is required to compute the values of tsp({0, 1}, 1)
and tsp({0, 2}, 2).

KIMURA et al.: ACCELERATING THE HELD-KARP ALGORITHM FOR THE SYMMETRIC TRAVELING SALESMAN PROBLEM

2331

Definition 3. Consider a path that has the weight oftsp(S , x); the path starts from the starting vertex n − 1, visits
all vertices in set S ⊆ V ′ , and reaches vertex x ∈ S . Letprev(S , x) ∈ V denote the vertex visited immediately before
vertex x on the path.
The vertex prev(S , x) then can be computed according
to Eq. (3) as follows.
prev(S , x) =
if |S | = 1,n − 1,argmin(tsp(S
\
{x},
y)
+
w(y,
x)),
otherwise, y∈S \{x}(3)where the function arg min(·) returns an argument that minimizes the value of the expression given to the function.
According to the above definitions, the optimal pathn − 1 → vn−1 → · · · → v1 → n − 1, where vk represents thek-th vertex visited before reaching the ending vertex n − 1,
can be computed as follows.
()(4)
v1 = arg min tsp(V ′ , x) + w(x, n − 1) ,
x∈V ′

vk =



prev(V ′ , v1 ),prev(V ′ \ {v , v , · · · , v
1

2

if k = 2,k−2 }, vk−1 ), otherwise.

Fig. 2 Data structure of the previous method [7]. Tables are stored in
two-dimensional arrays that have row s and column x. Arcs on the table
represent the execution order for filling up the table.

the subproblem that computes the values of tsp(S , x) andprev(S , x), i.e., fills in tables Dtsp and Dprev .
T(n) can be given byT(n) = T1 (n) + T2 (n),(8)

where T1 (n) is the time complexity for solving a set of all
subproblems, {P(S , x) | S ∈ P(V ′ ) \ ϕ, x ∈ S }, and T2 (n) is
the time complexity for obtaining the optimal value and the
optimal path from tables Dtsp and Dprev , respectively. T1 andT2 can be expressed as follows.

(5)
T1 (n) =n−1∑

k(k − 1)

)
(
n−1
,
k

(9)

That is, the optimal path can be identified from v1 by Eq. (4),then v2 , v3 , · · · , vn−1 by Eq. (5).
Note here that the value of prev(S , x) can be simultaneously computed with that of tsp(S , x). Similar to tsp(S , x),which we mentioned earlier, redundant computation can be
avoided by storing the value of prev(S , x) in table Dprev .

Thus, the worst-case time complexity of the Held-Karp algorithm is O(n2 2n ) because the binomial theorem rewrites
Eq. (9) into T1 (n) = (n − 1)(n − 2)2n−3 .

3.1 Complexity Analysis

3.2

We first analyze L(n), or the space complexity of the HeldKarp algorithm. We assume that a single value occupies
constant space. L(n) then can be given by

The previous method [7], which we denote hereafter by M0 ,deploys a bit string expression to represent set S ⊆ V ′(Fig. 2). That is, set S is expressed by bit string s, or a binary
sequence of n − 1 bits; the x-th bit of s, where 0 ≤ x ≤ n − 2,
is 1 if x ∈ S and 0 otherwise. Consequently, bit string s canbe interpreted as an integer value in the range [0, 2n−1 − 1].
This value can be directly used as an index of tables Dtsp andDprev , which can be stored in two-dimensional arrays. As
shown in Fig. 2, the tables are indexed with row s and column x to store the value for pair (S , x), where S ∈ P(V ′ ) \ ϕ
and x ∈ S . Thus, table values can be rapidly accessed inconstant time.
This indexing scheme has another benefit in obtaining acorrect execution order for solving subproblems. As shownin Fig. 1, some subproblems must be computed before others due to data dependence. Filling up the arrays from lower
to higher indices guarantees such a correct execution order(Fig. 2).

L(n) = Ltsp (n) + Lprev (n),(6)

where Ltsp (n) and Lprev (n) are the space complexities for tables Dtsp and Dprev , respectively. Since these tables have thesame indexing scheme (S , x), where S ∈ P(V ′ )\ϕ and x ∈ S ,we have
)
n−1 (
∑
n−1
.
(7)
Ltsp (n) = Lprev (n) =kk
k=1
Thus, the worst-case space complexity of the Held-Karp algorithm is O(n2n ) because the binomial theorem rewrites
Eq. (7) into Ltsp (n) = Lprev (n) = (n − 1)2n−2 .
We next investigate T(n), or the time complexity of theHeld-Karp algorithm, assuming a comparison as the elementary operation. We define subproblem P(S , x) as follows.
Definition 4. For S ⊆ V ′ and x ∈ S , let P(S , x) denote

k=2

T2 (n) =n − 1.

(10)Previous Method

4. Proposed Method
The Held-Karp algorithm produces not only data-dependent

IEICE TRANS. INF. & SYST., VOL.E102–D, NO.12 DECEMBER 2019

2332

Fig. 3 An example of data-independent subproblems that can be processed in parallel (n = 4). Subproblems {P(S , x)} of the same size |S | canbe processed in parallel, but those with diﬀerent |S | must be processed sequentially (from the top to bottom of the figure).

subproblems but also data-independent subproblems. The
proposed method M1 identifies such data-independent subproblems to process them in parallel. The method M2 further achieves acceleration by MITM, which exploits thesymmetric attribute available in the graph. The method M2
also reduces the memory usage.
4.1 Parallelization Technique
The proposed method M1 classifies subproblems into n − 1groups such that subproblems in each group can be processed in parallel. In other words, subproblems belonging
to diﬀerent groups must be processed sequentially with n−1
steps. The key idea for realizing this classification is togroup subproblems by their size k, where 1 ≤ k ≤ n − 1.

Fig. 4 Proposed data structure that facilitates locating data-independent
subproblems (n = 4). Subproblems are sorted by Hamming weight to store
data-independent subproblems in a contiguous region. Note here that weomit the 0-th row (i = 0) in this figure, because the 0-th row corresponds tothe subproblem of size zero (i.e., an empty set). However, this dummy row
is required to simplify index computation.

We introduce functions bstr and oﬀset to identify thecontiguous rows that correspond to data-independent subproblems.
Definition 7. Let bstr(i), where 0 ≤ i ≤ 2n−1 − 1, be thei-th bit string in a binary sequence of length n − 1 sorted byascending order of Hamming weight where bit sequences ofthe same Hamming weight are further sorted by ascendingorder of numerical value.
Definition 8. Let oﬀset(k) be the number of bit strings withHamming weight at most k − 1 given by
oﬀset(k) =

)
k−1 (
∑
n−1
,
l
l=0

(11)

Definition 5. Let |P(S , x)| = |S | be the size of subproblem
P(S , x).

where 1 ≤ k ≤ n − 1.

More specifically, Eqs. (2) and (3) indicate that subproblems of size k depend only on those of size k−1 (Fig. 3).
This implies that (1) at least n − 1 steps are required to process all 1 ≤ k ≤ n−1 and (2) subproblems of the same size k
are data-independent, allowing them to be solved in parallel.
Because the proposed method adopts the same indexing scheme as the previous method M0 [7], the subproblem
size |S | can be easily computed according to the Hamming
weight of bit string s, where s corresponds to set S ⊆ V ′ .

The number oﬀset(k) implies the minimum index forbit sequences of hamming weight of k. Consequently, sorted
subproblems of size k exist in a contiguous region from bit
strings bstr(oﬀset(k)) to bstr(oﬀset(k + 1) − 1). To locate
these rows rapidly, we precompute the function values ofbstr(i) and oﬀset(k) for all i and k, which are then stored intables Dbstr and Doﬀset , respectively.
We first analyze L′ (n), or the space complexity of M1 ,which can be expressed as follows.

Definition 6. Let d(s) be the number of 1s in bit string s.

L′ (n) = Lbstr (n) + Loﬀset (n) + Ltsp (n) + Lprev (n),Thus, we have |S | = d(s) for any S , which meansthat data-independent subproblems can be classified into thesame group by Hamming weight.
Notice here that (data-independent) subproblems thathave the same Hamming weight involve non-contiguous
rows in each table because the previous method M0 storessubproblems in an increasing order of s (Fig. 2). For example, subproblems of size 1 exist in the first, second, andfourth rows in Fig. 2. Thus, search operations are required
to locate subproblems that form the same group, eliminating the advantage of constant-time access. To deal with thisissue, our method sorts the subproblems, i.e., bit strings,
by Hamming weight so that data-independent subproblems
are stored in contiguous rows in each table (Fig. 4). This
operation incurs an overhead, but once sorted, the dataindependent subproblems can be located in constant time.

where Lbstr (n) and Loﬀset (n) are the space complexities forstoring the values of bstr(i) and oﬀset(k), for all 0 ≤ i ≤2n−1 − 1 and 1 ≤ k ≤ n − 1, respectively. Obviously, we
have Lbstr ∈ O(2n ) and Loﬀset ∈ O(n). Thus, Lbstr and Loﬀset
are asymptotically negligible in comparison to Ltsp , Lprev ∈O(n2n ). Therefore, the space complexity of M1 equals tothat of M0 .
Finally, T′ (n), or the time complexity of M1 , can begiven byT′ (n) = Tbstr (n) + Toﬀset (n) + T1 (n) + T2 (n),

(12)

(13)

where Tbstr (n) and Toﬀset (n) are the time complexities forcomputing the values of bstr(i) and oﬀset(k), for all 0 ≤ i ≤2n−1 − 1 and 1 ≤ k ≤ n − 1, respectively.
• Tbstr (n) ∈ O(n2n ) holds. To compute the values of

KIMURA et al.: ACCELERATING THE HELD-KARP ALGORITHM FOR THE SYMMETRIC TRAVELING SALESMAN PROBLEM

2333

Algorithm 1 Optimal path computation with MITM.
1: Compute set S A and midpoint m by Eq. (16);2: Compute path A by Eq. (5), which starts from prev(S A , m);
3: Compute path rev(B) by Eq. (5), which starts from prev(S A ∪ {m}, m);
4: Compute the optimal path A + rev(rev(B));
Fig. 5 MITM technique for computing the optimal path from both clockwise and counterclockwise directions.

bstr(i), we have to obtain the Hamming weights of all
bit strings. Since the number of bit strings of lengthn − 1 is 2n−1 , the worst-case time complexity for obtaining all Hamming weights is given by O(n2n ). As
for sorting subproblems, the merge sort algorithm has
the worst-case time complexity of O(n2n ).
• Toﬀset (n) ∈ O(n2n ) holds. The values of oﬀset(k) can
be computed by traversing row indices from the top oftable Dbstr ; oﬀset(k) is the first index that has Hamming
weight k.
Thus, the terms in Eq. (13) except T1 are asymptoticallynegligible in comparison to T1 . Since T1 is the time complexity for solving all subproblems, parallelizing this performance bottleneck leads to acceleration.
4.2 Meet in the Middle (MITM) Technique
The proposed method M2 integrates the MITM techniqueinto M1 . Given an undirected graph, MITM computes the
optimal path from both clockwise and counterclockwise directions, solving subproblems of at most size ⌈n/2⌉ insteadof those of at most n − 1 (Fig. 5). In comparison to methods
M0 and M1 , MITM roughly halves the number of subproblems that are to be solved.
We first explain how MITM obtains the optimal value.
To do this, we introduce the following operators for paths.
Definition 9. Let + be the binary operator that concatenates
two paths. That is, we obtain path R + Q = r → · · · → x →· · · → q from paths R = r → · · · → x and Q = x → · · · → q.
Definition 10. Let rev be the unary operator that returns thegiven path in reverse order. That is, we obtain path rev(R) =rn → rn−1 → · · · → r1 from path R = r1 → r2 → · · · → rn .
Theorem 1. Obviously, rev(rev(R)) = R holds.
Let m be the midpoint of the optimal path (Fig. 5); since
the optimal path contains n vertices, m is the ⌈n/2⌉-th vertex
from the starting vertex. Consequently, the optimal path canbe given by A+B, where A = n−1 → a1 → · · · → a⌈n/2⌉−1 →m and B = m → b1 → · · · → bn−⌈n/2⌉−1 → n − 1. Theorem 2
holds for the weights of paths A and B.
Theorem 2. Let S A be the set of vertices that constitute path
A. The weights of paths A and B are then given by tsp(S A , m)and tsp(S A ∪ {m}, m), respectively, where S A is the complementary set of S A in V ′ .
Proof. According to Definition 1, the weights of paths A

and rev(B) are tsp(S A , m) and tsp(S A ∪ {m}, m), respectively.
Since the graph consists of undirected edges, i.e., weights
are symmetric, the weight of path B equals to that of pathrev(B). Therefore, the weight of path B is given by tsp(S A ∪{m}, m).
Equation (14) gives the optimal value OPT.
()
OPT = min tsp(S , x) + tsp(S ∪ {x}, x) ,(S ,x)∈F(14)

where set F is given by
F = {(S , x) | S ∈ P(V ′ ) \ ϕ, x ∈ S such that |S | = ⌈n/2⌉}.
(15)
Once we set a pair (S , x) for the first term of Eq. (14), we
accordingly obtain a unique pair (S ∪ {x}, x) for the second
term. Furthermore, |S ∪ {x}| ≤ ⌈n/2⌉ holds. Therefore, we
can compute OPT from subproblems of at most size ⌈n/2⌉.
We next describe how the proposed method M2 obtains
the optimal path with MITM. Algorithm 1 returns the optimal path from set S A and midpoint m, which can be given
by()(S A , m) = arg min tsp(S , x) + tsp(S ∪ {x}, x) .
(16)
(S ,x)∈F

We next analyze T′′ (n), or the time complexity of M2 ,which is given by
T′′ (n) = Tbstr (n) + Toﬀset (n) + T′1 (n) + T′2 (n),(17)

where T′1 (n) is the time complexity for solving subproblemsof size k ≤ ⌈n/2⌉ and T′2 (n) is that for obtaining the optimalvalue and path. The former term is given by
T′1 (n) =

⌈n/2⌉
∑
k=2

k(k − 1)

)
(
n−1
.
k

(18)

Lemmas 1 and 2 then give Theorem 3 whose proofs
are presented in Appendix A, Appendix B and Appendix C,
respectively.
Lemma 1. The time complexity T′1 (n) for solving all subproblems of size k ≤ ⌈n/2⌉ is given by
T′1 (n)
T (n)/2,) if n is even,( 1=(n − 1)(n − 2) n − 3, otherwise.
T1 (n)/2 +(n − 3)/2
2(19)

IEICE TRANS. INF. & SYST., VOL.E102–D, NO.12 DECEMBER 2019

2334

(n)
√
Lemma 2. For any even number n, n/2∈ Θ(2n / n) holds.
Theorem 3. Given the time complexity T1 (n) for solving
subproblems for all 1 ≤ k ≤ n − 1 and that T′1 (n) for solving
subproblems for all k ≤ ⌈n/2⌉, we havelim T′1 (n)/T1 (n) = 1/2.

n→∞

(20)

As for the latter term of Eq. (17), we have
)
(
n−1
− 1,
(21)
T′2 (n) = ⌈n/2⌉⌈n/2⌉( n−1 )
because there are ⌈n/2⌉
sets whose size is |S | = ⌈n/2⌉. Ac√
cording to Lemma 2, T′2 (n) ∈ O( n2n ) holds. Therefore,
similar to Tbstr (n) and Toﬀset (n), T′2 (n) is asymptotically negligible in comparison to T1 in Eq. (17). In summary, MITM
is expected to double the performance according to Theorem 3.
Note that MITM can be integrated into the previous
method M0 instead of M1 ; this integration allows us to skip
subproblems of larger than ⌈n/2⌉. However, we avoid thisbecause MITM is able to reduce the memory usage if subproblems in tables Dtsp and Dprev are sorted by Hamming
weight. Thus, there exists a common procedure in methods
M1 and M2 .
4.3 Reducing Memory Usage
The proposed method M2 is allowed to store the solutions
for the subproblems of size k ≤ ⌈n/2⌉ whereas the previous
method M0 has to store the solutions for those of size k ≤n − 1. Therefore, the memory usage of M2 asymptoticallyreduces to 1/2 in comparison to that of M0 as follows.
Theorem 4. Let L′tsp (n) and L′prev (n) be the space complexities of tables Dtsp and Dprev in method M2 , respectively. We
then havelim L′tsp (n)/Ltsp (n) = 1/2,

(22)

lim L′prev (n)/Lprev (n) = 1/2.

(23)

n→∞
n→∞

The proof for Theorem 4 is presented in Appendix D.
The above discussion assumes that necessary rows intables Dtsp and Dprev exist in contiguous regions as shownin Fig. 6. In other words, tables in the previous method M0
fails to save the memory usage because necessary rows exist
in non-contiguous regions; the entire table must be allocatedfor execution.
Similar to the original method M0 , the tables mustbe accessed in constant time to achieve acceleration with
MITM. To realize this constant time access, we introducefunction rank as follows.
Definition 11. Consider the i-th bit string s, where 0 ≤ i ≤2n−1 −1 and s = bstr(i), in a sequence of bit strings sorted byascending order of Hamming weight where bit sequences ofthe same Hamming weight are further sorted by ascending

Fig. 6 Comparison of data structures deployed for table Dprev of the previous method M0 and that of the proposed method M2 (n = 6). Method M2
stores necessary rows in a contiguous region, allowing unnecessary rows
not to be allocated. Unnecessary rows are colored with light blue.

order of numerical value. Let rank(s) be the rank of string s
in the sequence of bit strings that have the same Hamming
weight. rank(s) is then given by rank(s) = i − oﬀset(d(s)).
Using the function rank, the value for (s, x) at the k-th
step can be stored in the element at (oﬀset(k) + rank(s), x) of
the two-dimensional array. As shown in Fig. 6, this indexing
scheme allows necessary rows to be stored in a contiguous
region and to be accessed in constant time. For example,
a bit sequence of 11010 has a hamming weight of 3. Accordingly, this sequence corresponds to the subproblems ofsize 3, which are processed at the third step (k = 3). We
then have oﬀset(3) = 16 according to Definition 8. Furthermore, we have rank(11010) = 8 because rank(11010) gives
the rank of 11010 in sorted sequences that have a hamming
weight of 3: 00111, 01011, 01101, 01110, 10011, 10101,
10110, 11001, 11010 and 11100, each corresponding to 7,
11, 13, 14, 19, 21, 22, 25, 26 and 28 in decimal.
The time complexity Trank (n) for computing rank(s) is
negligible in comparison to T′1 . In more detail, Trank (n) ∈O(2n ) holds because the values of rank(s) can be computedby traversing row indices from the top of table Dbstr . Similarly, the space complexity Lrank (n) for storing rank(s) is
negligible in comparison to Ltsp and Lprev ; Lrank (n) ∈ O(2n )holds if we store all values in table Drank .
The space complexity for table Dtsp can be further reduced by allocating rows only for two consecutive steps andoverwriting the allocated rows that are not to be accessed.
More specifically, tsp(S , x) at the (k + 1)-th step depends
only on that at the k-th step, which means that all valuesof tsp(S , x) except those at the last step can be overwritten. Accordingly, the proposed method M2 allocates only
the rows for the latest two steps rather than all rows in table
Dtsp . Furthermore, the method deploys a double buﬀering
scheme that stores the rows computed at odd and even steps
in tables Dtsp1 and Dtsp0 , respectively (Fig. 7). Each tableholds the value for bit string s in the rank(s)-th row to realize constant time access to a contiguous region.
Finally, we analyze the space complexities of tables
Dtsp0 and Dtsp1 . The number of values computed at the kth step is maximized when k = ⌈(n − 1)/2⌉. Each table

KIMURA et al.: ACCELERATING THE HELD-KARP ALGORITHM FOR THE SYMMETRIC TRAVELING SALESMAN PROBLEM

2335
Table 1

Specification of experimental machine.

Classification
CPU
Memory capacity (CPU)GPU
Memory capacity (GPU)GPU driverOS
Development environment
Compile option

Specification
Intel Xeon Silver 4110 (8 cores)
96 GB
NVIDIA Tesla V100 PCIe
16 GB
396.44
Ubuntu 16.045 LTS
GCC 5.4.0 and CUDA 9.2
O3 optimization

Fig. 7 Double buﬀering scheme that alternatively uses tables Dtsp0 and
Dtsp1 (n = 8). Instead of allocating all rows in Dtsp , limited rows are allocated only for two consecutive steps.

( n−1 )
then requires ⌈(n−1)/2⌉
rows and n − 1 columns to store allvalues generated at the k-th step. According to Lemma 2,
both the space √complexities for tables Dtsp0 and Dtsp1 canbe given by O( n2n ), which is negligible in comparison to
Lprev ∈ O(n2n ).
5.

Experiments

We evaluated the proposed methods M1 and M2 by comparing them with the previous method M0 . The evaluation
metrics were the execution time, the speedup ratio, and thememory usage, which demonstrated the impact of parallelization and MITM.
For the experiments, all methods were implemented
with the C language and OpenMP directives [12] to obtain performance results on a multi-core CPU. We adopted
bucket sort to implement sorting operations on the CPU.
Similarly, we implemented all methods except M0 with thecompute unified device architecture (CUDA) [13] to executethe methods on a GPU. Sort operations were implemented
with the Thrust library [14]. Table 1 shows the specification of the deployed machine. Hyper-threading was enabledon the 8-core CPU, and thus, at most 16 threads were executed on 16 logical cores. Our CPU-based implementation parallelizes the k-th step with all CPU threads. On theother hand, the(GPU-based
implementation parallelizes the)k-th step with n−1
threads;
each thread is responsible forkprocessing a bit sequence and each block consists of 1024
threads.
As for the problem size n, we varied the number ofvertices from 4 to 29, which was the maximum number that
could be solved on the GPU successfully; the execution for
n ≥ 30 failed due to exhaustion of GPU memory. Given acomplete weighted graph, all the methods returned the optimal path with its value OPT. Each edge in the graph had
a random weight specified as an integer value in the range[1, 1023]. For each size of n, we generated 10 graphs that
had random weights. For each graph, we executed 50 times(500 times for each n) and the minimum time was selected
for each n. We found that the maximum time was close to
the minimum time. For example, the execution times of M2
on the GPU ranged from 1.16 s to 1.17 s for n = 29. Thus,
the gap between the minimum and maximum times was negligible on our machine because the time complexity of theFig. 8 Comparison of the previous method M0 and the proposed methods M1 and M2 in terms of execution time.

Held-Karp algorithm depends only on n.
Figure 8 shows the execution time for diﬀerent problem sizes (4 ≤ n ≤ 29). The proposed methods M1 andM2 were faster than the previous method M0 when the problem size was suﬃciently large (n ≥ 13) where eﬃcient parallelization was achieved. In other words, small problem
sizes solved within 0.001 s could not be accelerated by parallelization, which involved a synchronization overhead; asynchronization operation is necessary at every step to increase the problem size k. We also found that M2 was fasterthan M1 because acceleration by MITM is independent fromthat by parallelization. That is, the former technique reduces
the number of subproblems that are to be solved whereas thelatter technique reduces the execution time for solving subproblems. In particular, M2 running on the GPU was 30–40
times faster than that on the CPU.
All of execution times in Fig. 8 include data transfer
time. The data transfer time of M2 on the GPU for n = 29
was 65.2 µs, which was negligible compared to the execution time of 1.16 s. Note that data transfer occurs only at
the program initialization and finalization to send the input
to the GPU and receive the output from the GPU, respectively. The input size is given by O(n2 ), which corresponds
to an adjacency matrix that stores graph information. On
the other hand, the output size is given by O(n), which corresponds to the optimal path and value. Thus, the amountof data transfer is asymptotically negligible to that of computation, which has the time complexity of O(n2 2n ). Notethat our GPU-based implementation failed to overcome theCPU-based implementation when n < 16, because the problem size was not suﬃcient large to maximize the eﬃciency
of GPU execution. The data transfer time for n = 4 occupied
12.6 % of the execution time.
Figure 9 shows the speedup ratio obtained on the CPU.

IEICE TRANS. INF. & SYST., VOL.E102–D, NO.12 DECEMBER 2019

2336

Fig. 9 Speedup ratio of the proposed methods M1 and M2 over the previous method M0 with diﬀerent problem sizes. All methods were executedon the CPU.

Fig. 10 Speedup ratio of the proposed method M1 over the previous
method M0 with diﬀerent numbers of CPU threads (n = 29). Executions
using more than 8 threads involved multiple threads to be executed on aphysical CPU core.

Our parallelization technique M1 achieved a 6.3× speedup
for large problems of n ≥ 25. More speedups were achievedby M2 , which integrated the MITM technique into M1 ; weobtained a 9.5× speedup when n was odd and a 10.5×
speedup otherwise. Thus, acceleration achieved by MITM
reached a factor of 1.5 and that of 1.7 when n was odd andeven, respectively. This zigzag behavior observed in Fig. 9,i.e., higher speedups for even numbers of vertices, perfectly
matches to our theoretical analysis presented in Sect. 4.2.
We next investigated how our parallelization technique
M1 increased the speedup with the number of CPU threads(Fig. 10). As shown in this figure, we observed a linear
speedup for n = 29, which increased in proportion to thenumber of threads. When n = 29, there were many subproblems suﬃcient to make parallel threads( busy.
) For exn−1
ample, the proposed method generated ⌈n/2⌉ ⌈n/2⌉
subproblems at the ⌈n/2⌉-th step, so that 561, 632, 400 subproblems
were assigned to 16 CPU threads when n = 29. In Fig. 10,
the speedup slightly decreased when we increased the number of threads from 8 to 9. This slight decrease was due
to hyper-threading, which dropped the eﬃciency of thread
execution when multiple threads were assigned to a single
physical core. A similar behavior can be found in othermemory-intensive applications [15]–[17].
Finally, we investigated the memory usage of the proposed and previous methods. The memory usage was computed assuming that tables Dprev , Dtsp , Dbstr , and Drank have
1 B, 2 B, 4 B, and 4 B values, respectively. Note that we

Fig. 11 Reduction rate γ of memory usage. A positive rate means thatthe proposed method M2 reduced the memory usage of the previous method
M0 . Execution of n = 30 failed due to exhaustion of GPU memory.

avoided considering variables that can be stored in polynomial space; In fact, these variables consume several KB
when n = 29, which is negligible to other data. Figure 11
shows the reduction rate of memory usage, γ = (α − β)/α,where α and β are the memory usage of M0 and that ofM2 , respectively. We observed that the rate γ had a negative value when n ≤ 7 but increased with n. The reason
for this negative value was due to the additional tables Dbstr
and Doﬀset . On the other hand, the increase was yielded by
the memory reduction technique presented in Sect. 4.3; theproposed method replaces table Dtsp with two smaller tables
Dtsp0 and Dtsp1√, which reduces the space complexity fromO(n2n ) to O( n2n ). Consequently, the gap from M0 to M2
increased with n.
Note here that the memory usage showed a zigzag behavior similar to that observed for the speedup (Fig. 9). Thisbehavior was also theoretically analyzed in Sect. 4.2. The
highest reduction rate of 48% was obtained when n = 28,
which was the largest even problem size processed successfully on the deployed GPU. In theoretical, the memory usage roughly doubles as the problem size increases by 1. Actually, method M2 consumed 5.2 GB and 11.2 GB of GPU
memory when n = 28 and n = 29, respectively. Thus, theproblem size of n = 30 required approximately 22 GB of
GPU memory for execution, which was beyond the capacity
available on the deployed GPU.
6. Conclusion
In this paper, we proposed an acceleration method for the
Held-Karp algorithm, which is an exact algorithm for thesTSP. The proposed method achieves acceleration with parallelization and MITM. We also showed theoretical analysis
for the time and space complexities of the proposed method.
Experimental results show that the proposed method
on an 8-core CPU was 9.5–10.5 times faster than the previous method on a single-core CPU. Moreover, the proposed method on the latest GPU was 30–40 times faster thanthat on the multi-core CPU. As a side eﬀect, the proposed
method reduced the memory usage by 48%, which matches
to the theoretical analysis.
Future work is to deal with large problem sizes thatcannot be naively stored in the GPU memory due to memory

KIMURA et al.: ACCELERATING THE HELD-KARP ALGORITHM FOR THE SYMMETRIC TRAVELING SALESMAN PROBLEM

2337

exhaustion. An out-of-core processing scheme [18], [19]may be useful for realizing this large-scale computation; thescheme divides data into small pieces and iteratively processes the pieces with overlapping GPU computation with
CPU–GPU data transfer.
Acknowledgments
This study was supported in part by the Japan Society
for the Promotion of Science KAKENHI Grant Numbers
JP15H01687 and JP16H02801. We are also grateful to theanonymous reviewers for their valuable comments.

monte carlo framework,” Proc. 5th Int’l Joint Conf. Biomedical Engineering Systems and Technologies (BIOSTEC’12), pp.222–235,
Feb. 2012.
[18] Y. Lu, F. Ino, and K. Hagihara, “Cache-aware GPU optimizationfor out-of-core cone beam CT reconstruction of high-resolution volumes,” IEICE Trans. Inf. & Syst., vol.E99-D, no.12, pp.3060–3071,
Dec. 2016.
[19] J. Shen, K. Shigeoka, F. Ino, and K. Hagihara, “GPU-based branch
and bound method to solve large 0-1 knapsack problems with datacentric strategies,” Concurrency and Computation: Practice and Experience, vol.31, no.4, e4954, Feb. 2019.
[20] J. Wallis, Arithmetica Infinitorum, Oxford, 1656.

Appendix A:

Proof for Lemma 1

References[1] X. Wang, B. Tang, Y. Yan, and X. Gu, “Time-optimal path planning for dual-welding robots based on intelligent optimization strategy,” Trans. Intelligent Welding Manufacturing, vol.1, no.2, pp.47–
59, Dec. 2017.
[2] M. Held and R.M. Karp, “A dynamic programming approach to sequencing problems,” J. Society for Industrial and Applied Mathematics, vol.10, no.1, pp.196–210, March 1962.
[3] W. Cook, “Concorde home,” 2016. http://www.math.uwaterloo.ca/
tsp/concorde.html.
[4] N. Christofides, “Worst-case analysis of a new heuristic for the travelling salesman problem,” Tech. Rep., Carnegie-Mellon University
Pittsburgh Pa Management Sciences Research Group, Feb. 1976.
[5] S. Lin and B.W. Kernighan, “An eﬀective heuristic algorithm forthe traveling-salesman problem,” Operations Research, vol.21, no.2,
pp.498–516, April 1973.
[6] W.J. Cook, In Pursuit of the Traveling Salesman: Mathematics at theLimits of Computation, Princeton University Press, Dec. 2011.
[7] M. Kubo, “An invitation to the traveling salesman problem ii,” Operations Research, vol.39, no.2, pp.91–96, Feb. 1994. (in Japanese).
[8] F. Ahammed and P. Moscato, “Evolving l-systems as an intelligent
design approach to find classes of diﬃcult-to-solve traveling salesman problem instances,” Proc. Int’l Conf. Applications of Evolutionary Computation (EvoApplications’11), Part I, pp.1–11, April
2011.
[9] A. Moﬀat, “Solving the TSP with WebGL and gpgpu.js,” 2017.
https://amoﬀat.github.io/held-karp-gpu-demo/.
[10] S. Sahni and T. Gonzalez, “P-complete approximation problems,” J.
ACM, vol.23, no.3, pp.555–565, July 1976.
[11] B. Chandra, H. Karloﬀ, and C. Tovey, “New results on the old kopt algorithm for the tsp,” Proc. 5th ACM-SIAM Symp. Discrete
Algorithms (SODA’94), pp.150–159, Jan. 1994.
[12] R. Chandra, L. Dagum, D. Kohr, D. Maydan, J. McDonald, and R.
Menon, Parallel Programming in OpenMP, Morgan Kaufmann, San
Mateo, CA, Oct. 2000.
[13] NVIDIA Corporation, “CUDA C Programming Guide Version 9.2,”
Aug. 2018.
[14] N. Bell and J. Hoberock, Thrust: A Productivity-Oriented Library
for CUDA, ch. 26, Morgan Kaufmann, San Mateo, CA, Jan. 2011.
http://thrust.github.io/.
[15] Y. Matsuzaki, N. Uchikoga, M. Ohue, T. Shimoda, T. Sato, T. Ishida,
and Y. Akiyama, “Megadock 3.0: a high-performance proteinprotein interaction prediction software using hybrid parallel computing for petascale supercomputing environments,” Source Codefor Biology and Medicine, vol.8, no.1, Sept. 2013.
[16] C. Couder-Castañeda, J.C. Ortiz-Alemán, M.G.O. del Castillo, and
M. Nava-Flores, “Forward modeling of gravitational fields on hybrid multi-threaded cluster,” Geofı́sica Internacional, vol.54, no.1,
pp.31–48, Jan. 2015.
[17] L.D. Antonov, C. Andreetta, and T. Hamelryck, “Parallel GPGPU
evaluation of small angle X-ray scattering profiles in a markov chain

Let N0 be a set of non-negative integers. We introduce Lemmas 3 and 4 to show a proof for Lemma 1.
Lemma 3. For any n ∈ N0 and k ∈ N0 , where k ≤ n,
)
(
( )
n−1
n
.
(A· 1)=nk
k−1
kProof. Lemma 3 obviously holds( ) according to the definitionof the binomial coeﬃcient: nk = n!/(k!(n − k)!).
Lemma 4. For any n ∈ N0 ,)n (∑2n + 1= 22n ,
k
k=0
( )
n ( )
∑
1 2n
2n
.
= 22n−1 +
k
2 n
k=0

(A· 2)(A· 3)

Proof. First, we show a proof for Eq. (A· 2). According to
the symmetry rule for binomial coeﬃcients, we have)) ((2n + 12n + 1,(A· 4)=
2n − k + 1
kfor all 0 ≤ k ≤ 2n + 1. Therefore, we have
) 2n+1n (∑ (2n + 1)∑
2n + 1
.
=
k
k
k=n+1
k=0

(A· 5)

Using this equivalence, the left-hand side of Eq. (A· 2) canbe rewritten as follows.
 n (
)
)
)
n (
n (
∑


1
2n + 1∑ 2n + 1 ∑ 2n + 1 
+
= 





k
k
k
2
k=0
k=0
k=0

 n (
) 2n+1∑ (2n + 1)
1
∑ 2n + 1
= 
+




k
k
2  k=0
k=n+1
)
(
2n+1
1 ∑ 2n + 1=(A· 6).
2 k=0
k
Applying the binomial theorem to Eq. (A· 6) gives))
2n+1 (n (∑1 ∑ 2n + 1
2n + 1= 22n .
=
2
k
k
k=0
k=0

(A· 7)

IEICE TRANS. INF. & SYST., VOL.E102–D, NO.12 DECEMBER 2019

2338

Thus, Eq. (A· 2) holds.
We next show a proof for Eq. (A· 3). Similar toEq. (A· 4), the symmetry rule for binomial coeﬃcients gives
)
( ) (
2n
2n
,
=
2n − kkfor all 0 ≤ k ≤ 2n. Therefore, we have2n ( )n−1 ( )
∑
∑
2n
2n
.
=
k
k
k=n+1
k=0

Therefore,
T′1 (n) = (n − 1)(n − 2)

= (n − 1)(n − 2)2= T1 (n)/2.

(A· 8)(A· 9)

Using this equivalence, the left-hand side of Eq. (A· 3) canbe rewritten as follows.
 n−1 ( ) ( ) n−1 ( ) ( )
n ( )


∑
2n 
2n ∑ 2n
1
2n

∑ 2n
+
+
+
= 





n
k
n
k
2
k
k=0
k=0
k=0
 n−1 ( ) ( )
)
)
(
(
2n
∑


1
2n 
2n
2n

∑ 2n
= 
+
+
+




n 
k
n
2 k=0 kk=n+1
 2n ( ) ( )


1
2n 

∑ 2n
= 
.
+
(A· 10)



n 
k
2
k=0

l=0
n−4

(A· 15)

Case 2. n is odd. Since n − 3 is an even number, there is an
integer m such that 2m = n − 3. Therefore, ⌈n/2⌉ − 2 =m. Using these equations, we have
)
) ∑
m (
2m
n−3
.
=
l
l
l=0

⌈n/2⌉−2∑ (l=0(A· 16)

Applying Lemma 4 to Eq. (A· 16), we obtain
)
( )
1 2m
n−3
2m−1
=2
+
l
2 m
)
(
1 n−3
n−4
.
=2 +
2 (n − 3)/2

⌈n/2⌉−2∑ (l=0(A· 17)

Therefore,
)
n−3
l
l=0
)
(
(n − 1)(n − 2) n − 3.
= T1 (n)/2 +(n − 3)/2
2(A· 18)

Applying the binomial theorem to Eq. (A· 10) gives
 2n ( ) ( )
n ( )


∑
2n 
1
2n

∑ 2n
+
= 





n
k
k
2
k=0
k=0
( )
1 2n
= 22n−1 +.
(A· 11)
2 nThus, Eq. (A· 3) holds.

)
n−3
l

⌈n/2⌉−2∑ (T′1 (n) = (n − 1)(n − 2)⌈n/2⌉−2∑ (Thus, by the two cases above, Lemma 1 holds.

Using Lemmas 3 and 4, we show a proof for Lemma 1.
Proof. By applying Lemma 3 twice to Eq. (18), we obtain
)
(
⌈n/2⌉
∑
n−1
T′1 (n) =k(k − 1)
k
k=2
⌈n/2⌉
∑ (n − 3)= (n − 1)(n − 2)k−2
k=2⌈n/2⌉−2∑ (n − 3)= (n − 1)(n − 2).
(A· 12)ll=0Appendix B:

Proof for Lemma 2

Proof. According to the Wallis product [20] for π, we have√
22n( ) = π.
n→∞ √ 2nnnlim

That is, using the (ε, δ)-definition of limit,∀ε > 0, ∃n0 ∈ N0 s.t. ∀n ∈ N0 ,We then consider two cases.
Case 1. n is even. Since n − 3 is an odd number, there isan integer m such that 2m + 1 = n − 3. Therefore,⌈n/2⌉ − 2 = m. Using these equations, we have)⌈n/2⌉−2m (∑ (n − 3) ∑2m + 1.
(A· 13)
=
l
l
l=0
l=0
Applying Lemma 4 to Eq. (A· 13), we obtain⌈n/2⌉−2∑ (n − 3)= 22m = 2n−4 .
ll=0

n > n0 ⇒√22n( ) − π < ε,
√ 2n
n
n

(A· 20)which can be rewritten as√∃ε0 s.t. 0 < ε0 < π, ∃n0 ∈ N0 s.t. ∀n ∈ N0 ,n > n0 ⇒(A· 14)(A· 19)√22n( ) − π < ε0 .
√ 2n
n
n

(A· 21)

KIMURA et al.: ACCELERATING THE HELD-KARP ALGORITHM FOR THE SYMMETRIC TRAVELING SALESMAN PROBLEM

2339

)
(
(n − 1)(n − 2) n − 3∈ Θ(n1.5 2n ).
(n − 3)/2
2

Then
√
22n
( ) − π < ε0
√ 2n
n
n
√
22n
( ) − π < ε0
⇒ − ε0 <
√ 2n
n
n
√
√
22n
( ) < π + ε0 .
⇒ π − ε0 <
√ 2n
n
n
√
√
Since 0 < ε0 < π gives π − ε0 > 0,
( )
√ 2n
n
n
1
1
< √
<
√
22n
π + ε0π − ε0
( )
1
22n
22n
1
2n
< √
⇒√
√ <
√ .
n
π + ε0 nπ − ε0 n√√
Using c1 = 1/( π + ε0 ) and c2 = 1/( π − ε0 ),(A· 28)Since T1 (n) = (n − 1)(n − 2)2n−3 ∈ Θ(n2 2n ), we obtainlim T′1 (n)/T1 (n) = 1/2.

n→∞(A· 29)Thus, by the two cases above, Theorem 3 holds.
Appendix D:(A· 22)

Proof of Theorem 4

Proof. We show a proof for Eq. (22). As for Eq. (23), weomit the proof but it can be presented in a similar manner.
First, L′tsp (n) can be expressed as follows.
L′tsp (n) =

)
(
n−1
.
k
k
k=1

⌈n/2⌉
∑

(A· 30)By applying Lemma 3 to Eq. (A· 30), we obtain(A· 23)∃n0 ∈ N0 s.t. ∀n ∈ N0 ,

( )
22n
22n
2n
< c2 √ . (A· 24)
n > n0 ⇒ c1 √ <nnn( 2n )( )22n∈ Θ √ . Assuming n be an even number, we
Thus,
n
n
obtain
( n )
( )
2
n
∈Θ √ .
(A· 25)
n/2
n

)
n−2
k−1
k=1
⌈n/2⌉−1∑ (n − 2).
= (n − 1)ll=0L′tsp (n) = (n − 1)

⌈n/2⌉
∑(

(A· 31)We then consider two cases.
Case 1. n is even. Since n − 2 is an even number, there is an
integer m such that 2m = n − 2. Therefore, ⌈n/2⌉ − 1 =m. Using these equations, we have
)
) ∑
m (
2m
n−2
.
=
l
l
l=0

⌈n/2⌉−1∑ (l=0(A· 32)

Applying Lemma 4 to Eq. (A· 32), we obtain
)
( )
1 2m
n−2
= 22m−1 +
l
2 m
)
(
1 n−2
.
= 2n−3 +2 (n − 2)/2

⌈n/2⌉−1∑ (Appendix C:

Proof of Theorem 3
l=0

Proof. We consider two cases.
Case 1. n is even. According to Lemma 1, Eq. (20) obviously holds.
Case 2. n is odd. According to Lemma 1,()(n − 1)(n − 2) n − 3(n − 3)/2
2.
T′1 (n)/T1 (n) = 1/2 +T1 (n)(A· 26)
According to Lemma 2, given an even number n − 3,
we have
)
( n )
(
2
n−3
∈Θ √ .
(A· 27)(n − 3)/2nThus,(A· 33)

Therefore,
)
n−2
l
l=0
)
(
n−1 n−2= (n − 1)2n−3 +2 (n − 2)/2)(
n−1 n−2.
= Ltsp (n)/2 +(A· 34)2 (n − 2)/2

L′tsp (n) = (n − 1)⌈n/2⌉−1∑ (According to Lemma 2, given an even number n − 2,
we have
)
( n )
(
2
n−2
∈Θ √ .
(A· 35)(n − 2)/2n

IEICE TRANS. INF. & SYST., VOL.E102–D, NO.12 DECEMBER 2019

2340

Since Ltsp (n) = (n − 1)2n−2 ∈ Θ(n2n ), we obtain

Thus,
)
(
√
n−1 n−2∈ Θ( n2n ).
2 (n − 2)/2

(A· 36)Since Ltsp (n) = (n − 1)2n−2 ∈ Θ(n2n ), we obtainlim L′tsp (n)/Ltsp (n))(

n−1 n−2 











2 (n − 2)/2 


1/2
+
= 1/2. (A· 37)= lim 





n→∞ 
L
(n)

tsp







lim L′tsp (n)/Ltsp (n)
)
(

n−3 






2(n
−
2)





(n − 3)/2 = 1/2.
1/2+= lim 




n→∞ 
L
(n)


tsp






n→∞

n→∞

Case 2. n is odd. Since n − 2 is an odd number, there is an
integer m such that 2m + 1 = n − 2. Therefore, ⌈n/2⌉ −1 = m + 1. Using these equations, we have
)
) ∑
m+1 (
2m + 1
n−2
=
l
l
l=0
(
m
∑ 2m + 1) (2m + 1)
.
+
=
m+1
l
l=0

⌈n/2⌉−1∑ (l=0(A· 44)Thus, by the two cases above, Lemma 4 holds.

Kazuro Kimura
received the B.E. degree
information and computer sciences from Osaka
University, Osaka, Japan, in 2019. He is currently an engineer at Forever Co., Ltd. His current research interests include computer graphics and high performance computing.

(A· 38)
Applying Lemma 4 to Eq. (A· 38), we obtain))(
2m + 1n−2= 22m +
m+1
l
)
(
n−2
.
= 2n−3 +(n − 1)/2⌈n/2⌉−1∑ (l=0

According to Lemma 3,)()(
2(n − 2) n − 3
n−2.
=
n − 1 (n − 3)/2
(n − 1)/2(A· 39)(A· 40)

Therefore,
L′tsp (n) = (n − 1)

)
n−2
l

⌈n/2⌉−1
∑ (
l=0

)}
(
{
2(n − 2) n − 3= (n − 1) 2n−3 +n − 1 (n − 3)/2)(
n−3
= (n − 1)2n−3 + 2(n − 2)(n − 3)/2)(
n−3. (A· 41)
= Ltsp (n)/2 + 2(n − 2)(n − 3)/2
According to Lemma 2, given an even number n − 3,
we have
)
( n )
(
2
n−3
∈Θ √ .
(A· 42)(n − 3)/2
n
Thus,
)
√
n−3
∈ Θ( n2n ).
2(n − 2)(n − 3)/2

Shinya Higa
received the B.E. degree information and computer sciences from Osaka
University, Osaka, Japan, in 2019. He is currently working toward the M.E. degree in information and computer sciences at Nagoya University. His current research interests include
computer graphics and high performance computing.

((A· 43)

Masao Okita
is an Associate Professor in
the Department of Computer Science, Graduate
School of Information Science and Technology,
Osaka University since 2009. He received the
B.E. and M.E., and Ph.D. degrees in information and computer sciences from Osaka University, Osaka, Japan, in 2001, 2003, and 2006, respectively. From 2006 to 2009, he worked forAcces Co., Ltd., as a section chief. His researchinterests include parallel processing, especially
physiological simulation on supercomputers andGPU-accelerated machine learning.

Fumihiko Ino
received the B.E., M.E.,
and Ph.D. degrees in information and computersciences from Osaka University, Osaka, Japan,
in 1998, 2000, and 2004, respectively. He iscurrently a Professor in the Graduate School ofInformation Science and Technology at Osaka
University. His research interests include parallel and distributed systems, software development tools, and performance evaluation.