1

我正在研究如何在 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 吃掉了。但是，通过使用带有 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 条件不计算大小写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's weave 、numba、pyrex 或 python 或者完全走 C 路径……</p>

提前感谢您的帮助！

编辑：

这就是数据样本的样子：

# 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()
4

3 回答 3

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])
于 2015-02-11T23:27:32.553 回答
1

那个双循环在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))
于 2015-02-12T01:00:57.320 回答
0

好的，非常感谢您的帮助。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，这是一个约 0 的加速，而对于我现在运行2*10^4的“真实数据”数组，这是一个惊人的加速约。非常感谢，IVlad 和 hpaulj，在最后一天学到了很多东西 :) PS 我很惊讶你们回答的速度之快让我花了半天时间才跟进；）N~6000125s0.0005s2.5*10^5

于 2015-02-12T20:19:24.477 回答