http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654
MATLAB Central Newsreader  Simple optimisation question
Feed for thread: Simple optimisation question
enus
©19942014 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Mon, 28 Jun 2010 13:36:23 +0000
Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758130
Frédéric Bergeron
Hi all,<br>
<br>
First, thanks for the answer I got in the past in this forum, hope I'll get some today! :)<br>
<br>
My situation:<br>
I got a set of 1dimensional data (data) and I want to find the scalar x that respect a previously defined ratio R:<br>
<br>
R=sum(data(over)x)/sum(xdata(under));<br>
where:<br>
data(over) represent all the data over x: data(over)=(data>=x);<br>
data(under) represent all the data under x: data(under)=(data<x);<br>
<br>
So R is the ratio of the sum of difference between data over x and x, over the sum of difference between data under x and x.<br>
<br>
My problem is that I need the value of x to define data(over) and data(under). Is there's a function that could help me? I read about lsqlin() that can solve a problem with inequality constraint but I didn't get success at using it.<br>
<br>
Thanks in advance!<br>
Best regards,

Mon, 28 Jun 2010 13:53:12 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758145
Alan Weiss
On 6/28/2010 9:36 AM, Frédéric Bergeron wrote:<br>
> Hi all,<br>
><br>
> First, thanks for the answer I got in the past in this forum, hope I'll<br>
> get some today! :)<br>
><br>
> My situation:<br>
> I got a set of 1dimensional data (data) and I want to find the scalar x<br>
> that respect a previously defined ratio R:<br>
><br>
> R=sum(data(over)x)/sum(xdata(under));<br>
> where:<br>
> data(over) represent all the data over x: data(over)=(data>=x);<br>
> data(under) represent all the data under x: data(under)=(data<x);<br>
><br>
> So R is the ratio of the sum of difference between data over x and x,<br>
> over the sum of difference between data under x and x.<br>
><br>
> My problem is that I need the value of x to define data(over) and<br>
> data(under). Is there's a function that could help me? I read about<br>
> lsqlin() that can solve a problem with inequality constraint but I<br>
> didn't get success at using it.<br>
><br>
> Thanks in advance!<br>
> Best regards,<br>
<br>
I think you would do better to use fzero.<br>
<br>
Alan Weiss<br>
MATLAB mathematical toolbox documentation

Mon, 28 Jun 2010 14:56:23 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758165
Matt J
"Frédéric Bergeron" <frederic.bergeron@logiag.com> wrote in message <i0a8gm$kju$1@fred.mathworks.com>...<br>
<br>
> My problem is that I need the value of x to define data(over) and data(under). Is there's a function that could help me? I read about lsqlin() that can solve a problem with inequality constraint but I didn't get success at using it.<br>
=======<br>
<br>
I could be wrong, but it seems like you can solve this by simple algebra. Let's make the following definitions<br>
<br>
N=length(data); %known <br>
S=sum(data); %known<br>
O=sum(data(over)x)<br>
U=sum(xdata(under))<br>
<br>
Then you have 2 equations that allow you to solve for O and U:<br>
<br>
OR*U=0<br>
O+U=S;<br>
<br>
Once you know O and U, you have the relationship<br>
<br>
OU=SN*x<br>
<br>
which can be used to solve explicitly for x<br>
<br>
x=(S+UO)/N

Mon, 28 Jun 2010 15:40:26 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758185
Matt J
"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i0ad6n$91i$1@fred.mathworks.com>...<br>
<br>
> <br>
> I could be wrong, but it seems like you can solve this by simple algebra. Let's make the following definitions<br>
> <br>
> N=length(data); %known <br>
> S=sum(data); %known<br>
> O=sum(data(over)x)<br>
> U=sum(xdata(under))<br>
> <br>
> Then you have 2 equations that allow you to solve for O and U:<br>
> <br>
> OR*U=0<br>
> O+U=S;<br>
========<br>
<br>
Forget it. I pulled the equation O+U=S from my imagination somehow.<br>
<br>
In any case, there are at most N different possibilities for the sets<br>
data(over) and data(under). For each one, x ranges over a known interval separating data(over) from data(under) and in that interval O R*U is a known linear function of x. You could simply loop over each of these N sets and test whether <br>
O(x)R*U(x) has a root in that interval.<br>
<br>

