Scipy的curve_fit/leastsq提供雅可比矩阵时变慢？

Question

我在阅读关于 curve_fit 的文档，具体内容可以在 这里 找到。文档里有一个例子：

import numpy as np
import scipy.optimize as so

def func(x, a,b,c ):
    return a * np.exp(-b * x) + c

a,b,c = 2.5, 1.3, 0.5
nx = 500
noiseAlpha = 0.5
#
xdata = np.linspace(0, 4, nx)
y = func(xdata, a,b,c)
ydata = y + noiseAlpha * np.random.normal(size=len(xdata))

如果我现在调用 curve_fit，它会因为我没有提供任何东西而自己估算导数。我们来测一下时间：

%%timeit
popt, pcov = so.curve_fit(func, xdata, ydata )
1000 loops, best of 3: 787 µs per loop

实际上，curve_fit 会调用 leastsq，你可以在 这里 查看文档，它接受一个 Dfun 参数来计算雅可比矩阵。所以我这样做了：

def myDfun( abc, xdata, ydata, f ) :
    a,b,c = abc
    ebx = np.exp(-b * xdata)
    res = np.vstack( ( ebx, a * -xdata * ebx, np.ones(len(xdata)) ) ).T
    return res

然后我又测了一次时间：

%%timeit
popt, pcov = so.curve_fit(func, xdata, ydata, Dfun=myDfun )
1000 loops, best of 3: 858 µs per loop

嗯？使用雅可比矩阵的速度比估算的还慢？我是不是做错什么了？

Answer 1

这不是一个真正的答案，但我觉得这要看问题的大小。对于小规模的问题（比如n=500），花额外的时间去计算雅可比矩阵（用提供的jac）可能最后得不到相应的回报。

当n=500，使用雅可比矩阵：

100 loops, best of 3: 1.50 ms per loop

不使用：

100 loops, best of 3: 1.57 ms per loop

当n=5000，使用雅可比矩阵：

100 loops, best of 3: 5.07 ms per loop

不使用：

100 loops, best of 3: 6.46 ms per loop

当n=50000，使用雅可比矩阵：

100 loops, best of 3: 49.1 ms per loop

不使用：

100 loops, best of 3: 59.2 ms per loop

另外，你可能还想考虑重写雅可比函数，比如去掉耗时的.T()，这样可以提高大约15%的速度：

def myDfun2( abc, xdata, ydata, f ) :
    a,b,c = abc
    ebx = np.exp(-b * xdata)
    res = np.hstack( ( ebx.reshape(-1,1), (a * -xdata * ebx).reshape(-1,1), np.ones((len(xdata),1)) ) )
    return res