Efficient uniform generation of random derangements with the expected distribution of cycle lengths

1 Introduction

Derangements are permutations $σ = σ_{1} \dots σ_{n}$ on $n \geq 2$ labels such that $σ_{i} \neq i$ for all $i = 1, \dots$ , $n$ . Derangements are useful in a number of applications like in the testing of software branch instructions and random paths and data randomization and experimental design bacher ; pd-rlg-sph ; sedgewick . A well known algorithm to generate random derangements is Sattolo’s algorithm, that outputs a random cyclic derangement in $O (n)$ time gries ; prodinger ; sattolo ; wilson . An $O (2 n)$ algorithm to generate random derangements in general (not only cyclic derangements) has been given in analco ; iran . Algorithms to generate all $n$ -derangements in lexicographic or Gray order have also been developed akl ; baril ; korsh .

In this paper we propose two procedures to generate random derangements with the expected distribution of cycle lengths: one based on the randomization of derangements by random restricted transpositions (a random walk in the set of derangements) and the other based on a simple sequential importance sampling scheme. The generation of restricted permutations by means of sequential importance sampling methods is closely related with the problem of estimating the permanent of a matrix, an important problem in, e. g., graph theory, statistical mechanics, and experimental design

beichl ; brualdi ; pd-rlg-sph . Simulations show that the randomization algorithm samples a derangement uniformly in

O (n^{a} log n^{2})

time, where

n

is the size of the derangement and

a ≃ \frac{1}{2}

, while the sequential importance sampling algorithm does it in

O (n)

time but with a small probability

O (1 / n)

of failing. The algorithms are straighforward to understand and implement and can be modified to perform related computations of interest in many areas.

2 Mathematical preliminaries

Let us briefly recapitulate some notation and terminology on permutations; for detailed accounts see arratia ; charalambides .

We denote a permutation of a set of $n \geq 2$ labels (an $n$ -permutation), formally a bijection of $[n] = {1, \dots, n} \subset N$ onto itself, by $σ = σ_{1} \dots σ_{n}$ , where $σ_{i} = σ (i)$ . If $σ$ and $π$ are two $n$ -permutations, their product is given by the composition $σ π = σ (π_{1}) \dots σ (π_{n})$ . A cycle of length $k \leq n$ in a $n$ -permutation $σ$ is a sequence of indices $i_{1}, \dots, i_{k}$ such that $σ (i_{1}) = i_{2}$ , …, $σ (i_{k - 1}) = i_{k}$ , and $σ (i_{k}) = i_{1}$ , completing the cycle. Fixed points are $1$ -cycles, transpositions are $2$ -cycles. An $n$ -permutation with $a_{k}$ cycles of length $k$ , $1 \leq k \leq n$ , is said to be of type $(a_{1}, \dots, a_{n})$ , with $\sum_{k} k a_{k} = n$ . For example, the $9$ -permutation $174326985 = (8) (6) (34) (2795) (1)$ has $5$ cycles and is of type $(3, 1, 0, 1)$ , where we have omitted the trailing $a_{5} = \dots = a_{9} = 0$ .

The number of $n$ -permutations with $k$ cycles is given by the unsigned Stirling number of the first kind $[\frac{n}{k}]$ . Useful formulae involving these numbers are $[\frac{0}{0}] = 1$ , $[\frac{n}{0}] = 0$ , and the recursion relation $[\frac{n + 1}{k}] = n [\frac{n}{k}] + [\frac{n}{k - 1}]$ . We have $[\frac{n}{n}] = 1$ , counting just the identity permutation $id = (1) (2) \dots (n)$ , $[\frac{n}{n - 1}] = (\frac{n}{2})$ , counting $n$ -permutations with $n - 2$ fixed points, that can be taken in $(\frac{n}{n - 2}) = (\frac{n}{2})$ different ways, plus a transposition of the remaining two labels, and $[\frac{n}{1}] = (n - 1)!$ , the number of cyclic $n$ -derangements. It can also be shown that $[\frac{n}{2}] = (n - 1)! H_{n - 1}$ , where $H_{k} = 1 + \frac{1}{2} + \dots + \frac{1}{k}$ is the $k$ th harmonic number. Obviously, $[\frac{n}{1}] + \dots + [\frac{n}{n}] = [\frac{n + 1}{1}] = n!$ , the total number of $n$ -permutations.

Let us denote the set of all $n$ -derangements by $D_{n}$ . It is well known that

d_{n} = | D_{n} | = n! (1 - \frac{1}{1!} + \dots + \frac{(- 1)^{n}}{} n!) = ⌊ \frac{n! + 1}{e} ⌋, n \geq 1,

(1)

the so-called rencontres (or subfactorial, $! n$ ) numbers. Let us also denote the set of $k$ -cycle $n$ -derangements, irrespective of their type, by $D_{n}^{(k)}$ . Note that $D_{n}^{(k)} = \emptyset$ for $k > ⌊ n / 2 ⌋$ . If we want to generate random $n$ -derangements uniformly over $D_{n} = D_{n}^{(1)} \cup \dots \cup D_{n}^{(⌊ n / 2 ⌋)}$ , we must be able to generate $k$ -cycle random $n$ -derangements with probabilities

P (σ \in D_{n}^{(k)}) = \frac{d_{n}^{(k)}}{d_{n}},

(2)

where $d_{n}^{(k)} = | D_{n}^{(k)} |$ . The following proposition establishes the (little known, hard-to-find) cardinality of the sets $D_{n}^{(k)}$ .

Proposition 2.1.

The cardinality of the set $D_{n}^{(k)}$ is given by

d_{n}^{(k)} = k \sum j = 0 (- 1)^{j} (\frac{n}{j}) [\frac{n - j}{k - j}] .

