Minimization with one residual term does not apply reduce_fcn in scalar_minimize

[hstaige]

First Time Issue Code

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

Description

When conducting a minimization with one residual term, the reduce function is sometimes not called. This can lead to the minimization converging to the wrong value or failing to converge entirely.

I believe this is due to minimizer.penalty checking if the array.size > 1 instead of >= 1.

A Minimal, Complete, and Verifiable example

import lmfit

def target1(params):
    return ((params['R0'] - 1) / 0.1)**2

def target2(params):
    return [(params['R0'] - 1) / 0.1]

params = lmfit.Parameters()
params.add('R0', 0)

mini_res1 = lmfit.minimize(target1, params, scale_covar=False, method='nelder') # first fit with direct squaring works fine
print(mini_res1.success)
mini_res2 = lmfit.minimize(target2, params, scale_covar=False, method='nelder') # second fit fails to converge

Fit report:

Not applicable, second fit crashes due to NaN value

Error message:

C:\Users\Hunter\AppData\Local\Temp\ipykernel_22232\4097992378.py:7: RuntimeWarning: overflow encountered in double_scalars
  return [(params['R0'] - 1) / 0.1]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [6], in <cell line: 14>()
     12 mini_res1 = lmfit.minimize(target1, params, scale_covar=False, method='nelder')
     13 print(mini_res1.success)
---> 14 mini_res2 = lmfit.minimize(target2, params, scale_covar=False, method='nelder')

File ~\anaconda3\envs\euv_test\lib\site-packages\lmfit\minimizer.py:2583, in minimize(fcn, params, method, args, kws, iter_cb, scale_covar, nan_policy, reduce_fcn, calc_covar, max_nfev, **fit_kws)
   2443 """Perform the minimization of the objective function.
   2444 
   2445 The minimize function takes an objective function to be minimized,
   (...)
   2577 
   2578 """
   2579 fitter = Minimizer(fcn, params, fcn_args=args, fcn_kws=kws,
   2580                    iter_cb=iter_cb, scale_covar=scale_covar,
   2581                    nan_policy=nan_policy, reduce_fcn=reduce_fcn,
   2582                    calc_covar=calc_covar, max_nfev=max_nfev, **fit_kws)
-> 2583 return fitter.minimize(method=method)

File ~\anaconda3\envs\euv_test\lib\site-packages\lmfit\minimizer.py:2352, in Minimizer.minimize(self, method, params, **kws)
   2349         if (key.lower().startswith(user_method) or
   2350                 val.lower().startswith(user_method)):
   2351             kwargs['method'] = val
-> 2352 return function(**kwargs)

File ~\anaconda3\envs\euv_test\lib\site-packages\lmfit\minimizer.py:1020, in Minimizer.scalar_minimize(self, method, params, max_nfev, **kws)
   1018 result.call_kws = fmin_kws
   1019 try:
-> 1020     ret = scipy_minimize(self.penalty, variables, **fmin_kws)
   1021 except AbortFitException:
   1022     pass

File ~\anaconda3\envs\euv_test\lib\site-packages\scipy\optimize\_minimize.py:680, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    677     bounds = standardize_bounds(bounds, x0, meth)
    679 if meth == 'nelder-mead':
--> 680     res = _minimize_neldermead(fun, x0, args, callback, bounds=bounds,
    681                                **options)
    682 elif meth == 'powell':
    683     res = _minimize_powell(fun, x0, args, callback, bounds, **options)

File ~\anaconda3\envs\euv_test\lib\site-packages\scipy\optimize\_optimize.py:825, in _minimize_neldermead(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, adaptive, bounds, **unknown_options)
    823 if bounds is not None:
    824     xe = np.clip(xe, lower_bound, upper_bound)
--> 825 fxe = func(xe)
    827 if fxe < fxr:
    828     sim[-1] = xe

File ~\anaconda3\envs\euv_test\lib\site-packages\scipy\optimize\_optimize.py:496, in _wrap_scalar_function_maxfun_validation.<locals>.function_wrapper(x, *wrapper_args)
    494 ncalls[0] += 1
    495 # A copy of x is sent to the user function (gh13740)
--> 496 fx = function(np.copy(x), *(wrapper_args + args))
    497 # Ideally, we'd like to a have a true scalar returned from f(x). For
    498 # backwards-compatibility, also allow np.array([1.3]),
    499 # np.array([[1.3]]) etc.
    500 if not np.isscalar(fx):

File ~\anaconda3\envs\euv_test\lib\site-packages\lmfit\minimizer.py:655, in Minimizer.penalty(self, fvars)
    652 else:
    653     apply_bounds_transformation = True
