FindMaximum ignores MaxIterations option

Question

I'm trying to find a maximum of interpolated function inside a particular mesh cell. However, for some cells the process is very slow and it returns a warning:

FindMinimum::eit: The algorithm does not converge to the tolerance of 1.`*^-6 in 500 iterations

The code is:

FF = Interpolation[q, InterpolationOrder -> 2];
dp = 
  DensityPlot[FF[x, y], {x, y} ∈ Disk[{0, 0}, 2], 
    PlotRange -> All, PlotPoints -> 40];
ms = DelaunayMesh[kekL, MeshCellStyle -> {2 -> Opacity[0]}];
pg = MeshPrimitives[ms, 2][[104]];
fm = 
  FindMaximum[{FF[x, y], {x, y} ∈ pg}, {{x, pg[[1, 1, 1]]}, {y, pg[[1, 1, 2]]}}, 
    MaxIterations -> 30]
Show[
  {dp, ms,
   Graphics[{White, Point[{x, y} /. fm[[2]]], Point[RegionCentroid[pg]]}]}, 
  ImageSize -> 600]

The definitions of q and kekL are on pastebin.

I don't want a particularly accurate solution (even 10^-2 precision is fine)m but I'd like it to be fast. I tried using option MaxIterations but no matter what the parameter was it still reported "500 iterations" warning. It seems like it's ignoring eveything I put like AccuracyGoal and PrecisionGoal as well. The version is 11.2

Answer 1

I suspect that the use of a region is calling an optimizer and options are not being passed (*see below). It should probably be reported to WRI as a potential bug to see what they say.

Here's a way to enforce the constraint, by constructing one's own penalty function. We also use a second starting point so that FindMaximum will avoid using the derivative of FF, which is not smooth.

dist[x_, y_] = RegionDistance[pg, {x, y}];
penalty[x_?NumericQ, y_?NumericQ] := Unitize[#] + # &@dist[x, y];
fm = FindMaximum[
  FF[x, y] - 
   penalty[x, y], {{x, y}, RegionCentroid@pg, pg[[1, 1]]}\[Transpose],
   MaxIterations -> 30]
(*  {0.200105, {x -> 0.789665, y -> -1.60565}}  *)

Show[
 Plot3D[{FF[x, y], FF[x, y] - penalty[x, y]}, {x, y} ∈ pg],
 Graphics3D[{Red, Sphere[{x, y, FF[x, y]} /. Last@fm, 0.01]}],
 PlotRange -> All
 ]

Update: Another workaround

The discussion in the addendum below suggests using RegionMember instead of Element as a workaround:

FindMaximum[{FF[x, y], 
  RegionMember[pg, {x, y}]}, {{x, pg[[1, 1, 1]]}, {y, pg[[1, 1, 2]]}},
  MaxIterations -> 2]

FindMaximum::cvmit: Failed to converge to the requested accuracy or precision within 2 iterations.

(* {0.195039, {x -> 0.742393, y -> -1.64621}} *)

With the option settings as follows, we get a result without error/warning messages:

FindMaximum[{FF[x, y], RegionMember[pg, {x, y}]},
 {{x, pg[[1, 1, 1]]}, {y, pg[[1, 1, 2]]}},
 PrecisionGoal -> 2, AccuracyGoal -> 2, MaxIterations -> 7]
(*  {0.200033, {x -> 0.788627, y -> -1.60954}}  *)

---

Addendum.

The following shows that FindMinimum is called with the default setting for the options. For some reason, after FindMinimum fails multiple times, FindMaximum seems to be run. I'm pretty sure it's an internal function being run as a proxy, but it's curious. A Trace[] show it is running FindMinimum on -FF[x,y] for two of the three line segments making up the triangle pg and then on the polygon; it seems it skips one of the line since it didn't show up. It would be nice to eliminate what seem to be extraneous FindMinimum runs, but they might be helping to narrow down a solution, I suppose. Someone with time to Trace[] the computation might be able to shed more light.

With[{opts = Options@FindMinimum},
 Internal`WithLocalSettings[
  SetOptions[FindMinimum, MaxIterations -> 1],
  ffm = FindMaximum[{FF[x, y], {x, y} ∈ pg}, {{x, 
      pg[[1, 1, 1]]}, {y, pg[[1, 1, 2]]}}, MaxIterations -> 30],
  SetOptions[FindMinimum, opts]
  ]
 ]

FindMinimum::cvmit: Failed to converge to the requested accuracy or precision within 1 iterations.
... [two more messages] ...
General::stop: Further output of FindMinimum::cvmit will be suppressed during this calculation.

FindMaximum::cvmit: Failed to converge to the requested accuracy or precision within 30 iterations.

(*  {0.200079, {x -> 0.794926, y -> -1.60892}} v *)

By reducing PrecisionGoal and AccuracyGoal, one can also reduce the number of iterations needed. Note the option values have to be set in both places:

With[{opts = Options@FindMinimum},
 Internal`WithLocalSettings[
  SetOptions[FindMinimum,
   {PrecisionGoal -> 2, AccuracyGoal -> 2, MaxIterations -> 30}],
  ffm = FindMaximum[{FF[x, y], {x, y} ∈ pg}, {{x, 
      pg[[1, 1, 1]]}, {y, pg[[1, 1, 2]]}},
    PrecisionGoal -> 2, AccuracyGoal -> 1, MaxIterations -> 30],
  SetOptions[FindMinimum, opts]
  ]
 ]

FindMinimum::eit: The algorithm does not converge to the tolerance of 0.01` in 30 iterations. The best estimated solution, with feasibility residual, KKT residual, or complementary residual of {2.78834*10^-11,0.0152532,1.39416*10^-11}, is returned.

(*  {0.200107, {x -> 0.790959, y -> -1.60642}}  *)

Here is the tracing code:

With[{opts = Options@FindMinimum},
 Internal`WithLocalSettings[
  SetOptions[FindMinimum, MaxIterations -> 1],
  Trace[
   FindMaximum[{FF[x, y], {x, y} ∈ pg}, {{x, 
      pg[[1, 1, 1]]}, {y, pg[[1, 1, 2]]}}, MaxIterations -> 30],
   _FindMinimum,
   TraceInternal -> True],
  SetOptions[FindMinimum, opts]
  ]
 ]

Answer 2

With method Spline in Interpolation works fine.

 FF = Interpolation[q, InterpolationOrder -> 2, Method -> "Spline"]
 dp = DensityPlot[FF[x, y], {x, y} \[Element] Disk[{0, 0}, 2], PlotRange -> All, PlotPoints -> 40];
 ms = DelaunayMesh[kekL, MeshCellStyle -> {2 -> Opacity[0]}];
 pg = MeshPrimitives[ms, 2][[104]];
 fm = FindMaximum[{FF[x, y], {x, y} \[Element] pg}, {{x, pg[[1, 1, 1]]}, {y, pg[[1, 1, 2]]}}, MaxIterations -> 30]

 (* {0.203597, {x -> 0.786015, y -> -1.64726}} *)

 Show[{dp, ms, Graphics[{White, Point[{x, y} /. fm[[2]]], 
 Point[RegionCentroid[pg]]}]}, ImageSize -> 600]