2
$\begingroup$

Define:

sums[x_List] := Module[{n = Length[x]},
  Flatten[Table[x[[i]] + x[[j]], {i, 1, n - 1}, {j, i + 1, n}]]
]

score[x_List] := CentralMoment[Differences[Sort[sums[x]]]/(Max[x] - Min[x]), 2]

f[x_, y_] := score[{0, x, y, 1}]

Then, f looks like this:

Plot3D[f[x, y], {x, 0, 1}, {y, 0, 1}, ImageSize -> Small]

Mathematica graphics

In the interior of the domain shown, the surface has 4 minima (all nicely-behaved):

ContourPlot[f[x, y], {x, 0, 1}, {y, 0, 1}, 
 Contours -> Table[2^-i, {i, 3, 14}],
 MaxRecursion -> 2, 
 PlotPoints -> 50,
 ImageSize -> Small]

Mathematica graphics

But I can't get any of the minimum-finding functions in Mathematica to find any of them. For example:

In[25]:= NMinimize[{f[x, y], 0 < x < 1 && 0 < y < 1}, {x, y}]

Out[25]= {0.514286, {x -> 0.714287, y -> 1.}}

In[26]:= FindMinimum[{f[x, y], 0 < x < 1 && 0 < y < 1}, {x, y}]

Out[26]= {0.514286, {x -> 0.714285, y -> 1.}}

Even when I give Mathematica one of the minima as starting point, it still fails to find it:

In[27]:= f[1/4, 1/2]

Out[27]= 0

In[28]:= FindMinimum[{f[x, y], 0 < x < 1 && 0 < y < 1}, {{x, 1/4}, {y, 1/2}}]

Out[28]= {0.514286, {x -> 0.714285, y -> 1.}}

(Note that score is non-negative, since it's a second central moment. Therefore, any point {x, y} for which f[x, y] equals $0$ is a minimum.)


What am I doing wrong?

$\endgroup$
3
  • 1
    $\begingroup$ Technically, that's "cheating". FindMinimum uses the slope to find the minimum point, so giving it a point with zero slope shouldn't work. $\endgroup$
    – user1722
    May 10, 2016 at 23:54
  • 2
    $\begingroup$ Specifying that the function only works for numbers helps: f[x_?NumberQ, y_?NumberQ] := .... $\endgroup$
    – wxffles
    May 11, 2016 at 0:00
  • $\begingroup$ @wxffles: Brilliant! That took care of it. Thank you. $\endgroup$
    – kjo
    May 11, 2016 at 1:10

2 Answers 2

3
$\begingroup$

In at least this case, one workaround avoiding Sort[] would be to use RankedMin[] instead, which is more amenable to symbolic processing, at the cost of being quite a bit slower to evaluate.

sums[x_?VectorQ] := Total[Subsets[x, {2}], {2}]

score[x_?VectorQ] := With[{s = sums[x]}, 
      CentralMoment[Differences[Table[RankedMin[s, k], {k, Length[s]}]],
                    2]/(Max[x] - Min[x])^2]

f[x_, y_] := score[{0, x, y, 1}]

Table[NMinimize[{f[x, y], 0 < x < 1 && 0 < y < 1}, {{x, 0, 1}, {y, 0, 1}}, 
                Method -> {"NelderMead", "RandomSeed" -> i}], {i, 5}] // Chop
   {{0, {x -> 0.5, y -> 0.75}}, {0, {x -> 0.5, y -> 0.25}}, {0, {x -> 0.5, y -> 0.25}},
    {0, {x -> 0.5, y -> 0.75}}, {0, {x -> 0.25, y -> 0.5}}}
$\endgroup$
1
  • $\begingroup$ @kjo, thanks for the acceptance, but I think wxffles's suggestion might be less expensive if your actual lists are long. $\endgroup$ May 11, 2016 at 1:10
3
$\begingroup$

Note that your expression for score[{0, x, y, 1}] seems to simplify considerably when $0<x,y<1$:

FullSimplify[score[{0, x, y, 1}], Thread[0 < {x, y} < 1]]
(* Out: 2/25 (7 x^2 - 2 x (2 + 3 y) + 2 (6 + (-2 + y) y)) *)

Minimize[Flatten@{%, Thread[0 < {x, y} < 1]}, {x, y}]
(* Out: {18/35, {x -> 5/7, y -> 1}} *)

The Minimize expression above also returns a warning:

Minimize::wksol: Warning: there is no minimum in the region in which the objective function is defined and the constraints are satisfied; returning a result on the boundary.

The symbolic result corresponds to the numerical one reported by FindMinimum and NMinimize.

$\endgroup$
2
  • $\begingroup$ I'm surprised that Mathematica can simplify the expression. There's a Sort step in the expression for score, and I did not think that was amenable to analytic reduction. Maybe therein lies the bug. $\endgroup$
    – kjo
    May 11, 2016 at 0:15
  • $\begingroup$ Afterthought: I think I did not get the point of your answer at first. Now (especially after reading wxffles's comment) I see that it identifies the source of the problem. $\endgroup$
    – kjo
    May 11, 2016 at 1:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.