Implement golden section search

optimization/objectives/samples.h

 #ifndef OPTIMIZATION_OBJECTIVES_SAMPLES_H
 #define OPTIMIZATION_OBJECTIVES_SAMPLES_H
 
+#include "utils/scalar.h"
+
 namespace optimization {
 
 //--------------------------------------------------------------------------------------------------
-template <typename Objective>
+template <typename Vector>
 struct Sample {
-    using Input = typename Objective::Input;
-
     Sample() {}
-    Sample(Input const& p, Scalar v)
+    Sample(Vector const& p, Scalar v)
         : point(p), value(v)
     {}
 
-    Input point;
+    Vector point;
     Scalar value;
 };
 
 //--------------------------------------------------------------------------------------------------
+// TODO: Templatize on the underlying data types, not the objective function.
 template <typename Objective>
 struct GradientSample {
     using Input = typename Objective::Input;

optimization/optimizers/line_search/golden_section.h

+#ifndef OPTIMIZATION_OPTIMIZERS_LINE_SEARCH_GOLDEN_SECTION_H
+#define OPTIMIZATION_OPTIMIZERS_LINE_SEARCH_GOLDEN_SECTION_H
+
+#include "bracket.h"
+#include "objectives/samples.h"
+#include <cmath>
+
+namespace optimization {
+
+//--------------------------------------------------------------------------------------------------
+class GoldenSection {
+public:
+    // Note: This tolerance default is suitable for double precision floats.
+    GoldenSection(Scalar tolerance = 3e-8): tolerance_(tolerance), n_evaluations_(0) {}
+
+    uint32_t n_evaluations() const { return n_evaluations_; }
+
+    template <typename Objective>
+    Sample<Scalar> optimize(Objective const& objective, Bracket const& bracket);
+
+private:
+    Scalar tolerance_;
+    uint32_t n_evaluations_;
+
+    static constexpr Scalar golden_ratio_big_ = 0.618034;
+    static constexpr Scalar golden_ratio_small_ = Scalar(1) - golden_ratio_big_;
+};
+
+//..................................................................................................
+template <typename Objective>
+Sample<Scalar> GoldenSection::optimize(Objective const& objective, Bracket const& bracket) {
+    // Invariants:
+    //   x_0 < x_1 < x_2 < x_3
+    //   y_1 < f(x_0) and y_1 < f(y_3)
+    //   y_2 < f(x_0) and y_2 < f(y_3)
+    Scalar x_0, x_1, x_2, x_3;
+    Scalar y_1, y_2;
+
+    // Copy in the edges of the bracket, ensuring order is respected.
+    if (bracket.x_1() < bracket.x_3()) {
+        x_0 = bracket.x_1();
+        x_3 = bracket.x_3();
+    } else {
+        x_0 = bracket.x_3();
+        x_3 = bracket.x_1();
+    }
+
+    // Copy in the middle point of the bracket, and perform the first golden section interpolation.
+    if (bracket.x_2() - x_0 <= x_3 - bracket.x_2()) {
+        // If the middle of the bracket is closer to the left edge, interpolate into the right half
+        // of the bracket.
+        x_1 = bracket.x_2();
+        y_1 = bracket.y_2();
+        x_2 = x_1 + golden_ratio_small_ * (x_3 - x_1);
+        objective.eval(x_2, y_2);
+    } else {
+        // If the middle of the bracket is closer to the right edge, interpolate into the left half
+        // of the bracket.
+        x_2 = bracket.x_2();
+        y_2 = bracket.y_2();
+        x_1 = x_2 - golden_ratio_small_ * (x_2 - x_0);
+        objective.eval(x_1, y_1);
+    }
+    n_evaluations_ = 1;
+
+    // Keep interpolating until our bracket is sufficiently tight.
+    // See Numerical Recipes for the thought behind this particular test.
+    while (x_3 - x_0 > tolerance_ * (std::abs(x_1) + std::abs(x_2))) {
+        if (y_2 < y_1) {
+            // x_1 is our new left edge; interpolate between x_2 and x_3
+            shift(x_0, x_1, x_2, golden_ratio_big_ * x_2 + golden_ratio_small_ * x_3);
+            y_1 = y_2;
+            objective.eval(x_2, y_2);
+            ++n_evaluations_;
+        } else {
+            // x_2 is our new left edge; interpolate between x_0 and x_1
+            shift(x_3, x_2, x_1, golden_ratio_small_ * x_0 + golden_ratio_big_ * x_1);
+            y_2 = y_1;
+            objective.eval(x_1, y_1);
+            ++n_evaluations_;
+        }
+    }
+
+    // Return the smaller of the two inner points.
+    return (y_1 < y_2) ? Sample<Scalar>(x_1, y_1) : Sample<Scalar>(x_2, y_2);
+}
+
+}
+
+#endif