TNC does not return optimal parameters

[yannikschaelte]

The truncated Newton (TNC) optimizer scipy.optimize.minimize(method='TNC') returns apparently not the optimal parameter vector and function value, but the last one, which often (depending on the dimension of the problem in 10-30% of cases in our applications) has a higher function value.

Expected behavior: It should return the optimal found parameters and function value.

This behavior was not observed for any other scipy optimizer considered.

Reproducing code example:

We encountered this problem in our parameter estimation tool github.com/icb-dcm/pypesto, which employs an objective function which records a history of all function evaluations and can thus report if the best value in the trace is better than the value reported by the optimizer. See here for a discussion: ICB-DCM/pyPESTO#327. The error can be reproduced in cell 5 of the notebook https://github.com/ICB-DCM/pyPESTO/blob/master/doc/example/rosenbrock.ipynb. The problem only has box constraints, which are not violated.

Error message:

Function values from history and optimizer do not match: 1.3168387678656086, 2.0179911928533514
Parameters obtained from history and optimizer do not match: [0.99553374 0.98995083 0.97544456 0.94768476 0.89550744 0.8019584
 0.64185456 0.40180156 0.14794723 0.00844711], [0.9940831  0.98497745 0.96635783 0.93780681 0.87380595 0.74153588
 0.53395112 0.25523609 0.03957006 0.00100739]
Function values from history and optimizer do not match: 1.6349167203078765, 2.336760892779954

Scipy/Numpy/Python version information:

1.4.1 / 1.18.4 / 3.7.4

[andyfaff]

Case that demonstrates:

import numpy as np
from scipy.optimize import rosen, minimize

class Func():
    def __init__(self):
        self.nfev = 0
        self.best_x = None
        self.best_f = np.inf

    def __call__(self, x, *args, **kwds):
        f = rosen(x)
        self.nfev += 1
        if f < self.best_f:
            self.best_f = f
            self.best_x = x
        return f

f = Func()
res = minimize(f, [1, -2, 3, 4, 5], method='TNC')

assert f.nfev == res.nfev
np.testing.assert_equal(f.best_f, res.fun)
np.testing.assert_equal(f.best_x, res.x)

gives

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-8-22d412083be5> in <module>
      2 res = minimize(f, [1, -2, 3, 4, 5], method='TNC')
      3 
----> 4 np.testing.assert_equal(res.fun, f.best_fun)

~/miniconda3/envs/dev3/lib/python3.8/site-packages/numpy/testing/_private/utils.py in assert_equal(actual, desired, err_msg, verbose)
    426         # Explicitly use __eq__ for comparison, gh-2552
    427         if not (desired == actual):
--> 428             raise AssertionError(msg)
    429 
    430     except (DeprecationWarning, FutureWarning) as e:

AssertionError: 
Items are not equal:
 ACTUAL: 2.7271265087711303
 DESIRED: 2.4253707839139977

[andyfaff]

In the case of TNC there is a big discrepancy. It's obviously a long way from the minimum ([1, 1, 1, 1, 1]). In this particular example 'TNC' reaches the maximum number of iterations. With more iterations than the default we get a lot closer to the True minimum and the discrepancy gets much less (9.625356727615683e-10 vs 9.606951571011453e-10). However, they're still not the same.
The best_fun is probably done by an intermediate function evaluation (e.g. during line search, or gradient evaluation), and is not the same as the best value being tracked by the iteration loop in the minimizer. Ideally the two would be the same. That particular kind of issue would probably pop up for a lot of minimizers. I'm not sure if one should classify that as a bug because those techniques typically don't keep track of the best function value in the way our example uses.

If, however, there is a stale value not being updated in a loop somewhere, then that's a more serious issue.

One possible solution in this situation is to add best_fun and best_x attributes to ScalarFunction. Pretty much all the minimizers go through this object, and it would be easy to cache the best solution evaluated so far (and check it against the final solution). However, I'm not sure how desirable this behaviour is. One would have to do checks against all the constraints, and whether the bounds were satisfied, and the returned jac and hess values might not sync exactly with the x location.

[yannikschaelte]

Hi @andyfaff , thanks for creating the case demonstration. This seems to be the exact problem. I agree that just caching the best solution inside ScalarFunction may lead to problems as some routines may evaluate points that do not satisfy the constraints. If one implements it this way, one would need to 1) check that all constraints are satisfied, and 2) check if the gradient etc. are right. One could perform one additional gradient evaluation at the end if in doubt or return None. Additional problems may come up if a Hessian approximation is used, which needs to be sequentially refined over the optimizer course and thus cannot be reproduced from a single point.

