Numerical anomalies in a minimax algorithm

[deberly (David Eberly)]

I am trying to compute error bounds for polynomial estimates to \sin(t\theta)/\sin(\theta) for t \in [0,1]. The polynomials are of the form p(x,t), where x = \cos(\theta). The polynomials have a constant u \in [0,1] that I want to choose to minimize the maximum error. The mathematical derivation is irrelevant, so I have skipped those details. I wrote code (Mathematica 11.3) to do this and plotted the minimax result as a function of u. (I omitted the NMinimize call for h[u] in this code sample.)

a[i_?NumericQ, t_?NumericQ] := If[i >= 1, a[i - 1, t]*(t^2 - i^2)/(i*(2*i + 1)), t]
p[n_?NumericQ, y_?NumericQ, t_?NumericQ] := (sum = a[n, t]; For[i = n - 1, i >= 0, i--, sum = a[i, t] + sum*y]; sum)
f[n_?NumericQ, x_?NumericQ, t_?NumericQ] := Sin[t*ArcCos[x]]/Sin[ArcCos[x]] - p[n, x - 1, t]
g[n_?NumericQ, u_?NumericQ, x_?NumericQ, t_?NumericQ] := Abs[f[n, x, t] - u*a[n, t]*(x - 1)^n]
h[n_?NumericQ, u_?NumericQ] := (result =  NMaximize[{g[n, u, x, t], 0 <= x <= 1 && 0 <= t <= 1}, {x, t}]; result[[1]])
Plot[h[8, u], {u, 0.7, 0.9}]

The output of Plot has some numerical anomalies.

When I program this in C++ using double precision, the function h(u) is smooth. Evaluating h[8,0.75], Mathematica produces 0.000058529. Evaluating h[8,0.751], Mathematica produces 9.13505e-06. I did not expect the sawtooth-like behavior of the graph. The valleys do not show up in my C++ computations, which shows effectively a V-shaped graph with vertex near (0.85352, 1.91558e-05). I tried to change the working precision, but the sawtooth behavior persisted.

I switched the method to “Simulated Annealing.” The output of the Plot function also has some anomalies.

The outputs at the two aforementioned locations are h[8,0.75] = 0.0000583938 and h[8,0.751] = 0.0000580131, but now the anomalies are in a different region of the graph.

Finally, I tried using “Differential Evolution” as the method. The output looks like what I expected.

I know how to debug numerical issues in C++ code using a debugger, but I am a novice at Mathematica and wish to know whether there is some standard approach or set of tools that allows me to diagnose such issues. Also, is there some general advice on choosing the method for minimizing or maximizing? Or is this simply something one has to use trial-and-error to determine? Thank you.

[trounev (Alexander Trounev)]

You do not need to use ?NumericQ and Plot in this code. Note that my code works 25 times faster and gives the same result

a[i_, t_] := If[i >= 1, a[i - 1, t]*(t^2 - i^2)/(i*(2*i + 1)), t]
p[n_, y_, t_] := (sum = a[n, t]; 
  For[i = n - 1, i >= 0, i--, sum = a[i, t] + sum*y]; sum)
f[n_, x_, t_] := Sin[t*ArcCos[x]]/Sin[ArcCos[x]] - p[n, x - 1, t]
g[n_, u_, x_, t_] := Abs[f[n, x, t] - u*a[n, t]*(x - 1)^n]
h[n_, u_] := (result = 
   NMaximize[{g[n, u, x, t], 0 <= x <= 1 && 0 <= t <= 1}, {x, t}, 
    Method -> "DifferentialEvolution"]; result[[1]])
ListLinePlot[
  Table[{u, h[8, u]}, {u, 0.7, 0.9, .005}]] // AbsoluteTiming

Your code on my computer

[deberly (David Eberly)]

I mentioned in my post that I omitted the final line of the minimax operation. After the definition of h[n_,u_] my actual code has

hmin = NMinimize[{h[8, u], 0 <= u <= 1}, u]

Without the NumericQ, this generates several lines of errors “… NMaximize: The function value 0 is not a number at (t,x) = {0.918621,0.831757}” followed by “…General: Further output of NMaximize::nnum will be suppressed during this calculation.” and the NMinimize fails (but the ListLinePlot succeeds). My goal is to obtain the minimum and where it occurs. Is there some way to avoid these errors? With NumericQ and the NMinimize line included, the program ran for over 2 hours (on a high-end desktop computer, but Mathematica appeared to be using 20% CPU [about 2 logical processors]) . Also, your code shows the same numerical anomalies when NMaximize uses its default method.

Thank you.

[trounev (Alexander Trounev)]

Can you publish the final result, which is obtained after 2 hours?

[deberly (David Eberly)]

With the NMinimize function included, the reported value is u=0.852.

That said, I tracked down the problem. To compare with my C++ results, where I use a dense sampling of (x,t) for each u in a set of 1024 and select the maximum, I used

ListLinePlot[Table[{u, h[8, u]}, {u, 0.8519, 0.8521, 0.0002/1024}]]

The output image is

which makes it appear that there is a discontinuity at u = 0.852. The call

NMaximize[{g[8, 0.8519, x, t], 0 <= x <= 1 && 0 <= t <= 1}, {x, t}, Method -> "DifferentialEvolution"]