(3)

Proof.

The number of $n$ -permutations with $k$ cycles is $[\frac{n}{k}]$ . Of these, $n [\frac{n - 1}{k - 1}]$ have at least one fixed point, $(\frac{n}{2}) [\frac{n - 2}{k - 2}]$ have at least two fixed points, and so on. Perusal of the inclusion-exclusion principle then furnishes the result. ∎

Proposition 2.2.

The numbers $d_{n}^{(k)}$ obey the recursion relation

d_{n + 1}^{(k)} = n (d_{n}^{(k)} + d_{n - 1}^{(k - 1)})

(4)

with $d_{0}^{(0)} = 1$ and $d_{n}^{(0)} = 0$ , $n \geq 1$ .

Proposition 2.2 follows from (3) by index manipulations and the properties of the binomial coefficients and Stirling numbers of the first kind. The numbers $d_{n}^{(k)}$ are sometimes called associated Stirling number of the first kind. Equation (4) generalizes the recursion relation $d_{n + 1} = n (d_{n} + d_{n - 1})$ for the rencontres numbers. Equation (3) recovers $d_{n}^{(0)} = 0$ and $d_{n}^{(1)} = [\frac{n}{1}] = (n - 1)!$ for $n \geq 1$ , while we find that $d_{n}^{(2)} = (n - 1)! (H_{n - 2} - 1)$ for $n \geq 2$ . Another noteworthy identity, valid for $n$ even, is $d_{n}^{(n / 2)} = (n / 2 - 1)!!$ , the number of derangements $σ$ with the property that $σ^{2} = id$ , a. k. a. fixed-point-free involutions. We see that already for small $n$ we have $P (σ \in D_{n}^{(1)}) ≃ e / n$ and $P (σ \in D_{n}^{(2)}) ≃ (H_{n - 2} - 1) e / n$ .

Remark 2.3.

One can consider the distribution of $n$ -derangements over possible cycle types (instead of cycle lengths) for a “finer” view of the distribution. The number of $n$ -permutations of type $(a_{1}, \dots, a_{n})$ is given by Cauchy’s formula

k_{n} (a_{1}, \dots, a_{n}) = \frac{n!}{1^{a_{1}} a_{1}! \dots n^{a_{n}} a_{n}!} .

(5)

The analogue of (2) is given by $P (σ \in K_{n} (0, a_{2}, \dots, a_{n})) = k_{n} (0, a_{2}, \dots, a_{n}) / d_{n}$ , where $K_{n} (0, a_{2}, \dots, a_{n})$ is the conjugacy class formed by all $n$ -permutations of type $(0, a_{2}, \dots, a_{n})$ . For example, for cyclic derangements $a_{n} = 1$ and all other $a_{k} = 0$ , and we obtain $P (σ \in K_{n} (0, \dots, 1)) = k_{n} (0, \dots, 1) / d_{n} ≃ e / n$ , as before.

3 Generating random derangements by random transpositions

Our first approach to generate random $n$

-derangements uniformly distributed over

D_{n}

consists in taking an initial

n

-derangement and to scramble it by random restricted transpositions enough to obtain a sample from the uniform distribution over

D_{n}

. By restricted transpositions we mean swaps

σ_{i} \leftrightarrow σ_{j}

avoiding pairs for which

σ_{i} = j

σ_{j} = i

. Algorithm T describes the generation of random

n

-derangements according to this idea, where

m i x

is a constant establishing the amount of random restricted transpositions to be attempted and

r n d

is a computer generated pseudorandom uniform deviate in

(0, 1)

Remark 3.1.

Algorithm T is applicable only for $n \geq 4$ , since it is not possible to connect the even permutations $231$ and $312$ by a single transposition.

A good choice for the initial derangement in Algorithm T is any cyclic derangement (cycle length $k = 1$ ), for example, $σ = (23 \dots n 1)$ . A particularly bad choice would be an involution ( $n$ even, cycle length $k = n / 2$ ), for example, $σ = (n n - 1) \dots (21)$ , because then the algorithm would never be able to generate derangements with $k \neq n / 2$ . (Incidentally, this suggests the use of Algorithm T to generate random fixed-point-free involutions.) To avoid this problem we hardcoded the requirement that Algorithm T starts with a cyclic derangement. If several parallel streams of random derangements are needed, one can set different initial random cyclic derangements from a one-line implementation of Sattolo’s algorithm.

Remark 3.2.

The minimum number of transpositions necessary to take a cyclic $n$ -derangement into a $k$ -cycle $n$ -derangement is $k - 1$ , $1 \leq k \leq ⌊ n / 2 ⌋$ , since transpositions of labels that belong to the same cycle split it into two cycles,

(a b) (i_{1} \dots i_{a - 1} i_{a} i_{a + 1} \dots i_{b - 1} i_{b} i_{b + 1} \dots i_{k}) = (i_{1} \dots i_{a - 1} i_{b} i_{b + 1} \dots i_{k}) (i_{a + 1} \dots i_{b - 1} i_{a})

(6)

and, conversely, transpositions involving labels of different cycles join them into a single cycle. If Algorithm T is started with a cyclic derangement then one must set $m i x \geq n / 2$ .

0: Initial cyclic

n

-derangement

σ_{1} σ_{2} \dots σ_{n}

m i x \leftarrow

number of restricted transpositions to attempt

2: for

m = 1

m i x

i \leftarrow ⌈ r n d \cdot n ⌉

j \leftarrow ⌈ r n d \cdot n ⌉

4: if

(σ_{i} \neq j) \land (σ_{j} \neq i)

then

5: swap

σ_{i} \leftrightarrow σ_{j}

6: end if

