Golden Rule#

For simple situations a straightforward approach may be all that is required. Using different ratios we may be able to improve on the interval halving method. The golden ratio creates a very robust search method. This ratio can be thought as dividing a line into unequal parts such that the ratio of the whole to the larger part is the same as the ratio of the larger to the smaller part.

The golden search method reduces the interval by first finding out in which section contains a minimum / maximum.

\[\frac {a + b}{a} = \frac a b = \phi\]

This ratio translates into

\[\frac {\sqrt 5 +1}{2} = \phi\]

To use this rule first bracket the point of interest (maximum/minimum), by points along the x-axis at a and b such that a < b, then create two intermediate points (x1, x2) such that

\[\begin{split}x1 &= a + d \\ x2 &= b -d \\ d &= \phi \cdot (b - a)\end{split}\]

Initialise by bracketing the point of interest, insert two intermediate positions x1, x2

1. Evaluate the function at the two intermediate points f(x1), f(x2)

1.1. If f(x1) is less than f(x2) then x1 is closer to the minimum than x2.

1.2. The upper limit is discarded, which is beyond x2 (in this case b).

1.3. The lower limit a, x1, x2 become the new a, x2, b and a new x1 is calculated, x2 is the old x1.

2. else

2.1. The lower limit is discarded

2.2. x1, x2, b become the new a, x1, and b, a new x2 is calculated

2.3. Test whether the absolute difference between the new upper and lower values reach the required limit or not.

Hint

When searching for a maximum the same algorithm is used, only when comparing values f(x1) is more than f(x2) then x1 is closer to the maximum than x2

Apart from initialisation, when 2 additional points and their evaluations are required, iteration only requires one new point position and its evaluation. Three positions and their evaluation are re-used in the new iteration.

../_images/max_gold.png — Finding the mminimum, initialisation in magenta, first iteration in blue#

A plot with the initial limits and intermediate points divided according to the golden ratio

After the lower and upper limits are chosen the first two calculated points(x1 and x2) are evaluated. Since f(x1) is less than f(x2), the magenta dashed lines, at the first iteration x2 becomes the new upper limit xun, the old lower limit xl becomes the new lower limit xln, the old point x1 became x2n and x1n has to be recalculated. At the next iteration as f(x2n) is less than f(x1n), the blue dashed lines, the lower limit will be discarded and a new x2 will have to be calculated.

Look carefully at the two sets of arrows, each set of two outer arrows are equal in length, the middle overlap being about a third of either the other two arrows.

In the example below a minimum is being found, more evaluations than strictly necessary are being made in order to show the progress of the script.

from math import pi, sqrt
from prettytable import PrettyTable

x = PrettyTable(['step', 'a', 'x1', 'x2', 'b', 'f(a)', 'f(x1)', 'f(x2)', 'f(b)', 'emax'])

f = lambda r: 2 * (pi*r**2 + 50/r)

def golden(f, a, b, tol):
    step = 0
    d = (sqrt(5) - 1)/2 * (b - a)
    x2 = a + d
    x1 = b - d
    F1 = f(x1)
    F2 = f(x2)
    Fa = f(a)
    Fb = f(b)
    x.add_row([step, round(a,5), round(x1,5), round(x2,5), round(b,5), \
            round(Fa,5), round(F1,5), round(F2,5), round(Fb,5), round(abs(a-b),5)])

    while abs(F1 - F2) >= tol and step < 19:

        if F1 < F2:
            b, x2 = x2, x1
            Fb, F2 = F2, F1
            d = b - x2
            x1 = a + d
            F1 = f(x1)
        else:
            a, x1 = x1, x2
            Fa, F1 = F1, F2
            d = x1 - a
            x2 = b - d
            F2 = f(x2)

        step += 1
        x.add_row([step, round(a,5), round(x1,5), round(x2,5), round(b,5), \
            round(Fa,5), round(F1,5), round(F2,5), round(Fb,5), round(abs(a-b),5)])

    return(x1, F1, step)

A = 1
B = 5

Tol = 1e-4

x1, F1, step = golden(f, A, B, Tol)

print(x)
print(f'Minimum found : at {x1:.7f} value {F1:.7f} in {step} steps')

This produced:

