Почему scipy.optimize.curve_fit повторно оценивает исходное предположение (и, вероятно, дорого)?

При использовании curve_fit функция модели неоднократно и бесполезно (и, вероятно, дорого) оценивается без изменения параметров. Почему это происходит?

Рассмотрим следующий пример для определения параметров квадратичной функции:

from scipy.optimize import curve_fit

i=0

def f(x, a, b):
    global i; i += 1
    print('run: {:2}, p: {:<11.10}, {:<11.10}'.format(i, a, b))
    return(x**a+b)

popt, pcov = curve_fit(f, [0, 1, 2, 3], [1, 2, 5, 10], p0 = [1.5, 0.5])
print('\npopt:', popt)

Выход:

run:  1, p: 1.5        , 0.5        
run:  2, p: 1.5        , 0.5        
run:  3, p: 1.5        , 0.5        
run:  4, p: 1.500000022, 0.5        
run:  5, p: 1.5        , 0.5000000075
run:  6, p: 2.166073038, 0.9668143807
run:  7, p: 2.166073071, 0.9668143807
run:  8, p: 2.166073038, 0.9668143951
run:  9, p: 2.014374939, 0.9956632744
run: 10, p: 2.014374969, 0.9956632744
run: 11, p: 2.014374939, 0.9956632892
run: 12, p: 2.000113621, 0.9999686478
run: 13, p: 2.000113651, 0.9999686478
run: 14, p: 2.000113621, 0.9999686627
run: 15, p: 2.000000007, 0.999999998
run: 16, p: 2.000000037, 0.999999998
run: 17, p: 2.000000007, 1.000000013
run: 18, p: 2.0        , 1.0        

popt: [2. 1.]

Первая оценка вычисляет значения для начального предположения, но затем это делается снова во втором и третьем прогоне. Только в четвертой и пятой оценке оптимизация начинается с вычисления производных. Если оценки функций являются дорогостоящими и допускается некоторый допуск в более крупном масштабе, эти избыточные оценки могут занимать значительное количество времени.


python scipy scipy-optimize
person theBridge    schedule 28.01.2020    source источник

Ответы (1)


arrow_upward
1
arrow_downward

С full_output (см. leastsq документы):

In [125]: curve_fit(f, [0, 1, 2, 3], [1, 2, 5, 10], p0 = [1.5, 0.5],full_output=True)            
run:  1, p: 1.5        , 0.5        
run:  2, p: 1.5        , 0.5        
run:  3, p: 1.5        , 0.5        
run:  4, p: 1.500000022, 0.5        
run:  5, p: 1.5        , 0.5000000075
run:  6, p: 2.166073038, 0.9668143807
run:  7, p: 2.166073071, 0.9668143807
run:  8, p: 2.166073038, 0.9668143951
run:  9, p: 2.014374939, 0.9956632744
run: 10, p: 2.014374969, 0.9956632744
run: 11, p: 2.014374939, 0.9956632892
run: 12, p: 2.000113621, 0.9999686478
run: 13, p: 2.000113651, 0.9999686478
run: 14, p: 2.000113621, 0.9999686627
run: 15, p: 2.000000007, 0.999999998
run: 16, p: 2.000000037, 0.999999998
run: 17, p: 2.000000007, 1.000000013
run: 18, p: 2.0        , 1.0        
Out[125]: 
(array([2., 1.]), array([[ 0., -0.],
        [-0.,  0.]]), {'fvec': array([0., 0., 0., 0.]),
  'nfev': 16,
  'fjac': array([[-10.2688908 ,   0.        ,   0.26999885,   0.96286064],
         [ -1.2328595 ,  -1.5748198 ,   0.2521752 ,  -0.73019944]]),
  'ipvt': array([1, 2], dtype=int32),
  'qtf': array([-7.08708997e-08,  3.07498129e-09])}, 'The relative error between two consecutive iterates is at most 0.000000', 2)

Где возвращенная информация:

infodict : dict
a dictionary of optional outputs with the keys:

``nfev``
    The number of function calls
``fvec``
    The function evaluated at the output
``fjac``
    A permutation of the R matrix of a QR
    factorization of the final approximate
    Jacobian matrix, stored column wise.
    Together with ipvt, the covariance of the
    estimate can be approximated.
