Improving LSQCcurvefit solutions through improving x0 through loops
Show older comments
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.
Answers (2)
Matt J
on 12 Mar 2019
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)]
17 Comments
Sam Mahdi
on 12 Mar 2019
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
on 12 Mar 2019
Matt J
on 12 Mar 2019
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
on 13 Mar 2019
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.
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). 
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
on 13 Mar 2019
Matt J
on 13 Mar 2019
You are getting negative values because you have lb=[].
Sam Mahdi
on 13 Mar 2019
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
on 13 Mar 2019
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).
x2As 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
/x2That 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
on 13 Mar 2019
The equations look fine, but they're not what you've implemented in code.
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
on 13 Mar 2019
Note that x.^1/2 is not the square root of x. Compare:
>> 2.^1/2
ans =
1
>> sqrt(2)
ans =
1.4142
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
];
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
1 Comment
Sam Mahdi
on 13 Mar 2019
Categories
Find more on Mathematics and Optimization in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
