операции Оптимизировать с массивами в NumPy

Question

я должен применить некоторую математическую формулу, которую я написал в питона, как:

for s in range(tdim): 
     sum1 = 0.0 
     for i in range(dim): 
      for j in range(dim): 
       sum1+=0.5*np.cos(theta[s]*(i-j))* 
       eig1[i]*eig1[j]+eig2[i]+eig2[j])-0.5*np.sin(theta[s]*(i-j))*eig1[j]*eig2[i]-eig1[i]*eig2[j]) 

     PHi2.append(sum1) 

Теперь, это правильно, но явно неэффективна, наоборот, чтобы сделать :

for i in range(dim): 
      for j in range(dim): 
       PHi2 = 0.5*np.cos(theta*(i-j))*(eig1[i]*eig1[j]+eig2[i]+eig2[j])-0.5*np.sin(theta*(i-j))*(eig1[j]*eig2[i]-eig1[i]*eig2[j]) 

Однако второй пример дает мне такое же число во всех элементах phi2, так что это быстрее, но ответ неверен. Как вы можете сделать это правильно и более эффективно?

ПРИМЕЧАНИЕ: eig1 и eig2 имеют одинаковую размерность d, тета и PHi2 имеют одинаковую размерность D, BUT d! = D.

Answer 1

В первом разряде кода у вас есть 0.5*np.cos(theta[s]*(i-j))..., а во втором - 0.5*np.cos(theta*(i-j)).... Если у вас нет тета, определенной по-разному для второго бита кода, это может быть причиной проблемы.

Answer 2

Это работает вещанием.
Для tdim = 200 and dim = 100.
14 секунд с оригиналом.
120 мс с версией.

s = np.array(range(tdim))[:, None, None] 
i = np.array(range(dim))[:, None] 
j = np.array(range(dim)) 
PHi2 =(0.5*np.cos(theta[s]*(i-j))*(eig1[i]*eig1[j]+eig2[i]+eig2[j])-0.5*np.sin(theta[s]*(i-j))*(eig1[j]*eig2[i]-eig1[i]*eig2[j])).sum(axis=2).sum(axis=1) 

Answer 3

Вы можете использовать грубую силу подхода вещания, но вы создаете промежуточный массив формы (D, d, d), которая может выйти из под контролем, если ваши массивы даже умеренно большие. Кроме того, при использовании трансляции без каких-либо уточнений вы повторно вычисляете множество вычислений из самого внутреннего цикла, который вам нужно сделать только один раз. Если вы сначала вычислить необходимые параметры для всех возможных значений i - j и добавить их вместе, вы можете использовать эти значения на внешнем контуре, например:

def fast_ops(eig1, eig2, theta): 
    d = len(eig1) 
    d_arr = np.arange(d) 
    i_j = d_arr[:, None] - d_arr[None, :] 
    reidx = i_j + d - 1 
    mult1 = eig1[:, None] * eig1[ None, :] + eig2[:, None] + eig2[None, :] 
    mult2 = eig1[None, :] * eig2[:, None] - eig1[:, None] * eig2[None, :] 
    mult1_reidx = np.bincount(reidx.ravel(), weights=mult1.ravel()) 
    mult2_reidx = np.bincount(reidx.ravel(), weights=mult2.ravel()) 

    angles = theta[:, None] * np.arange(1 - d, d) 

    return 0.5 * (np.einsum('ij,j->i', np.cos(angles), mult1_reidx) - 
        np.einsum('ij,j->i', np.sin(angles), mult2_reidx)) 

, если переписать код M4rtini как функция для сравнения:

def fast_ops1(eig1, eig2, theta): 
    d = len(eig1) 
    D = len(theta) 
    s = np.array(range(D))[:, None, None] 
    i = np.array(range(d))[:, None] 
    j = np.array(range(d)) 
    ret = 0.5 * (np.cos(theta[s]*(i-j))*(eig1[i]*eig1[j]+eig2[i]+eig2[j]) - 
       np.sin(theta[s]*(i-j))*(eig1[j]*eig2[i]-eig1[i]*eig2[j])) 
    return ret.sum(axis=(-1, -2)) 

И мы составляем некоторые данные:

d, D = 100, 200 
eig1 = np.random.rand(d) 
eig2 = np.random.rand(d) 
theta = np.random.rand(D) 

улучшения скорости очень заметно, 80x на вершине 115x над исходным кодом, что приводит к ускоренному 9000 разному ускорению:

In [22]: np.allclose(fast_ops1(eig1, eig2, theta), fast_ops(eig1, eig2, theta)) 
Out[22]: True 

In [23]: %timeit fast_ops1(eig1, eig2, theta) 
10 loops, best of 3: 145 ms per loop 

In [24]: %timeit fast_ops(eig1, eig2, theta) 
1000 loops, best of 3: 1.85 ms per loop