2

我正在尝试使用scipy.optimize.leastsq使用预分配内存来存储残差的 fit 函数。我有数百万次非线性拟合，时间很关键。我在 中编写了 fit 函数C，我注意到scipy.optimize.leastsq当 fit 函数只是返回对作为输入获得的缓冲存储器的引用时，它不能正常工作。Py_INCREF当我注意到问题可以在 pure 中重现时，我以为我搞砸了s 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提供了带有内存的用户定义函数 in *fvec，结果将存储在其中。这个指针可能指向暂存内存，因为LMDIF需要缓存几个函数评估的结果来估计雅可比。

由于从 调用用户定义的函数，不能直接使用Pythonin 的内存，因此包装器的工作方式是先评估函数，然后将结果处理到。在进入之前，会调用一次用户自定义函数来找出输出数组的形状。*fvecraw_multipack_lm_function()Python*fvecLMDIF

出现问题是因为第一次函数评估的内存然后LMDIF作为原始传递给*fvec，就好像它是新的一次性内存一样。LMDIF继续使用它来存储第一个函数评估，然后再次使用不同的*fvec. 但是在带有 的示例中dy()，用户函数在将结果复制到所需的内存之前，会覆盖上一次函数调用的内存LMDIF。这愚蠢LMDIF地认为结果永远不会改变，因此退出拟合参数对拟合质量没有影响的情况。

我认为这是minpack_lmdif()from中的一个错误scipy/optimize/__minpack.h，因为它假定用户定义的函数总是返回新的一次性内存，但情况并非如此dy()（这似乎是编写 fit 函数的一种完全合法的方式）。下面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

那是东西。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来自摆脱元组构建代码和memcpys。最后，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 回答