7: end for

7: For sufficiently large

m i x

σ_{1} \dots σ_{n}

is a random derangement uniformly distributed in

D_{n}

Algorithm T Random derangements by random restricted transpositions

We run Algorithm T for $n = 64$ and different values of $m i x \geq n$ and collect data. Simulations were performed on Intel Xeon E5-1650 v3 processors running -O3 compiler-optimized C code (GCC v. $5.4.0$ ) over Linux kernel $4.15.0$ at $3.50$ GHz, while the numbers (3) were calculated on the software package Mathematica 11.3 wolfram . We draw our pseudo-random numbers from Vigna’s superb xoshiro256+ generator xoshiro . Our results appear in Table 1. We see from that table that with $m i x = n$ random restricted transpositions there is a slight excess of probability mass in the lower $k$ -cycle sets with $k = 1, 2$ , and $3$ . Trying to scramble the initial $n$ -derrangement by $2 n$ restricted transpositions performs better. The difference between attempting $2 n$ and $n log n$ random restricted transpositions is much less pronounced. Figures for derangements of higher cycle number fluctuate more due to the finite size of the sample. The data suggest that Algorithm T can generate a random $n$ -derangement uniformly distributed on $D_{n}$ with $2 n$ random restricted transpositions, employing $4 n$ pseudorandom numbers in the process. This is further discussed in Section 5.

Remark 3.3.

It is a classic result that $O (n log n)$ transpositions are needed before a shuffle by unrestricted transpositions becomes “sufficiently random” aldous ; persi . A similar analysis for random restricted transpositions over derangements is complicated by the fact that derangements do not form a group. Recently, the analysis of the spectral gap of the Markov transition kernel of the process provided the upper bound $m i x < C n + a n log n^{2}$ , with $a > 0$ and $C \geq 0$ a decreasing function of $n$ aaron . This bound results from involved estimations and approximations and may not be very accurate. We are not aware of other rigorous results on this particular problem. Related results for the mixing time of the random transposition walk on permutations with “one-sided restrictions” ( $σ_{i} \geq b_{i}$ for given $b_{i} \geq 1$ ) appear in hanlon .

Cycles	Algorithm T ( $m i x$ )			Algorithm S	Exact
$k$	$n$	$2 n$	$n log n$	—	Eqs. (1)–(3)
$1$	$0.042 933$	$0.042 479$	$0.042 473$	$0.042 475$	$0.042 473$
$2$	$0.158 395$	$0.157 691$	$0.157 679$	$0.157 684$	$0.157 677$
$3$	$0.260 129$	$0.258 787$	$0.258 765$	$0.258 788$	$0.258 772$
$4$	$0.252 739$	$0.253 304$	$0.253 305$	$0.253 306$	$0.253 301$
$5$	$0.167 189$	$0.167 621$	$0.167 639$	$0.167 622$	$0.167 635$
$6$	$0.079 498$	$0.080 390$	$0.080 402$	$0.080 389$	$0.080 400$
$7$	$0.028 825$	$0.029 192$	$0.029 195$	$0.029 196$	$0.029 200$
$8$	$0.008 087$	$0.008 269$	$0.008 274$	$0.008 272$	$0.008 274$
$9$	$0.001 821$	$0.001 868$	$0.001 869$	$0.001 868$	$0.001 869$
$10$	${3.292}_{- 4}$	${3.416}_{- 4}$	${3.418}_{- 4}$	${3.412}_{- 4}$	${3.417}_{- 4}$
$11$	${4.914}_{- 5}$	${5.109}_{- 5}$	${5.120}_{- 5}$	${5.103}_{- 5}$	${5.116}_{- 5}$
$12$	${5.997}_{- 6}$	${6.322}_{- 6}$	${6.301}_{- 6}$	${6.354}_{- 6}$	${6.326}_{- 6}$
$13$	${6.215}_{- 7}$	${6.493}_{- 7}$	${6.301}_{- 7}$	${6.507}_{- 7}$	${6.499}_{- 7}$
$14$	${4.83}_{- 8}$	${5.40}_{- 8}$	${5.57}_{- 8}$	${5.44}_{- 8}$	${5.569}_{- 8}$
$15$	${4.6}_{- 9}$	${3.1}_{- 9}$	${3.0}_{- 9}$	${4.1}_{- 9}$	${3.989}_{- 9}$
$16$	$4_{- 10}$	$1_{- 10}$	$3_{- 10}$	$1_{- 10}$	${2.390}_{- 10}$
runtime (sec)	$11298$	$15068$	$25238$	$18954$	—

Table 1: Proportion of

n

-derangements in

D_{n}^{(k)}

measured in

10^{10}

samples generated by Algorithms T and S for

n = 64

. The notation

x_{- a}

reads

x \times 10^{- a}

. Data for Algorithm S are based on a run with a ratio of completed/attempted derangements of

0.985472

. Runtimes are given for comparison.

4 Sequential importance sampling of derangements

4.1 The SIS algorithm

Sequential importance sampling (SIS) is an importance sampling scheme with the sampling weights built up sequentially. The idea is particularly suited to sample composite objects $X = X_{1} \dots X_{n}$ from a complicated sample space $X$ with many restrictions and for which the high-dimensional volume $| X |$ , from which the uniform distribution $P (X) = | X |^{- 1}$ follows, may not be easily calculable. However, since we can always write

P (X_{1} \dots X_{n}) = P (X_{1}) P (X_{2} ∣ X_{1}) \dots P (X_{n} ∣ X_{1} \dots X_{n - 1}) .

(7)

