Random permutations without double rises (avoiding consecutive pattern $\underline{123}$)

Question

A permutation avoiding a consecutive pattern $\underline{123}$ is permutation $\pi = \pi_1 \pi_2 \ldots \pi_n$ with the property that there does not exists $i \in [1, n-2]$ such that $\pi_i < \pi_{i+1} < \pi_{i+2}$. Example: $53241$ is $\underline{123}$-avoiding; while $314562$ is not $\underline{123}$ avoiding, as it contains $145$.

Let $a_n$ be the number of $\underline{123}$ avoiding permutations, it is known, in corresponding OEIS entry A049774, that $$ \frac{a_n}{n!} \sim e^{\frac{\pi}{3\sqrt{3}}} \left(\frac{3\sqrt{3}}{{2\pi}}\right)^{n+1}.$$

So, the simple rejection sampling is inefficient even for moderate values of $n$. How does one generate such permutations uniformly at random?

I'm aware of Boltzmann sampling. But maybe you are aware of more simple and/or faster algorithms as in the case of alternating permutations.

Answer 1

I have a sampling algorithm for you written in sage / python.

Instead of $\underline{123}$ avoiding permutations of length $n$ I consider words $w_1 w_2 \dots w_n$ satisfying the following properties:

• $0\le w_k \le n-k$
• the word is weakly $\underline{123}$ avoiding, which means there is no index $i\in [1,n-2]$ such that $w_i\le w_{i+1} \le w_{i+2}$.

I call these words $\underline{123}$ avoiding indexwords. They are in bijection with $\underline{123}$ avoiding permutations using the following bijection: Let $S=\{1,2,\dots,n\}$ then for a given indexword $w_1 w_2 \dots w_n$ we define $\pi_i = w_i\text{-th smallest entry of } S\setminus\{\pi_1,\dots,\pi_{i-1}\}$, where I start indexing the entries of the set with $0$. Put differently: If we keep track of the unused letters to build up the permutation in an ordered list, then $\pi_i$ is precisely the $w_i\text{-th}$ entry in that list (explaining my name for these words).

These words will simplify the following algorithm for random uniform sampling. The algorithm comes in two steps: First count the number of $\underline{123}$ avoiding indexwords with a given prefix. Then create a random word iteratively by appending letters to the already constructed word using certain probabilities using the precomputed values.

For this let $\alpha_n^{w_1 w_2 \dots w_l}$ be the number of $\underline{123}$ avoiding indexwords of length $n$ starting with $w_1 w_2 \dots w_l$. For calculating these values, we have two rules: Let $w_1 w_2 \dots w_l$ be weakly $\underline{123}$ avoiding, then

1. Expansion rule $$ \alpha_n^{w_1 w_2 \dots w_l} = \sum_{0\le x\le n-l-1 \text{ and } \neg(w_{l-1}\le w_l\le x)} \alpha_n^{w_1 w_2 \dots w_l x} $$
2. Reduction rule $$ \alpha_n^{w_1 w_2 \dots w_l} = \alpha_{n-l+2}^{w_{l-1}w_l} $$
3. Start cases $$ \alpha_2^{00}=\alpha_2^{10}=1 $$

Combining the rules we obtain $$ \alpha_n^{xy}=\sum_{0\le z\le n-3 \text{ and } \neg(x\le y\le z)}\alpha_{n-1}^{yz} $$

Note that you can obtain $a_n$ by summing $\alpha_n^{xy}$ over all pairs $(x,y)\in [0,n-1]\times [0,n-2]$

Now for the sampling part. Randomly choose a pair $(x,y)\in [0,n-1]\times [0,n-2]$ with probability $\frac{\alpha_n^{xy}}{a_n}$. Assume now, we have already created a sequence $w_1 \dots w_l$ with $l<n$, we append an appendable letter $x\in [0,n-l-1]$ with probabilty $\frac{\alpha_{n+1-l}^{w_l x}}{\alpha_{n+2-l}^{w_{l-1} w_{l}}}$. (A letter is appendable, if it does not create the $\underline{123}$ pattern). Repeat this until you have a $\underline{123}$ avoiding indexword of length $n$.

Here is the sage code:

@cached_function
def alpha(n,x,y):
    assert n>1, "n needs to be at least 2"
    if n == 2:
        return 1
    return sum(alpha(n-1,y,z) for z in range(n-2) if not (x<=y<=z))

def random_indexword(n):
    assert n>=2, "Algorithm only works for length n>=2"
    pairs = [(x,y) for x in range(n) for y in range(n-1)]
    P = [alpha(n,x,y) for x,y in pairs]
    s1,s2 = pairs[GeneralDiscreteDistribution(P).get_random_element()]
    indexword = [s1,s2]
    while len(indexword)<n:
        appendable = [x for x in range(n-len(indexword)) if not (indexword[-2] <= indexword[-1] <= x)]
        P = [alpha(n+1-len(indexword), indexword[-1],x) for x in appendable]
        indexword.append(appendable[GeneralDiscreteDistribution(P).get_random_element()])

    return indexword

def indexword_to_permutation(word):
    letters = [i for i in range(1,len(word)+1)]
    p = []
    for i in word:
        p.append(letters[i])
        letters.remove(p[-1])
    return p

# for testing
def is_123_avoiding(perm):
    return not any(perm[i]<perm[i+1]<perm[i+2] for i in range(len(perm)-2))

w = random_indexword(10)
perm = indexword_to_permutation(w)
perm,is_123_avoiding(perm)

Edit: Updated typo in bijection between index words and permutations.

Answer 2

In Igor Pak's answer about alternating permutations that is linked to, the key is to precompute a triangle of numbers, and then you can use the triangle to inform the probabilities of your selections.

One way to break apart your search space (that is, the $\underline{123}$-avoiding permutations) would be via Doron Zeilberger's "enumeration schemes". Andrew Baxter and Lara Pudwell have looked at enumeration schemes for permutations avoiding patterns such as $\underline{123}$ (which can be viewed as a dashed pattern or a vincular pattern). Here are slides from a talk Pudwell gave about their work. Here is the arXiv preprint. Here is the published version (subscription required).

They don't mention the pattern $\underline{123}$ explicitly from what I can see (it would be denoted simply $123$ in their notation), but they do mention the pattern $\underline{1234}$, so I have to imagine there is a finite enumeration scheme for the pattern $\underline{123}$.

Assuming you went down this route, you would precompute the enumeration scheme for $\underline{123}$-avoiding permutations, precompute all the numbers involved up to the length you are interested in, and then use those numbers to inform random selections (as in the alternating permutations example). Assuming I'm correct that the $\underline{123}$-avoiding permutations have a finite enumeration scheme, this would all be polynomial time in the length of the permutation you want to generate.