了解上升法

距离我上次写博客文章已经很久了，比我预期的要长。正如我在最初的 GSoC 提案中所预期的那样，Held-Karp 松弛证明很难实现。

我和我的导师一致认为，Held 和 Karp 在 1970 年的论文《旅行推销员问题与最小生成树》中讨论的分支定界法，首先需要实现上升法，因为它被用在分支定界法中。在过去的一周半时间里，我一直致力于实现和调试上升法，并且想花一些时间反思我所学到的东西。

首先我要说，截至撰写本文时，我的上升法版本并没有给出我期望的最佳解决方案。为了进行测试，我使用了 Held 和 Karp 在分支定界法示例中使用的图，一个加权的 $\mathcal{K}_6$，并将其转换为下面的邻接矩阵中给出的有向但对称的版本。

$$ \begin{bmatrix} 0 & 97 & 60 & 73 & 17 & 52 \\\ 97 & 0 & 41 & 52 & 90 & 30 \\\ 60 & 41 & 0 & 21 & 35 & 41 \\\ 73 & 52 & 21 & 0 & 95 & 46 \\\ 17 & 90 & 35 & 95 & 0 & 81 \\\ 52 & 30 & 41 & 46 & 81 & 0 \end{bmatrix} $$

原始解是一个无向的巡回，但在有向版本中，期望的解取决于它们的遍历方式。这两个循环的总权重均为 207。

这是程序返回的循环，其总权重为 246。

所有这些代码都进入了 NetworkX 中 traveling_saleaman.py 文件内的 _held_karp() 函数，我试图尽可能地遵循论文中概述的算法。_held_karp() 函数本身有三个内部函数，k_pi()、direction_of_ascent() 和 find_epsilon()，它们代表了上升法每次迭代中使用的三个主要步骤。

k_pi()#

k_pi() 使用我在编码第一周为 Google Summer of Code 实现的 ArborescenceIterator 来查找图中的所有最小 1-树。我对创建 1-树的最初评估略有错误。我曾说过

为了连接顶点 1，我们将选择具有最小成本的出边和具有最小成本的入边。

实际上，这种方法将生成几乎是树的图，仅仅是因为出边几乎肯定会创建一个具有两个入边的顶点。相反，我们需要使用最低成本的入边连接顶点 1，以及连接到节点 ${2, 3, \dots, n}$ 上树的根节点的边，这样就不会违反入度约束。

对于测试图，在上升法的第一次迭代中，k_pi() 返回了 10 个 1-树，但它们的成本并不都相同。请注意，因为我们在选择顶点 1 的出边时没有自由，所以 1-树的总成本将因连接到根节点的最便宜节点与最昂贵的节点之间的差异而异。我最初编写的这个函数效率不高，它从所有最小生成树中创建了 1-树，然后对其进行迭代以删除所有非最小树。

昨天，我重写了这个函数，以便一旦找到权重较低的 1-树，它就会删除所有当前的最小树，以支持新的树，并且不会将权重更大的任何 1-树添加到最小 1-树集中。

我重写该方法的真正原因是尝试一些新方法，希望将程序从次优解推向最优解。正如我之前提到的，强制选择连接到根节点的方式会创建不同权重的 1-树。因此我怀疑，选择不同的顶点 1 可能会创建比仅使用 next(G.__iter__()) 返回的顶点所创建的 1-树权重更低的 1-树。所以我用一个遍历图中顶点的 for 循环将整个 k_pi() 函数包装起来，发现顶点 1 的选择会产生影响。

Excluded node: 0, Total Weight: 161.0
Chosen incoming edge for node 0: (4, 0), chosen outgoing edge for node 0: (0, 4)
(2, 3, 21)
(2, 5, 41)
(4, 2, 35)
(4, 0, 17)
(5, 1, 30)
(0, 4, 17)

Excluded node: 0, Total Weight: 161.0
Chosen incoming edge for node 0: (4, 0), chosen outgoing edge for node 0: (0, 4)
(1, 5, 30)
(2, 1, 41)
(2, 3, 21)
(4, 2, 35)
(4, 0, 17)
(0, 4, 17)

Excluded node: 1, Total Weight: 174.0
Chosen incoming edge for node 1: (5, 1), chosen outgoing edge for node 1: (1, 5)
(2, 3, 21)
(2, 4, 35)
(4, 0, 17)
(5, 2, 41)
(5, 1, 30)
(1, 5, 30)

Excluded node: 2, Total Weight: 187.0
Chosen incoming edge for node 2: (3, 2), chosen outgoing edge for node 2: (2, 3)
(0, 4, 17)
(3, 5, 46)
(3, 2, 21)
(5, 0, 52)
(5, 1, 30)
(2, 3, 21)

