殷勤昨夜三更雨 又得浮生一日凉

二维随机游动(二):一个随机的完美迷宫分别有多少个死角、拐角、岔路和十字路口?

赵亮 / 2020-10-06


此为二维随机游动系列的第二篇文章,但本文介绍的理论也适用于一般的图。

问题: 在 \(n\times n\) 的正方形网格图 \(G_n\) 的所有生成树中等概率地随机任选一个,记这个随机生成树为 \(T\)\(T\) 叫做 \(G_n\) 的一个均匀生成树。对 \(G_n\) 中任一顶点 \(v\)\(v\)\(T\) 的叶节点的概率是多少?

这个问题也可以表述为,"对一个完全随机的 \(n\times n\) 的完美迷宫,求它包含的死角的比例"。这里一个迷宫称作是完美的,如果迷宫中的任何两个房间之间都有且仅有唯一的道路相连 (这正是生成树的等价描述)。迷宫中的一个房间是称作是"死角",当且仅当它只有一条道路与其它房间相通。换句话说,一旦进去了就必须原路返回,没有其它出路 (这正是叶节点的等价描述)。

下图显示了三个不同的均匀生成树,它们分别来自大小为 \(80\times 80\)\(120\times120\)\(200\times200\) 的三个网格图,这三个生成树的叶节点 (用蓝色标出) 占全体顶点的比例分别为 \(1884/6400=0.294375\)\(4234/14400\approx0.294028\)\(11776/40000=0.2944\)。咦?看起来好像是在围着一个固定的值波动喔?

绘制上述插图的 Python 代码在这里,你可以自己试一试。

后面我们会证明对 \(G_n\) 的一个均匀生成树 \(T\)\(T\) 包含的叶节点的比例随着 \(n\) 趋于无穷收敛到常数 \(8(\pi-2)/\pi^3\approx 0.2945449\) (想不到 \(\pi\) 会出现在这里吧?),可见上面的三个例子与真实的极限值已经非常接近。

实际上对于度数是 \(2,3,4\) 的顶点,它们的比例也随着 \(n\to\infty\) 趋于固定的值,详细结果见下表:

度数 极限概率 近似值
1:死角 \[\dfrac{8}{\pi^2}(1-\dfrac{2}{\pi})\] 0.294
2:拐角 \[\dfrac{4}{\pi}(2-\dfrac{9}{\pi}+\dfrac{12}{\pi^2})\] 0.447
3:岔路 \[2(1-\dfrac{2}{\pi})(1-\dfrac{6}{\pi}+\dfrac{12}{\pi^2})\] 0.222
4:十字路口 \[(\dfrac{4}{\pi}-1)(1-\dfrac{2}{\pi})^2\] 0.036

所以在一个随机的完美迷宫中,最多的其实是"拐角" (45%),其次是"死角" (30%),再次是"岔路" (22%),而"十字路口"是非常少的 (3.6%)。

这些概率是怎样计算出来的呢?它和我们这个系列的主题二维随机游动又有何关系呢?这里面的故事相当精彩。粗略地说,图的随机生成树中的概率量可以翻译为该图上随机游动中某些事件发生的概率,利用随机游动与电网络的联系,我们又可以把这些事件发生的概率翻译为电网络的物理量,而这些物理量是可以用各种方法来计算的。在我们这个问题中,一个房间是死角/拐角/岔路/十字路口的概率都可以从一个矩阵出发通过求其子矩阵的行列式得到 (或者稍微加点变化),这个矩阵叫做转移电流矩阵 (transfer current matrix),是本文内容的核心。转移电流矩阵的物理意义很简单:将 \(G_n\) 看作一个电网络,每条边的电阻都是单位 1。在其中一条边 \(e\) 的两端加上电极,使得流入整个网络的电流总量为 1 个单位,记此时流过边 \(e'\) 的电流总量为 \(Y(e,e')\),则矩阵 \(Y=Y(e,e')_{e,e'\in E}\) 叫做网络 \(G_n\) 的转移电流矩阵。

转移电流矩阵有许多有趣的性质,比如:

所以解决问题的关键就是求解 \(\mathbb{Z}^2\) 上的转移电流矩阵,本文接下来就来介绍具体的理论。这个故事并无特别的难点,但确实比较长,把它勉强压缩在一篇博客文章中很难不影响可读性,所以我尽量只写与上面问题相关的部分,至于 Rayleigh 定律、等效电阻、常返与暂态的流判定等其它电网络中的重要内容就只好略去了。在写作过程中主要参考了 Lyons 等人关于图随机游动的专著 1,Stanley 的代数组合学教材 2 以及网上一些公开的讲义。

1 图上的线性代数

一个网络 \(\mathcal{N}=(G, c)\) 由一个有限图 \(G=(V,E)\) 和一个权值函数 \(c:E\to\mathbb{R}_{\geq0}\) 组成。这里 \(G\) 是连通的无向图,其每条边 \(e\in E\) 有一个权值 \(c(e)\in[0, +\infty)\),叫做 \(e\)电导,其倒数 \(r(e)=1/c(e)\) 叫做 \(e\)电阻

\(x,y\in V\) 是两个顶点,我们用 \(x\sim y\) 表示 \(x,y\) 是相邻的。由于本文主要考察边上的流函数,而流函数是有方向的,所以我们约定边集 \(E\) 同时包含有序对 \((x,y)\)\((y,x)\),即一条边作为有序对在 \(E\) 中出现两次。当 \(e=(x,y)\) 时,我们记 \(-e=(y,x)\)\(e\)反向边。并用记号 \(E_{1/2}\) 表示对每个 \(e\),在 \(e\)\(-e\) 中任选一个组成的集合。此外我们称 \(e^-=x\)\(e\)尾部\(e^+=y\)\(e\)头部。注意在反向边 \(-e\)\(x\) 变为头部,\(y\) 变为尾部。

