Combinatorics : Derangements

This whole tutorial starts with a simple example problem

Suppose we have 99 empty papers and we wrote numbers from 1 to 99 (using each number) on one side of the papers randomly. Then we mixed all the papers randomly and started to write numbers from 1 to 99 on the other sides (empty sides) of each paper randomly. What is the probability of having the same number on both sides of at least one paper?

https://math.stackexchange.com/questions/394007/a-tricky-probability-question

Give it a go and quickly you will see how tricky it is. The answer requires an understanding of derangements and we will go through this idea in detail in this post.

What is a derangement?

So what is a derangement? Consider a set of integers

1 2 3 4 5

A permutation is any ordering of these elements. For a list of $n$ integers; there are $n!$ permutations. A permutation may also be called an arrangement. A derangement is a permutation/arrangement in which none of the integers appear in their correct position.

2 3 5 1 4  # Example derangement

From this description you may be able to see that the problem posed at the start is closely related to the number of derangements of the sequence of all integers up to 100.

We will take a concrete example that you can reference throughout (see the appendix for the code used to produce this) involving a set with 5 elements. Each column refers to the number of elements in their correct position.

   0      1      2      3      4      5      
===========================================
|BADEC |ACBED |ABDEC |ABCED |      |ABCDE |
|BAECD |ACDEB |ABECD |ABDCE |      |      |
|BCAED |ACEBD |ACDBE |ABEDC |      |      |
|BCDEA |ADBEC |ACEDB |ACBDE |      |      |
|BCEAD |ADEBC |ADBCE |ADCBE |      |      |
|BDAEC |ADECB |ADCEB |AECDB |      |      |
|BDEAC |AEBCD |AEBDC |BACDE |      |      |
|BDECA |AEDBC |AECBD |CBADE |      |      |
|BEACD |AEDCB |BCADE |DBCAE |      |      |
|BEDAC |BACED |BDCAE |EBCDA |      |      |
|BEDCA |BADCE |BECDA |      |      |      |
|CABED |BAEDC |CABDE |      |      |      |
|CADEB |BCDAE |CBDAE |      |      |      |
|CAEBD |BCEDA |CBEDA |      |      |      |
|CDAEB |BDACE |DACBE |      |      |      |
|CDBEA |BDCEA |DBACE |      |      |      |
|CDEAB |BEADC |DBCEA |      |      |      |
|CDEBA |BECAD |EACDB |      |      |      |
|CEABD |CADBE |EBADC |      |      |      |
|CEBAD |CAEDB |EBCAD |      |      |      |
|CEDAB |CBAED |      |      |      |      |
|CEDBA |CBDEA |      |      |      |      |
|DABEC |CBEAD |      |      |      |      |
|DAEBC |CDABE |      |      |      |      |
|DAECB |CDBAE |      |      |      |      |
|DCAEB |CEADB |      |      |      |      |
|DCBEA |CEBDA |      |      |      |      |
|DCEAB |DABCE |      |      |      |      |
|DCEBA |DACEB |      |      |      |      |
|DEABC |DBAEC |      |      |      |      |
|DEACB |DBEAC |      |      |      |      |
|DEBAC |DBECA |      |      |      |      |
|DEBCA |DCABE |      |      |      |      |
|EABCD |DCBAE |      |      |      |      |
|EADBC |DECAB |      |      |      |      |
|EADCB |DECBA |      |      |      |      |
|ECABD |EABDC |      |      |      |      |
|ECBAD |EACBD |      |      |      |      |
|ECDAB |EBACD |      |      |      |      |
|ECDBA |EBDAC |      |      |      |      |
|EDABC |EBDCA |      |      |      |      |
|EDACB |ECADB |      |      |      |      |
|EDBAC |ECBDA |      |      |      |      |
|EDBCA |EDCAB |      |      |      |      |
|      |EDCBA |      |      |      |      |
===========================================
   44     45     20     10     0     1     

Getting it Wrong (1)

Getting it Right (1)

Getting it Right (2)

Proof using PIE

Example Problems using Derangements

Derangements using Python

This was a maths-heavy post so far, which is good for Olympiad style maths problems; but there are also potentially occasions where derangements are needed in code. If nothing else; having the code allows you to check solutions to problems.

So, let’s write some code to calculate the total number of derangements.

The first option is to just “brute force” it; let’s enumerate every possible permutation and check for each one whether the permutation is a derangement or not.

Below is the code

import numpy as np
import math
from itertools import permutations

def brute_force_derangement(n):
    
    a = np.arange(1,n,1)

    deranged = 0
    for perm in permutations(a):
        if not np.any(perm == a):
            deranged += 1
            
    return deranged / math.factorial(n)

