Neighbors sum of numpy array with mask

Question

I have two large arrays, one containing values, and one being a mask basically. The code below shows the function I want to implement.

from scipy.signal import convolve2d
import numpy as np

sample = np.array([[6, 4, 5, 5, 5],
                   [7, 1, 0, 8, 3],
                   [2, 5, 4, 8, 4],
                   [2, 0, 2, 6, 0],
                   [5, 7, 2, 3, 2]])

mask = np.array([[1, 0, 1, 1, 0],
                 [0, 0, 1, 0, 1],
                 [0, 1, 0, 0, 0],
                 [0, 0, 0, 1, 0],
                 [1, 1, 0, 0, 1]])

neighbors_sum = convolve2d(sample, np.ones((3,3), dtype=int), mode='same', boundary='wrap')

# neighbors_sum = np.array([[40, 37, 35, 33, 44],
#                           [37, 34, 40, 42, 48],
#                           [24, 23, 34, 35, 40],
#                           [27, 29, 37, 31, 32],
#                           [31, 33, 34, 30, 34]])

result = np.where(mask, neighbors_sum, 0)

print(result)

This code works, and gets me what I expects:

np.array([[40,  0, 35, 33,  0],
          [ 0,  0, 40,  0, 48],
          [ 0, 23,  0,  0,  0],
          [ 0,  0,  0, 31,  0],
          [31, 33,  0,  0, 34]])

So far, so good. However, where I'm encountering some large issue is when I increase the size of the arrays. In my case, instead of a 5x5 input and a 3x3 summing mask, I need a 50,000x20,000 input and a 100x100 summing mask. And when I move to that, the convolve2d function is in all kinds of trouble and the calculation is extremely long.

Given that I only care about the masked result, and thus only care about the summation from convolve2d at those points, can anyone think of a smart approach to take here? Going to a for loop and selecting only the points of interest would lose the speed advantage of the vectorization so I'm not convinced this would be worth it.

Any suggestion welcome!