Generate Random Latin Squares

I'm looking for algorithms to generate randomized instances of Latin squares.

I found only one paper:

M. T. Jacobson and P. Matthews, Generating uniformly distributed random Latin squares, J. Combinatorial Design 4 (1996), 405-437

Is there any other method known which generates truly random instances not just isomorphic instances of other instances?


In the combinatorics community, the Jacobson and Matthews approach is widely considered the best approach to obtain random Latin squares (approximately) uniformly at random from the set of all Latin squares. A practical non-MCMC approach that samples uniformly would be extremely well-received (but seems beyond current techniques).

Other uniform sampling methods are:

  • Generating Latin squares row-by-row by appending random permutations and restarting whenever their is a clash gives the uniform distribution. [Or equivalently, uniformly sampling from the set of row-Latin squares, then restarting if there is a clash.]
  • Generating a list of all Latin squares, and picking one at random. Storage requirements could be reduced by (a) storing only a list of normalised Latin squares [i.e. first row in order] then randomly permuting the columns after sampling, and (b) storing only the differences between subsequent Latin squares (i.e. Latin trades) [although, this makes the algorithm more complicated].

In either case, for not particularly large n [maybe n>5], these approaches are either impractically slow, or requires impractically large storage space. However, for some applications, we don't need to sample $n \times n$ Latin squares for $n>5$, in which case, this is not a problem.

Moreover, for statistical applications, sampling from all possible Latin squares is often not necessary, in which case we just apply a random isotopism to any given Latin square (i.e. we pick a Latin square, then permute its rows, columns and symbols randomly).

Any attempt to sample via extending a Latin rectangle (or partial Latin square) to a Latin square without restarting from scratch after a clash occurs will almost certainly result in a non-uniform distribution. [I suppose theoretically you could add weights to your intermediate choices.] Different Latin rectangles admit a different number of completions, so if we don't restart from scratch, we will favour Latin squares that have Latin rectangles that admit fewer completions (i.e. there's less competition for those Latin squares).

This non-uniformity might seem like a subtle difference, but consider a $(n-2) \times n$ Latin rectangle. The number of completions is always a power of 2. If the number of completions is 1 (which can happen: take a cyclic group's Cayley table and delete the last two rows), then its completion is guaranteed to be generated from that point on. If the number of completions is $2^{n/2}$ (which can also happen: take the elementary abelian 2-group's Cayley table and delete the last two rows), then the probability of it being generated from that point on could be $2^{-n/2}$ (depending on how things are implemented). So, the difference in probabilities can be at least exponential in n.

Even if you don't care too much about the uniform distribution, the Jacobson and Matthews is still reasonable: it is quite fast and simple to implement (there's also implementations for GAP ("loops") and SAGE available, and probably others I'm unaware of).


Any partial $n\times n$ Latin square with the first $r$ rows filled in ($0 \le r < n$) can be extended to a partial Latin square with the first $r+1$ rows filled in (and thus to a complete $n\times n$ Latin square). Do you know how to do this? If so, it's fairly easy to randomise the choices that have to be made. If not, let me know, and I'll try to remember how I did it 17 years ago.

Edit Suppose we have an $n \times n$ Latin square $A$ that is partially filled in: the first $r-1$ rows have been filled in, and the first $c-1$ columns in row $r$. We want to extend the Latin square to column $c$ in row $r$, i.e. we want to fill in $A_{rc}$.

If there exist any values that satisfy the Latin square constraints for $A_{rc}$, choose one of them at random. Otherwise we have to backtrack, as follows:

  1. Construct a directed graph (the 'replacement graph') which specifies which elements in row $r$ may be replaced by which elements. Its vertices are:

    • a start vertex $S$, representing $A_{rc}$
    • an end vertex $E$, representing values unused in the current row
    • for each element $A_{ri}$ in the current row with $i<c$, a vertex $V_i$
  2. Add edges as follows:

    • for each vertex $V_i$, and for each value $v$ unused in column $i$, add an edge from $V_i$ to $V_j$ if $A_{rj} = v$ for some $j$, else add an edge from $V_i$ to $E$.
    • for each value $v$ unused in column $c$, add an edge from $S$ to $V_j$ if $A_{rj} = v$.
  3. Select a path $P$ from $S$ to $E$ at random.

  4. Replace the values in row $r$ as follows:

    • for the first edge in $P$, from $S$ to $V_i$: set $A_{rc} = A_{ri}$
    • for internal edges in $P$, from $V_i$ to $V_j$: set $A_{ri} = A_{rj}$
    • for the last edge in $P$, from $V_i$ to $E$: set $A_{ri}=$any unused value, chosen at random

This algorithm runs quickly, and you won't have any trouble generating $12 \times 12$ Latin squares with it. The problem remains that it probably doesn't generate all Latin squares with equal probability. In particular, it's not clear how to select the path $P$ at random.

Perhaps I should hand-run this algorithm on my counterexample, to show you how it works. We have the partial $4\times4$ Latin square

0 1 2 3
1 3 0 2
3 0 1

The replacement graph has vertices $S$, $E$, $V_1$, $V_2$, $V_3$. Its edges are:

  • $V_1 \rightarrow E$ (because $2$ is unused)
  • $V_2 \rightarrow E$ (because $2$ is unused)
  • $V_3 \rightarrow V_1$ (because $A_{r1}=3$)
  • $S \rightarrow V_1$ (because $A_{r1}=3$)
  • $S \rightarrow V_2$ (because $A_{r2}=0$)

Pick the path $P = S \rightarrow V_2 \rightarrow V_0 \rightarrow E$. Then we set:

  • $A_{3,4} = A_{3,2}$
  • $A_{3,2} = A_{3,0}$
  • $A_{3,0} = $ the only unused value, $2$

And we get:

0 1 2 3
1 3 0 2
2 0 3 1