Re: Re: FindMinimum and the minimum-radius circle
- To: mathgroup at smc.vnet.net
- Subject: [mg50357] Re: [mg50348] Re: FindMinimum and the minimum-radius circle
- From: DrBob <drbob at bigfoot.com>
- Date: Sun, 29 Aug 2004 03:54:05 -0400 (EDT)
- References: <cfuq6g$653$1@smc.vnet.net> <200408180834.EAA08732@smc.vnet.net> <cg204q$msq$1@smc.vnet.net> <200408200857.EAA12448@smc.vnet.net> <200408210704.DAA24766@smc.vnet.net> <cg97e9$a5p$1@smc.vnet.net> <200408261051.GAA16446@smc.vnet.net> <cgmmqi$ts$1@smc.vnet.net> <200408280838.EAA19082@smc.vnet.net>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
For both data sets, I can plainly SEE that the answer is wrong on a plot, as the circle clearly doesn't contain three of the data points and clearly doesn't contain two of them on a diameter. I'm not talking about sixth-decimal-place errors. The other methods that survived this far don't make visible errors, and Fred Simons' method is also faster. Bobby On Sat, 28 Aug 2004 04:38:01 -0400 (EDT), Maxim <ab_def at prontomail.com> wrote: > I just have a small objection to using the word 'fails'; after all, > we're talking about numerical optimization. For example, here's what > NMinimize gives for my first example: > > In[1]:= > SeedRandom[1]; > points = Partition[Table[Random[], {1000}], 2]; > > goal[x_, y_] = Max[Norm[{x, y} - #]& /@ points]; > > grad = Compile[ > {x, y, {points, _Real, 2}}, > Module[ > {n = Length@points, Ld = Norm[{x, y} - #]& /@ points, > Ldd, i1, i2}, > {i2, i1} = Ordering[Ld, -2]; > Ldd = ({x, y} - points[[#]])/Ld[[#]]& /@ {i1, i2}; > If[Chop[Ld[[i1]] - Ld[[i2]]] != 0, > Ldd[[1]], > (Ldd[[1]] + Ldd[[2]])/2 > ]], > {{_Ordering, _Integer, 1}} > ]; > > NMinimize[goal[x, y], {x, y}] > > Out[5]= > {0.66281538, {x -> 0.48683054, y -> 0.4614583}} > > It turns out that this is not the exact minimum either. Does it mean > that NMinimize fails as well? The advantage of the FindMinimum > approach is that FindMinimum can be called repeatedly to obtain finer > approximations: > > In[6]:= > FixedPoint[ > FindMinimum[goal[x, y], > {x, x /. #[[2]]}, {y, y /. #[[2]]}, > Method -> Gradient, Gradient :> grad[x, y, points]]&, > {Infinity, {x -> Mean[points][[1]], y -> Mean[points][[2]]}}] > > Out[6]= > {0.66281504, {x -> 0.48684176, y -> 0.46146944}} > > If you want speed increase, replace Norm in the definitions of goal > and grad with Sqrt[#.#]&; Norm is slow. Certainly some more > fine-tuning can be done; for example, replace Chop[Ld[[i1]]-Ld[[i2]]] > with Chop[Ld[[i1]]-Ld[[i2]],10^-6], this seems to give better results > in most cases. > > Maxim Rytin > m.r at inbox.ru > > > DrBob <drbob at bigfoot.com> wrote in message news:<cgmmqi$ts$1 at smc.vnet.net>... >> That frequently works, and it's almost as fast as the Fred Simons method. >> >> But it fails on this data: >> >> {{-0.694112,-1.47331},{-2.26344,0.220021},{-0.8867,-1.85288},{0.527656,-2.\ >> 81915},{-1.75879,-0.28898},{ >> 0.728263,1.20724},{-1.34788,-0.799066},{-0.676629,-1.03418},{1.45111,0.\ >> 961224},{0.914209,-0.539339}} >> >> and this data: >> >> {{1.60034, 0.0176019}, {-0.355331, -1.62872}, {-0.2654, -0.000635396}, \ >> {-0.777536, 0.890349}} >> >> In both cases, the result circle is too big and passes through two points when it should be three. >> >> I'm using Mean@data as the starting point (as you did), but I'm using the ConvexHull to define goal[x_,y_] and as the third argument of grad. >> >> Bobby >> >> On Thu, 26 Aug 2004 06:51:34 -0400 (EDT), Maxim <ab_def at prontomail.com> wrote: >> >> > Daniel Lichtblau <danl at wolfram.com> wrote in message news:<cg97e9$a5p$1 at smc.vnet.net>... >> >> Janos D. Pinter wrote: >> >> > Hello mathgroup, >> >> > >> >> > as stated in previous messages, the minimum-radius circle problem is >> >> > convex, but constrained: hence, an unconstrained optimization approach may >> >> > or may not work... >> >> >> >> But it can certainly be made to work on an unconstrained formulation. >> >> This takes some nondefault option settings to work around problems that >> >> arise because such formulation is not everywhere smooth. >> >> >> >> >> >> > Using the convex hull as a constraint set will work, of >> >> > course. NMInimize and other constrained optimizers should then also work, >> >> > even in local search mode if they have that option. >> >> > >> >> > Instead of further philosophical discussions on the subject, why don't we >> >> > solve a standardized test model(s), using the same numerical example. (I >> >> > assume here that SeedRandom is platform and version independent... as it >> >> > should be.) >> >> > >> >> > Here is a simple model-instance, for your perusal and experiments. >> >> > >> >> > SeedRandom[1]; >> >> > points = Partition[Table[Random[], {1000}], 2]; >> >> > constraints = (#[[1]] - x0)^2 + (#[[2]] - y0)^2 <= r^2 & /@ points; >> >> > soln = NMinimize[{r, constraints}, {x0, y0, {r, 0, 10}}] >> >> > {0.6628150350002375, {r -> 0.6628150350002375, x0 -> 0.48684175822621106, >> >> > y0 -> 0.46146943870460916}} >> >> > >> >> > timedsoln = NMinimize[{r, constraints}, {x0, y0, {r, 0, 10}}] // AbsoluteTiming >> >> > {30.694136`8.938600406559628*Second, {0.6628150350002375, >> >> > {r -> 0.6628150350002375, x0 -> 0.48684175822621106, y0 -> >> >> > 0.46146943870460916}}} >> >> > >> >> > That is (on a P4 1.6 WinXP machine) it takes about 30 seconds to generate >> >> > this solution. >> >> > >> >> > Next, I solve the same model using the MathOptimizer Professional package, >> >> > in local search mode: >> >> > >> >> > Needs["MathOptimizerPro`callLGO`"] >> >> > callLGO[r, constraints, {{x0, 0, 1}, {y0, 0, 1}, {r, 0, 1}}, Method -> LS] >> >> > // AbsoluteTiming >> >> > {6.158856`8.24104504348061*Second, >> >> > {0.662815036, {x0 -> 0.4868417589, y0 -> 0.4614694394, r -> 0.662815036}, >> >> > 4.770200900949817*^-11}} >> >> > >> >> > The two solutions are fairly close (the solution methods are different). >> >> > The MathOptimizer Professional (LGO) solver time is ~ 6 seconds. The result >> >> > also shows that the max. constraint error of this solution is ~ >> >> > 4.770200900949817*^-11. (This error could be reduced, if necessary by >> >> > changing the model and some LGO options.) >> >> > >> >> > Obviously, one could use also more sophisticated models and other point >> >> > sets etc., as long as we all use the same one(s), for objective >> >> > comparisons. (The proof of the pudding principle.) >> >> > >> >> > Regards, >> >> > Janos >> >> > [...] >> >> >> >> What about the unconstrained formulation? >> >> >> >> In[47]:= >> >> Timing[NMinimize[Sqrt[Max[Map[(x-#[[1]])^2+(y-#[[2]])^2&,points]]],{x,y}]] >> >> >> >> Out[47]= >> >> {3.54 Second,{0.662815,{x\[Rule]0.486831,y\[Rule]0.461458}}} >> >> >> >> I note that this can be made faster still by preprocessing to extract >> >> the convex hull of the data points. Also one can attain a faster result, >> >> with but slight loss of accuracy, using >> >> >> >> FindMinimum[...,Gradient->{"FiniteDifference","DifferenceOrder" -> 2}] >> >> >> >> You can do this with arbitrary reasonable starting points but clearly >> >> the mean or median values would be good choices. >> >> >> >> In[52]:= >> >> Timing[FindMinimum[ >> >> Sqrt[Max[ >> >> Map[(x-#[[1]])^2+(y-#[[2]])^2&,points]]],{x,#[[1]]},{y,#[[2]]}, >> >> Gradient->{"FiniteDifference","DifferenceOrder"->2}]&[ >> >> Mean[points]]] >> >> >> >> Out[52]= >> >> {1.34 Second,{0.663207,{x->0.487394,y->0.461453}}} >> >> >> >> Experimentation reveals that by taking DifferenceOrder->4 I can recover >> >> the NMinimize/MOP result for this particular example, albeit at some >> >> loss of speed. >> >> >> >> Again, this will perform much better if you work with the convex hull. >> >> Here is how fast this becomes. >> >> >> >> Needs["DiscreteMath`ComputationalGeometry`"] >> >> >> >> In[69]:= >> >> Timing[With[{hull=points[[ConvexHull[points]]]}, >> >> FindMinimum[ >> >> Sqrt[Max[ >> >> Map[(x-#[[1]])^2+(y-#[[2]])^2&,hull]]],{x,#[[1]]},{y,#[[2]]}, >> >> Gradient->{"FiniteDifference","DifferenceOrder"->2}]&[ >> >> Mean[hull]]]] >> >> >> >> Out[69]= >> >> {0.04 Second,{0.66286,{x->0.486905,y->0.461467}}} >> >> >> >> >> >> My thanks to Bobby Treat, Rob Knapp, and Tom Burton for bringing to my >> >> attention some of the finer points regarding speed and avoidance of >> >> premature convergence by FindMinimum. >> >> >> >> >> >> Daniel Lichtblau >> >> Wolfram Research >> > >> > >> > Another method would be to construct the gradient explicitly; I used a >> > simple heuristic formula >> > >> > (D[Max[f1,...,fn], x] /. x->x0) == (D[f1, x] /. x->x0) /; >> > f1 > f2,...,f1 > fn, >> > (D[Max[f1,...,fn], x] /. x->x0) == (D[(f1+f2)/2, x] /. x->x0) /; >> > f1 == f2, f2 > f3,...,f2 > fn. >> > >> > Here's what I get for the same set of points: >> > >> > In[1]:= >> > SeedRandom[1]; >> > points = Partition[Table[Random[], {1000}], 2]; >> > >> > goal[x_, y_] = Max[Norm[{x, y} - #]& /@ points]; >> > >> > grad = Compile[ >> > {x, y, {points, _Real, 2}}, >> > Module[ >> > {n = Length@points, Ld = Norm[{x, y} - #]& /@ points, >> > Ldd, i1, i2}, >> > {i2, i1} = Ordering[Ld, -2]; >> > Ldd = ({x, y} - points[[#]])/Ld[[#]]& /@ {i1, i2}; >> > If[Chop[Ld[[i1]] - Ld[[i2]]] != 0, >> > Ldd[[1]], >> > (Ldd[[1]] + Ldd[[2]])/2 >> > ]], >> > {{_Ordering, _Integer, 1}} >> > ]; >> > >> > FindMinimum[goal[x, y], >> > {x, Mean[points][[1]]}, {y, Mean[points][[2]]}, >> > Method -> Gradient, Gradient :> grad[x, y, points]] >> > >> > Out[5]= >> > {0.662817, {x -> 0.486845, y -> 0.461469}} >> > >> > This seems to work quite well for any starting point, but might fail >> > if there are triple points where f1[x0]==f2[x0]==f3[x0]. >> > >> > The effect of DifferenceOrder is probably version specific; I tried >> > all combinations with different methods like >> > FindMinimum[...,Method->Gradient,Gradient->{FiniteDifference,DifferenceOrder->4}] >> > but changing DifferenceOrder didn't seem to have any effect. >> > >> > Maxim Rytin >> > m.r at inbox.ru >> > >> > >> > > > > -- DrBob at bigfoot.com www.eclecticdreams.net
- References:
- Re: FindMinimum and the minimum-radius circle
- From: Thomas Burton <tburton@brahea.com>
- Re: Re: FindMinimum and the minimum-radius circle
- From: Thomas Burton <tburton@brahea.com>
- Re: Re: Re: FindMinimum and the minimum-radius circle
- From: "Janos D. Pinter" <jdpinter@hfx.eastlink.ca>
- Re: FindMinimum and the minimum-radius circle
- From: ab_def@prontomail.com (Maxim)
- Re: FindMinimum and the minimum-radius circle
- From: ab_def@prontomail.com (Maxim)
- Re: FindMinimum and the minimum-radius circle