def vectorized_approach(allData):
    M, N, R = allData.shape
valid_mask = np.mod(np.arange(M * M), M + 1) != 0
out = np.einsum('ijkl,ijkl->lij', allData[:, (None), :, :], allData[(None),
    :, :, :])
return out.reshape(R, -1)[:, (valid_mask)]
numpy cross-correlation - vectorizing