ENH: Function to generate all permutations of an array

[timhoffm]

Proposed new feature or change:

I would like to be able to generate all permutations of an array; e.g.

>>> np.permutations([0, 1, 2])
array([[0, 1, 2],
       [0, 2, 1],
       [1, 0, 2],
       [1, 2, 0],
       [2, 0, 1],
       [2, 1, 0]])

AFAICS, there's neither a dedicated function for that in numpy already nor is it easily expressable through a combination of existing functions. I've only found the choice and permutation functions from the random generators, but these only create samples not a complete set of permutations.

Would such a function be of interest for numpy? If so, I could start with a PR proposal.

Note: The standard library provides itertools.permutations(...) but that's significantly slower that a pure numpy implementation.

[mattip]

How do you know itertools.permutations is slower? What are you comparing it too?

[timhoffm]

I've used the following to create an array using itertools:

indices = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
permutations = np.array(list(itertools.permutations(indices)), dtype=np.uint8)

Note: Not instantiating the intermediate list would save some time, but numpy.fromiter only creates 1D arrays. Is there a way to create a 2D array from an iterable yielding sequences?

I compared to the numpy solutions proposed at https://stackoverflow.com/questions/64291076/generating-all-permutations-efficiently. They are about a factor 50 faster, which is not an uncommon Pure-Python -> numpy speedup given that permutations must construct the sequence object for each yielded permutation.

Edit: This reference may be interesting for a fast algorithm as well: https://en.wikipedia.org/wiki/Heap%27s_algorithm#cite_note-3

[seberg]

FWIW, you can use fromiter if you use a structured dtype, but that only gives you a ~30% speedup or so (might be slower also due to the structured part being slow to deal with):

def permutations(nvals, dtype=np.intp):
    dt = np.dtype([("", dtype)] * nvals)
    arr = np.fromiter(itertools.permutations(range(nvals)), dtype=dt)
    return arr.view(dtype)

I think @stefanv also wanted this always and I could see adding it. There is maybe also a bit of a question of whether you need to have some option to gracefully deal with a huge number of permutations (since this can quickly explode in size).

[timhoffm]

There is maybe also a bit of a question of whether you need to have some option to gracefully deal with a huge number of permutations (since this can quickly explode in size).

This is a good idea since people may not be aware that this becomes a significant issue with increasing N. One option would be a default size limit (in N or memory usage), and if people want bigger arrays, they have to consciously opt in. An advanced feature would be some batch mode where you receive a generator for a sequence of subarrays of the permutations. But to keep things simple that's likely not for the first shot.

[OyiboRivers]

@seberg would you consider including other combinatorial indexing functions as well?

Some options would be:

• permutations
• permutations_with_replacement
• combinations
• combinations_with_replacement
• cartesian_product
• ...

@timhoffm, your feature request is very useful.
It would be great and very convenient if numpy could offer combinatorial indexing functions.

I'm skeptical about the expectations for a substantial increase in performance.
Most benchmarks should be taken with a grain of salt. On stackoverflow the apparent speed difference between functions, where one is claimed to be 50 times faster than itertools, may stem from differences in functionality:

• The supposedly faster function does not enumerate permutations lexicographically, potentially requiring additional sorting.
• It does not generate permutations of arbitrary arrays, you have to use fancy indexing additionally.
• It is using a dtype of uint8 instead of int64, which might be sufficient in many cases but maybe not all.
• Additionally, it does not support selecting k items from a set of n items.

Here's a benchmark that aims to make a fairer comparison between faster_permutations and itertools.permutations. The uint8 version is approximately 5 times faster than itertools, while the int64 version offers similar performance on reduced functionality.

import math
import itertools
from timeit import timeit
import numpy as np

def faster_permutations(n, dtype=np.uint8):
    """https://stackoverflow.com/questions/64291076/generating-all-permutations-efficiently"""
    perms = np.empty((math.factorial(n), n), dtype=dtype, order='F')
    perms[0, 0] = 0
    rows_to_copy = 1
    for i in range(1, n):
        perms[:rows_to_copy, i] = i
        for j in range(1, i + 1):
            start_row = rows_to_copy * j
            end_row = rows_to_copy * (j + 1)
            splitter = i - j
            perms[start_row: end_row, splitter] = i
            perms[start_row: end_row, :splitter] = perms[:rows_to_copy, :splitter]
            perms[start_row: end_row, splitter + 1:i + 1] = perms[:rows_to_copy, splitter:i]
        rows_to_copy *= i + 1
    return perms

def apply_faster_permutations(arr, dtype=np.uint8):
    """
    Apply faster_permutations does some additional work which itertools takes care of, too.
    - get index permutations from faster_permutations
    - Make sure enumeration of permutations is done lexicographically
    - Use fancy indexing to generate permutations of the actual array
    """
    n = arr.size
    perms = faster_permutations(n, dtype)
    idx = np.lexsort(perms.T[::-1])
    perms_sorted = perms[idx]
    return arr[perms_sorted]

def permutations_itertools(arr: np.ndarray, k: int = None) -> np.ndarray:
    k = arr.size if k is None else k
    dtype = np.dtype((arr.dtype, k))
    return np.fromiter(itertools.permutations(arr, k), dtype=dtype, count=-1)

# setup
n = 10  # The number of items within the set
arr = np.arange(1, n+1)   # The set

