2

scipy.optimize.leastsq事前に割り当てられたメモリを使用して残差を格納するフィット関数を使用しようとしています。何百万もの非線形フィットがあり、時間が重要です。で fit 関数をコーディングしましたが、fit 関数が、入力として取得されたバッファ メモリへの参照を返すだけの場合、正しく動作しないCことに気付きました。scipy.optimize.leastsq問題がPy_INCREF純粋なPython. この問題を示すモックアップ コードを次に示します (実際の問題にはさまざまな適合関数があり、はるかに複雑です)。

import numpy as np
import scipy.optimize as opt

# the mock-up data to be fitted
x = np.arange(5.)
noise = np.array([0.0, -0.02, 0.01, -0.03, 0.01])
y = np.exp(-x) + noise

# preallocate the buffer memory
buf = np.empty_like(x)

# fit function writing residuals to preallocated memory and returning a reference 
def dy(p, x, buf, y):
    buf[:] = p[0]*np.exp(-p[1]*x) - y
    return buf

# starting values (true values are [1., 1.])
p0 = np.array([1.2, 0.8])

# leastsq with dy(): DOESN'T WORK CORRECTLY
opt.leastsq(dy, p0, args=(x, buf, y))
# -> (array([ 1.2,  0.8]), 4)

正しく動作させるには、コピーを作成する関数に fit 関数をラップする必要があります。

# wrapper that calls dy() and returns a copy of the array
def dy2(*args): 
    return dy(*args).copy()

# leastsq with dy2(): WORKS CORRECTLY
opt.leastsq(dy2, p0, args=(x, buf, y))
# -> (array([ 0.99917134,  1.03603201]), 1)

...しかし、そもそもコピーを作成することは明らかにバッファ メモリの使用に不利です! 補足として、次のopt.fminようにバッファメモリでも動作します（ただし、実際には私のアプリケーションには遅すぎます）：

def sum2(p, x, buf, y):
    dy(p, x, buf, y)
    return buf.dot(buf)

opt.fmin(sum2, p0, args=(x, buf, y))
# Optimization terminated successfully.
#     Current function value: 0.001200
#     Iterations: 32
#     Function evaluations: 63
# -> array([ 0.99915812,  1.03600273])

上記の例で正しく動作し、scipy.optimize.leastsq正しく動作しdy2()ない理由はありますか?dy()

4

1 に答える 1

2

わかりました、これがここで起こっていることだと思います: 基礎となるFORTRANルーチンは、結果が格納されるLMDIFメモリ内のユーザー定義関数を提示します。このポインターは、ヤコビアンを推定するためにいくつかの関数評価の結果をキャッシュする必要があるため*fvec、スクラッチ メモリを指す場合があります。LMDIF

ユーザー定義関数は から呼び出されるPythonため、 のメモリを*fvec直接使用することはできません。そのため、ラッパーは、最初に関数raw_multipack_lm_function()を評価してから結果を にコピーすることで機能します。に入る前に、ユーザー定義関数が 1 回呼び出されて、出力配列の形状が検出されます。Python*fvecLMDIF

最初の関数評価からのメモリが、あたかも新しい使い捨てメモリであるかのようにLMDIF、元の として渡されるため、問題が発生します。これを使用して最初の関数評価を に格納し、別の. しかし、 の例では、ユーザー関数は、結果が必要なメモリにコピーされる前に、前の関数呼び出しからメモリを上書きします。これは、結果が決して変化しないと思い込ませるため、適合パラメータが適合の質に影響を与えないというケースで終了します。*fvecLMDIF*fvecdy()LMDIFLMDIF

ユーザー定義関数が常に新しい使い捨てメモリを返すと仮定しているため、これはminpack_lmdif()fromのバグであると考えていますが、これは当てはまりません(fit 関数をコーディングする完全に正当な方法のようです)。以下は、簡単な修正方法を示しています。scipy/optimize/__minpack.hdy()git diff