we can think of “telescoping” the sampling of $X$ by first sampling $X_{1}$ , then use the updated information brought by the knowledge of $X_{1}$ to sample $X_{2}$ and so on. In Monte Carlo simulations, the right-hand side of (7) actually becomes $P_{1} (X_{1}) P_{2} (X_{2} ∣ X_{1}) \dots$ $P_{n} (X_{n} ∣ X_{1} \dots X_{n - 1})$ , with the distributions $P_{i} (\cdot)$ estimated or inferred incrementally based on approximate weighting functions for the partial objects $X_{1} \dots X_{i - 1}$ . Expositions of the SIS framework of interest to what follows appear in cdhl ; pd-rlg-sph .

Algorithm S

describes a SIS algorithm to generate random derangements inspired by the analogous problem of sampling contingency tables with restrictions

cdhl ; pd-rlg-sph as well as by the problem of estimating the permanent of a matrix kuznetsov ; liang ; rasmussen . Our presentation of Algorithm S is not the most efficient for implementation; the auxiliary sets

J_{i}

, for instance, are not actually needed and were included only to facilitate the analysis of the algorithm, and the

n

tests in line 4 can be reduced to a single test in the last pass, since all

J_{i} \neq \emptyset

except perhaps

J_{n}

The distribution of cycle lengths in $10^{10}$ derangements generated by Algorithm S is presented in Table 1; there is excellent agreement between the data and the exact values.

J \leftarrow [n]

2: for

i = 1

n

J_{i} \leftarrow J ∖ {i}

4: if

J_{i} \neq \emptyset

then

5: choose

j_{i} \in J_{i}

uniformly at random

σ_{i} \leftarrow j_{i}

J \leftarrow J ∖ {j_{i}}

8: else

9: fail

10: end if

11: end for

11: If completed,

σ_{1} \dots σ_{n}

is a random derangement uniformly distributed in

D_{n}

Algorithm S Random derangements by sequential importance sampling

4.2 Failure probability of the SIS algorithm

In the $i$ th pass of the loop in Algorithm S, $σ_{i}$ can pick (lines 5–6) one of either $n - i$ or $n - i + 1$ labels, depending on whether label $i$ has already been picked. This guarantees the construction of the $n$ -derangement till the last but one element $σ_{n - 1}$ . The $n$ -derangement is completed only if the last remaining label is different from $n$ , such that $σ_{n}$ does not pick $n$ . The probability that Algorithm S fails is thus given by

P (σ_{n} = n ∣ σ_{1} \dots σ_{n - 1}) = P (σ_{1} \neq n) P (σ_{2} \neq n ∣ σ_{1}) \dots P (σ_{n - 1} \neq n ∣ σ_{1} \dots σ_{n - 2}) .

(8)

Now, according to Algorithm S, line 5, we have

P (σ_{i} \neq n ∣ σ_{1} \dots σ_{i - 1}) = 1 - P (σ_{i} = n ∣ σ_{1} \dots σ_{i - 1}) = 1 - \frac{1}{E (J_{i} ∣ σ_{1} \dots σ_{i - 1})},

(9)

where we write $J_{i}$ for $| J_{i} (σ_{1} \dots σ_{i - 1}) |$ in the argument of the expectation to improve readability, and the failure probability becomes

(10)

where $E_{i}$ stand for $E (J_{i} ∣ σ_{1} \dots σ_{i - 1})$ . Expression (10) is but a probabilistic version of the inclusion-exclusion principle known as Poincaré’s formula in elementary probability.

The explicit computation of (10) is a cumbersome business, and we will not pursue it here. Otherwise, the following theorem establishes an easy upper bound on the failure probability of Algorithm S.

Theorem 4.1.

Algorithm S fails with probability $O (1 / n)$ .

Proof.

Recall that in the $i$ th pass of the loop in Algorithm S we have $J_{i} = | J_{i} (σ_{1} \dots σ_{i - 1}) | = n - i$ or $n - i + 1$ , such that the expectation $E_{i} = E (J_{i} ∣ σ_{1} \dots σ_{i - 1})$ obeys

1 - \frac{1}{n - i} < 1 - \frac{1}{E_{i}} < 1 - \frac{1}{n - i + 1}

(11)

and it immediately follows that

(12)

∎

Remark 4.2.

In A we provide an improved upper bound to $P (σ_{n} = n ∣ σ_{1} \dots σ_{n - 1})$ based on an “independent sets” approximation that works remarkably well.

The measured failure rate for the SIS data in Table 1 is $1 - 0.985472 = 0.014528$ , not far from $1 / 64 = 0.015625$ . A sample of $10^{4}$ runs of Algorithm S of $10^{6}$ derangements each with $n = 64$ gives an average failure rate of $0.01453 (12)$

with a sample minimum of 0.014130 and maximum of 0.014991, where the digits within parentheses indicate the uncertainty at one standard deviation in the corresponding last digits of the datum. Figure

1 depicts Monte Carlo data for the failure probability (10) against the upper bound

1 / n

. Each data point was obtained as an average over

10^{4}

runs of Algorithm S of

10^{6}

derangements each except for

n = 512

, for which the runs are of

2 \times 10^{5}

derangements each.

Figure 1: Measured failure rate for Algorithm S against the upper bound $1 / n$ . Error bars in the data are much smaller than the symbols shown.

4.3 Uniformity of the SIS algorithm

In the SIS approach, the ensuing sampling probabilities may deviate considerably from the uniform distribution. An important aspect of Algorithm S is that it generates each $σ \in D_{n}$ with probability $1 / d_{n}$ , as the data in Table 1 suggest. We prove that this is indeed the case based on related results in cdhl ; kuznetsov ; rasmussen .

Definition 4.3.

We call the $n \times n$ matrix $A_{n} = (a_{i j})$ with $0$ on the diagonal and $1$ elsewhere the derangement matrix of order $n$ .

