Google Summer of Code Logo with NetworkX logo

欢迎!这篇文章不会讨论我的 Google Summer of Code 项目的技术实现细节或理论工作,而是作为我今年夏天所做工作的总结和回顾。

我对能够完成的工作感到非常高兴,并且相信我成功地完成了我的项目。

概述#

我的项目标题为 NetworkX:实现 Asadpour 非对称旅行商问题算法。在 Summer of Code 项目项目页面上给出的更新摘要如下。

该项目似乎实现了由 Asadpour 等人开发的非对称旅行商问题,最初发表于 2010 年,并在 2017 年进行了修订。该项目分为多种方法,每种方法在项目期间都有一个设定的时间表。我们首先使用 Held 和 Karp 原论文中的 Ascent 方法求解 Held-Karp 松弛。假设结果是分数的,我们继续进入 Asadpour 算法(整数解根据定义是最优的,并立即返回)。我们使用最大熵舍入方法近似 Held Karp 解的无向支撑上的生成树分布,以构建树的分布。粗略地说,采样任何给定树的概率与所有边 lambda 值的乘积成正比。我们使用 V. G. Kulkarni 开发的迭代方法从分布中采样 2 log n 棵树,并在将方向返回到弧后选择成本最小的树。最后,使用最小网络流算法增强最小树,并将其缩短到最小哈密顿循环的 O(log n / log log n) 近似值。

我 2021 年夏季代码的提案 PDF 可以在这里找到

我对 NetworkX 的所有更改和添加都是此拉取请求的一部分,也可以在我的 GitHub 存储库的分支此分支上找到,但我将在稍后更详细地讨论更改和提交。另请注意,对于我在每个部分中列出的提交,这只是一个不完整的列表,仅涉及对该函数或其测试的集中提交。有关完整列表,请参阅拉取请求或我 NetworkX 分支上的 bothTSP GitHub 分支。

我今年夏天对 NetworkX 的贡献主要包括以下函数和类,我将在本文的各个部分中讨论每个函数和类。面向用户的函数和类也链接到下面列表中的 NetworkX 的开发者文档以及它们的部分标题。

这些函数也进行了单元测试,一旦拉取请求合并,这些测试将集成到 NetworkX 中。

参考文献#

以下论文是所有这些算法的起源,当然它们对项目的完成至关重要。

[1] A. Asadpour、M. X. Goemans、A. Madry、S. O. Gharan 和 A. Saberi,非对称旅行商问题的 O (log n / log log n) 近似算法,SODA '10,工业与应用数学学会,2010 年,第 379 - 389 页https://dl.acm.org/doi/abs/10.5555/1873601.1873633

[2] J. Edmonds,最优分支,国家标准局研究杂志,1967 年,第 71B 卷,第 233-240 页,https://archive.org/details/jresv71Bn4p233

[3] M. Held、R.M. Karp,旅行商问题和最小生成树。运筹学,1970 年 11 月 1 日,第 18 卷(6),第 1138-1162 页。https://www.jstor.org/stable/169411

[4] G.K. Janssens、K. Sörensen,按成本递增顺序生成所有生成树的算法,Pesquisa Operacional,2005 年 8 月,第 25 卷(2),第 219-229 页,https://www.scielo.br/j/pope/a/XHswBwRwJyrfL88dmMwYNWp/?lang=en

[5] V. G. Kulkarni,生成随机组合对象,算法杂志,11(1990),第 185-207 页。

SpanningTreeIterator#

SpanningTreeIterator 是我作为 GSoC 项目完成的第一个贡献。此类接受一个图并按成本递增顺序返回其中的每个生成树,这使其成为 [4] 的直接实现。

这个迭代器有趣的一点是它不用于 Asadpour 算法的一部分,而是作为中间步骤,以便我可以开发 Held Karp 松弛所需的 ArborescenceIterator。它的工作原理是将图的边划分为包含、排除或打开,然后找到尊重图边上分区数据的最小生成树。为了使其工作,我创建了一个名为 kruskal_mst_edges_partition 的新最小生成树函数,它正是这样做的。为了防止冗余,所有 Kruskal 最小生成树现在都使用此函数(原始 kruskal_mst_edges 函数现在只是分隔版本的包装器)。一旦从迭代器返回一个生成树,该树的分区数据就会被拆分,以便新生成的区间的并集是分区中除返回的最小生成树之外的所有生成树的集合。