Mon, 28 Jun 2010 16:14:05 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758200
Matt J
"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i0afpa$3is$1@fred.mathworks.com>...<br>
<br>
> In any case, there are at most N different possibilities for the sets<br>
> data(over) and data(under). For each one, x ranges over a known interval separating data(over) from data(under) and in that interval O R*U is a known linear function of x. You could simply loop over each of these N sets and test whether <br>
> O(x)R*U(x) has a root in that interval.<br>
===========<br>
<br>
The following code implements this, assuming that all data points are unique. I leave the generalization as an exercise.<br>
<br>
%data<br>
R=0.75;<br>
data=rand(1000,1);<br>
<br>
<br>
%engine<br>
data=sort(data);<br>
S=sum(data);<br>
N=length(data);<br>
<br>
Sunder=cumsum(data(1:end1));<br>
Sover=SSunder;<br>
Nunder=(1:N1).';<br>
Nover=NNunder;<br>
<br>
<br>
X=(Sover+R*Sunder)./(Nover+R*Nunder);<br>
<br>
x=X(X>=data(1:end1) & X<=data(2:end))<br>
<br>
%Check<br>
over=(data>=x);<br>
under=~over;<br>
Rcheck=sum(data(over)x)/sum(xdata(under))

Mon, 28 Jun 2010 16:25:30 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758205
Frédéric Bergeron
"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i0afpa$3is$1@fred.mathworks.com>...<br>
> "Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i0ad6n$91i$1@fred.mathworks.com>...<br>
> <br>
> > <br>
> > I could be wrong, but it seems like you can solve this by simple algebra. Let's make the following definitions<br>
> > <br>
> > N=length(data); %known <br>
> > S=sum(data); %known<br>
> > O=sum(data(over)x)<br>
> > U=sum(xdata(under))<br>
> > <br>
> > Then you have 2 equations that allow you to solve for O and U:<br>
> > <br>
> > OR*U=0<br>
> > O+U=S;<br>
> ========<br>
> <br>
> Forget it. I pulled the equation O+U=S from my imagination somehow.<br>
> <br>
> In any case, there are at most N different possibilities for the sets<br>
> data(over) and data(under). For each one, x ranges over a known interval separating data(over) from data(under) and in that interval O R*U is a known linear function of x. You could simply loop over each of these N sets and test whether <br>
> O(x)R*U(x) has a root in that interval.<br>
<br>
<br>
Hey thank you,<br>
<br>
<br>
The O+U=S equation is fine for me, I just don't get the OU=SN*x.<br>
Anyway, I decided to do a loop for a lot of values in the interval [min(data) max(data)], than calulate the ratio R for each of them and then take the value associated with the ratio R the closest than the previously defined one.<br>
<br>
So I started by making a for loop, than I tought a while loop will be better and more efficient to locate the value with the least error on R. I redefine the upper or lower value of the interval for each while loop, and it's way faster. Here's my code if someone interested:<br>
<br>
*****************************<br>
clear all;<br>
<br>
%Input<br>
R=2; %Cut/fill ratio<br>
datamin=55;<br>
datamax=65;<br>
data=datamin+(datamaxdatamin)*rand(1,50);<br>
<br>
%% First method<br>
tic;<br>
x=min(data):(max(data)min(data))/200:max(data);<br>
for i=1:length(x)<br>
y(i)=abs(R*sum(x(i)data(data<x(i)))sum(data(data>=x(i))x(i)));<br>
end<br>
x_best1=x(y==min(y));<br>
R2_1=sum(data(data>=x_best1)x_best1)/sum(x_best1data(data<x_best1));<br>
erreur_ratio1=abs(R2_1R)/R;<br>
t1=toc;<br>
<br>
%% Second method<br>
tic;<br>
erreur_ratio2=1; %error on R<br>
x1=min(data); %Lower value of the interval around the x_best<br>
x2=max(data); %Upper value<br>
<br>
while erreur_ratio2>0.01<br>
<br>
%Calculs<br>
x_best2=x1+(x2x1)/2; %Taking the middle value of the interval<br>
<br>
R2_2=sum(data(data>=x_best2)x_best2)/sum(x_best2data(data<x_best2));<br>
if R2_2==R;<br>
break %If the good ratio is found<br>
elseif R2_2<R<br>
x2=x_best2;<br>
else %R2>R<br>
x1=x_best2;<br>
end<br>
erreur_ratio2=abs(R2_2R)/R;<br>
end<br>
t2=toc;<br>
<br>
%% Methods comparison<br>
fprintf('\nt1= %g\nerreur_ratio1= %g\nx_best1= %g\n',t1,erreur_ratio1,x_best1);<br>
fprintf('\nt2= %g\nerreur_ratio2= %g\nx_best2= %g\n\n',t2,erreur_ratio2,x_best2);<br>
<br>
****************<br>
and I get:<br>
<br>
t1= 0.030075<br>
erreur_ratio1= 0.0113545<br>
x_best1= 59.3732<br>
<br>
t2= 0.0003564<br>
erreur_ratio2= 0.00205718<br>
x_best2= 59.3845<br>
<br>
Thanks again for the answers!!<br>
Fred

