numpy数组中每对唯一列的Hadamard乘积

Question

我使用Python（3.7.7）和numpy（1.17.4）处理中等大小的2d numpy阵列（从5000x80到200000x120）。对于给定的数组，我想计算该数组的所有列向量对之间的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

另一种与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]]

Answer 2

方法#1

最简单的带^{}-

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

用于大型阵列的高效内存-

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扩展为numba1-

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

基准测试

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

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的人则是裸体人