lmfit.minimize with method='lbfgsb' doesn't obey setting max_nfev > 15000

[rowlesmr]

First Time Issue Code

Yes, I read the instructions and I am sure this is a GitHub Issue.

Description

When minimising a residual using lmfit.minimize with method='lbfgsb', setting max_nfev to a value greater than the default value of maxiter (15000) of scipy.optimize.fmin_l_bfgs_b results in the minimisation still terminating at nfev=15000. Setting the keyword argument maxiter=50000 gives RuntimeWarning: ignoring 'maxiter' argument to 'Minimizer()'. Use 'max_nfev' instead..

A Minimal, Complete, and Verifiable example

import lmfit
import numpy as np


def G(a, b, x) -> float:
    return a * np.exp(-b * x * x)


def residual(params, x, data):
    a_s = [params['a1'], params['a2'], params['a3'], params['a4'], params['a5']]
    b_s = [params['b1'], params['b2'], params['b3'], params['b4'], params['b5']]
    c = params['c']

    model = c * np.ones(x.size)
    for a, b in zip(a_s, b_s):
        model += G(a, b, x)

    return data - model


X = np.asarray([0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.22, 0.24, 0.25, 0.26, 0.28, 0.30, 0.32, 0.34, 0.35, 0.36, 0.38, 0.40, 0.42, 0.44, 0.45, 0.46, 0.48, 0.50, 0.55, 0.60, 0.65, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.50, 3.00, 3.50, 4.00, 5.00, 6.00])
Y = np.asarray([57.000, 56.926, 56.708, 56.360, 55.900, 55.351, 54.736, 54.076, 53.388, 52.687, 51.982, 51.278, 50.580, 49.888, 49.202, 48.523, 47.849, 47.182, 46.519, 45.862, 45.212, 43.932, 42.686, 42.078, 41.481, 40.321, 39.212, 38.153, 37.145, 36.659, 36.185, 35.270, 34.397, 33.562, 32.760, 32.370, 31.988, 31.243, 30.523, 28.817, 27.231, 25.759, 24.401, 22.031, 20.106, 18.561, 17.300, 16.227, 15.265, 14.362, 13.489, 12.636, 11.807, 11.009, 10.253, 9.550, 6.917, 5.550, 4.820, 4.270, 3.240, 2.410])

params = lmfit.Parameters()

params.add('a1', value=11.4)
params.add('a2', value=21.66)
params.add('a3', value=14.82)
params.add('a4', value=6.84)
params.add('a5', value=1.71)

params.add('b1', value=0.2379567949423352)
params.add('b2', value=1.5687819528065554)
params.add('b3', value=16.957819441754367)
params.add('b4', value=37.086721911270914)
params.add('b5', value=149.49345053809822)

params.add('c', value=3.213258144928004)

out = lmfit.minimize(residual, params, args=(X, Y), method='lbfgsb', max_nfev=50000)
# out.nfev=15012 out.success=False out.message='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT'
print(f"Doesn't follow max_nfev > maxiter default:\n{out.nfev=} {out.success=} {out.message=}")

out = lmfit.minimize(residual, params, args=(X, Y), method='lbfgsb', max_nfev=5)
# out.nfev=6 out.success=False out.message='Fit aborted: number of function evaluations > 5'
print(f"\nDoes follow max_nfev < maxiter default:\n{out.nfev=} {out.success=} {out.message=}")

out = lmfit.minimize(residual, params, args=(X, Y), method='lbfgsb', maxiter=50000)
# out.nfev=15012 out.success=False out.message='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT'
print(f"\nRuntimeWarning because I used maxiter:\n{out.nfev=} {out.success=} {out.message=}")

Error message:

C:\Users\184277j\AppData\Local\Programs\Python\Python39\lib\site-packages\lmfit\minimizer.py:497: RuntimeWarning: ignoring `maxiter` argument to `Minimizer()`. Use `max_nfev` instead.
  warnings.warn(maxeval_warning.format(maxnfev_alias, 'Minimizer'),

Version information

Python: 3.9.4 (tags/v3.9.4:1f2e308, Apr 6 2021, 13:40:21) [MSC v.1928 64 bit (AMD64)]

lmfit: 1.2.0, scipy: 1.6.3, numpy: 1.21.4,asteval: 0.9.29, uncertainties: 3.1.7

Link(s)

#865

[newville]

@rowlesmr Thanks. It looks like the scipy l-bfgs-b code has keywords for both maxfun and maxiter, both defaulting to 15000. Neither the docs nor the scipy.optimize.show_options state that value, nor do they clarify what "an iteration" is.

Currently, when you say

lmfit.minimize(residual, params, args=(X, Y), method='lbfgsb', max_nfev=50000)

we actually convert that to

scipy.optimize.minimize(Minimizer.penalty, inital_values, method='L-BFGS-B', options={'maxiter': 100000})

where Minimizer.penalty will wrap the conversion of variable values to Parameters, call your residual function, and also actually count the number of function evaluations and halt if you exceed your value for max_nfev. In fact, if you had set max_nfev to 4000, that would have correctly aborted at 4000 (+ a few) function calls.

But we are not also setting maxfun, so that currently remains at the default value of 15000. OK, we can blame that as "very weird, probably buggy" behavior for the l-bfgs-b solver in scipy.optimize. But compared to all the other API oddities of scipy.optimize, that seems like a relatively minor annoyance. We should probably also set maxfun to 2*max_nfev for the l-bfgs-b method (as far as I can tell, no other method uses both), probably adding to

lmfit-py/lmfit/minimizer.py

Line 967 in a03ac0f

fmin_kws = dict(method=method, options={'maxiter': 2*self.max_nfev})

to read

        fmin_kws = dict(method=method, options={'maxiter': 2*self.max_nfev})
        if method == 'L-BFGS-B':
            fmin_kws['options']['maxfun'] = 2*self.max_nfev

I'm not 100% sure that is the best place or way to do this: it is effectively applying a bandage to a confusing API, but that is sort of the story on max_nfev anyway.

FWIW, as a workaround that I would be hesitant to recommend as a viable long-term solution, you could do:

   out = lmfit.minimize(residual, params, args=(X, Y), method='lbfgsb', options={'maxfun': 50000, 'maxiter': 50000})

to set both of these.

[rowlesmr]

The paper for algorithm 778 says (somewhere in the middle)

If the line search is unable to find a point with a sufficiently lower value of the objective after 20 evaluations of the objective function, we conclude that the current direction is not useful.

Maybe that is the (at least the original) meaning between function evaluations and iterations. Iterations are when parameter values are updated and function evaluations takes into account the line search per iteration. I think you mean (from only using this package since yesterday) max_nfev to mean iterations.

Maybe.

[newville]

@rowlesmr It would be sensible to guess that an iteration is between N_varys and 3*N_varys function calls (perhaps to first compute or update a derivative, then take a step, though maybe after the first step that is not needed .... ). But the documentation nor the paper (( think) are explicit and clear about that.

And the code (in scipy.optimize, I have not looked at the Fortran in detail) sets a default value for both maxiter and maxfun of 15000, which is honestly sort of weird. It is hard to imagine how one could have 15000 "iterations" with fewer than 15000 function calls. Also, most solvers have a default that scales with n_varys, and avoid two separate ways of limiting the number of function calls made.

But, I think we can work around this wart in scipy.optimize now that we're aware of it.

[rowlesmr]

No worries. I'm in the middle of a calculation right now, so I'll probably be able to test tomorrow sometime.

.

Also just noticed that it looks like the message I'm getting corresponds to exceeding the maximum function evaluations, and not the maximum iterations.

https://github.com/scipy/scipy/blob/c1ed5ece8ffbf05356a22a8106affcd11bd3aee0/scipy/optimize/_lbfgsb_py.py#L360-L370

Anyhoo, I'll leave it to you; I tried to follow what you're doing, but got lost in the wrapping and indirection. Thanks for the package!