fminunc throws an error; can't access source of error

When I use an optimization routine (fminunc) I get the error:

Error using lineSearch Search direction is not a descent direction; roundoff errors may be affecting convergence.

If there is going to be an error with one of the millions of iterations, I would prefer it just give an NAN for that value instead of exiting the entire code. The fact that it doesn't do this strikes me as a bug. I traced the source of the error to Matlab 2019\toolbox\optim\optim\private in a function called lineSearch.p

By design this function cannot be opened or modified. Is there any open source alternative to this code?

I reproduced the error by running the following lines of code (MatLab 2019a)

options = optimset('Display','off','TolFun',4e-16,'LargeScale','off');
par12=[4.45873897125075124848,2.45448132035209054536,398.23544583453281120455];
psfY12=[332.20478974188495158160,105.81515908458436570072,0,217.29169454128577854135,47.84122489713877257600,398.23544583453281120455,186.32064258820469149214,69.24275387165039319370,33.13754887984555352887];
y12=1:9;
fp1D = fminunc(@errfun,par12,options,psfY12,y12);
function [z] = errfun(p,v,x);
%cx = p(1);
%wx = p(2);
%amp = p(3);
zx = p(3)*exp(-0.5*(x-p(1)).^2./(p(2)^2)) - v;
z = sum(zx.^2);
end

[Cris LaPierre]

A better place to start might be having you share your code and data. It's hard to offer suggestions if we don't have anything to work with.

[Benjamin Brenner]

Here is a code that recreates the error

options = optimset('Display','off','TolFun',4e-16,'LargeScale','off');
par12=[4.45873897125075124848,2.45448132035209054536,398.23544583453281120455];
psfY12=[332.20478974188495158160,105.81515908458436570072,0,217.29169454128577854135,47.84122489713877257600,398.23544583453281120455,186.32064258820469149214,69.24275387165039319370,33.13754887984555352887];
y12=1:9;
fp1D = fminunc(@errfun,par12,options,psfY12,y12);


function [z] = errfun(p,v,x);

%cx = p(1);
%wx = p(2);
%amp = p(3);

zx = p(3)*exp(-0.5*(x-p(1)).^2./(p(2)^2)) - v;

z = sum(zx.^2);
end

Answer 1

[Matt J]

Since you are apparently trying to fit a Gaussian lobe to your data, you might find it helpful to use gaussfitn (which you must Download). However, the psfY12 data that you provided looks a very poor match to your Gaussian model. On the chance that the data was corrupted somehow, I demonstrate with simulated data below.

Fun=@(p,x) p(3)*exp(-0.5*(x-p(1)).^2./(p(2)^2)); %Gaussian model

par12=[4.45873897125075124848,2.45448132035209054536,398.23544583453281120455];
y12=1:9;
psfY12=Fun(par12,y12)+30*randn(size(y12)); %Simulate data

%%%% Do the fit
params=gaussfitn(y12(:),psfY12(:),{0,max(psfY12),mean(y12),[]},{0,0,0},{0,[],[]},...
                                  'Display','off','TolFun',1e-8,'MaxIter',1e6);

[p(3),p(1),p(2)]=deal(params{2:end}); 
p(2)=sqrt(p(2));


%%% Plot
fun=@(x) Fun(p,x);
t=linspace(min(y12),max(y12),1000);
plot(y12, psfY12,'--x',t, fun(t));
xlabel 'y12'; ylabel 'psfY12'; legend('Data samples','Gaussian fit')