Simple optimisation question
in [Matlab]
|
From: Matt J on 29 Jun 2010 14:49 "Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <i0da7l$e7j$1(a)fred.mathworks.com>... > 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 > > x=X(X>=data(1:end-1) & X<=data(2:end)) > > 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. ====================== 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)). The idea being that this numerical problem happens when the computed X(i) values are barely distinguishable from one of the data points.
From: Roger Stafford on 29 Jun 2010 15:22 "Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <i0df7h$br3$1(a)fred.mathworks.com>... > 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)). > > The idea being that this numerical problem happens when the computed X(i) values are barely distinguishable from one of the data points. - - - - - - - - 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 round-off error estimates. Also you don't have to make a special case of x = [] beforehand. [ignore,ix] = min((max(data(1:end-1)-X,0)+max(X-data(2:end),0)); x = X(ix) Roger Stafford
From: Matt J on 29 Jun 2010 16:55 "Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <i0dh57$hnq$1(a)fred.mathworks.com>... > - - - - - - - - > 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 round-off error estimates. Also you don't have to make a special case of x = [] beforehand. > > [ignore,ix] = min((max(data(1:end-1)-X,0)+max(X-data(2:end),0)); > x = X(ix) > ============ 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 if R==0 x=max(data); elseif R==inf x=min(data); else % min(data) < x < max(data) data=sort(data(:)); S=sum(data); N=length(data); Xl=data(1:end-1); Sunder=cumsum(Xl); Sover=S-Sunder; Nunder=(1:N-1).'; Nover=N-Nunder; Rgrid=(Sover-Xl.*Nover)./(Xl.*Nunder-Sunder); idx=find(R-Rgrid>0,1); x=(Sover(idx)+R*Sunder(idx))./(Nover(idx)+R*Nunder(idx)); end
From: Roger Stafford on 29 Jun 2010 17:26 "Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <i0dmju$g1k$1(a)fred.mathworks.com>... > 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 > ........ - - - - - - - - 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. Roger Stafford
From: Matt J on 29 Jun 2010 18:10
"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <i0dode$eg3$1(a)fred.mathworks.com>... > 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. ============ Well, just goes to show. You shouldn't give up so quickly... |