Эффективный способ вычисления скалярного произведения в двойной сумме в Python3

Question

Я ищу в том, чтобы вычислить как можно более эффективные в Python3 продукте точки внутри двойной суммы формы:

import cmath 
for j in range(0,N): 
    for k in range(0,N): 
     sum_p += cmath.exp(-1j * sum(a*b for a,b in zip(x, [l - m for l, m in zip(r_p[j], r_p[k])]))) 

где r_np является массив из нескольких тысяч троек, а xa - постоянный тройной. Сроки для длины N=1000 троек составляет около 2.4s. То же самое с помощью NumPy:

import numpy as np 
for j in range(0,N): 
    for k in range(0,N): 
     sum_np = np.add(sum_np, np.exp(-1j * np.inner(x_np,(r_np[j] - r_np[k])))) 

на самом деле медленнее с выполнения около 4.0s. Я предполагаю, что это связано с отсутствием большого преимущества векторизации, только короткая 3 точка 3 является np.dot, которая съедается путем запуска N^2 из циклов. Однако скромное убыстрение по первому примеру, я мог бы получить с помощью простого python3 с картой и мулом:

from operator import mul 
for j in range(0,N): 
    for k in range(0,N): 
     sum_p += cmath.exp(-1j * sum(map(mul,x, [l - m for l, m in zip(r_p[j], r_p[k])]))) 

с выполнением около 2.0s

Попытки использовать либо если условие не рассчитал случай j=k, где

r_np[j] - r_np[k] = 0 

и, таким образом, скалярное произведение также становится равным 0, или разделив сумму в два, чтобы достигнуть того же самого

for j in range(0,N): 
     for k in range(j+1,N): 
    ... 
for k in range(0,N): 
     for j in range(k+1,N): 
    ... 

оба сделали это еще медленнее. Итак, все это масштабируется с помощью O (N^2), и мне интересно, можно ли с помощью некоторых методов сортировки или других вещей избавиться от циклов и сделать их масштабируемыми с помощью O (N logN). Проблема в том, что мне нужно однократное второе время автономной работы для набора из N~6000 троек, поскольку у меня есть тысячи этих сумм для вычисления. В противном случае мне нужно попробовать плетение scipy, numba, pyrex или python или полностью спуститься по пути C ...

Заранее благодарим за любую помощь!

Edit:

это как образец данных будет выглядеть следующим образом:

# numpy arrays 
x_np = np.array([0,0,1], dtype=np.float64) 
N=1000 
xy = np.multiply(np.subtract(np.random.rand(N,2),0.5),8) 
z = np.linspace(0,40,N).reshape(N,1) 
r_np = np.hstack((xy,z)) 

# in python format 
x = (0,0,1) 
r_p = r_np.tolist() 

Answer 1

Я использовал это, чтобы генерировать тестовые данные:

x = (1, 2, 3) 
r_p = [(i, j, k) for i in range(10) for j in range(10) for k in range(10)] 

На моей машине это произошло 2.7 секунд с вашим алгоритмом.

Тогда я избавилась от zip с и sum:

for j in range(0,N): 
    for k in range(0,N): 
     s = 0 
     for t in range(3): 
      s += x[t] * (r_p[j][t] - r_p[k][t]) 
     sum_p += cmath.exp(-1j * s) 

Это привело его до 2.4 секунд.

Тогда я заметил, что x постоянен так:

x * (p - q) = x1*p1 - x1*q1 + x2*p2 - x2*q2 - ... 

Так что я изменил код поколения к:

x = (1, 2, 3) 
r_p = [(x[0] * i, x[1] * j, x[2] * k) for i in range(10) for j in range(10) for k in range(10)] 

И алгоритм:

for j in range(0,N): 
    for k in range(0,N): 
     s = 0 
     for t in range(3): 
      s += r_p[j][t] - r_p[k][t] 
     sum_p += cmath.exp(-1j * s) 

Какие у меня до 2.0 секунд.

Тогда я понял, что мы можем переписать как:

for j in range(0,N): 
    for k in range(0,N): 
     sum_p += cmath.exp(-1j * (sum(r_p[j]) - sum(r_p[k]))) 

Что удивительно, у меня 1.1 секунд, которые я не могу объяснить - возможно, некоторые кэширование происходит?