正如我之前提到的,SpanningTreeIterator 在我的 GSoC 项目中没有直接使用,但我仍然决定实现它以了解分区过程并能够直接使用 [4] 中的示例,然后再继续使用 ArborescenceIterator。我相信这个类对 NetworkX 的其他用户将非常有用,并且为构建 ArborescenceIterator 提供了坚实的基础。

关于 SpanningTreeIterator 的博客文章

2021 年 6 月 5 日 - 查找所有最小树

2021 年 6 月 10 日 - 实现迭代器

关于 SpanningTreeIterator 的提交

现在,在这个项目开始时,我的提交信息不是很好……在我不小心提交到错误的分支后,我遇到了一些合并冲突的问题,这是我第一次使用预提交钩子。

我没有在这里更改提交信息,以便您可以被我毫无帮助的消息所迷惑,但对其进行了注释以提供更准确的提交描述。

测试 - 重写了 Kruskal 算法以尊重分区并在单独的文件中存根迭代器时对其进行了测试

我不确定提交钩子是如何工作的…… - 添加了测试用例并完成了 Spanning Tree Iterator 在错误文件中的实现

将迭代器移动到正确的文件以维护正确的代码库可见性 - 意识到迭代器分别需要在 mst.pybranchings.py 中,以保持私有函数隐藏

迭代器的文档更新 - 无需解释

更新 mst.py 以接受建议 - 接受了代码审查中的文档字符串编辑

来自 dshult 的审查建议 - 实施了我的导师之一的代码审查建议

清理代码,如果可能,合并函数并向所有人开放分区功能

实施来自 boothby 的建议

ArborescenceIterator#

ArborescenceIterator 是 [4] 中讨论的算法的修改版本,以便它迭代生成树。

这个迭代器实现起来稍微困难一些,但这是由于最小生成树算法的结构方式,而不是分区方案不适用于有向图。事实上,分区方案与无向 SpanningTreeIterator 相同,但 Edmonds 算法更复杂,并且关于节点如何收缩以及尊重分区数据意味着什么的几个边缘情况。为了充分理解 NetworkX 实现,我不得不阅读原始的 Edmonds 论文 [2]。

最显著的变化是,当迭代器在执行 Edmonds 算法之前将下一个分区写入图的边缘时,如果任何传入边被标记为包含,则所有其他边都被标记为排除。这是 SpanningTreeIterator 的一个隐式部分,但需要在这里显式地执行,以便如果相关顶点在 Edmonds 算法期间被合并,我们不能在合并被反转后选择两条到同一顶点的传入边。

最后需要说明的是,ArborescenceIteratorSpanningTreeIterator 多一个初始参数,即能够为它提供一个初始分区并迭代所有成本大于初始分区的生成树。这在分支定界方法中被使用,但不再是我 Asadpour 算法实现的一部分。

关于 ArborescenceIterator 的博文

2021 年 6 月 5 日 - 查找所有最小树

2021 年 6 月 10 日 - 实现迭代器

关于 ArborescenceIterator 的提交

此处列出的我的提交仍然有注释,并且大部分工作是在同一时间完成的。

测试 - 重写了 Kruskal 算法以尊重分区并在单独的文件中存根迭代器时对其进行了测试

将迭代器移动到正确的文件以维护正确的代码库可见性 - 意识到迭代器分别需要在 mst.pybranchings.py 中,以保持私有函数隐藏

包含 Black 格式化 - 修改了 Edmonds 算法以尊重分区

修改了 ArborescenceIterator 以接受初始分区 - 无需解释

迭代器的文档更新 - 无需解释

更新 branchings.py 接受文档字符串编辑 - 无需解释

来自 dshult 的审查建议 - 实施了我的导师之一的代码审查建议

清理代码,如果可能,合并函数并向所有人开放分区功能

实现了来自 rossbar 的审查建议

实施来自 boothby 的建议

held_karp_ascent#

Held Karp 松弛是我 GSoC 项目中最困难的部分,也是我今年 5 月份最担心的一部分。

我关于如何解决松弛的计划在整个夏天也在不断发展,最终形成了 held_karp_ascent。在我的 GSoC 提案中,我讨论了使用 scipy 来解决松弛,但 Held Karp 松弛是一个半无限线性问题(即,它是有限的但指数级的),因此我很快就会超过几乎任何运行代码的计算机的能力。幸运的是,我意识到这一点是在我还在撰写提案的时候,并且能够进行修改。接下来,我想使用椭球算法,因为这是 Asadpour 论文 [1] 中建议的方法。