Since we did not observe this behavior for any other optimizers (we had one possible case with L-BFGS, but I cannot reproduce that, so I'm not sure), it might be specific to TNC and there indeed a weird behavior of the algorithm, or a bug, as you point out.

[andyfaff]

So lets have a look at all the minimizers:

import numpy as np
from scipy.optimize import rosen, minimize, rosen_der, rosen_hess

class F():
    def __init__(self):
        self.best_fun = np.inf
        self.best_x = None

    def __call__(self, x):
        fun = rosen(x)
        if fun < self.best_fun:
            self.best_fun = fun
            self.best_x = x
        return fun

meths = ['nelder-mead', 'cg', 'bfgs', 'l-bfgs-b', 'tnc', 'cobyla', 'slsqp', 'trust-constr',
 'dogleg', 'trust-ncg', 'trust-exact', 'trust-krylov']

for method in meths:
    jac = hess = None
    if method in ('newton-cg', 'trust-krylov', 'trust-exact',
                  'trust-ncg', 'dogleg'):
        jac = rosen_der
    if method in ('trust-krylov', 'trust-exact', 'trust-ncg',
                  'dogleg'):
        hess = rosen_hess
            
    f = F()
    res = minimize(f, [1, -2, 3, 4, 5], method=method, jac=jac, hess=hess)
    print(method, np.equal(res.fun, f.best_fun), res.fun, f.best_fun, res.success, res.message)

gives

nelder-mead True 0.6841550038341955 0.6841550038341955 False Maximum number of function evaluations has been exceeded.
cg False 1.944273508185588e-10 1.9431177708679375e-10 True Optimization terminated successfully.
bfgs False 4.609168005831311e-11 4.437134484491673e-11 False Desired error not necessarily achieved due to precision loss.
l-bfgs-b False 2.2490199818718852e-11 2.2321289184993123e-11 True b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
tnc False 2.7271265087711303 2.4253707839139977 False Max. number of function evaluations reached
cobyla True 134.13260788002674 134.13260788002674 False Maximum number of function evaluations has been exceeded.
slsqp False 1.563858402399705e-08 1.5589201215269426e-08 True Optimization terminated successfully.
trust-constr False 5.143020768594074e-11 5.1369046800844246e-11 True `xtol` termination condition is satisfied.
dogleg True 11.718733542698956 11.718733542698956 False A linalg error occurred, such as a non-psd Hessian.
trust-ncg True 2.190618908717441e-15 2.190618908717441e-15 True Optimization terminated successfully.
trust-exact True 9.449221937588392e-16 9.449221937588392e-16 True Optimization terminated successfully.
trust-krylov True 1.549679356983471e-17 1.549679356983471e-17 True Optimization terminated successfully.

You'll see that the minimisation with TNC does not result in OptimizeResult.success being True. i.e. the minimisation failed. Ostensibly several other methods failed as well, mainly because iteration limits were reached. You definitely need to be checking the success attribute. I suspect in your examples you weren't.

Now let's revisit the minimum function value being visited by the objective. Imagine a solver whose mode of operation is fun(x) --> jac(x) --> fun(x).
Lets say on the second fun(x) evaluation (fval2) that the amount of fval decrease over the first fun(x) (fval1) has reached the criterion for stopping. At this point you'd expect the OptimizeResult.fun value to be the same as the second fval2. However, calculation of jac(x) uses finite differences (with many fun(x) evaluations). It may be that some of those FD evaluations may be at slightly lower fval than the second fval2. In that case the right thing to do is still return fval2, it's not a bug, rather the way that the solvers work. One such example is the l-bfgs-b above, for which OptimizerResult == True, but for which OptimizerResult.fun is slightly larger than that observed in the trace. The way to get to lower function values is to set stopping tolerances much tighter.

What would be a bug would be the optimizer halting based on the process I've just described, but with the OptimizeResult.fun being equal to fval1.

As such, I don't think the OP is evidence for an actual bug in the solvers.

[yannikschaelte]

Hi, thanks for the profound analysis.

You'll see that the minimisation with TNC does not result in OptimizeResult.success being True. i.e. the minimisation failed. Ostensibly several other methods failed as well, mainly because iteration limits were reached. You definitely need to be checking the success attribute. I suspect in your examples you weren't.

Yes, we did not check success as we employ a multi-start optimization where we usually have to live with whatever each of 100s of single starts returns. Of course makes sense to check, maybe we'll implement that for promising points, but just to clarify.

However, calculation of jac(x) uses finite differences (with many fun(x) evaluations). It may be that some of those FD evaluations may be at slightly lower fval than the second fval2. In that case the right thing to do is still return fval2, it's not a bug, rather the way that the solvers work.

I see your point about FD. Since we provided analytic gradients (which you did not for all optimizers above), I suspect far less optimizers showed this behavior than above, except some which evaluated points which were then not used to update the trace history, like, potentially, TNC. I agree that returning fval2 is the "right thing to do" from the search algorithm perspective. If one is just interested in overall best parameter values, it is of course always possible to manually self-record them, as in Func above, which we also do.

As such, I don't think the OP is evidence for an actual bug in the solvers.

Agreed. It was striking that TNC showed this behavior quite often compared to the others (providing gradients), and seems to always return the last evaluated function value, not the best, the difference being beyond any tolerance threshold. Whether this is a flaw of the algorithm (corresponding to returning fval1, not fval2) or not this use case fails to reveal, this I guess only a code analysis could clarify.

[yannikschaelte]

Here's for completion an illustration of how only TNC fails in gradient-based mode:

import numpy as np
from scipy.optimize import rosen, minimize, rosen_der, rosen_hess

class F():
    def __init__(self):
        self.best_fun = np.inf
        self.best_x = None

    def __call__(self, x):
        fun = rosen(x)
        if fun < self.best_fun:
            self.best_fun = fun
            self.best_x = x
        return fun

meths = ['nelder-mead', 'cg', 'bfgs', 'l-bfgs-b', 'tnc', 'cobyla', 'slsqp', 'trust-constr',
 'dogleg', 'trust-ncg', 'trust-exact', 'trust-krylov']

for method in meths:
    jac = hess = None
    #if method in ('newton-cg', 'trust-krylov', 'trust-exact',
    #              'trust-ncg', 'dogleg'):
    jac = rosen_der
    #if method in ('trust-krylov', 'trust-exact', 'trust-ncg',
    #              'dogleg'):
    hess = rosen_hess
    x0 = [-0.16584467,  1.31897148, -3.26098005,  3.45927947,  0.29453825,
        1.31031958, -1.75957634, -1.83706183,  3.42357846, -0.78687309]
    f = F()
    res = minimize(f, x0, method=method, jac=jac, hess=hess, options={'maxiter': 100})
    print(method, np.equal(res.fun, f.best_fun), res.fun, f.best_fun, res.success, res.message)

gives

nelder-mead True 2120.269071102054 2120.269071102054 False Maximum number of iterations has been exceeded.
cg True 0.26533809379430673 0.26533809379430673 False Maximum number of iterations has been exceeded.
bfgs True 2.5370442739135357e-14 2.5370442739135357e-14 False Maximum number of iterations has been exceeded.
l-bfgs-b True 2.7822653573972946e-09 2.7822653573972946e-09 True b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
tnc False 4.755755190559484 3.987300664638279 False Max. number of function evaluations reached
cobyla True 69.73757558350557 69.73757558350557 False Maximum number of function evaluations has been exceeded.
slsqp True 2.6276887439165257 2.6276887439165257 False Iteration limit exceeded
trust-constr True 8.597660362768444e-21 8.597660362768444e-21 True `gtol` termination condition is satisfied.
dogleg True 40923.23016898274 40923.23016898274 False A linalg error occurred, such as a non-psd Hessian.
trust-ncg True 7.202214306674385e-21 7.202214306674385e-21 True Optimization terminated successfully.
trust-exact True 6.761705476594788e-16 6.761705476594788e-16 True Optimization terminated successfully.
trust-krylov True 4.911506820708559e-13 4.911506820708559e-13 True Optimization terminated successfully.

Indeed, setting maxiter=1000, the found values match (True) for all optimizers, indicating this may only occur when convergence was not reached.

[andyfaff]

Closing because I think the following is True:

However, calculation of jac(x) uses finite differences (with many fun(x) evaluations). It may be that some of those FD evaluations (or evaluations in linesearch) may be at slightly lower fval than the second fval2.

Note that some of the examples had res.success=False because of hitting the maxfun limit. If this is happening it's not a True demonstration of the issue; set maxfun higher and repeat.

Furthermore, there is an issue when using TNC with a callback function, the callback mechanism in the TNC code was mutating the output. This is being fixed in #14882.