SIAM 2002 challenge #5

5. Let f(z) = 1/G(z), where G(z) is the gamma function, and let p(z) be the cubic polynomial that best approximates f(z) on the unit disk in the supremum norm ||·||oo What is || f - p||oo ?

Translation (thanks to BM): find the maximum of | f(z) - p(z)| on the unit disk |z| <= 1 for that cubic p(z) for which that maximum is smallest.

This solution was supplied by Marijke van Gans in HTML form

Above: Plot of f = 1/G(z) over the unit disk. Re part cyan, Im part yellow, z abscissa grid grey. We're looking approx. in the direction of the +ve real axis (marked "1"). The plot only looks asymmetric because of perspective.

Above: Plot of  f - p (colors and view as before), for a cubic p that has the best fit in terms of the problem.

Below: Plot of | f - p| for that p.

Preliminary analysis:

So we're left with a search over 4-D space with the 4 (real) coordinates being a, b, c, d, the cubic's coefficients, and we only have to look at function values along the unit circle rather than disk.

Why this is a deviously hard problem <g>

Consider. Suppose that at one "point" p(z) in abcd space, the corresponding | f(z) - p(z)| has k local maxima (on the unit circle), m0, m1, m2 ... say, attained at z0 and z0*, at z1 and z1*, at z2 and z2*, ... say. Also suppose that over some region of abcd space the corresponding function continues to have the "same" k local maxima (i.e. the values of the mi may vary, their location zi may vary, but identity is maintained, no local maxima arise, disappear, or merge in this region). Now what do we actually minimise?

The global maximum m of | f - p| for any one function p is max(m0, m1, m2 ...).

Let us plot the height of m0 at each point of our 4D region in abcd space. Again, let us plot m1. Again, plot m2 ... Each of those plots forms an analytically smooth 4-D sheet on which the usual minimisation techniques would work well. But what would a sheet for m look like?

Visualise the problem with just a 2-D control space rather than our 4-D control space abcd. Now we can use the third dimension to plot the height of m0, and (superimposed on the same 3-D plot) m1, m2 and so on. Each mi plots as a landscape, sloping hillsides. m however is everywhere the highest of these topographies. So the sheet for m can have bad discontinuities. In one region we may find that m0 and m1 form the twin banks of a "gully" that has a V crosssection, with m being m0 on the one "bank" and m1 on the other "bank".

There are numeric techniques that work with narrow steep sided valleys, but they all assume that gradient is meaningful, and doesn't become a singularity just where you need it most. We can (and indeed do) find a gully with steep banks whereas the gully itself has only a very shallow descent downstream. In this situation, almost all directions are up. There's only a very narrow wedge of directions in which you can go down. Once in the gully, a cursory check in all ± coordinate axis directions gives the misleading impression we're at the minimum because "all" directions are up.

Conjugate gradient methods, direction set methods, simplex methods, they all fail if, as here, the gully is truly relentlessly undifferentiably V-shaped rather than just a rather sharp U at small enough scale. Every trick in the book gets stuck in the mud of this gully.

Time for new tricks

But first, consider that in 3-D or 4-D there can be tunnely gullies with 3 banks, and in 4-D even 4-bank ones. I didn't come across a 4-bank one but did wind my way along a 3-bank gully.

The key of the method i used is to look at the gradients of the individual mi sheets, and concoct a step vector that has a nonnegative dot product with all the k individual gradients (then, of course, step in the opposite direction).

How to find local and global maxima for any one cubic (a, b, c, d)? Scan, say, 512 points along unit circle and then, starting from local maxima of that coarse set, home in further by bisection until the difference is not worth bothering. Keep up to 3 maxima (coupled with visual inspection of plot with color coding to highlight maxima, which showed 3 was enough).

How to find gradients? Let current step size be s, and current location (a, b, c, d). Then look at the i-th local maximum of (a ± s, b, c, d), subtract the one for a - s from that for a + s, that gives you the a-component. Repeat for b ± s, c ± s, d ± s. Four components give you the gradient of the mi sheet. Repeat for all k different values of i.