This certainly does the job and is very explicit; let’s see how long it takes for n=10 (on my machine, MacBook with a core i5 processor)

%timeit brute_force_derangement(10)

4.55 s ± 244 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Ok, a few seconds is reasonable for a quick “back of the envelope” code check. But this won’t cut it for a real application; and it isn’t scalable at all.

Lets dig into the scalability a bit more…

First of all, there is a factorial growth of executions the for loop. That’s going to kill the algorithm for n a little above 10. Even the algorithm in the loop itself is inefficient; it requires an element by element check of a list, which grows as O(n).

We could use our analytical solution derived in Eq.(1). In principle this is a good idea; but you will find that it requires the code to deal with some large numbers (via the factorial) which may slow things down. Let’s code it and check

def derangement_formula(n):
    
    tot= 0
    for i in range(n):
        tot += (-1)**i / math.factorial(i)
    return int(math.factorial(n)*tot)

%timeit derangement_formula(10)

6.4 µs ± 261 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Ok, pretty good! a million times faster. That is a victory for the maths; analytical derivation of a formula has prevented our code requiring a exhaustive brute force solution. Now let’s grab a victory for the code.

So how do we improve this? Although the recursion formula we derived earlier was of limited use in practical problems (our brains don’t really do recursion very well) aside from its use in derivation. Let’s try the first recursive formula Eq.(4)

def derangement_recursive_1(n):
    current_n   = 1
    current_d_n = 1
    while current_n < n:
        current_d_n = current_n * current_d_n + (-1)**current_n
        current_n += 1
        
    return current_n

%timeit derangement_recursive_1(10)

5.01 µs ± 606 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Ok this has sped things up a bit; but not substantially. What is holding us up? Well, generally the power operation has a high computational cost. Multiplication and addition are lower on computation cost.

So let's try to remove the power; in this case it's very easy since all that's happening is that each term alternates between positive and negative. This can be done with multiplication, but it requires one extra variable to be held in memory. In this case it will benefit us...

def derangement_recursive_2(n):
    current_n   = 1
    current_d_n = 1
    prev = 1
    while current_n < n:
        prev *= -1
        current_d_n = current_n * current_d_n + prev
        current_n += 1
        
    return current_n

%timeit derangement_recursive_2(10)

1.76 µs ± 86.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

More than twice as fast!

Can we do even better? Within each iteration of the loop there are 2 multiplies and 2 adds. We can reduce this to one multiply and three add by using the other recursion relation Eq.(3)

def derangement_recursive_3(n):
    current_n             = 1
    current_d_n           = 1
    current_d_n_minus_one = 1 
    while current_n < n:
        tmp = current_d_n
        current_d_n = (current_n-1) * (current_d_n + current_d_n_minus_one)
        current_d_n_minus_one = tmp
        current_n += 1
        
    return current_n

%timeit derangement_recursive_3(10)

1.71 µs ± 42.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

This doesn't really result in much gain. However we can re-organise to reduce the number of adds operations to two

def derangement_recursive_4(n):
    current_n             = 1
    current_d_n           = 1
    current_d_n_minus_one = 1 
    while current_n < n:
        tmp = current_d_n
        current_n += 1
        current_d_n = current_n * (current_d_n + current_d_n_minus_one)
        current_d_n_minus_one = tmp
        
        
    return current_n

%timeit derangement_recursive_4(10)

1.63 µs ± 95.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Now the gains we are seeing are very modest and we should probably stop.

To tie it all off, I have run these functions through 10 repeats of the timeit function for n between 2 and 20

Appendix : Code

Here is the code used to produce the ASCII style example of all permutations of set at the start of this post

import numpy as np

def pretty_print_permutations(a):
    n = len(a)
    
    count    = 0 
    deranged = []
    corr = []

    for d in permutations(a):
        corr     += [sum(list(d) == a)]
        deranged += ["".join(d)]

    picks = dict()
    for i in range(n+1):
        examples = []
        for index in np.argwhere(np.array(corr) == i):
            examples += [deranged[index[0]]]
        picks[i] = examples

    columns = n+1

    text = " "*(int(n/2)+1) + "".join(["{}".format(i) + " "*(n+1) for i in range(n+1)]) + "\n"
    linelen = len(text)
    text += "=" * (linelen-3) + "\n"

    max_i = max([len(x) for x in picks.values()])
    for i in range(max_i):
        line = ""
        for j in range(columns):
            try:
                line += "|" + str(picks[j][i]) + " "
            except:
                line += "|" + " "*(n+1)   
        line += "|\n"
        text += line

    text += "=" * (linelen-3) + "\n"
    text += " "*(int(n/2)+1) + "".join(["{}".format(len(i)) + " "*n for i in picks.values()]) + "\n"
    
    return text