The following well known property of $A_{n}$ clarifies its denomination and will be useful later.

Lemma 4.4.

The permanent of $A_{n}$ is given by the $n$ th rencontre number $d_{n}$ .

Proof.

The result is classic and we only outline its proof; for details consult brualdi . The permanent of $A_{n}$ is given by

per A_{n} = \sum σ n \prod i = 1 a_{i σ_{i}} = n \sum i = 1 a_{i j} per A_{n - 1}^{(i j)},

(13)

where the first sum is over all $n!$ permutations $σ$ of $[n]$ and the other sum is the Laplace (row) development of the permanent in terms of the submatrices $A_{n - 1}^{(i j)}$ obtained from $A_{n}$ by deleting its $i$ th row and $j$ th column. Perusal of the Laplace development of $per A_{n}$ twice eventually leads to the recursion relation $per A_{n} = (n - 1) (per A_{n - 1} + per A_{n - 2})$ with $per A_{1} = 0$ and $per A_{2} = 1$ . But this is exactly the recursion relation for the rencontre numbers. ∎

The next lemma elucidates the close relationship between the sets $J_{i}$ of Algorithm S and $A_{n}$ .

Lemma 4.5.

The product of the $E_{i} = E (J_{i} ∣ σ_{1} \dots σ_{i - 1})$ , $1 \leq i \leq n$ , generated by Algorithm S

is a one-sample unbiased estimator for

per A_{n}

Proof.

We have to prove that $E (E_{1} \dots E_{n}) = E (J_{1}) E (J_{2} ∣ σ_{1}) \dots$ $E (J_{n} ∣ σ_{1} \dots σ_{n - 1}) = per A_{n}$ . From Algorithm S, lines 3 and 7, it follows, with $J_{i} = [n] ∖ {i, j_{1}, \dots, j_{i - 1}}$ for some realization of ${j_{1}, \dots, j_{i - 1}}$ , that

E (J_{i} ∣ σ_{1} \dots σ_{i - 1}) = n \sum j_{i} = 1 1{1}% {j_{i} \in J_{i}} = \sum j_{i} \neq j_{1}, \dots, j_{i - 1} 1{1} {j_{i} \in [n] ∖ {i}},

(14)

where $1{1} {P}$ is the indicator function that equals $1$ if $P$ is true and $0$ otherwise. Now if we write $a_{i j_{i}}$ for $1{1} {j_{i} \in [n] ∖ {i}}$ and multiply the expectations (14) we find that

E (E_{1} \dots E_{n}) = \sum j_{1} a_{1 j_{1}} \sum j_{2} \neq j_{1} a_{2 j_{2}} \dots \sum j_{n} \neq j_{1}, \dots, j_{n - 1} a_{n j_{n}} .

(15)

Expression (15) is nothing but the expression for the permanent of a matrix with elements $a_{i j} = 0$ if $j = i$ and $1$ otherwise, that is, of the derangement matrix $A_{n}$ , and by Lemma 4.4 it follows that $E (E_{1} \dots E_{n}) = per A_{n}$ . ∎

Theorem 4.6.

Algorithm S samples $D_{n}$ uniformly.

Proof.

Algorithm S generates permutations $σ$ with probability $P (σ) = P (σ_{1}) P (σ_{2} ∣ σ_{1}) \dots$ $P (σ_{n} ∣ σ_{1} \dots σ_{n - 1})$ . Since

P (σ_{i} ∣ σ_{1} \dots σ_{i - 1}) = \frac{1}{E (J_{i} ∣ σ_{1} \dots σ_{i - 1})}

(16)

we have

P (σ) = \frac{1}{E (J_{1})} \frac{1}{E (J_{2} ∣ σ_{1})} \dots \frac{1}{E (J_{n} ∣ σ_{1} \dots σ_{n - 1})} .

(17)

By Lemmas 4.4 and 4.5 the r. h. s. of (17) equals $1 / per A_{n} = 1 / d_{n}$ , and we are golden. ∎

5 Mixing time of the restricted transpositions shuffle

To shed some light on the question of how many random restricted transpositions are necessary to generate random derangements uniformly over $D_{n}$ , we investigate the convergence of Algorithm T numerically. This can be done by monitoring the evolution of the empirical probabilities along the run of the algorithm towards the exact probabilities given by (2).

Let $ν$ be the measure that puts mass $ν (k) = d_{n}^{(k)} / d_{n}$ on the set $D_{n}^{(k)}$ and $μ_{t}$ be the empirical measure

μ_{t} (k) = \frac{1}{t} t \sum s = 1 1{1} {σ_{s} \in D_{n}^{(k)}},

(18)

where $σ_{s}$ is the derangement obtained after attempting $s$ restricted transpositions by Algorithm T on a given initial derangement $σ_{0}$ . The total variance distance between $μ_{t}$ and $ν$ is given by aldous ; persi

d_{T V} (t) = ∥ μ_{t} - ν ∥_{T V} = \frac{1}{2} ⌊ n / 2 ⌋ \sum k = 1 | μ_{t} (k) - ν (k) | .

(19)

The right-hand side of (19) can be seen as the “histogram distance” between $μ_{t}$ and $ν$ in the $ℓ_{1}$ norm. Clearly, $0 \leq d_{T V} (t) \leq 1$ . This distance allows us to define $t_{m i x} (ϵ)$ as the time it takes for $μ_{t}$ to fall within distance $ϵ$ of $ν$ ,

t_{m i x} (ϵ) = min {t \geq 0 : d_{T V} (t) < ϵ} .

(20)

It is usual to define the mixing time $t_{m i x}$ by setting $ϵ = \frac{1}{4}$ or $ϵ = \frac{1}{2} e^{- 1} ≃ 0.184$

