Improving LSQCcurvefit solutions through improving x0 through loops

Hello everyone,

I had a quick question in regards to attempting to improve my lsqccurvefit values by improving x0 through multiple datasets.

Data1=[0.596421471	0.06
0.5859375	0.119284294
0.566037736	0.29296875
0.530035336	0.622641509
0.418994413	1.219081272
0.388619855	3.058659218
];
Data2=[5.00E+04	3.82E+04	3.45E+04	3.42E+04	3.74E+04	3.21E+04	2.81E+04	2.88E+04
5.00E+04	3.82E+04	3.45E+04	3.42E+04	3.74E+04	3.21E+04	2.81E+04	2.88E+04
3.08E+04	3.07E+04	2.19E+04	2.23E+04	2.53E+04	2.05E+04	1.98E+04	1.89E+04
2.05E+04	1.87E+04	1.30E+04	1.10E+04	1.62E+04	1.31E+04	1.05E+04	1.05E+04
8963	1.11E+04	6243	3504	6454	4346	4448	4357
0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00
];
Pt=Data1(:,1);
Lt=Data1(:,2);
for n = 1:size( Data2, 2 )
   Int=Data2(:,n);
  ...

plot(Lt,Int, 'ro');
xdata=[Pt,Lt]
F=@(x,xdata)((xdata(:,1)+xdata(:,2)+x(2))-((xdata(:,1)+xdata(:,2)+x(2)).^2-4.*xdata(:,1).*xdata(:,2).^1/2).*x(1)/2.*xdata(:,1))
x0=[1 1]
options=optimoptions('lsqcurvefit','MaxFunEvals',5e2)
lb=[];
up=[]
[x,resnorm,~,exitflag,output]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
lb = [];
ub = [];
[x2,resnorm2,~,exitflag2,output2]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)
M=mean(x) 
M=mean(x2)
hold on 
plot (Lt, F(x,xdata))
plot (Lt, F(x2,xdata))
hold off
end

In a nutshell, I have a function F that is trying to fit my Data 2 using 2 different algorithms. My end plot currently has a poor fit with my initial plot. I've seen that the initial starting values [x0] can have a signifigant impact on your results. Is there a way to have the later iterations in my loop change the x0 value and sample various x0 values and try and get closer to the real value that fits my data well (this can be done by doing an RMSD or something like that with my initial plot and vlues). The simplest way I was thinking of having this done is using having x0 change to the results of the previous run. I.E. The first run using the x0 I gave it, the 2nd x0 that is used in the 2nd run uses the results of the previous run as it's own x0.

I don't quite know how to do this, whether this is a modification through x0 (like my example), or through using a loop within my current loop that runs the same dataset n=1 repeatedly through multiple iterations using the RMSD example I gave above to try and improve it. Then do the same thing for the 2nd dataset (n=2), although I can see this method taking quite a while and maybe potentially creating an endless loop.

Answer 1

[Matt J]

For such a small problem, it is a simple matter to do a vectorized exhaustive search, rather than looping

[X1,X2]=meshgrid(linspace(-50000,0,500), linspace(-60,60,100));
[mm,nn]=size(X1);
XX1=reshape(X1,1,1,[]);
XX2=reshape(X2,1,1,[]);

Y=reshape(vecnorm(F([XX1,XX2],xdata)-Int,inf,1),mm,nn);

[minval,idx]=min(Y(:)),
x0=[X1(idx),X2(idx)]

[Sam Mahdi]

Thank you for the quick answer!

Could you point me in the proper direction for learning what a vectorized exhaustive search is? And why it would be better than looping? I'm a little bit unfamiliar with both matlab and the mathematical functions/operations.

My main source has just been using doc "item", but there doesn't appear to be anything on exhaustive searches for that.

[Matt J]

