Randomize List with guaranteed Index Change

Question

I usually post on Stackoverflow, but I got stuck with a problem which I think fits better here and I hope someone can help me. The task ist pretty simple:

Shuffle/Randomize an array/list

conditions:

1. each Index has to change; an Item which was at index 0 before the shuffle has to be at an index != 0 after the shuffle
2. all possible resulting permutations which satisfy the first condition should be distributed equally

I am unable to find a satisfying solution. Brute force (shuffling the array with Fisher-Yates until first condition is met) obviously works but is not a solution I am looking for. Brute Force with an array with four elements resulted in following number of permutations (which is roughly distributed equally) for one million tries: [111467, 111172, 111497, 111047, 111114, 111053, 110699, 110655, 111296]

Another thing I tried was something like this:

1. iterate through the array object by object and randomly choose a new index for the object
2. the randomly chosen index must not have been chosen before and must not be same as the original index
3. if the last element would fall on its original index, swap the last element with a random other element

Using this approach, the distribution looks really bad: [74166, 110734, 83581, 129486, 83210, 101904, 111305, 111235, 194379].

In theory, there must be a better solution than brute force but unfortunately everybody I asked so far (including google) didn't have a solution, or just had another non uniform solution comparable to my second approach.

Can anyone help me please? I would appreciate it a lot.

Answer 1

Algorithm

There is an idea based on subfactorial formula !n = (n-1)(!(n-1) + !(n-2)).

There are two ways to generate derangement of size n:

1. Generate derangement of size n-1, append n-th item, swap n-th and random item in range [1, n-1].

2. Generate derangement of size n-2, get i = random[1, n-1], increment all values >= i in the derangement, insert value i to i-th place. The modified permutation is not derangement any more, but the only item (i) is on his place. Append n-th item, swap n-th and i-th items. The result is derangement.

To generate derangement we select first of second way with probabilities !(n-1) / (!(n-1) + !(n-2)) and !(n-2) / (!(n-1) + !(n-2)).

The only issue is that this algorithm is not linear of n in time. I think it is O(n^2) due to item insertion and increment in the second way.

Evaluation

Below is test in Python. It runs permutation generator 1000000 times and shows number of calls to random generator, expected number of permutations of particular kind, distribution of permutations.

MotiNK has been added in the last moment to check it. Unfortunately it does not generate all permutations.