, this last figure being reminiscent of the spectral analysis of Markov chains. We set

ϵ = \frac{1}{2} e^{- 1}

. This choice is motivated by the following pragmatic reasons:

We want the derangements output by Algorithm T to be as uniformly distributed over $D_{n}$ as possible, so the smaller the $ϵ$ the better the assessment of the algorithm and the choice of the constant $m i x$ ;
Most of the probability mass is concentrated on a few cycle numbers (see Table 1 and Remark 5.1 below), such that even relatively small differences between $μ_{t}$ and $ν$ are likely to induce noticeable biases in the output of Algorithm T;
With $ϵ = \frac{1}{4}$ we found that $t_{m i x} < n / 2$ , meaning that not even every possible derangement had chance to be generated if the initial derangement is cyclic (see Remark 3.2).

Remark 5.1.

It is well known that the number of $k$ -cycles of random $n$

-permutations is Poisson distributed with mean

1 / k

, such that as

n ↗ \infty

the CLT implies that the length of the cycles of random permutations follow a normal ditribution with mean

log n

and variance

log n

; see, e. g., arratia and the references there in. soria proved that the same holds for permutations with no cycles of length less than a given

ℓ > 1

, using complex asymptotics of exponential generating functions; analco and iran provide the analysis for the particular case of derangements. Figure 2 displays the exact distribution of

k

-cycles for derangements with

n = 32768

together with the normal density

N (log n, \sqrt{log n})

. For

n = 32768

we obtain from equations (2)–(3) that

⟨ k ⟩ = 9.967 \dots

and

\sqrt{⟨ k^{2} ⟩ - ⟨ k ⟩^{2}} = 2.872 \dots

, while

log n = 10.397 \dots

and

\sqrt{log n} = 3.224 \dots

. The distribution of cycle lengths in Figure 2 indeed looks close to a normal

N (log n, \sqrt{log n})

, albeit slightly skewed. We did not go to greater

n

because Stirling numbers of the first kind are notoriously hard to compute even by computer algebra systems running on modern workstations. Recently, the cycle structure of certain types of restricted permutations (with

σ_{i} \geq i - 1

) was also shown to be asymptotically normal ozel .

Figure 2: Distribution of cycle lengths of random $n$ -derangements for $n = 32768$ together with the normal densities $N (log n, \sqrt{log n})$ (shorter, in blue) and $N (m, s)$ (taller, in orange) with $m = ⟨ k ⟩ = 9.967 \dots$ and $s = \sqrt{⟨ k^{2} ⟩ - ⟨ k ⟩^{2}} = 2.872 \dots$ .

Starting with a cyclic derangement, i. e., with $μ_{0} (1) = 1$ and all other $μ_{0} (k) = 0$ , we run Algorithm T and collect statistics on $d_{T V} (t)$ . Figure 3 displays the average $⟨ d_{T V} (t) ⟩$ over $10^{6}$ runs for $n = 128$ . The behavior of $⟨ d_{T V} (t) ⟩$ does not show sign of the cutoff phenomenon—a sharp transition from unmixed state ( $d_{T V} (t_{m i x} - δ) \approx 1$ ) to mixed state ( $d_{T V} (t_{m i x} + δ) \approx 0$ ) over a small window of time $δ ≪ t_{m i x}$ . Table 2 lists the average $⟨ t_{m i x} ⟩$ obtained over $10^{6}$ samples for larger derangements at $ϵ = \frac{1}{2} e^{- 1}$ . An adjustment of the data to the form

t_{m i x} = c n^{a} log n^{2}

(21)

furnishes $a = 0.527 (2)$ and $c = 0.90 (1)$ . Our data thus suggest that $t_{m i x} \sim O (n^{a} log n^{2})$ with $a ≃ \frac{1}{2}$ , roughly an $O (\sqrt{n})$ lower than the upper bound given by aaron . It is tempting to conjecture that $a = \frac{1}{2}$ (and, perhaps, that $c = 1$ ) exactly, cf. last two lines of Table 2, although our data do not support the case unequivocally.

Figure 3: Total variance distance $⟨ d_{T V} (t) ⟩$ (averaged over $10^{6}$ runs) between the empirical measure $μ_{t}$ (with $μ_{0} (1) = 1$ ) and the stationary measure $ν$ of the process defined by Algorithm T for $n = 128$ . The dashed line indicates the level $ϵ = \frac{1}{2} e^{- 1}$ .

$n$	$64$	$128$	$192$	$256$	$320$	$384$	$448$	$512$
$⟨ t_{m i x} ⟩$	$67$	$112$	$150$	$184$	$216$	$245$	$274$	$301$
$\sqrt{n} log n^{2}$	$67$	$110$	$146$	$177$	$206$	$233$	$258$	$282$
$a$ in $n^{a} log n^{2}$	$0.502$	$0.504$	$0.505$	$0.507$	$0.508$	$0.508$	$0.510$	$0.510$

Table 2: Mixing time

t_{m i x}

evaluated at

ϵ = \frac{1}{2} e^{- 1}

obtained from an average trajectory

⟨ d_{T V} (t) ⟩

over

10^{6}

samples; see Figure 3. The second line displays the best guess to

n^{a} log n^{b}

involving only integer and semi-integer exponents. The last line displays the adjusted value of

a

supposing a dependence like in (21) with

c = 1

6 Summary and conclusions

