FindMaxValue specifics

Question

I'm using FindMaxValue to study the distribution of maxima of Abs[RiemannSiegelZ[t]] between consecutive values of ZetaZero,

Some background: RiemannSiegelZ[t] is Zeta[1/2+I* t] 'with the phase taken out'; i.e. real for real $t$ and with the same absolute value as Zeta. Outside of Mathematica this is known as Hardy's function $Z(t)$. On the Riemann Hypothesis, $Z(t)$ has no positive local minima or negative local maxima 'for $t$ sufficiently large' (see Edwards' book "Riemann's Zeta Function".) In fact $t>5$ is sufficient. Thus on the Riemann Hypothesis $Z(t)$ is positive (resp. negative) and concave down (resp. up) with a unique local maximum (resp. minimum) for $t$ between ZetaZero[2k-1] and ZetaZero[2k] (resp ZetaZero[2k] and ZetaZero[2k+1]).

To determine the maxima, I'm less concerned with speed or precision than I am with robustness. The documentation says that

FindMaxValue[f,{x, x0, x1}] searches for a local maximum in f using x0 and x1 as the first two values of x, avoiding the use of derivatives.

Is this routine guaranteed to return a value between x0 and x1 if there is in fact a unique local maximum there (i.e., guaranteed to not jump outside the interval)? Should I be using a Method other than Automatic, and if so, which one?

The documentation also says:

FindMaxValue[f,{x, x0, xmin, xmax}] searches for a local maximum, stopping the search if x ever gets outside the range xmin to xmax.

It does not say what the parameter x0 determines but from earlier documentation one would suppose that it is the starting point for the search. What does it mean if the search 'stops'? Does Mathematica return an error message?

What is the difference between

FindMaxValue[f,{x, x0, xmin, xmax}]

and

FindMaxValue[{f, xmin < x < xmax}, {x, x0}]

Finally, I'm having difficulties with (I think) the fact that FindMaxValue has the attribute HoldAll, for example

In[3]:= t1 = N[Im[ZetaZero[1]]]; t2 = 
 N[Im[ZetaZero[2]]]; tmid = (t1 + t2)/2; FindMaxValue[
 RiemannSiegelZ[t], {t, tmid, t1, t2}]

FindMaxValue::nrgnum: The gradient is not a vector of real numbers at {t} = {17.8742}. >>

Out[3]= 2.31595

Answer 1

Is this routine guaranteed to return a value between x0 and x1 if there is in fact a unique local maximum there (i.e., guaranteed to not jump outside the interval)?

Nope. Consider for instance a simple example:

FindArgMax[Sin[x], {x, π/2 - π/6, π/2 - π/12}]
   {1.570796326793869}

where we see that the known maximum at $x=\pi/2$ was found, which is of course outside the bounds given. On the other hand, look at the following:

FindArgMax[Sin[x], {x, 3 π/8, π/2 - π/6, π/2 - π/12}] // Quiet
   {1.3089969389957472}

(I had to suppress the warning message that the optimum was found at the interval's edge for this, though.) On the other hand, I would admit that I don't know if the local optimization functions retain the possibility of returning a point outside a given interval if the function being optimized looks "parabolic" within the given interval.

---

Since the derivative of RiemannSiegelZ looks a bit nasty (try executing RiemannSiegelZ'!), one can either elect to use finite differences for approximating the gradient, or using a gradient-free method altogether.

To illustrate the first approach,

t1 = Im[N[ZetaZero[1]]]; t2 = Im[N[ZetaZero[2]]];

FindMaximum[RiemannSiegelZ[t], {t, (t1 + t2)/2, t1, t2}, 
            Gradient :> {"FiniteDifference", "DifferenceOrder" -> 2}]
   {2.340551029908819, {t -> 17.882582076641626}}

For the second approach, the only genuinely derivative-free method that seems to be available is the "principal axis" method of Richard Brent:

FindMaximum[RiemannSiegelZ[t], {t, (t1 + t2)/2, t1, t2},
            Method -> "PrincipalAxis"]
   {2.3405510299088164, {t -> 17.88258207729408}}

(Mathematica apparently has an implementation of the genuinely univariate optimization method by Brent, but alas, it is only available for line searches. A pity.)

I have used FindMaximum[] here, but the same should apply to its cognate functions FindArgMax[] and FindMaxValue[].