diff --git a/scipy/optimize/__minpack.h b/scipy/optimize/__minpack.h
index 2c0ea33..465724b 100644
--- a/scipy/optimize/__minpack.h
+++ b/scipy/optimize/__minpack.h
@@ -483,7 +483,7 @@ static PyObject *minpack_lmdif(PyObject *dummy, PyObject *args) {
   qtf = (double *) ap_qtf->data;
   fjac = (double *) ap_fjac->data;
   ldfjac = dims[1];
-  wa = (double *)malloc((3*n + m)* sizeof(double));
+  wa = (double *)malloc((3*n + 2*m)* sizeof(double));
   if (wa == NULL) {
     PyErr_NoMemory();
     goto fail;
@@ -492,12 +492,15 @@ static PyObject *minpack_lmdif(PyObject *dummy, PyObject *args) {

   /* Call the underlying FORTRAN routines. */
   n_int = n; /* to provide int*-pointed storage for int argument of LMDIF */
-  LMDIF(raw_multipack_lm_function, &m, &n_int, x, fvec, &ftol, &xtol, &gtol, &maxfev, &epsfcn, diag,
-    
+  LMDIF(raw_multipack_lm_function, &m, &n_int, x, wa+3*n+m, &ftol, &xtol, &gtol, &maxfev, &epsfcn, d
+
   RESTORE_FUNC();

   if (info < 0) goto fail;           /* Python error */

+  /* Copy final residuals back to initial array */
+  memcpy(fvec, wa+3*n+m, m*sizeof(double));
+
   free(wa);
   Py_DECREF(extra_args); 
   Py_DECREF(ap_diag);

mそのため、より多くのスクラッチ スペースの要素を に割り当てLMDIF、追加のメモリを初期の として使用します*fvec。これにより、ユーザー関数が呼び出されたときのメモリ衝突が回避されます。正しい最終残差を返すにmemcpy()は、元の配列に最終結果を格納する追加のものが必要です。

dy() を使用した元の例のように、これにより fit 関数をコーディングしてメモリ割り当てから解放することができます。内側のループ全体LMDIFがメモリ割り当てから解放されるため、パフォーマンスの向上が期待できます。

アップデート

ここにいくつかのタイミング結果があります。明らかに、テストの問題は非常に小さく、収束が速いため、実際のアプリケーションを代表するものではない可能性があります。これは、パッチを適用したバージョンのscipy.optimize.leastsq:

In [1]: def dy0(p, x, y): return p[0]*np.exp(-p[1]*x) - y

In [2]: %timeit p = opt.leastsq(dy2, p0, args=(x, buf, y))
1000 loops, best of 3: 399 us per loop

In [3]: %timeit p = opt.leastsq(dy, p0, args=(x, buf, y))
1000 loops, best of 3: 363 us per loop

In [4]: %timeit p = opt.leastsq(dy0, p0, args=(x, y))
1000 loops, best of 3: 341 us per loop

したがって、事前に割り当てられたメモリに書き込むことによって得られるものは何もありません。 の直接的な実装がdy0()最も高速です。LMDIFしかし、事前に割り当てられたメモリをより有効に活用するためのより効率的なラッパーを作成するとどうなるでしょうか? ここに私が書いたものがあります：

In [5]: %timeit p = mp.leastsq(dy, (p0.copy(), x, buf, y))
1000 loops, best of 3: 279 us per loop

それは何かです。1 番目の引数が結果で上書きされ、3 番目の引数がバッファ メモリであるという制限付きで、mp.leastsq一般的な関数を評価します。しかし、次のようにコーディングするとPythonどうなるか見てみましょう:dy()C

In [6]: %timeit p = opt.leastsq(fitfun.e2_diff, p0, args=(x, buf, y))
10000 loops, best of 3: 48.2 us per loop

痛い！numpy完全にベクトル化されたコードの のパフォーマンス (少なくとも短い配列の場合) については、これで終わりです。改善されたラッパーを使用しましょう。

In [7]: %timeit p = mp.leastsq(fitfun.e2_diff, (p0.copy(), x, buf, y))
100000 loops, best of 3: 6.94 us per loop

opt.leastsqとの間の追加の高速化はmp.leastsq、タプル構築コードと を取り除くことによるものmemcpyです。最後に、LMDIFコールバックなしの生のパフォーマンスPython:

In [8]: %timeit p = fitfun.e2_fit(p0.copy(), x, buf, y)
100000 loops, best of 3: 6.13 us per loop

それほど違いはありません。そのため、Python へのコールバックはあまりコストがかからないように見えますnumpyが、計算を任せてはいけません! その後の多くのフィット (私のアプリケーション) のさらなる高速化はwa、LMDIF.

最も重要なことは、すべての計算が正しい結果を返すようになったことです!

于 2013-01-20T14:55:45.567 に答える