Mon, 28 Jun 2010 16:37:05 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758214
Frédéric Bergeron
"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <br>
> The following code implements this, assuming that all data points are unique. I leave the generalization as an exercise.<br>
> ...<br>
> %data<br>
> R=0.75;<br>
> data=rand(1000,1);<br>
> <br>
> <br>
> %engine<br>
> data=sort(data);<br>
> S=sum(data);<br>
> N=length(data);<br>
> <br>
> Sunder=cumsum(data(1:end1));<br>
> Sover=SSunder;<br>
> Nunder=(1:N1).';<br>
> Nover=NNunder;<br>
> <br>
> <br>
> X=(Sover+R*Sunder)./(Nover+R*Nunder);<br>
> <br>
> x=X(X>=data(1:end1) & X<=data(2:end))<br>
> <br>
> %Check<br>
> over=(data>=x);<br>
> under=~over;<br>
> Rcheck=sum(data(over)x)/sum(xdata(under))<br>
<br>
=========<br>
<br>
Hey,<br>
<br>
Thanks for this code, but since I got a working code in my last message, I think I'll use the code with the while loop. It's maybe longer than your code but since I don't understand all the things you do, I'll use my code, while referring to yours if I get any problem. Thanks!

Mon, 28 Jun 2010 17:11:08 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758236
Matt J
"Frédéric Bergeron" <frederic.bergeron@logiag.com> wrote in message <i0aj3h$e0m$1@fred.mathworks.com>...<br>
<br>
> Thanks for this code, but since I got a working code in my last message, I think I'll use the code with the while loop. It's maybe longer than your code but since I don't understand all the things you do, I'll use my code, while referring to yours if I get any problem. Thanks!<br>
==============<br>
<br>
That's fine, as long as you only need a modestly accurate result. Just for your reference though, I reran your test, but substituting your Method #1 with my closedform method. As you can see, the closedform method accuracy achieves machine precision accuracy in approximately the same amount of time.<br>
<br>
%Closed form<br>
t1= 0.000141917<br>
erreur_ratio1= 2.96059e016<br>
x_best1= 0.524563<br>
<br>
%Iterative approx.<br>
t2= 0.000346971<br>
erreur_ratio2= 0.00396906<br>
x_best2= 0.52504

