How to supply user gradient in simultaneous fitting with lsqcurvefit

xdata={[1 3 5 7]' [5 7 9 11]' [9 11 13 15]'};
y={[10 11 12 13]' [10 11 12 13]' [10 11 13]'};
yy = vertcat(y{:});
%variables v(1:4) are m, c1, c2, c3
fn = @(v,xdata) [v(1).*xdata{1} + v(2); v(1).*xdata{2} + v(3); v(1).*xdata{3} + v(4)];
x0 = [1; 10; 9; 8];
fitvars = lsqcurvefit(fn,x0,xdata,yy);

Based on the above code snippet example is there any way that we can tell lsqcurvefit to use user supplied gradients for each function?

Answer 1

[John D'Errico]

This problem is linear in the parameters. Supplying the gradient is a waste of time. In fact, using lsqcurvefit is a waste of time too, since you could solve the problem using backslash.

But since lsqcurvefit will differentiate the problem itself, and since the derivatives of a linear function are rather trivial to estimate exactly, lsqcurvefit will work adequately well. Just far more slowly than you could have done it using backslash.

If you want to supply the gradients for some reason, like you think lsqcurvefit is not doing sufficiently well, don't bother. If you think it too slow, then you need to learn to formulate the problem for solution using backslash. The funny thing is, that the gradient that you wanted to learn to supply to lsqcurvefit is the same array that you need to generate to use backslash.

As far as actually writing any code, I would point out that you have 11 values in xdata, and 12 in y. So nothing you do will really work here.

[John D'Errico]

Please don't answer just to make a comment. Moved your answer to a comment:

"Dear John first of all thank you very much for your time answering my question. You are absolutely right this problem besides the missing value may be formulated and solved with the backslash operator easily. In my actual problem the parameters enter non-linearly and I would like to utilize my analytic gradients in lsqcurvefit for speed and accuracy. This code snippet serves as an analogy, a bad one though since it's linear in the parameters.

The main characteristics though are similar, I have an anonymous vector of functions which I would like to apply for different Xdata-Ydata values while the fitting parameters are shared. The lsqcurvefit works perfectly with central finite differences but because I want to call the solver from 10^6 different starting points, (via multistart) I am looking for a way to implement each functions' gradient in order to signifficantly reduce the total function evaluations and thus overall runtime.

The error I get is "Too many output arguments. etc.." Is there any work around?

Thanks again in advance"

[John D'Errico]

It is often true that analytical gradients do not improve either speed OR accuracy. People seem to deceive themselves into thinking this will be an improvement, yet it rarely really helps. Very often, the gradient computation itself takes more time to compute than multiple calls to the objective, so there is little speed gain. As well, tiny errors in the gradient from a well chosen finite difference approximation, as long as they do not result in a non-descent direction, are not a problem.

The error of too many output arguments is something I cannot resolve, since I do not have your code. You have not written the code so it runs porperly. The fix is to write proper code. Sorry. That may seem trite, but how can I say it better?

[Pas]

I am sorry for my mistake thank you form moving it appropriately. The reason i believe it will offer a benefit is that in simple cases when a single function is used with the gradient the total function evaluations are an order of magnitude less than when central finite differences, (CFD) are used.Specifically for the simplified version of my code shown below when I use the Jacobian I get (1132 Fcounts vs 10179 Fcounts with CFD and a speed up for a single core of ~20%. But when i run problem2 with 'SpecifyObjectiveGradient',true I get the following error:

Error using vertcat, Too many output arguments.

Error in @(m,xdata)[fun1(m,xdata(1:19,1));fun1(m,xdata(20:end,1))]

Error in lsqcurvefit (line 213)

[initVals.F,initVals.J] =

feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in fmultistart Error in MultiStart/run (line 268)

fmultistart(problem,startPointSets,msoptions);

Caused by:

Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

Failure in call to the local solver with supplied problem structure.

Any ideas?

% Clear temporary variables & Workspace
clear all;clc
%
% Xdata
xk=[-0.0449531768;-0.0615182117;-0.0717112888;-0.0790843544;-0.0848645537;
     -0.0896207495;-0.0936631056;-0.0971793402;-0.100291694;-0.103084217;
     -0.121658488;-0.132708653;-0.140652425;-0.146885302;-0.152031629;
     -0.15642554;-0.160267022;-0.163685376;-0.166769074];
    	
% Ydata
yk= [(-0.0001:-0.0001:-0.0009),(-0.001:-0.001:-0.01)]';
%
% For this example i am using the same data twice
xdata=[xk];
ydata=[yk];

xdata2=[xk;xk];
ydata2=[yk;yk];

% Lower and Upper bounds
lb=[-6 -12 -12 -12];
ub=[6 12 12 12];

l1 = lb(1);
u1 = ub(1);
l2 = lb(2);
u2 = ub(2);
l3 = lb(3);
u3 = ub(3);
l4 = lb(4);
u4 = ub(4);

% Create a handle for the objective(s) function 
fun1=@icellT;

% Create an array of the data respective functions 
fun= @(m,xdata)[fun1(m,xdata(1:19,1))];

% This works only with Central Differences no User Gradients
fun2= @(m,xdata2)[fun1(m,xdata2(1:19,1));fun1(m,xdata2(20:end,1))];

runs=100;  %24*41667 = 1M random Quadruples (my rig has 24 cores)

rng(9) % for reproducibility 
M = zeros(runs,4);
m0=zeros(runs,4);
m01=zeros(1,runs);
m02=zeros(1,runs);
m03=zeros(1,runs);
m04=zeros(1,runs);

% Create uniformly distributed bounded random starting points 
m01 =(unifrnd(l1,u1,1,runs));
m02 =(unifrnd(l2,u2,1,runs));
m03 =(unifrnd(l3,u3,1,runs));
m04 =(unifrnd(l4,u4,1,runs));

m0=[m01' m02' m03' m04'];
m0start=m0(1,:);

pts = m0;
tpoints = CustomStartPointSet(pts);
pts = list(tpoints);

% Provide Options and Create the Multistart Problem
options = optimoptions('lsqcurvefit','Algorithm','trust-region-reflective',...
  'OptimalityTolerance',1e-6,'MaxFunctionEvaluations',5e3,...
   'MaxIterations',5e3,'display','off','FiniteDifferenceType','central',...
 'TolFun',1e-6,'TolX',1e-6,'SpecifyObjectiveGradient',true,'CheckGradients',false);

problem = createOptimProblem('lsqcurvefit','x0',m0start,'objective',fun,...
          'lb',lb,'ub',ub,'xdata',xdata,'ydata',ydata,'options',options);

problem2 = createOptimProblem('lsqcurvefit','x0',m0start,'objective',fun2,...
          'lb',lb,'ub',ub,'xdata',xdata2,'ydata',ydata2,'options',options);
      
ms = MultiStart('StartPointsToRun','all','UseParallel',false,'Display','off');

% Call the Optimization and Store the possible Minima and respective Function Values

tic 
[mbest,errorbest,exitflag,Fcount,soln] = run(ms,problem2,tpoints);
TIME_Fitting=toc
Fcount
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [icell,J]=icellT(m,nin)

m1=m(1);
m2=m(2);
m3=m(3);
m4=m(4);
n=nin;

icell=1.*((5753448039979137964463181271138304.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./35184372088832) - 1))./((3060513257434037./1125899906842624).^((1802042270621195.*n)./140737488355328).*(5244825130281784884355015968292864.*(3060513257434037./1125899906842624).^((1802042270621195.*n)./70368744177664).*(10.^m3...
    + (2.*10.^m1.*10.^m3.*3853.^(1./2).*20265.^(1./2))./20265 + (2.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2))./(20265.*10.^m3)) + 5244825130281784884355015968292864.*10.^m4 + (7764359926397905084167307132928.*10.^m1.*3853.^(1./2).*20265.^(1./2).*(10.^m2 + 10.^m4))./15)));

if nargout>1
J(:,1)=1.*((13247803690171774509280495290233856.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./35184372088832)...
    - 1))./((3060513257434037./1125899906842624).^((1802042270621195.*n)./140737488355328).*(5244825130281784884355015968292864.*(3060513257434037./1125899906842624).^((1802042270621195.*n)./70368744177664).*(10.^m3 ...
    + (2.*10.^m1.*10.^m3.*3853.^(1./2).*20265.^(1./2))./20265 + (2.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2))./(20265.*10.^m3)) + 5244825130281784884355015968292864.*10.^m4 ...
    + (7764359926397905084167307132928.*10.^m1.*3853.^(1./2).*20265.^(1./2).*(10.^m2 + 10.^m4))./15)) ...
    - (4990320691382279.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./70368744177664).*((864160113899737.*10.^m1.*10.^m3.*3853.^(1./2).*20265.^(1./2))./3802726935360962560 ...
    + (864160113899737.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2))./(3802726935360962560.*10.^m3)) + (864160113899737.*10.^m1.*3853.^(1./2).*20265.^(1./2).*(10.^m2...
    + 10.^m4))./3802726935360962560).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./35184372088832) - 1))./(4549160640446464.*(3060513257434037./1125899906842624).^((1802042270621195.*n)./140737488355328).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./70368744177664).*(10.^m3...
    + (2.*10.^m1.*10.^m3.*3853.^(1./2).*20265.^(1./2))./20265 + (2.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2))./(20265.*10.^m3)) + 10.^m4 + (2.*10.^m1.*3853.^(1./2).*20265.^(1./2).*(10.^m2 + 10.^m4))./20265).^2));

J(:,2)=1.*((13247803690171774509280495290233856.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./35184372088832)...
    - 1))./((3060513257434037./1125899906842624).^((1802042270621195.*n)./140737488355328).*(5244825130281784884355015968292864.*(3060513257434037./1125899906842624).^((1802042270621195.*n)./70368744177664).*(10.^m3 ...
    + (2.*10.^m1.*10.^m3.*3853.^(1./2).*20265.^(1./2))./20265 + (2.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2))./(20265.*10.^m3)) + 5244825130281784884355015968292864.*10.^m4 ...
    + (7764359926397905084167307132928.*10.^m1.*3853.^(1./2).*20265.^(1./2).*(10.^m2 + 10.^m4))./15)) - (4990320691382279.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./35184372088832)...
    - 1).*((864160113899737.*10.^m1.*10.^m2.*3853.^(1./2).*20265.^(1./2))./3802726935360962560 + (864160113899737.*(3060513257434037./1125899906842624).^((1802042270621195.*n)./70368744177664).*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2))./(3802726935360962560.*10.^m3))...
    )./(4549160640446464.*(3060513257434037./1125899906842624).^((1802042270621195.*n)./140737488355328).*((3060513257434037./1125899906842624).^(...
    (1802042270621195.*n)./70368744177664).*(10.^m3 + (2.*10.^m1.*10.^m3.*3853.^(1./2).*20265.^(1./2))./20265 + (2.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2))./(20265.*10.^m3))...
    + 10.^m4 + (2.*10.^m1.*3853.^(1./2).*20265.^(1./2).*(10.^m2 + 10.^m4))./20265).^2));

J(:,3)=1.*(-(4990320691382279.*(3060513257434037./1125899906842624).^((1802042270621195.*n)./70368744177664).*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./35184372088832) ...
    - 1).*((2592480341699211.*10.^m3)./1125899906842624 + (864160113899737.*10.^m1.*10.^m3.*3853.^(1./2).*20265.^(1./2))./3802726935360962560 ...
    - (864160113899737.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2))./(3802726935360962560.*10.^m3)))./(4549160640446464.*(3060513257434037./1125899906842624).^(...
    (1802042270621195.*n)./140737488355328).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./70368744177664).*(10.^m3 + (2.*10.^m1.*10.^m3.*3853.^(1./2).*20265.^(1./2))./20265 ...
    + (2.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2))./(20265.*10.^m3)) + 10.^m4 + (2.*10.^m1.*3853.^(1./2).*20265.^(1./2).*(10.^m2 + 10.^m4))./20265).^2));

J(:,4)=1.*((13247803690171774509280495290233856.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./35184372088832) ...
    - 1))./((3060513257434037./1125899906842624).^((1802042270621195.*n)./140737488355328).*(5244825130281784884355015968292864.*(3060513257434037./1125899906842624).^((1802042270621195.*n)./70368744177664).*(10.^m3 ...
    + (2.*10.^m1.*10.^m3.*3853.^(1./2).*20265.^(1./2))./20265 + (2.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2))./(20265.*10.^m3)) ...
    + 5244825130281784884355015968292864.*10.^m4 + (7764359926397905084167307132928.*10.^m1.*3853.^(1./2).*20265.^(1./2).*(10.^m2 + 10.^m4))./15)) ...
    - (4990320691382279.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2).*((3060513257434037./1125899906842624).^((1802042270621195.*n)./35184372088832) ...
    - 1).*((2592480341699211.*10.^m4)./1125899906842624 + (864160113899737.*10.^m1.*10.^m4.*3853.^(1./2).*20265.^(1./2))./3802726935360962560 ...
    + (864160113899737.*(3060513257434037./1125899906842624).^((1802042270621195.*n)./70368744177664).*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2)...
    )./(3802726935360962560.*10.^m3)))./(4549160640446464.*(3060513257434037./1125899906842624).^((1802042270621195.*n)./140737488355328).*((3060513257434037./1125899906842624).^(...
    (1802042270621195.*n)./70368744177664).*(10.^m3 + (2.*10.^m1.*10.^m3.*3853.^(1./2).*20265.^(1./2))./20265 + (2.*10.^m1.*10.^m2.*10.^m4.*3853.^(1./2).*20265.^(1./2)...
    )./(20265.*10.^m3)) + 10.^m4 + (2.*10.^m1.*3853.^(1./2).*20265.^(1./2).*(10.^m2 + 10.^m4))./20265).^2));