对每个 \(e\in E\),当 \(e=(x,y)\) 时我们也记 \(c(x,y)=c(e)\)注意我们要求电导函数 \(c(x,y)\) 关于 \(x,y\) 是对称的:\(c(x,y)=c(y,x)\)。对每个顶点 \(x\in V\),记 \(c(x)=\sum_{y\sim x}c(x,y)\)

所有 \(V\) 上的实值函数 \(f:V\to\mathbb{R}\) 构成一个内积空间 \(\ell^2(V)\)\[(f,\, g)_V = \sum_{x\in V}c(x)f(x)g(x).\]

所有 \(E\) 上的反对称实值函数 \(\theta: E\to\mathbb{R}\) 也构成一个内积空间 \(\ell^2_{-}(E)\)\[(\theta,\,\theta')_E = \frac{1}{2}\sum_{e\in E}r(e)\theta(e)\theta'(e)=\sum_{e\in E_{1/2}}r(e)\theta(e)\theta'(e).\] 这里反对称的意思是对任何 \(e=\{x,y\}\)\(\theta(x, y)=-\theta(y,x)\)

以上定义基本上与 Lyons 等人的书保持一致,但是接下来的叙述会有比较大的差别。

定义 1.1

  1. 定义梯度算子 \(\nabla: \ell^2(V)\to\ell^{2}_{-}(E)\)\[ (\nabla f)(x,y) = c(x, y)\big[f(x)-f(y)\big]=c(e)\big[f(e^-)-f(e^+)\big]. \]

  2. 定义散度算子 \(\mathrm{div}:\ell^{2}_{-}(E)\to\ell^2(V)\)\[ (\mathrm{div}\,\theta)(x) = \frac{1}{c(x)}\sum_{y\sim x}\theta(x, y)=\frac{1}{c(x)}\sum_{e^-=x}\theta(e). \]

  3. 定义拉普拉斯算子 \(\Delta:\ell^2(V)\to\ell^{2}(V)\)\[ (\Delta f)(x) = (\mathrm{div}\,\nabla f)(x) = f(x)-\sum_{y\sim x}\frac{c(x,y)}{c(x)}f(y). \] 如果 \((\Delta f)(x)=0\)\(x\in V\) 成立,我们就称 \(f\)\(x\) 处是调和的。

定理 1.1:梯度算子和散度算子是一对共轭算子:对任何 \(f\in\ell^2(V)\)\(\theta\in\ell_{-}^2(E)\)\[(f,\,\mathrm{div}\,\theta)_V = (\nabla f,\,\theta)_E.\]

证明:直接验证即可,注意中间交换一次求和顺序: \[ \begin{align*} (\nabla f,\,\theta)_E&=\frac{1}{2}\sum_{e\in E}\theta(e)r(e)\nabla f(e)\\ &=\frac{1}{2}\sum_{e\in E}\theta(e)\bigg[f(e^-)-f(e^+)\bigg]\\ &\overset{*}{=}\sum_{e\in E}\theta(e)f(e^-)\\ &=\sum_{e\in E}\theta(e)\sum_{x\in V}1_{\{e^-=x\}}f(x)\\ &=\sum_{x\in V}f(x)\sum_{e\in E}\theta(e)1_{\{e^-=x\}}\\ &=\sum_{x\in V}f(x)\sum_{e^-=x}\theta(e)\\ &=\sum_{x\in V}f(x)c(x)(\mathrm{div}\,\theta)(x)\\ &=(f,\,\mathrm{div}\,\theta)_V. \end{align*} \] 其中 \((\ast)\) 式是因为当 \(e\) 跑遍 \(E\) 时,\(-e\) 也跑遍 \(E\),而 \(\theta(e)f(e^-)\)\(\theta(e)f(e^+)\)\(e\)\(-e\) 的取值是相反的。

定义 1.2:我们称映射 \(\mathrm{div}\) 的核 \(\mathrm{Ker\,div}\) 为环空间 (cycle space),记作 \(\lozenge\)。称映射 \(\nabla\) 的像 \(\mathrm{Im}\nabla\) 为星空间 (star space),记作 \(\bigstar\)。环空间和星空间都是 \(\ell_{-}^2(E)\) 的子空间。

定理 1.2\(\ell_{-}^2(E)\) 是环空间和星空间的正交直和:\(\ell_{-}^2(E)=\lozenge\oplus\bigstar\)

证明:由于 \(\{1_{x},x\in V\}\) 构成 \(\ell^2(V)\) 的一组基,所以它们在 \(\nabla\) 下的像构成 \(\bigstar\) 的一组基。对任何 \(\theta\in\ell^2_{-}(E)\)\[(\nabla 1_{x},\,\theta)_E=(1_x,\,\mathrm{div}\,\theta)_V=c(x)(\mathrm{div}\,\theta)(x).\] 于是 \(\theta\) 与任何 \(\nabla 1_x\) 正交当且仅当 \((\mathrm{div}\,\theta)(x)=0\) 对任何 \(x\in V\) 成立,即当且仅当 \(\theta\in\lozenge\),从而 \(\lozenge\)\(\bigstar\) 互为正交补空间,证毕。


我们解释一下环空间这个名字的来历。对任一 \(e\in E\),记 \[ \chi^e=1_{e} - 1_{-e}\in\ell^2_{-}(E) \] 为边 \(e\) 的示性函数。设 \(e_1=(v_1,v_2), e_2=(v_2,v_3),\ldots,e_k=(v_k,v_1)\) 首尾相接构成一条回路 \(C\),即 \[ v_1\xrightarrow{e_1}v_2\xrightarrow{e_2}\cdots\xrightarrow{e_{k-1}}v_k\xrightarrow{e_k}v_1. \] 定义 \(C\) 对应的环函数 \(f^C = \sum_{i=1}^k\chi^{e_i}\),则不难验证 \(f^C\in\lozenge\)。我们断言 \(\lozenge\) 可以由所有的 \(f^C\) 张成,这就是为什么 \(\lozenge\) 叫做环空间。

\(\theta\in\lozenge\),定义 \(\theta\) 的支集为 \[ \mathrm{supp}(\theta)=\{e\in E\,|\,\theta(e)\ne0\}. \] 注意如果 \(\mathrm{supp}(\theta)\) 非空的话,则它必然包含一个回路,否则 \(\mathrm{supp}(\theta)\) 构成的子图中必然有一个叶顶点 (度数恰好为 1),从而这一点处的散度不可能为 0,矛盾!于是 \(\mathrm{supp}(\theta)\) 包含一个回路 \(C\),我们可以选择合适的常数 \(\alpha\) 使得 \(\widetilde{\theta}=\theta-\alpha f^C\) 在某个 \(e\in C\) 上为 0,则 \(\widetilde{\theta}\) 的支集严格小于 \(\theta\)。对 \(\widetilde{\theta}\) 重复此操作一直下去必然可以在有限次后将支集变为空集,这就证明了 \(\theta\) 可以表示为形如 \(f^C\) 的线性组合。

2 网络中的流

仍然设 \(\mathcal{N}=(G,c)\) 是一个网络,\(Z\)\(V\) 的一个子集,\(a\notin Z\) 是一个顶点。

定义 2.1:我们称 \(\theta\in\ell_{-}^2(E)\) 是一个从 \(a\)\(Z\) 的流,如果以下两个条件成立:

  1. \((\mathrm{div}\,\theta)(a)\geq 0\)
  2. 对任何 \(x\notin \{a\}\cup Z\)\((\mathrm{div}\,\theta)(x)=0\)

\(a\) 叫做 \(\theta\) 的 "源点",\(Z\) 叫做 \(\theta\) 的 "汇点"。

定义流 \(\theta\) 的强度为从源点 \(a\) 流入网络的量: \[\mathrm{strength}(\theta)=\sum_{x\sim a}\theta(a, x).\] 我们也把 \(\mathrm{strength}(\theta)\) 记作 \(\|\theta\|\),如果有 \(\|\theta\|=1\) 成立,就称 \(\theta\) 是一个单位流

此外定义 \(\theta\) 的能量为 \[ \mathcal{E}(\theta)=(\theta,\,\theta)_E. \]

此定义自动蕴含了 \(\sum_{z\in Z}(\mathrm{div}\,\theta)(z)=-(\mathrm{div}\,\theta)(a)\leq=0\):考虑 \(V\) 上的常函数 \(1\),则 \(\nabla 1=0\),于是 \[0=(\nabla 1,\,\theta)_E=(1,\,\mathrm{div}\,\theta)_V=(\mathrm{div}\,\theta)(a)+\sum_{z\in Z}(\mathrm{div}\,\theta)(z).\] 即对整个网络来说,从 \(a\) 流出的量等于流入 \(Z\) 的量,中间没有任何"淤积"。

定义 2.2:设 \(\theta\) 是一个从 \(a\)\(Z\) 的流,如果还有 \(\theta=\nabla g\in\bigstar\) 成立,就称 \(\theta\) 是一个电流\(g\) 叫做 \(\theta\) 的电势。

我们总结一下流和电流的特点:

  1. 你可以把 \(a\) 理解为电源正极,\(Z\) 理解为网络接地的部分。
  2. 流函数要求对任何 \(x\notin \{a\}\cup Z\)\((\mathrm{div}\,\theta)(x)=0\) 等价于物理中的基尔霍夫定律:任何中间节点的流量守恒。
  3. 电流函数要求 \(\theta=\nabla g\in\bigstar\) 等价于物理中的欧姆定律,即电流 \(\theta\) 等于电势差 \(\nabla g\) 除以电阻 \(r(e)\)

在本文中,我们最关心的电流函数是当将一条给定边 \(e\) 的一个端点 \(x\) 通入单位电流,另一端点 \(y\) 接地时网络中对应的电流函数,此函数记作 \(i^e\)。于是对任一边 \(e'\)\(i^e(e')\) 就是此时流过 \(e'\) 的电流量。

推论 2.1\(\chi^e\)\(\bigstar\) 上的正交投影即为 \(i^e\)

证明:显然 \(\sum_{y\sim e^-}\chi^e(e^-, y)=1\),以及 \(\mathrm{div}\,\chi^e(x)=0\) 对任何 \(x\notin e\) 成立。而到 \(\bigstar\) 的投影不改变任何顶点处的散度值,所以此投影分量满足基尔霍夫定律和欧姆定律,且是从 \(e^-\)\(e^+\) 的单位流,从而必然是电流 \(i^e\)

推论 2.2:设 \(e,e'\) 是两条边,则 \[i^e(e')=c(e')(i^e, i^{e'})_E.\]

