查找数组元素的所有可能组合。不是乘积。

Question

我研究了这个问题，人们一直建议使用 np.meshgrid() 来找到一个数组的所有可能组合。但是问题在于，np.meshgrid() 产生的不是组合，而是乘积（类似于itertools.product()）。

在组合中，元素不重复且无序。

arr = np.arange(5)
r = 3

以下是组合的样子

np.array(
    list(itertools.combinations(arr, r))
)

>>>   [[0, 1, 2],
       [0, 1, 3],
       [0, 1, 4],
       [0, 2, 3],
       [0, 2, 4],
       [0, 3, 4],
       [1, 2, 3],
       [1, 2, 4],
       [1, 3, 4],
       [2, 3, 4]]

以下不是组合：

np.array(
    list(itertools.product(arr, arr, arr))
)

>>>   [[0, 0, 0],
       [0, 0, 1],
       [0, 0, 2],
       [0, 0, 3],
       [0, 0, 4],
       [0, 1, 0],
       [0, 1, 1],
       [0, 1, 2],
       ....,
       [4, 3, 2],
       [4, 3, 3],
       [4, 3, 4],
       [4, 4, 0],
       [4, 4, 1],
       [4, 4, 2],
       [4, 4, 3],
       [4, 4, 4]])

np.array(
    np.meshgrid(arr, arr, arr)
).transpose([2, 1, 3, 0]).reshape(-1, r)

>>>   [[0, 0, 0],
       [0, 0, 1],
       [0, 0, 2],
       [0, 0, 3],
       [0, 0, 4],
       [0, 1, 0],
       [0, 1, 1],
       [0, 1, 2],
       ....,
       [4, 3, 2],
       [4, 3, 3],
       [4, 3, 4],
       [4, 4, 0],
       [4, 4, 1],
       [4, 4, 2],
       [4, 4, 3],
       [4, 4, 4]])

对于r = 2，我找到了一种简洁的组合方法。

np.array(
    np.triu_indices(len(arr), 1)
).T

>>>   [[0, 1],
       [0, 2],
       [0, 3],
       [0, 4],
       [1, 2],
       [1, 3],
       [1, 4],
       [2, 3],
       [2, 4],
       [3, 4]]

但我很难找到任何适用于 r > 2 的向量化方法。

注意： 即使我的数组不是[0, 1, 2, 3, 4]，我也可以使用上面的答案作为索引。

如果有助于想象,

当 r = 2 时，所需答案是一个大小为 len(arr) 的正方形矩阵的右上角三角形的索引，忽略对角线。

当 r = 3 时，所需答案是一个大小为 len(arr) 的三维数组的顶部右上角上四面体（在图像中为中间位置），忽略三维等效的对角线。

- Hammad

---

"排列组合"? from itertools import permutations list(permutations(arr, 3)) 可能。 - Scott Boston

排列与组合是不同的，我正在寻找一种向量化实现，所以 itertools 不可行。 - hammi

@Hammad 你需要向量化它的原因是什么？ - Ehsan

@Ehsan 我有一个大数据集，所以我想避免循环和生成器。 - hammi

Answer 1

这与您关于3D的想法相似:

n = 5
r = 3
a = np.argwhere(np.triu(np.ones((n,)*r),1))
a[a[:,0]<a[:,1]]

输出：

[[0 1 2]
 [0 1 3]
 [0 1 4]
 [0 2 3]
 [0 2 4]
 [0 3 4]
 [1 2 3]
 [1 2 4]
 [1 3 4]
 [2 3 4]]

对于4-D（以此类推）的情况，您可以像这样扩展（性能不确定）：

n = 5
r = 4
a = np.argwhere(np.triu(np.ones((n,)*r),1))
a[(a[:,0]<a[:,1]) & (a[:,1]<a[:,2])]

输出：

[[0 1 2 3]
 [0 1 2 4]
 [0 1 3 4]
 [0 2 3 4]
 [1 2 3 4]]

如果你的目标是速度，那么Itertools似乎更快：

def m1(n):
  r = 3
  a = np.argwhere(np.triu(np.ones((n,)*r),1))
  return a[a[:,0]<a[:,1]]

def m2(n):
  r = 3
  return combinations(np.arange(n),r)

in_ = [5,10,100,200]

- Ehsan

---

谢谢你的回答，itertools实现之所以更快，是因为你只返回生成器，这些生成器稍后需要转换成数组并花费大量时间。 - hammi

