Список случайных подмножеств указанного размера

Question

Я хочу вычислить список случайных подмножеств массива, где порядок в подмножествах является случайным, и каждый элемент в их соответствующем подмножестве уникален, и я хочу сделать это эффективно (с точки зрения времени и пространства).

Например, массив [1,2,3,4,5] с некоторыми параметрами k=3 и n=5 должен создавать матрицу

[[4,3,1],
 [1,2,5],
 [2,4,5],
 [3,2,5],
 [2,3,1]]

, то есть мы получаем n=5 списки с k=3 random уникальных элементов из массива.

Мне бы хотелось, чтобы это было как можно быстрее, без создания гигантской таблицы поиска возможных комбинаций, поскольку время и пространство имеют существенное значение.

Я пытался сделать это с numpy.random.choice(array,n*k).reshape((n,k)), который почти делает то, что я хочу, за исключением части уникальности. Я остановился на следующем

subsets = numpy.zeros(n).reshape((n,1))
subsets = numpy.apply_along_axis(lambda x: numpy.random.choice(array, k, replace=False),1, subsets)

Однако, так как это не просто тупица, это медленно. Слишком медленно для моего приложения. Есть ли какой-нибудь способ улучшить время выполнения и, возможно, добиться этого с помощью чисто командных команд?

Буду признателен за любую помощь.

Answer 1

при условии, что k будет намного меньше, чем n, вы можете реализовать выборку без замены самостоятельно:

idx = np.random.randint([5,4,3], size=[5,3])
idx[:,2:] += idx[:,2:] >= idx[:,1,None]
idx[:,1:] += idx[:,1:] >= idx[:,0,None]
np.array([1,2,3,4,5])[arr]

, вставив в функцию, это может быть:

def sample_noreplace(arr, n, k):
    assert k <= len(arr)
    idx = np.random.randint(len(arr) - np.arange(k), size=[n, k])
    for i in range(k-1, 0, -1):
        idx[:,i:] += idx[:,i:] >= idx[:,i-1,None]
    return np.array(arr)[idx]

для запуска sample_noreplace([1,2,3,4,5], 10000, 3) требуется ~ 1 мс, в то время как для кода OP требуется ~ 800 мс

, чтобы показать, что выборка выполняется равномерно, и в таблицу были бы выведены большие объемы:

arr = sample_noreplace(range(5), 1_000_000, 3)

# summarise
dict(zip(*np.unique(arr, return_counts=True)))

который дает мне: {0: 600288, 1: 599656, 2: 600494, 3: 599233, 4: 600329}, который выглядит довольно равномерно. затем мы можем проверить в каждой позиции:

out = np.zeros([3, 5])
for i in range(3):
    n, c = np.unique(arr[:,i], return_counts=True)
    out[i, n] = c

, что дает мне таблицу, похожую на:

array([[199936., 199701., 200106., 199843., 200414.],
       [200227., 200044., 200345., 199897., 199487.],
       [200125., 199911., 200043., 199493., 200428.]])

, которая также выглядит равномерно.

Answer 2

Учитывая ваши числа, один хороший подход - просто рисовать с заменой и отбрасывать розыгрыши с дубликатами.

Ниже приведены временные интервалы для рисования 5000 выборок по 10% или 10, в зависимости от того, что меньше из массивов различного размера. pp_overdraw рисует слишком много сэмплов и сбрасывает, pp_fillin обращается к точному числу, перерисовывает плохие сэмплы, пока не останется ни одного. Pythran - это скомпилированное решение. Поскольку OP запрашивает чистый numpy, он здесь только для справки.

Код:

import pyflip
import numpy as np
from simple_benchmark import BenchmarkBuilder, MultiArgument

B = BenchmarkBuilder()

@B.add_function()
def pythran(A,k,n):
    assert len(A) >= k
    if A.dtype not in (float,int) or A.ndim > 1:
        return A[pyflip.without_repl(np.arange(len(A)),k,n)]
    else:
        return pyflip.without_repl(A,k,n)

from scipy import stats

@B.add_function()
def pp_overdraw(A,k,n):
    N = len(A)
    p = np.linspace(1-(k-1)/N,1-1/N,k-1).prod()
    nn = int(n/p + 4*np.sqrt(n-n*p)) + 1
    out = np.random.randint(0,N,(nn,k))
    os = np.sort(out,1)
    valid = (os[:,:-1] != os[:,1:]).all(1)
    validx, = valid.nonzero()
    while len(validx)<n: # very unlikely
        replace = np.random.randint(0,N,(nn-len(validx),k))
        rs = np.sort(replace,1)
        rvalid = (rs[:,:-1] != rs[:,1:]).all(1)
        out[~valid,:] = replace
        valid[~valid] = rvalid
        validx, = valid.nonzero()
    return A[out[validx[:n]]]

@B.add_function()
def pp_fillin(A,k,n):
    N = len(A)
    out = np.random.randint(0,N,(n,k))
    os = np.sort(out,1)
    valid = (os[:,:-1] != os[:,1:]).all(1)
    validx, = valid.nonzero()
    while len(validx)<n:
        replace = np.random.randint(0,N,(n-len(validx),k))
        rs = np.sort(replace,1)
        rvalid = (rs[:,:-1] != rs[:,1:]).all(1)
        out[~valid,:] = replace
        valid[~valid] = rvalid
        validx, = valid.nonzero()
    return A[out]


@B.add_function()
def OP(A,k,n):
    subsets = np.zeros(n).reshape((n,1))
    subsets = np.apply_along_axis(lambda x: np.random.choice(A, k, replace=False),1, subsets)

@B.add_function()
def Sam_Mason(A,k,n):
    assert k <= len(A)
    idx = np.random.randint(len(A) - np.arange(k), size=[n, k])
    for i in range(k-1, 0, -1):
        idx[:,i:] += idx[:,i:] >= idx[:,i-1,None]
    return np.array(A)[idx]

@B.add_arguments('array size')
def argument_provider():
    for exp in range(4, 15):
        sz = 2**exp
        A = np.arange(sz)
        yield sz, MultiArgument([A,min(10,sz//10+1),5000])

r = B.run()
r.plot()

import pylab
pylab.savefig('norepl.png')

модуль пифрана. В принципе, он работает на Python, но это будет мучительно медленно. Лучше скомпилировать с

pythran -O3 pyflip.py

file <pyflip.py>

import numpy as np

#pythran export without_repl(int[:],int,int)
#pythran export without_repl(float[:],int,int)

def without_repl(A,k,n):
    N = A.size
    out = np.empty((n,k),A.dtype)
    A = A.copy()
    for i in range(n):
        for j in range(k):
            l = np.random.randint(j,N)
            out[i,j] = A[l]
            A[l] = A[j]
            A[j] = out[i,j]
    return out

Answer 3

N = [1, 2, 3, 4, 5]
k = 3
n = 5

arr = np.array([N] * n)

arr
array([[1, 2, 3, 4, 5],
       [1, 2, 3, 4, 5],
       [1, 2, 3, 4, 5],
       [1, 2, 3, 4, 5],
       [1, 2, 3, 4, 5]])

[np.random.shuffle(i) for i in arr]

arr
array([[3, 5, 1, 2, 4],
       [5, 1, 3, 2, 4],
       [3, 5, 2, 4, 1],
       [4, 5, 1, 2, 3],
       [1, 5, 3, 2, 4]])

arr[:, :3]
array([[3, 5, 1],
       [5, 1, 3],
       [3, 5, 2],
       [4, 5, 1],
       [1, 5, 3]])