распараллелить цикл над iter

Во-первых, давайте посмотрим на простую векторизацию вашей проблемы. У меня такое ощущение, что вы хотите, чтобы ваши array_a и array_b были точно такими же, то есть координаты частиц, но я держу их здесь отдельно.

Я превратил ваш код в функцию, чтобы упростить синхронизацию:

def IIII(array_a, array_b, bins) :
    hist_data=np.zeros([len(bins)-1])
    for i in product(array_a,array_b):
        (result,bins) = np.histogram(newdist(i), bins=bins)
        hist_data+=result
    hist_data = np.array(map(float, hist_data))
    return hist_data

Кстати, вы можете сгенерировать свои образцы данных менее запутанным способом следующим образом:

n = 100
array_a = np.random.uniform(0, 10, size=(n, 2))
array_b = np.random.uniform(0, 10, size=(n, 2))

Итак, сначала нам нужно векторизовать ваш func. Я сделал так, что он может принимать любую array формы (4, ...). Для экономии памяти он делает расчет на месте и возвращает первую плоскость, т.е. array[0].

def func_vectorized(a) :
    a[1] **= 22
    np.sin(a[2], out=a[2])
    a[0] *= a[1]
    a[0] += a[2]
    a[0] += a[3]
    return a[0]

Имея эту функцию, мы можем написать векторизованную версию IIII:

def IIII_vec(array_a, array_b, bins) :
    array_ab = np.empty((4, len(array_a), len(array_b)))
    a = np.swapaxes(array_a, 0 ,1)
    b = np.swapaxes(array_b, 0 ,1)
    array_ab[[0, 1]] = a[:, :, None]
    array_ab[[2, 3]] = b[:, None, :]
    newdist = func_vectorized(array_ab)
    hist, _ = np.histogram(newdist, bins=bins)
    return hist

С n = 100 точками они оба возвращают одно и то же:

In [2]: h1 = IIII(array_a, array_b, x_edges)

In [3]: h2 = IIII_bis(array_a, array_b, x_edges)

In [4]: np.testing.assert_almost_equal(h1, h2)

А вот временные различия уже очень актуальны:

In [5]: %timeit IIII(array_a, array_b, x_edges)
1 loops, best of 3: 654 ms per loop

In [6]: %timeit IIII_vec(array_a, array_b, x_edges)
100 loops, best of 3: 2.08 ms per loop

Ускорение в 300 раз!. Если вы попробуете еще раз с более длинными выборочными данными, n = 1000, вы увидите, что они оба одинаково плохо масштабируются, как n**2, поэтому 300x остается там:

In [10]: %timeit IIII(array_a, array_b, x_edges)
1 loops, best of 3: 68.2 s per loop

In [11]: %timeit IIII_bis(array_a, array_b, x_edges)
1 loops, best of 3: 229 ms per loop

Итак, вы все еще смотрите на хорошие 10 мин. обработки, что на самом деле не так много по сравнению с более чем 2 днями, которые потребуются вашему текущему решению.

Конечно, чтобы все было так хорошо, вам нужно поместить в память массив (4, 50000, 50000) чисел с плавающей запятой, что моя система не может обработать. Но вы по-прежнему можете делать вещи относительно быстро, обрабатывая их по частям. Следующая версия IIII_vec делит каждый массив на d фрагментов. Как написано, длина массива должна быть кратна d. Было бы несложно преодолеть это ограничение, но оно запутало бы истинную цель:

def IIII_vec_bis(array_a, array_b, bins, d=1) :
    a = np.swapaxes(array_a, 0 ,1)
    b = np.swapaxes(array_b, 0 ,1)
    a_chunk = len(array_a) // d
    b_chunk = len(array_b) // d
    array_ab = np.empty((4, a_chunk, b_chunk))
    hist_data = np.zeros((len(bins) - 1,))
    for j in xrange(d) :
        for k in xrange(d) :
            array_ab[[0, 1]] = a[:, a_chunk * j:a_chunk * (j + 1), None]
            array_ab[[2, 3]] = b[:, None, b_chunk * k:b_chunk * (k + 1)]
            newdist = func_vectorized(array_ab)
            hist, _ = np.histogram(newdist, bins=bins)
            hist_data += hist
    return hist_data

Во-первых, давайте проверим, что это действительно работает:

