Hadamard product for each unique pair of columns in numpy array

Question

Using Python (3.7.7) and numpy (1.17.4), I am working with medium sized 2d numpy arrays (from 5000x80 up to 200,000x120). For a given array, I want to calculate the Hadamard product between all possbible uniqe pairs of column-vectors of that array.

I have:

    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]

and I want to get:

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

I already have a solution from a colleague using np.kron which I adapated a bit:

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())

Note that hadamard_A is what I want, but transposed (which is also what I want for further processing). Also, ind1_ (ind2_) gives the indices for the objects which feature as the first (second) element in the pair for which the hadamard product is calculated. I need those as well.

However, I feel this code is too inefficient: it takes to long and since I call this function several times during my algorithm, I was wondering whether there is a cleverer solution? Am I overlooking some numpy/scipy tools I could cleverly combine for this task?

Thanks all! :)

Answer 1

Approach #1

Simplest one with 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]])

Approach #2

Memory-efficient one for large arrays -

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:]

Approach #3

Masking based one -

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)

Approach #4

Extend approach #2 for a numba one -

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

Then, leverage parallel processing (as pointed out in comments by @max9111), like so -

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

Benchmarking

Using benchit package (few benchmarking tools packaged together; disclaimer: I am its author) to benchmark proposed solutions.

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')

Conclusion : Numba ones seem to be doing pretty well and app2 among NumPy ones.