Efficiently compute pairwise equal for NumPy arrays

Question

Given two NumPy arrays, say:

import numpy as np
import numpy.random as rand

n = 1000
x = rand.binomial(n=1, p=.5, size=(n, 10))
y = rand.binomial(n=1, p=.5, size=(n, 10))

Is there a more efficient way to compute X in the following:

X = np.zeros((n, n))
for i in range(n):
    for j in range(n):
        X[i, j] = 1 * np.all(x[i] == y[j])

Yes. Look into [broadcasting](https://docs.scipy.org/doc/numpy-1.13.0/user/basics.broadcasting.html), and then try `(x[:, None, :] == y).all(axis=-1)`. However, this is time efficient but not particularly space efficient (though in this particular case we're only talking about needing on the order of 10MB of space). Are you primarily worried about space efficiency or time efficiency? How large will `n` get in practice? — Mark Dickinson, Jan 08 '18 at 20:26
FYI: The default data type created by `X = np.zeros((n, n))` is floating point. Do you really want `X` to be an array of floating point values? If all you care about is equality, you could make it an array of small integers (`X = np.zeros((n, n), dtype=np.int8)`) or an array of boolean values (`X = np.zeros((n, n), dtype=np.bool_)`). — Warren Weckesser, Jan 08 '18 at 21:18
@MarkDickinson I think `n` will be at most `10^6` in practice. — p-value, Jan 15 '18 at 17:36
Hmm. Then you're going to have trouble getting `X` into memory, unless you have a terabyte or more or RAM. What do you do with `X` *after* you've created it? — Mark Dickinson, Jan 15 '18 at 20:14
That's a good point... I'm hoping that `X` will be a sparse matrix. Afterwards I would need to perform a number of mat-vec products using `X`. — p-value, Jan 15 '18 at 20:43

Divakar · Accepted Answer · 2018-01-08T20:48:19.203

Approach #1 : Input arrays with 0s & 1s

For input arrays with 0s and 1s only, we can reduce each of their rows to scalars and hence the input arrays to 1D and then leverage broadcasting, like so -

n = x.shape[1]        
s = 2**np.arange(n)
x1D = x.dot(s)
y1D = y.dot(s)
Xout = (x1D[:,None] == y1D).astype(float)

Approach #2 : Generic case

For a generic case, we can use views -

# https://stackoverflow.com/a/45313353/ @Divakar
def view1D(a, b): # a, b are arrays
    a = np.ascontiguousarray(a)
    b = np.ascontiguousarray(b)
    void_dt = np.dtype((np.void, a.dtype.itemsize * a.shape[1]))
    return a.view(void_dt).ravel(),  b.view(void_dt).ravel()

x1D, y1D = view1D(x, y)
Xout = (x1D[:,None] == y1D).astype(float)

Runtime test

# Setup
In [287]: np.random.seed(0)
     ...: n = 1000
     ...: x = rand.binomial(n=1, p=.5, size=(n, 10))
     ...: y = rand.binomial(n=1, p=.5, size=(n, 10))

# Original approach
In [288]: %%timeit
     ...: X = np.zeros((n, n))
     ...: for i in range(n):
     ...:     for j in range(n):
     ...:         X[i, j] = 1 * np.all(x[i] == y[j])
1 loop, best of 3: 4.69 s per loop

# Approach #1
In [290]: %%timeit
     ...: n = x.shape[1]        
     ...: s = 2**np.arange(n)
     ...: x1D = x.dot(s)
     ...: y1D = y.dot(s)
     ...: Xout = (x1D[:,None] == y1D).astype(float)
1000 loops, best of 3: 1.42 ms per loop

# Approach #2
In [291]: %%timeit
     ...: x1D, y1D = view1D(x, y)
     ...: Xout = (x1D[:,None] == y1D).astype(float)
100 loops, best of 3: 18.5 ms per loop

Thanks for the comprehensive answer! I'm new to the concept of views, and I'm having a little trouble understanding Approach #2. Do you mind providing some more explanation for the approach and the reason why it's so efficient? — p-value, Jan 09 '18 at 05:06
@Hilbert The idea has been explored before and has some comments in this post - https://stackoverflow.com/a/43481447/. — Divakar, Jan 09 '18 at 06:55
Sorry I'm new to using views... I looked through the other posts but I still don't understand why we need `void_dt` and what `a.view(void_dt)` does (also, what doe an `np.void` data type do). Could you please elaborate a bit? — p-value, Jan 16 '18 at 05:17
@Hilbert `void_dt` is the datatype that we use to view each row as a scalar, basically packing more data into one element. Going the other way would be a case of unpacking, as shown in this example case - https://stackoverflow.com/a/48163550/3293881, as we unpack one element into three elements. — Divakar, Jan 16 '18 at 06:11

Efficiently compute pairwise equal for NumPy arrays

1 Answers1