How to sample a uniform random polyomino?

Question

A polyomino is formed by joining finitely many unit squares edge to edge. It may be regarded as a finite subset of the regular square tiling with a connected interior. In particular, for us, polyominos may have holes.

(Update: these are sometimes also called two-dimensional lattice animals. For another model of random lattice animal, see this earlier MO question and discussion and references there.)

I would like to be able to sample a uniform random polyomino with $n$ squares, where $n$ is of the order $100$ or $1000$.

I am thinking of fixed polyominos, where two polyominos are considered equivalent only if they agree after a translation, but where rotations and reflections are in general distinct. But I might also be okay with sampling one-sided or fixed polyominos if that is for some reason easier or more natural.

The point is that I don't want to have to enumerate all the polyominos of this size. For $n=100$, enumerating all the polyominoes would already be far beyond the reach of modern computing. I am hoping that there is a random construction, or else a Markov chain that converges to the uniform distribution (i.e. every polyomino having equal probability). I would be interested in such a Markov chain, if it could be proved to converge to the uniform distribution, even without guarantees on the mixing time.

My first attempt at a random construction was to apply a randomized version of Redelmeier's method, or at least Wikipedia's description of this algorithm. But I am not sure if Wikipedia's description is correct. I copy it here, for the sake of completeness.

The basic idea is that we begin with a single square, and from there, recursively add squares. Depending on the details, it may count each n-omino n times, once from starting from each of its n squares, or may be arranged to count each once only.

The simplest implementation involves adding one square at a time. Beginning with an initial square, number the adjacent squares, clockwise from the top, 1, 2, 3, and 4. Now pick a number between 1 and 4, and add a square at that location. Number the unnumbered adjacent squares, starting with 5. Then, pick a number larger than the previously picked number, and add that square. Continue picking a number larger than the number of the current square, adding that square, and then numbering the new adjacent squares. When n squares have been created, an n-omino has been created.

This method ensures that each fixed polyomino is counted exactly n times, once for each starting square.

I was hoping that by making a uniform random legal choice at each step, I could generate a large random polyomino, but unfortunately the description of the algorithm above seems to sometimes get stuck. See, for example, the figure below. We have consructed a polyomino with 16 tiles, but no further moves are possible.

Answer 1

The following Markov chain Monte Carlo algorithm would suffice, as it satisfies detailed balance and is ergodic. It may suffer from some sort of "critical slowing down" as the proportion of moves which are legal might decay as a power of the system size, but I don't have clear intuition about this.

I don't know if there's anything better in the literature. This is a "local" move; maybe there is a "global" move which is akin to the pivot move but which still gives a correct algorithm (e.g. rotating/reflecting/translating a piece of the polyomino and re-attaching it).

I'll just refer to the objects as "polyominoes", but they could be referred to as lattice animals, or for dimensions greater than 2 as $d$-dimensional polycubes.

1. Initialise system, perhaps starting with a "straight rod" of length $N$, starting at the origin and extending in the positive x-axis direction, or (better) starting with a $d$-dimensional cube with a few cells missing, or (better still) starting with a configuration obtained by attaching cells at random until the correct number of cells have been added (following suggestion of Brendan McKay).
2. Choose a random occupied cell, and move it to a randomly chosen unoccupied neighbour of the remaining $N-1$ cells. (A list of the unoccupied neighbours needs to be maintained to do this efficiently. This data structure could be based on hash table but must allow insertions and deletions.)
3. If the resulting configuration is a polyomino, accept the move, otherwise reject the move and keep the original configuration. (The condition that needs to be checked is that the cells are connected. It suffices to check that the occupied neighbours of the original cell are connected, and this could be efficiently done via a greedy graph search algorithm.)
4. Return to step 2.

The algorithm generates translates of polyominoes, but that's trivial to partial out.

The Markov chain can generate any polyomino starting from a seed state, and is therefore ergodic. (Start with the rod seed state, and assume that the cell at the origin is the rightmost cell of an arbitrary polyomino that you want to construct. You can do this by taking the rightmost cell of the rod and building up the required configuration cell-by-cell to the left of the origin.)