produces {0.0000195003, {x → 0., t → 0.476764}} and the call

NMaximize[{g[8, 0.8521, x, t], 0 <= x <= 1 && 0 <= t <= 1}, {x, t}, Method -> "DifferentialEvolution"]

produces {0.0000189766, {x → 0.129624, t → 0.521673}}. The result for u = 0.8521 is what I believe to be incorrect (and in fact incorrect for any u > 0.852 in the line-list plot). I called

Plot3D[g[8, 0.8521, x, t], {x, 0, 1}, {t, 0, 1}, PlotRange -> All, PlotPoints -> 50]

which generated the output

so the graph is bimodal. For u=0.8251, NMaximize reports the local maximum at (x,t)=(0.129624,0.521673). The local maximum on the boundary x=0 occurs at t=0.476541 and the value is 0.000019424, which is larger than the local maximum in the interior. I obtained the boundary computation from

NMaximize[{g[8, 0.8521, 0, t], 0 <= t <= 1}, t, Method -> "DifferentialEvolution"]

[trounev (Alexander Trounev)]

Thank you. I got a similar result {{0.852001, 0.0000189641}} after a few minutes without using NMinimize. I created Table and found minimum using MinimalBy

[deberly (David Eberly)]

After some detailed experiments with C++ code, the bimodal behavior and switching between the two local maxima (in terms of reporting the maximum) has the flavor of Chebyshev equioscilation that is obtained via minimax algorithms. In fact I am trying to do some error balancing, and the u-parameter is what I vary to do the balancing. When you remove the Abs[] from g[…] and call Plot3D as I did, you get a local minimum (on boundary) and a local maximum (in the interior) that have same magnitude but opposite sign (within numerical rounding errors in the code). The difference between the C++ results and the Mathematica results is that the C++ code makes it clear that when one local maxima becomes smaller than the other as u varies past the switching point, the magnitudes are nearly the same.

Thanks for your input. Perhaps someone will eventually help me with those Mathematica errors when I don’t use NumericQ.

[danl (Daniel Lichtblau)]

the maximization step seems to be key. If x is restricted more tightly, say 0<=x<=.25, then NMaximize appears less likely to get trapped on the wrong side of the nearly equal divide, and the subsequent NMinimize does better. Below is a recoding, using this tighter restriction.

a[0, t_] = t;
a[n_Integer /; n >= 1, t_] := a[n - 1, t]*(t^2 - n^2)/(n*(2*n + 1))
p[n_Integer /; n >= 1, x_, t_] := 
 Table[a[j, t], {j, 0, n}].x^Range[0, n]
f[n_Integer /; n >= 1, x_, t_] := 
 Sin[t*ArcCos[x]]/Sin[ArcCos[x]] - p[n, x - 1, t]
g[n_Integer /; n >= 1, u_, x_, t_] := 
 Abs[f[n, x, t] - u*a[n, t]*(x - 1)^n]

obj8[u_] = g[8, u, x, t]

(* Out[129]= Abs[-t - 1/3 t (-1 + t^2) (-1 + x) - 
  1/30 t (-4 + t^2) (-1 + t^2) (-1 + x)^2 - 
  1/630 t (-9 + t^2) (-4 + t^2) (-1 + t^2) (-1 + x)^3 - (
  t (-16 + t^2) (-9 + t^2) (-4 + t^2) (-1 + t^2) (-1 + x)^4)/22680 - (
  t (-25 + t^2) (-16 + t^2) (-9 + t^2) (-4 + t^2) (-1 + t^2) (-1 + 
     x)^5)/1247400 - (
  t (-36 + t^2) (-25 + t^2) (-16 + t^2) (-9 + t^2) (-4 + t^2) (-1 + 
     t^2) (-1 + x)^6)/97297200 - (
  t (-49 + t^2) (-36 + t^2) (-25 + t^2) (-16 + t^2) (-9 + t^2) (-4 + 
     t^2) (-1 + t^2) (-1 + x)^7)/10216206000 - (
  t (-64 + t^2) (-49 + t^2) (-36 + t^2) (-25 + t^2) (-16 + t^2) (-9 + 
     t^2) (-4 + t^2) (-1 + t^2) (-1 + x)^8)/1389404016000 - (1/
  1389404016000)
  t (-64 + t^2) (-49 + t^2) (-36 + t^2) (-25 + t^2) (-16 + t^2) (-9 + 
     t^2) (-4 + t^2) (-1 + t^2) u (-1 + x)^8 + Sin[t ArcCos[x]]/Sqrt[
  1 - x^2]] *)

h[u_?NumericQ] := 
 NMaxValue[{obj8[u], 0 <= x <= .3, 0 <= t <= 1}, {x, t}, 
  Method -> "DifferentialEvolution"]

This next minmax might or might not be correct, but it is clearly closer than before. When I check the max of g[8, 0.8521, x, t] it no longer gets the (slightly) incorrect value, so that at least is a plus.

AbsoluteTiming[hmin = NMinimize[{h[u], 0 <= u <= 1}, u]]

(* Out[184]= {289.592236, {0.0000190048812285, {u -> 0.852323872136}}} *)