В любом случае, кэширование или нет, вы можете предварительно скопировать суммы ваших троек, а затем вам не придется полагаться на механизм кеширования. Я сделал это:

sums = [sum(a) for a in r_p] 

sum_p = 0 
N = len(r_p) 
start = time.clock() 
for j in range(0,N): 
    for k in range(0,N): 
     sum_p += cmath.exp(-1j * (sums[j] - sums[k])) 

Который доставил меня до 0.73 секунд.

Надеюсь, это будет достаточно!

Update:

Вот один около 0.01 секунд с одной цикл. Это кажется математически здоровым, но он дает несколько разные результаты, которые, как я предполагаю, обусловлены проблемами точности. Я не уверен, как это исправить, но я думал, что опубликую его, если вы сможете жить с проблемами точности или кто-то знает, как их исправить.

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

sums = [sum(a) for a in r_p] 
e_denom = sum([cmath.exp(1j * p) for p in sums]) 
sum_p = 0 
N = len(r_p) 
start = time.clock() 
for j in range(0,N): 
    sum_p += e_denom * cmath.exp(-1j * sums[j]) 

print(sum_p) 
end = time.clock() 
print(end - start) 

Обновление 2:

То же самое, за исключением того, с меньшим количеством умножений и вызовите функцию sum:

sum_p = e_denom * sum([np.exp(-1j * p) for p in sums]) 

Answer 2

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

In [1764]: sum_np=0 

In [1765]: for j in range(0,N): 
    for k in range(0,N): 
     sum_np += np.exp(-1j * np.inner(x_np,(r_np[j] - r_np[k]))) 

In [1766]: sum_np 
Out[1766]: (2116.3316526447466-1.0796252780664872e-11j) 

In [1767]: np.exp(-1j * np.inner(x_np, (r_np[:N,None,:]-r_np[None,:N,:]))).sum((0,1)) 
Out[1767]: (2116.3316526447466-1.0796252780664872e-11j) 

Тайминги:

In [1768]: timeit np.exp(-1j * np.inner(x_np, (r_np[:N,None,:]-r_np[None,:N,:]))).sum((0,1)) 
1 loops, best of 3: 506 ms per loop 

In [1769]: %%timeit 
sum_np=0 
for j in range(0,N): 
    for k in range(0,N): 
     sum_np += np.exp(-1j * np.inner(x_np,(r_np[j] - r_np[k]))) 
1 loops, best of 3: 12.9 s per loop 

замещающие np.inner с np.einsum бреет 20% от времени

np.exp(-1j * np.einsum('k,ijk', x_np, r_np[:N,None,:]-r_np[None,:N,:])).sum((0,1)) 

Answer 3

Ok ребята, спасибо большое за помощь. Последнее значение имеет код IVlads, который использует идентификатор sum_j sum_k a[j]*a[k] = sum_j a[j] * sum_k a[k]. Это теперь масштабируется также с меньшим, чем O (N^2). предварительное вычисления скалярного произведения, прежде чем сумма составляет hpaulj в NumPy предложение точно так же быстро:

sum_np = 0 
dotprods = np.inner(q_np,r_np) 
sum_rkexp = np.exp(1j * dotprods).sum() 
sum_np = sum_rkexp * np.exp(-1j * dotprods).sum() 

как с выполнением около 0.0003s. Тем не менее, я нашел еще одну вещь, которая дает еще ~ 50% -ное увеличение, вместо вычисления экспоненты дважды, я беру комплексное сопряжение в сумме:

sum_np = 0 
dotprods = np.inner(q_np,r_np) 
rkexp = np.exp(1j * dotprods) 
sum_rkexp = rkexp.sum() 
sum_np = sum_rkexp * np.conj(rkexp).sum() 

который проходит около 0.0002s. За моими первые попытки с не векторизованной NumPy, которые имели ~4s, это убыстрение о 2*10^4, и для массивов моих реальных данных "и N~6000, которые бегают 125s я теперь получить 0.0005s, который является удивительным убыстрением о 2.5*10^5. Большое спасибо, IVlad и hpaulj, многому научились в последний день :) P.S. Я поражен тем, как быстро вы, ребята, отвечаете с вещами, которые потребовали мне полдня, чтобы просто следить;)