Exhaustive search just means to search for the best x by doing a sweep over all parameter combinations. In the code I've shown you the objective function is calculated at all lattice points on some 2D grid search grid (except I do it without using explicit loops - instead I just use Matlab's matrix arithmetic) and I can just take the minimum over all the sample points. Exhaustive is a bit of a misnomer here - it's not really exhaustive because you have a continuum of feasible solutions, and this search sweeps through a discretization of that. Also, you still have to guess the grid box size big enough to contain the optimal solution. However, the search can give you an x0 as close to the global optimum as the fineness of the search grid allows, and also tells you approximately the best fitting accuracy (minval) that you can achieve.

You can plot the objective function at all the search points by adding a few lines to the code,

[X1,X2]=meshgrid(linspace(-1e5,1e5,100), linspace(-600,600,100));
[mm,nn]=size(X1);
XX1=reshape(X1,1,1,[]);
XX2=reshape(X2,1,1,[]);

Y=reshape(vecnorm(F([XX1,XX2],xdata)-Int,inf,1),mm,nn);

[minval,idx]=min(Y(:)),
x0=[X1(idx),X2(idx)]

surf(X1,X2,Y);scf
xlabel x1, ylabel x2

[Sam Mahdi]

Thank you for the explanation. So from my understanding, your code defines x0 as being 2 variables that are searching a grid with the size you provided for the local minimum. So I included this into my loop, since I'd desirably want this done to every single datapoint that I have

for n = 1:size( Data2, 2 )
   Int=Data2(:,n);
  ...
hold on
plot(Lt,Int, 'ro');
xdata=[Pt,Lt]
F=@(x,xdata)((xdata(:,1)+xdata(:,2)+x(2))-((xdata(:,1)+xdata(:,2)+x(2)).^2-4.*xdata(:,1).*xdata(:,2).^1/2).*x(1)/2.*xdata(:,1))
options=optimoptions('lsqcurvefit','MaxFunEvals',5e2)
lb=[];
up=[]
[x,resnorm,~,exitflag,output]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
lb = [];
ub = [];
[x2,resnorm2,~,exitflag2,output2]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)
[x1,x2]=meshgrid(linspace(-50000,0,500), linspace(-60,60,100));
[mm,nn]=size(x1);
xx1=reshape(x1,1,1,[]);
xx2=reshape(x2,1,1,[]);
Y=reshape(vecnorm(F([xx1,xx2],xdata)-Int,inf,1),mm,nn);
[minval,idx]=min(Y(:)),
x0=[x1(idx),x2(idx)]
M=mean(x) 
M=mean(x2) 
plot (Lt, F(x,xdata))
plot (Lt, F(x2,xdata))
hold off
end

And I got this error.

Error using reshape
To RESHAPE the number of elements must not
change.

Error in Testscript4loop (line 36)
Y=reshape(vecnorm(F([xx1,xx2],xdata)-Int,inf,1),mm,nn); 
>> 

Initially you had everything as capital x's, so I changed them all to lower case (since in my function, it was a lower case x, although I don't know if it's cap sensitive or not. Then I thought it might be an issue since you have it written as x1 and x2, but I don't technically have an x1 in y script (unless x is automatically assumed to be x1 by the program), but I still got the exact same error. I did the reverse as well, changing all the x1s to just x, and that also did not fix it. So it doesn't appear to be an issue regarding syntax.

[Matt J]

You should have the following:

for n = 1:size( Data2, 2 )
   Int=Data2(:,n);
  ...
hold on
plot(Lt,Int, 'ro');
xdata=[Pt,Lt]

F=@(x,xdata)((xdata(:,1)+xdata(:,2)+x(2))-((xdata(:,1)+xdata(1,2,:)+x(2)).^2-4.*xdata(:,1).*xdata(:,2).^1/2).*x(1,1,:)/2.*xdata(:,1))

[x1,x2]=meshgrid(linspace(-50000,0,500), linspace(-60,60,100));
[mm,nn]=size(x1);
xx1=reshape(x1,1,1,[]);
xx2=reshape(x2,1,1,[]);
Y=reshape(vecnorm(F([xx1,xx2],xdata)-Int,inf,1),mm,nn);
[minval,idx]=min(Y(:)),
x0=[x1(idx),x2(idx)]


