numpy 数组中每对唯一的列的 Hadamard 产品

Question

使用 Python (3.7.7) 和 numpy (1.17.4)，我正在处理中等大小的 2d numpy 数组（从 5000x80 到 200,000x120）。对于给定的数组，我想计算该数组的所有可能的唯一列向量对之间的 Hadamard 乘积。

我有：

    A           A
[a,b,c,d]   [a,b,c,d]
[1,2,3,4]   [1,2,3,4]
[4,5,6,7] * [4,5,6,7]
[7,8,9,1]   [7,8,9,1]

我想得到：

[a*b, ac,  ad,  bc,  bd,  cd]
[ 2.,  3.,  4.,  6.,  8., 12.]
[20., 24., 28., 30., 35., 42.]
[56., 63.,  7., 72.,  8.,  9.]

我已经从同事那里得到了一个使用 np.kron 的解决方案，我对其进行了一些修改：

def hadamard_kron(A: np.ndarray) -> :
    """Returns the hadamard products of all unique pairs of all columns, 

    and return indices signifying which columns constitute a given pair.
    """

    n = raw_inputs.shape[0]
    ind1 = (np.kron(np.arange(0, n).reshape((n, 1)), np.ones((n, 1)))).squeeze().astype(int)
    ind2 = (np.kron(np.ones((n, 1)), np.arange(0, n).reshape((n, 1)))).squeeze().astype(int)
    xmat2 = np.kron(raw_inputs, np.ones((n, 1))) * np.kron(np.ones((n, 1)), raw_inputs)

    hadamard_inputs =  xmat2[ind2 > ind1, :]
    ind1_ = ind1[ind1 < ind2]
    ind2_ = ind2[ind1 < ind2]
    return hadamard_A, ind1_, ind2_

hadamard_A, first_pair_members, second_pair_members = hadamard_kron(a.transpose())

注意 hadamard_A 是我想要的，但是转置了（这也是我想要进一步处理的）。此外，ind1_ (ind2_) 给出了对象的索引，这些对象的特征是计算 hadamard 乘积的对中的第一个（第二个）元素。我也需要这些。

但是，我觉得这段代码效率太低：需要很长时间，而且我在算法中多次调用这个函数，我想知道是否有更聪明的解决方案？我是否忽略了一些 numpy/scipy 可以巧妙组合用于此任务的工具？

谢谢大家！ :)

Answer 1

方法 #1

最简单的 np.triu_indices -

In [45]: a
Out[45]: 
array([[1, 2, 3, 4],
       [4, 5, 6, 7],
       [7, 8, 9, 1]])

In [46]: r,c = np.triu_indices(a.shape[1],1)

In [47]: a[:,c]*a[:,r]
Out[47]: 
array([[ 2,  3,  4,  6,  8, 12],
       [20, 24, 28, 30, 35, 42],
       [56, 63,  7, 72,  8,  9]])

方法 #2

Memory-efficient 一个用于大数组 -

m,n = a.shape
s = np.r_[0,np.arange(n-1,-1,-1).cumsum()]
out = np.empty((m, n*(n-1)//2), dtype=a.dtype)
for i,(s0,s1) in enumerate(zip(s[:-1], s[1:])):
    out[:,s0:s1] = a[:,i,None] * a[:,i+1:]

方法 #3

基于掩蔽的一个 -

m,n = a.shape
mask = ~np.tri(n,dtype=bool)
m3D = np.broadcast_to(mask, (m,n,n))

b1 = np.broadcast_to(a[...,None], (m,n,n))
b2 = np.broadcast_to(a[:,None,:], (m,n,n))
out = (b1[m3D]* b2[m3D]).reshape(m,-1)

方法 #4

将方法 #2 扩展为 numba 一个 -

from numba import njit

def numba_app(a):
    m,n = a.shape
    out = np.empty((m, n*(n-1)//2), dtype=a.dtype)
    return numba_func(a,out,m,n)

@njit
def numba_func(a,out,m,n):
    for p in range(m):
        I = 0
        for i in range(n):
            for j in range(i+1,n):
                out[p,I] = a[p,i] * a[p,j]
                I += 1
    return out

然后，利用 parallel 处理（正如@max9111 在评论中指出的那样），就像这样 -

from numba import prange

def numba_app_parallel(a):
    m,n = a.shape
    out = np.empty((m, n*(n-1)//2), dtype=a.dtype)
    return numba_func_parallel(a,out,m,n)

@njit(parallel=True)
def numba_func_parallel(a,out,m,n):
    for p in prange(m):
        I = 0
        for i in range(n):
            for j in range(i+1,n):
                out[p,I] = a[p,i] * a[p,j]
                I += 1
    return out

基准测试

使用 benchit 包（几个基准测试工具打包在一起；免责声明：我是它的作者）对建议的解决方案进行基准测试。

import benchit
in_ = [np.random.rand(5000, 80), np.random.rand(10000, 100), np.random.rand(20000, 120)]
funcs = [ehsan, app1, app2, app3, numba_app, numba_app_parallel]
t = benchit.timings(funcs, in_, indexby='shape')
t.rank()
t.plot(logx=False, save='timings.png')

结论：Numba 个似乎做得很好，app2 个在 NumPy 中。

Answer 2

Divakar 第一种方法的另一种等效方法：

r,c = np.triu_indices(A.shape[1],1)
np.einsum('ij,ik->ijk',A,A)[:,r,c]

输出：

[[ 2  3  4  6  8 12]
 [20 24 28 30 35 42]
 [56 63  7 72  8  9]]

numpy 数组中每对唯一的列的 Hadamard 产品

Hadamard product for each unique pair of columns in numpy array

python

arrays

numpy

unique

基准测试