I want calculate the matrix product of a boolean array with its transpose: import numpy as np a = np.array([[1, 0, 1], [1, 1, 0]], dtype=bool) What's the best/fastest way of doing this? I made several attemps: Attempt 1 out1 = np.matmul(a, a.T) print(out1) [[ True True] [ True True]] Just applying np.matmul doesn't work because of overflow. Attempt 2 So I came up with the idea of passing a np.uint64 output array to the function: Doesn't work either. Seems like the sum is calculated with bools, then casted to uint64. out2 = np.zeros((2, 2), dtype=np.uint64) np.matmul(a, a.T, out=out2) print(out2) [[1 1] [1 1]] Attempt 3 Manually calculating the row*column product then calling sum works: But it takes way too much memory, as it stores all product results before summing them up. out3 = (a[None,:, :] * a[:, None, :]).sum(axis=-1) print(out3) [[2 1] [1 2]] Attempt 4 So I tried casting the array before calculation: Does work too. But casting the input array takes 64 times of the memory and I guess that the multiplication takes more time if conducted with uint64 compared to two booleans. out4 = np.matmul(a.astype(np.uint64), a.T) print(out4) [[2 1] [1 2]] Is there any better solution for this problem that I might have overseen? Continue reading...