And here is the code used to create the speed test plot;

import timeit
import numpy as np
import math
from itertools import permutations
import matplotlib.pyplot as plt

"""
A speed test of various ways to calculate the derangement probability 
"""


def brute_force_derangement(n):
    """
    Find the probability of a derangement being found
    from a random permutation of n elements by a brute
    force search
    """
    a = np.arange(1, n + 1, 1)

    deranged = 0
    for perm in permutations(a):
        if not np.any(perm == a):
            deranged += 1

    return deranged / math.factorial(n)


def derangement_formula(n):
    """
    Find the probability of a derangement being found
    from a random permutation of n elements using the
    standard analytic form
    """
    tot = 0
    for i in range(n):
        tot += (-1) ** i / math.factorial(i)
    return int(math.factorial(n) * tot)


def derangement_recursive_1(n):
    """
    Find the probability of a derangement being found
    from a random permutation of n elements using a
    recursive formulae
    """
    current_n = 1
    current_d_n = 1
    while current_n < n:
        current_d_n = current_n * current_d_n + (-1) ** current_n
        current_n += 1

    return current_n


def derangement_recursive_2(n):
    """
    Find the probability of a derangement being found
    from a random permutation of n elements using a
    recursive formulae
    """
    current_n = 1
    current_d_n = 1
    prev = 1
    while current_n < n:
        prev *= -1
        current_d_n = current_n * current_d_n + prev
        current_n += 1

    return current_n


def derangement_recursive_3(n):
    """
    Find the probability of a derangement being found
    from a random permutation of n elements using a
    recursive formulae
    """
    current_n = 1
    current_d_n = 1
    current_d_n_minus_one = 1
    while current_n < n:
        tmp = current_d_n
        current_d_n = (current_n - 1) * (current_d_n + current_d_n_minus_one)
        current_d_n_minus_one = tmp
        current_n += 1

    return current_n


def derangement_recursive_4(n):
    """
    Find the probability of a derangement being found
    from a random permutation of n elements using a
    recursive formulae
    """
    current_n = 1
    current_d_n = 1
    current_d_n_minus_one = 1
    while current_n < n:
        tmp = current_d_n
        current_n += 1
        current_d_n = current_n * (current_d_n + current_d_n_minus_one)
        current_d_n_minus_one = tmp

    return current_n


if __name__ == "__main__":

    def get_times(func, n_values, repeats=1):
        res = np.zeros(len(n_values))
        for r in range(repeats):
            res += [timeit.timeit(lambda: func(n)) for n in n_values]
        return np.vstack((n_values, res / repeats))


    n_values = np.arange(2, 21, 1)

    d1 = get_times(derangement_formula, n_values, repeats=10)
    d2 = get_times(derangement_recursive_1, n_values, repeats=10)
    d3 = get_times(derangement_recursive_2, n_values, repeats=10)
    d4 = get_times(derangement_recursive_3, n_values, repeats=10)
    d5 = get_times(derangement_recursive_4, n_values, repeats=10)

    np.savetxt("speed_test.csv", np.hstack((n_values, d1[0], d2[0], d3[0], d4[0], d5[0])), delimiter=",")

    fig, ax = plt.subplots(1, 1, figsize=(8, 8))

    ax.plot(d1[0], d1[1], label="Standard Formula", color="black")
    ax.plot(d2[0], d2[1], label="Recursive 1", color="b")
    ax.plot(d3[0], d3[1], label="Recursive 2", color="r")
    ax.plot(d4[0], d4[1], label="Recursive 3", color="g")
    ax.plot(d5[0], d5[1], label="Recursive 4", color="purple")

    ax.scatter(d1[0], d1[1], color="black")
    ax.scatter(d2[0], d2[1], color="b")
    ax.scatter(d3[0], d3[1], color="r")
    ax.scatter(d4[0], d4[1], color="g")
    ax.scatter(d5[0], d5[1], color="purple")

    ax.set_title(r"$\rm Timings \ for \ derangement \ implementations $", fontsize=18)
    ax.set_xlabel(r"$n$", fontsize=18)
    ax.set_ylabel(r"$\rm Time \ (\mu s)$", fontsize=18)

    ax.set_xlim([1, 20])
    ax.set_ylim([0, 25])

    ax.legend(frameon=False, fontsize=16)
    ax.minorticks_on()
    ax.set_xticks(n_values)

    fig.savefig("derangement_speed_test.pdf")

    plt.show()