--> 655 r = self.__residual(fvars, apply_bounds_transformation)
    656 if isinstance(r, np.ndarray) and r.size > 1:
    657     r = self.reduce_fcn(r)

File ~\anaconda3\envs\euv_test\lib\site-packages\lmfit\minimizer.py:601, in Minimizer.__residual(self, fvars, apply_bounds_transformation)
    599     raise AbortFitException("fit aborted by user.")
    600 else:
--> 601     return _nan_policy(np.asarray(out).ravel(),
    602                        nan_policy=self.nan_policy)

File ~\anaconda3\envs\euv_test\lib\site-packages\lmfit\minimizer.py:2436, in _nan_policy(arr, nan_policy, handle_inf)
   2431     if contains_nan:
   2432         msg = ('NaN values detected in your input data or the output of '
   2433                'your objective/model function - fitting algorithms cannot '
   2434                'handle this! Please read https://lmfit.github.io/lmfit-py/faq.html#i-get-errors-from-nan-in-my-fit-what-can-i-do '
   2435                'for more information.')
-> 2436         raise ValueError(msg)
   2437 return arr

ValueError: NaN values detected in your input data or the output of your objective/model function - fitting algorithms cannot handle this! Please read https://lmfit.github.io/lmfit-py/faq.html#i-get-errors-from-nan-in-my-fit-what-can-i-do for more information.`


Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  ...

Version information

Python: 3.13.0 | packaged by Anaconda, Inc. | (main, Oct 7 2024, 21:21:52) [MSC v.1929 64 bit (AMD64)]

lmfit: 1.3.2, scipy: 1.14.1, numpy: 2.1.3,asteval: 1.0.5, uncertainties: 3.2.2

Link(s)

Relevant code block

[reneeotten]

The first problem is that NaN values are detected - you will need to make sure that that cannot happen.

Also, this doesn't look like a complete example as you're not actually fitting any data.

[hstaige]

One change I can test is adjusting this check from

if isinstance(r, np.ndarray) and r.size > 1:

to

if instance(r, np.ndarray) and r.size >= 1 and r.shape != ():

This returns true for np.asarray([1]) (userfcn returns a list with one element) and false for np.asarray(1) (userfcn returns a scalar).

I am not too familiar with GitHub, but I'll try for an hour and see if this passes all existing tests or if it breaks things fundamentally.

[newville]

@hstaige In my opinion, we do not need a Pull Request, because no changes are needed. I think this is the documented behavior -- the scalar minimizers minimize the value returned by the objective function, not necessarily in the least squares sense.

If the objective returns an array of length > 1, lmfit Minimize will reduce that array to a scalar, with sum-of-squares as the default array-to-scalar reduction method, and others available. In this way, the same objective function can be used with multiple methods.

When a scalar is returned, we do not reduce it to a scalar. A one-element array is ambiguous... we treat it as a single value, allowing us to convert all outputs to ndarrays.

Your function target2() returns scalar, so I think the current behavior is correct. Note that if you had

def target2(params):
    return [0.0, (params['R0'] - 1) / 0.1]

The value for R0 would be 1, as the need to reduce the list/array would be clear. That is not to say that this is the recommended approach, though it would allow non-scalar algorithms. Using your target1 seems reasonable, and other "cost functions" could be explored too.

But, to put this as a clear statement: if you are using a scalar minimizer and your objective function returns a scalar, that scalar will be minimized.

And, that minimum could go to -Inf.

But, I'll also ask @reneeotten and anyone else interested if that seems right to them.

[hstaige]

I tried a few things, but, as you stated above, it is already assumed in a few locations that a function that returns a single element array should be interpreted as a scalar. I can modify my code to respect this.

It is odd that different fitting methods minimize different quantities, but I think that is an issue with scipy, not lmfit.

[newville]

@hstaige The (many) scalar minimizers (Nelder-Mead, Powell, etc) do act the same to minimize a scalar value, and scipy.optimize.minimize() wraps those, as does lmfit.minimize().
The multi-variate minimizers (leastsq, least_squares) want an array to minimize in the least squares sense (with some options to use other cost functions).

Lmfit tries to give a uniform interface so that if your objective returns an ndarray (or list of numbers) instead of a scalar value, the scalar minimizers will not panic, or do something silly like use only the first value - it will make a sum-of-squares (by default). That means you do not have to change your objective function to switch solvers.

But if your solver returns a scalar, there is nothing to sum, and that objective would not work with least squares anyway (as the return array must have more elements than variables). So there is no need to square the returned scalar value.

Hope that makes sense.

[hstaige]

I understand the logic and agree that subtly changing the behavior of lmfit would likely break more code than it fixes. I am fine with closing this discussion then.

Thank you both for your time here and in developing lmfit!