What is the fastest way to compare patches of an array?

I think dependant on what your some_metric function does, it's possible there is already a highly-optimized implementation out there. For example if it is just a convolution then take a look at Theano which will even let you leverage GPU quite easily. Even if your function isn't quite as simple as a straightforward convolution it is likely there'll be building blocks within Theano you can use rather than trying to go really low-level with Cython.

initialy posted another answer based on pattern matching (carried away by the title), post another answer

use one loop to loop over both arrays (big and small) and apply partial correlation metric (or whatever) without slicing lists all the time:

as a sidenote, observe the fact (depending on metric of course) that once one patch has been completed, the next patch (either one row down or one column right) shares much with the previous patch, only its initial and final rows (or columns) change thus the new correlation can be computed faster from the previous correlation by subtracting previous row and adding new row (or vice-versa). Since metric function is not given is added as observation.

def get_best_fit(big_array, small_array):    # we assume the small array is square    patch_size = small_array.shape[0]    x = 0     y = 0    x2 = 0     y2 = 0    W = big_array.shape[0]    H = big_array.shape[1]    W2 = patch_size    H2 = patch_size    min_value = np.inf    tmp = 0    min_patch = None    start = 0     end = H*W    start2 = 0    end2 = W2*H2    while start < end:        tmp = 0        start2 = 0        x2 = 0         y2 = 0        valid = True        while start2 < end2:            if x+x2 >= W or y+y2 >= H:                 valid = False                break            # !!compute partial metric!!            # no need to slice lists all the time            tmp += some_metric_partial(x, y, x2, y2, big_array, small_array)            x2 += 1             if x2>=W2:                 x2 = 0                 y2 += 1            start2 += 1        # one patch covered        # try next patch        if valid and min_value > tmp:            min_value = tmp            min_patch = [x, y, W2, H2]        x += 1         if x>=W:             x = 0             y += 1        start += 1    return min_patch, min_value

The other issue with your time measurement of the parallel function is that you call reshape on your array object after you run your parallel function. It could be the case that the parallel function is faster, but then reshape is adding extra time (although reshape should be quite fast, but I'm not sure how fast).

The other issue is we don't know what your some_metric term is, and it doesn't appear that you are using that in parallel. The way I see your parallel code working is that you are getting the patches in parallel and then putting them in a queue to be analyzed by the some_metric function, hence defeating the purpose of the parallelization of your code.

As John Greenhall said, using FFTs and convolutions can be quite fast. However, to take advantage of parallel processing, you would still need to do the analysis of the patch and small array in parallel as well.

How big are the arrays? Is memory an issue here as well?