Perfectly random sampling (1): domino shuffling algorithm on Aztec diamonds

Example:

Requirements:

  1. cairocffi and matplotlib for drawing random tilings.
  2. ImageMagick for making GIF animations.

Usage

To sample a random tiling of \(AZ(n)\), run

python random_tiling.py -size imgsize -order n -prog matplotlib

To make GIF animation of domino shuffling algorithm up to \(AZ(n)\), run

python domino_shuffling_animation.py -size imgsize -order n

Windows users need to manually set the variable CONVERTER to be the path to your convert.exe.

What is domino shuffling algorithm

The following picture shows an Aztec diamond of order 10:

You can see from the top row to the bottom row there are \[2, 4, \ldots, 20, 20, \ldots, 4, 2\] unit squares piled up layers by layers.

In general an Aztec diamond of order \(n\) (denoted as \(AZ(n)\) for short) is obtained via this way by piling up \[2, 4, \ldots, 2n, 2n, \ldots, 4, 2\] unit squares.

Now we consider domino tilings of \(AZ(n)\), i.e. tilings of \(AZ(n)\) with 1x2 dominoes. The following image shows one possible tiling of the above \(AZ(10)\):

Our problem is:

Problem 1: How many domino tilings of \(AZ(n)\) are there?

This problem seems quite innocent at the first glance: anyone with a high school level can understand it, but unfortunately all solutions known today (there are more than one dozen solutions) are not elementary, they either use quite sophisticated math or requre genius insights.

The answer is \(2^{n(n+1)/2}\), a short and neat expression which suggests the problem itself should also have a neat and beautiful solution, and it does. This solution is called the domino shuffling algorithm we implement here.

Before we state the algorithm, let's step further to another problem:

Problem 2: How to choose a random tiling of \(AZ(n)\) among all its \(2^{n(n+1)/2}\) tilings with uniform probability?

Compare this with the Wilson's uniform spanning tree algorithm and the coupling from the past algorithm in this repo, basically all of them are devised to do the same thing: sample with uniform probability from a very very large set (in fact they can be generalized to sample under any given probability distribution). For \(n=100\), \(AZ(100)\) has \(2^{5050}\) different tilings, far more larger than the number of partices in the universe!

Domino Shuffling Algorithm

Think the plane as an infinite chessboard, put \(AZ(n)\) at the origin, the chessboard are colored in the fashion that the leftmost cell in the top row of \(AZ(n)\) is white.

A 2x2 block is called "black" if its topleft cell is colored black.

Assume \(\mathcal{T}\) is a domino tiling of \(AZ(n)\), a 2x2 block is called "bad" (with respect to \(\mathcal{T}\)) if it's black and it contains a pair of parellel dominoes in \(\mathcal{T}\).

The algorithm runs as follows:

  1. Deletion: find all bad blocks and remove all pairs of parallel dominoes in them.

  2. Sliding: for each domino that remains in \(\mathcal{T}\), move it one step according to its orientation. The orientation of a domino \(d\) is determined by the following rule: find the unique 2x2 black block \(B\) that contains \(d\), move \(d\) to the oppsite position in \(B\). After sliding these dominoes are scattered in the region bounded by \(AZ(n+1)\).

  3. Creation: it can be proven that after the moves the holes in \(AZ(n+1)\) can be filled up by disjoint 2x2 black blocks in a unique way (this assertion is the crux of the algorithm), use either a pair of horizontal or vertical dominoes to tile each block, then one gets a domino tiling \(\mathcal{T}'\) of \(AZ(n+1)\). Further more, if \(\mathcal{T}\) is sampled from tilings of \(AZ(n)\) with uniform probability, then \(\mathcal{T}'\) is also sampled from tilings of \(AZ(n+1)\) with uniform probability.

So to sample a random tiling of \(AZ(n)\) with uniform probability, one just starts from a uniformly random tiling of \(AZ(1)\) (which is a 2x2 square, use a coin to choose one of its two tilings), and repeat the procedure deletion --> sliding --> creation --> deletion --> ... until one gets a tiling of \(AZ(n)\), then this tiling must be a perfectly random sampled one.

The proof of this algorithm is quite non-trivial and one may refer to the two papers by Elkies, etc. and Propp for more details.

About the code

I don't have much to say about the code, once you understand the algorithm the code becomes quite straight-forward. The only trap one should be aware of is, you must start the search from the boundary in the deletion and creation steps.

当前网速较慢或者你使用的浏览器不支持博客特定功能,请尝试刷新或换用Chrome、Firefox等现代浏览器