11

私はPython/numpy / scipyを使用して、小さなレイトレーサーを作成しています。サーフェスは、法平面より上の高さを与える2次元関数としてモデル化されます。光線とサーフェスの交点を見つける問題を、1つの変数を持つ関数のルートを見つける問題に減らしました。機能は継続的であり、継続的に微分可能です。

scipyルートファインダーを使用して（そしておそらく複数のプロセスを使用して）、単にすべての関数をループするよりも効率的にこれを行う方法はありますか？

編集：関数は、光線を表す線形関数と、交点の平面に制約された表面関数との違いです。

4

2 に答える 2

4

次の例は、二分法を使用して、関数x **（a + 1）-b（すべて異なるaとb）の100万コピーの根を並列に計算する方法を示しています。ここでは約12秒かかります。

import numpy

def F(x, a, b):
    return numpy.power(x, a+1.0) - b

N = 1000000

a = numpy.random.rand(N)
b = numpy.random.rand(N)

x0 = numpy.zeros(N)
x1 = numpy.ones(N) * 1000.0

max_step = 100
for step in range(max_step):
    x_mid = (x0 + x1)/2.0
    F0 = F(x0, a, b)
    F1 = F(x1, a, b)
    F_mid = F(x_mid, a, b)
    x0 = numpy.where( numpy.sign(F_mid) == numpy.sign(F0), x_mid, x0 )
    x1 = numpy.where( numpy.sign(F_mid) == numpy.sign(F1), x_mid, x1 )
    error_max = numpy.amax(numpy.abs(x1 - x0))
    print "step=%d error max=%f" % (step, error_max)
    if error_max < 1e-6: break

基本的な考え方は、変数のベクトルと個々のコンポーネント関数を定義するパラメーターの同等のベクトルで評価できる関数を使用して、変数のベクトルでルートファインダーの通常のすべてのステップを並列に実行することです。条件文は、マスクとnumpy.where（）の組み合わせに置き換えられます。これは、すべてのルートが必要な精度で見つかるまで、または問題からそれらを削除してそれらのルートを除外する小さな問題を続行する価値がある十分なルートが見つかるまで続けることができます。

私が解決するために選択した関数は任意ですが、関数が正常に動作している場合に役立ちます。この場合、ファミリ内のすべての関数は単調であり、正の根が1つだけあります。さらに、二分法では、関数のさまざまな符号を与える変数を推測する必要があります。これらは、ここでも非常に簡単に思い付くことができます（x0とx1の初期値）。

上記のコードはおそらく最も単純な求根アルゴリズム（二分法）を使用していますが、同じ手法をニュートンラプソン法やリダー法などに簡単に適用できます。求根法の条件が少ないほど、これに適しています。ただし、必要なアルゴリズムを再実装する必要があります。既存のライブラリルートファインダー関数を直接使用する方法はありません。

上記のコードスニペットは、速度ではなく明確さを念頭に置いて書かれています。一部の計算の繰り返しを回避し、特に関数を3回ではなく、反復ごとに1回だけ評価すると、次のように最大9秒高速化されます。

...
F0 = F(x0, a, b)
F1 = F(x1, a, b)

max_step = 100
for step in range(max_step):
    x_mid = (x0 + x1)/2.0
    F_mid = F(x_mid, a, b)
    mask0 = numpy.sign(F_mid) == numpy.sign(F0)
    mask1 = numpy.sign(F_mid) == numpy.sign(F1)
    x0 = numpy.where( mask0, x_mid, x0 )
    x1 = numpy.where( mask1, x_mid, x1 )
    F0 = numpy.where( mask0, F_mid, F0 )
    F1 = numpy.where( mask1, F_mid, F1 )
...

比較のために、scipy.bisect（）を使用して一度に1つのルートを見つけるには、約94秒かかります。

for i in range(N):
    x_root = scipy.optimize.bisect(lambda x: F(x, a[i], b[i]), x0[i], x1[i], xtol=1e-6)
于 2012-11-27T03:16:37.217 に答える
1

過去数年のいつか、scipy.optimize.newtonベクトル化のサポートを得ました。他の回答の例を使用すると、次のようになります。

import numpy as np
from scipy import optimize

def F(x, a, b):
    return np.power(x, a+1.0) - b

N = 1000000

a = np.random.rand(N)
b = np.random.rand(N)

optimize.newton(F, np.zeros(N), args=(a, b))

これは、他の回答のベクトル化された二分法と同じくらい速く実行されます。

于 2021-11-24T06:48:44.383 に答える