# Check if results are the same
res_itertools = permutations_itertools(arr)
res_faster = apply_faster_permutations(arr)
print('Is result equal? => ', np.array_equal(res_itertools, res_faster))
# Is result equal? =>  True

# measure performance
time_itertools = timeit("permutations_itertools(arr)", globals=globals(), number=10)
time_faster = timeit("apply_faster_permutations(arr)", globals=globals(), number=10)
time_faster_int64 = timeit("apply_faster_permutations(arr, dtype=np.int64)", globals=globals(), number=10)

print(f'N={n}: faster_permutations vs itertools is {round(time_itertools/time_faster, 1)} times faster')
print(f'N={n}: faster_permutations_int64 vs itertools is {round(time_itertools/time_faster_int64, 1)} times faster')
# N=10: faster_permutations vs itertools is 5.4 times faster
# N=10: faster_permutations_int64 vs itertools is 1.1 times faster

[charris]

I would be inclined to make all those functions part of a separate package. Knuth has covered the topic in volumes 4A, 4B, 4C, 4D, 4E, and 4F of TAOCP, of which 4A and 4B have been published :) I think 4A covers what you want.

[timhoffm]

I feel permutations() would be a good addition, but the others are too special.

This is opinionated, but I regard permutations() more as anarray creation routine, rather than an algorithm.

The first priority is ease of use and expressiveness. Better performance than itertools is nice but only secondary. - It's likely also only relevant for a small range of elments (N~10-14). For lower numbers, every approach is fast enough (I don't expect people to call this many times in a loop). The combinatorical explosion limits towards the high numbers because of memory usage.
If one wanted to keep it really simple, I find even just wrapping itertools would be valuable because of a concise API in the numpy world.

[ngoldbaum]

I agree this definitely makes sense in an external package. Starting with an external implementation that we could think about importing into numpy would also motivate the discussion about including it in numpy.

Note that a cython implementation would probably be less appealing than a C or C++ implementation - cython produces very large code sizes and I don't know if there's a ton of appetite for adding more cython code to numpy outside of the random module.

[timhoffm]

I agree this definitely makes sense in an external package

Are we talking about permutations() or a broader set of combinatorial functions like proposed above?

I would be willing to contribute a permutations() function to numpy, but I don't have motivation or capacity to start a new package.

Note that a cython implementation would probably be less appealing than a C or C++ implementation

As mentioned above, I see the first priority in ease of use and expressiveness. Getting the API and tests right would be the first step for me: e.g. to be discussed whether one wants one or both of permutations of the first N integers (permutations(N)) and permutations of a given array (permutations([2, 4, 6])). The initial implementation may even be in terms of itertools.permutation. Improving performance is a second step. Whether that works in pure numpy reasonably well or needs C/C++/cython is an optimization detail.

[OyiboRivers]

@seberg , is it possible that using itertools.permutations for large sets with n>10 causes a bottleneck due to the conversion of generated tuples into a numpy array?

The following code calculates permutations for sets of size n.

• permutations_list() converts the generated tuples into a python list.
• permutations_list_asarray() converts the generated tuples into a python list and then into a numpy array.
• permutations_fromiter() converts the generated tuples into a numpy array using np.fromiter.

import timeit
import itertools as it
import numpy as np

def permutations_list(n: int, k: int = None) -> np.ndarray:
    """Generate permutations using itertools and convert to list."""
    return list(it.permutations(np.arange(n), k))

def permutations_list_asarray(n: int, k: int = None) -> np.ndarray:
    """Generate permutations using itertools and convert to numpy array."""
    return np.asarray(list(it.permutations(np.arange(n), k)))

def permutations_fromiter(n: int, k: int = None) -> np.ndarray:
    """Generate permutations using itertools and convert to numpy array using fromiter."""
    dtype = np.intp
    return np.fromiter(it.permutations(np.arange(n, dtype=dtype), k), dtype=np.dtype((dtype, k)))

for n, k in [(8,8), (9,9), (10,10), (11,11)]:
    time_list = timeit.timeit(stmt="permutations_list(n, k)", globals=globals(), number=1)
    time_list_asarray = timeit.timeit(stmt="permutations_list_asarray(n, k)", globals=globals(), number=1)
    time_from_iter = timeit.timeit(stmt="permutations_fromiter(n, k)", globals=globals(), number=1)
    print(f'Permutations using itertools for n={n}, k={k} => rows={len(permutations_list(n, k)):,} columns={k}')
    print("list:                    {:.3f} seconds".format(time_list))
    print("np.array from list:      {:.3f} seconds".format(time_list_asarray))
    print("np.array using fromiter: {:.3f} seconds".format(time_from_iter))
    print()

# Permutations using itertools for n=8, k=8 => rows=40,320 columns=8
# list:          0.005 seconds
# list to array: 0.034 seconds
# fromiter:      0.029 seconds

# Permutations using itertools for n=9, k=9 => rows=362,880 columns=9
# list:          0.051 seconds
# list to array: 0.337 seconds
# fromiter:      0.256 seconds

# Permutations using itertools for n=10, k=10 => rows=3,628,800 columns=10
# list:          0.538 seconds
# list to array: 3.676 seconds
# fromiter:      2.654 seconds

# Permutations using itertools for n=11, k=11 => rows=39,916,800 columns=11
# list:          5.912 seconds
# list to array: 44.436 seconds
# fromiter:      30.511 seconds

The conversion of generated tuples into a Python list remains stable, while conversions into numpy arrays deteriorate beyond set size n>10.
Can you confirm this observed behavior?