While simple rejection-sampling generates random derangements with an acceptance rate of $\sim e^{- 1}$ $≃ 0.368$ , thus being $O (e \cdot n)$ (plus the cost of verifying if the permutation generated is a derangement, which does not impact the complexity of the algorithm but impacts its runtime), Sattolo’s $O (n)$ algorithm only generates cyclic derangements, and Martínez-Panholzer-Prodinger algorithm, with guaranteed uniformity, is $2 n + O ({log}^{2} n)$ , we described two procedures, Algorithms T and S, that are competitive for the efficient uniform generation of random derangements. We note in passing that Algorithm T can easily be used also to generate random fixed-point-free involutions.

We found, empirically, that $O (n^{a} log n^{2})$ random restricted transpositions with $a ≃ 0.5$ suffice to spread an initial $n$ -derangement over $D_{n}$ uniformly. The fact that $2 n > c n^{a} log n^{2}$ for all $n \geq 1$ as long as $a \leq 0.63$ and $c \leq 1$ explains the good statistics (see Table 1) displayed by Algorithm T with $m i x = 2 n$ . There are currently very few analytical results on the mixing time of the random restricted transposition walk implemented by Algorithm T; the upper bound $O (n log n^{2})$ obtained in aaron is roughly $O (\sqrt{n})$ above our empirical findings.

Algorithm T employs $2 m i x$ pseudorandom numbers and Algorithm S employs $O (n)$ pseudorandom numbers to generate an $n$ -derangement uniformly distributed over $D_{n}$ . In this way, even if we set $m i x = c \sqrt{n} log n^{2}$ with some $1 < c \sim O (1)$ , both algorithms perform better than currently known methods, with comparable runtime performances between them.

Acknowledgments

The author thanks Aaron Smith (U. Ottawa) for useful correspondence and suggestions on a previous version of the manuscript, the LPTMS for kind hospitality during a sabbatical leave in France, and FAPESP (Brazil) for partial support through grants nr. 2015/21580-0 and 2017/22166-9.

Appendix A Improved upper bound for the failure probability of Algorithm S

In the $i$ th pass of the loop in Algorithm S we have

| J_{i} (σ_{1} \dots σ_{i - 1}) | = n - i + i - 1 \sum j = 1 % 1{1} {(σ_{j} = i ∣ σ_{1}, \dots, σ_{j - 1})} .

(22)

The difficulty in the calculation of $E_{i} = E (J_{i} ∣ σ_{1} \dots σ_{i - 1})$ resides in the calculation of $E (1{1} {(σ_{j} = i ∣ σ_{1}, \dots, σ_{j - 1})})$ . We can simplify this calculation by ignoring the conditioning of the event $(σ_{j} = i)$ on the event $(σ_{1} \dots σ_{j - 1})$ , i. e., by ignoring correlations between the values assumed by the many $σ_{j}$ along a “path” in the algorithm. This approximation is clearly better in the beginning of the construction of $σ$ , when $j$ is small, than later. In doing so we obtain

E (i - 1 \sum j = 1 1{1} {(σ_{j} = i ∣ σ_{1}, \dots, σ_{j - 1})}) = i - 1 \sum j = 1 E (1{1} {(σ_{j} = i ∣ σ_{1}, \dots, σ_{j - 1})}) \approx i - 1 \sum j = 1 E (1{1} {σ_{j} = i}) = \frac{i - 1}{n - 1},

(23)

such that $E_{i} \approx n - i + \frac{i - 1}{n - 1}$ . This approximate $E_{i}$ is greater than the true $E_{i}$ , because conditioning $J_{i}$ on $σ_{1} \dots σ_{i - 1}$ can only restrict, not enlarge, the set of indices available to $σ_{i}$ . Numerical evidence supports this claim. We thus have that the approximate value of $1 - 1 / E_{i}$ is greater than its true value, and the failure probability (10) becomes bounded by

P (σ_{n} = n ∣ σ_{1} \dots σ_{n - 1}) < n - 1 \prod i = 1 (1 - \frac{1}{n - i + \frac{i - 1}{n - 1}}) = \frac{1}{n - 1} n - 1 \prod i = 1 [1 + \frac{1}{(n - 2) (n - i)}]^{- 1} .

(24)

The surprisingly excellent agreement between this upper bound and the Monte Carlo data can be appreciated in Figure 4; compare with Figure 1. Note how the data lie only very slightly below the upper bound (24).

Figure 4: Measured failure rate for Algorithm S against the upper bound (24). Error bars in the data are much smaller than the symbols shown.

References

