给定雅可比行列式时，Scipy的curve_fit / leastsq会变慢吗？答案

【问题标题】：Scipy's curve_fit / leastsq become slower when given the Jacobian?给定雅可比行列式时，Scipy的curve_fit / leastsq会变慢吗？
【发布时间】：2014-08-14 23:42:22
【问题描述】：

所以我开始阅读有关 curve_fit here 的文档。它包含以下示例：

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 (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

我又重新计时了：

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

嗯？使用雅可比比近似它慢？我是不是做错了什么？

【问题讨论】：

标签： python scipy curve-fitting least-squares

【解决方案1】：

不是一个真正的答案，但我的感觉是它取决于问题的大小。对于小尺寸（n=500），花费在评估 jacobian（使用提供的 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

你也可以考虑重写 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

【讨论】：