Derangement implements the idea described above. It works OK but it is quadratic. :(

FisherYatesDerangement1 is Fisher-Yates "try and reject".

FisherYatesDerangement2 is optimized slightly. It is possible to reject permutation early while Fisher-Yates is generating permutation.

FisherYates is added as a baseline. It does not generate derangements but all permutations.

import collections
import random


def is_derangement(a):
    for i, v in enumerate(a):
        if i == v:
            return False
    return True


# suggested by [MotiNK](https://cs.stackexchange.com/users/37884/motink)
#
# for i=0..n-2
#   swapIdx = random(i+1,n) // random(start,end) --> start <= i < end
#   swapElement(array, i, swapIdx)


class MotiNK:
    def __init__(self, n):
        self._n = n
        self._c = 0

    def next(self):
        a = list(range(self._n))
        for i in range(self._n - 1):
            j = random.randrange(i + 1, self._n)
            self._c += 1
            a[i], a[j] = a[j], a[i]
        assert is_derangement(a)
        return tuple(a)

    def count(self):
        return self._c


class FisherYates:
    def __init__(self, n):
        self._n = n
        self._c = 0

    def next(self):
        a = list(range(self._n))
        for i in range(self._n):
            j = random.randrange(i, self._n)
            self._c += 1
            a[i], a[j] = a[j], a[i]
        return tuple(a)

    def count(self):
        return self._c


class FisherYatesDerangement1:
    def __init__(self, n):
        self._fy = FisherYates(n)

    def next(self):
        while True:
            a = self._fy.next()
            if is_derangement(a):
                return a

    def count(self):
        return self._fy.count()


class FisherYatesDerangement2:
    def __init__(self, n):
        self._n = n
        self._c = 0

    def next(self):
        while True:
            a = list(range(self._n))
            if self._shuffle(a):
                assert is_derangement(a)
                return tuple(a)

    def count(self):
        return self._c

    def _shuffle(self, a):
        for i in range(self._n):
            j = random.randrange(i, self._n)
            self._c += 1
            if i == a[j]:
                return False
            a[i], a[j] = a[j], a[i]
        return True


def subfactorial(n):
    f1 = 0
    f2 = 1
    for i in range(n):
        f1, f2 = f2, i * (f1 + f2)
    return f2


class Derangement:
    def __init__(self, n):
        self._n = n
        self._c = 0

    def next(self):
        a = self._derangement(self._n)
        assert(is_derangement(a))
        return tuple(a)

    def count(self):
        return self._c

    def _derangement(self, n):
        assert n != 1
        if n == 0:
            return []
        n1 = subfactorial(n - 1)
        n2 = subfactorial(n - 2)
        k = random.randrange(n1 + n2)
        self._c += 1
        i = random.randrange(n - 1)
        self._c += 1
        if k < n1:
            a = self._derangement(n - 1)
        else:
            a = self._derangement(n - 2)
            for j in range(n - 2):
                if a[j] >= i:
                    a[j] += 1
            a.insert(i, i)
        a.append(n - 1)
        a[i], a[n - 1] = a[n - 1], a[i]
        return a


def test_generator(k, label, gen):
    c = collections.Counter(gen.next() for _ in range(k))
    m = len(c)
    print(
        label,
        f'{gen.count():10}',
        f'{k / m:10.2f}',
        sorted(c.values())
    )


for n in range(3, 6):
    print('n =', n)
    m = 1000000
    test_generator(m, 'MotiNK                 ', MotiNK(n))
    test_generator(m, 'Derangement            ', Derangement(n))
    test_generator(m, 'FisherYatesDerangement1', FisherYatesDerangement1(n))
    test_generator(m, 'FisherYatesDerangement2', FisherYatesDerangement2(n))
    test_generator(m, 'FisherYates            ', FisherYates(n))

Results:

n = 3
MotiNK                     2000000  500000.00 [499557, 500443]
Derangement                4000000  500000.00 [499959, 500041]
FisherYatesDerangement1    8986053  500000.00 [499855, 500145]
FisherYatesDerangement2    6493519  500000.00 [498982, 501018]
FisherYates                3000000  166666.67 [166124, 166291, 166399, 166842, 166971, 167373]
n = 4
MotiNK                     3000000  166666.67 [165922, 166499, 166550, 166680, 166937, 167412]
Derangement                5332446  111111.11 [110578, 110775, 111087, 111098, 111105, 111165, 111184, 111501, 111507]
FisherYatesDerangement1   10663028  111111.11 [110525, 110740, 110906, 111113, 111193, 111203, 111421, 111436, 111463]
FisherYatesDerangement2    7434896  111111.11 [110378, 110400, 111099, 111141, 111186, 111194, 111300, 111543, 111759]
FisherYates                4000000   41666.67 [41325, ..., 41911]
n = 5
MotiNK                     4000000   41666.67 [41338, ..., 42069]
Derangement                7091636   22727.27 [22449, ..., 23012]
FisherYatesDerangement1   13626165   22727.27 [22403, ..., 23004]
FisherYatesDerangement2    9344157   22727.27 [22392, ..., 23014]
FisherYates                5000000    8333.33 [8118, ..., 8632]

Answer 2

NOTE: updated (before it had low=i+1 always, this missed some valid results, and missed the extra check needed for the last index) **NOTE: update 2 - randomize order to ensure uniformity **

Consider the following approach:

loopIndex[0..n-1] = permute([0..n-1])
index[0..n-1] = [0..n-1] // track indices
for i=0..n-2
   idx = loopIndex[i]
   if index[idx] == idx
      low = i+1
   else low = i
   swapIdx = random(low,n) // random(start,end) --> start <= i < end
   swapElement(index, idx, loopIndex[swapIdx])
   swapElement(array, idx, loopIndex[swapIdx])
finalIdx = loopIndex[n-1]
if index[finalIdx] == finalIdx
   finalSwapIdx = random(0,n-1)
   swapElement(index, finalSwapIdx, finalIdx)
   swapElement(array, finalSwapIdx, finalIdx)

PROOF NOT YET UPDATED Note the following:

1. loopIndex is used to give a random ordering on which we process the elements. This ensures that any non-uniformity is "smoothed out".
2. If an index $j$ is selected as swapIdx then that element is swapped to a place where it's not processed again in the loop.
3. If an index $j < n-1$ is never selected as swapIdx then when $i==j$ the element will be swapped out from $j$ and never can return to $j$.
4. The final check handles the case where the final element was never selected. In this case, we swap between $n-1$ and some other random element.

Proof of uniformity by induction on the index - Let $p(j,k)$ be the probability that after completing $i=j$, the element at index $j$ is the element that originally had index $k \neq j$ (clearly $p(j,j)=0$). Then $p(0,k)=\frac{1}{n-1}$ clearly.

Assume then that the statement holds for all $j \leq t$ and we now prove for $j=t+1$.

$p(t+1,k)$ is the probability the element with original index $k$ is not in the first $t$ elements at step $i=t+1$, times the probability that $k$ is swapIdx now. By the assumption, this is then: $$ \left(\frac{n-2}{n-1}\right)\left(\frac{n-3}{n-2}\right)\dots\left(\frac{n-t-1}{n-t}\right)\left(\frac{1}{n-t-1}\right) = \frac{1}{n-1} $$ (Note that there is a minor difference between $k \leq t$ and $k > t+1$. For the former, it must not have been selected as a swapIdx for $j \neq k \leq t$, and it must not have selected $t+1$ as its swapIdx. For the latter, it must not have been selected as a swapIdx. The result though is the same in terms of the math).