给予雅可比人时，Scipy的curve_fit / leastsq会变慢？

So I wad reading the documentation about curve_fit here. It contains the following example:

所以我在这里阅读有关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))

If I call curve_fit now, it will approximate the derivatives since I didn't provide anything. Let's time it:

如果我现在调用curve_fit，它将接近导数，因为我没有提供任何东西。我们来说吧：

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

In fact curve_fit calls leastsq (doc here), which accepts a Dfun argument for computing the Jacobian. So I did this:

事实上，curve_fit调用leastsq（doc here），它接受用于计算雅可比行列式的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

And I timed it again:

我再次计时：

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

Uh? Using the Jacobian is slower than approximating it? Did i do something wrong?

呃？使用雅可比行列式比接近它更慢？我做错什么了吗？

1 个解决方案

#1

Not really an answer, but my feeling is it dependents on the size of the problem. For small size (n=500), the extra time spend in evaluating jacobian (with the supplied jac) probably don't pay off in the end.

不是一个真正的答案，但我的感觉是它依赖于问题的大小。对于小尺寸（n = 500），在评估雅可比（使用提供的jac）时花费的额外时间可能最终没有得到回报。

n=500, with jab:

n = 500，用jab：

100 loops, best of 3: 1.50 ms per loop

Without:

无：

100 loops, best of 3: 1.57 ms per loop

n=5000, with jab:

n = 5000，用刺戳：

100 loops, best of 3: 5.07 ms per loop

Without:

无：

100 loops, best of 3: 6.46 ms per loop

n=50000, with jab:

n = 50000，用jab：

100 loops, best of 3: 49.1 ms per loop

Without:

无：

100 loops, best of 3: 59.2 ms per loop

Also you may want to consider rewrite the jacobian function, for example getting rid of the expensive .T() can bring about ~15% speed up:

你也可以考虑重写jacobian函数，例如摆脱昂贵的.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