``ipvt``
    An integer array of length N which defines
    a permutation matrix, p, such that
    fjac*p = q*r, where r is upper triangular
    with diagonal elements of nonincreasing
    magnitude. Column j of p is column ipvt(j)
    of the identity matrix.
``qtf``
    The vector (transpose(q) * fvec).

Таким образом, он утверждает, что оценивает 16 раз, а не ваши 18. Таким образом, одна или несколько первоначальных оценок могут быть получены из проверки параметров curve_fit (или что-то в этом роде).

С другими методами:

In [134]: i=0                                                                                    
In [135]: curve_fit(f, [0, 1, 2, 3], [1, 2, 5, 10], p0 = [1.5, 0.5],method='trf')                
run:  1, p: 1.5        , 0.5        
run:  2, p: 1.500000022, 0.5        
run:  3, p: 1.5        , 0.5000000149
run:  4, p: 2.166073038, 0.9668143807
run:  5, p: 2.166073071, 0.9668143807
run:  6, p: 2.166073038, 0.9668143956
run:  7, p: 2.01437494 , 0.9956632758
run:  8, p: 2.01437497 , 0.9956632758
run:  9, p: 2.01437494 , 0.9956632907
run: 10, p: 2.000113621, 0.9999686478
run: 11, p: 2.000113651, 0.9999686478
run: 12, p: 2.000113621, 0.9999686627
run: 13, p: 2.000000007, 0.999999998
run: 14, p: 2.000000037, 0.999999998
run: 15, p: 2.000000007, 1.000000013
run: 16, p: 2.0        , 1.0        
run: 17, p: 2.00000003 , 1.0        
run: 18, p: 2.0        , 1.000000015
Out[135]: 
(array([2., 1.]), array([[ 3.77052335e-34, -1.19338002e-33],
        [-1.19338002e-33,  9.94005329e-33]]))

а также

In [136]: i=0                                                                                    
In [137]: curve_fit(f, [0, 1, 2, 3], [1, 2, 5, 10], p0 = [1.5, 0.5],method='dogbox')             
run:  1, p: 1.5        , 0.5        
run:  2, p: 1.500000022, 0.5        
run:  3, p: 1.5        , 0.5000000149
run:  4, p: 2.166073038, 0.9668143807
run:  5, p: 2.166073071, 0.9668143807
run:  6, p: 2.166073038, 0.9668143956
run:  7, p: 2.01437494 , 0.9956632758
run:  8, p: 2.01437497 , 0.9956632758
run:  9, p: 2.01437494 , 0.9956632907
run: 10, p: 2.000113621, 0.9999686478
run: 11, p: 2.000113651, 0.9999686478
run: 12, p: 2.000113621, 0.9999686627
run: 13, p: 2.000000007, 0.999999998
run: 14, p: 2.000000037, 0.999999998
run: 15, p: 2.000000007, 1.000000013
run: 16, p: 2.0        , 1.0        
run: 17, p: 2.00000003 , 1.0        
run: 18, p: 2.0        , 1.000000015
Out[137]: 
(array([2., 1.]), array([[ 3.77052335e-34, -1.19338002e-33],
        [-1.19338002e-33,  9.94005329e-33]]))
person hpaulj    schedule 28.01.2020
comment
Спасибо за вашу помощь! Чтобы разобраться в проблеме, я попытался напрямую использовать leastsq: хотя он утверждает (через nfev), что функция оценивается только 16 раз, на самом деле функция по-прежнему оценивается 18 раз. Мне также интересно, почему curve_fit использует устаревшую функцию leastsq, а не более новую least_squares. Поскольку на мой вопрос на самом деле нет ответа, я думаю, мне придется отредактировать свой вопрос? - person theBridge; 31.01.2020
comment
nfev происходит от retval, созданного вызовом _minpack. Я вижу в leastsq как минимум один вызов _check_func (вызывает функцию для проверки формы и типа) до этого. - person hpaulj; 31.01.2020
comment
Интересно, что least_squares не проверяет функцию, а сразу вычисляет производные. Однако в целом требует такого же количества вызовов, как leastsq (хотя я пытался прописать ту же точность). - person theBridge; 02.02.2020