碰巧的是,椭球算法没有在 numpyscipy 中实现,在讨论将该算法作为本项目的一部分实现的可行性后,我们决定一个健壮的椭球求解器本身就是一个 GSoC 项目,超出了 Asadpour 算法的范围。需要另一种方法,并且找到了。在 Held 和 Karp [3] 的原始论文中,他们提出了三种不同的算法来解决松弛问题,即列生成技术、上升法和分支定界法。在阅读了论文并比较了所有方法后,我决定分支定界法在性能方面是最好的,并希望实现它。

分支定界法是上升法的一个修改版本,所以我首先实现了上升法,然后在其基础上实现了分支定界法。这额外的好处是允许我比较两者并确定哪个实际上更好。

实现上升法证明是困难的。在寻找最小 1-树和寻找 epsilon 的值时,存在一些细微的错误,因为没有意识到图中所有有效的边替换。有关这些问题的更多信息,请参阅我题为“理解上升法”的文章。即使在之后,上升法也没有正常工作,但我决定转向分支定界法,希望能更多地了解这个过程,以便我可以修复上升法。

这正是发生的事情!在调试分支定界法时,我意识到我用于查找最小 1-树集的函数会过早停止搜索,并可能错过最小 1-树。一旦我修复了这个错误,上升法和分支定界法都开始产生正确的结果。

但是哪个会在最终项目中使用呢?

好吧,这取决于哪个输出更兼容 Asadpour 算法的其余部分。上升法可以找到一个分数解,其中边不是完全在解中或不在解中,而分支定界法将花费时间来确保解是整数解。碰巧的是,Asadpour 算法期望 Held Karp 松弛的分数解,因此最终上升法胜出,分支定界法从项目中删除。

所有这些都详细记录在我为这个主题撰写的(许多)博文中,这些博文列在下面。

关于 Held Karp 松弛的博文

我的前两篇文章是关于 scipy 解和椭球算法的。

2021 年 4 月 11 日 - Held Karp 松弛

2021 年 5 月 8 日 - Held Karp 分离预言机

下一篇文章讨论了 Held 和 Karp 原始论文 [3] 中提出的每种算法的优缺点。

2021 年 6 月 3 日 - 深入了解 Held Karp

最后,最后三篇与 Held Karp 相关的文章是关于我所实现的算法的调试。

2021 年 6 月 22 日 - 理解上升法

2021 年 6 月 28 日 - 实现 Held Karp 松弛

2021 年 7 月 7 日 - 完成 Held Karp

关于 Held Karp 松弛的提交

仅在需要时提供注释。

获取 Black 格式化 - 初始上升法实现

正在调试上升法并进行 Black 格式化

上升法终止,但在非最优解处终止

少量编辑 - 删除了一些调试语句

修复了终止条件,但仍然给出非最优结果

小错误修复,仍然是非最优结果 - 确保报告的答案是多个选项中的循环

修复了 find_epsilon() 中的细微错误 - 修复了不正确的替换检测错误

清理代码并尝试了一些不起作用的东西

Black 格式化 - 初始分支定界实现

分支定界返回最优解

black 格式化更改 - 将上升法和分支定界法拆分为不同的函数

性能调整和测试分数答案

修复了测试错误,希望如此

上升法的 Asadpour 输出

删除了分支定界法。一个单元测试行为异常

为上升法添加了非对称分数测试

删除了 printn 语句并调整了最终测试使其更加非对称

更改 HK 以仅报告答案的支持

文档更新

spanning_tree_distribution#

一旦我们有了 Held Karp 松弛的支持,我们就计算支持的边权重 $\gamma$,以便任何树被采样的概率与 $e^{\gamma}$ 在其边上的乘积成正比。这在 Asadpour 论文中称为最大熵分布。此过程包含在 Asadpour 论文 [1] 第 386 页中。

  1. 设置 $\gamma = \vec{0}$。
  2. 当存在一条边 $e$ 使得 $q_e(\gamma) > (1 + \epsilon)z_e$ 时
  • 计算 $\delta$,使得如果我们定义 $\gamma’$ 为 $\gamma_e’ = \gamma_e - \delta$ 且 $\gamma_f’ = \gamma_e$ 对于所有 $f \in E \backslash {e}$,则 $q_e(\gamma’) = (1 + \epsilon / 2)z_e$
  • 设置 $\gamma \leftarrow \gamma'$
  1. 输出 $\tilde{\gamma} := \gamma$。

