Google Summer of Code Logo with NetworkX logo

为了测试我使用spanning_tree_distribution生成的指数分布，我需要能够从该分布中抽样一个树。Asadpour论文中使用的主要参考文献是V. G. Kulkarni (1989) 的《生成随机组合对象》。虽然我无法找到这篇论文的在线版本，但密歇根理工大学图书馆有一个副本，我能够阅读它。

Kulkarni算法是否适用于Asadpour？#

Kulkarni在第2节中对该算法进行了概述，但第5节的标题为“随机生成树”，从第200页开始。首先，让我们检查一下Kulkarni论文第200页的预备知识是否与Asadpour算法匹配。

令$G = (V, E)$为一个具有$M$个节点和$N$条弧的无向网络……令$\mathfrak{B}$为$G$中所有生成树的集合。令$\alpha_i$为弧$i \in E$的正权重。定义生成树$B \in \mathfrak{B}$的权重$w(B)$为

$$w(B) = \prod_{i \in B} \alpha_i$$

同样定义

$$n(G) = \sum_{B \in \mathfrak{B}} w(B)$$

在本节中，我们描述了一种生成$B \in \mathfrak{B}$的算法，使得

$$P\{B \text{ is generated}\} = \frac{w(B)}{n(G)}$$

我们可以立即看到$\mathfrak{B}$与Asadpour论文中的$\mathcal{T}$相同，即所有生成树的集合。每条边的权重对于Kulkarni是$\alpha_i$，对于Asadpour是$\lambda_e$。至于图的权重的乘积作为概率，Asadpour论文在第382页指出

给定$\lambda*e \geq 0$对于$e \in E$，$G$的$\lambda$*-随机树_ $T$是从$G$的所有生成树集合中选择的树，其概率与$\prod_{e \in T} \lambda_e$成正比。

所以这不是问题。最后，$n(G)$可以写成

$$\sum_{T \in \mathcal{T}} \prod_{e \in T} \lambda_e$$

这在Asadpour论文中多次出现。因此，Kulkarni和Asadpour论文之间的预备知识是一致的。

Kulkarni算法#

Kulkarni给出的通用算法的专门版本是第202页上的算法A8。

$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$是所需的生成树。

现在我们必须理解这个算法，以便我们能够为它创建伪代码。首先作为符号解释，“生成$Z \sim U[0, 1]$”意味着在区间$[0, 1]$上选择一个均匀随机变量，该变量独立于之前生成的所有随机变量（有关更多信息，请参见Kulkarni第188页）。内置的python模块random可以在这里使用。查看实值分布，我相信使用random.uniform(0, 1)优于random.random()，因为后者没有生成“1”的概率，而这在Kulkarni论文中讨论的区间中明确包含。

另一个符号上的奇特之处将是类似于$G(U, V)$的语句，在这种情况下，它并不指以$U$作为顶点集和$V$作为边集的图，因为$U$和$V$都是完整边集$E$的子集。

$G(U, V)$在Kulkarni论文第201页中定义为

令$G(U, V)$为$G$的子图，通过删除不在$V$中的弧，以及收缩在$U$中的弧（即，识别$U$中弧的端节点）并删除由此产生的所有自环和收缩来获得。

这种语言似乎有点……笨拙，特别是对于$U$中的边。在这种情况下，“收缩在$U$中的弧”将是收缩这些边而不产生自环。幸运的是，此功能是NetworkX的一部分，使用networkx.algorithms.minors.contracted_edge并将self_loops关键字参数设置为False。

至于$E - V$中的边，这可以通过使用networkx.MultiGraph.remove_edges_from轻松实现。