The Markov chain satisfies detailed balance by inspection.

Therefore the Markov chain will sample polyominoes uniformly at random. It would take a bit of work to estimate the mixing time by running computer experiments and then coming up with a heuristic argument. I think it's unlikely that you will be able to get an exact result for the mixing time, but it may be possible to derive an upper bound.

One significant difficulty in implementation will be the efficient checking of connectivity for each attempted move. Naively, this could take CPU time $O(N)$, but in practice a greedy search algorithm should be much faster than that in the average case. It should be possible to do better by using a dynamic bounding volume hierarchy (like my SAW-tree implementation of the pivot algorithm for self-avoiding walks) or a dynamic data structure which maintains information about connectivity (can't remember the name of the appropriate data structure for fully dynamic insertions and deletions of vertices, but they are related to link/cut trees).

Better algorithms are surely possible - e.g. the idea sketched above for a global move - but I think that this will suffice to get random polyominoes with say 1000 cells, even for naive connectivity checking. If you want to do this for a million cells then we might have to think a bit harder!

Answer 2

Here is an approach to estimating parameters of uniformly random polyominoes of $n$ squares without generating random ones uniformly.

First define a rooted tree whose root is the single square and whose nodes at distance $m$ from the root are the $m$-square polyominoes, with no nodes further than $n$ from the root. This is fairly easy: the parent of polyomino $P$ can be $P-s$ where $s$ is the most-northern most-western square such that $P-s$ is a polyomino.

Now you need a procedure that, given a polyomino $P$ tells you how many children it has and chooses a random child. If nothing better comes to mind, just try adding a new square in each possible place and reject it if $P$ is not the parent.

Now start at the root and follow a random path down until you come to a leaf $L$ (with $m$ squares, say). At each node, choose a child uniformly at random. If $d_0,\ldots,d_{m-1}$ are the degrees (numbers of children) of the nodes along that path (excluding $L$), then $p(L)=(d_0\cdots d_{m-1})^{-1}$ is the probability that a random path ends at $L$.

Now, if $\chi(L)$ is the indicator function of "$L$ has $n$ squares", then the expectation of $\chi(L)/p(L)$ is exactly equal to the number of $n$-square polyominoes, so you can estimate that number by sampling repeatedly. In fact, for any function $f$, the sum of $f(P)$ over all $n$-vertex polyominoes is exactly the expectation of $f(P)\chi(L)/p(L)$ for a random path. This lets you estimate the expectation of $f$ by repeated sampling.

I've considered the simplest case, where all the children of a node are chosen with equal probability as a random path is traced. The theory is valid for any distribution on the children. To get good estimates you would like the probability of choosing a child to be roughly proportional to the number of $n$-square polyominoes descended from that child. Then each $n$-vertex polyomino comes out with similar probability. If that is far from the truth, the estimates may have large variance even though they have the correct expectation.

This is an example of "sequential importance sampling". Searching on that phrase will bring up some papers; look for Diaconis in particular. I've seen this work extremely well and extremely badly on different problems.

Answer 3

Consider reverse engineering the problem. Every sufficiently large polyomino has a decomposition into triominoes and tetraminoes (with perhaps a vanishingly small set of exceptions). Enumerate the polyominoes up to some feasible number (say 20), and see how one can decompose each animal into tetraminoes, and get a frequency of how often each tetramino is used. Then you can start with an appropriately small polyomino which agrees mod 4 with the target size, and randomly add tetraminoes, and see if this gives you the desired distribution. It may be that you can't just add tetraminoes, but need a mix of triominoes and tetraminoes.

It may also be that most animals can be bisected in two or fewer ways into animals with half the number of squares. If so, then most of the animals can be built up by merging two animals of almost the same size.

Gerhard "Merge And Conquer Your Menagerie" Paseman, 2016.07.17.