# How to sample a uniform random polyomino?

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.

• I wonder if it would help to look at the literature on lattice animals, as cited in this earlier MO question, "Random polycube shapes." Jul 17 '16 at 20:49
• I think that at this stage and the previous ones you are next allowed to pick any empty numbered square Jul 17 '16 at 22:41
• Aaron Meyerowitz, it seems that this keeps you from getting stuck, but then not every polyomino is equally likely. For example, the square tetromino is more likely than the straight tetromino. Jul 17 '16 at 23:12
• Igor Pak. Thanks for the reference. I have been looking at this book, but still haven't found the answer to my question. Maybe it is right in front of me, and I am missing it. As far as I understand it, self-avoiding random walks can be "nearly uniformly" sampled via relatively easy to implement Markov chain Monte Carlo methods such as the pivot algorithm, for example: arxiv.org/abs/cond-mat/0109308 . I guess I am asking if there are similar methods for simulating uniform polyominoes. I am guessing that proving things about mixing time is hard, and I am happy to just do experiments. Jul 18 '16 at 2:30
• I'll just add a comment on your description of Redelmeier's method: although it may get stuck this doesn't mean the algorithm is incorrect, instead you just have to repeatedly perform the procedure until you have successfully generated a polyomino of the correct size. Although this probability will decay exponentially with size N, it may be good enough in practice for sizes you are interested in. Self-avoiding walks built up step-by-step tend to get trapped after tens (for $d=2$) or hundreds ($d=3$) of steps. N.B. dimerization is a far more effective, but still exact, method for sampling SAWs. Jul 18 '16 at 3:54

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.)

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!

• Given that a close-to-reality bound on the mixing time is unlikely to available, I'd recommend starting the chain at some sort of pseudo-random polyomino rather than a fixed one. For example, make one by adding new squares at random starting with a single square, then start the chain with that. Then the chain is just fixing up an approximation that spans the whole space rather than trying to reach everywhere from one point. Jul 18 '16 at 6:05
• I agree with Brendan's point, with the caveat that when trying to infer the convergence properties of the Markov chain one may obtain a clearer signal by starting at a configuration that is far from equilibrium (e.g. straight rod or a solid ball). As an aside, the "kinetic growth" algorithm (as I would think of it) suggested can generate any polyomino, but overweights dense configurations and underweights sparse ones. Jul 18 '16 at 7:15
• Does it really satisfy detailed balance? If I am understanding the algorithm correctly, it seems like the probability of an L-tetromino transforming to square-tetromino is twice as large as a square-tetromino transforming to L-tetromino. Jul 19 '16 at 14:22
• Matthew, thanks for pointing this out, as previously described the algorithm was incorrect, now I think it should be fine. Jul 20 '16 at 6:41

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.

• Very closely related to this, but from the physics literature, is the PERM (pruned-enriched Rosenbluth method) introduced by Peter Grassberger. See: arxiv.org/abs/cond-mat/0201313 for an overview, and references therein. This will work with the kinetic growth algorithm. The unweighted method described above corresponds to the "Rosenbluth method" which typically has large variance, while PERM generates many configurations at once which allows you to automatically optimise the weights and so reduce the variance. Jul 18 '16 at 7:40
• Like sequential importance sampling, it can fail for some problems. My sense is that this happens when there is a global constraint that needs to be satisfied which cannot be easily decomposed into a sequence of local constraints as the object is built up. (I think it will work pretty well for sampling polyominoes.) Jul 18 '16 at 7:43

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.