在 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 吃掉了。但是，通过将普通 python3 与 map 和 mul:

结合使用，我可以获得比第一个示例适度的加速

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 条件不计算大小写 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 循环。这在数学上似乎是合理的，但它给出的结果略有不同，我猜这是由于精度问题。我不确定如何解决这些问题，但我想我会 post 它以防您可以忍受精度问题或有人知道如何解决这些问题。

然而，考虑到我使用的 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)

更新二：

相同，除了乘法较少和 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 的加速，对于我的 'real data' 个 N~6000 数组运行 125s 我现在得到 0.0005s，这是大约 2.5*10^5 的惊人加速。非常感谢，IVlad 和 hpaulj，在最后一天学到了很多东西:) P.S。我很惊讶你们回答的速度如此之快，我花了半天时间才跟进；)

在 python3 中计算双和内的点积的有效方法

Efficient way of computing dot product inside double sum in python3

algorithm

python-3.x

numpy

sum

dot-product