Compute operations of all rows in matrix with itself efficiently

Question

Let say I have a 2D matrix M of some arbitrary size (n, m). Now I want to do some vector operations of all rows in M with all other rows in M efficiently. In my case, I want to compute a bitwise xnor between every row.

What we do know is that there's no need to compute xnor(x, y) and xnor(y, x), it's enough to compute one of them. Also, in my case, there's no need to compute xnor(x, x). I know two solutions to this problem:

def xnor(p, q):
    return abs(p-1) * abs(q-1) + (p * q)

def solution_one(M):
    return xor(M[:, None], M).all(axis=2)

def solution_two(M):
    n = M.shape[0]
    results = np.zeros((n, n))

    # Upper triangular indices of 
    ii, jj = np.triu_indices(n, 1)
    for i, j in zip(ii,jj):
        results[i,j] = xnor(M[i], M[j]).all(axis=0)

    return results

Running the two with some arbitrary large matrices we get for solution one

# 100x100
%timeit solution_one(M)
# 3.55 ms ± 73.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# 500x500
%timeit solution_one(M)
# 1.36 s ± 59.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

# 1000x1000
%timeit solution_one(M)
# 27.4 s ± 603 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

and for solution two

# 100x100
%timeit solution_two(M)
# 29.5 ms ± 535 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

# 500x500
%timeit solution_two(M)
# 1.06 s ± 41.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

# 1000x1000
%timeit solution_two(M)
# 5.58 s ± 95.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

As one see, for smaller instances solution_one is faster than solution_two, but for bigger it is the opposite.

We can see that solution_one is doing n x n calculations whereas solution_two will do (as mentioned by Naphat Amundsen in the comments) ((n x n) - n) / 2. Is there better solutions to this? Can one somehow construct a new matrix of M such that solution one will do closer to ((n x n) - n) / 2 calculations?

user_13091658 · Answer

I believe the lowest number of computations you can do in this case is namely ((n x n) - n) / 2. If that is the case, then in terms of complexity your solution_two is already optimal. However, I believe can you make solution_two more vectorized, which will probably improve runtime performance. Here is my attempt:

import numpy as np

def solution_two(M):
    n = M.shape[0]
    results = np.zeros((n, n))

    # Upper triangular indices of
    ii, jj = np.triu_indices(n, 1)
    for i, j in zip(ii,jj):
        results[i,j] = xnor(M[i], M[j]).all(axis=0)

    return results

def solution_two_vectorized(M):
    n = M.shape[0]
    results = np.zeros((n, n))
    ii, jj = np.triu_indices(n, 1)
    results[ii,jj] = xnor(M[ii], M[jj]).all(-1)
    return results

# Quick test checking that your outputs are same as mine on some random matrices
import time
# Generator of random test matrices
Ms = (np.random.randint(-10,10,(np.random.randint(100,200),np.random.randint(100,200))) for i in
      range(100))

vectortimes = []
loopytimes = []

for M in Ms:
    t0 = time.time()
    out_vectorized = solution_two_vectorized(M)
    vectortimes.append(time.time() - t0)

    t0 = time.time()
    out_loopy = solution_two(M)
    loopytimes.append(time.time() - t0)

    if not np.allclose(out_vectorized, out_loopy):
        print("Output mismatch!")
        break
else:
    print("Success!")
# Success!

mean_vectortimes = np.mean(vectortimes)
mean_loopytimes = np.mean(loopytimes)

print("Vector times: ", mean_vectortimes)
# Vector times:  0.017428524494171142
print("Loopy times: ", mean_loopytimes)
# Loopy times:  0.07040918111801148
print(f"Performance improvement: {mean_loopytimes/mean_vectortimes} times")
# Performance improvement: 4.039881926984661 times

Compute operations of all rows in matrix with itself efficiently

1 Answers