Scipy curve_fit / leastsq становятся медленнее при задании якобиана?

#python #scipy #curve-fitting #least-squares

#python #scipy #подгонка кривой #метод наименьших квадратов

Вопрос:

Итак, я прочитал документацию о 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 (документ здесь), который принимает 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
 

А? Использование якобиана медленнее, чем его аппроксимация? Я сделал что-то не так?

Ответ №1:

На самом деле это не ответ, но я чувствую, что это зависит от размера проблемы. При небольшом размере (n = 500) дополнительное время, затрачиваемое на оценку якобиана (с предоставленным jac), вероятно, в конечном итоге не окупится.

n = 500, с jab:

 100 loops, best of 3: 1.50 ms per loop
 

Без:

 100 loops, best of 3: 1.57 ms per loop
 

n = 5000, с jab:

 100 loops, best of 3: 5.07 ms per loop
 

Без:

 100 loops, best of 3: 6.46 ms per loop
 

n = 50000, с jab:

 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