(1) D. Aldous, J. A. Fill, Reversible Markov Chains and Random Walks on Graphs. Unfinished monograph (recompiled 2014) available at http://www.stat.berkeley.edu/~aldous/RWG/book.html.
(2) S. G. Akl, A new algorithm for generating derangements, BIT Numer. Math. 20 (1) (1989) 2–7.
(3) R. Arratia, A. D. Barbour, S. Tavaré, Logarithmic Combinatorial Structures: A Probabilistic Approach. EMS, Zürich (2003).
(4) A. Bacher, O. Bodini, H.-K. Hwang, T.-H. Tsai, Generating random permutations by coin tossing: Classical algorithms, new analysis, and modern implementation, ACM Trans. Algorithms 13 (2) (2017) 24.
(5) J. L. Baril, V. Vajnovszki, Gray code for derangements, Discrete Appl. Math. 140 (1–3) (2004) 207–221.
(6) I. Beichl, F. Sullivan, Approximating the permanent via importance sampling with application to the dimer covering problem, J. Comput. Phys. 149 (1) (1999) 128–147.
(7) R. A. Brualdi, R. J. Ryser, Combinatorial Matrix Theory, Cambridge University Press, Cambridge, 1991.
(8) C. A. Charalambides, Enumerative Combinatorics, Chapman & Hall/CRC, Boca Raton, 2002.
(9) Y. Chen, P. Diaconis, S. P. Holmes, J. S. Liu, Sequential Monte Carlo methods for statistical analysis of tables, J. Am. Stat. Assoc. 100 (469) (2005) 109–120.
(10) P. Diaconis, Group Representations in Probability and Statistics, IMS, Hayward (1988).
(11) P. Diaconis, R. L. Graham, S. P. Holmes, Statistical problems involving permutations with restricted positions, in: M. de Gunst, C. Klaassen, A. Van der Vaart (Eds.), State of the Art in Probability and Statistics: Festschrift for Willem R. van Zwet, IMS, Beachwood, 2001, pp. 195–222.
(12) P. Flajolet, M. Soria, Gaussian limiting distributions for the number of components in combinatorial structures, J. Comb. Theor. Ser. A 53 (2) (1990) 165–182.
(13) D. Gries, J. Xue, Generating a random cyclic permutation, BIT Numer. Math. 28 (3) (1988) 569–572.
(14) P. Hanlon, A random walk on the rook placements on a Ferrer’s board, Electron. J. Comb. 3 (2) (1996) 26.
(15) J. F. Korsh, P. S. LaFollette, Constant time generation of derangements, Inf. Process. Lett. 90 (4) (2004) 181–186.
(16) N. Y. Kuznetsov, Computing the permanent by importance sampling method, Cybern. Syst. Anal. 32 (6) (1996) 749–755.
(17) H. Liang, L. Shi, F. Bai, X. Liu, Random path method with pivoting for computing permanents of matrices, Appl. Math. Comput. 185 (1) (2007) 59–71.
(18) C. Martínez, A. Panholzer, H. Prodinger, Generating random derangements, in: R. Sedgewick, W. Szpankowski (Eds.), 2008 Proc. Fifth Workshop on Analytic Algorithmics and Combinatorics – ANALCO, SIAM, Philadelphia, 2008, pp. 234–240.
(19) E. Ozel, The number of $k$ -cycles in a family of restricted permutations, arXiv:1710.07885 (2017).
(20) A. Panholzer, H. Prodinger, M. Riedel, Measuring post-quickselect disorder, J. Iran. Stat. Soc. 3 (2) (2004) 219–249.
(21) H. Prodinger, On the analysis of an algorithm to generate a random cyclic permutation, Ars Comb. 65 (2002) 75–78.
(22) L. E. Rasmussen, Approximating the permanent: A simple approach, Random Struct. Algor. 5 (2) (1994) 349–361.
(23) S. Sattolo, An algorithm to generate a random cyclic permutation, Inf. Process. Lett. 22 (6) (1986) 315–317.
(24) R. Sedgewick, Permutation generation methods, Comput. Surv. 9 (2) (1977) 137–164.
(25) A. Smith, Comparison theory for Markov chains on different state spaces and application to random walk on derangements, J. Theor. Probab. 28 (4) (2015) 1406–1430.
(26) S. Vigna, xoshiro/xoroshiro generators and the PRNG shootout. Available at http://xoshiro.di.unimi.it/.
(27) M. C. Wilson, Random and exhaustive generation of permutations and cycles, Ann. Comb. 12 (4) (2009) 509–520.
(28) Wolfram Research, Inc., Mathematica, Version 11.3, Champaign, IL (2018).

Efficient uniform generation of random derangements with the expected distribution of cycle lengths

On the uniform generation of random derangements

Symmetrized importance samplers for stochastic differential equations

Stein Variational Adaptive Importance Sampling

An Importance Sampling Algorithm Based on Evidence Pre-propagation

Random Sampling using k-vector

Randomized sequential importance sampling for estimating the number of perfect matchings in bipartite graphs

Random derangements and the Ewens Sampling Formula

1 Introduction

2 Mathematical preliminaries

Proposition 2.1.

Proof.

Proposition 2.2.

Remark 2.3.

3 Generating random derangements by random transpositions

Remark 3.1.

Remark 3.2.

Remark 3.3.

4 Sequential importance sampling of derangements

4.1 The SIS algorithm

4.2 Failure probability of the SIS algorithm

Theorem 4.1.

Proof.

Remark 4.2.

4.3 Uniformity of the SIS algorithm

Definition 4.3.

Lemma 4.4.

Proof.

Lemma 4.5.

Proof.

Theorem 4.6.

Proof.

5 Mixing time of the restricted transpositions shuffle

Remark 5.1.

6 Summary and conclusions

Acknowledgments

Appendix A Improved upper bound for the failure probability of Algorithm S

References

Efficient uniform generation of random derangements with the expected distribution of cycle lengths

Related Research

On the uniform generation of random derangements

Symmetrized importance samplers for stochastic differential equations

Stein Variational Adaptive Importance Sampling

An Importance Sampling Algorithm Based on Evidence Pre-propagation

Random Sampling using k-vector

Randomized sequential importance sampling for estimating the number of perfect matchings in bipartite graphs

Random derangements and the Ewens Sampling Formula

1 Introduction

2 Mathematical preliminaries

Proposition 2.1.

Proof.

Proposition 2.2.

Remark 2.3.

3 Generating random derangements by random transpositions

Remark 3.1.

Remark 3.2.

Remark 3.3.

4 Sequential importance sampling of derangements

4.1 The SIS algorithm

4.2 Failure probability of the SIS algorithm

Theorem 4.1.

Proof.

Remark 4.2.

4.3 Uniformity of the SIS algorithm

Definition 4.3.

Lemma 4.4.

Proof.

Lemma 4.5.

Proof.

Theorem 4.6.

Proof.

5 Mixing time of the restricted transpositions shuffle

Remark 5.1.

6 Summary and conclusions

Acknowledgments

Appendix A Improved upper bound for the failure probability of Algorithm S

References