其中 $q_e(\gamma)$ 是任何给定边 $e$ 出现在以与 $\exp(\gamma(T))$ 成正比的概率选择的采样生成树中的概率。$\delta$ 也给出为

$$ \delta = \frac{q_e(\gamma)(1-(1+\epsilon/2)z_e)}{(1-q_e(\gamma))(1+\epsilon/2)z_e} $$

因此,Asadpour 论文几乎完成了此函数的所有繁重工作。然而,他们并没有非常清楚如何计算 $q_e(\gamma)$,除了可以使用基尔霍夫树矩阵定理之外。

我最初计算 $q_e(\gamma)$ 的方法是将基尔霍夫定理应用于原始拉普拉斯矩阵和从图中收缩边 $e$ 后产生的拉普拉斯矩阵。测试很快表明,一旦从图中收缩了边,它就不会影响拉普拉斯矩阵的值,因此在减去 $\delta$ 后,该边的概率会增加而不是减少。将我的原始 $q_e(\gamma)$ 值乘以 $\exp(\gamma_e)$ 被证明是这里的解决方案,原因在我的博文“熵分布”以及特别是“更新!(2021 年 7 月 28 日)”部分中进行了广泛的讨论。

关于 spanning_tree_distribution 的博文

2021 年 7 月 13 日 - 熵分布设置

2021 年 7 月 20 日 - 熵分布

关于 spanning_tree_distribution 的提交

spanning_tree_distribution 的草稿

更改 HK 以仅报告答案的支持 - 需要将 $\gamma$ 限制在 Held Karp 松弛的支持范围内是导致此更改的原因

通过更改为 MultiGraph 修复了收缩错误。概率 > 1 的问题 - 因为概率仅与边权重的乘积成正比,所以这实际上不是问题

Black 格式化 - 重写了测试并清理了代码

修复了 pypi 测试错误 - pypi 测试没有 numpyscipy,我忘记将测试标记为如果它们不可用则跳过

进一步测试 dist 修复 - 修复了将 $q_e(\gamma)$ 乘以 $\exp(\gamma_e)$ 的函数,并在 $\delta$ 出现异常行为时实现了异常

可以采样生成树 - 使用新的辅助函数简化了查找 $q_e(\gamma)$ 的过程

文档更新

来自 dshult 的审查建议 - 实施了我的导师之一的代码审查建议

实施来自 boothby 的建议

sample_spanning_tree#

如果我们无法从中采样,那么生成树分布有什么用呢?

虽然 Asadpour 论文 [1] 对采样过程进行了粗略的概述,但其大部分方法都来自 Kulkarni 论文,即 *生成随机组合对象* [5]。该论文给出了更详细的解释,甚至提供了第 202 页的伪代码。

$U = \emptyset,$ $V = E$
执行 $i = 1$ 到 $N$;
$\qquad$令 $a = n(G(U, V))$
$\qquad\qquad a’$ $= n(G(U \cup {i}, V))$
$\qquad$生成 $Z \sim U[0, 1]$
$\qquad$如果 $Z \leq \alpha_i \times \left(a’ / a\right)$
$\qquad\qquad$则 $U = U \cup {i}$,
$\qquad\qquad$否则 $V = V - {i}$
$\qquad$结束。
停止。$U$ 是所需的生成树。

这里唯一的真正困难是跟踪节点是如何收缩的。我的第一次尝试是一团糟的 if 语句等等,但切换到合并查找数据结构(或不相交集数据结构)被证明是一个明智的决定。

当然,能够采样生成树是一回事,而知道采样技术是否与预期分布相符则是另一回事。我针对 sample_spanning_tree 的第一次迭代测试只是采样了大量树木(50000),并打印了与生成树归一化分布的百分比误差。在 50000 的样本量下,所有误差都在 10% 以内,但我仍然希望找到一个更好的测试。

从高中时期的 AP 统计学课程中,我记住了 $X^2$(卡方)检验,并意识到它在这里非常完美。scipy 甚至能够进行卡方检验。通过转换为卡方检验,我能够将样本量减少到 1200(接近进行有效卡方检验所需的最小样本量),并在 $\alpha = 0.01$ 的显著性水平下使用适当的假设检验。不幸的是,在为 sample_spanning_tree 添加 @py_random_state 装饰器之前,该测试仍然会在 1% 的时间内失败,然后测试可以在 Random 对象中传递以生成可重复的结果。

