Local minimum at initial point when fitting gaussian with `lsqcurvefit`

Question

I am writing MATLAB code with the intention to do some fittings. I simulated a plot using a second-order Gaussian (see my code below) and tried fitting using the lsqcurvefit function. Unfortunately, MATLAB returns the same guess values as 'optimized' parameters and apparently gets stuck in a local minimum.

Can someone please provide some advice as to what might be wrong here? I know that if the guess is far from the 'true' values then this can happen, but I expected MATLAB to return an answer which was closer to true value. Improving the initial guess to [29,0] (which is much closer to the actual value) results in the same output: that the initial value is a local minimum.

%%%%%%%%%%

function z= testfunction(x, xdata);

sigma=x(1)/(2*sqrt(2*log(2)));

z=((xdata.^2-2*x(2)*xdata-sigma.^2+x(2)^2)./(sigma^5*sqrt(2*pi))).*exp(-(xdata-x(2)).^2/(2.*sigma.^2));

end
%%%%%%%%


% Simulate Data

xdata= -50:1:50;

ydata = testfunction([30,0],xdata);

% Fit Data 

xfit = lsqcurvefit(@testfunction,[19,-4],xdata, ydata );

xfit(1)
xfit(2)

yfit=testfunction([xfit(1),xfit(2)],xdata);

% Plot Data;
plot(xdata,yfit);
hold on;
plot(xdata,ydata,'o')

Output:

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point 
is less than the default value of the optimality tolerance.

<stopping criteria details>


ans =

    19


ans =

    -4

Answer 1

Short answer: check the stopping criteria details and change the stopping criteria accordingly:

options = optimoptions('lsqcurvefit','OptimalityTolerance', 1e-16, 'FunctionTolerance', 1e-16);
xfit = lsqcurvefit(@testfunction,[19,-4],xdata, ydata, [], [], options);

What is the problem?

lsqcurvefit is a numerical solver and therefore uses a stopping criteria to determine if the local minimum is sufficiently reached. In general, you never reach the exact solution. Therefore, the solution to your problem is to change the stopping criteria to request a more accurate solution.

How to solve it?

By clicking on stopping criteria details, you receive the following explanation:

Optimization completed: The final point is the initial point. The first-order optimality measure, 7.254593e-07, is less than options.OptimalityTolerance = 1.000000e-06.

Optimization Metric                                       Options 
relative first-order optimality =   7.25e-07  OptimalityTolerance =   1e-06 (default)

So, you should decrease the OptimalityTolerance (f.e. to 1e-16):

options = optimoptions('lsqcurvefit','OptimalityTolerance', 1e-16);
xfit = lsqcurvefit(@testfunction,[19,-4],xdata, ydata, [], [], options);

The above image visualizes the new result, which is better than the previous one, but not yet very good. By inspecting the stopping criteria again, you will see that you also need to change FunctionTolerance:

options = optimoptions('lsqcurvefit','OptimalityTolerance', 1e-16, 'FunctionTolerance', 1e-16);

Why are the default options that bad?

Note that you need to tweak the stopping criteria, because your function returns relative small values. Multiplying z with a factor of 1000, without any option specification would also result in a nice fit: