Random Alternating Permutations

Question

An alternating permutation of {1, ..., n} is one were π(1) > π(2) < π(3) > π(4) < ... For example: (24153) is an alternating permutation of length 5.

If $E_n$ is the number of alternating permutations of length n, the $\sec x + \tan x = \sum_{n \geq 0} E_n \frac{x^n}{n!}$ is the exponential generating function.

How does one sample the alternating permutations uniformly at random? Rejection sampling would quickly become inefficient since $\frac{E_n}{n!} \approx \frac{4}{\pi} \left(\frac{2}{\pi} \right)^n$ decays exponentially in n.

Answer 1

This is very easy - I often teach this in my combinatorics classes. You start with this Pascal type triangle which you need to precompute: http://mathworld.wolfram.com/Seidel-Entringer-ArnoldTriangle.html (read Arnold's paper referenced there if the pattern is unclear or my survey paper with Postnikov "Increasing Trees and Alternating Permutations" - see my webpage).

Now note that $a_{n,i}$, the i-th number in n-th row is the number of alternating permutations with the first letter $\ge i$ (you need to alternate directions). Then you use dynamic programming or simply choose the first letter to be i with probability $a_{n,i} - a_{n,i+1}$, which reduces the problem to a alternating permutation of $[n-1]$. I am skipping (easy) details, but I trust you can figure out the rest from here.

UPDATE (May 10, 2015):
Apparently, this algorithm requires $O^\ast(n^3)$ time to compute the numbers $a_{n,i}$ which is not bad. However, if one only wants to compute the number of alternating permutations, there is a faster way: see this talk by Richard Brent and the slides here.

Answer 2

An $O(n \log n)$ algorithm to uniformly sample alternating permutation is described by P. Marchal in *Generating alternating permutations in time $O(n \log n)$." (by contrast, Pak's answer is $O(n^3),$ Boltzmann sampling is $O(n^2).$

Answer 3

This is covered in some detail in the paper "Boltzmann sampling of ordered structures".

Answer 4

Below is Mathematica code based on Igor Pak's answer. To get a random downup permutation on $[n]$, start by choosing the first entry $p_1$ with the appropriate probability; then randomly choose an updown permutation of size $n-1$ from the updown permutations with first entry $< p_1$; then join them together (incrementing entries $\ge p_1$ in the updown permutation).

To implement this method, we actually need code to generate a random downup permutation with first entry $\ge$ a specified number $k$ and the code below does so. It uses the ComplementPermutation operation to interchange updown and downup permutations.

(* e[n,k] is the Entringer number *)

e[0,0] = 1;
e[n_,0]/;n>=1 := 0;
e[n_,k_]/;k>n || k<0 := 0
e[n_,k_] := e[n,k] = e[n,k-1] + e[n-1,n-k]

ComplementPermutation[perm_] := Module[{n=Length[perm]}, n+1-perm];

incrementSpecifiedAndUp[perm_,k_]:=perm/.{i_/;i>=k :> i+1};

partialSums[list_] := Drop[FoldList[Plus,0,list],1];

RandomUpDownPermFirstEntryAtMostk[n_,k_]/;k==n :=
RandomUpDownPermFirstEntryAtMostk[n,n-1];
RandomUpDownPermFirstEntryAtMostk[n_,k_]/;1<=k<n :=
ComplementPermutation[RandomDownUpPermFirstEntryAtLeastk[n,n+1-k]]

RandomDownUpPermFirstEntryAtLeastk[1,1]={1}; RandomDownUpPermFirstEntryAtLeastk[2,2]={2,1};

RandomDownUpPermFirstEntryAtLeastk[n_,k_]/; n>=3 && 2<=k<=n := Module[{keys,m,i,firstEntry,restOfPerm},

(* pick first entry using the Entringer distribution *)
keys=partialSums[Table[e[n-1,j],{j,k-1,n-1}]];
m=Random[Integer,{1,Last[keys]}];
i=1;
While[Not[ m<=keys[[i]] ],i=i+1];
firstEntry=k-1+i;

(* choose restOfPerm uniformly from updowns with their first entry < firstEntry *)
restOfPerm=RandomUpDownPermFirstEntryAtMostk[n-1,k-2+i];

(* amalgamate firstEntry and restOfPerm *)
Join[{firstEntry},incrementSpecifiedAndUp[restOfPerm,firstEntry]] ]

RandomDownUpPerm[1]={1};
RandomDownUpPerm[n_]/;n>=2 := RandomDownUpPermFirstEntryAtLeastk[n,2]

Sample output:
In[264]:=RandomDownUpPerm[15]
Out[264]=
{8, 2, 4, 1, 15, 6, 7, 3, 10, 9, 13, 11, 14, 5, 12}

Answer 5

One observation is that (as mentioned in Stanley's Survey ), the Alternating Permutations satisfy a recurrence $$E_{n+1}=\sum_{\textrm{odd } j \leq n}^n \binom{n}{j} E_j E_{n-j}$$ Here $j$ on the right hand side corresponds to the number of elements appearing before $1$ in the permutation.

This is probably not be the most efficient way, but it seems like you could compute $E_k$ for all $k$ up to $n$ and then build up your permutation recursively by successively choosing where the smallest element goes and dividing the remaining variables into two alternating permutations of the appropriate size based on that.

Answer 6

Edit: in previous versions of this post I claimed that the algorithm below generated random alternating permutations with uniform probability. That's not true, I apologize for not checking before posting. I don't know yet a simple way of generating them with uniform probability.

Start with an empty list L. Places in the list are numbered 1 thru n, and inserting a new element at the j-th place means the elements at places j thru the end of the list get pushed one place to the right. Then at the end of the following algorithm L will contain a random alternating permutation in L, but not with uniform probability. Input is n, the length of the permutation.

set j <- 0
set k <- 1
loop
  set j <- random integer in the interval [j+1, k]
  insert k into the list L at the j-th place
  if k = n exit loop
  k <- k+1
  set j <- random integer in the  interval [1, j]
  insert k into the list L at the j-th place
  if k = n exit loop
  k <- k+1
end loop

For the list you can either use a linked list, or a balanced binary search tree with each node containing the size of its subtree. The first is easier to code but runs in O(n^2), the second more complicated to code, but runs in O(n log n).

Edit: Here's a more concrete implementation in C++, using a linked list (I'm using the fancy random number generator, you may have to use gcc 4.2 or later)

#include <list>
#include <tr1/random>

using namespace std;
using namespace std::tr1;

list<int> random_alt_perm(int n) {
    variate_generator<mt19937, uniform_real<> > urng(mt19937(time(NULL)), uniform_real<>());
    list<int> l;
    for (int k=1, j=0; k<=n; ++k) {
        double rnd=urng();
        if (k&1) j = j + 1 + (k-j)*rnd;
        else j = 1 + j*rnd;
        list<int>::iterator it=l.begin();
        advance(it, j-1);
        l.insert(it, k);
    }
    return l;
}