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 条件を使用して casej=kを計算しないようにします。

r_np[j] - r_np[k] = 0

したがって、ドット積も0になるか、合計を2つに分割して同じ結果を達成します

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問題は、数千の合計を計算する必要があるため、トリプルのセットに対して 1 桁の秒ランタイムが必要なことです。それ以外の場合は、scipy の 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数秒になりました。

これで十分だと思います！

アップデート：

これは、単一の for ループで約 10.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])
于 2015-02-11T23:27:32.553 に答える
1

その二重ループは、 のタイム キラーですnumpy。ベクトル化された配列操作を使用すると、評価は 1 秒未満に短縮されます。

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との交換20% オフnp.einsum

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 に答える