【问题标题】:Numpy list of random subsets of specified size指定大小的随机子集的 Numpy 列表
【发布时间】:2025-11-23 14:40:02
【问题描述】:

我想计算一个数组的随机子集列表,其中子集中的顺序是随机的,并且它们各自子集中的每个元素都是唯一的,我想有效地这样做(在时间和空间方面)。

例如数组[1,2,3,4,5] 带有一些参数k=3n=5 应该产生一个矩阵

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

即,我们从数组中获得带有k=3 随机唯一元素的n=5 列表。

我希望这尽可能快,而不是创建一个巨大的可能组合查找表,因为时间和空间都是至关重要的。

我尝试使用 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)

但是,由于这不是纯粹的 numpy,所以这很慢。我的申请太慢了。有什么方法可以改进运行时,也许可以用纯 numpy 命令来实现?

任何帮助将不胜感激。

【问题讨论】:

  • 不是一个很好的解决方案,但使用地图更快。 np.array(list(map(lambda x: np.random.choice(array, 3, replace=False),range(5)))).
  • 您期望nk 以及输入数组的大小是多少?
  • @PaulPanzer 目前,len(array) 在 500 左右,n 在 1000 到 5000 之间,k 在 2 到 10 之间
  • @PaulPanzer 好点,也应该问这个问题!请注意 OP,如果 k 远低于 sqrt(len(array)) 则使用不同的算法会变得有用,因为选择重复的机会相对较低。参见例如birthday problem 表示 sqrt 的来源
  • @SamMason 如果您有兴趣,我已经添加了一个解决方案的代码,该解决方案适用于任何参数(O(len(array)+nk)),但依赖于 pythran(numba 也可以)来提高速度,因为它是循环的。 - - 顺便提一句。您认为您的算法适用于哪个政权?它在 k 中是二次的,所以它也必须在 k

标签: python arrays numpy random


【解决方案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]

这需要 ~1ms 来运行 sample_noreplace([1,2,3,4,5], 10000, 3) 而 OPs 代码需要 ~800ms

为了表明这是均匀采样,我将大量输出制成表格:

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.]])

看起来也很统一。

【讨论】:

  • 非常干净的解决方案,谢谢。我也很欣赏概率分析!
【解决方案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')

pythran 模块。原则上,它在 Python 中运行,但这会非常缓慢。更好地编译

pythran -O3 pyflip.py

文件&lt;pyflip.py&gt;

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

【讨论】:

    【解决方案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]])
    

    【讨论】:

    • 输出矩阵中的第二个数组具有非唯一成员。 1 in [1,5,1]
    • 如果 N 很大但 k 很小,这个新版本在空间方面可能真的很昂贵。
    最近更新 更多