end

Answer 2

[Torsten]

From the documentation:

If the Jacobian can also be computed and the Jacobian option is 'on', set by

options = optimoptions('lsqcurvefit','SpecifyObjectiveGradient',true)

then the function fun must return a second output argument with the Jacobian 
value J (a matrix) at x. 
By checking the value of nargout, the function can avoid computing J when fun 
is called with only one output argument (in the case where the optimization 
algorithm only needs the value of F but not J).

function [F,J] = myfun(x,xdata)
F = ...          % objective function values at x
if nargout > 1   % two output arguments
   J = ...   % Jacobian of the function evaluated at x
end


If fun returns a vector (matrix) of m components and x has n elements, 
where n is the number of elements of x0, the Jacobian J is an m-by-n matrix 
where J(i,j) is the partial derivative of F(i) with respect to x(j). 
(The Jacobian J is the transpose of the gradient of F.) 

For more information, see Writing Vector and Matrix Objective Functions.

[Pas]

Dear Torsten thanks for the feedback. I have uploaded an example of my code to a comment under John D'Errico's answer. Any feedback would be appreciatied

[Torsten]

fun2= @(m,xdata2)[fun1(m,xdata2(1:19,1));fun1(m,xdata2(20:end,1))];

instead of

fun2= @(m,xdata)[fun1(m,xdata(1:19,1));fun1(m,xdata(20:end,1))];

[Pas]

Thanks for pointing out this typing error. I corrected this in the above code as well. This doesn't change the outpout error of "Too many output arguments".

The point of what I would like to do is to instruct if possible lsqcurvefit to use the gradient for each function,(in the function vector list), i.e. for example for xdata2(1:19) use objective function fun1 and its corresponding gradient and for xdata2(20:end) do the same with function fun1(in this example it's the same objective function) with it's respective gradient.

If I use deal like this :

fun2= @(m,xdata2)deal([fun1(m,xdata2(1:19,1));fun1(m,xdata2(20:end,1))]);

I get the following error:

"Error using lsqncommon (line 45)

Supplied Jacobian is not the correct size:

the Jacobian matrix should be 38-by-4.

Error in lsqcurvefit (line 257)

lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,caller,...

Error in fmultistart

Error in MultiStart/run (line 268)

fmultistart(problem,startPointSets,msoptions);

Caused by:

Failure in call to the local solver with supplied problem structure."

Is there any way I can supply lsqqcurvefit with the jacobian in separate function?

[Torsten]

Here is a simple code that reproduces your error.

function main
  p(1) = 2.0; 
  p(2) = 3.0;
  xdata(1,1) = 4.0;
  xdata(2,1) = 6.0;  
  fun = @(p,xdata)[fun1(p,xdata(1:1,1));fun1(p,xdata(2:end,1))];
  I = fun(p,xdata)
  size(I)
  [I,J] = fun(p,xdata)
  size(I)
  size(J)
end

function [I,J] = fun1(p,xdata)
 I = xdata;
 if nargout > 1
   J = xdata;
 end  
end 

For one single output argument, MATLAB seems to be able to "vertcat" the output as requested by your function definition in fun2.

I don't know of a way to tell MATLAB that it should do the same with the second output argument if nargout > 1.

Maybe somebody in the forum can help here.

Best wishes

Torsten.

[Pas]

Thank you very much for the Input Torsten

[Torsten]

Of course you can pass a function "fun" to lsqcurvefit as

fun = @(p,xdata2) fun1(p,xdata2(1:19,1),xdata2(20:end,1));

and do the concatenation on your own ...

[Pas]

So if I understand your suggestion correctly I should create a single function which contains all listed sub-objectivefunctions (fun1, fun2 fun3 etc, here is the same function ), as well as corresponding gradients and set for example either an if or switch conditition depended on the xdata index such as to calculate only the corresponding (to xdata(i:j) ) objective (scalar) and jacobian (19,4) for i:j?

[Torsten]

Pass all the information to "fun" needed to build I and J for your specific purpose.

In the case above, the call to "driver" could be e.g.

fun = @(p,xdata) driver(p,xdata(1:19,1),fun1,xdata(20:end),fun2)

which would indicate that xdata(1:19,1) is to be evaluated with fun1 (a function handle) and xdata(20:end) is to be evaluated with fun2 (maybe a different function handle).

[Pas]

Thank you very much Torsten, I will try it.