7

I have a boolean matrix with 1.5E6 rows and 20E3 columns, similar to this example:

M = [[ True,  True, False,  True, ...],
     [False,  True,  True,  True, ...],
     [False, False, False, False, ...],
     [False,  True, False, False, ...],
     ...
     [ True,  True, False, False, ...]
     ]

Also, I have another matrix N ( 1.5E6 rows, 1 column):

 N = [[ True],
      [False],
      [ True],
      [ True],
      ...
      [ True]
      ]

What I need to do, is to go through each column pair from matrix M (1&1, 1&2, 1&3, 1&N, 2&1, 2&2 etc) combined by the AND operator, and count how many overlaps there are between the result and matrix N.

My Python/Numpy code would look like this:

for i in range(M.shape[1]):
  for j in range(M.shape[1]):
    result = M[:,i] & M[:,j] # Combine the columns with AND operator
    count = np.sum(result & N.ravel()) # Counts the True occurrences
    ... # Save the count for the i and j pair

The problem is, going through 20E3 x 20E3 combinations with two for loops is computationally expensive (takes around 5-10 days to compute). A better option I tried is comparing each column to the whole matrix M:

for i in range(M.shape[1]):
  result = M[:,i]*M.shape[1] & M # np.tile or np.repeat is used to horizontally repeat the column
  counts = np.sum(result & N*M.shape[1], axis=0)
  ... # Save the counts

This reduces overhead and calculation time to around 10%, but it's still taking 1 day or so to compute.

My question would be :
what is the fastest way (non Python maybe?) to make these calculations (basically just AND and SUM)?

I was thinking about low level languages, GPU processing, quantum computing etc.. but I don't know much about any of these so any advice regarding the direction is appreciated!

Additional thoughts: Currently thinking if there is a fast way using the dot product (as Davikar proposed) for computing triplets of combinations:

def compute(M, N):
    out = np.zeros((M.shape[1], M.shape[1], M.shape[1]), np.int32)
    for i in range(M.shape[1]):
        for j in range(M.shape[1]):
            for k in range(M.shape[1]):
                result = M[:, i] & M[:, j] & M[:, k]
                out[i, j, k] = np.sum(result & N.ravel())
    return out
