I have an expensive function f which is unimodal and I want to find its minimum. However f is only defined at integer values. I read that golden section search is the right thing to do. My implementation which ignores the integer restriction is:
def golden_section_search(f, a, b, tolerance=1e-6):
"""
Perform the Golden Section Search for finding the minimum of a unimodal function f.
Parameters:
f (function): The unimodal function to minimize.
a (float): Left boundary of the interval.
b (float): Right boundary of the interval.
tolerance (float, optional): Tolerance for stopping criterion. Default is 1e-6.
Returns:
float: Argument that minimizes the function f within the interval [a, b].
"""
phi = (math.sqrt(5) - 1) / 2 # Golden ratio constant
# Initial points
x1 = a + (1 - phi) * (b - a)
x2 = a + phi * (b - a)
# Initial function evaluations
f_x1 = f(x1)
f_x2 = f(x2)
while (b - a) > tolerance:
print(a, b)
if f_x1 < f_x2:
b = x2 # Discard the interval [x2, b]
x2 = x1 # Move x2 to the left
f_x2 = f_x1 # Reuse the already computed value
# Recompute x1 and its function value
x1 = a + (1 - phi) * (b - a)
f_x1 = f(x1)
else:
a = x1 # Discard the interval [a, x1]
x1 = x2 # Move x1 to the right
f_x1 = f_x2 # Reuse the already computed value
# Recompute x2 and its function value
x2 = a + phi * (b - a)
f_x2 = f(x2)
# The minimum is at the midpoint of the final interval
return (a + b) / 2
As a test function let's use this (this is of course defined everywhere but we pretend it isn't):
def f(x):
return (x - 3)**2 + 5
To make sure the search only evaluates at integers I tried just doing round(x1) and round(x2) every time x1 and x2 are changed but that made a version that either gives the wrong answer or iterates forever. It gives the wrong answer if I set tolerance = 1 and iterates forever if I set tolerance = 0.9. I am using a = 0 and b = 5 in the function call.
The true minimum is at x = 3.
This is the code which only evaluates f at integer values but does not work:
def golden_section_search_integer(f, a, b, tolerance=1e-6):
"""
Perform the Golden Section Search for finding the minimum of a unimodal function f,
but only evaluate f at integer values.
Parameters:
f (function): The unimodal function to minimize.
a (int): Left boundary of the interval (integer).
b (int): Right boundary of the interval (integer).
tolerance (float, optional): Tolerance for stopping criterion. Default is 1e-6.
Returns:
int: Integer value that minimizes the function f within the interval [a, b].
"""
phi = (math.sqrt(5) - 1) / 2 # Golden ratio constant
# Initial points, rounded to integers
x1 = round(a + (1 - phi) * (b - a))
x2 = round(a + phi * (b - a))
# Initial function evaluations
f_x1 = f(x1)
f_x2 = f(x2)
while (b - a) > tolerance:
print(a, b, x1, x2)
if f_x1 < f_x2:
b = x2 # Discard the interval [x2, b]
x2 = x1 # Move x2 to the left
f_x2 = f_x1 # Reuse the already computed value
# Recompute x1, ensuring it's an integer, and its function value
x1 = round(a + (1 - phi) * (b - a))
if x1 != x2: # Only evaluate if x1 is different from x2
f_x1 = f(x1)
else:
a = x1 # Discard the interval [a, x1]
x1 = x2 # Move x1 to the right
f_x1 = f_x2 # Reuse the already computed value
# Recompute x2, ensuring it's an integer, and its function value
x2 = round(a + phi * (b - a))
if x2 != x1: # Only evaluate if x2 is different from x1
f_x2 = f(x2)
# The minimum is at the midpoint of the final interval, rounded to the nearest integer
return round((a + b) / 2)
x1 and x2 can and should be float values otherwise your convergence values will be quantized and do weird things around 1.0. Only the input to f() needs to be integral which can be done without modifying x{1,2}:
def golden_section_search_integer(f, a, b, tolerance=1e-6):
# ...
# Initial function evaluations
f_x1 = f(round(x1)) # <-- round
f_x2 = f(round(x2)) # <-- round
while (b - a) > tolerance:
print(a, b)
if f_x1 < f_x2:
# ...
f_x1 = f(round(x1)) # <-- round
else:
# ...
f_x2 = f(round(x2)) # <-- round
# ...
return (a + b) / 2 # <-- round here if necessary
>>> golden_section_search_integer(f, 0, 5, 0.9)
0 5
1.9098300562505255 5
1.9098300562505255 3.819660112501052
2.639320225002103 3.819660112501052
3.454915028125263