options=optimoptions('lsqcurvefit','MaxFunEvals',5e2)
lb=[];
up=[]
[x,resnorm,~,exitflag,output]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
lb = [];
ub = [];
[x2,resnorm2,~,exitflag2,output2]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)

M=mean(x) 
M=mean(x2) 
plot (Lt, F(x,xdata))
plot (Lt, F(x2,xdata))
hold off
end

[Sam Mahdi]

Thank you, this definitely helped my function fit my data better. However, a few questions again.

1) I noticed the main difference here was the change in the my function (so xdata 2nd column [xdata :,2] is now xdata 1,2, ) and my x(1) variable is now x(1,1,:). I'm curious what this is doing exactly, and why you only changed one of the xdata (:,2) and not the other ones. These changes seem to make distinct changes to my output, but I don't exactly understand what's going on. As well as why xdata (:,1) isn't touched.

F=@(x,xdata)((xdata(:,1)+xdata(:,2) [why this isn't changed]+x(2))-((xdata(:,1)+xdata(1,2,:)[but this is changed]+x(2)).^2-4.*xdata(:,1).*xdata(:,2).^1/2).*x(1,1,:)/2.*xdata(:,1))

2) While my fit is much better now, my results are still very poor. The issue with this is I get negative values; however these results are actually indicating concentration, so negative values don't actually make sense. I have attempted to resolve this by inputing a lower bound of 0

lb=[0];

for my functions but when I do, I get this error.

Warning: Length of lower bounds is <
length(x); filling in missing lower bounds
with -Inf. 

And I can't use lower bounds for the levenberg algorithm since it doesn't seem to accept it (which is a shame, since this is the one that's actually fitting my data decently). I believe this is the main other issue I'm now at, simply that the function should approach zero (and this will give it the nice logarithmic curve), but as it stands, the function goes straight through zero, thus giving me the negative values for x.

[Matt J]

I'm curious what this is doing exactly, and why you only changed one of the xdata (:,2) and not the other ones.

I made a mistake. The redefinition should look like

F=@(x,xdata)((xdata(:,1)+xdata(:,2)+x(1,2,:))-((xdata(:,1)+xdata(:,2)+x(1,2,:)).^2-4.*xdata(:,1).*xdata(:,2).^1/2).*x(1,1,:)/2.*xdata(:,1));

In other words, only the indexing of x should change. This makes it so that you can evaluate a stack of pairs [x1,x2] of dimension 1x2xN simultaneously.

I have attempted to resolve this by inputing a lower bound of 0

You have 2 unknowns, so it should be lb=[0,0]. You should also choose a search grid that excludes negative values, perhaps something like

[x1,x2]=meshgrid(linspace(0,1e4,500), linspace(0,600,100));

While my fit is much better now, my results are still very poor.

That doesn't necessarily indicate a failure in lsqcurvefit. Your model equation and your data could simply be badly matched, i.e., this might be the best fit you can hope for. The minval variable computed by my code will give you the approximate magnitude of the fitting errors you should expect, assuming the search grid is chosen large enough to contain the optimal solution.

It might also be worth mentioning that the surface plot that I presented to you for your first data set shows a very suspicious shape. It is flat looking over a very large region. Essentially, every [x1,x2] in that region looks like it would give nearly the same curve. This makes me suspicious even moreso about the validity of the model equation.

And I can't use lower bounds for the levenberg algorithm since it doesn't seem to accept it (which is a shame, since this is the one that's actually fitting my data decently)

Your other technique limits the number of iterations greatly. This is usually a bad idea - it stopes the solver sooner than it wants to and may be the reason the resulting fit is poorer.

[Sam Mahdi]

With these modifications, now my first model (x) just gives a flat line as an output, and my 2nd model gives a decent curve, but even limiting the available starting solutions to only be positive (by limited the grid to only positive values), I still am getting negative values for x2. And the plot is also dipping in the negatives as well. Also, one other thing I noticed, why are the linear grids different for x1 and x2 (by 10 orders of magnitude).

[Matt J]

I still am getting negative values for x2.

Did you correct the lb input? In any case, I would need to see the code as you are running it now.

Also, one other thing I noticed, why are the linear grids different for x1 and x2 (by 10 orders of magnitude).

You can adjust the search grid as you please, but experimentation with the surf plots showed the minima at points where x2 was much closer to zero than x1, so I chose the limits based on that.

[Sam Mahdi]

Here is my script so far.

Data1=[0.596421471	0.06
0.5859375	0.119284294
0.566037736	0.29296875
0.530035336	0.622641509
0.418994413	1.219081272
0.388619855	3.058659218
];
Data2=[5.00E+04	3.82E+04	3.45E+04	3.42E+04	3.74E+04	3.21E+04	2.81E+04	2.88E+04
5.00E+04	3.82E+04	3.45E+04	3.42E+04	3.74E+04	3.21E+04	2.81E+04	2.88E+04
3.08E+04	3.07E+04	2.19E+04	2.23E+04	2.53E+04	2.05E+04	1.98E+04	1.89E+04
2.05E+04	1.87E+04	1.30E+04	1.10E+04	1.62E+04	1.31E+04	1.05E+04	1.05E+04
8963	1.11E+04	6243	3504	6454	4346	4448	4357
0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00	0.00E+00
];
Pt=Data1(:,1);
Lt=Data1(:,2);

for n = 1:size( Data2, 2 )
   Int=Data2(:,n);
  ...
hold on
plot(Lt,Int, 'ro');
xdata=[Pt,Lt]
F=@(x,xdata)((xdata(:,1)+xdata(:,2)+x(1,2,:))-((xdata(:,1)+xdata(:,2)+x(1,2,:)).^2-4.*xdata(:,1).*xdata(:,2).^1/2).*x(1,1,:)/2.*xdata(:,1))
[x1,x2]=meshgrid(linspace(0,1e4,5e4), linspace(0,60,100));
[mm,nn]=size(x1);
xx1=reshape(x1,1,1,[]);
xx2=reshape(x2,1,1,[]);
Y=reshape(vecnorm(F([xx1,xx2],xdata)-Int,inf,1),mm,nn);
[minval,idx]=min(Y(:)),
x0=[x1(idx),x2(idx)]
options=optimoptions('lsqcurvefit','MaxFunEvals',5e2)
lb=[0,0];
up=[]
[x,resnorm,~,exitflag,output]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
lb = [];
ub = [];
[x2,resnorm2,~,exitflag2,output2]=lsqcurvefit(F,x0,xdata,Int,lb,up,options)
M=mean(x) 
M=mean(x2)
plot (Lt, F(x,xdata))
plot (Lt, F(x2,xdata))
hold off
end

And yes, I do expect the x2 values to be close to zero (looking at hopefully somewhere in the 0.1-0.001 range.

[Matt J]

You are getting negative values because you have lb=[].

[Sam Mahdi]

Yes, but that's for the levenberg equation, and I can't set a lower bound for that. In the graph I posted, the levenberg algorithm are the curved lines, and the trust region algorithms are the flat lines. I can't set lb for levenberg because then I get errors.

[Matt J]

And the plot is also dipping in the negatives as well.

Your model function has a quadratic term -Lt.^2 in it so I can't see why it would be expected to go to 0 as Lt-->Inf. Again, I think your model equation is the problem. It is not in agreement with the data.

[Sam Mahdi]

I'm a bit confused as to what you mean. The equation should approach zero. If we simplify the equation to just xs, we have.

In this case, we can take x2 or x3 (Pt or Lt), and eliminate all other variables for simplicity (assume it's a constant).

x2

As x2 approaches Inf, the 4 and 2 is irrelevant, so those can also be removed. And since the squared version of x2 will be much bigger than 4*x2, we can remove that as well. Leaving us with

/x2

That simplifies to

y=x2-x2 /x2

y=0

This can be applied to any variable. This equation should reach zero if any of these variables approach infinity (outside of x1).

[Matt J]

The equations look fine, but they're not what you've implemented in code.

[Sam Mahdi]

I found my error, I hadn't closed of the interior function, and I had an extra parenthesis at the end of my function. I have thus modified my function :

F=@(x,xdata)((xdata(:,1)+xdata(:,2)+x(1,2,:))-((xdata(:,1)+xdata(:,2)+x(1,2,:)).^2-4.*xdata(:,1).*xdata(:,2).^1/2)).*x(1,1,:)/2.*xdata(:,1)

Which should now properly represent the equations I have posted in my previous post. However, this provides even stranger results. I'm still getting negative values, and now my graph looks funky.

[Matt J]

Note that x.^1/2 is not the square root of x. Compare:

>> 2.^1/2

ans =

     1

>> sqrt(2)

ans =

    1.4142

[Sam Mahdi]

Oh, well that makes a huge difference. With that modification, my graph is no longer in the negative. Although, now it looks even weirder! I should note, I had forgotten to normalize my data prior, it won't change the plot, but does make the values fit better to the quadratic function you have below.

Data2=[0.914458052	0.698775361	0.630597697	0.625479803	0.684335588	0.587461159	0.512703345	0.526777554
0.914458052	0.698775361	0.630597697	0.625479803	0.684335588	0.587461159	0.512703345	0.526777554
0.535527691	0.534656914	0.38174852	0.388714734	0.440961338	0.356844305	0.34430512	0.328631139
0.401724476	0.36664707	0.253968254	0.215755438	0.316872428	0.257103665	0.206153243	0.204977464
0.182397232	0.224867725	0.127045177	0.071306471	0.131339031	0.088441188	0.090516891	0.088665039
0	0	0	0	0	0	0	0
];

Answer 2

[Matt J]

I can see that your model function derives in some way from the quadratic formula, but I can't quite see where all the pieces are supposed to go. In any case, using the implementation below, I get a fit that approaches sensible. Maybe you can tweak it to perfection.

for n = 1:size( Data2, 2 )
   Int=Data2(:,n);
 
xdata=[Pt,Lt]
%F=@(x,xdata)((xdata(:,1)+xdata(:,2)+x(1,2,:))-((xdata(:,1)+xdata(:,2)+x(1,2,:)).^2-4.*xdata(:,1).*xdata(:,2).^1/2).*x(1,1,:)/2.*xdata(:,1))
F=@modelfun;
[X1,X2]=meshgrid(linspace(0,1e4,500), linspace(0,60,100));
[mm,nn]=size(X1);
xx1=reshape(X1,1,1,[]);
xx2=reshape(X2,1,1,[]);
Y=reshape(vecnorm(F([xx1,xx2],xdata)-Int,inf,1),mm,nn);
[minval,idx]=min(Y(:)),
x0=[X1(idx),X2(idx)];

options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
[x2,resnorm2,~,exitflag2,output2]=lsqcurvefit(F,x0,xdata,Int,[],[],options)

plot (Lt,Int,'ro',Lt, F(x2,xdata))
xlabel 'Lt'

end

function out=modelfun(x,xdata)

     Pt=xdata(:,1);
     Lt=xdata(:,2);
     x1=x(1,1,:);
     x2=x(1,2,:);
     
     A=Lt.*x1;
     B=-(Pt+Lt+x2);
     C=Pt;


     out = ( -B + sqrt(B.^2-4.*A.*C) )./(2.*A);
 

end

[Sam Mahdi]

It is the quadratic equation. The only addition is the quadratic is being multiplied by a variable (x1). You have the right set up in your script in terms of what everything stands for, except x1 is in the numerator, and you currently have it in the denominator. This is why I'm confused as to why my graphs look funky, because if you plot it out using it the way you have it written, you get a nice curve. While my function is essentially the same thing, but provides the above funky results.