在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 是一个包含数千个三元组的数组，而 x 是一个常数三元组。长度为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 吃掉了。 但是，通过使用带有 map 和 mul 的普通 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

尝试使用 if 条件不计算 case 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 的 weave、numba、pyrex 或 python，或者完全走 C 路径……

提前感谢您的帮助！

编辑：

这就是数据样本的样子：

# 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 秒。

然后我摆脱了zips 和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 秒。

我希望这已经足够好了！

更新：

这里是 0.01 秒左右，带有一个 for 循环。它在数学上似乎是合理的，但它给出的结果略有不同，我猜这是由于精度问题。我不知道如何解决这些问题，但我想我会发布它，以防你可以忍受精度问题或者有人知道如何解决它们。

但考虑到我使用的 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

好的，非常感谢您的帮助。 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 附近运行。在我第一次尝试使用 ~4s 的非矢量化 numpy 时，这是一个大约 2*10^4 的加速，对于我的 N~6000 的“真实数据”数组，它运行大约 125s 我现在得到 0.0005s，它是大约2.5*10^5 的惊人加速。非常感谢，IVlad 和 hpaulj，在最后一天学到了很多东西 :) 附言我很惊讶你们回答的速度之快让我花了半天的时间才跟进；）