C[omp]ute

From: andrew cooke <andrew@...>

Date: Fri, 14 Dec 2018 18:04:10 -0300

Maybe this is already known - I guess it must be - but I just dreamt
it up myself.  I need to find the max (could be min, obvs) of a noise
function in 1D.

Since it's noisy I can't assume it's smooth, hope it has a global
minimum, and bisect.  I need something more robust.

At the same time, I don't want to be doing repeated evaluations of the
function where they're not needed, so any iterative deepening of a
search should re-use old values when possible.

So here's the solution:

  * Evaluate the function at 5 equally spaced points across the range
    (including the two extremes).

  * Throw away two end points.

  * Expand the remaining 3 points back to 5 by inserting new points at
    the mid-points of the existing points.

  * Repeat until you're happy.

The trick is to know which end points to discard.  It could be one
either end, or two from one end.

What I do is calculate the average of the x values at the points,
weighted by the function values there.  Then I compare this to the
unweighted average.  If it's higher, then the max is towards the high
end, so discard a low point.  Or vice-versa.  And repeat for a second
point.

Here's an example.  The lines are (x, f(x)) pairs:

[(0.0, 0), (0.25, 5), (0.5, 10), (0.75, 1), (1.0, 1)]
[(0.25, 5), (0.375, 12), (0.5, 10), (0.625, 2), (0.75, 1)]
[(0.25, 5), (0.3125, 11), (0.375, 12), (0.4375, 10), (0.5, 10)]
[(0.3125, 11), (0.34375, 12), (0.375, 12), (0.40625, 10), (0.4375, 10)]
[(0.3125, 11), (0.328125, 11), (0.34375, 12), (0.359375, 13), (0.375, 12)]

On the first step, both ends were discarded.  On the second, two from
the "high" side, etc.  The maximum is 13 at 0.359375, roughly.

And here's the code:

def expand_max(lo, hi, n, f):
    data = [(x, f(x)) for x in (lo + i * (hi - lo) / 4 for i in range(5))]
    for _ in range(n):
        print(data)
        while len(data) > 3:
            w = sum(x*fx for (i, (x, fx)) in enumerate(data)) / sum(fx for (x, fx) in data)
            m = sum(x for (x, fx) in data) / len(data)
            if w > m:
                del data[0]
            else:
                del data[-1]
        x = (data[0][0] + data[1][0]) / 2
        data.insert(1, (x, f(x)))
        x = (data[-2][0] + data[-1][0]) / 2
        data.insert(-1, (x, f(x)))
    x_max, fx_max = None, None
    for (x, fx) in data:
        if x_max is None or fx > fx_max:
            x_max, fx_max = x, fx
    return x_max, fx_max

Andrew