In [4]: h1 = IIII_vec(array_a, array_b, x_edges)

In [5]: h2 = IIII_vec_bis(array_a, array_b, x_edges, d=10)

In [6]: np.testing.assert_almost_equal(h1, h2)

А теперь немного таймингов. С n = 100:

In [7]: %timeit IIII_vec(array_a, array_b, x_edges)
100 loops, best of 3: 2.02 ms per loop

In [8]: %timeit IIII_vec_bis(array_a, array_b, x_edges, d=10)
100 loops, best of 3: 12 ms per loop

Но по мере того, как вы начинаете иметь все больший и больший массив в памяти, выполнение его по частям начинает окупаться. С n = 1000:

In [12]: %timeit IIII_vec(array_a, array_b, x_edges)
1 loops, best of 3: 223 ms per loop

In [13]: %timeit IIII_vec_bis(array_a, array_b, x_edges, d=10)
1 loops, best of 3: 208 ms per loop

С n = 10000 я больше не могу вызывать IIII_vec без ошибки array is too big, но краткая версия все еще работает:

In [18]: %timeit IIII_vec_bis(array_a, array_b, x_edges, d=10)
1 loops, best of 3: 21.8 s per loop

И просто чтобы показать, что это можно сделать, я запустил его один раз с n = 50000:

In [23]: %timeit -n1 -r1 IIII_vec_bis(array_a, array_b, x_edges, d=50)
1 loops, best of 1: 543 s per loop

Итак, хорошие 9 минут обработки чисел, что не так уж и плохо, учитывая, что было вычислено 2,5 миллиарда взаимодействий.

Используйте векторизованные операции numpy. Замените цикл for поверх product() одним вызовом newdist(), создав аргументы с помощью meshgrid().

Чтобы парализовать задачу, вычислите newdist() на срезах array_a, array_b, которые соответствуют подблокам meshgrid(). Вот пример использования срезов и многопроцессорной обработки.

Вот еще один пример для демонстрации шагов: цикл python -> векторизованная версия numpy -> parallel:

#!/usr/bin/env python
from __future__ import division
import math
import multiprocessing as mp
import numpy as np

try:
    from itertools import izip as zip
except ImportError:
    zip = zip # Python 3

def pi_loop(x, y, npoints):
    """Compute pi using Monte-Carlo method."""
    #  note: the method converges to pi very slowly.
    return 4 * sum(1 for xx, yy in zip(x, y) if (xx**2 + yy**2) < 1) / npoints

def pi_vectorized(x, y, npoints):
    return 4 * ((x**2 + y**2) < 1).sum() / npoints # or just .mean()

def mp_init(x_shared, y_shared):
    global mp_x, mp_y
    mp_x, mp_y = map(np.frombuffer, [x_shared, y_shared]) # no copy

def mp_pi(args):
    # perform computations on slices of mp_x, mp_y
    start, end = args
    x = mp_x[start:end] # no copy
    y = mp_y[start:end]
    return ((x**2 + y**2) < 1).sum()

def pi_parallel(x, y, npoints):
    # compute pi using multiple processes
    pool = mp.Pool(initializer=mp_init, initargs=[x, y])
    step = 100000
    slices = ((start, start + step) for start in range(0, npoints, step))
    return 4 * sum(pool.imap_unordered(mp_pi, slices)) / npoints

def main():
    npoints = 1000000

    # create shared arrays
    x_sh, y_sh = [mp.RawArray('d', npoints) for _ in range(2)]

    # initialize arrays
    x, y = map(np.frombuffer, [x_sh, y_sh])
    x[:] = np.random.uniform(size=npoints)
    y[:] = np.random.uniform(size=npoints)

    for f, a, b in [(pi_loop, x, y), 
                    (pi_vectorized, x, y), 
                    (pi_parallel, x_sh, y_sh)]:
        pi = f(a, b, npoints)
        precision = int(math.floor(math.log10(npoints)) / 2 - 1 + 0.5)
        print("%.*f %.1e" % (precision + 1, pi, abs(pi - math.pi)))

if __name__=="__main__":
    main()

Показатели времени для npoints = 10_000_000:

pi_loop pi_vectorized pi_parallel 
   32.6         0.159       0.069 # seconds

Он показывает, что основное преимущество в производительности заключается в преобразовании цикла python в его векторизованный аналог numpy.