I have a simple optimization problem that, with some specific data, makes scipy.optimize.minimize ignore the tol argument. From the documentation, tol determines the "tolerance for termination", that is, the maximum error accepted for the objective function, in my understanding (am I wrong?). However in the next working example, when tol is set to 0.1 for example, or other small numbers, the optimizations finishes with a "Optimization terminated successfully" message even when the objective function > tol. Is this a bug in Scipy's method or am I misunderstanding something here?
The optimization problem: I need to make a linear combination of var1 and var2, which are two time series, scaling them by parameters Btd and Bta. I need that the mean of the linear combination approximates to a target value Target, a scalar. So I simply minimize the absolute difference between np.mean(Btd*var1 + Bta*var2) and Target. The constraints are that the scaling coefficients must be >0 and that the ratio of means np.mean(Btd*var1)/np.mean(Bta*var2) should approximate to the function gi/(1-gi), where gi is a scalar in the interval [0,1].
Reproducible code:
import numpy as np
import scipy.optimize as opt
# The data that exactly reproduce the error:
time = np.arange(1979,2011)
var2=np.array([ 88.95705521, 74.5398773 , 72.08588957, 65.64417178,
50. , 72.39263804, 77.3006135 , 72.08588957,
64.41717791, 96.62576687, 69.93865031, 84.96932515,
86.50306748, 82.20858896, 80.98159509, 73.00613497,
66.25766871, 67.48466258, 79.75460123, 65.64417178,
70.24539877, 84.66257669, 76.3803681 , 83.74233129,
83.74233129, 78.2208589 , 88.03680982, 87.73006135,
100. , 71.16564417, 73.6196319 , 85.58282209])
var1=np.array([300. , 420.89552239, 333.58208955, 355.97014925,
376.11940299, 510.44776119, 420.89552239, 434.32835821,
333.58208955, 394.02985075, 523.88059701, 411.94029851,
353.73134328, 434.32835821, 355.97014925, 398.50746269,
476.86567164, 371.64179104, 445.52238806, 544.02985075,
416.41791045, 427.6119403 , 541.79104478, 579.85074627,
429.85074627, 414.17910448, 420.89552239, 528.35820896,
577.6119403 , 490.29850746, 600. , 454.47761194])
X=np.transpose([var1, var2])
# Global parameters
Target = 3.0
gi = 0.7
# This model is a simple linear combination of the two time series.
def MyModel(modelparams, X, gi):
Bta, Btd = modelparams
Eta = Bta*X[:,0]
Etd = Btd*X[:,1]
Etot = Eta + Etd
return Etot, Eta, Etd
# Objective function
def Obj(modelparams):
Bta, Bdt = modelparams
Etot, Eta, Etd = MyModel([Bta, Bdt], X, gi)
return abs(np.mean(Etot)-Target)
# Ratio constraint
def Ratio(modelparams):
import numpy as np
Bta, Btd = modelparams
Etot, Eta, Etd = MyModel([Bta, Btd], X, gi)
A = np.mean(Etd)/np.mean(Eta)
B = gi/(1-gi)
# The epsilon comes in to loosen a bit only this constraint
epsilon = 0.1
return -abs(abs(A-B)-epsilon)
# This is my solution to make the parameters different from zero.
# The ineq-type constraint makes them >=0.
def TDPos(modelparams):
Bta, Btd = modelparams
return Btd - 10**(-5)
def TAPos(modelparams):
Bta, Btd = modelparams
return Bta - 10**(-5)
constraints=[{'type': 'ineq', 'fun': Ratio},
{'type': 'ineq', 'fun': TDPos},
{'type': 'ineq', 'fun': TAPos}]
# Bounds or Model Parameters
bounds=((0, None), (0, None))
# Minimize
modelparams0=[Target/np.nanmean(var1), Target/np.nanmean(var2)]
result = opt.minimize(Obj, modelparams0,
tol=0.1,
method='SLSQP',
options={'maxiter': 40000 }, #,'ftol': 0.1},
bounds=bounds,
constraints=constraints)
print(result)
Prints out:
fun: 3.0
jac: array([439.92537314, 77.31019938])
message: 'Optimization terminated successfully.'
nfev: 20
nit: 4
njev: 4
status: 0
success: True
x: array([0., 0.])
My problem: fun: 3.0 > tol: 0.1 which is not desired.
TL;DR: scipy.optimize.minimize ignores the stop argument tol. Why?
EDIT: Moreover, the optimal solution [0, 0] ignores two of the ineq constraints, designed to make this couple of parameters > 10**(-5). Is this part of the same problem?