4
  • 1
    Is there any reason for the tensorflow tag in the question? When you say "1.5mm", is that 1.5 million rows? Just to clarify. Also, what do you do with count? Are you accumulating the total sum, or each count is stored individually, or something else?
    – javidcf
    Commented Feb 4, 2020 at 13:49
  • None of the tags are too relevant since it's just a theoretical question. Let me know if there are better tags. 1.5 million rows, yes. I am storing the count of True occurrences, np.sum is returning the count for a bool array. Commented Feb 4, 2020 at 13:56
  • Yes but I meant, all of the count values should be summed into a global total, or the count for each pair of columns should be stored separately? (in your loops, you are just replacing the value of count in each iteration, so I'm not sure if you do something else with it later or just accumulate it).
    – javidcf
    Commented Feb 4, 2020 at 13:58
  • So i put dot dot dot, this is exactly what I'm doing after, basically storing the count per i and j pair. Commented Feb 4, 2020 at 13:59

3 Answers 3

9

Simply use np.einsum to get all the counts -

np.einsum('ij,ik,i->jk',M,M.astype(int),N.ravel())

Feel free to play around with optimize flag with np.einsum. Also, feel free to play around with different dtypes conversion.

To leverage GPU, we can use tensorflow package that also supports einsum.

Faster alternatives with np.dot :

(M&N).T.dot(M.astype(int))
(M&N).T.dot(M.astype(np.float32))

Timings -

In [110]: np.random.seed(0)
     ...: M = np.random.rand(500,300)>0.5
     ...: N = np.random.rand(500,1)>0.5

In [111]: %timeit np.einsum('ij,ik,i->jk',M,M.astype(int),N.ravel())
     ...: %timeit (M&N).T.dot(M.astype(int))
     ...: %timeit (M&N).T.dot(M.astype(np.float32))
227 ms ± 191 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
66.8 ms ± 198 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
3.26 ms ± 753 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

And take it a bit further with float32 conversions for both of the boolean arrays -

In [122]: %%timeit
     ...: p1 = (M&N).astype(np.float32)
     ...: p2 = M.astype(np.float32)
     ...: out = p1.T.dot(p2)
2.7 ms ± 34.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
9
  • This is very much appreciated! Seems to work in theory, but somehow counts are higher than expected: repl.it/repls/PlainUnwieldyCurrency Commented Feb 4, 2020 at 15:06
  • 1
    @FrancWeser When you say results are higher, do you mean the counts are larger than the ones that you are getting from your loop based solution(s)?
    – Divakar
    Commented Feb 4, 2020 at 15:08
  • 1
    @FrancWeser Well theoretically this should work. I have also tested with a random data of (500,300) and that worked too. Are you sure you have checked it properly. How are you checking again?
    – Divakar
    Commented Feb 4, 2020 at 15:24
  • 1
    @FrancWeser Yeah, just add one more arg for M in einsum : np.einsum('ij,ik,il,i->jkl',M,M,M.astype(int),N.ravel()).
    – Divakar
    Commented Feb 4, 2020 at 16:52
  • 1
    @FrancWeser Well sadly np.dot accepts two arguments only. So, we could work with two matrices as was the case before. With three versions of M as needed for the edited code, we don't have a good way. Think that einsum is the best you can do there.
    – Divakar
    Commented Feb 4, 2020 at 17:16
6

EDIT: To fix the code below to the fit the corrected question, just a couple of minor changes are required in compute:

def compute(m, n):
    m = np.asarray(m)
    n = np.asarray(n)
    # Apply mask N in advance
    m2 = m & n
    # Pack booleans into uint8 for more efficient bitwise operations
    # Also transpose for better caching (maybe?)
    mb = np.packbits(m2.T, axis=1)
    # Table with number of ones in each uint8
    num_bits = (np.arange(256)[:, np.newaxis] & (1 << np.arange(8))).astype(bool).sum(1)
    # Allocate output array
    out = np.zeros((m2.shape[1], m2.shape[1]), np.int32)
    # Do the counting with Numba
    _compute_nb(mb, num_bits, out)
    # Make output symmetric
    out = out + out.T
    # Add values in diagonal
    out[np.diag_indices_from(out)] = m2.sum(0)
    # Scale by number of ones in n
    return out

I would do this with Numba, using a few tricks. First, you can do only half of the column-wise operations, since the other half is repeated. Second, you can pack the boolean values into bytes so with each & you are operating over eight values instead of one. Third, you can use multiprocessing to parallelize it. In total, you could do it like this:

import numpy as np
import numba as nb

def compute(m, n):
    m = np.asarray(m)
    n = np.asarray(n)
    # Pack booleans into uint8 for more efficient bitwise operations
    # Also transpose for better caching (maybe?)
    mb = np.packbits(m.T, axis=1)
    # Table with number of ones in each uint8
    num_bits = (np.arange(256)[:, np.newaxis] & (1 << np.arange(8))).astype(bool).sum(1)
    # Allocate output array
    out = np.zeros((m.shape[1], m.shape[1]), np.int32)
    # Do the counting with Numba
    _compute_nb(mb, num_bits, out)
    # Make output symmetric
    out = out + out.T
    # Add values in diagonal
    out[np.diag_indices_from(out)] = m.sum(0)
    # Scale by number of ones in n
    out *= n.sum()
    return out

@nb.njit(parallel=True)
def _compute_nb(mb, num_bits, out):
    # Go through each pair of columns without repetitions
    for i in nb.prange(mb.shape[0] - 1):
        for j in nb.prange(1, mb.shape[0]):
            # Count common bits
            v = 0
            for k in range(mb.shape[1]):
                v += num_bits[mb[i, k] & mb[j, k]]
            out[i, j] = v

# Test
m = np.array([[ True,  True, False,  True],
              [False,  True,  True,  True],
              [False, False, False, False],
              [False,  True, False, False],
              [ True,  True, False, False]])
n = np.array([[ True],
              [False],
              [ True],
              [ True],
              [ True]])
out = compute(m, n)
print(out)
# [[ 8  8  0  4]
#  [ 8 16  4  8]
#  [ 0  4  4  4]
#  [ 4  8  4  8]]

As a quick comparison, here is a small benchmark against the original loop and NumPy-only methods (I am pretty sure the proposals by Divakar are the best you can get from NumPy):

import numpy as np

# Original loop

def compute_loop(m, n):
    out = np.zeros((m.shape[1], m.shape[1]), np.int32)
    for i in range(m.shape[1]):
        for j in range(m.shape[1]):
            result = m[:, i] & m[:, j]
            out[i, j] = np.sum(result & n)
    return out

# Divakar methods

def compute2(m, n):
    return np.einsum('ij,ik,lm->jk', m, m.astype(int), n)

def compute3(m, n):
    return np.einsum('ij,ik->jk',m, m.astype(int)) * n.sum()

def compute4(m, n):
    return np.tensordot(m, m.astype(int),axes=((0,0))) * n.sum()

def compute5(m, n):
    return m.T.dot(m.astype(int))*n.sum()

# Make random data
np.random.seed(0)
m = np.random.rand(1000, 100) > .5
n = np.random.rand(1000, 1) > .5
print(compute(m, n).shape)
# (100, 100)

%timeit compute(m, n)
# 768 µs ± 17.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit compute_loop(m, n)
# 11 s ± 1.23 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit compute2(m, n)
# 7.65 s ± 1.06 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit compute3(m, n)
# 23.5 ms ± 1.53 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit compute4(m, n)
# 8.96 ms ± 194 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit compute5(m, n)
# 8.35 ms ± 266 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
2
  • This is awesome, thanks, going to check it! As mentioned in comments to Divakar's post, I made the mistake of writing "np.sum(result & n)" instead of "np.sum(result & n.ravel())" which wrongly results in a 1D*2D array sum. Intended count is of 1D*1D array. I'll see if I can adjust your code for that. Commented Feb 4, 2020 at 15:40
  • 1
    @FrancWeser I suspected as much, as it seemed a bit of a strange operation. In any case, I edited the answer with the fix for that.
    – javidcf
    Commented Feb 4, 2020 at 15:47
0

I'd suggest trying to the get Python out of the way: convert your columns to bit fields, convert your N to a bit field as well, & together each triplet, then use (bin(num).count('1')) (or a proper popcnt if numpy has one).

1
  • Thanks, best language/resource to do this? Commented Feb 4, 2020 at 14:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.