Behaviour of least_squares fit for aborted fits

[eendebakpt]

In one of my applications I was optimizing on the reduced chi squared value. For one of the bad fits, the reduced chi squared reported by lmfit was negative infinity, which causes my application to pick this bad solution.

For the case at hand (minimal example below), the fitting is aborted with message Fit aborted: number of function evaluations > 14000. The values of ModelResult.residual and ModelResult.redchi are both -inf.

• The value -inf compares to other values such as 0.0012. A better value in case of an aborted fit could be NaN
• The fit has a valid best_values, so the residual can be calculated
• Other fitting procedures that abort early (Nelder-Mead, see minimal example below) do report a valid residual and redchi.

The change that introduced this behaviour is 5c15004 @newville.

Can we adapt the least-squares fitting code to:

• Use NaN instead of -inf to signal a bad fit?
• Calculate the residual, even if the fit was aborted (for some reason)?

Minimal example reproducing the behavior:

from lmfit import Model
import numpy as np

def gaussian(
    x, mean: float, std: float, amplitude: float = 1.0, offset: float = 0.0):
    x0 = x - mean
    y = offset + amplitude * np.exp(-x0 * x0 / (2.0 * std * std))
    return y


def double_gaussian(x, mean_left: float, mean_right: float,std_left: float, std_right: float, 
                        amplitude_left: float, amplitude_right: float ):
    gauss_left = gaussian(x, mean_left, std_left, amplitude_left)
    gauss_right = gaussian(x, mean_right, std_right, amplitude_right)
    y = gauss_left + gauss_right
    return y

x_data, y_data = (np.array([18.025, 18.075, 18.125, 18.175, 18.225, 18.275, 18.325, 18.375,
        18.425, 18.475, 18.525, 18.575, 18.625, 18.675, 18.725, 18.775,
        18.825, 18.875, 18.925, 18.975, 19.025, 19.075, 19.125, 19.175,
        19.225, 19.275, 19.325, 19.375, 19.425, 19.475, 19.525, 19.575,
        19.625, 19.675, 19.725, 19.775, 19.825, 19.875, 19.925, 19.975,
        20.025, 20.075, 20.125, 20.175, 20.225, 20.275, 20.325, 20.375,
        20.425, 20.475, 20.525, 20.575, 20.625, 20.675, 20.725, 20.775,
        20.825, 20.875, 20.925, 20.975]), np.array([ 43,  27,  14,  16,  20,  10,  15,   8,   3,   9,   2,   4,   8,
          3,   3,   1,   2,   4,   5,   0,   1,   3,   1,   0,   0,   1,
          0,   0,   0,   0,   0,   0,   1,  23, 202, 824, 344, 279, 276,
         23,   2,   7,   8,  17,  22,   3,  25,  24,  17,  99, 288, 172,
        252, 103,  24,  20,   4,   1,   3,   1], dtype=int))



initial_parameters = {'mean_left': 18.286083743842365,
 'mean_right': 20.07653394255875,
 'std_left': 0.0941689336296408,
 'std_right': 0.0741722783464804,
 'amplitude_left': 43,
 'amplitude_right': 824}

m=Model(double_gaussian)
m.set_param_hint("amplitude_left", min=0)
m.set_param_hint("amplitude_right", min=0)
r=m.fit(y_data, **initial_parameters, x=x_data,  method="least_squares")
print(r.success, r.message)
print(r.redchi) # is -inf, which seems odd


r=m.fit(y_data, **initial_parameters, x=x_data,  method="nelder_mead", max_nfev=2)
print(r.success,r.message)
print(r.redchi) # has a valid value

[newville]

@eendebakpt Thanks for the report. It looks like we are not handling that "abort because of too many function evals" with least_squares very well. I think we can (and should) fix that - as you suggest, calculate the residual and fit statistics with the last set of parameter values. In fact, we are doing that for other solvers ("if aborted, use last values, calculate all the stats"), just not "least_squares".

I'm also sort of surprised that "least_squares/trf" is not converging for this fitting problem. It doesn't look that hard to me, and "leastsq" seems to handle it OK. It is certainly possible that we could be using least_squares/trf better, but I think we are using basically the default settings....