How to find out how many banks you're on? Keep the values rather than just the differences. If the mi sheet is so far above the mj sheet that you can not reach the intersection in one step then mj does not play a rôle, this step. Look only at the mi that are near enough to the highest local maximum.

How to behave when on one bank? Just slide down the gradient.

How to behave when on two banks? Let u0 and u1 be the relevant gradients. Now normalise u0 and u1 to unit length, then add. This has nonnegative dot product with both.

How to behave when on three banks? Let u0 and u1 and u2 be the relevant gradients. We must (thanks to BM) find vectors v0 perpendicular to both u1 and u2 (and in positive u0 sense); v1 perp. to both u2 and u0; and v2 perp. to both u0 and u1. Then xv0 + yv1 + zv2 with positive x, y, z works. For a 3-D control space, cross products of the u's give the right v's. For 4-D, i used a similar but longer expression (6 terms per component).

Above: at an earlier stage i believed d should be zero. This screenshot is from a survey in abc space (9×9×9 cubic grid, tasting the global max of | f - p| for p ranging up to ±4 steps away from current p in each direction), and shows well how thin the wedge of directions is in which the global maximum is lower (green and yellow), once you're in the gully.

Above: from a survey in full abcd space (5×5×5×5 grid, tasting the global max of | f - p| for p ranging up to ±2 steps away from current p in each direction).

Below: same view, but in stead of goodness of global max we display which of 3 banks we're on (i.e. which of the 3 local maxima forms the global one).

My search for p started at the Taylor series for f, truncated to cubic, and the program got tweaked numerous times as the terrain got harder, each time continuing from the best coefficients found by the previous version. After ironing out many bugs the program finally did the steps right, sort of. I never sorted out the problem of disappearing or merging local maxima, i just kept my finger on the key to do a step at a time, having assigned the backspace to an undo function, saving often, and restarting from last saved whenever the program got lost. Decreasing step size whenever the gully became too shallow and steps would overshoot. The program did not try to adjust the relative weight of the v's but just alternated a few 3-bank steps with a 2- or 1-bank step automatically whenever it had climbed too far up some of the banks.

Finite accuracy of C doubles means that most likely not all the digits of the result are those of the true minimum maximum [as it happens, that was too pessimistic: the method did get the required ten digits right :) -- MvG]. The program has oscillated between incarnations in UBasic (for virtually unlimited precision) and ones in C for Windows (for the nice visuals and effortless incorporation of lots of menu options to tweak what to show and how to step). It is in the latter form that it finally became relatively bugfree enough to reach

      0.2143352346

somewhere in the lower reaches of a 3-bank gully for coinciding local maxima at z = -1 and two pairs of conjugate z:

plot of |f-p| for p = -0.603343177993 z^3 + 0.625212051899 z^2 + 1.019761988440 z + 0.005541993177

Above: | f - p| again for possibly one of its lowest maxima (the program just reloaded the final coefficients from where i called it a day). This display option has color-coded hatching to highlight the five local maxima, as height differences this shallow wouldn't show up clearly on just the 3-D plot.

PS

Cornelius Lanczos' celebrated approximation method for complex Gamma, as sum of a finite number of terms, actually consists of an infinite family of such algorithms, for different integer parameter g and different number of terms. The tricky bit is finding good coefficients for the terms. One of his solutions, with g = 5 and seven terms up to N = 6, is well known as a canned method (see e.g. section 6.1 of Numerical Recipes), it has accuracy good to 2 parts in 10-10 across half the complex plane (the other half being accessible by Gamma's reflexion formula).

It is worth noting that Paul Godfrey offers not only a ready-made solution with g = 9 and eleven terms up to N = 10, good to 1 part in 10-13, but also a method to roll your own coefficients for any precision.

As it happens, i did code up Godfrey's N = 10 Gamma in my UBasic programs but the C version that ended up finding the result to ten places still used Lanczos' N = 6 Gamma. Both programs used a version of the reflexion formula that mirrors about Re z = +1, and i checked that within the relevant region the two Gammas were both accurate to enough decimal places.

Brian Medley has been mentioned above several times already for his valuable help with this problem. He also found Godfrey's web page as soon as i expressed concern about the accuracy of the Gamma i was using.

--marijke

<back