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