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)\) 的一种可能的铺砌：
可以看到图中出现了四种不同颜色的骨牌，为什么这样染色后面会解释。
我们的问题是：
问题 1：\(\mathrm{AZ}(n)\) 总共有多少种不同的多米诺骨牌铺砌？
这里不考虑铺砌的对称性，比如全部用水平的骨牌铺砌和全部用竖直的骨牌铺砌是两种不同的铺砌。
问题 2：如何在 \(\mathrm{AZ}(n)\) 的所有铺砌中等概率地随机任选一种？
第一个问题的答案是 \(2^{n(n+1)/2}\)，这个表达式很简洁，这暗示这个问题也应该有一个简洁的解答（确实如此）。
第二个问题的答案叫做多米诺洗牌算法，正是本文要介绍的。
Aztec 钻石图多米诺铺砌的计数问题最早由 Elkies、Kuperberg、Propp 和 Larson 在 1992 年的论文 Alternating sign matrices and domino tilings 中进行了深入的研究，在这篇文章中他们一共给出了 4 种不同的解答。时至今日人们发现的解法已经超过一打，不幸的是没有一种可以算是“简单的”，但其中最精彩的仍然要数 Elkies 等人论文中的第四个解法：多米诺洗牌算法。后来 Propp 在另一篇文章 Generalized domino-shuffling 中用图变换的方式重新表述了这个算法，并同时解决了在给定边的权重情况下对铺砌随机取样的问题。本文就来介绍 Propp 的结果。
本文的代码放在 Github 上。
什么是多米诺洗牌算法
理解多米诺洗牌算法的关键是理解多米诺骨牌的定向。
把 \(\mathbb{Z}^2\) 看做一张无穷大的国际象棋棋盘，每个方格染成黑白两色之一，相邻方格的颜色是不同的，于是任何 1x2 的骨牌都恰好覆盖一个黑方格和一个白方格。
乍看起来，骨牌只有水平或者竖直两种不同的类型，其实不然，骨牌有
N、S、W、E（北南西东）四种不同的类型，认识到这一点非常重要。下面介绍怎样定义和区分骨牌的类型。
如果棋盘上一个 \(2\times2\) 的正方形区域的左上角的方格是黑色的，我们就称这是一个黑方块，如下图所示：
类似地可以定义白方块为左上角方格为白色的 2x2 正方形区域。
对任意一张骨牌 \(d\)，它必然落在唯一的一个黑方块 \(B\) 中。规定 \(d\) 的移动规则为：把 \(d\) 移动到它在 \(B\) 中的对面的位置上，根据 \(d\)
的移动方向上、下、左、右规定其定向分别为
N、S、W、E。如下图所示：
注意到在移动以后，定向为 N 的骨牌变成了定向为
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)\) 的背景区域。
移走 \(T\) 的所有坏方块中的骨牌，这一操作叫做删除 (deletion)。上图中的铺砌在删除后的结果如下图所示：
将 \(T\) 中剩下的骨牌按照其定向各自移动一步，这一操作叫做移动 (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 等人的原证明很简单，但是并未从本质上解释多米诺算法为什么是成立的。这个算法用 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\) 的某个局部如下图左边所示：
这里位于中心处的胞腔的四条边的权值分别为 \(x,y,z,w\)，胞腔与周围四个顶点 \(A,B,C,D\) 相连的四条边的权值都是 1。
现在我们删去中心的胞腔，把 \(A,B,C,D\) 连接为一个新的胞腔，并且规定边的权值为（如图右边所示） \[x'=\frac{z}{\Delta},\quad y'=\frac{w}{\Delta},\quad z'=\frac{x}{\Delta},\quad w'=\frac{y}{\Delta}.\] 其中 \(\Delta = xz+yw\)，图 \(G\) 其它的部分保持保持不变，则得到的新图 \(G'\) 与原图 \(G\) 的权函数的关系为 \[w(G')=w(G)/\Delta.\]
证明：对 \(G\) 的任一匹配 \(\pi\)，我们考察 \(\pi\) 在这个局部的匹配情况。这有三种可能：
- \(A,B,C,D\) 无一与胞腔的四个顶点匹配；
- \(A,B,C,D\) 中有两个和胞腔中的顶点匹配，另外两个与外部的顶点匹配；
- \(A,B,C,D\) 与胞腔的四个顶点一一匹配。
设 \(Q\) 是 \(\pi\) 在除去这个局部外的权。依次讨论这三种可能：
- 若 \(\pi\) 属于第一类情形，则存在同属于第一类情形的另一个匹配 \(\pi'\)，\(\pi\) 和 \(\pi'\) 在外部是相同的，仅在胞腔处不同。两者在蜘蛛变换下对应于 \(G'\) 的同一个匹配（保持外部的匹配不变），如下图所示：
\(\pi\) 和 \(\pi'\) 的权一个是 \(xzQ\)，一个是 \(ywQ\)，二者的权和是 \(\Delta Q\)，变换后对应的匹配的权值是 \(Q\)，变换后的权值是变换前的 \(1/\Delta\)。
- 若 \(\pi\) 属于第二类情形，则不妨设顶点 \(\{B,C\}\) 与胞腔中的顶点匹配（另外三种可能的情形是 \(\{A,B\}\)，\(\{C,D\}\)，\(\{A,D\}\)，分析是类似的，而 \(\{A,C\}\) 和 \(\{B,D\}\) 是不可能出现的），则 \(\pi\) 可以唯一地映射为 \(G'\) 的一个匹配，如下图所示：
这里 \(\pi\) 在外部的部分仍然保持不变。\(\pi\) 的权值为 \(wQ\)，变换后的匹配的权值为 \(wQ/\Delta\)，变换后的权值仍然是变换前的 \(1/\Delta\)。
- 若 \(\pi\) 属于第三类情形，则 \(\pi\) 可以映射为 \(G'\) 中的两个匹配：
由于 \(\pi\) 在这个局部的四条边的权值都是 1，因此 \(\pi\) 的权值是 \(Q\)，变换后对应的两个匹配的权值分别是 \(y'w'Q\) 和 \(x'z'Q\)，它们的和 \((y'w'+x'z')Q=Q/\Delta\) 仍然是变换前的 \(1/\Delta\)。
于是 \(G\) 的匹配，根据属于情形 1-3，可以划分为不相交的三个子集 \(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\) 中的两个匹配。在这个对应下 \(w(Y_i)=w(X_i)/\Delta\)，从而 \(w(G')=w(G)/\Delta\)。
最后是一个简单的引理，在下一节要用到。
顶点分裂：设顶点 \(v\in G\) 的邻居 \(N(v)\) 可以分成两个不相交的集合 \(N(v)=H\cup K\)，这里 \(H\) 和 \(K\) 中允许有一个是空集。我们对图 \(G\) 作一个简单的变换：从 \(v\) 中分裂出两个新的顶点 \(v',v''\)，用 \(v'\) 代替 \(v\) 去和 \(H\) 中的顶点相连，用 \(v''\) 代替 \(v\) 去和 \(K\) 中的顶点相连，\(v\) 与 \(v'\) 和 \(v''\) 之间的边权值都定义为 1，则得到的新图 \(G'\) 与原图有同样的权函数：\(w(G)=w(G')\)。
证明是非常简单的，这里不再赘述。
铺砌计数的另一种解法
我们用蜘蛛移动技巧给出 Aztec 钻石图多米诺铺砌的另一种解答。将 \(G_n\) 倾斜 45 度放置，这样清楚地显示了 \(n^2\) 个胞腔，每个胞腔的边的权值与蜘蛛引理中一样，依顺时针方向为 \(x,y,z,w\)。设 \(G_n\) 的权函数为 \(w(G_n)\)，我们来求出序列 \(\{w(G_n)\}_{n\geq1}\) 所满足的递归关系：
首先对 \(G_n\) 的每一个顶点都使用顶点分裂，我们得到图 \(G_n'\)，新增加的边的权值都是 1：
由于顶点分裂不改变权函数，因此 \(G_n'\) 的权函数仍然是 \(w(G_n)\)。
然后对 \(G_n'\) 在各个胞腔处分别使用蜘蛛移动，得到图 \(G_n^{'\phantom{}'}\)：
\(n^2\) 个胞腔总共使用了 \(n^2\) 次蜘蛛移动，因此 \(G_n^{'\phantom{}'}\) 的权函数是 \(\frac{w(G_n)}{\Delta^{n^2}}\)。
我们再换另一个角度来计算 \(G_n^{'\phantom{}'}\) 的权函数，这也是整个证明峰回路转，柳暗花明之处！注意 \(G_n^{'\phantom{}'}\) 的匹配在边界上是限制死的，真正自由的部分是内部的一个低一阶的图 \(G_{n-1}\)，所以可以剥去外面的这层顶点！
而且这个 \(G_{n-1}\) 图的边的权值的标记规则是 \(z/\Delta,w/\Delta,x/\Delta,y/\Delta\)，因此其权函数为 \(\frac{w(G_{n-1})}{\Delta^{n(n-1)}}\)（因为 \(G_{n-1}\) 的每个匹配包含 \(n(n-1)\) 条边，是 \(G_{n-1}\) 的顶点个数的一半）。从而我们用两种方法得到了 \(G_n^{'\phantom{}'}\) 的权函数，即 \[\frac{w(G_n)}{\Delta^{n^2}}=\frac{w(G_{n-1})}{\Delta^{n(n-1)}},\] 也就是 \(w(G_n)=\Delta^nw(G_{n-1})\)。
注意 \(n=1\) 时 \(w(G_1)=\Delta\)，从而 \(w(G_n)=\Delta^{\frac{n(n+1)}{2}}\)。特别令 \(w=x=y=z=1\) 我们就得到 \(n\) 阶 Aztec 钻石图 \(\mathrm{AZ}(n)\) 的所有多米诺骨牌覆盖的个数为 \(2^{\frac{n(n+1)}{2}}\)。
多米诺洗牌算法和蜘蛛移动是等价的
最后我们来解释为什么多米诺洗牌算法和蜘蛛移动是等价的，特别还要补上前面多米诺洗牌算法表述中留下的坑：在 slide 操作后 \(\mathrm{AZ}(n+1)\) 区域内所有的洞可以唯一地表示为不相交的黑方块的并。下面将会看到虽然在多米诺洗牌算法中三个操作是有先后顺序的，但它们不过是蜘蛛移动作用在不同胞腔上的表现，即多米诺洗牌算法的删除、移动、填充三步分别对应于蜘蛛移动中的三种情形。
注意到对一个图 \(G\) 的所有胞腔都进行顶点分裂 + 蜘蛛移动后得到的是一个新图 \(G'\)，\(G'\) 的形状依赖于 \(G\) 的边界，因此不是每个胞腔都能恢复原状，可能某些胞腔会有悬在外面的边 （正如上面证明中看到的），但是当 \(G=\mathbb{Z}^2\) 是网格图时我们得到的仍然是 \(\mathbb{Z}^2\)。我们来考察 \(\mathbb{Z}^2\) 的一个局部的胞腔在变换下的结果。为了表述方便，我们使用先蜘蛛移动然后进行顶点合并的操作，这个操作是先顶点分裂再蜘蛛移动的逆操作，因此仍然把 \(\mathbb{Z}^2\) 变为自身。
delete：
其中第一个箭头表示蜘蛛移动，第二个箭头表示顶点合并后又回到一个正方形胞腔。slide：
同上，第一个箭头表示蜘蛛移动，第二个箭头表示顶点合并后又回到一个正方形胞腔。create：
同样地第一个箭头表示蜘蛛移动，第二个箭头表示顶点合并后又回到一个正方形胞腔。
你可以看到变换前后匹配的变化情况确实与洗牌算法的三步一一对应。注意到新匹配的边要么来自 slide 中移动的原匹配中的边，要么来自 creation 中补充的边，从而我们证明了如下结论：
推论：将 \(\mathbb{Z}^2\) 看做无穷大的国际象棋棋盘，设 \(\mathcal{T}\) 是其一个多米诺铺砌，则在 delete 和 slide 操作后平面上的空白区域可以表示为不相交的黑方块的并，因而可以用 create 操作补全为新的匹配。
那么怎么说明在 slide 操作后 \(\mathrm{AZ}(n+1)\) 区域内的洞可以唯一地表示为不相交的黑方块的并呢？
为此我们先将最初的 \(\mathrm{AZ}(n)\)
的铺砌扩展为整个平面的铺砌：对 \(\mathrm{AZ}(n)\) 外部的区域，位于 \(x\) 轴上方的方格一律用形如
白 | 黑 的水平骨牌铺砌，位于 \(x\) 轴下方的方格一律用形如
黑 | 白 的水平骨牌铺砌，容易看出在这个扩展后的铺砌中 \(\mathrm{AZ}(n)\)
外部没有坏方块，所有的坏方块都位于区域 \(\mathrm{AZ}(n)\) 中，所以 delete 操作后
\(\mathrm{AZ}(n)\) 外部区域不变；在
slide 操作中 \(\mathrm{AZ}(n)\)
的外部位于 \(x\)
轴上方的骨牌整体向上移动一步，位于 \(x\)
轴下方的骨牌整体向下移动一步，因此在得到的 \(\mathrm{AZ}(n+1)\)
的两侧出现了两个半无穷长的宽度为 2
的带状空白。整个棋盘的空白区域由这两个半无穷带状区域和 \(\mathrm{AZ}(n+1)\)
内部的洞组成。根据上面的推论，这些空白区域可以表示为不相交的黑方块的并，但是容易观察看出这些黑方块不可能跨越
\(\mathrm{AZ}(n+1)\)
的边界，从而两个半无穷带状区域可以表示为不相交的黑方块的并
（易见这个表示是唯一的），\(\mathrm{AZ}(n+1)\)
内部也可以表示为不相交黑方块的并
（易见也是唯一的，考虑最右上方的黑方块，这个必然是唯一确定的，拿走以后继续考虑剩下的黑方块中最右上方的，如此继续下去可得表示的唯一性）。
至此我们就从图变换的角度论证了多米诺洗牌算法的正确性。
结语
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 的第一个正式程序，从开始看语法到写出最初的粗糙版本花了一周左右的时间。