实现 spanning_tree_distribution 证明存在一些 NetworkX 问题和一个算法问题。回想一下，创建分布的算法在 Asadpour 论文中给出如下所示
- 设置 $\gamma = \vec{0}$。
- 当存在一个边 $e$ 使得 $q_e(\gamma) > (1 + \epsilon) z_e$ 时
- 计算 $\delta$，使得如果我们定义 $\gamma’$ 为 $\gamma_e’ = \gamma_e - \delta$，并且对于所有 $f \in E\ \backslash {e}$，$\gamma_f’ = \gamma_f$，那么 $q_e(\gamma’) = (1 + \epsilon/2)z_e$。
- 设置 $\gamma \leftarrow \gamma’$。
- 输出 $\tilde{\gamma} := \gamma$。
现在，我在上一篇题为 熵分布设置 的博客中提出的流程在 while 循环部分运行良好。我在 NetworkX API 中遇到的所有困难都发生在 q 内部函数中。
在我编写完函数后，我当然需要运行它，最初我只是打印 gamma 字典，以便查看每条边的值。我的第一个测试使用对称分数 Held Karp 解，令我惊讶的是，所有 $\gamma$ 的值都返回为 0。我不认为这是预期的行为，因为如果是这样，那么在整个 Asadpour 算法中就没有理由包含这一步，因此我开始使用 PyCharm 的调试器来深入研究代码。结果正如我所料，并不正确。我正在对原始图运行 Krichhoff 树矩阵定理，因此返回的概率比我比较的 $z_e$ 值小一个数量级。此外，所有值都相同，所以我知道这是一个问题，而不是我检查的第一条边的概率异常小。
所以，我回到了 Asadpour 论文，并开始问自己一些问题，比如
- 我是否需要以某种方式对 Held Karp 答案进行归一化？
- 我是否需要考虑 $E$（Held Karp 松弛解的无向支撑）之外的边，或者只处理 $E$ 中的边？
第一个问题很容易被驳回，如果需要归一化，Asadpour 论文中会提到，并且没有对如何进行归一化的描述，我找到“正确”归一化方法的可能性几乎为零。第二个问题确实需要一些深入研究。Asadpour 论文中关于使用 Krichhoff 定理的部分都讨论了在图 $G$ 中使用它，这就是我最初使用 $G$ 中所有边而不是 $E$ 中的边的原因。一些提示表明我只需要考虑 $E$ 中的边，第一个是算法概述，它指出
查找权重 ${\tilde{\gamma}}_{e \in E}$
特别是 $e \in E$ 语句表明我不需要考虑不在 $E$ 中的边。其次，引理 7.2 以以下语句开头
设 $G = (V, E)$ 为一个图，其中对于 $e \in E$，边权为 $\gamma_e$。
根据函数的当前状态和这些提示，我决定将输入图简化为 spanning_tree_distribution，只保留 $z_e > 0$ 的边。现在对对称分数解进行测试，它仍然返回 $\gamma = \vec{0}$，但在第一次迭代中，它比较的概率要接近得多。由于我没有示例图和分布可以用来进行测试，这可能是正确的答案，但所有值都相同这一事实仍然让我感到困惑。
我的下一步是确定在 $\gamma = \vec{0}$ 的第一次迭代中，边出现在生成树中的实际概率。这可以使用我的 SpanningTreeIterator 很容易地完成，并利用了 $\gamma = \vec{0} \equiv \lambda_e = 1\ \forall\ e \in \gamma$ 的事实，因此我们可以遍历生成树并统计每条边出现的次数。
该脚本列在下面
import networkx as nx
edges = [
(0, 1),
(0, 2),
(0, 5),
(1, 2),
(1, 4),
(2, 3),
(3, 4),
(3, 5),
(4, 5),
]
G = nx.from_edgelist(edges, create_using=nx.Graph)
edge_frequency = {}
sp_count = 0
for tree in nx.SpanningTreeIterator(G):
sp_count += 1
for e in tree.edges:
if e in edge_frequency:
edge_frequency[e] += 1
else:
edge_frequency[e] = 1
for u, v in edge_frequency:
print(
f"({u}, {v}): {edge_frequency[(u, v)]} / {sp_count} = {edge_frequency[(u, v)] / sp_count}"
)此输出表明 q 返回的概率应在边之间有所不同，并且 $\gamma$ 的正确解当然不是 $\vec{0}$。
(networkx-dev) mjs@mjs-ubuntu:~/Workspace$ python3 spanning_tree_frequency.py
(0, 1): 40 / 75 = 0.5333333333333333
(0, 2): 40 / 75 = 0.5333333333333333
(0, 5): 45 / 75 = 0.6
(1, 4): 45 / 75 = 0.6
(2, 3): 45 / 75 = 0.6
(1, 2): 40 / 75 = 0.5333333333333333
(5, 3): 40 / 75 = 0.5333333333333333
(5, 4): 40 / 75 = 0.5333333333333333
(4, 3): 40 / 75 = 0.5333333333333333
让我们关注第一条边，$(0, 1)$。我的暴力脚本表明它出现在以下图中 75 棵生成树中的 40 棵中，其中每条边都标记了其 $z_e$ 值。
然而，q 却说这条边出现在 75 棵生成树中的 24 棵中。由于分母是正确的，我决定专注于分子，即 $G\ \backslash\ \{(0, 1)\}$ 中的生成树数量。该图如下所示。
可以说，此图应在顶点 0 上有一个自环，但这不会影响拉普拉斯矩阵，因此此处省略。基本上，邻接矩阵的 $[0, 0]$ 元素将为 1，顶点 0 的度数将为 5，而 $5 - 1 = 4$ 将是该元素在没有自环的情况下所具有的值。
发生的情况是，我向 nx.contracted_edge 提供了一个 Graph 类图（不是有向图，因为 $E$ 是无向的），并且得到了一个 Graph 类的图。Graph 类不支持两个节点之间的多条边，因此返回的图在节点 0 和节点 2 之间只有一条边，这会影响整体拉普拉斯矩阵，从而影响生成树的数量。从 Graph 切换到 MultiGraph 解决了这个问题，但这种细微的更改应在 NetworkX 文档中提到，该文档链接 在此。我绝对相信，如果收缩一条边，输出应该自动包含两条 $(0, 2)$ 边。可以说，可以更改默认行为以匹配这一点，但至少文档应解释此问题。
现在，q 函数返回了 $(0, 1)$ 的正确答案 $40 / 75$，以及所有其他边的正确值，只要所有 $\gamma_e$ 的值为 0。但是当我尝试计算 $\delta$ 时，测试以 ValueError 错误退出。q 返回了一条边出现在采样生成树中的概率大于 1，这显然是不可能的，但也导致 $\delta$ 的分母变为负数，违反了自然对数的定义域。
在我调查此问题期间，我注意到，在计算 $\delta$ 并将其从 $\gamma_e$ 中减去后，它没有对 $q_e$ 产生预期效果。回想一下，我们定义 $\delta$，使得 $\gamma_e - \delta$ 产生 $(1 + \epsilon / 2) z_e$ 的 $q_e$。换句话说，$\delta$ 的作用是降低过高的边概率，但在我的当前实现中，它产生了相反的效果。$q_{(0, 1)}$ 的值从 0.5333 变为略大于 0.6。如果让这种趋势继续下去，程序最终将遇到 $q_e \geq 1$ 的情况之一，并导致程序崩溃。
这里我可以使用边 $(0, 1)$ 作为示例来展示问题。具有 $\gamma = \vec{0}$ 的 $G$ 的原始拉普拉斯矩阵为
$$ \begin{bmatrix} 3 & -1 & -1 & 0 & 0 & -1 \\\ -1 & 3 & -1 & 0 & -1 & 0 \\\ -1 & -1 & 3 & -1 & 0 & 0 \\\ 0 & 0 & -1 & 3 & -1 & -1 \\\ 0 & -1 & 0 & -1 & 3 & -1 \\\ -1 & 0 & 0 & -1 & -1 & 3 \end{bmatrix} $$
而 $G\ \backslash\ \{(0, 1)\}$ 的拉普拉斯矩阵为
$$ \begin{bmatrix} 4 & -2 & -1 & -1 & 0 \\\ -2 & 3 & 0 & 0 & -1 \\\ -1 & 0 & 3 & -1 & -1 \\\ -1 & 0 & -1 & 3 & -1 \\\ 0 & -1 & -1 & -1 & 3 \end{bmatrix} $$
第一个余子式的行列式是 $40 / 75$ 的来源。现在，考虑我们第一次更新 $\gamma_{(0, 1)}$ 后的拉普拉斯矩阵。$G$ 的拉普拉斯矩阵变为
$$ \begin{bmatrix} 2.74 & -0.74 & -1 & 0 & 0 & -1 \\\ -0.74 & 2.74 & -1 & 0 & -1 & 0 \\\ -1 & -1 & 3 & -1 & 0 & 0 \\\ 0 & 0 & -1 & 3 & -1 & -1 \\\ 0 & -1 & 0 & -1 & 3 & -1 \\\ -1 & 0 & 0 & -1 & -1 & 3 \end{bmatrix} $$
其第一个余子式行列式从 75 减少到 61.6。我们期望 $G\ \backslash\ \{(0, 1)\}$ 的矩阵值为多少？好吧，我们知道 $q_e$ 的最终值需要为 $(1 + \epsilon / 2) z_e$ 或 $1.1 \times 0.41\overline{6}$，即 $0.458\overline{3}$。所以
$$ \begin{array}{r c l} \displaystyle\frac{x}{61.6} &=& 0.458\overline{3} \\\ x &=& 28.2\overline{3} \end{array} $$
并且第一个余子式行列式的值应为 $28.2\overline{3}$。但是，更新 $\gamma_e$ 值后，$(0, 1)$ 的收缩拉普拉斯矩阵为
$$ \begin{bmatrix} 4 & -2 & -1 & -1 & 0 \\\ -2 & 3 & 0 & 0 & -1 \\\ -1 & 0 & 3 & -1 & -1 \\\ -1 & 0 & -1 & 3 & -1 \\\ 0 & -1 & -1 & -1 & 3 \end{bmatrix} $$
与之前相同！ 唯一一个与之前不同的 $\gamma_e$ 的边是 (0, 1)，但由于它是收缩边，所以它不再在图中，因此不会影响第一个余子式的行列式值！
但是，如果我们将算法改为对 $\gamma_e$ 加 $\delta$ 而不是减它，那么 $G\ \backslash\ \{e\}$ 的拉普拉斯算子的第一个余子式的行列式不会改变，但 $G$ 的第一个余子式的拉普拉斯算子的行列式会增加。这会降低在生成树中选取 $e$ 的整体概率。而且，如果我们恰好对 (0, 1) 的例子使用与之前相同的 $\delta$ 公式，那么 $q_{(0, 1)}$ 会变成 0.449307。回想一下我们目标值为 0.458$\overline{3}$。这个答案的误差为 -1.96%。
$$ \begin{array}{r c l} \text{误差} &=& \frac{0.449307 - 0.458333}{0.458333} \times 100 \\\ &=& \frac{-0.009026}{0.458333} \times 100 \\\ &=& -0.019693 \times 100 \\\ &=& -1.9693% \end{array} $$
此外，测试现在可以顺利完成。
更新！ (2021 年 7 月 28 日)#
进一步的研究和与导师的讨论揭示了我最初的分析存在多么大的缺陷。在下一步，对生成树进行采样时，对 $\gamma$ 添加任何值都会直接增加采样该边的概率。也就是说，我发现的原始问题仍然存在。
回到我们在每个生成树映射到包含所需边的所有生成树的图上的概念，这仍然是让我们使用基尔霍夫树矩阵定理的关键思想。而且，收缩边仍然会得到一个图，在这个图中每个生成树都可以映射到包含 $e$ 的对应生成树。然而，$G \backslash \{e\}$ 中这些生成树的权重在两个图之间并没有完全映射。
回想一下，我们正在处理一个乘法权重函数，所以树的最终权重是其所有边上的 $\lambda$ 的乘积。
$$ c(T) = \prod_{e \in E} \lambda_e $$
上述语句可以扩展为
$$ c(T) = \lambda_1 \times \lambda_2 \times \dots \times \lambda_{|E|} $$
其中边 $1, 2, \dots |E|$ 的顺序是任意的。由于边的顺序是任意的，并且由于乘法的结合律，我们可以假设在不失一般性的前提下，所需边 $e$ 是序列中的最后一条边。
$G \backslash \{e\}$ 中的任何生成树都不能包含最后的 $\lambda$，因为该边在图中不存在。因此，为了将权重从 $G \backslash \{e\}$ 中的树转换为权重，我们需要将 $\lambda_e$ 乘回收缩树的权重。所以，我们现在可以声明
$$ c(T \in \mathcal{T}: T \ni e) = \lambda_e \prod_{f \in E} \lambda_f\ \forall\ T \in G \backslash \{e\} $$
或者对于 $G \backslash \{e\}$ 中的所有树，$G$ 中对应树的成本是其边 $\lambda$ 的乘积乘以所需边的权重。现在回想一下 $q_e(\gamma)$ 是
$$ \frac{\sum_{T \ni e} \exp(\gamma(T))}{\sum_{T \in \mathcal{T}} \exp(\gamma(T))} $$
特别地，我们正在处理上述分数的分子，并使用 $\lambda_e = \exp(\gamma_e)$，我们可以将其改写为
$$ \sum_{T \ni e} \exp(\gamma(T)) = \sum_{T \ni e} \prod_{e \in T} \lambda_e $$
既然我们现在知道缺少 $\lambda_e$ 项，我们可以将其添加到表达式中。
$$ \sum_{T \ni e} \lambda_e \times \prod_{f \in T, f \not= e} \lambda_f $$
使用求和规则，我们可以将 $\lambda_e$ 因子从求和中提出来得到
$$ \lambda_e \times \sum_{T \ni e} \prod_{f \in T, f \not= e} \lambda_f $$
由于我们使用的是对 $G \backslash \{e\}$ 应用基尔霍夫定理将产生除 $\lambda_e$ 因子之外的所有内容，我们可以手动将其乘回来。这将使 q 的伪代码变为
def q
input: e, the edge of interest
# Create the laplacian matrices
write lambda = exp(gamma) into the edges of G
G_laplace = laplacian(G, lambda)
G_e = nx.contracted_edge(G, e)
G_e_laplace = laplacian(G, lambda)
# Delete a row and column from each matrix to made a cofactor matrix
G_laplace.delete((0, 0))
G_e_laplace.delete((0, 0))
# Calculate the determinant of the cofactor matrices
det_G_laplace = G_laplace.det
det_G_e_laplace = G_e_laplace.det
# return q_e
return lambda_e * det_G_e_laplace / det_G_laplace
对 q 进行这个小小的更改非常有效。我能够将 $\delta$ 的减法改回 Asadpour 论文的做法，甚至在代码中添加了检查，以便每次我们更新 $\gamma$ 中的值时，我们都知道 $\delta$ 产生了正确的影响。
# Check that delta had the desired effect
new_q_e = q(e)
desired_q_e = (1 + EPSILON / 2) * z_e
if round(new_q_e, 8) != round(desired_q_e, 8):
raise Exception测试通过，没有任何错误！
下一步#
从技术上讲，在开始从中采样之前，我不知道这个分布是否正确。我已经将我一直在使用的测试编写成一个正式的测试，但由于我的预言机是程序本身，所以它唯一的失败方式是我在不知情的情况下改变了函数的行为。
因此，我必须继续编写 sample_spanning_tree，并为这两个函数获取更好的测试。
至于 spanning_tree_distribution 的测试，我当然希望添加更多测试用例。但是，如果 Held Karp 松弛返回循环作为答案，那么将有 $n - 1$ 个路径生成树，并且创建此分布的概念从一开始就不成立，因为我们已经找到了 ATSP 的解决方案。我真的需要更多真正的分数 Held Karp 解来扩展这两个函数的测试。
参考文献#
A. Asadpour, M. X. Goemans, A. Mardry, S. O. Ghran 和 A. Saberi，《不对称旅行商问题的 O(log n / log log n) 近似算法》，运筹学，65 (2017)，第 1043-1061 页。