BUG: Regression in 1.11.3 can still fail for optimize.least_squares with method='trf', maybe due to 0 lower bounds?

[rhugonnet]

Describe your issue.

Hi @tylerjereddy (thanks for all the maintenance work!)

This issue follows up on #19103, maybe also echoes #19309 (but the error message is not linked to division by zero).

Optimization with optimize.least_squares and method='trf' fails since 1.11.2, and for some cases it was not fixed in 1.11.3 (but some were). We noticed this from our CI test suite in SciKit-GStat (https://github.com/mmaelicke/scikit-gstat).
I'm not sure what the origin is, maybe related to lower bounds equal to zero?
Reproducible example below.

Reproducing Code Example

Here's an example that:

• Passes with SciPy 1.11.1,
• Fails with SciPy 1.11.2,
• Still fails with SciPy 1.11.3 despite the fix in BUG: Fixes #19103 by adding back make_strictly_feasible to lsq trf. #19111,

import numpy as np
from scipy.optimize import curve_fit
from functools import wraps
import math

def variogram(func):
    @wraps(func)
    def wrapper(*args, **kwargs):
        if hasattr(args[0], '__iter__'):
            new_args = args[1:]
            mapping = map(lambda h: func(h, *new_args, **kwargs), args[0])
            return np.fromiter(mapping, dtype=float)
        else:
            return func(*args, **kwargs)
    return wrapper

def spherical(h, r, c0, b=0):
    a = r / 1.

    if h <= r:
        return b + c0 * ((1.5 * (h / a)) - (0.5 * ((h / a) ** 3.0)))
    else:
        return b + c0

def gaussian(h, r, c0, b=0):
    # prepare parameters
    a = r / 2.

    return b + c0 * (1. - math.exp(- (h ** 2 / a ** 2)))

@variogram
def model_sum(h, r1, c1, r2, c2):

    return spherical(h, r1, c1) + gaussian(h, r2, c2)

f = model_sum

xdata = np.array([ 12.64307503,  25.28615006,  37.9292251 ,  50.57230013,
                   63.21537516,  75.85845019,  88.50152522, 101.14460026,
                   113.78767529, 126.43075032])
ydata = np.array([10.33885898, 20.33271591, 17.6010188 , 18.12232903, 19.60947486,
                  16.96256666, 17.08413617, 18.35997835, 13.34964295, 16.91325055])
p0 = np.array([126.43075032,  20.33271591, 126.43075032,  20.33271591])
bounds = (0, [126.43075032075254, 20.33271591029816, 126.43075032075254, 20.33271591029816])

curve_fit(f, xdata, ydata, p0=p0, bounds=bounds)

Error message

As for previous issues on this topic, the error message is (in both 1.11.2 and 1.11.3):

RuntimeError: Optimal parameters not found: The maximum number of function evaluations is exceeded.

SciPy/NumPy/Python version and system information

SciPy 1.11.1 to 1.11.3.

Thanks for looking into this!

[tylerjereddy]

I had a pretty good suspicion that the problem wasn't entirely fixed as you can see on the release discussion here: #19279 (comment), and obviously from the issue you linked, which we knew about before I did the 1.11.3 release.

To get these issues squashed more decisively I probably need some domain expertise to chime in. I think we usually ping @nmayorov for this corner of optimize.

@mdhaber also has more knowledge of optimization algorithms than I do in general (pretty sure at least). Any thoughts on the sequence of patches below and the issues that remain open?

Patch 1: #18896
Patch 2: #19111
(both sets of regression tests passed when combined)

New issues that remain and/or popped up as a result of merging those patches:
#19309
#19351

[mdhaber]

Hmm thanks but unfortunately I don't know much about these algorithms, either. I would defer to @nmayorov, too.

[nmayorov]

I think the wisest decision is to restore the code prior to #18896

The original complaint (about x not being exactly 0) I believe is not substantial or at least not critical. It is rather a small issue or a special property of the algorithm, certainly not a bug.

[tylerjereddy]

Maybe the docs could briefly explain the special property so folks don't think it is a bug?

[nmayorov]

I plan to create a PR with the revert + the test adjustment + small clarification of what to expect from the algorithms.