6

numpyを使って効率的に解きたい小さな線形方程式系がたくさんあります。基本的に、与えられたA[:,:,:]と、私は与えられたb[:,:]を見つけたいと思います。したがって、速度を気にしない場合は、次のように解決できます。x[:,:]A[i,:,:].dot(x[i,:]) = b[i,:]

for i in range(n):
    x[i,:] = np.linalg.solve(A[i,:,:],b[i,:])

しかし、これにはPythonでの明示的なループが含まれ、A通常はのような形をしているため(1000000,3,3)、このようなソリューションは非常に遅くなります。numpyがこれに達していない場合は、Fortranでこのループを実行できます（つまり、f2pyを使用）が、可能であればPythonのままにしておきたいと思います。

4

4 に答える 4

4

今この質問を読みに戻ってきた人たちのために、私は他の人の時間を節約し、numpyが今放送を使ってこれを処理していると言いたいと思いました。

したがって、numpy 1.8.0以降では、以下を使用してN個の一次方程式を解くことができます。

x = np.linalg.solve(A,b)
于 2015-12-01T18:42:44.203 に答える
3

自分自身に答えるのはちょっとしたことだと思いますが、これは私が今持っているFortranソリューションです。つまり、速度と簡潔さの両方で、他のソリューションが効果的に競合しているものです。

function pixsolve(A, b) result(x)
    implicit none
    real*8    :: A(:,:,:), b(:,:), x(size(b,1),size(b,2))
    integer*4 :: i, n, m, piv(size(b,1)), err
    n = size(A,3); m = size(A,1)
    x = b
    do i = 1, n
        call dgesv(m, 1, A(:,:,i), m, piv, x(:,i), m, err)
    end do
end function

これは次のようにコンパイルされます。

f2py -c -m foo{,.f90} -llapack -lblas

そしてPythonから次のように呼び出されます

x = foo.pixsolve(A.T, b.T).T

（.Tf2pyの設計上の選択が不十分なため、sが必要になります。これにより、不要なコピー、非効率的なメモリアクセスパターン、および.Tsが省略された場合のFortranインデックスの不自然な検索の両方が発生します。）

これにより、setup.pyなども回避されます。Fortranで選択するボーンはありません（文字列が含まれていない限り）が、numpyに同じことができる短くてエレガントなものがあることを期待していました。

于 2012-11-17T17:50:19.680 に答える
3

明示的なループが問題になるのは間違っていると思います。通常、それは最適化する価値のある最も内側のループにすぎません。これはここでも当てはまると思います。たとえば、オーバーヘッドのコードと実際の計算のコストを比較できます。

import numpy as np

n = 10**6
A = np.random.random(size=(n, 3, 3))
b = np.random.random(size=(n, 3))
x = b*0

def f():
    for i in xrange(n):
        x[i,:] = np.linalg.solve(A[i,:,:],b[i,:])

np.linalg.pseudosolve = lambda a,b: b

def g():
    for i in xrange(n):
        x[i,:] = np.linalg.pseudosolve(A[i,:,:],b[i,:])

それは私に

In [66]: time f()
CPU times: user 54.83 s, sys: 0.12 s, total: 54.94 s
Wall time: 55.62 s

In [67]: time g()
CPU times: user 5.37 s, sys: 0.01 s, total: 5.38 s
Wall time: 5.40 s

IOW、実際に問題を解決する以外のことに費やしている時間は10％にすぎません。さて、Fortranから得られるものに比べて、それ自体が遅すぎると完全に信じることnp.linalg.solveができたので、あなたは何か他のことをしたいと思っています。これはおそらく小さな問題に特に当てはまります。考えてみてください。IIRC私はかつて、特定の小さなソリューションを手作業で展開する方が速いことがわかりましたが、それはしばらく前のことです。

しかし、それ自体では、最初のインデックスで明示的なループを使用すると、ソリューション全体が非常に遅くなるというのは事実ではありません。十分に速い場合np.linalg.solve、ループはここではあまり変更しません。

于 2012-11-18T15:21:01.187 に答える
0

対角線の周りの3x3ブロックで構成される（3x100000,3x100000）行列を使用して、一度に実行できると思います。

未検証 ：

b_new = np.vstack([ b[i,:] for i in range(len(i)) ])
x_new = np.zeros(shape=(3x10000,3) )

A_new = np.zeros(shape=(3x10000,3x10000) )
n,m = A.shape
for i in range(n):
   A_new[3*i:3*(i+1),3*i:3*(i+1)] = A[i,:,:]

x = np.linalg.solve(A_new,b_new)
于 2012-11-17T16:02:55.913 に答える