fmin_slsqp not converging to sensible solution

[astrofrog]

I am trying to use the fmin_slsqp function while fitting a Gaussian to fake data:

import numpy as np
from scipy.optimize import fmin_slsqp

w = np.linspace(1, 2, 100)
f = 1.2e-8 * np.exp(-0.5 * (w-1.4)**2 / 0.11**2)

def gaussian(p, x):
    return p[0] * np.exp(-0.5 * (x - p[1])**2 / p[2]**2)

def errfunc(p, x, y):
    return np.sum((gaussian(p, x) - y)**2)

p0 = [1.5e-8, 1.6, 0.2]
p_new = fmin_slsqp(errfunc, x0=p0, args=(w, f))

print('p0   ', p0)
print('p_new', p_new)

However this does not converge to the expected solution despite the reasonable initial guess:

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 4.793585605039844e-15
            Iterations: 1
            Function evaluations: 5
            Gradient evaluations: 1
p0    [1.5e-08, 1.6, 0.2]
p_new [  1.50000000e-08   1.60000000e+00   2.00000000e-01]

I suspected that this was due to the acc tolerance which may be absolute, so I changed it to:

p_new = fmin_slsqp(errfunc, x0=p0, args=(w, f), iter=10000, acc=1.e-40)

However while this now does many iterations instead of one, it still does not converge:

Iteration limit exceeded    (Exit mode 9)
            Current function value: 2.7710800492819493e-15
            Iterations: 10001
            Function evaluations: 149997
            Gradient evaluations: 10001
p0    [1.5e-08, 1.6, 0.2]
p_new [  2.15505669e-11   1.60000000e+00   2.00000000e-01]

It seems like this could be a bug?

I'm using Python 3.6.1 with Scipy 0.19.0 and Numpy 1.11.3.

(cc @ibusko - just so you are in the loop)

[jaimefrio]

I have been playing with your code, and it seems to fail for all optimizers, not just 'slsqp'. It's easier to try various if you use optimize.minimize rather than the algorithm-specific methods, e.g.:

p_new = optimize.minimize(errfunc, x0=p0, args=(w, f), method='slsqp')

The other interesting thing is that if you multiply your whole function by 1e+8, and scale accordingly the first starting point, then it converges to the right solution without any problems. You may be able to work around this by playing with the eps option, that controls the step used for gradient approximation, but my quick tests didn't take me anywhere...

So I don't really think it is a bug in scipy.optimize, but more of a numerically ill-conditioned problem being thrown at it.

[astrofrog]

@jaimefrio - I haven't had this issue with all other optimizers. For example, if I just use fmin, it works:

import numpy as np
from scipy.optimize import fmin

w = np.linspace(1, 2, 100)
f = 1.2e-8 * np.exp(-0.5 * (w-1.4)**2 / 0.11**2)

def gaussian(p, x):
    return p[0] * np.exp(-0.5 * (x - p[1])**2 / p[2]**2)

def errfunc(p, x, y):
    return np.sum((gaussian(p, x) - y)**2)

p0 = [1.5e-8, 1.6, 0.2]
p_new = fmin(errfunc, x0=p0, args=(w, f))

print('p0   ', p0)
print('p_new', p_new)

with the following output:

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 76
         Function evaluations: 139
p0    [1.5e-08, 1.6, 0.2]
p_new [  1.20061578e-08   1.39997215e+00   1.09967149e-01]

The p_new parameters match the input ones.

I don't think this is a numerically ill-conditioned problem but more that some algorithms may have built-in absolute tolerances or other issues. I will try and look into it further.

[pv]

Several of the optimizers may have assumptions the problem is well scaled. This is not necessarily only a question of absolute tolerances.

[pv]

Namely, here it is the first variable that has a very different scale than the rest. An algorithm using e.g. spherical trust regions will encounter problems. I don't recall outright what SLSQP exactly does, but maybe it does just that.

[astrofrog]

@pv - I'm reaching the same conclusion. I can get better (though still not ideal) results if I modify slsqp.py so that epsilon can be a vector - at the moment epsilon is used as an additive difference on all parameters, which doesn't work well if some of the parameters have a very different scale.

[pv]

The scaling of the problem is pretty much extra information that you pass to the solver --- you can pass it to the solver by ensuring typical parameters and values are of order 1. . It is in principle possible to add some logic that tries to rescale the problem, for example by looking at the jacobian etc. But this is probably not always reliable, especially if the initial point is not representative of the solution.

[astrofrog]

@pv - ok, that makes sense, thanks! I have a way to fix the scaling in an automated way, but I'm still running into issues even when the scales are appropriate and the initial parameters are a bit off. Consider the following example:

import numpy as np
from scipy.optimize import fmin_slsqp

w = np.linspace(1, 2, 100)
f = np.exp(-0.5 * (w-1.4)**2 / 0.11**2)

def gaussian(p, x):
    return p[0] * np.exp(-0.5 * (x - p[1])**2 / p[2]**2)

def errfunc(p, x, y):
    return np.sum((gaussian(p, x) - y)**2)

p0 = [3, 0.85, 0.1]
p_new = fmin_slsqp(errfunc, x0=p0, args=(w, f))

print('p0   ', p0)
print('p_new', p_new)

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(w, f, 'k.')
ax.plot(w, gaussian(p0, w), 'r-')
ax.plot(w, gaussian(p_new, w), 'g-')
fig.savefig('fit.png')

Is this also expected? Other optimizers (e.g. fmin) don't have any issue with such cases, but is this just a peculiarity of SLSQP that it is very sensitive to the initial conditions?

[ibusko]

Re-running my test script with the data normalized to ~1., I found that SLSQP improves tremendously. It’s able to find the solution, and also appears to be relatively insensitive to the initial guesses for the parameters. Other methods seem to not care at all for the small values of the dependent variable (~1.e-15), with the exception of L-M (the MINPACK implementation) which changes its behavior in a funky way depending on the data normalization. I still have to work more on that topic. Ivo

…

On Jun 23, 2017, at 9:50 AM, Thomas Robitaille ***@***.***> wrote: @pv - ok, that makes sense, thanks! I have a way to fix the scaling in an automated way, but I'm still running into issues even when the scales are appropriate and the initial parameters are a bit off. Consider the following example: import numpy as np from scipy.optimize import fmin_slsqp w = np.linspace(1, 2, 100 ) f = np.exp(-0.5 * (w-1.4)**2 / 0.11**2 ) def gaussian(p, x ): return p[0] * np.exp(-0.5 * (x - p[1])**2 / p[2]**2 ) def errfunc(p, x, y ): return np.sum((gaussian(p, x) - y)**2 ) p0 = [3, 0.85, 0.1 ] p_new = fmin_slsqp(errfunc, x0=p0, args= (w, f)) print('p0 ' , p0) print('p_new' , p_new) import matplotlib.pyplot as plt fig = plt.figure() ax = fig.add_subplot(1,1,1 ) ax.plot(w, f, 'k.' ) ax.plot(w, gaussian(p0, w), 'r-' ) ax.plot(w, gaussian(p_new, w), 'g-' ) fig.savefig( 'fit.png') Is this also expected? Other optimizers (e.g. fmin) don't have any issue with such cases, but is this just a peculiarity of SLSQP that it is very sensitive to the initial conditions? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.