+------+---------+---------+---------+---------+-----------+----------+-----------+-----------+---------+
| step |    a    |    x1   |    x2   |    b    |    f(a)   |  f(x1)   |   f(x2)   |    f(b)   |   emax  |
+------+---------+---------+---------+---------+-----------+----------+-----------+-----------+---------+
|  0   |    1    | 2.52786 | 3.47214 |    5    | 106.28319 | 79.70925 | 104.54909 | 177.07963 |    4    |
|  1   |    1    | 1.94427 | 2.52786 | 3.47214 | 106.28319 | 75.18479 |  79.70925 | 104.54909 | 2.47214 |
|  2   |    1    | 1.58359 | 1.94427 | 2.52786 | 106.28319 | 78.90432 |  75.18479 |  79.70925 | 1.52786 |
|  3   | 1.58359 | 1.94427 | 2.16718 | 2.52786 |  78.90432 | 75.18479 |  75.65298 |  79.70925 | 0.94427 |
|  4   | 1.58359 |  1.8065 | 1.94427 | 2.16718 |  78.90432 | 75.86044 |  75.18479 |  75.65298 | 0.58359 |
|  5   |  1.8065 | 1.94427 | 2.02942 | 2.16718 |  75.86044 | 75.18479 |  75.15274 |  75.65298 | 0.36068 |
|  6   | 1.94427 | 2.02942 | 2.08204 | 2.16718 |  75.18479 | 75.15274 |  75.26674 |  75.65298 | 0.22291 |
|  7   | 1.94427 | 1.99689 | 2.02942 | 2.08204 |  75.18479 | 75.13251 |  75.15274 |  75.26674 | 0.13777 |
|  8   | 1.94427 | 1.97679 | 1.99689 | 2.02942 |  75.18479 | 75.13985 |  75.13251 |  75.15274 | 0.08514 |
|  9   | 1.97679 | 1.99689 | 2.00932 | 2.02942 |  75.13985 | 75.13251 |  75.1356  |  75.15274 | 0.05262 |
|  10  | 1.97679 | 1.98922 | 1.99689 | 2.00932 |  75.13985 | 75.1335  |  75.13251 |  75.1356  | 0.03252 |
|  11  | 1.98922 | 1.99689 | 2.00164 | 2.00932 |  75.1335  | 75.13251 |  75.13301 |  75.1356  |  0.0201 |
|  12  | 1.98922 | 1.99396 | 1.99689 | 2.00164 |  75.1335  | 75.13263 |  75.13251 |  75.13301 | 0.01242 |
|  13  | 1.99396 | 1.99689 | 1.99871 | 2.00164 |  75.13263 | 75.13251 |  75.1326  |  75.13301 | 0.00768 |
+------+---------+---------+---------+---------+-----------+----------+-----------+-----------+---------+
Minimum found : at 1.9968944 value 75.1325103 in 13 steps

It was comparable with the interval halving.

This example is then revised to become recursive.

import math


invphi = (math.sqrt(5) - 1) / 2  # 1 / phi
invphi2 = (3 - math.sqrt(5)) / 2  # 1 / phi^2


def gssrec(f, a, b, tol=1e-5, h=None, x1=None, x2=None, F1=None, F2=None):
    """Golden-section search, recursive.

    Given a function f with a single local minimum in
    the interval [a,b], gss returns a subset interval
    [c,d] that contains the minimum with d-c <= tol.

    Example:
    >>> f = lambda x: (x-2)**2
    >>> a = 1
    >>> b = 5
    >>> tol = 1e-5
    >>> (c,d) = gssrec(f, a, b, tol)
    >>> print (c, d)
    1.9999959837979107 2.0000050911830893
    """

    (a, b) = (min(a, b), max(a, b))
    if h is None:
        h = b - a
    if h <= tol:
        return (a, b)
    if x1 is None:
        x1 = a + invphi2 * h
    if x2 is None:
        x2 = a + invphi * h
    if F1 is None:
        F1 = f(x1)
    if F2 is None:
        F2 = f(x2)
    if F1 < F2:  # F1 > F2 to find the maximum
        return gssrec(f, a, x2, tol, h * invphi, x1=None, F1=None, x2=x1, F2=F1)
    else:
        return gssrec(f, x1, b, tol, h * invphi, x1=x2, F1=F2, x2=None, F2=None)

f = lambda x: 2 * (math.pi*x**2 + 50/x) # (x-2)**2
a = 1
b = 5
tol = 1e-5
(c,d) = gssrec(f, a, b, tol)
print ("Minimum value of f(x) is", f((c+d)/2), "at x =", (c + d)/2)

Golden Rule#

This Page