TSP问题从DP算法到深度学习1： 递归DP方法 AC AIZU TSP问题

旅行商问题（TSP）是计算机算法中经典的NP hard 问题。 在本系列文章中，我们将首先使用动态规划 AC aizu中的TSP问题，然后再利用深度学习求大规模下的近似解。深度学习应用解决问题时先以PyTorch实现监督学习算法 Pointer Network，进而结合强化学习来无监督学习，提高数据使用效率。 本系列完整列表如下：

• 第一篇: 递归DP方法 AC AIZU TSP问题

• 第二篇: 二维空间TSP数据集及其DP解法

• 第三篇: 深度学习 Pointer Networks 的 Pytorch实现

• 第四篇: 搜寻最有可能路径：Viterbi算法和其他

• 第五篇: 深度强化学习无监督算法的 Pytorch实现

TSP 问题回顾

TSP可以用图模型来表达，无论有向图或无向图，无论全连通图或者部分连通的图都可以作为TSP问题。 Wikipedia TSP 中举了一个无向全连通的TSP例子。如下图所示，四个顶点A，B，C，D构成无向全连通图。TSP问题要求在所有遍历所有点后返回初始点的回路中找到最短的回路。例如，\(A \rightarrow B \rightarrow C \rightarrow D \rightarrow A\) 和 \(A \rightarrow C \rightarrow B \rightarrow D \rightarrow A\) 都是有效的回路，但是TSP需要返回这些回路中的最短回路（注意，最短回路可能会有多条）。

Wikipedia 4个顶点组成的图

无论是哪种类型的图，我们都能用邻接矩阵表示出一个图。上面的Wikipedia中的图可以用下面的矩阵来描述。

\[ \begin{matrix} & \begin{matrix}A&B&C&D\end{matrix} \\\\ \begin{matrix}A\\\\B\\\\C\\\\D\end{matrix} & \begin{bmatrix}-&20&42&35\\\\20&-&30&34\\\\42&30&-&12\\\\35&34&12&-\end{bmatrix}\\\\ \end{matrix} \]

当然，大多数情况下，TSP问题会被限定在欧氏空间，即二维地图中的全连通无向图。因为，如果将顶点表示一个地理位置，一般来说它可以和其他所有顶点连通，回来的距离相同，由此构成无向图。

AIZU TSP 问题

AIZU在线题库 有一道有向不完全连通图的TSP问题。给定V个顶点和E条边，输出最小回路值。例如，题目里的例子如下所示，由4个顶点和6条单向边构成。

这个示例的答案是16，对应的回路是 \(0\rightarrow1\rightarrow3\rightarrow2\rightarrow0\)，由下图的红色边构成。注意，这个题目可能不存在合法解，原因是无回路存在，此时返回-1，可以合理地理解成无穷大。

暴力解法

一种暴力方法是枚举所有可能的从某一顶点的回路，取其中的最小值即可。下面的 Python 示例如何枚举4个顶点构成的图中从顶点0出发的所有回路。