Excluded node: 3, Total Weight: 165.0
Chosen incoming edge for node 3: (2, 3), chosen outgoing edge for node 3: (3, 2)
(1, 5, 30)
(2, 1, 41)
(2, 4, 35)
(2, 3, 21)
(4, 0, 17)
(3, 2, 21)

Excluded node: 3, Total Weight: 165.0
Chosen incoming edge for node 3: (2, 3), chosen outgoing edge for node 3: (3, 2)
(2, 4, 35)
(2, 5, 41)
(2, 3, 21)
(4, 0, 17)
(5, 1, 30)
(3, 2, 21)

Excluded node: 4, Total Weight: 178.0
Chosen incoming edge for node 4: (0, 4), chosen outgoing edge for node 4: (4, 0)
(0, 5, 52)
(0, 4, 17)
(1, 2, 41)
(2, 3, 21)
(5, 1, 30)
(4, 0, 17)

Excluded node: 4, Total Weight: 178.0
Chosen incoming edge for node 4: (0, 4), chosen outgoing edge for node 4: (4, 0)
(0, 5, 52)
(0, 4, 17)
(2, 3, 21)
(5, 1, 30)
(5, 2, 41)
(4, 0, 17)

Excluded node: 5, Total Weight: 174.0
Chosen incoming edge for node 5: (1, 5), chosen outgoing edge for node 5: (5, 1)
(1, 2, 41)
(1, 5, 30)
(2, 3, 21)
(2, 4, 35)
(4, 0, 17)
(5, 1, 30)

请注意，由于我的测试图是对称的，它喜欢创建只有两个节点的循环。这些 1-树的权重范围从 161 到 178，所以我尝试运行使用新方法进行的测试，该测试大约需要 300 毫秒... 并且程序没有终止。我在上升法迭代 200 次后在 PyCharm 中创建了断点，发现程序陷入了一个循环中，它在两个不同的最小 1-树之间交替。这只是一个尝试，但没有奏效，所以我将代码恢复为始终为顶点 1 选择相同的顶点。

无论如何，我几乎完全重写了这个函数，而输出没有变化，这表明这个函数不是问题的根源。

direction_of_ascent()#

这是 Held 和 Karp 论文中包含伪代码的一个函数

1. 将 $d$ 设置为零 $n$ 维向量。
2. 找到一个 1-树 $T^k$，使得 $k \in K(\pi, d)$。[执行步骤 2 的方法来自第 6 节的结果（贪婪算法）。]
3. 如果 $\sum_{i=1}^{i=n} d_i v_{i k} > 0$，则停止。
4. $d_i \rightarrow d_i + v_{i k}$，对于 $i = 2, 3, \dots, n$
5. 转到 2。

以它为指导，这个函数的实现很简单，直到我遇到了终止条件，这是一个线性规划问题，在第 1149 页进行了讨论，如下所示

因此，当怀疑无法终止时，有必要检查是否存在上升方向；根据 Minkowski-Farkas 引理，这等同于存在非负系数 $\alpha_k$，使得

$ \sum_{k \in K(\pi)} \alpha_kv_{i k} = 0, \quad i = 1, 2, \dots, n $

可以通过线性规划来检查这一点。

虽然我能够在没有太多问题的情况下实现这一点，但线性规划的一个非常重要的约束条件在这里没有提到，而是在前一页的证明中提到了。该约束条件是

$$ \sum_{k \in K(\pi)} \alpha_k = 1 $$

在我花了好几个小时试图调试原始线性规划并注意到缺少的约束条件后。线性规划开始正常工作，在找到巡回时终止程序。

find_epsilon()#

与 Held 和 Karp 论文中描述的方法相比，这个函数需要完全不同的实现。

我的有向图实现和无向图描述中的基本思想是找到相互替代的边，或者找到不在 1-树中的边，它可以替换树中的边，并且将导致 1-树。

无向版本使用树中的基本循环的概念来找到替代边，我试图在 NetworkX 库中使用 find_cycle() 函数来实现这一想法。我手动执行了上升法的第一次迭代，发现我计算的所有可能的 $\epsilon$ 值与程序找到的值不匹配。我发现了一些程序遗漏的值，它也发现了一些我遗漏的值。对于示例图，我发现以下边对是替代边，其中第一条边不在 1-树中，第二条边是在 1-树中可以被替换的边，使用下面的最小 1-树。

