7

私が解決したい実際の問題は、N個の単位ベクトルのセットと別のM個のベクトルのセットが与えられた場合、各単位ベクトルについて、 M個のベクトルのすべてとの内積の絶対値の平均を計算することです。基本的に、これは 2 つの行列の外積を計算し、絶対値を間に挟んで合計および平均化します。

NとMが大きすぎない場合、これは難しくなく、続行するには多くの方法があります (以下を参照)。問題は、NとMが大きい場合、作成される一時ファイルが巨大になり、提供されたアプローチに実際的な制限が生じることです。この計算は、一時的に作成せずに実行できますか? 私が抱えている主な困難は、絶対値の存在によるものです。そのような計算を「スレッド化」するための一般的な手法はありますか?

例として、次のコードを考えてみましょう

N = 7
M = 5

# Create the unit vectors, just so we have some examples,
# this is not meant to be elegant
phi = np.random.rand(N)*2*np.pi
ctheta = np.random.rand(N)*2 - 1
stheta = np.sqrt(1-ctheta**2)
nhat = np.array([stheta*np.cos(phi), stheta*np.sin(phi), ctheta]).T

# Create the other vectors
m = np.random.rand(M,3)

# Calculate the quantity we desire, here using broadcasting.
S = np.average(np.abs(np.sum(nhat*m[:,np.newaxis,:], axis=-1)), axis=0)

これは素晴らしいことです。S は長さNの配列になり、目的の結果が含まれています。残念ながら、その過程で巨大な可能性のあるアレイをいくつか作成してしまいました。結果として

np.sum(nhat*m[:,np.newaxis,:], axis=-1)

はM X N配列です。もちろん、最終結果はサイズNのみです。NとMのサイズを増やし始めると、すぐにメモリ エラーが発生します。

上記のように、絶対値が必要ない場合は、次のように進めることができます。einsum()

T = np.einsum('ik,jk,j', nhat, m, np.ones(M)) / M

これは、非常に大きなNおよびMに対しても機能し、迅速に機能します。特定の問題については、を含める必要がありますabs()が、より一般的な解決策 (おそらくより一般的な ufunc) も興味深いでしょう。

4

3 に答える 3

3

いくつかのコメントに基づいて、cython を使用するのが最善の方法のようです。私はばかげて cython の使用を検討したことがありません。動作するコードを作成するのは比較的簡単であることがわかりました。

いくつか検索した後、次の cython コードをまとめました。これは最も一般的なコードではなく、おそらく最適な記述方法ではなく、おそらくより効率的にすることができます。それでもeinsum()、元の質問のコードよりも約 25% 遅いだけなので、それほど悪くはありません! 元の質問で行われたように作成された配列で明示的に機能するように書かれています（したがって、入力配列の想定モード）。
警告にもかかわらず、元の問題に対して合理的に効率的な解決策を提供し、同様の状況での出発点として役立ちます。

import numpy as np
cimport numpy as np
import cython
DTYPE = np.float64
ctypedef np.float64_t DTYPE_t
cdef inline double d_abs (double a) : return a if a >= 0 else -a

@cython.boundscheck(False)
@cython.wraparound(False)
def process_vectors (np.ndarray[DTYPE_t, ndim=2, mode="fortran"] nhat not None,
                     np.ndarray[DTYPE_t, ndim=2, mode="c"] m not None) :
    if nhat.shape[1] != m.shape[1] :
        raise ValueError ("Arrays must contain vectors of the same dimension")
    cdef Py_ssize_t imax = nhat.shape[0]
    cdef Py_ssize_t jmax = m.shape[0]
    cdef Py_ssize_t kmax = nhat.shape[1] # same as m.shape[1]
    cdef np.ndarray[DTYPE_t, ndim=1] S = np.zeros(imax, dtype=DTYPE)
    cdef Py_ssize_t i, j, k
    cdef DTYPE_t val, tmp
    for i in range(imax) :
        val = 0
        for j in range(jmax) :
            tmp = 0
            for k in range(kmax) :
                tmp += nhat[i,k] * m[j,k]
            val += d_abs(tmp)
        S[i] = val / jmax
    return S
于 2013-07-14T03:02:43.537 に答える
1

正確な操作を高速化する簡単な方法（Cythonなど以外）はないと思います。しかし、計算しているものを本当に計算する必要があるかどうかを検討したい場合があります。絶対値の平均の代わりに二乗平均平方根を使用できる場合でも、何らかの方法で内積の大きさを平均化できますが、次のように 1 回のショットで取得できます。

rms = np.sqrt(np.einsum('ij,il,kj,kl,k->i', nhat, nhat, m, m, np.ones(M)/M))

これは、次のことと同じです。

rms_2 = np.sqrt(np.average(np.einsum('ij,kj->ik', nhat, m)**2, axis=-1))

はい、それはまさにあなたが求めていたものではありませんが、ベクトル化されたアプローチで得られるものと同じくらい近いと思います. この道を進むことにした場合は、大規模なandnp.einsumのパフォーマンスを確認してください。あまりにも多くのパラメーターとインデックスを渡すと、動きが鈍くなる傾向があります。NM

于 2013-07-13T05:44:57.113 に答える
0

これはかなり遅くなりますが、大きな中間行列は作成されません。

vals = np.zeros(N)
for i in xrange(N):
    u = nhat[i]
    for v in m:
        vals[i]+=abs(np.dot(u,v))
    vals[i]=vals[i]/M

編集: for ループの外で M で除算するように移動しました。

edit2: 新しいアイデア。後世と関連するコメントのために古いアイデアを保持します。

m2 = np.average(m,0)
vals = np.zeros(N)
for i in xrange(N):
    u=nhat[i]
    vals[i]=abs(np.dot(u,m2))

これは高速ですが、異なる値が返されることがあります。その理由については現在取り組んでいますが、当面は役立つかもしれません。

編集 3: ああ、それは絶対値のことです。うーん

>>> S
array([ 0.28620962,  0.65337876,  0.37470707,  0.46500913,  0.49579837,
        0.29348924,  0.27444208,  0.74586928,  0.35789315,  0.3079964 ,
        0.298353  ,  0.42571445,  0.32535728,  0.87505053,  0.25547394,
        0.23964505,  0.44773271,  0.25235646,  0.4722281 ,  0.33003338])
>>> vals
array([ 0.2099343 ,  0.6532155 ,  0.33039334,  0.45366889,  0.48921527,
        0.20467291,  0.16585856,  0.74586928,  0.31234917,  0.22198642,
        0.21013519,  0.41422894,  0.26020981,  0.87505053,  0.1199069 ,
        0.06542492,  0.44145805,  0.08455833,  0.46824704,  0.28483342])

time to compute S: 0.000342130661011 seconds
time to compute vals: 7.29560852051e-05 seconds

編集 4: 単位ベクトルのほとんどが正の値である場合、m のベクトルがダミー データのように常に正であると仮定すると、これはより速く実行されるはずです。

m2 = np.average(m,0)
vals = np.zeros(N)
for i in xrange(N):
    u=nhat[i]
    if u[0] >= 0 and u[1] >= 0 and u[2] >= 0:
        vals[i] = abs(np.dot(u,m2))
    else:
        for j in xrange(M):
            vals[i]+=abs(np.dot(u,m[j]))
        vals[i]/=M
于 2013-07-12T21:39:32.063 に答える