{linenos

from itertools import permutations
v = [1,2,3]
p = permutations(v)
for t in list(p):
  print([0] + list(t) + [0])

所有从顶点0出发的回路如下：

{linenos

[0, 1, 2, 3, 0]
[0, 1, 3, 2, 0]
[0, 2, 1, 3, 0]
[0, 2, 3, 1, 0]
[0, 3, 1, 2, 0]
[0, 3, 2, 1, 0]

很显然，这种方式的时间复杂度是 O(\(n!\))，无法通过AIZU。

动态规划求解

我们可以使用位状态压缩的动态规划来AC这道题。 首先，需要将回路过程中的状态编码成二进制的表示。例如，在四顶点的例子中，如果顶点2和1都被访问过，并且此时停留在顶点1。将已经访问的顶点对应的位置1，那么编码成0110，此外，还需要保存当前顶点的位置，因此我们将代表状态的数组扩展成二维，第一维是位状态，第二维是顶点所在位置，即 \(dp[bitstate][v]\)。这个例子的状态表示就是 \(dp["0110"][1]\)。

状态转移方程如下： \[ dp[bitstate][v] = \min ( dp[bitstate \cup \{u\}][u] + dist(v,u) \mid u \notin bitstate ) \] 这种方法对应的时间复杂度是 O(\(n^2*2^n\) )，因为总共有 \(2^n * n\) 个状态，而每个状态又需要一次遍历。虽然都是指数级复杂度，但是它们的巨大区别由下面可以看出区别。

	\(n!\)	\(n^2*2^n\)
n=8	40320	16384
n=10	3628800	102400
n=12	479001600	589824
n=14	87178291200	3211264

暂停思考一下为什么状态压缩DP能工作。注意到之前暴力解法中其实是有很多重复计算，下面红圈表示重复的计算节点。

在本篇中，我们将会用Python 3和Java 8 实现自顶向下的DP 缓存版本。这种方式比较符合直觉，因为我们不需要预先考虑计算节点的依赖关系。在Java中我们使用了一个小技巧，dp数组初始化成Integer.MAX_VALUE，如此只需要一条语句就能完成更新dp值。

res = Math.min(res, s + g.edges[v][u]);

当然，为了AC 这道题，我们需要区分出真正无法到达的情况并返回-1。 在Python实现中，也可以使用同样的技巧，但是这次示例一般的实现方法：将dp数组初始化成-1并通过 if-else 来区分不同情况。

INT_INF = -1

if s != INT_INF and edges[v][u] != INT_INF:
    if ret == INT_INF:
        ret = s + edges[v][u]
    else:
        ret = min(ret, s + edges[v][u])

下面附完整的Python 3和Java 8的AC代码，同步在 github。

AIZU Java 8 递归DP版本

{linenos

// passed http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=DPL_2_A
import java.util.Arrays;
import java.util.Scanner;

public class Main {
    public static class Graph {
        public final int V_NUM;
        public final int[][] edges;

        public Graph(int V_NUM) {
            this.V_NUM = V_NUM;
            this.edges = new int[V_NUM][V_NUM];
            for (int i = 0; i < V_NUM; i++) {
                Arrays.fill(this.edges[i], Integer.MAX_VALUE);
            }
        }
    
        public void setDist(int src, int dest, int dist) {
            this.edges[src][dest] = dist;
        }
    
    }
    
    public static class TSP {
        public final Graph g;
        long[][] dp;
    
        public TSP(Graph g) {
            this.g = g;
        }
    
        public long solve() {
            int N = g.V_NUM;
            dp = new long[1 << N][N];
            for (int i = 0; i < dp.length; i++) {
                Arrays.fill(dp[i], -1);
            }
    
            long ret = recurse(0, 0);
            return ret == Integer.MAX_VALUE ? -1 : ret;
        }
    
        private long recurse(int state, int v) {
            int ALL = (1 << g.V_NUM) - 1;
            if (dp[state][v] >= 0) {
                return dp[state][v];
            }
            if (state == ALL && v == 0) {
                dp[state][v] = 0;
                return 0;
            }
            long res = Integer.MAX_VALUE;
            for (int u = 0; u < g.V_NUM; u++) {
                if ((state & (1 << u)) == 0) {
                    long s = recurse(state | 1 << u, u);
                    res = Math.min(res, s + g.edges[v][u]);
                }
            }
            dp[state][v] = res;
            return res;
    
        }
    
    }
    
    public static void main(String[] args) {
    
        Scanner in = new Scanner(System.in);
        int V = in.nextInt();
        int E = in.nextInt();
        Graph g = new Graph(V);
        while (E > 0) {
            int src = in.nextInt();
            int dest = in.nextInt();
            int dist = in.nextInt();
            g.setDist(src, dest, dist);
            E--;
        }
        System.out.println(new TSP(g).solve());
    }
}

AIZU Python 3 递归DP版本

{linenos

from typing import List

INT_INF = -1

class Graph:
    v_num: int
    edges: List[List[int]]

    def __init__(self, v_num: int):
        self.v_num = v_num
        self.edges = [[INT_INF for c in range(v_num)] for r in range(v_num)]
    
    def setDist(self, src: int, dest: int, dist: int):
        self.edges[src][dest] = dist


class TSPSolver:
    g: Graph
    dp: List[List[int]]

    def __init__(self, g: Graph):
        self.g = g
        self.dp = [[None for c in range(g.v_num)] for r in range(1 << g.v_num)]
    
    def solve(self) -> int:
        return self._recurse(0, 0)
    
    def _recurse(self, v: int, state: int) -> int:
        """
    
        :param v:
        :param state:
        :return: -1 means INF
        """
        dp = self.dp
        edges = self.g.edges
    
        if dp[state][v] is not None:
            return dp[state][v]
    
        if (state == (1 << self.g.v_num) - 1) and (v == 0):
            dp[state][v] = 0
            return dp[state][v]
    
        ret: int = INT_INF
        for u in range(self.g.v_num):
            if (state & (1 << u)) == 0:
                s: int = self._recurse(u, state | 1 << u)
                if s != INT_INF and edges[v][u] != INT_INF:
                    if ret == INT_INF:
                        ret = s + edges[v][u]
                    else:
                        ret = min(ret, s + edges[v][u])
        dp[state][v] = ret
        return ret


def main():
    V, E = map(int, input().split())
    g: Graph = Graph(V)
    for _ in range(E):
        src, dest, dist = map(int, input().split())
        g.setDist(src, dest, dist)

    tsp: TSPSolver = TSPSolver(g)
    print(tsp.solve())


if __name__ == "__main__":
    main()