在与我的 GSoC 导师讨论我们都认为 Asadpour 算法中最困难的部分——Held-Karp 松弛后，我们得出了一些结论。
- Asadpour 论文建议使用椭球法，以便他们的算法在多项式时间内运行。我们不需要多项式时间，只需要一个具有合理执行时间的算法。例如，椭球算法与单纯形算法。虽然单纯形算法是指数级的，但在实践中它几乎总是比椭球算法快。
- 我们对椭球算法的兴趣并非基于性能，而是基于椭球算法能够处理具有指数数量约束的线性规划的能力。这是通过分离预言机实现的，有关预言机的更多信息，请参阅我的文章此处。
- 实现一个健壮的椭球算法求解器（科学 Python 生态系统中缺少的显著内容）本身就是一个 GSoC 项目，并且超出了 NetworkX 此项目的范围。
因此，需要研究解决 Held-Karp 松弛的替代方法。为此，我们转向了 Held 和 Karp 在 1970 年的原始论文《旅行商问题和最小生成树》，以了解他们如何提出解决松弛的方法（请注意，这篇论文是在 1979 年将椭球算法应用于线性规划之前发表的）。Held 和 Karp 的论文讨论了三种解决松弛的方法
- **列生成：**一种较旧的解决非常大的线性规划的方法，其中只需要检查影响最优解的变量。
- **上升法：**一种基于最大化线性规划对偶的方法，其最佳描述是在类似于多元微积分中梯度概念的方式下寻找目标函数的上升方向。
- **分支定界：**这种方法具有最多的理论优势，并试图增强上升法以避免引入分数权重，分数权重是导致收敛速度慢的主要因素。
但在我们探索 Held 和 Karp 讨论的方法之前，我们需要确保这些方法仍然适用于在 Asadpour 论文的背景下解决 Held-Karp 松弛。我一直在本博客中使用的 Held-Karp 松弛的定义来自 Asadpour 论文的第 3 节，如下所示。
$$ \begin{array}{c l l} \text{min} & \sum_{a} c(a)x_a \\\ \text{s.t.} & x(\delta^+(U)) \geqslant 1 & \forall\ U \subset V \text{ and } U \not= \emptyset \\\ & x(\delta^+(v)) = x(\delta^-(v)) = 1 & \forall\ v \in V \\\ & x_a \geqslant 0 & \forall\ a \end{array} $$
Held Karp 论文中与此程序最接近的匹配是他们的线性规划 3，它是整个旅行商问题的线性规划表示，而不仅仅是松弛版本。请注意，Held 和 Karp 正在处理对称 TSP（STSP），而 Asadpour 则处理不对称或有向 TSP（ATSP）。
$$ \begin{array}{c l l} \text{min} & \sum_{1 \leq i < j \leq n} c_{i j}x_{i j} \\\ \text{s.t.} & \sum_{j > i} x_{i j} + \sum_{j < i} x_{j i} = 2 & (i = 1, 2, \dots, n) \\\ & \sum_{i \in S\\\ j \in S\\\ i < j} x_{i j} \leq |S| - 1 & \text{for any proper subset } S \subset {2, 3, \dots, n} \\\ & 0 \leq x_{i j} \leq 1 & (1 \leq i < j \leq n) \\\ & x_{i j} \text{integer} \\\ \end{array} $$
第二个线性规划中的最后两个约束是正确约束的，并且符合原始问题的范围，而前两个约束在找到 TSP 巡回路方面做了大部分工作。此外，将最后两个约束更改为 $x_{i j} \geq 0$ 就是 Held Karp 松弛。第一个约束，$\sum_{j > i} x_{i j} + \sum_{j < i} x_{j i} = 2$，确保在生成的巡回路中的每个顶点处，都有一条边到达那里，还有一条边离开。这与 Asadpour ATSP 松弛中的第二个约束相匹配。Held Karp 公式中的第二个约束是 Asadpour 线性规划中看到的子巡回消除约束的另一种形式。
Held 和 Karp 还指出
在本节中，我们证明最小化差距 $f(\pi)$ 等效于无需整数约束地求解此程序。
在第 1141 页，因此，似乎求解 Held 和 Karp 公式化的等效程序之一应该在这里起作用。
列生成技术#
列生成技术旨在求解 Held 和 Karp 论文中的线性规划 2，表示为
$$ \begin{array}{c l} \text{min} & \sum_{k} c_ky_k \\\ \text{s.t.} & y_k \geq 0 \\\ & \sum_k y_k = 1 \\\ & \sum_{i = 2}^{n - 1} (-v_{i k})y_k = 0 \\\ \end{array} $$
其中 $v_{i k}$ 是 1-树 $k$ 中顶点 $i$ 的度数减去 2，即 $v_{i k} = d_{i k} - 2$，每个变量 $y_k$ 对应于一个 1-树 $T^k$。每棵树的相关成本 $c_k$ 是 $T^k$ 的权重。
此方法的其余部分使用单纯形算法来求解线性规划。我们只关注每个 1-树中的边，使每一列的形式为
$$ \begin{bmatrix} 1 & -v_{2k} & -v_{3k} & \dots & -v_{n-1,k} \end{bmatrix}^T $$
以及进入解的列，该列位于 $c_k + \theta + \sum_{j=2}^{n-1} \pi_jv_{j k}$ 最小的 1-树中，其中 $\theta$ 和 $\pi_j$ 来自由 $(\theta, \pi_2, \pi_3, \dots, \pi_{n-1})$ 给出的“影子价格”向量。现在基是 $(n - 1) \times (n - 1)$，我们可以使用最小 1-树算法找到要添加到基中的 1-树，Held 和 Karp 说这可以在 $O(n^2)$ 步内完成。
我已经熟悉单纯形方法，因此我不会在这里详细介绍它的实现。
列生成技术的性能#
这种技术收敛速度很慢。Held 和 Karp 在 IBM/360 上对其进行了编程，并且能够解决最多 $n = 12$ 的一致问题。现在，在现代计算机上，时钟频率快了 210 到 101,500 倍（取决于使用的 IBM/360 模型），因此我们预计会有更好的性能，但目前还无法说明改进的程度。
他们还谈到了一个启发式过程，其中每当其相邻顶点的选择“明显”时，就会从程序中消除一个顶点。启发式的技术细节基本上不存在，但是
该程序在最多 $n = 48$ 的示例中显示出希望，但没有进行系统地探索。
上升法#
Held 和 Karp 的这篇论文是关于最小化 $f(\pi)$，其中 $f(\pi)$ 是排列的 1-树与 TSP 巡回路之间的差距。一种方法是最大化 $f(\pi)$ 的对偶，写成 $\text{max}_{\pi}\ w(\pi)$，其中
$$ w(\pi) = \text{min}_k\ (c_k + \sum_{i=1}^{i=n} \pi_iv_{i k}) $$
此方法使用一组 1-树的索引，这些 1-树相对于权重 $\overline{c}_{i j} = c_{i j} + \pi_i + \pi_j$ 具有最小权重。
$$ K(\pi) = {k\ |\ w(\pi) = c_k + \sum_{i=1}^{i=n} \pi_i v_{i k}} $$
如果 $\pi$ 不是 $w$ 的最大点，则将存在一个称为 $\pi$ 处上升方向的向量 $d$。这是定理 3，并在第 1148 页给出了证明。令函数 $\Delta(\pi, d)$ 和 $K(\pi, d)$ 定义如下。
$$ \Delta(\pi, d) = \text{min}_{k \in K(\pi)}\ \sum_{i=1}^{i=n} d_iv_{i k} \\\ K(\pi, d) = {k\ |\ k \in K(\pi) \text{ and } \sum_{i=1}^{i=n} d_iv_{i k} = \Delta(\pi, d)} $$
现在，对于足够小的 $\epsilon$，$K(\pi + \epsilon d) = K(\pi, d)$ 且 $w(\pi + \epsilon d) = w(\pi) + \epsilon \Delta(\pi, d)$，或者 $w(\pi)$ 的值增加，并且最小 1-树的增长率最小，因此我们保持低权重 1-树并进一步朝着最优值前进。最后，令 $\epsilon(\pi, d)$ 为以下数量
$$ \epsilon(\pi, d) = \text{max}\ {\epsilon\ |\text{ for } \epsilon’ < \epsilon,\ K(\pi + \epsilon’d = K(\pi, d)} $$
所以换句话说，$\epsilon(\pi, d)$ 是我们可以沿 $d$ 方向行进的最大距离，以保持所需的特性。
如果我们能找到 $d$ 和 $\epsilon$，那么我们就可以设置 $\pi = \pi + \epsilon d$ 并移动到上升法的下一轮迭代。Held 和 Karp 在第 1149 页给出了一个查找 $d$ 的协议。
- 将 $d$ 设置为零 $n$ 向量。
- 找到一个 1-树 $T^k$，使得 $k \in K(\pi, d)$。
- 如果 $\sum_{i=1}^{i=n} d_iv_{i k} > 0$ 停止。
- $d_i \leftarrow d_i + v_{i k},$ 对于 $i = 2, 3, \dots, n$
- 转到 2。
为了使该过程在 Python 中可实现，必须对其进行两方面的改进。
- 我们如何在步骤 2 中找到提到的 1-树？
- 我们如何知道何时不存在上升方向？（即，我们如何知道何时达到 $w(\pi)$ 的最大值？）
Held 和 Karp 针对这两个问题都提供了指导。在关于拟阵的第 6 节中，我们被告知要使用 Dijkstra 在《关于图的两个问题的注记》中开发的一种方法，但在这种特殊情况下，该方法并不是最有帮助的。我找到了这份文档，但 NetworkX 中已经有一个名为 minimum_spanning_arborescence 的函数，我们可以用它来创建一个最小 1-树形图。该过程是在 ${2, 3, \dots, n}$ 中的顶点上找到一个最小生成树形图，然后连接顶点 1 以创建循环。为了连接顶点 1，我们将选择成本最小的出边和成本最小的入边。
最后，在 $w(\pi)$ 的最大值处，没有上升方向，并且 Held 和 Karp 提出的程序将无法终止。他们的文章在第 1149 页指出
因此，当怀疑无法终止时，有必要检查是否存在上升方向；根据 Minkowski-Farkas 引理，这等价于存在非负系数 $\alpha_k$，使得
$ \sum_{k \in K(\pi)} \alpha_kv_{i k} = 0, \quad i = 1, 2, \dots, n $
这可以通过线性规划来检查。
虽然他们给出了这个求和式，但线性规划的其余部分也很有用。整个线性规划可以写成如下形式
$$ \begin{array}{c l l} \text{max} & \sum_k \alpha_k \\\ \text{s.t.} & \sum_{k \in K(\pi)} \alpha_k v_{i k} = 0 & \forall\ i \in {1, 2, \dots n} \\\ & \alpha_k \geq 0 & \forall\ k \\\ \end{array} $$
此线性规划不是标准形式，但将其转换并不困难。首先，通过最小化负值将最大化更改为最小化。
$$ \begin{array}{c l l} \text{min} & \sum_k -\alpha_k \\\ \text{s.t.} & \sum_{k \in K(\pi)} \alpha_k v_{i k} = 0 & \forall\ i \in {1, 2, \dots n} \\\ & \alpha_k \geq 0 & \forall\ k \\\ \end{array} $$
虽然约束直观上不是标准形式，但仔细观察就会发现它确实是。矩阵形式中的每一列都对应于 $\alpha_k$ 的一个条目，每一行都代表一个不同的 $i$ 值，或一个不同的顶点。一个约束实际上是一组非常相似的约束，可以写成
$$ \begin{array}{c l} \text{min} & \sum_k -\alpha_k \\\ \text{s.t.} & \sum_{k \in K(\pi)} \alpha_k v_{1 k} = 0 \\\ & \sum_{k \in K(\pi)} \alpha_k v_{2 k} = 0 \\\ & \vdots \\\ & \sum_{k \in K(\pi)} \alpha_k v_{n k} = 0 \\\ & \alpha_k \geq 0 & \forall\ k \\\ \end{array} $$
因为所有求和必须等于零，所以不需要松弛变量和剩余变量，因此该程序的约束矩阵为 $n \times k$。$n$ 显然具有线性增长率，但我不知道 $k$ 会增长到多大。$k$ 是最小 1-树的集合，所以我认为它是可以管理的。此线性规划可以使用 SciPy 库中内置的 linprog 函数求解。
作为实现说明，首先我可能会在每次迭代中检查终止条件，但最终我们可以找到它必须执行的迭代次数，然后才能开始检查终止条件以节省计算能力。
终止条件的一个可能困难是，我们需要使用来自每个最小 1-树或 1-树形图的数据运行线性规划，这意味着我们需要能够生成所有最小 1-树。目前，在 NetworkX 中似乎没有简单的方法可以做到这一点。查看此处 树算法，它们似乎只专注于找到所需类型的一个最小分支，而不是所有这些分支。
现在我们必须找到 $\epsilon$。第 1150 页上的定理 4 指出
令 $k$ 为 $K(\pi, d)$ 的任何元素，其中 $d$ 是 $\pi$ 处的上升方向。则 $\epsilon(\pi, d) = \text{min}{\epsilon\ |\text{ 对于某些对 } (e, e’),\ e’ \text{ 是 } e \text{ 在 } T^k \text{ 中的替代 } \\\ \text{ 并且 } e \text{ 和 } e’ \text{ 在 } \epsilon \text{ 处交叉 }}$
第一步是确定 $e$ 和 $e’$ 是否为替代。如果对于 1-树 $T^k$，$(T^k - {e}) \cup {e’}$ 也是 1-树，则 $e’$ 为替代。如果对 $(\overline{c}_{i j}, d_i + d_j)$ 和 $(\overline{c}_{r s}, d_r + d_s)$ 不同，但边 $e = {r, s}$ 和 $e’ = {i, j}$ 在 $\epsilon$ 处交叉
$$ \overline{c}_{i j} + \epsilon(d_i + d_j) = \overline{c}_{r s} + \epsilon(d_r + d_s) $$
从该方程，我们可以推导出 $\epsilon$ 的公式。
$$ \begin{array}{r c l} \overline{c}_{i j} + \epsilon(d_i + d_j) &=& \overline{c}_{r s} + \epsilon(d_r + d_s) \\\ \epsilon(d_i + d_j) &=& \overline{c}_{r s} + \epsilon(d_r + d_s) - \overline{c}_{i j} \\\ \epsilon(d_i + d_j) - \epsilon(d_r + d_s) &=& \overline{c}_{r s} - \overline{c}_{i j} \\\ \epsilon\left((d_i + d_j) - (d_r + d_s)\right) &=& \overline{c}_{r s} - \overline{c}_{i j} \\\ \epsilon(d_i + d_j - d_r - d_s) &=& \overline{c}_{r s} - \overline{c}_{i j} \\\ \epsilon &=& \displaystyle \frac{\overline{c}_{r s} - \overline{c}_{i j}}{d_i + d_j - d_r - d_s} \end{array} $$
因此，我们现在可以找到任何两对彼此替代的边的 $\epsilon$，但我们需要能够在 1-树中找到替代。我们知道当且仅当 $e$ 和 $e’$ 都与顶点 1 相邻或 $e$ 位于不通过顶点 1 的 $T^k \cup {e’}$ 的循环中时，$e’$ 为 $e$ 的替代。更正式地说，我们试图在与 $e’$ 相同的基本循环中找到边。当任何不在生成树中的边添加到该生成树时，就会创建一个基本循环。因为此边的端点通过一条唯一的路径连接，所以会创建一个唯一的循环。为了找到这个循环，我们将利用 NetworkX 库中的 find_cycle。
下面是我草拟的利用定理 4 查找 $\epsilon(\pi, d)$ 的伪代码过程。它没有经过很好的优化，但会找到 $\epsilon(\pi, d)$。
# Input: An element k of K(pi, d), the vector pi and the vector d.
# Output: epsilon(pi, d) using Theorem 4 on page 1150.
for each edge e in the graph G
if e is in k:
continue
else:
add e to k
let v be the terminating end of e
c = find_cycle(k, v)
for each edge a in c not e:
if a[cost] = e[cost] and d[i] + d[j] = d[r] + d[s]:
continue
epsilon = (a[cost] - e[cost])/(d[i] + d[j] - d[r] - d[s])
min_epsilon = min(min_epsilon, epsilon)
remove e from k
return min_epsilon
上升方法的性能#
上升方法也很慢，但在现代计算机上会更好。当 Held 和 Karp 对其进行编程时，他们在一些最多 25 个顶点的小问题上对其进行了测试，虽然每次迭代的时间很短，但迭代次数却迅速增加。他们没有评论这种方法是否比列生成技术更好，但确实指出他们没有确定这种方法是否始终收敛到 $w(\pi)$ 的最大点。
分支定界方法#
在与我的 GSoC 导师交谈后，我们认为这是我们可以为 Held-Karp 松弛（Asadpour 算法所需）实现的最佳方法。上升方法嵌套在此方法中，因此需要深入探讨之前的方法才能实现此方法。此方法中的大多数符号都从上升方法中重复使用。
分支定界方法利用了顶点可能出现不平衡的概念。如果
$$ \forall\ k \in K(\pi),\ v_{i k} \geq 1 $$
则顶点 $i$ 为高不平衡。类似地，如果
$$ \forall\ k \in K(\pi),\ v_{i k} = -1 $$
则顶点 $i$ 为低不平衡。请记住，$v_{i k}$ 是顶点的度数减 2。我们知道所有顶点的度数至少为 1，否则 1-树 $T^k$ 将不会连通。高不平衡顶点在每个最小 1-树中的度数为 3 或更高，而低不平衡顶点在所有最小 1-树中的度数仅为 1。我们的目标是每个顶点的度数都为 2 的最小 1-树。
如果我们知道某个顶点在任一方向上都是不平衡的，那么我们就知道上升方向，并且该方向是单位向量。令 $u_i$ 为一个 $n$ 维单位向量，在第 $i$ 个坐标处为 1。如果顶点 $i$ 为高不平衡，则 $u_i$ 为上升方向；如果顶点 $i$ 为低不平衡，则 $-u_i$ 为上升方向。
第 1151 页上的推论 3 和 4 也表明，当顶点不平衡时，查找 $\epsilon(\pi, d)$ 也更简单。
推论 3. 假设顶点 $i$ 为低不平衡，并令 $k$ 为 $K(\pi, -u_i)$ 的元素。则 $\epsilon(\pi, -u_i) = \text{min} (\overline{c}_{i j} - \overline{c}_{r s})$，使得 ${i, j}$ 是 $T^k$ 中 ${r, s}$ 的替代，并且 $i \not\in {r, s}$。
推论 4. 假设顶点 $r$ 为高不平衡。则 $\epsilon(\pi, u_r) = \text{min} (\overline{c}_{i j} - \overline{c}_{r s})$，使得 ${i, j}$ 是 $T^k$ 中 ${r, s}$ 的替代，并且 $r \not\in {i, j}$。
这些推论可以使用上面在上升方法部分查找 $\epsilon$ 的伪代码列表的修改版本来实现。
一旦不再存在不平衡顶点，上升方向就不是单位向量，并且会引入分数权重。这是上升方法收敛到最优解速度大幅下降的原因，因此应尽可能避免这种情况。
在我们讨论实现细节之前，仍有一些主要内容需要回顾。令 $X$ 和 $Y$ 为图中不相交的边集。然后令 $\mathsf{T}(X, Y)$ 表示包含 $X$ 中的所有边但不包含 $Y$ 中的任何边的 1-树集。最后，定义 $w_{X, Y}(\pi)$ 和 $K_{X, Y}(\pi)$ 如下。
$$ w_{X, Y}(\pi) = \text{min}_{k \in \mathsf{T}(X, Y)} (c_k + \sum_{i=1}^{i=n} \pi_i v_{i k}) \\\ K_{X, Y}(\pi) = {k\ |\ c_k + \sum \pi_i v_{i k} = w_{X, Y}(\pi)} $$
从这些函数中，出现了高不平衡和低不平衡的修订定义，允许顶点相对于 $X$ 和 $Y$ 不平衡。
在分支定界方法完成过程中，分支在一个列表中跟踪，其中每个条目具有以下格式。
$$[X, Y, \pi, w_{X, Y}(\pi)]$$
其中 $X$ 和 $Y$ 是前面讨论的不相交集，$\pi$ 是我们用来扰动边权重的向量，$w_{X, Y}(\pi)$ 是条目的界限。
在方法的每次迭代中，我们考虑具有最小界限的列表条目，并尝试查找不平衡顶点。如果我们找到一个，我们将使用简化的单位向量作为上升方向，应用一次上升方法迭代。在这里，如果存在整数权重，我们可以利用它们。也许 NetworkX 中 Asadpour 实现的文档应该说明整数边权重会表现更好，但该说法需要由我们的测试来支持。
如果不存在失衡顶点，我们仍然需要找到上升方向以确定是否处于 $w(\pi)$ 的最大值。如果存在上升方向，则进行分支。如果没有上升方向，则在 $K_{X, Y}(\pi)$ 中搜索一个巡回路线，如果未找到，也进行分支。
分支过程如下。从条目 $[X, Y, \pi, w_{X, Y}(\pi)]$ 中选择一条边 $e \not\in X \cup Y$（Held 和 Karp 没有给出任何分支标准，因此我相信选择可以是任意的），并将父条目替换为两种形式的其他两个条目
$$ [X \cup {e}, Y^*, \pi, w_{X \cup {e}, Y^*}(\pi)] \quad \text{和} \quad [X^*, Y \cup {e}, \pi, w_{X^*, Y \cup {e}}(\pi)] $$
在 Held 和 Karp 论文的第 1153 页到第 1156 页给出了分支定界方法的一个例子。
为了实现此方法，除了修改上升方法的一些细节外，我们还需要能够确定。
- 如果一个顶点相对于 $X$ 和 $Y$ 在任一方向上都处于失衡状态。
- 在 $K_{X, Y}(\pi)$ 中搜索一个巡回路线。
Held 和 Karp 论文指出，为了找到失衡顶点，我们只需要测试单位向量即可。如果对于 $K(\pi, u_i)$ 的任意成员 $k$，$v_{i k} \geq 1$ 并且对于失衡低点成立相应的逆关系。通过此过程，我们可以通过按顺序检查 $u_i$ 以 $O(n^2)$ 的过程找到失衡顶点。
如果我们可以枚举该集合中的最小 1-树，那么在 $K_{X, Y}(\pi)$ 中搜索一个巡回路线将很容易。虽然我知道如何找到其中一个最小 1-树或 $K(\pi)$ 的一个成员，但我不知道如何找到 $K(\pi, d)$ 中的元素，甚至不知道如何找到 $K(\pi)$ 的所有成员。使用 Held 和 Karp 论文中的属性，我确实知道如何将 $K(\pi)$ 细化为 $K(\pi, d)$ 以及将 $K(\pi)$ 细化为 $K_{X, Y}(\pi)$。这将是以后博客文章的主题。
我能够找到的关于此问题的最有希望的研究论文是 Sörensen 和 Janssens 在 2005 年发表的题为《An Algorithm to Generate all Spanning Trees of a Graph in Order of Increasing Cost》的论文 这里。从这里，我们生成生成树或树枝直到成本上升，此时我们已经找到了 $K(\pi)$ 的所有元素。
分支定界方法的性能#
Held 和 Karp 没有编写此方法的程序。我们有一些理由相信此方法的性能将是最佳的，因为它是为了改进上升方法而设计的，而上升方法经过了（某种程度上的）测试，直到 $n = 25$，这仍然优于列生成技术，列生成技术只能持续解决到 $n = 12$。
参考文献#
A. Asadpour、M. X. Goemans、A. Mardry、S. O. Ghran 和 A. Saberi，《An o(log n / log log n)-approximation algorithm for the asymmetric traveling salesman problem》，运筹学，65（2017），第 1043-1061 页，https://homes.cs.washington.edu/~shayan/atsp.pdf。
Held, M., Karp, R.M. 《The traveling-salesman problem and minimum spanning trees》。运筹学，1970 年 11 月 1 日，第 18 卷（第 6 期），第 1138-1162 页。 https://www.jstor.org/stable/169411