关于 sample_spanning_tree 的博文

2021 年 7 月 21 日 - 采样生成树的准备工作

2021 年 7 月 28 日 - 采样生成树

关于 sample_spanning_tree 的提交

可以采样生成树

开发采样生成树的测试

sample_spanning_tree 测试更改为卡方检验

添加测试用例 - *实现了 @py_random_state 装饰器*

文档更新

来自 dshult 的审查建议 - 实施了我的导师之一的代码审查建议

asadpour_atsp#

此函数是拼图的最后一块,将所有其他部分连接在一起并产生最终结果!

此函数的实现实际上相当顺利。我遇到的唯一技术困难是读取 flow_dict 的支持,而理论上的困难是调整 min_cost_flow 函数以解决最小循环问题。哦,如果流量大于 1,我需要向图中添加平行边,以便它仍然是欧拉图。

下面简要概述了整个算法

  1. 求解 Held Karp 松弛并使其对称以使其无向。
  2. 计算 Held Karp 支撑图上的最大熵生成树分布。
  3. 采样 $2 \lceil \ln n \rceil$ 生成树,并在重新引入边的方向之前记录最小的权重树。
  4. 找到最小成本循环以创建包含采样树的欧拉图。
  5. 获取该图的欧拉路径并对其进行捷径处理。
  6. 返回捷径答案。

关于 asadpour_atsp 的博文

2021 年 7 月 29 日 - 纵观全局

2021 年 8 月 10 日 - 完成 Asadpour 算法

关于 asadpour_atsp 的提交

asadpour_tsp 的未经测试的实现

修复了读取 flow_dict 的问题

修复了 asadpour_tsp 中的运行时错误 - *一般旅行商问题函数假设图是无向的。这与 atsp 算法不兼容*

black 格式化 - *修复了来自流支持错误的平行边*

添加测试用例

文档更新

一个新的测试和检查

修复了测试中的舍入错误

来自 dshult 的审查建议 - 实施了我的导师之一的代码审查建议

实现了来自 rossbar 的审查建议

未来参与 NetworkX#

总的来说,我真的很享受这次暑期代码编写活动。我能够扩展分支,继续学习 Python 和更多关于图和图算法的知识,这是一个我感兴趣的领域。

假设我在即将到来的秋季学期有任何空闲时间,我都很乐意继续参与 NetworkX。事实上,即使我目前的代码可以正常工作,我心中也有一些想法。

  • sample_spanning_tree 移动到 mst.py 并将其重命名为 random_spanning_tree。采样随机生成树的功能不是 NetworkX 库的主要部分,但可能对其他人有用。我的导师之一提到它与 斯坦纳树 相关,如果我能帮助其他开发人员和用户,我就会去做。

  • 调整 sample_spanning_tree,使其可以使用加性和乘性权重函数。Asadpour 算法只需要乘性权重,但 Kulkarni 论文 [5] 确实讨论了使用加性权重函数,这可能对其他 NetworkX 用户更有用。

  • 将我的克里奇霍夫树矩阵定理辅助函数移动到 laplacian_matrix.py,以便其他 NetworkX 用户可以访问它。

  • 研究以下关于 Held Karp 松弛的文章。虽然我没有这方面的确凿证据,但我确实认为 Held Karp 松弛是我实现 Asadpour 算法中最慢的部分,因此是改进它的最佳位置。我正在使用的上升方法来自最初的 Held 和 Karp 论文 [3],但他们确实发布了第二部分,其中可能包含更好的算法。引用如下。

    M. Held, R.M. Karp,*旅行商问题和最小生成树:第二部分*。数学规划,1971,1(1),第 6-25 页。https://doi.org/10.1007/BF01584070

  • 重构 branchings.py 中的 Edmonds 类。该类是 Edmonds 分支算法的实现,但使用迭代方法而不是 Edmonds 论文 [2] 中讨论的递归方法。我还同意与另一个人 lkora 合作,帮助重做此类,并可能添加一个 minimum_maximal_branching 函数来查找仍然连接尽可能多节点的最小分支。这类似于无向图中的生成森林。目前,我们两人都没有时间开始此类工作。有关更多信息,请参考问题 #4836

虽然我可以在这个问题的某些方面进行改进,但重要的是要记住,这个项目仍然是完全成功的。NetworkX 现在有一个算法可以近似非对称或有向图中的旅行商问题。