# Inconsistency in optimize.minimize

I am trying to fit a time-dependent curve at each time step. I do so in minimizing along $x_c$ the quadratic error between the curve and a reference solution $1/(1 + \exp\left(\sqrt{S}(x-x_c)\right)$:

$$\frac{1}{N+2} \sum_{i=0}^{N+1} \left(y_i - \frac{1}{1 + \exp\left(\sqrt{S}(x_i-x_c)\right)}\right)^2$$

Because the calculations are easy, I calculated the gradient so as to use a conjugate gradient method implemented in optimize.minimize in scipy.

So: fairly easy stuff.

Below is my MWE and here is my data file:

import numpy as np
from matplotlib import pyplot as plt
from scipy import optimize

def energie(centre, x, y, S):
erreur = 0.
for i in range(0, y.size):
erreur += (y[i] - 1./(1.+np.exp(np.sqrt(S)*(x[i]-centre))))**2.
return erreur/y.size

for i in range(0, y.size):
grad += - (y[i] - 1./(1.+np.exp(np.sqrt(S)*(x[i]-centre)))) * np.sqrt(S) * \
np.exp(np.sqrt(S)*(x[i]-centre)) * (1.+np.exp(np.sqrt(S)*(x[i]-centre)))**-2.

T = np.arange(0., 36., 0.5/np.sqrt(2))
x = np.arange(301) * 0.5
centres = np.zeros(q.shape[0])

for n in range(q.shape[0]):
# initial guess
centreq = centres[n-1]
# call to optimize.minimize
resultatq = optimize.minimize(energie, centreq, method = 'CG', jac = gradenergie, args = (x, q[n , : ], 0.1))
# write to centres
centres[n] = resultatq.x

plt.plot(dcentresqdt[2 : ])
plt.show()


which produces the following graph: with obvious irregularities. Any idea why?

--

Edit:

in the code centre or centreq stands for $x_c$, each q[n, :] stands for the $(y_i)$, and the $(x_i)$ are the same at each time step, they are in the array x.

Indeed I plot the derivative of centres with respect to time, but if I plot centres itself I have the same kind of problem.

--

Edit2:

I realized my problem is so standard that there is actually a scipy function that does the job in no time (101 seconds for Nelder Mead vs 0.39 for the version below):

from scipy.optimize import curve_fit

def kink(x, centre):
S = 0.1
return 1./(1.+np.exp(np.sqrt(S)*(x-centre)))


and later on in the code:

for n in range(q.shape[0]):
centreq = centres[n-1]
popt, pcov = curve_fit(kink, x, q[n , : ], p0 = centreq)
centres[n] = popt

• I find it hard to tell from this question what the relationship is between $y_i,x_i,x_c$ in the formula and centres,q in the code. What exactly is being plotted here? It looks like there are multiple sets of data $\{(y^{(k)}_{0:n-1},x^{(k)}_{0:n-1})\}_k$ each leading to an $x_c^{(k)}$ and then $x_c^{(k+1)}-x_c^{(k)}$ is plotted as a function of $k$? It probably makes sense to you, but could you add a mathematical description that more precisely matches the code, please? – Kirill Jul 13 '18 at 16:14
• I just changed the method of minimization to Nelder-Mead and obtained a smooth curve. Although, I agree with @Kirill that it's not clear what the plot is about. – nicoguaro Jul 13 '18 at 17:39
• @Kirill you are right, I'm adding details. – bela83 Jul 14 '18 at 21:05