Answer 2

看一下代码：

def triu_indices(n, k=0, m=None):
    return nonzero(~tri(n, m, k=k-1, dtype=bool))

def tri(N, M=None, k=0, dtype=float):
    m = greater_equal.outer(arange(N, dtype=_min_int(0, N)),
                        arange(-k, M-k, dtype=_min_int(-k, M - k)))

所以对于一个正方形：

In [75]: np.greater_equal.outer(np.arange(5),np.arange(5))                                           
Out[75]: 
array([[ True, False, False, False, False],
       [ True,  True, False, False, False],
       [ True,  True,  True, False, False],
       [ True,  True,  True,  True, False],
       [ True,  True,  True,  True,  True]])
In [76]: np.argwhere(~np.greater_equal.outer(np.arange(5),np.arange(5)))                             
Out[76]: 
array([[0, 1],
       [0, 2],
       [0, 3],
       [0, 4],
       [1, 2],
       [1, 3],
       [1, 4],
       [2, 3],
       [2, 4],
       [3, 4]])

或使用ix_创建可广播的范围：

In [77]: I,J = np.ix_(np.arange(5),np.arange(5))                                                     
In [78]: I,J                                                                                         
Out[78]: 
(array([[0],
        [1],
        [2],
        [3],
        [4]]),
 array([[0, 1, 2, 3, 4]]))
In [79]: I>=J                                                                                        
Out[79]: 
array([[ True, False, False, False, False],
       [ True,  True, False, False, False],
       [ True,  True,  True, False, False],
       [ True,  True,  True,  True, False],
       [ True,  True,  True,  True,  True]])
In [80]: I<J                                                                                         
Out[80]: 
array([[False,  True,  True,  True,  True],
       [False, False,  True,  True,  True],
       [False, False, False,  True,  True],
       [False, False, False, False,  True],
       [False, False, False, False, False]])

将其概括为三维情况：

In [90]: I,J,K = np.ix_(np.arange(5),np.arange(5),np.arange(5))                                      
In [91]: np.argwhere((I<J) & (J<K))                                                                  
Out[91]: 
array([[0, 1, 2],
       [0, 1, 3],
       [0, 1, 4],
       [0, 2, 3],
       [0, 2, 4],
       [0, 3, 4],
       [1, 2, 3],
       [1, 2, 4],
       [1, 3, 4],
       [2, 3, 4]])

这创建了一个(5,5,5)的数组，然后查找True。因此，它的时间效率不如itertools出乎意料。

In [92]: %%timeit 
    ...: I,J,K = np.ix_(np.arange(5),np.arange(5),np.arange(5))   
    ...: np.argwhere((I<J) & (J<K))                                                                                             
60.4 µs ± 97 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [99]: timeit np.array(list(itertools.combinations(np.arange(5), 3)))                              
20.9 µs ± 1.74 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

- hpaulj

Answer 3

正在运行一些针对答案的基准测试

import numpy as np
import itertools
import perfplot

itertools 实现

def func0(n, r):
    
    r = np.array(
        list(itertools.combinations(range(n), r))
    )
    
    return r

@hpaulj的实现

def func1(n, r):
    
    prod = np.argwhere(np.ones(r*[n], dtype= bool))
    mask = True
    
    for i in range(r-1): mask = mask & (prod[:, i] < prod[:, i+1])
        
    return prod[mask]

@Ehsan的实现

def func2(n, r):
    
    rngs = (np.arange(n), )*r
    indx = np.ix_(*rngs)
    mask = True
    
    for i in range(r-1): mask = mask & (indx[i] < indx[i+1])
        
    return np.argwhere(mask)

性能与 n 相关

bench_n = perfplot.bench(
    n_range= range(4, 50),
    setup  = lambda n: n,
    kernels= [
        lambda n: func0(n, 4),
        lambda n: func1(n, 4),
        lambda n: func2(n, 4)
    ],
    labels = ['func0', 'func1', 'func2']
)

bench_n.show()

性能与 r 相关

bench_r = perfplot.bench(
    n_range= range(2, 8),
    setup  = lambda r: r,
    kernels= [
        lambda r: func0(10, r),
        lambda r: func1(10, r),
        lambda r: func2(10, r)
    ],
    labels = ['func0', 'func1', 'func2']
)

bench_r.show()

- hammi