在scipy.optimize.leastsq中使用預分配的緩衝區內存函數擬合函數

Question

我正在嘗試使用scipy.optimize.leastsq，該函數使用預分配的內存來存儲殘差。我有數以百萬計的非線性擬合，時間非常關鍵。我在C中編寫了擬合函數，並且我注意到，當fit函數簡單地返回對作爲輸入獲取的緩衝區內存的引用時，scipy.optimize.leastsq無法正常工作。當我發現問題可以在純粹的Python中複製時，我以爲我搞砸了Py_INCREF。下面是一些說明問題的實體模型代碼（實際問題有不同的適合的功能和複雜得多）：

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) 

爲了使它正常工作，我必須包裝擬合函數爲一體，使一個副本：

# 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()在上面的例子？

Answer 1

好吧，我認爲這是在這裏發生的事情：底層FORTRAN例程LMDIF呈現用戶定義的函數與內存*fvec，其中結果將被存儲。該指針可能指向暫存內存，因爲LMDIF需要緩存幾個函數評估的結果來估計雅可比行列式。

由於用戶定義函數由Python調用，在*fvec存儲器不能直接使用，所以raw_multipack_lm_function()由第一和然後評估Python函數的結果應對到*fvec作品的包裝。在輸入LMDIF之前，用戶定義的函數被調用一次以找出輸出數組的形狀。

問題出現是因爲來自第一次功能評估的內存然後作爲原始*fvec傳遞到LMDIF，就好像它是新的一次性內存一樣。 LMDIF繼續使用它來存儲第一個功能評估，然後用不同的*fvec再次調用用戶功能。但在dy()的示例中，用戶函數在將結果複製到LMDIF需要的內存之前，覆蓋前一個函數調用的內存。這使得LMDIF認爲結果從不改變，因此退出的情況是擬合參數對擬合質量沒有影響。

我認爲這是一個錯誤的minpack_lmdif()從scipy/optimize/__minpack.h，因爲它假定用戶定義的函數總是返回新的一次性內存，這是不符合dy()的情況下（這似乎是編寫擬合函數完全合法的方式）。下面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（）的原始示例中一樣，這允許將適配函數編碼爲不受內存分配限制。由於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 

這是事。 mp.leastsq仍然評估一般的Python函數，限制條件是第一個參數將被結果覆蓋，第三個參數是緩衝區內存。但是，讓我們看看會發生什麼，如果我們在C代碼dy()：

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的完美矢量化代碼的性能（至少對於短陣列）來說非常重要。讓我們使用改進包裝：和mp.leastsq

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來自除暴安良元組建築規範和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爲你做計算！對於許多後續配合（我的應用程序）的進一步加速可以來自對LMDIF重新使用暫存存儲器wa。

最重要的是，所有計算現在都返​​回正確的結果！