Mon, 28 Jun 2010 18:33:05 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758276
Frédéric Bergeron
> That's fine, as long as you only need a modestly accurate result. Just for your reference though, I reran your test, but substituting your Method #1 with my closedform method. As you can see, the closedform method accuracy achieves machine precision accuracy in approximately the same amount of time.<br>
> <br>
> %Closed form<br>
> t1= 0.000141917<br>
> erreur_ratio1= 2.96059e016<br>
> x_best1= 0.524563<br>
> <br>
> %Iterative approx.<br>
> t2= 0.000346971<br>
> erreur_ratio2= 0.00396906<br>
> x_best2= 0.52504<br>
<br>
Cool thanks again... Yes I only need modestly accurate result. This was for a GUI leveling agricultural fields. So I have just finished transforming my program in 3D and instead of finding a single scalar, it finds the constant of the plane equation, i.e. c in z=ax+by+c, that respect the ratio R.

Tue, 29 Jun 2010 17:24:05 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758594
Roger Stafford
"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i0ahod$f3u$1@fred.mathworks.com>...<br>
> The following code implements this, assuming that all data points are unique. I leave the generalization as an exercise.<br>
> <br>
> R=0.75;<br>
> data=rand(1000,1);<br>
> %engine<br>
> data=sort(data);<br>
> S=sum(data);<br>
> N=length(data);<br>
> Sunder=cumsum(data(1:end1));<br>
> Sover=SSunder;<br>
> Nunder=(1:N1).';<br>
> Nover=NNunder;<br>
> X=(Sover+R*Sunder)./(Nover+R*Nunder);<br>
> x=X(X>=data(1:end1) & X<=data(2:end))<br>
         <br>
Hello Matt J. I want to congratulate you on that solution you came up with so quickly. I worked for quite a while on the problem before looking at yours and came up with something that was based on the same principles but a lot more complicated. That is, I directly solved for all the R values that correspond to boundaries between the discrete (and sorted) data values, and then checked to see which pair of these the given R lay in. Then I solved for x using the appropriate Sunder/Sover pair. Your way is more efficient.<br>
<br>
However, I notice that both our methods are subject to a possible error which, with an unfortunate choice of R, can actually occur. Because of round off errors it is possible that the x you obtain in<br>
<br>
x=X(X>=data(1:end1) & X<=data(2:end))<br>
<br>
can come up with either two answers which are essentially the same or no answer. I have checked and found this to happen when the given R value is very near one of the R boundary values mentioned above. The case of two values is easily dealt with by just choosing the first one since they are essentially equal. However when no answer is provided, that requires a little more effort. Perhaps the two inequalities above can be slightly widened by a sufficient "epsilon" amount based on possible round off errors to always provide at least one value. Or possibly using a more robust method of defining an "interval distance" function which is zero within an interval but equal to the distance from the closer end otherwise, and then choosing the nearest interval with the min function.<br>
<br>
Hello Fred. My advice to you is to make use of Matt's solution. It is a very efficient and accurate algorithm. The rare difficulty I have discussed here is quite easily overcome.<br>
<br>
Roger Stafford

Tue, 29 Jun 2010 18:49:21 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758626
Matt J
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i0da7l$e7j$1@fred.mathworks.com>...<br>
<br>
> However, I notice that both our methods are subject to a possible error which, with an unfortunate choice of R, can actually occur. Because of round off errors it is possible that the x you obtain in<br>
> <br>
> x=X(X>=data(1:end1) & X<=data(2:end))<br>
> <br>
> can come up with either two answers which are essentially the same or no answer. I have checked and found this to happen when the given R value is very near one of the R boundary values mentioned above. The case of two values is easily dealt with by just choosing the first one since they are essentially equal. However when no answer is provided, that requires a little more effort. Perhaps the two inequalities above can be slightly widened by a sufficient "epsilon" amount based on possible round off errors to always provide at least one value. Or possibly using a more robust method of defining an "interval distance" function which is zero within an interval but equal to the distance from the closer end otherwise, and then choosing the nearest interval with the min function.<br>
======================<br>
<br>
Roger  I see what you mean. In the case where x=[], my solution would probably be to assign x to the X(i) that minimizes norm(X(i)data(i)).<br>
<br>
The idea being that this numerical problem happens when the computed X(i) values are barely distinguishable from one of the data points.