证明:记 \(P_\bigstar\) 为从 \(\ell^2_{-}(E)\)\(\bigstar\) 的投影算子,则 \(P_\bigstar\) 是自共轭算子以及 \(P_\bigstar^2=P_\bigstar\),从而 \[(i^e, i^{e'})_E=(P_\bigstar\chi^e,P_\bigstar\chi^{e'})_E=(P_\bigstar\chi^e,\chi^{e'})_E=(i^e,\chi^{e'})=r(e')i^e(e').\] 即得结论。

定义 2.3:记 \(Y(e,e')=i^e(e')\),则矩阵 \(\big(Y(e,e')\big)_{e,e'\in E}\) 叫做转移电流矩阵。\(Y(e,e')\) 满足互反律 \[Y(e,e')r(e')=Y(e',e)r(e)=(i^e,i^{e'})_E.\]

3 电网络物理量的概率解释

在这一节中,我们会给出电势、电流对应的概率解释。假设 \(a\) 是一个顶点,\(Z\) 是一个顶点集合且 \(a\notin Z\)

电压的概率解释 (第一种)

电压有一个很简单的解释,它是直线上随机游动中赌徒获胜/破产概率的类比。

考虑从 \(x\) 出发的随机游动,设 \(\tau_a\) 是这一随机游动首次到达 \(a\) 点的时间 (类比赌徒获胜),\(\tau_Z\) 是这一随机游动首次到达集合 \(Z\) 的时间 (类比赌徒输光),并考虑概率 \[h(x) = \mathbb{P}_x(\tau_a<\tau_Z).\] 其中我们规定 \(h(a)=1, h(Z)=0\),则 \(h(x)\) 在任何 \(x\notin \{a\}\cup Z\) 上是调和的: \[h(x) = \sum_{y\sim x}P(x, y)h(y).\] 此即为全概率公式,下一步到某个 \(y\sim x\) 的概率是 \(P(x,y)\),从 \(y\) 继续出发获胜的概率是 \(h(y)\),求和即可。

现在假设 \(a\) 点加上的电压大小是 1,则其电势 \(v(a)=1\),由于 \(Z\) 接地所以 \(v(Z)=0\)。我们断言电势函数 \(v(x)\) 在任何一点 \(x\notin\{a\}\cup Z\) 处也是调和的: \[v(x) = \sum_{y\sim x}P(x, y)v(y).\] 这是因为设 \(v\) 对应的电流是 \(i\),则 \[ \begin{align*} v(x) - \sum_{y\sim x}P(x, y)v(y)&= \sum_{y\sim x}P(x, y)(v(x)-v(y))\\&= \sum_{y\sim x}\frac{c(x,y)}{c(x)}i(x,y)r(x,y)\\&= \frac{1}{c(x)}\sum_{y\sim x}i(x,y)\\&= 0. \end{align*} \] 其中最后一个等号是根据基尔霍夫定律对任何中间节点 \(x\)\((\mathrm{div}\,i)(x)=0\)

现在获胜概率 \(\mathbb{P}_x(\tau_a<\tau_Z)\) 和电势函数 \(v(x)\)\(\{a\}\cup Z\) 上有相同的边界条件,在 \((\{a\}\cup Z)^c\) 上都是调和的,则它俩必然是同一个函数!(提示:考虑两个函数之差并利用调和函数必然在边界上取得最大最小值,注意这是在有限图上) 于是我们得到

电压的概率解释:在 \(a\) 加上单位 1 的电压,将 \(Z\) 接地,则任一顶点 \(x\) 处的电势等于从 \(x\) 出发的随机游动在到达 \(Z\) 之前先到达 \(a\) 的概率:\(v(x)=\mathbb{P}_x(\tau_a<\tau_Z)\)

如果 \(a\) 处不是单位电压的话,则结论中的 \(v(x)\) 应替换为 \(v(x)/v(a)\)

电流的概率解释

\(\mathbb{P}(a\to Z)=\mathbb{P}_a(\tau_Z < \tau_a^+)\) 为从 \(a\) 出发的随机游动在返回 \(a\) 之前先到达 \(Z\) 的概率。此概率称作从 \(a\)\(Z\)逃逸概率

仍然在 \(a\) 处加上电压,\(Z\) 接地,考虑此网络的电势 \(v\) 和电流 \(i\)

引理 3.1\[\frac{v(a)}{\| i\|} = \frac{1}{c(a)\mathbb{P}(a\to Z)}.\]

证明是常用的单步转移概率: \[ \begin{align*} \mathbb{P}(a\to Z) &= \sum_{x\sim a}P(a, x)\mathbb{P}_x(\tau_Z < \tau_a)\\&= \sum_{x\sim a}\frac{c(a,x)}{c(a)}\left(1-\frac{v(x)}{v(a)}\right)\\&= \sum_{x\sim a}\frac{c(a,x)}{c(a)}\frac{v(a)-v(x)}{v(a)}\\&= \sum_{x\sim a}\frac{i(a,x)}{c(a)v(a)}\\&= \frac{\|i\|}{c(a)v(a)}. \end{align*} \]

对任何边 \(x\sim y\),定义随机变量 \(N_{x\to y}^Z\) 为从 \(a\) 出发的随机游动在访问 \(Z\) 之前从 \(x\) 走到 \(y\) 的次数,则我们有如下定理:

定理 3.1:设电流是单位电流:\(\|i\|=1\),此时的电势函数为 \(v\),则

  1. \(v(x) = \mathcal{G}_Z(a, x)/ c(x)\)
  2. \(i(x,y)=\mathbb{E}_a\left[N_{x\to y}^Z - N_{y\to x}^Z\right]\)

其中 \(\mathcal{G}_Z(a, x)\) 是随机游动的 Green 函数,其定义为从 \(a\) 出发的随机游动在到达 \(Z\) 之前访问 \(x\) 的期望次数: \[\mathcal{G}_Z(a, x) = \sum_{k=0}^{\tau_Z}1_{\{X_k=x\}}.\]

证明

  1. 首先边界条件给出了 \(\mathcal{G}_Z(a, Z)=0\)。此外 \(\mathcal{G}_Z(a, a)\) 是一个每次试验成功概率为 \(\mathbb{P}(a\to Z)\) 的几何分布的随机变量的期望,故 \[\mathcal{G}_Z(a, a)=\frac{1}{\mathbb{P}(a\to Z)},\] 再由上面的引理,有 \[\frac{\mathcal{G}_Z(a, a)}{c(a)}=\frac{1}{c(a)\mathbb{P}(a\to Z)}=v(a).\] 所以 \(\mathcal{G}_Z(a, x)/c(x)\) 在边界上与电势函数是一致的,只要再证明在任意中间节点 \(x\) 的调和性质。注意随机游动必然是从与 \(x\) 相邻的顶点走过去访问 \(x\),所以 \[\frac{\mathcal{G}_Z(a, x)}{c(x)}=\frac{1}{c(x)}\sum_{y\sim x}\mathbb{E}_a[N_{y\to x}^Z].\] 另一方面 \[ \begin{align*} \mathbb{E}_a[N_{y\to x}^Z]&= \mathbb{E}_a\left[\sum_{t=0}^{\tau_Z}1_{\{X_t=y,X_{t+1}=x\}}\right]\\&= \sum_{t=0}^{\tau_Z}\mathbb{P}\left[X_t=y,X_{t+1}=x\right]\\&= \sum_{t=0}^{\tau_Z}\mathbb{P}(X_t=y)P(y,x)\\&= P(y,x)\mathbb{E}_a\left[\sum_{t=0}^{\tau_Z}1_{\{X_t=y\}}\right]\\&= P(y,x)\mathcal{G}_Z(a,y). \end{align*} \label{eq:cross} \] 所以代入求和后有 \[ \frac{1}{c(x)}\sum_{y\sim x}\mathbb{E}_a[N_{y\to x}^Z]=\frac{1}{c(x)}\sum_{y\sim x}P(y,x)\mathcal{G}_Z(a,y)=\sum_{y\sim x}P(x,y)\frac{\mathcal{G}_Z(a,y)}{c(y)}. \] 其中我们使用了 Markov 链的可逆性:\(c(x)P(x,y)=c(y)P(y,x)\)。这就证明了 \(\mathcal{G}_Z(a, x)/c(x)\) 在中间节点的调和性,从而它就是电势函数,这证明了 1。
  2. \[ \begin{align*} \mathbb{E}_a\left[N_{x\to y}^Z - N_{y\to x}^Z\right]&=P(x,y)\mathcal{G}_Z(a,x)-P(y,x)\mathcal{G}_Z(a,y)\\&= P(x,y)c(x)v(x) - P(y,x)c(y)v(y)\\&= c(x,y)(v(x)-v(y))\\&= i(x,y). \end{align*} \]

4 Wilson 均匀生成树算法

这一节的目的是介绍怎样将一条给定边 \(e\) 属于一个随机生成树 \(T\) 的概率 \(\mathbb{P}(e\in T)\) 表示为随机游动中的概率。这里的关键是 Wilson 均匀生成树算法,它通过擦圈的随机游动按照一定的权重分布在全体生成树中随机地取样,这里一个生成树 \(T\) 的权重 \(w(T)=\prod_{e\in T}c(e)\),任一生成树 \(T_0\) 被选中的概率为 \(w(T_0)/\sum_{T}w(T)\)

在下文中我们"随机生成树"均指服从上述权重分布的随机生成树。

Wilson 算法的步骤如下:

  1. 任取一个根节点 \(r\),维护一个树 \(T\),初始时 \(T=\{r\}\)
  2. 任取一个顶点 \(v\notin T\),从 \(v\) 出发作 \(G\) 上的擦圈随机游动,直到随机游动撞到 \(T\) 为止,然后将随机游动的路径 \(p\) 加入 \(T\) 中:\(T=T\cup p\)
  3. 重复步骤 2 直到 \(T\) 是一个生成树为止。

下面这个动画演示了 Wilson 算法的过程,你可以随时单击画布来重启动画。

Wilson 算法正确性的证明比较繁琐,这里不再赘述,你可以参考我之前的一篇文章。初次遇到这个算法的话可以先跳过证明,只记住结论。这个算法的重点就是初始时的根节点是可以任意选取的,每次进行擦圈的随机游动时选择的出发顶点也可以是任意的,这都不影响得到的生成树 \(T\) 的分布。

定理 4.1:设 \(\mathcal{N}=(G,c)\) 是一个网络,给定边 \(e=(x,y)\in E\),在 \(x\) 处加上电源并将 \(y\) 接地,使得流入网络的电流是单位流 \(i^e\)。则对 \(G\) 的随机生成树 \(T\)\(e\) 属于 \(T\) 的概率即为此时 \(e\) 中的电流量 \(i^e(e)\),它也等于从 \(x\) 出发的随机游动首次访问 \(y\) 是从边 \(e\) 走过去的这一事件的概率: \[\mathbb{P}(e \in T) = i^e(e) = \mathbb{P}_x\left[\text{the first visit to $y$ is across $(x,y)$}\right].\]

一个非常出人意料的结论,它把生成树、随机游动和电网络统一在一起。

证明:在 Wilson 算法中将初始根节点选为 \(y\),并将第一次出发作擦圈随机游动的顶点选为 \(x\),则此随机游动首次到达 \(y\) 有两种可能:

  1. 是从边 \((x,y)\) 走到 \(y\) 的,此随机游动的路径属于 \(T\),自然也有 \(e\in T\)
  2. 是从其它边走到 \(y\) 的,则此随机游动的路径属于 \(T\) 但不可能有 \(e\in T\),否则 \(T\) 中会出现回路。

于是我们证明了 \(e\in T\) 当且仅当从 \(x\) 出发的随机游动首次到达 \(y\) 是从边 \(e\) 走过去的。

另一方面在定理 3.1 中取 \(Z=\{y\}\)\[ i^e(x,y)=\mathbb{E}_x\left[N^{\{y\}}_{x\to y} - N^{\{y\}}_{y\to x}\right]. \] 注意到从 \(x\) 出发的随机游动不可能在到达 \(y\) 之前从 \(y\) 走到 \(x\),所以 \(N^{\{y\}}_{y\to x}=0\)。同时它在到达 \(y\) 之前从 \(x\) 走到 \(y\) 的次数只能是 0 或者 1,所以 \(N^{\{y\}}_{x\to y}\in\{0,1\}\),从而 \[\mathbb{E}_xN^{\{y\}}_{x\to y}=\mathbb{P}_x(N^{\{y\}}_{x\to y}=1)=\mathbb{P}_x\left[\text{the first visit to $y$ is across $(x,y)$}\right]. \] 这就完成了证明。

5 图的收缩

\(\mathcal{N}=(G,c)\) 是一个网络,\(F\subset E\) 是边的子集且 \(F\) 不含回路 (不考虑边的方向),\(G\) 关于 \(F\) 的图收缩是一个图变换,它把每个 \(F\) 中的每条边 \(f\) 的两个端点合并为一个新顶点,并将 \(f\) 变成新顶点的一条自边 (loop)。这样得到的新图记作 \(G/F\)。注意当一条边 \(e\notin F\) 如果和 \(F\) 中的边构成回路的话,它在 \(G/F\) 中也会变成自边。

在收缩变换下,\(G\) 的边集 \(E\)\(G/F\) 的边集 \(\widehat{E}\) 一一对应。不难验证如下结论:

引理 5.1:当 \(F\) 不含回路时,\(G\) 的所有包含 \(F\) 的生成树 \(\{T_G, F\subset T_G\}\)\(G/F\) 的所有生成树 \(\{T_{G/F}\}\) 一一对应,且 \[ w(T_G)=\prod_{f\in F}c(f)\cdot w(T_{G/F}). \] 于是对任何 \(e\in E\)\[ \mathbb{P}(e\in T_G \,\big| \, F\subset T_G) = \mathbb{P}(e\in T_{G/F})=\widehat{i^e}(e).\] 其中 \(\widehat{i^e}\)\(G/F\) 上的单位流 \(\chi^e\) 对应的电流函数。

我们想知道图 \(G\) 中的电流函数 \(i^e\) 和图 \(G/F\) 中的电流函数 \(\widehat{i^e}\) 之间的关系。这个关系乍看起来不知道怎么下手才好,但其实它有非常直观的解释:\(F\) 中的边在 \(G/F\) 中都变成了自边,这些自边中的电流都是 0,这时 \(G/F\) 中的电流函数 \(\widehat{i^e}\) 相当于在 \(G\) 中我们强制令 \(F\) 中的边的端点为等势点,并求解这时 \(G\) 上的电流函数 \(I^e\)。从物理上讲,这相当于我们在 \(e\) 的一端通入单位电流并将 \(e\) 的另一端接地的同时,还在每个 \(f\in F\) 的两端加上电源,使得 \(f\) 的两端电势相等,从而 \(f\) 中电流为 0。你可能要担心了:这不会在网络中引入新的电流吗?流过网络的还是单位流吗?答案是不会的,我们还要求这时的电流在任何 \(z\notin e\) 处满足基尔霍夫定律,这就不会有额外的电流进入网络了。

于是 \(I^e\) 必须满足以下三个条件:

  1. \(I^e\in \bigstar\)。(Ohm 定律)
  2. \(I^e\downharpoonright F=0\)。(任何 \(F\) 中的边两端等势)
  3. 流量限制:
    • \(F\) 的任何连通成分 \(F_C\),若 \(F_C\) 不包含 \(e\) 的任何端点,则从 \(F_C\) 净流出的电流量为 0。
    • \(F\) 的任何连通成分 \(F_C\),若 \(e^-\in F_C (e^+\in F_C)\) 则从 \(F_C\) 净流出的电流为 1 (-1)。
    • \(I^e\)\(G\) 的其它顶点 \(x\notin e\cup F\) 处流量守恒。

注意因为 \(F\) 的任何连通成分 \(F_C\) 会在 \(G/F\) 中 "坍缩" 为同一个顶点 \(x_C\),所以 \(x_C\) 处净的流出量对应于在 \(G\) 中将 \(F_C\) 作为一个整体净的流出量。由于 \(I^e\) 在每个 \(f\in F\) 上都是 0,所以 \(F_C\) 内部没有 "淤积" 的电流。

这个函数 \(I^e\) 不难找到:我们注意一个简单的线性代数事实:

Fact:若 \(U,V\) 是内积空间的子空间且 \(U\supseteq V^\bot\),则 \[U = V^\bot\oplus (U\cap V) = V^\bot\oplus P_VU.\] 其中 \(P_V\) 是从全空间到 \(V\) 的正交投影算子。

\(\chi^F=\{\chi^f,f\in F\}\) 为所有 \(f\in F\) 对应的单位流 \(\chi^F\) 张成的子空间,\(Z=\{i^f,f\in F\}=P_\bigstar\chi^F\),则 \(Z\subset\bigstar\)。对 \(U=\bigstar\)\(V=Z^\bot\) 使用上述事实得到 \[ \bigstar = Z \oplus P_{Z}^\bot\bigstar. \]推论 2.2 可知 \(I^e\downharpoonright F=0\) 等价于 \((I^e,i^f)_E=0\) 对任何 \(f\in F\) 成立,即 \(I^e\in Z^\bot\),而 \(I^e\in\bigstar\),从而 \(I^e\in \bigstar\cap Z^{\bot}=P_Z^\bot\bigstar\),因此这个 \(I^e\) 是什么就很清楚了,它正是电流函数 \(i^e\)\(Z^\bot\) 上的正交投影分量: \[I^e = P_Z^\bot i^e.\] 稍微麻烦一点的是验证 \(I^e\) 满足条件 3。记 \[I^e = i^e - \sum_{f\in F}\alpha_f f.\] 则在情形 1 中,每个 \(i^f\)\(f\) 两端的净流出量抵消,\(F_C\) 总的净流出量为 0;在情形 2 中边 \(e\) 只有一端在 \(F\) 中,其贡献的净流出量为 \(+1\) 或者 \(-1\),其它边 \(i^f\)\(f\) 两端的净流出量抵消,\(F_C\) 总的净流出量为 \(\pm1\);在情形 3 中 \(i^e\) 和所有 \(i^f\) 在顶点处的净流出量都是 0,\(F_C\) 总的净流出量为 0。

于是结合推论 2.2 我们有如下定理:

定理 5.1:设 \(F\subset E\) 不含回路,\(e\notin F\) 且与 \(F\) 中的边不构成回路,\(T_G\)\(G\) 的随机生成树,则 \(G/F\) 中的单位电流 \(\widehat{i^e}\)\(G\) 中单位电流 \(i^e\) 在子空间 \(Z^\bot\) 上的正交投影: \(\widehat{i^e}=P_{Z}^\bot i^e\)。于是 \[\mathbb{P}(e\in T_G \,\big| \, F\subset T_G) = \widehat{i^e}(e)=c(e)(\widehat{i^e},\widehat{i^e})_{\widehat{E}}=c(e)\big\|P_Z^\bot i^e\big\|^2.\] 其中 \(Z\) 是由 \(\{i^f,f\in F\}\) 张成的子空间,\(P_Z^\bot\) 是到 \(Z^\bot\) 的正交投影。

6 转移电流定理

定理 6.1 [转移电流定理]:设 \(\{e_1,\ldots,e_k\}\) 一组边的集合,\(T\) 是一个随机生成树,\(Y\)如前定义的转移电流矩阵,定义 \(k\times k\) 矩阵 \(Y^{(k)}(i,j)=Y(e_i,e_j)\),则 \[\mathbb{P}(e_1,\ldots,e_k\in T) = \det Y^{(k)}.\]

推论:考虑其中 \(r\) 条边 \(\{e_1,\ldots,e_r\}\) 不属于 \(T\),另外 \(k-r\) 条边 \(\{e_{r+1},\ldots,e_k\}\) 属于 \(T\) 的概率 \(\mathbb{P}(e_1,\ldots,e_r\notin T, e_{r+1},\ldots,e_k\in T)\),令矩阵 \(M^{(r)}\)\[ M^{(r)}(i, j) = \begin{cases} \delta_{ij} - Y^{(k)}(i, j), & i\leq r\\ Y^{(k)}(i,j), & r+1\leq i\leq k. \end{cases} \]\[ \mathbb{P}(e_1,\ldots,e_r\notin T, e_{r+1},\ldots,e_k\in T)=\det M^{(r)}. \]

我们首先说明怎样用转移电流定理得出推论的结论。对 "排除" 的边的个数 \(r\) 归纳,\(r=0\) 的情形对应的就是转移电流定理因而成立。设推论对任何小于 \(r\) 的情形成立,在 \(r\) 的情形考虑矩阵 \[ P(i, j) = \begin{cases} \delta_{ij} - Y^{(k)}(i,j), & i < r \\ \delta_{ij}, & i = r \\ Y^{(k)}(i,j), & r+1\leq i\leq k. \end{cases} \]\(M^{(r)},M^{(r-1)},P\) 在所有 \(i\ne r\) 的行都相等,但是 \(M^{(r)}\)\(M^{(r-1)}\) 的第 \(r\) 行之和等于 \(P\)\(r\) 行,于是根据行列式的多重线性性质, \[\det M^{(r)}+\det M^{(r-1)}=\det P.\]\(\det P\) 按照 \(P\) 的第 \(r\) 行展开并对 \(\{e_1,\ldots,e_{r-1},e_{r+1},\ldots,e_k\}\) 利用归纳假设,得到 \[ \det P=\mathbb{P}(e_1,\ldots,e_{r-1}\notin T, e_{r+1},\ldots,e_k\in T). \] 又对 \(\{e_1,\ldots,e_{r-1},e_r,\ldots,e_k\}\) 利用归纳假设得到 \[ \det M^{(r-1)} = \mathbb{P}(e_1,\ldots,e_{r-1}\notin T, e_{r},\ldots,e_k\in T) \] 相减即得结论。

回到转移电流定理的证明。分两种情形:

  1. 如果 \(\{e_1,\ldots,e_k\}\) 中包含某个回路 \(C\) 的话,我们可以适当重排它们并选取 \(a_i\in\{+1,-1,0\}\) 使得 \(\sum_i a_i\chi^{e_i}=f_C\in\lozenge\),这里对不属于 \(C\)\(e_i\)\(a_i=0\),对属于 \(C\)\(e_i\) 通过选择 \(a_i=\pm1\) 使得它们的方向为依次首尾相接。则根据推论 2.2 \(Y^{(k)}\) 的行是线性相关的,从而 \(\mathbb{P}(e_1,\ldots,e_k\in T)=\det Y^{(k)}=0\),结论成立。

  2. 如果 \(\{e_1,\ldots,e_k\}\) 中不包含回路的话,对 \(k\) 归纳,定理 4.1 已经给出了 \(k=1\) 的情形,所以可以假设 \(k>1\) 且结论对所有小于 \(k\) 的情形成立。 我们注意到 \(Y^{(k)}\) "差不多" 是一个 Gram 矩阵 (见之前的推论): \[ Y^{(k)}(i,j)=Y(e_i,e_j)=c(e_j)(i^{e_i},i^{e_j})_E. \]\(Y^{(k)}\) 不过是对每个 \(j\),在 Gram 矩阵 \(I^{(k)}=\big((i^{e_i},i^{e_j})\big)_{1\leq i,j\leq k}\) 的第 \(j\) 列上乘以 \(c(e_j)\) 得到的,所以 \[ \det Y^{(k)} = \left(\prod_{i=1}^k c(e_i)\right)\det I^{(k)}. \] 回忆内积空间中 Gram 矩阵的行列式的几何意义是其张成的平行多面体的体积的平方,它等于 \(\{i^{e_1},\ldots,i^{e_{k-1}}\}\) 张成的底面 \(Z\) 面积平方乘以 \(i^{e_k}\) 到这个底面距离的平方: \[ \det I^{(k)} = \big\|P_Z^\bot i^{e_k}\big\|^2\det I^{(k-1)}. \] 其中 \(Z=\mathrm{span}\{i^{e_1},\ldots,i^{e_{k-1}}\}\)\(P_Z^\bot\) 是到 \(Z\) 的正交补空间 \(Z^\bot\) 的投影算子。把 \(\prod_{i=1}^kc(e_i)\) 乘回去就得到 \[ \det Y^{(k)} = c(e_k) \big\|P_Z^\bot i^{e_k}\big\|^2 \det Y^{(k-1)}. \] 另一方面根据条件概率公式有 \[ \mathbb{P}(e_1,\ldots,e_k\in T)= \mathbb{P}(e_k\in T\,\big|\, e_1,\ldots,e_{k-1}\in T)\cdot \mathbb{P}(e_1,\ldots,e_{k-1}\in T). \] 而根据定理 5.1 的结论 \[ \mathbb{P}(e_k\in T\,\big|\, e_1,\ldots,e_{k-1}\in T) = c(e_k) \big\|P_Z^\bot i^{e_k}\big\|^2, \] 从而 \(\det Y^{(k)}\)\(\mathbb{P}(e_1,\ldots,e_k\in T)\) 有相同的初始条件和递推关系,从而它们相等,转移电流定理得证。

7 回答本文开头的问题

如本文开头所述,解决 \(\mathbb{Z}^2\) 中顶点在随机迷宫中度数分布的问题归结为求解 \(\mathbb{Z}^2\) 上的转移电流矩阵,即令所有边的电导值均为 1,在一条边 \(e=(x,y)\) 两端通入单位电流时,流过另一条边 \(e'\) 的电流流量,这等价于求解对应的电势函数 \(v\),也等价于求解 \(\mathbb{Z}^2\) 上拉普拉斯算子的解 \[ \Delta v (x) = 1_{x=e^-} - 1_{x=e^+}. \] 这一步可以通过 Fourier 分析的方法来做,比如 Lyons 等人的教材就是采用的这种方法。但是假定你看过本系列的第一篇文章的话,我们可以省很多事情,直接用二维随机游动的势核函数 \(a(x)\) 来构造电势函数 \(v\)。注意到 \(a(x)\) 在除去原点之外的任何点是调和的,并且在原点处的值为 0。记 \[ g(u) = \frac{a(u-y)-a(u-x)}{4}. \]\(g\) 在除去 \(x,y\) 外任何点处是调和的,从而 \(i=\nabla g\) 满足欧姆定律、在除 \(x,y\) 之外满足基尔霍夫定律,且 \(g(x)=a(x-y)/4 = -g(y)\)。又 \[ \begin{align*} \|i\| &= \sum_{z\sim x}i(x, z)=\sum_{z\sim x}[g(x)-g(z)]\\ &=a(x-y)-\sum_{z\sim x}\frac{a(z-y)-a(z-x)}{4}\\ &=\left(a(x-y)-\sum_{z\sim x}\frac{a(z-y)}{4}\right) + \sum_{z\in\{\pm(1,0),\pm(0,1)\}}\frac{a(z)}{4}\\ &=1. \end{align*} \] 其中最后一个等号是因为势核在除去原点之外是调和的,且在原点的四个邻居处值为 1。

于是 \(v(u)=g(u)-g(y)\) 即为当 \(x\) 通入单位电流,\(y\) 接地时 \(\mathbb{Z}^2\) 上的电势函数,从而其在任何一边 \(e'=(z,w)\) 上的电流为 \[ i(z,w) = g(z) - g(w)=\frac{1}{4}\big[a(z-y)-a(z-x)-a(w-y)+a(w-x)\big]. \]\(e=(x,y)\)\(e'=(z,w)\) 都取自原点相邻的四条边时,根据下图中的势核值代入上式即得本文开头的转移电流矩阵。

举个例子验证一下:取 \(e_1=(0, 0)\to(1,0)\)\(e_2=(0,0)\to (0,1)\),则 \(x=z=(0,0)\)\(y=(1,0)\)\(w=(0,1)\),于是 \[i^{e_1}(e_2)=\frac{a(-1,0)-a(0, 0)-a(-1,1)+a(0,1)}{4}=\frac{1}{2}-\frac{1}{\pi}.\] Bingo!


  1. Lyons, Russell, and Peres, Yuval. 2016. Probability on Trees and Networks. Cambridge Series in Statistical and Probabilistic Mathematics, vol. 42.

  2. Richard Stanley. 2013. Algebraic Combinatorics: Walks, Trees, Tableaux, and More. Undergraduate Texts in Mathematics. Springer-Verlag New York.