$$ \begin{array}{l} (0, 1) \rightarrow (2, 1) \text{ 有效: } \epsilon = 56 \\\ (0, 2) \rightarrow (4, 2) \text{ 有效: } \epsilon = 25 \\\ (0, 3) \rightarrow (2, 3) \text{ 有效: } \epsilon = 52 \\\ (0, 5) \rightarrow (1, 5) \text{ 有效: } \epsilon = \frac{30 - 52}{0 - 0} \text{，无效} \\\ (1, 3) \rightarrow (2, 3) \text{ 有效: } \epsilon = 15.5 \\\ (2, 5) \rightarrow (1, 5) \text{ 有效: } \epsilon = 5.5 \\\ (3, 1) \rightarrow (2, 1) \text{ 有效: } \epsilon = 5.5 \\\ (3, 5) \rightarrow (1, 5) \text{ 有效: } \epsilon = \frac{30 - 46}{-1 + 1} \text{，无效} \\\ (4, 1) \rightarrow (2, 1) \text{ 有效: } \epsilon = \frac{41 - 90}{1 - 1} \text{，无效} \\\ (4, 3) \rightarrow (2, 3) \text{ 有效: } \epsilon = \frac{30 - 95}{1 - 1} \text{，无效} \\\ (4, 5) \rightarrow (1, 5) \text{ 有效: } \epsilon = -25.5 \text{，无效（负的 }\epsilon) \\\ (5, 3) \rightarrow (2, 3) \text{ 有效: } \epsilon = 25 \\\ \end{array} $$

我遗漏了程序找到的以下替代边。

$$ \begin{array}{l} (1, 0) \rightarrow (4, 0) \text{ 有效: } \epsilon = 80 \\\ (1, 4) \rightarrow (0, 4) \text{ 有效: } \epsilon = 73 \\\ (2, 0) \rightarrow (4, 0) \text{ 有效: } \epsilon = \frac{17 - 60}{1 - 1} \text{，无效} \\\ (2, 4) \rightarrow (0, 4) \text{ 有效: } \epsilon = -18 \text{，无效（负的 }\epsilon) \\\ (3, 0) \rightarrow (4, 0) \text{ 有效: } \epsilon = 28 \\\ (3, 4) \rightarrow (0, 4) \text{ 有效: } \epsilon = 78 \\\ (5, 0) \rightarrow (4, 0) \text{ 有效: } \epsilon = 35 \\\ (5, 4) \rightarrow (0, 4) \text{ 有效: } \epsilon = \frac{17 - 81}{0 - 0} \text{，无效} \\\ \end{array} $$

请注意，如果我们在上升方向上移动，有些替代边不会交叉，这些边对的分母为零。此外，$\epsilon$ 是距离，负距离的概念没有意义。将负距离解释为相反方向上的正距离，如果我们需要朝那个方向移动，则上升方向向量将指向相反的方向。

我的列表与程序列表不匹配的原因是，find_cycle() 并不总是返回包含新边的基本循环。如果我在图中另一个循环中的顶点（在本例中为 ${(0, 4), (4, 0)}$）上调用 find_cycle()，它会返回该循环而不是真正的基本循环。

这促使我思考，是什么真正决定了 1-树形图中的边是否可以互相替代。在所有有效的替代情况下，这两条边都指向同一个顶点。如果它们没有指向同一个顶点，那么树形图的度数约束就会被破坏，因为我们没有用指向同一个节点的另一条边替换指向该节点的边。这无论边是否属于同一个基本循环都是成立的。

因此，find_epsilon() 现在会遍历图中的所有边，但不包括所选的 1-树形图 $k \in K(\pi, d)$ 中的边，并找到 $k$ 中指向同一个顶点的另一条边，交换它们，然后检查度数约束是否未被破坏，它是否有正确数量的边，以及它是否仍然是连通的。这是一个更有效的方法，它也找到了更多有效的替换，因此我希望能最终将返回的解降至最优解，也许是因为它在其中一次迭代中甚至只缺少了 $\epsilon$ 的正确值。

它没有。

下一步#

此时，我没有任何明确的解决方案，但有两个不令人满意的选择。

• 我通过手动执行上升法的第一次迭代找到了 find_epsilon() 的问题。这花了大约 90 分钟。我可以尝试继续这个过程，并希望在迭代 1 正确执行的同时找到代码中的其他错误，但我怀疑我永远无法达到程序找到错误解所需的 9 次迭代。
• 继续进行 Held-Karp 松弛的枝剪界定部分。我希望，由于 Held 和 Karp 对枝剪界定方法进行了完整的执行，我能够利用它来追踪松弛的完整执行，并通过这种方式找到上升法中的缺陷。

我很快就会与我的 GSoC 指导老师讨论下一步行动。

参考文献#

Held, M., Karp, R.M. The traveling-salesman problem and minimum spanning trees. Operations research, 1970-11-01, Vol.18 (6), p.1138-1162. https://www.jstor.org/stable/169411