Aztec 钻石图的完美匹配与多米诺洗牌算法
2021/03/15 更新：刚得知 Youtube 上的博主 mathologer 制作了一期非常精彩的节目，介绍 Aztec 钻石图与多米诺洗牌算法，非常值得一看！
2021/01/01 更新：2021 年的第一天，有人在 Shadertoy 上放了一个精彩的动画，演示多米诺洗牌算法的步骤：
Aigner 和 Ziegler 所著的 Proofs from the book （中文译版《数学天书中的证明》） 是一本非常精彩的数学读物，其中包含了 40 余个著名的数学问题和它们的巧妙解答。这些问题并不深奥，但也绝非没有受过严格数学训练的人所能欣赏，其中往往包含了相当的洞察力和聪明才智，读起来让人神清气爽，大叹数学之妙。
然而读完这本书的人恐怕都会有意犹未尽的感觉：这就没啦？我还没看够呢！
有一个问题我想是非常适合放在这本书里的，我也非常期待能在未来的版本中看到它，这就是 Aztec 钻石图的多米诺铺砌的计数问题。这个问题完美符合该书选题的标准：
- 表述初等，不需要太多的背景知识就能能理解。
- 内涵丰富。Aztec 钻石图是当前代数组合学中一个热点问题，它与交错符号矩阵、表示论、概率论、统计力学都有着深刻而奇妙的联系，有许多悬而未决的问题有待解决。
- 有多种令人拍案叫绝的解答，每个解答都不平凡，要么需要深刻的数学知识，要么需要开很大的脑洞。
我来介绍一下这个问题：
依次把 \(2,4,\ldots,2n\) 个单位正方形对称地摞在一起，放在 \(x\) 轴上方，然后关于 \(x\) 轴对称地反射过去，得到的图形叫做 \(n\) 阶的 Aztec 钻石图，记作 \(\mathrm{AZ}(n)\)。下图显示的是 \(\mathrm{AZ}(10)\)：
用 1x2 的多米诺骨牌不重叠不遗漏地盖住这些方格，可以得到 Aztec 钻石图的一个多米诺骨牌铺砌。下图显示的是 \(\mathrm{AZ}(10)\) 的一种可能的铺砌：
可以看到图中出现了四种不同颜色的骨牌，为什么这样染色后面会解释。
我们的问题是：
问题： \(\mathrm{AZ}(n)\) 总共有多少种不同的多米诺骨牌铺砌？
这里不考虑铺砌的对称性，比如全部用水平的骨牌铺砌和全部用竖直的骨牌铺砌是两种不同的铺砌。
问题： 如何在 \(\mathrm{AZ}(n)\) 的所有铺砌中等概率地随机任选一种？
第一个问题的答案是 \(2^{n(n+1)/2}\)，这个表达式很简洁，这暗示这个问题也应该有一个简洁的解答（确实如此）。
第二个问题的答案叫做多米诺洗牌算法，正是本文要介绍的。
Aztec 钻石图多米诺铺砌的计数问题最早由 Elkies、Kuperberg、Propp 和 Larson 在论文 (Elkies et al. 1991) 中进行了深入的研究，在这篇文章中他们一共给出了 4 种不同的解答。时至今日人们发现的解法已经超过一打，不幸的是没有一种可以算是「简单的」，但其中最精彩的仍然要数 Elkies 等人论文中的第四个解法：洗牌算法。后来 Propp 在另一篇文章 (Propp 2001) 中用图变换的方式重新表述了这个算法，并同时解决了在给定边的权重情况下对铺砌随机取样的问题。本文就来介绍 Propp 的方法。
洗牌算法
理解洗牌算法的第一个关键是理解骨牌的定向。
把 \(\mathbb{Z}^2\) 看做一张无穷大的国际象棋棋盘，每个方格染成黑白两色之一，相邻方格的颜色是不同的，于是每个骨牌都恰好覆盖一个黑方格和一个白方格。
乍看起来，骨牌只有水平或者竖直两种不同的类型，其实不然，骨牌有
N、S、W、E（北南西东）四种不同的类型，认识到这一点非常重要。下面介绍怎样定义和区分骨牌的类型。
如果棋盘上一个 \(2\times2\) 的正方形区域的左上角的方格是黑色的，我们就称这是一个黑方块，如下图所示：
类似地可以定义白方块为左上角方格为白色的 2x2 正方形区域。
对一张骨牌 \(d\)，它必然落在唯一的一个黑方块 \(B\) 中。规定 \(d\) 的移动规则为：把 \(d\) 移动到它在 \(B\) 中对面的位置上，根据 \(d\) 的移动方向规定其定向分别为
N、S、W、E。如下图所示：
注意到在移动以后，定向为 N/S 的骨牌变成了定向为
S 的骨牌，定向为 S 的骨牌变成了定向为
N 的骨牌，对 E, W
类型的骨牌亦然。
设 \(T\) 是 \(\mathrm{AZ}(n)\) 的一个多米诺铺砌，如果一个黑方块 \(B\) 恰好包含 \(T\) 中一对平行放置的骨牌，就称 \(B\) 是 \(T\) 的一个坏方块。在一个坏方块中，两个骨牌的定向是互相朝着对方的。
把 \(\mathrm{AZ}(n)\)
放在棋盘上，使得其最上方一行的最左边的方格是白色，规定 \(\mathrm{AZ}(n)\) 的对称中心为原点。设 \(T\) 是 \(\mathrm{AZ}(n)\)
的任一铺砌，如下图所示。图中根据骨牌的定向
N、S、W、E
将其染成了红、黄、绿、蓝四种颜色。我还画出了一个大小为 \(\mathrm{AZ}(n+1)\) 的背景区域。
移走所有坏方块中的骨牌，这一操作叫做「删除 」(deletion)。上图中的铺砌在删除后的结果如下图所示：
将剩下的骨牌按照其定向各自移动一步，这一操作叫做「移动」 (sliding)：
可以看到，这些剩下的骨牌分布在更大一些的 \(\mathrm{AZ}(n+1)\) 的区域内，而且不会出现骨牌重叠的情况。
可以证明在移动结束后 \(\mathrm{AZ}(n+1)\) 中的空白部分可以唯一地表示为若干不相交的黑方块的并，对每个这样的黑方块，任意用一对水平的或者是竖直的骨牌将其填充，就得到了 \(\mathrm{AZ}(n+1)\) 的一个铺砌。这一操作叫做「生成 」(creation)。上图中的空白区域可以表示为 14 个不相交的黑方块的并，因此一共有 \(2^{14}\) 种不同的生成方式，其中一种如下：
于是从 \(\mathrm{AZ}(n)\) 的任一铺砌出发，经过删除、移动、生成三次操作后，可以得到 \(\mathrm{AZ}(n+1)\) 的一个铺砌。这个步骤就叫做洗牌。
注意： 算法开始时要求 \(\mathrm{AZ}(n)\) 最左上边的方格是白色，但是经过一次洗牌操作后，得到的 \(\mathrm{AZ}(n+1)\) 的铺砌最左上边的方格是黑色。如果我们想从这个 \(\mathrm{AZ}(n+1)\) 的铺砌继续洗牌得到 \(\mathrm{AZ}(n+2)\) 的铺砌的话，就要翻转棋盘的染色。否则我们又会回到一个 \(\mathrm{AZ}(n)\) 的铺砌。
用 GIF 动图演示这个过程：
上图演示的是从 1 阶 Aztec 钻石图的一个随机铺砌出发，反复地执行删除 -> 移动 -> 生成 -> 删除 -> 移动 -> … 的步骤，最终生成 30 阶 Aztec 钻石图的随机铺砌。其中每次生成之后都翻转棋盘的染色。
算法中最难的部分是证明在第二步移动结束后，\(\mathrm{AZ}(n+1)\) 区域中的空白部分可以表示为不相交的黑方块的并。Elkies 等人的原证明很简单，但是给人的感觉不够优雅和本质。我们先承认这个结论是对的，然后可以给出计数问题的一个快速解答：假设一个 \(\mathrm{AZ}(n)\) 的铺砌 \(T\) 包含 \(k\) 个坏方块，则 \(\mathrm{AZ}(n)\) 有 \(2^k\) 个不同的铺砌，它们包含的坏方块恰好就是这 \(k\) 个。在删除和移动操作以后这 \(2^k\) 个铺砌给出相同的结果。\(\mathrm{AZ}(n+1)\) 的空白部分由 \(k+n+1\) 个黑方块组成，这里多出来的 \(n+1\) 是因为 \(\mathrm{AZ}(n+1)\) 比 \(\mathrm{AZ}(n)\) 多了 \(4(n+1)\) 个方格，所以要多 \(n+1\) 个方块。所以有 \(2^{k+n+1}\) 种不同的生成方式得到 \(\mathrm{AZ}(n+1)\) 的铺砌，数目正好是原来的 \(2^{n+1}\) 倍。这就是递推关系。
这个算法用 Propp 论文中介绍的图变换的角度更容易看清楚。
蜘蛛移动
多米诺铺砌实质上是图的完美匹配：考虑图 \(G_n\)，\(G_n\) 的顶点与 \(\mathrm{AZ}(n)\) 中的方格一一对应，两个顶点相邻当且仅当它们对应的方格在 \(\mathrm{AZ}(n)\) 中有公共的边，于是 \(\mathrm{AZ}(n)\) 的多米诺铺砌与 \(G_n\) 的完美匹配 (perfect matching) 一一对应。
现在问题转化为求 \(G_n\) 的完美匹配的个数，基本的想法是权函数。
设 \(G\) 是一个简单平面图，\(G\) 的每条边 \(e\) 有一个权值 \(w(e)\)，\(w(e)\) 是一个变量，比如 \(x,y,z,w\)。对 \(G\) 的一个完美匹配 \(\pi\)，定义 \(\pi\) 的权值 \(w(\pi)\) 为 \(\pi\) 中所有边的权值的乘积：
\[w(\pi)=\prod_{e\in \pi}w(e).\]
并定义 \(G\) 的权函数 \(w(G)\) 为 \(G\) 的所有完美匹配的权值的和： \[w(G)=\sum_{\pi}w(\pi).\]
如果 \(G\) 的每条边的权值都是 1，那么 \(w(G)\) 就是 \(G\) 的完美匹配的个数。但是一般情况下边的权值是未定元，所以 \(w(G)\) 是一个包含未定元的多元函数。但是只要求出了 \(w(G)\) 的表达式，把其中所有变元都赋值为 1，就得到了 \(G\) 的完美匹配的个数。
能求出权函数来当然是一件好事，因为权函数里面包含了图的非常多的信息，可以帮助我们计算出许多其它感兴趣的量来。比如说指定一条边 \(e\)，问 \(G\) 有多少个匹配不包含 \(e\)？为此只要把边 \(e\) 的权值设置为 0，其它的边权值保持为 1，代入权函数中即可。
看起来求出图的权函数是一个比直接计算其完美匹配的个数更复杂的问题，那为什么我们要舍近求远呢？权函数方法的奥秘在哪里呢？
这就是关键所在：对于一个复杂的图，我们想通过一些「手术」或者「变换」把它变成简单一些的图，比如删去一些顶点和边，或者把其中的一部分用一个新图去替换。这些变换通常会改变图的完美匹配的个数，但是权函数却在变换前后保持某种递推关系从而可以求解出来。这种操作我们其实都见到过，高中物理中电网络的各种等效电路替换（如 \(Y-\Delta\) 变换）就是图变换。
Propp 使用的图变换基于下面的蜘蛛引理：
蜘蛛移动. 假设图 \(G\) 的某个局部是一个正方形，如下图左边所示：
这里正方形的四个顶点记作 \(ABCD\)，四条边权值分别为 \(x,y,z,w\)。虚线表示它们可能与 \(G\) 的其它部分相连。
现在我们保持 \(G\) 其它的部分不动，将正方形 \(ABCD\) 替换为一个新的局部图，如上图右边所示。新局部的正方形部分四条边的权值分别是 \(z/\Delta,w/\Delta,x/\Delta,y/\Delta\)，其中 \(\Delta = xz+yw\)。多出来的四条边的权值都是 1。记替换后则得到的新图为 \(G'\)，则 \(G'\) 与原图 \(G\) 的权函数的关系为 \[w(G)=w(G')\cdot\Delta.\]
证明：对 \(G\) 的任一匹配 \(\pi\)，我们考察 \(\pi\) 在正方形 \(ABCD\) 处的匹配情况。根据 \(ABCD\) 中包含的 \(\pi\) 中边的个数，有三种可能：
- \(ABCD\) 正好匹配成两对，即它们包含 \(\pi\) 的两条边；
- \(ABCD\) 中有两个匹配在一起，另外两个与外部匹配，即它们包含 \(\pi\) 的一条边。
- \(ABCD\) 互不匹配，即它们不含 \(\pi\) 的边。
设 \(Q\) 是 \(\pi\) 在除去正方形 \(ABCD\) 后剩下部分的权。依次讨论这三种可能：
若 \(\pi\) 属于第一类情形，则有两种可能：\(\{A,B\}, \{C,D\}\) 或者 \(\{A,D\}, \{B,C\}\)。这两种可能对应于 \(G'\) 的一个匹配 \(\pi^\ast\)，如下图所示：
\(\pi\) 的两种可能性的权分别是 \(xzQ\) 和 \(ywQ\)，和是 \(\Delta Q\)；\(\pi^\ast\) 的权是 \(Q\)，变换后的权值是变换前的 \(1/\Delta\)。
若 \(\pi\) 属于第二类情形，则不妨设顶点 \(\{A,B\}\) 匹配（另外三种可能的情形是 \(\{B,C\}\)，\(\{C,D\}\)，\(\{A,D\}\)，分析是类似的），则 \(\pi\) 可以唯一地对应到 \(G'\) 的一个匹配 \(\pi^\ast\)，如下图所示：
\(\pi\) 的权是 \(xQ\)，\(\pi^\ast\) 的权是 \(xQ/\Delta\)，变换后的权值仍然是变换前的 \(1/\Delta\)。
若 \(\pi\) 属于第三类情形，则 \(\pi\) 可以映射为 \(G'\) 中的两个匹配 \(\pi^\ast_1,\pi^\ast_2\)：
\(\pi\) 的权是 \(Q\)，\(\pi^\ast_1,\pi^\ast_2\) 的权分别是 \(xzQ/\Delta^2\) 和 \(ywQ/\Delta^2\)，它们的和是 \(Q/\Delta\)，仍然是变换前的 \(1/\Delta\)。
于是 \(G\) 的匹配可以划分为不相交的三个子集 \(X_1,X_2,X_3\)，它们分别对应上面列出的三种情形。相应地 \(G'\) 的匹配也可以分为三个不相交的子集 \(Y_1,Y_2,Y_3\)。\(X_1\) 中每两个对应 \(Y_1\) 中的一个，\(X_2\) 和 \(Y_2\) 一一对应，\(X_3\) 中每一个对应 \(Y_3\) 中的两个。在这个对应下 \(Y_i\) 中所有匹配的权之和 \(w(Y_i)=w(X_i)/\Delta\)，从而 \(w(G')=\sum_{i=1}^3w(Y_i)=w(G)/\Delta\)。
下面这个引理在下节会用到，它的证明非常简单，所以我省略它的证明。
顶点收缩. 如果某个顶点只有两个邻居，并且它和这两个邻居之间的边的权值都是 1，则我们可以将这个顶点以及它的两个邻居收缩为一个顶点。这样得到的新图与原图有同样的权函数。如下图所示：
用图变换求解计数问题
我们用蜘蛛移动的技巧来计算 \(G_n\) 的权函数 \(w(G_n)\) 满足的递推关系。
把 \(\mathbb{Z}^2\) 的顶点间隔染成黑白两色。如果一个 \(2\times 2\) 的方块的左上角的顶点是黑色，就称它是一个「胞腔」。胞腔就是前面洗牌算法中的「黑方块」。对每个胞腔的四条边，从右侧开始顺时针依次标记其权重为 \(x,y,z,w\)。
我们从 \(G_{n-1}\) 开始。下图是 \(G_2\,(n=3)\) 的例子：
然后我们给 \(G_{n-1}\) 外侧装饰上一些边，把它变成一个更大的图 \(G_n^{(1)}\)。\(G_n^{(1)}\) 包含 \(G_n\) 作为子图，但是又比 \(G_n\) 多出一些蓝色的边。这些蓝色的边权值都是 1。\(G_n^{(1)}\) 的匹配在这些装饰的边处是限制死的，必须包含蓝色的边，所以 \(G_n^{(1)}\) 和 \(G_{n-1}\) 的权函数完全相同：
对 \(G_n^{(1)}\) 在各个胞腔处分别使用 蜘蛛移动，得到图 \(G_n^{(2)}\)：
由于我们总共对 \(n^2\) 个胞腔使用了蜘蛛移动，因此 \(G_n^{(2)}\) 的权函数是 \[w(G_n^{(2)})=\frac{w(G_n^{(1)})}{\Delta^{n^2}} = \frac{w(G_{n-1})}{\Delta^{n^2}}.\]
对 \(G_n^{(2)}\) 使用 顶点收缩，我们得到了 \(G_n^{(3)}\)：
这一步不改变权函数，所以 \[w(G_n^{(3)}) = w(G_n^{(2)}) = \frac{w(G_{n-1})}{\Delta^{n^2}}.\]
\(G_n^{(3)}\) 作为图和 \(G_n\) 是一样的，但是每个胞腔四条边的权值变成了 \(x/\Delta,y/\Delta,z/\Delta,w/\Delta\)，因此其权函数为 \(w(G_n)/\Delta^{n(n+1)}\)（因为 \(G_n\) 的每个匹配都包含 \(n(n+1)\) 条边，是 \(G_n\) 的顶点个数的一半）。从而我们用两种方法得到了 \(G_n^{(3)}\) 的权函数，即 \[\frac{w(G_{n-1})}{\Delta^{n^2}}=\frac{w(G_n)}{\Delta^{n(n+1)}}.\] 也就是 \(w(G_n)=\Delta^nw(G_{n-1})\)。这就是 \(\{w(G_n)\}\) 满足的递推关系。
由初始条件 \(w(G_1)=\Delta\) 可得 \(w(G_n)=\Delta^{n(n+1)/2}\)。特别令 \(w=x=y=z=1\) 我们就得到 \(\mathrm{AZ}(n)\) 的多米诺铺砌的个数为 \(2^{n(n+1)/2}\)。
从图变换的角度看洗牌算法
最后我们从图变换的角度解释为什么洗牌算法是正确的。你会看到，洗牌算法的三个操作步骤恰好对应 蜘蛛移动 中的三种情形。
在下图中，胞腔是深色，匹配的边是红色。三个胞腔 I, II, III 分别对应 蜘蛛移动 的三种情形：I 包含匹配的两条边；II 包含匹配的一条边；III 不含匹配的边。
蜘蛛移动以后的结果如下。注意初始 \(\mathbb{Z}^2\) 中的边不属于新图。
 |
 |
收缩后我们得到了 \(\mathbb{Z}^2\) 的一个新匹配。如下图所示：
与原匹配相比，I 中的边「删除」了，II 中的边「移动」到了对面，III 中「生成」了两条新边，正好与洗牌算法一致。同时胞腔的染色发生了翻转，即白方块和黑方块互换了颜色。
结语
Aztec 钻石图背后有许多有趣而深刻的数学，即便是计数这样看似初等的问题也有许多奥妙在里面，比如还有基于 graph condensation 的证明、转化为不相交路径组计数的证明、利用 \({\rm GL}_n(\mathbb{C})\) 的表示等精彩的证明，它们背后的思想已经广泛应用在计数组合学的许多问题中。借用丘吉尔的一句话作为结尾：“Now this is not the end. It is not even the beginning of the end. But it is, perhaps, the end of the beginning”。
顺便一提，演示多米诺算法的 GIF 动图是我学习 Python 的第一个正式程序，从开始看语法到写出最初的粗糙版本花了一周左右的时间。