一旦我们生成了$G(U, V)$，我们就需要找到$n(G(U, V)$。这可以使用我们已经熟悉的东西来完成：基尔霍夫树矩阵定理。我们只需要创建拉普拉斯矩阵，然后找到第一个余子式的行列式即可。此代码可能会直接从spanning_tree_distribution函数中获取。实际上，这是一个创建更广泛的辅助函数krichhoffs的地方，它将接收一个图并返回其中加权生成树的数量，然后将其用作spanning_tree_distribution中的q和sample_spanning_tree的一部分。

从这里，我们将$Z$与$\alpha_i \left(a’ / a\right)$进行比较，以查看该边是添加到图中还是被丢弃。理解算法的过程为$U$和$V$的含义提供了背景。$U$是我们已决定包含在生成树中的边的集合，而$V$是尚未考虑用于$U$的边的集合（粗略地说）。

现在，Kulkarni给出的算法中仍然存在一些模糊性，主要是关于$i$的。在循环条件中，$i$是从1到$N$的整数，即图中弧的数量，但稍后将其添加到$U$中，因此它必须是一条边。参考Asadpour论文，它在第383页开始描述抽样$\lambda$ -随机树，说“思路是任意对$G$的边$e_1, \dots, e_m$进行排序，并逐个处理它们”。所以我认为边的解释是正确的，Kulkarni使用的整数表示法假设已经发生了将边映射到${1, 2, \dots, N}$的操作。

sample_spanning_tree伪代码#

是时候编写一些伪代码了！从函数签名开始

def sample_spanning_tree
    Input: A multigraph G whose edges contain a lambda value stored at lambda_key
    Output: A new graph which is a spanning tree of G

接下来是初始化

    U = set()
    V = set(G.edges)
    shuffled_edges = shuffle(G.edges)

现在U和V的定义直接来自算法A8，但shuffled_edges是新的。我的想法是这将是我们用于$i$的东西。我们对图的边进行洗牌，然后在循环中，我们迭代shuffled_edges中的边。接下来我们有循环。

    for edge e in shuffled_edges
        G_total_tree_weight = kirchhoffs(prepare_graph(G, U, V))
        G_i_total_tree_weight = kirchhoffs(prepare_graph(G, U.add(e), V))
        z = uniform(0, 1)
        if z <= e[lambda_key] * G_i_total_tree_weight / G_total_tree_weight
            U = U.add(e)
            if len(U) == G.number_of_edges - 1
                # Spanning tree complete, no need to continue to consider edges.
                spanning_tree = nx.Graph
                spanning_tree.add_edges_from(U)
                return spanning_tree
        else
            V = V.remove(e)

主循环体确实使用了两个不在标准NetworkX库中的其他函数，krichhoffs和prepare_graph。正如我之前提到的，krichhoffs将对图应用基尔霍夫定理。此函数的伪代码如下，并严格基于spanning_tree_distribution的q中现有代码，该代码将更新为使用此新的辅助函数。

def krichhoffs
    Input: A multigraph G and weight key, weight
    Output: The total weight of the graph's spanning trees

    G_laplacian = laplacian_matrix(G, weight=weight)
    G_laplacian = G_laplacian.delete(0, 0)
    G_laplacian = G_laplacian.delete(0, 1)

    return det(G_laplacian)

另一个辅助函数prepare_graph的过程也给出了。

def prepare_graph
    Input: A graph G, set of contracted edges U and edges which are not removed V
    Output: A subgraph of G in which all vertices in U are contracted and edges not in V are
			removed

    result = G.copy
    edges_to_remove = set(result.edges).difference(V)
    result.remove_edges_from(edges_to_remove)

    for edge e in U
        nx.contracted_edge(e)

    return result

我想对NetworkX API进行的另一个更改是。目前，networkx.algorithms.minors.contracted_edge被编程为始终返回图的副本。由于我需要一次收缩多条边，因此就地执行收缩更有意义。我想向contracted_edge添加一个可选的关键字参数copy，它默认为True，以便整体功能不会改变，但我将能够执行就地收缩。

后续步骤#

最明显的一个是实现我在伪代码步骤中列出的函数，但测试仍然是一个令人担忧的领域。我最好的办法是抽样1000棵树，并检查每棵树的概率是否等于其边上所有lambda的乘积。

这实际上让我想到了一种新的spanning_tree_distribution测试。如果我生成分布，然后使用SpanningTreeIterator迭代所有生成树，我可以对每棵树被抽样的总概率求和，如果该和不等于1（或非常接近1），那么我没有一个有效的生成树分布。

参考文献#

A. Asadpour、M. X. Goemans、A. Mardry、S. O. Ghran 和 A. Saberi，《非对称旅行商问题的o(log n / log log n)近似算法》，SODA ’10，工业与应用数学学会，2010 年，第 379-389 页，https://dl.acm.org/doi/abs/10.5555/1873601.1873633。

V. G. Kulkarni，《生成随机组合对象》，算法杂志，11 (1990)，第 185–207 页。