CMA-ES

doc/math.qbk

@@ -728,6 +728,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [include optimization/differential_evolution.qbk]
 [include optimization/jso.qbk]
 [include optimization/random_search.qbk]
+[include optimization/cma_es.qbk]
 [endmathpart] [/mathpart optimization Optimization]
 
 [mathpart poly Polynomials and Rational Functions]

doc/optimization/cma_es.qbk

@@ -0,0 +1,89 @@
+[/
+Copyright (c) 2024 Nick Thompson
+Use, modification and distribution are subject to the
+Boost Software License, Version 1.0. (See accompanying file
+LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+]
+
+[section:cma_es Evolution Strategy with Covariance Matrix Adaptation]
+
+[heading Synopsis]
+
+``
+    #include <boost/math/optimization/cma_es.hpp>
+
+    namespace boost::math::optimization {
+
+    template <typename ArgumentContainer> struct cma_es_parameters {
+        using Real = typename ArgumentContainer::value_type;
+        ArgumentContainer lower_bounds;
+        ArgumentContainer upper_bounds;
+        size_t max_generations = 1000;
+        ArgumentContainer const * initial_guess = nullptr;
+        // If zero, choose the default from the reference.
+        size_t population_size = 0;
+        Real learning_rate = 1;
+    };
+
+    template <typename ArgumentContainer, class Func, class URBG>
+    ArgumentContainer cma_es(
+        const Func cost_function,
+        random_search_parameters<ArgumentContainer> const &params,
+        URBG &gen,
+        std::invoke_result_t<Func, ArgumentContainer> target_value = std::numeric_limits<std::invoke_result_t<Func, ArgumentContainer>>::quiet_NaN(),
+        std::atomic<bool> *cancellation = nullptr,
+        std::atomic<std::invoke_result_t<Func, ArgumentContainer>> *current_minimum_cost = nullptr,
+        std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr);
+
+    } // namespaces
+``
+
+The `cma_es` optimizer searches for a global minimum of a function.
+Our implementation attempts to follow "The CMA evolution strategy: A tutorial" exactly.
+
+
+[heading Parameters]
+
+    `lower_bounds`: A container representing the lower bounds of the optimization space along each dimension. The `.size()` of the bounds should return the dimension of the problem.
+    `upper_bounds`: A container representing the upper bounds of the optimization space along each dimension. It should have the same size of `lower_bounds`, and each element should be >= the corresponding element of `lower_bounds`.
+    `initial_guess`: An optional guess for where we should start looking for solutions. This is provided for consistency with other optimization functions-it's not particularly useful for this function.
+    `max_generations`: The maximum number of generations before giving up.
+    `population_size`: The number of function calls per generation. Defaults to 4 + 4log(n), where n is the dimension.
+    `learning_rate`: Unlikely to be necessary-refer to the reference for when this should be changed.
+
+[heading The function]
+
+``
+template <typename ArgumentContainer, class Func, class URBG>
+ArgumentContainer cma_es(const Func cost_function,
+    cma_es_parameters<ArgumentContainer> const &params,
+    URBG &gen,
+    std::invoke_result_t<Func, ArgumentContainer> value_to_reach = std::numeric_limits<std::invoke_result_t<Func, ArgumentContainer>>::quiet_NaN(),
+    std::atomic<bool> *cancellation = nullptr,
+    std::atomic<std::invoke_result_t<Func, ArgumentContainer>> *current_minimum_cost = nullptr,
+    std::vector<std::pair<ArgumentContainer, std::invoke_result_t<Func, ArgumentContainer>>> *queries = nullptr)
+``
+
+Parameters:
+
+    `cost_function`: The cost function to be minimized.
+    `params`: The parameters to the algorithm as described above.
+    `gen`: A uniform random bit generator, like `std::mt19937_64`.
+    `value_to_reach`: An optional value that, if reached, stops the optimization. This is the most robust way to terminate the calculation, but in most cases the optimal value of the cost function is unknown. If it is, use it! Physical considerations can often be used to find optimal values for cost functions.
+    `cancellation`: An optional atomic boolean to allow the user to stop the computation and gracefully return the best result found up to that point. N.B.: Cancellation is not immediate; the in-progress generation finishes.
+    `current_minimum_cost`: An optional atomic variable to store the current minimum cost during optimization. This allows developers to (e.g.) plot the progress of the minimization over time and in conjunction with the `cancellation` argument allow halting the computation when the progress stagnates.
+    `queries`: An optional vector to store intermediate results during optimization. This is useful for debugging and perhaps volume rendering of the objective function after the calculation is complete.
+
+Returns:
+
+The argument vector corresponding to the minimum cost found by the optimization.
+
+[h4:examples Examples]
+
+An example exhibiting graceful cancellation and progress observability can be studied in [@../../example/cma_es_example.cpp cma_es_example.cpp].
+
+[h4:references References]
+
+Hansen, N. (2016). The CMA evolution strategy: A tutorial. arXiv preprint arXiv:1604.00772.
+
+[endsect] [/section:cma_es]

example/cma_es_example.cpp

@@ -0,0 +1,76 @@
+/*
+ * Copyright Nick Thompson, 2024
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+#if __APPLE__ || __linux__
+#include <signal.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <future>
+#include <chrono>
+#include <iostream>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/optimization/cma_es.hpp>
+
+using boost::math::optimization::cma_es_parameters;
+using boost::math::optimization::cma_es;
+using namespace std::chrono_literals;
+
+template <class Real> Real rastrigin(std::vector<Real> const &v) {
+  using std::cos;
+  using boost::math::constants::two_pi;
+  Real A = 10;
+  Real y = 10 * v.size();
+  for (auto x : v) {
+    y += x * x - A * cos(two_pi<Real>() * x);
+  }
+  return y;
+}
+
+std::atomic<bool> cancel = false;
+
+void ctrl_c_handler(int){
+    cancel = true;
+    std::cout << "Cancellation requested-this could take a second . . ." << std::endl;
+}
+
+int main() {
+  std::cout << "Running random search on Rastrigin function (global minimum = (0,0,...,0))\n";
+  signal(SIGINT, ctrl_c_handler);
+  using ArgType = std::vector<double>;
+  auto params = cma_es_parameters<ArgType>();
+  params.lower_bounds.resize(50, -5.12);
+  params.upper_bounds.resize(50, 5.12);
+  std::random_device rd;
+  std::mt19937_64 gen(rd());
+
+  // By definition, the value of the function which a target value is provided must be <= target_value.
+  double target_value = 1e-3;
+  std::atomic<double> current_minimum_cost;
+  std::cout << "Hit ctrl-C to gracefully terminate the optimization." << std::endl;
+  auto f = [&]() {
+    return cma_es(rastrigin<double>, params, gen, target_value, &cancel, &current_minimum_cost);
+  };
+  auto future = std::async(std::launch::async, f);
+  std::future_status status = future.wait_for(3ms);
+  while (!cancel && (status != std::future_status::ready)) {
+    status = future.wait_for(3ms);
+    std::cout << "Current cost is " << current_minimum_cost << "\r";
+  }
+
+  auto local_minima = future.get();
+  std::cout << "Local minimum is {";
+  for (size_t i = 0; i < local_minima.size() - 1; ++i) {
+    std::cout << local_minima[i] << ", ";
+  }
+  std::cout << local_minima.back() << "}.\n";
+  std::cout << "Final cost: " << current_minimum_cost << "\n";
+}
+#else
+#warning "Signal handling for the random search example only works on Linux and Mac."
+int main() {
+    return 0;
+}
+#endif