Tue, 29 Jun 2010 19:22:15 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758633
Roger Stafford
"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i0df7h$br3$1@fred.mathworks.com>...<br>
> Roger  I see what you mean. In the case where x=[], my solution would probably be to assign x to the X(i) that minimizes norm(X(i)data(i)).<br>
> <br>
> The idea being that this numerical problem happens when the computed X(i) values are barely distinguishable from one of the data points. <br>
       <br>
How about this for a fix, Matt? It seems the most robust. There is no way it can obtain other than a single value and we don't have to go through messy roundoff error estimates. Also you don't have to make a special case of x = [] beforehand.<br>
<br>
[ignore,ix] = min((max(data(1:end1)X,0)+max(Xdata(2:end),0));<br>
x = X(ix)<br>
<br>
Roger Stafford

Tue, 29 Jun 2010 20:55:26 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758656
Matt J
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i0dh57$hnq$1@fred.mathworks.com>...<br>
<br>
>        <br>
> How about this for a fix, Matt? It seems the most robust. There is no way it can obtain other than a single value and we don't have to go through messy roundoff error estimates. Also you don't have to make a special case of x = [] beforehand.<br>
> <br>
> [ignore,ix] = min((max(data(1:end1)X,0)+max(Xdata(2:end),0));<br>
> x = X(ix)<br>
> <br>
============<br>
<br>
I guess that'll do it. I'm assuming the approach you started out with was something like the following? It works with the boundary values of R like you were describing. If so, I'm not sure why you abandoned it, since it seems a lot faster and also numerically well behaved as long as R is not too close to 0<br>
<br>
<br>
if R==0<br>
<br>
x=max(data); <br>
<br>
elseif R==inf<br>
<br>
x=min(data);<br>
<br>
else % min(data) < x < max(data)<br>
<br>
data=sort(data(:));<br>
S=sum(data);<br>
N=length(data);<br>
<br>
Xl=data(1:end1);<br>
<br>
Sunder=cumsum(Xl);<br>
Sover=SSunder;<br>
Nunder=(1:N1).';<br>
Nover=NNunder;<br>
<br>
<br>
Rgrid=(SoverXl.*Nover)./(Xl.*NunderSunder);<br>
<br>
idx=find(RRgrid>0,1);<br>
<br>
<br>
x=(Sover(idx)+R*Sunder(idx))./(Nover(idx)+R*Nunder(idx));<br>
<br>
end

Tue, 29 Jun 2010 21:26:06 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758665
Roger Stafford
"Matt J " <mattjacREMOVE@THISieee.spam> wrote in message <i0dmju$g1k$1@fred.mathworks.com>...<br>
> I guess that'll do it. I'm assuming the approach you started out with was something like the following? It works with the boundary values of R like you were describing. If so, I'm not sure why you abandoned it, since it seems a lot faster and also numerically well behaved as long as R is not too close to 0<br>
> ........<br>
       <br>
Well, it was something like that except that I fooled around with transforming R to W by w = 1/(R+1) which maps R onto [0,1], thinking that might simplify things. As it turned out it was no particular advantage, and about that time I saw your solution and said, "Aw heck, Matt J has beaten me to it again, and stopped work." So it's all your fault for being so quick  no, I'm kidding.<br>
<br>
Roger Stafford

Tue, 29 Jun 2010 22:10:23 +0000
Re: Simple optimisation question
http://www.mathworks.com/matlabcentral/newsreader/view_thread/285654#758675
Matt J
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <i0dode$eg3$1@fred.mathworks.com>...<br>
<br>
> Well, it was something like that except that I fooled around with transforming R to W by w = 1/(R+1) which maps R onto [0,1], thinking that might simplify things. As it turned out it was no particular advantage, and about that time I saw your solution and said, "Aw heck, Matt J has beaten me to it again, and stopped work." So it's all your fault for being so quick  no, I'm kidding.<br>
============<br>
<br>
Well, just goes to show. You shouldn't give up so quickly...