Traveling Salesman Problem（TSP），旅行商问题

在Hamiltonion Cycle的问题上，再增加一个条件，即要求cost最小，就是经典的TSP（Traveling Salesman Problem）问题。经过所有的点仅一次，最后回到原点，所经过的路程cost最小。

下面演示回溯算法，时间复杂度$O(N!)$，不过空间复杂度很低，只有一些调用栈的开销，$O(N)$。假设一个Complete Graph，有n个节点，用代码随机产生edge cost。

import random
import time


N = 4
edge = [[0]*N for _ in range(N)]
for i in range(N):
    for j in range(1,N):
        if i == j:
            edge[i][j] = edge[j][i] = 0  # simple graph, no self-loop
        else:
            edge[i][j] = edge[j][i] = random.randint(10,99)

print('node number', N, ', edge matrix:')
for i in range(N):
    for j in range(N):
        print(edge[i][j], end=' ')
    print()


def tsp_backtracking(V):
    """ TSP by backtracking, O(N!), https://cs.pynote.net """
    assert V < N
    mcost = None
    mpath = []

    def __tbt(stage, path, pcost):
        nonlocal mcost
        nonlocal mpath
        if stage == N:
            print('#', path, pcost)
            # time.sleep(2)
            if mcost is None or mcost>pcost:
                mcost = pcost
                mpath = path[:]
            return
        if len(path) == N:  # back to start
            endv = path[-1]
            __tbt(stage+1, path+[V], pcost+edge[endv][V])
            return
        for i in range(N):
            endv = path[-1]
            if i not in path:
                path.append(i)
                __tbt(stage+1, path, pcost+edge[endv][i])
                path.pop()

    __tbt(0, [V], 0)
    #c = 0
    #for i in range(len(mpath)-1):
    #    c += edge[mpath[i]][mpath[i+1]]
    #assert c == mcost
    return mcost, mpath


for i in range(N):
    print(tsp_backtracking(i))

node number 4 , edge matrix:
0 19 77 12
19 0 19 25
77 19 0 40
12 25 40 0
# [0, 1, 2, 3, 0] 90
# [0, 1, 3, 2, 0] 161
# [0, 2, 1, 3, 0] 133
# [0, 2, 3, 1, 0] 161
# [0, 3, 1, 2, 0] 133
# [0, 3, 2, 1, 0] 90
(90, [0, 1, 2, 3, 0])
# [1, 0, 2, 3, 1] 161
# [1, 0, 3, 2, 1] 90
# [1, 2, 0, 3, 1] 133
# [1, 2, 3, 0, 1] 90
# [1, 3, 0, 2, 1] 133
# [1, 3, 2, 0, 1] 161
(90, [1, 0, 3, 2, 1])
# [2, 0, 1, 3, 2] 161
# [2, 0, 3, 1, 2] 133
# [2, 1, 0, 3, 2] 90
# [2, 1, 3, 0, 2] 133
# [2, 3, 0, 1, 2] 90
# [2, 3, 1, 0, 2] 161
(90, [2, 1, 0, 3, 2])
# [3, 0, 1, 2, 3] 90
# [3, 0, 2, 1, 3] 133
# [3, 1, 0, 2, 3] 161
# [3, 1, 2, 0, 3] 133
# [3, 2, 0, 1, 3] 161
# [3, 2, 1, 0, 3] 90
(90, [3, 0, 1, 2, 3])

动态规划

如果使用动态规划，可以把复杂度降下来一些，但还是指数级的增长！

$TSP(i,V_i)$表示从vertex i开始，经过$V_i$这个vertex集合，最后到达起点的cost。假设1为起点，TSP问题可以表达为：

$$\begin{cases} TSP(1,V_1)=min_{v\in Nbr(1)}\{w(1,v)+TSP(v,V_1-v)\} \\ 1\notin V_1 \\ TSP(u,\empty)=w(u,1) \\ \end{cases}$$

对于Complete Graph而言，可以认为subproblem的数量为$2^n$，即顶点的任意数量组合数，每一个顶点组合最多有$n$个前缀顶点，min操作再产生一个$n$，因此最后复杂度为：$O(N^2\cdot 2^N)$。

s is start vertex
V is a set of all vertice without s

TSP(p,V):
    if V is empty:
        return w(p,s)
    cost = []
    for each v in V:
        cost.append(w(p,v)+TSP(v,V-{v}))
    return min(cost)

TSP(s,V) is TSP

Top-Down

class tsp_dp:

    def __init__(self, start, remain):
        self.s = start
        self.rem = remain
        self.rec = {}

    def _tsp(self, v, rem):
        """ TSP by Top-Down DP, O(N^2 2^N), https://cs.pynote.net """
        if len(rem) == 0:
            return edge[v][self.s]
        assert v not in rem
        key = (v, frozenset(rem))
        if key in self.rec:
            return self.rec[key]
        mcost = None
        for i in rem:
            cost = edge[v][i] + self._tsp(i, rem-set({i}))
            if mcost is None or mcost>cost:
                mcost = cost
            rem.add(i)
        self.rec[key] = mcost
        return mcost

    def go(self):
        return self._tsp(self.s, self.rem)


for i in range(N):
    tsp = tsp_dp(i, set(j for j in range(N) if j!=i))
    print(tsp.go())

动态规划的核心是，想清楚subproblem是什么，以及如何利用subproblem。TSP问题的subproblem，是一个典型的permutation转combination后避免了一些重复计算的典型。比如当某个节点与k节点连接，以k节点开始的各种其它节点的组合，就是subproblem。这些组合对应的值可以被重复使用。这些组合是不同size的。如何选择组合，完全根据当前状态下剩余还没有使用过的节点。组合实际上隐藏了内部结果，只暴露出k作为头结点。组合数为0就是base case，此时需要回到最开始的节点。$O(N\cdot 2^N)$就是组合数，另外一个$N$是在选择minimum时产生的。从$O(N!)$到$O(N^2\cdot 2^N)$，快了许多倍，代价是存储空间需求的暴涨。

Bottom-Up

def tsp_bu(V, rem):
    """ TSP by Bottom-up DP, O(N^2 2^N), https://cs.pynote.net """
    assert V not in rem
    rec = {(i,frozenset({})):edge[i][V] for i in rem}
    for i in range(1,N-1):
        trec = {}
        for k in rem:
            for r,c in rec.items():
                a = (r[0],) + tuple(r[1])
                if k in a:
                    continue
                cost = edge[k][r[0]] + c
                key = (k, frozenset(a))
                if key not in trec or trec[key]>cost:
                    trec[key] = cost
        rec = trec
    mcost = None
    for r,c in rec.items():
        cost = edge[V][r[0]] + c
        if mcost is None or mcost>cost:
            mcost = cost
    return mcost