## Efficient, Simple, Directed Maximisation of Noisy Function

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