diff --git a/.gitignore b/.gitignore index 26f2766687..b498904cf4 100644 --- a/.gitignore +++ b/.gitignore @@ -37,3 +37,4 @@ a.out /test/test_* !/test/test_*.* !/test/test_*/ +__pycache__/ diff --git a/doc/math.qbk b/doc/math.qbk index 4182fe8149..fd79cdd91b 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -739,6 +739,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/multilevel_coordinate_search.qbk] [include optimization/cma_es.qbk] [include optimization/gradient_optimizers.qbk] [endmathpart] [/mathpart optimization Optimization] diff --git a/doc/optimization/multilevel_coordinate_search.qbk b/doc/optimization/multilevel_coordinate_search.qbk new file mode 100644 index 0000000000..84913f2311 --- /dev/null +++ b/doc/optimization/multilevel_coordinate_search.qbk @@ -0,0 +1,153 @@ +[/ +Copyright (c) 2026 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:multilevel_coordinate_search Multilevel Coordinate Search] + +[heading Synopsis] + +`` + #include + + namespace boost::math::optimization { + + template struct multilevel_coordinate_search_parameters { + using Real = typename ArgumentContainer::value_type; + ArgumentContainer lower_bounds; + ArgumentContainer upper_bounds; + ArgumentContainer const * initial_guess = nullptr; + size_t max_function_calls = 0; + size_t max_level = 0; + bool local_search = true; + size_t max_local_search_steps = 50; + size_t line_search_points = 15; + }; + + template + void validate_multilevel_coordinate_search_parameters(multilevel_coordinate_search_parameters const ¶ms); + + template + ArgumentContainer multilevel_coordinate_search(const Func cost_function, + multilevel_coordinate_search_parameters const ¶ms, + std::invoke_result_t target_value + = std::numeric_limits>::quiet_NaN(), + std::atomic *cancellation = nullptr, + std::atomic> *current_minimum_cost = nullptr, + std::vector>> *queries = nullptr); + + } // namespaces +`` + +The `multilevel_coordinate_search` function searches for a global minimum of a function on a box, +using the multilevel coordinate search (MCS) algorithm of Huyer and Neumaier. +MCS splits the box into smaller boxes along one coordinate at a time, each with a base point at which the function has been evaluated. +Every box has a level, which measures how often it has been split; a sweep through the levels splits the box with the smallest function value on each level, +so that large unexplored boxes (global search) and boxes with small function values (local search) are both refined. +Boxes that reach the maximum level start local searches, which build quadratic models of the function by coordinate searches and minimize them over the box. + +Unlike differential evolution, jSO, random search, and CMA-ES, MCS is deterministic: it needs no random number generator and always makes the same function calls. +It is single-threaded, and is at its best on smooth functions in a modest number of dimensions, where it often finds the global minimum in a few hundred function calls. + + +[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 as `lower_bounds`, and each element must be strictly greater than the corresponding element of `lower_bounds`. +* `initial_guess`: An optional point in the box which is used as the base point of the initial box. + It is the first point evaluated. Defaults to the midpoint of the box. +* `max_function_calls`: Defaults to 1000 times the dimension. +* `max_level`: The number of levels /s/[sub max]; boxes that reach this level are no longer split. + Defaults to 5/n/ + 10, where /n/ is the dimension, as in Section 7.3 of the paper. It must not be 1. +* `local_search`: Whether to start local searches from the boxes that reach the maximum level. Defaults to true. +* `max_local_search_steps`: The maximum number of iterations of each local search. Defaults to 50. +* `line_search_points`: The maximum number of function values in a single line search of the local search. Defaults to 15, and must be at least 3. + +`validate_multilevel_coordinate_search_parameters` throws a `std::domain_error` if the bounds or initial guess are invalid, +and a `std::invalid_argument` if `max_level` or `line_search_points` is out of range. +`multilevel_coordinate_search` calls it before the first function evaluation. + +[heading The function] + +`` +template +ArgumentContainer multilevel_coordinate_search(const Func cost_function, + multilevel_coordinate_search_parameters const ¶ms, + std::invoke_result_t target_value + = std::numeric_limits>::quiet_NaN(), + std::atomic *cancellation = nullptr, + std::atomic> *current_minimum_cost = nullptr, + std::vector>> *queries = nullptr) +`` + +Parameters: + +* `cost_function`: The cost function to be minimized. A NaN return value is treated as +[infin]. +* `params`: The parameters to the algorithm as described above. +* `target_value`: 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. + It is checked before every function call. +* `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 `ArgumentContainer` corresponding to the minimum cost found by the optimization. + +[heading Termination] + +As in the paper, there is no convergence criterion. +The computation stops when the `target_value` is reached, when `max_function_calls` calls have been made, when `cancellation` is set, +or when every box has reached `max_level`. +In practice the budget is usually what stops it, so if the minimum value is known, pass it as `target_value`. + +[heading Implementation notes] + +The global part (the initialization procedure, splitting by rank and by expected gain, and the sweeps through the levels) follows Sections 3 and 4 of the paper, +and the local search and the "shopping basket" of local minimizers follow Section 5. +There are two departures from the paper in the local search: +the univariate line searches are our own (extrapolation until the best point is bracketed, then safeguarded parabolic interpolation) rather than the authors' `gls`, +and the quadratic models are minimized over the box by a small bound-constrained quadratic minimizer +(sweeps of exact coordinate minimization, each followed by a Newton step on the free variables via a Cholesky factorization, +or a step along a direction of negative curvature when the factorization fails) rather than MINQ. +Like MINQ, it finds a local minimizer of a possibly indefinite quadratic. + +[h4:examples Examples] + +The six-hump camel function, one of the Dixon-Szego test functions of the paper, has two global minima with value -1.0316284534898774 on the box \[-3, 3\]×\[-2, 2\]: + +`` +#include +#include +#include + +int main() { + using boost::math::optimization::multilevel_coordinate_search; + using boost::math::optimization::multilevel_coordinate_search_parameters; + using ArgType = std::array; + auto camel = [](ArgType const & v) { + double x = v[0]; + double y = v[1]; + return (4 - 2.1*x*x + x*x*x*x/3)*x*x + x*y + (-4 + 4*y*y)*y*y; + }; + multilevel_coordinate_search_parameters params; + params.lower_bounds = {-3, -2}; + params.upper_bounds = {3, 2}; + auto x = multilevel_coordinate_search(camel, params); + std::cout << "f(" << x[0] << ", " << x[1] << ") = " << camel(x) << "\n"; +} +`` + +[h4:references References] + +* W. Huyer and A. Neumaier, ['Global optimization by multilevel coordinate search.] Journal of Global Optimization, vol. 14, pp. 331-355, 1999, doi: 10.1023/A:1008382309369. [@https://arnold-neumaier.at/ms/mcs.pdf Preprint]. + +[endsect] [/section:multilevel_coordinate_search] diff --git a/include/boost/math/optimization/multilevel_coordinate_search.hpp b/include/boost/math/optimization/multilevel_coordinate_search.hpp new file mode 100644 index 0000000000..7115df4085 --- /dev/null +++ b/include/boost/math/optimization/multilevel_coordinate_search.hpp @@ -0,0 +1,1343 @@ +/* + * Copyright Nick Thompson, 2026 + * 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) + */ +#ifndef BOOST_MATH_OPTIMIZATION_MULTILEVEL_COORDINATE_SEARCH_HPP +#define BOOST_MATH_OPTIMIZATION_MULTILEVEL_COORDINATE_SEARCH_HPP + +#include +#ifndef BOOST_MATH_BUILD_MODULE +#include +#endif +#include +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +BOOST_MATH_NAMESPACE_BEGIN namespace optimization { + +// Huyer, W., Neumaier, A. (1999). Global optimization by multilevel coordinate search. +// Journal of Global Optimization, 14, 331-355. +// See: https://arnold-neumaier.at/ms/mcs.pdf +// Section, step and equation numbers in the comments below refer to this paper. + +BOOST_MATH_EXPORT template struct multilevel_coordinate_search_parameters { + using Real = typename ArgumentContainer::value_type; + ArgumentContainer lower_bounds; + ArgumentContainer upper_bounds; + // Used as the base point x^0 of the root box; the default is the midpoint of the box. + ArgumentContainer const *initial_guess = nullptr; + // 0 means 1000 times the dimension. + size_t max_function_calls = 0; + // The number of levels s_max; 0 means 5n + 10, the value used in Section 7.3. + size_t max_level = 0; + // Start local searches from the base points of the boxes that reach level s_max (Section 5). + bool local_search = true; + // The limit on visits to Step 3 per local search (Section 7.3). + size_t max_local_search_steps = 50; + // smaxls in Section 7.3: the most function values in one line search of the local search. + size_t line_search_points = 15; +}; + +BOOST_MATH_EXPORT template +void validate_multilevel_coordinate_search_parameters(multilevel_coordinate_search_parameters const ¶ms) { + std::ostringstream oss; + detail::validate_bounds(params.lower_bounds, params.upper_bounds); + for (size_t i = 0; i < params.lower_bounds.size(); ++i) { + if (!(params.lower_bounds[i] < params.upper_bounds[i])) { + oss << __FILE__ << ":" << __LINE__ << ":" << __func__; + oss << ": MCS requires lower_bounds[i] < upper_bounds[i], but at i = " << i << " the bounds are " + << params.lower_bounds[i] << " and " << params.upper_bounds[i] << "."; + throw std::domain_error(oss.str()); + } + } + if (params.initial_guess) { + detail::validate_initial_guess(*params.initial_guess, params.lower_bounds, params.upper_bounds); + } + if (params.max_level == 1) { + oss << __FILE__ << ":" << __LINE__ << ":" << __func__; + oss << ": There must be at least two levels."; + throw std::invalid_argument(oss.str()); + } + if (params.line_search_points < 3) { + oss << __FILE__ << ":" << __LINE__ << ":" << __func__; + oss << ": A line search needs at least 3 points, but line_search_points = " << params.line_search_points << "."; + throw std::invalid_argument(oss.str()); + } +} + +namespace detail { + +// Minimizes g^T h + h^T G h/2 over lo <= h <= hi, where lo <= 0 <= hi and G is symmetric (row-major n x n) but may be +// indefinite. This stands in for MINQ: sweeps of exact coordinate minimization, each followed by a Newton step on the +// free variables when the Hessian restricted to them is positive definite. It finds a local minimizer, as MINQ does. +template +std::vector minimize_quadratic_on_box(std::vector const &g, std::vector const &G, + std::vector const &lo, std::vector const &hi) { + using std::abs; + using std::isfinite; + using std::sqrt; + size_t const n = g.size(); + std::vector h(n, Real(0)); + for (size_t i = 0; i < n; ++i) { + if (!isfinite(g[i])) { + return h; + } + } + for (auto const &Gij : G) { + if (!isfinite(Gij)) { + return h; + } + } + // r = g + G h, the gradient at h: + std::vector r = g; + auto add_column = [&](size_t j, Real t) { + for (size_t i = 0; i < n; ++i) { + r[i] += t * G[i * n + j]; + } + }; + // Moves h_j by t within [lo_j, hi_j], keeping r = g + G h: + auto move = [&](size_t j, Real t) { + Real const old = h[j]; + h[j] = (std::min)((std::max)(h[j] + t, lo[j]), hi[j]); + add_column(j, h[j] - old); + }; + auto value = [&]() { + Real q = 0; + for (size_t i = 0; i < n; ++i) { + q += h[i] * (g[i] + r[i]) / 2; + } + return q; + }; + std::vector free_set; + std::vector L, p; + Real q_old = 0; + for (size_t sweep = 0; sweep < 10 * n + 10; ++sweep) { + for (size_t i = 0; i < n; ++i) { + // Minimize r_i t + G_ii t^2/2 over h_i + t in [lo_i, hi_i]: + Real const a = G[i * n + i]; + Real const tlo = lo[i] - h[i]; + Real const thi = hi[i] - h[i]; + Real t = 0; + if (a > 0) { + t = (std::min)((std::max)(-r[i] / a, tlo), thi); + } + else { + // Not convex along coordinate i, so the minimum is at an end of the interval (or t = 0 if neither is lower): + Real const qlo = r[i] * tlo + a * tlo * tlo / 2; + Real const qhi = r[i] * thi + a * thi * thi / 2; + if (qlo < 0 && qlo <= qhi) { + t = tlo; + } + else if (qhi < 0) { + t = thi; + } + } + if (t != 0) { + move(i, t); + } + } + // Newton step on the free variables, by a hand-rolled Cholesky factorization: + free_set.clear(); + for (size_t i = 0; i < n; ++i) { + if (lo[i] < h[i] && h[i] < hi[i]) { + free_set.push_back(i); + } + } + size_t const m = free_set.size(); + bool positive_definite = m > 0; + L.assign(m * m, Real(0)); + for (size_t j = 0; j < m && positive_definite; ++j) { + Real d = G[free_set[j] * n + free_set[j]]; + for (size_t k = 0; k < j; ++k) { + d -= L[j * m + k] * L[j * m + k]; + } + if (!(d > 0)) { + positive_definite = false; + // w = (-G_11^{-1} G_1j, 1, 0, ...) has curvature w^T G w = d <= 0, where G_11 = L_11 L_11^T is the leading + // block already factored. Go downhill along it as far as the box allows. + p.assign(m, Real(0)); + p[j] = 1; + for (size_t k = j; k-- > 0;) { + Real s = L[j * m + k]; + for (size_t q = k + 1; q < j; ++q) { + s += L[q * m + k] * p[q]; + } + p[k] = -s / L[k * m + k]; + } + Real slope = 0; + for (size_t i = 0; i <= j; ++i) { + slope += r[free_set[i]] * p[i]; + } + if (slope > 0) { + for (auto &pi : p) { + pi = -pi; + } + } + if (d < 0 || slope != 0) { + Real t = std::numeric_limits::infinity(); + for (size_t i = 0; i <= j; ++i) { + size_t const k = free_set[i]; + if (p[i] > 0) { + t = (std::min)(t, (hi[k] - h[k]) / p[i]); + } + else if (p[i] < 0) { + t = (std::min)(t, (lo[k] - h[k]) / p[i]); + } + } + if (t > 0 && t < std::numeric_limits::infinity()) { + for (size_t i = 0; i <= j; ++i) { + size_t const k = free_set[i]; + move(k, t * p[i]); + } + } + } + break; + } + L[j * m + j] = sqrt(d); + for (size_t i = j + 1; i < m; ++i) { + Real s = G[free_set[i] * n + free_set[j]]; + for (size_t k = 0; k < j; ++k) { + s -= L[i * m + k] * L[j * m + k]; + } + L[i * m + j] = s / L[j * m + j]; + } + } + if (positive_definite) { + // Solve L L^T p = -r_F: + p.assign(m, Real(0)); + for (size_t i = 0; i < m; ++i) { + Real s = -r[free_set[i]]; + for (size_t k = 0; k < i; ++k) { + s -= L[i * m + k] * p[k]; + } + p[i] = s / L[i * m + i]; + } + for (size_t i = m; i-- > 0;) { + Real s = p[i]; + for (size_t k = i + 1; k < m; ++k) { + s -= L[k * m + i] * p[k]; + } + p[i] = s / L[i * m + i]; + } + // The largest step t <= 1 along p that stays in the box: + Real t = 1; + for (size_t i = 0; i < m; ++i) { + size_t const j = free_set[i]; + if (p[i] > 0) { + t = (std::min)(t, (hi[j] - h[j]) / p[i]); + } + else if (p[i] < 0) { + t = (std::min)(t, (lo[j] - h[j]) / p[i]); + } + } + if (t > 0) { + for (size_t i = 0; i < m; ++i) { + size_t const j = free_set[i]; + move(j, t * p[i]); + } + } + } + Real const q = value(); + if (!(q < q_old - 4 * std::numeric_limits::epsilon() * abs(q))) { + break; + } + q_old = q; + } + return h; +} + +// The state of an MCS run. multilevel_coordinate_search() below drives it. It is separate so that the boxes can be inspected between sweeps. +template class multilevel_coordinate_search_engine { +public: + using Real = typename ArgumentContainer::value_type; + using ResultType = std::invoke_result_t; + + // A box B[x, y] (Section 3) with base point x and opposite point y. + struct box { + std::vector x; + std::vector y; + Real f; + // 0 for a box that has been split: + size_t level; + // The number of times each coordinate has been split in the history of the box: + std::vector splits; + // For each coordinate, two more points (the closest in the history of the box) and the differences of their + // function values from that of the base point they were split from, for the separable model of Section 4.2: + std::vector xi1, df1, xi2, df2; + }; + + multilevel_coordinate_search_engine(Func const &cost_function, multilevel_coordinate_search_parameters const ¶ms, ResultType target_value, + std::atomic *cancellation, std::atomic *current_minimum_cost, + std::vector> *queries) + : cost_function_(cost_function), params_(params), target_value_(target_value), cancellation_(cancellation), + current_minimum_cost_(current_minimum_cost), queries_(queries) { + validate_multilevel_coordinate_search_parameters(params); + n_ = params.lower_bounds.size(); + u_.assign(params.lower_bounds.begin(), params.lower_bounds.end()); + v_.assign(params.upper_bounds.begin(), params.upper_bounds.end()); + smax_ = params.max_level ? params.max_level : 5 * n_ + 10; + max_function_calls_ = params.max_function_calls ? params.max_function_calls : 1000 * n_; + heaps_.resize(smax_); + } + + // The initialization procedure of Section 3, with the initialization list {u_i, x^0_i, v_i}. + void initialize() { + std::vector x0(n_); + for (size_t i = 0; i < n_; ++i) { + x0[i] = params_.initial_guess ? Real((*params_.initial_guess)[i]) : (u_[i] + v_[i]) / 2; + } + init_list_.resize(n_); + init_f_.resize(n_); + l0_.resize(n_); + for (size_t i = 0; i < n_; ++i) { + auto &list = init_list_[i]; + list = {u_[i], x0[i], v_[i]}; + if (x0[i] == u_[i] || x0[i] == v_[i]) { + list = {u_[i], (u_[i] + v_[i]) / 2, v_[i]}; + } + l0_[i] = static_cast(std::find(list.begin(), list.end(), x0[i]) - list.begin()); + } + box root; + root.x = x0; + // The best point stays x0 if the run is over before f is ever called: + xbest_ = x0; + root.f = evaluate(x0); + root.y.resize(n_); + for (size_t i = 0; i < n_; ++i) { + root.y[i] = (x0[i] - u_[i] > v_[i] - x0[i]) ? u_[i] : v_[i]; + } + root.level = 1; + root.splits.assign(n_, 0); + root.xi1.assign(n_, std::numeric_limits::quiet_NaN()); + root.df1 = root.xi2 = root.df2 = root.xi1; + boxes_.push_back(std::move(root)); + size_t current = 0; + for (size_t i = 0; i < n_ && !done(); ++i) { + init_f_[i] = split_by_list(current, i, ¤t); + } + if (!done()) { + push(current); + } + f0_ = fbest_; + // The variability of f along each coordinate, from the quadratic interpolants through consecutive triples: + std::vector variability(n_, Real(0)); + for (size_t i = 0; i < n_ && !done(); ++i) { + auto const &a = init_list_[i]; + auto const &f = init_f_[i]; + Real lo = std::numeric_limits::infinity(); + Real hi = -lo; + for (size_t l = 0; l + 2 < a.size(); ++l) { + Real c[3]; + quadratic_through(a[l], f[l], a[l + 1], f[l + 1], a[l + 2], f[l + 2], c); + // The range over [a_l, a_{l+2}]: the ends, and the stationary point when it lies between them. + Real const stationary = c[2] != 0 ? -c[1] / (2 * c[2]) : a[l]; + for (Real t : {a[l], a[l + 2], std::clamp(stationary, a[l], a[l + 2])}) { + Real const q = c[0] + t * (c[1] + t * c[2]); + lo = (std::min)(lo, q); + hi = (std::max)(hi, q); + } + } + using std::isnan; + variability[i] = isnan(hi - lo) ? std::numeric_limits::infinity() : hi - lo; + } + rank_.resize(n_); + std::vector order(n_); + for (size_t i = 0; i < n_; ++i) { + order[i] = i; + } + std::stable_sort(order.begin(), order.end(), [&](size_t a, size_t b) { return variability[a] > variability[b]; }); + for (size_t r = 0; r < n_; ++r) { + rank_[order[r]] = r; + } + } + + // One sweep through the levels (Section 3), then the local searches of Section 5.4. + // Returns false once the run is over. + bool sweep() { + if (done()) { + return false; + } + bool any = false; + for (size_t s = 1; s < smax_ && !done(); ++s) { + if (heaps_[s].empty()) { + continue; + } + any = true; + std::pop_heap(heaps_[s].begin(), heaps_[s].end(), box_greater{this}); + size_t const b = heaps_[s].back(); + heaps_[s].pop_back(); + process(b); + } + if (params_.local_search && !done()) { + shopping_basket(); + } + candidates_.clear(); + return any && !done(); + } + + std::vector const &boxes() const { return boxes_; } + // The local minimizers found by local search (the shopping basket of Section 5.4): + std::vector, Real>> const &basket() const { return basket_; } + std::vector const &best_point() const { return xbest_; } + Real best_value() const { return fbest_; } + size_t function_calls() const { return function_calls_; } + size_t max_level() const { return smax_; } + + ArgumentContainer best_argument() const { return to_argument(xbest_); } + +private: + struct box_greater { + multilevel_coordinate_search_engine const *engine; + bool operator()(size_t a, size_t b) const { + Real const fa = engine->boxes_[a].f; + Real const fb = engine->boxes_[b].f; + return fa > fb || (fa == fb && a > b); + } + }; + + static Real golden() { return Real(0.6180339887498948482045868343656381177203L); } + + bool done() const { + using std::isnan; + if (cancellation_ && *cancellation_) { + return true; + } + if (!isnan(target_value_) && fbest_ <= Real(target_value_)) { + return true; + } + return function_calls_ >= max_function_calls_; + } + + ArgumentContainer to_argument(std::vector const &x) const { + ArgumentContainer arg; + if constexpr (detail::has_resize_v) { + arg.resize(n_); + } + for (size_t i = 0; i < n_; ++i) { + arg[i] = x[i]; + } + return arg; + } + + // Calls the cost function (NaN counts as +infinity), or returns +infinity once the run is over. + Real evaluate(std::vector const &x) { + using std::isnan; + if (done()) { + return std::numeric_limits::infinity(); + } + auto const arg = to_argument(x); + ResultType const value = cost_function_(arg); + ++function_calls_; + if (queries_) { + queries_->push_back(std::make_pair(arg, value)); + } + Real f = static_cast(value); + if (isnan(f)) { + f = std::numeric_limits::infinity(); + } + if (f < fbest_ || function_calls_ == 1) { + fbest_ = f; + xbest_ = x; + if (current_minimum_cost_) { + *current_minimum_cost_ = static_cast(f); + } + } + return f; + } + + // The coefficients c of the quadratic c0 + c1 t + c2 t^2 through three points with distinct abscissae. + static void quadratic_through(Real a, Real fa, Real b, Real fb, Real c, Real fc, Real *coef) { + Real const dab = (fb - fa) / (b - a); + Real const dabc = ((fc - fb) / (c - b) - dab) / (c - a); + coef[2] = dabc; + coef[1] = dab - dabc * (a + b); + coef[0] = fa - a * (dab - dabc * b); + } + + // The minimizer over [lo, hi] (in either order) of the quadratic with coefficients coef. + static Real vertex(Real const *coef, Real lo, Real hi) { + if (lo > hi) { + std::swap(lo, hi); + } + auto q = [&](Real t) { return coef[0] + t * (coef[1] + t * coef[2]); }; + Real best = q(lo) <= q(hi) ? lo : hi; + if (coef[2] > 0) { + Real const t = -coef[1] / (2 * coef[2]); + if (lo < t && t < hi && q(t) < q(best)) { + best = t; + } + } + return best; + } + + // subint of (4), as in the authors' Matlab code: the second case is for 1000|x| >= 1. + static Real subint(Real x, Real y) { + using std::abs; + Real const s = y < 0 ? Real(-1) : Real(1); + if (1000 * abs(x) < 1) { + return abs(y) > 1000 ? s : y; + } + return abs(y) > 1000 * abs(x) ? 10 * s * abs(x) : y; + } + + void push(size_t b) { + if (boxes_[b].level >= smax_) { + boxes_[b].level = smax_; + candidates_.push_back(b); + return; + } + heaps_[boxes_[b].level].push_back(b); + std::push_heap(heaps_[boxes_[b].level].begin(), heaps_[boxes_[b].level].end(), box_greater{this}); + } + + // A child of box `parent` split along coordinate i, with base point coordinate xi (value f), the other end `end` of + // its ith interval, and the points `others` (coordinate, value) of this split, nearest first. The differences kept + // for the other coordinates stay relative to the base points they were measured from: the model is separable. + size_t add_child(size_t parent, size_t i, Real xi, Real f, Real end, size_t level, + std::vector> others) { + box child = boxes_[parent]; + child.x[i] = xi; + child.f = f; + child.y[i] = end; + child.level = (std::min)(level, smax_); + ++child.splits[i]; + // The two points closest in the history, other than the base point: + using std::abs; + using std::isnan; + if (!isnan(boxes_[parent].xi1[i])) { + others.emplace_back(boxes_[parent].xi1[i], boxes_[parent].f + boxes_[parent].df1[i]); + } + if (!isnan(boxes_[parent].xi2[i])) { + others.emplace_back(boxes_[parent].xi2[i], boxes_[parent].f + boxes_[parent].df2[i]); + } + std::stable_sort(others.begin(), others.end(), + [&](auto const &a, auto const &b) { return abs(a.first - xi) < abs(b.first - xi); }); + size_t k = 0; + for (auto const &o : others) { + if (o.first == xi || (k == 1 && o.first == child.xi1[i])) { + continue; + } + (k == 0 ? child.xi1[i] : child.xi2[i]) = o.first; + (k == 0 ? child.df1[i] : child.df2[i]) = o.second - f; + if (++k == 2) { + break; + } + } + boxes_.push_back(std::move(child)); + return boxes_.size() - 1; + } + + // Splits box b along coordinate i (never split before) at the initialization list and the golden section points + // between them (Section 3). Returns the function values at the list. When `current` is given (the initialization + // procedure), it receives the child to split next, and that child is not pushed. + std::vector split_by_list(size_t b, size_t i, size_t *current = nullptr) { + auto const &a = init_list_[i]; + size_t const L = a.size(); + std::vector f(L); + std::vector x = boxes_[b].x; + for (size_t l = 0; l < L; ++l) { + if (l == l0_[i]) { + f[l] = boxes_[b].f; + } + else { + x[i] = a[l]; + f[l] = evaluate(x); + } + } + if (done()) { + return f; + } + size_t const s = boxes_[b].level; + boxes_[b].level = 0; + auto others = [&](size_t l) { + std::vector> o; + std::vector idx; + for (size_t k = 0; k < L; ++k) { + if (k != l) { + idx.push_back(k); + } + } + std::stable_sort(idx.begin(), idx.end(), [&](size_t p, size_t q) { + using std::abs; + return abs(a[p] - a[l]) < abs(a[q] - a[l]); + }); + for (size_t k : idx) { + o.emplace_back(a[k], f[k]); + } + return o; + }; + size_t lbest = 0; + for (size_t l = 1; l < L; ++l) { + if (f[l] < f[lbest]) { + lbest = l; + } + } + // The children whose base point is the best point, left and right of it: + size_t left = 0, right = 0; + bool has_left = false, has_right = false; + std::vector children; + if (a[0] > u_[i]) { + children.push_back(add_child(b, i, a[0], f[0], u_[i], s + 1, others(0))); + if (lbest == 0) { + left = children.back(); + has_left = true; + } + } + for (size_t l = 1; l < L; ++l) { + // The part with the smaller function value gets the larger fraction: + bool const left_larger = f[l - 1] <= f[l]; + Real const z = a[l - 1] + (left_larger ? golden() : golden() * golden()) * (a[l] - a[l - 1]); + children.push_back(add_child(b, i, a[l - 1], f[l - 1], z, left_larger ? s + 1 : s + 2, others(l - 1))); + if (l - 1 == lbest) { + right = children.back(); + has_right = true; + } + children.push_back(add_child(b, i, a[l], f[l], z, left_larger ? s + 2 : s + 1, others(l))); + if (l == lbest) { + left = children.back(); + has_left = true; + } + } + if (a[L - 1] < v_[i]) { + children.push_back(add_child(b, i, a[L - 1], f[L - 1], v_[i], s + 1, others(L - 1))); + if (lbest == L - 1) { + right = children.back(); + has_right = true; + } + } + size_t next = has_left ? left : right; + if (has_left && has_right) { + // The one containing the minimizer of the quadratic model through the best point and its neighbors: + size_t const l = (std::min)((std::max)(lbest, size_t(1)), L - 2); + Real c[3]; + quadratic_through(a[l - 1], f[l - 1], a[l], f[l], a[l + 1], f[l + 1], c); + next = vertex(c, a[0], a[L - 1]) < a[lbest] ? left : right; + } + for (size_t c : children) { + if (current && c == next) { + continue; + } + push(c); + } + if (current) { + *current = next; + } + return f; + } + + // Splits box b along coordinate i at z (new value fz) and at the golden section point between x_i and z, giving + // two parts, or three when z != y_i (Sections 4.1 and 4.2). + void split_at(size_t b, size_t i, Real z, Real fz, bool by_rank) { + using std::abs; + box const &B = boxes_[b]; + Real const xi = B.x[i]; + Real const yi = B.y[i]; + Real const f = B.f; + size_t const s = B.level; + bool const near_larger = f <= fz; + Real const fraction = near_larger ? golden() : golden() * golden(); + Real const g = xi + fraction * (z - xi); + if (g == xi || g == z) { + // Too narrow to split in floating point; let the box move up a level instead. + boxes_[b].level = s + 1; + push(b); + return; + } + boxes_[b].level = 0; + std::vector> from_x{{z, fz}}; + std::vector> from_z{{xi, f}}; + push(add_child(b, i, xi, f, g, near_larger ? s + 1 : s + 2, from_x)); + push(add_child(b, i, z, fz, g, near_larger ? s + 2 : s + 1, from_z)); + if (z != yi) { + size_t level = s + 1; + if (!by_rank && !(abs(yi - z) > (std::min)(abs(g - xi), abs(z - g)))) { + level = s + 2; + } + push(add_child(b, i, z, fz, yi, level, from_z)); + } + } + + // Step 2 of a sweep (Section 3) for box b: split it (Section 4), or move it up a level. + void process(size_t b) { + using std::abs; + using std::isfinite; + box const &B = boxes_[b]; + size_t const s = B.level; + size_t const m = *std::min_element(B.splits.begin(), B.splits.end()); + if (s > 2 * n_ * (m + 1)) { + // Splitting by rank (Section 4.1): + size_t i = n_; + for (size_t j = 0; j < n_; ++j) { + if (B.splits[j] == m && (i == n_ || rank_[j] < rank_[i])) { + i = j; + } + } + if (m == 0) { + split_by_list(b, i); + return; + } + Real const z = B.x[i] + 2 * (subint(B.x[i], B.y[i]) - B.x[i]) / 3; + std::vector x = B.x; + x[i] = z; + Real const fz = evaluate(x); + if (!done()) { + split_at(b, i, z, fz, true); + } + return; + } + // Splitting by expected gain (Section 4.2): + Real best_gain = std::numeric_limits::infinity(); + size_t i = n_; + Real z = 0; + for (size_t j = 0; j < n_; ++j) { + Real gain; + Real zj = 0; + if (B.splits[j] == 0) { + gain = *std::min_element(init_f_[j].begin(), init_f_[j].end()) - init_f_[j][l0_[j]]; + } + else { + Real const xj = B.x[j]; + Real const outer = subint(xj, B.y[j]); + Real const inner = xj + (outer - xj) / 10; + // e_j(t) = alpha (t - x_j) + beta (t - x_j)^2 through the two more points: + Real const h1 = B.xi1[j] - xj; + Real const h2 = B.xi2[j] - xj; + Real const d1 = B.df1[j] / h1; + Real const d2 = B.df2[j] / h2; + Real const beta = (d1 - d2) / (h1 - h2); + Real const alpha = d1 - beta * h1; + Real const c[3] = {0, alpha, beta}; + Real const t = vertex(c, inner - xj, outer - xj); + gain = alpha * t + beta * t * t; + zj = xj + t; + if (!isfinite(gain) || outer == xj) { + gain = std::numeric_limits::infinity(); + } + } + if (gain < best_gain) { + best_gain = gain; + i = j; + z = zj; + } + } + // (5): split only if the model expects to beat the best value so far. + if (i < n_ && B.f + best_gain < fbest_) { + if (B.splits[i] == 0) { + split_by_list(b, i); + return; + } + std::vector x = B.x; + x[i] = z; + Real const fz = evaluate(x); + if (!done()) { + split_at(b, i, z, fz, false); + } + return; + } + boxes_[b].level = s + 1; + push(b); + } + + // --------------------------------------------------------------------------------------------------------------- + // Local search (Section 5). + + std::vector along(std::vector const &x, std::vector const &p, Real alpha) const { + std::vector y(n_); + for (size_t i = 0; i < n_; ++i) { + y[i] = std::clamp(x[i] + alpha * p[i], u_[i], v_[i]); + } + return y; + } + + // A univariate line search along x + alpha p, alpha in [amin, amax], which stands in for gls. `samples` holds the + // (alpha, f) pairs known so far, including alpha = 0, and receives the new ones, sorted by alpha. With only one + // sample, the first step is `step`. It extrapolates until the best point is bracketed, then interpolates + // parabolically (golden section steps as a safeguard), using at most `max_points` points in all. + void line_search(std::vector const &x, std::vector const &p, Real amin, Real amax, + std::vector> &samples, size_t max_points, Real step) { + using std::abs; + using std::sqrt; + Real const tol = sqrt(std::numeric_limits::epsilon()); + auto by_alpha = [](auto const &a, auto const &b) { return a.first < b.first; }; + std::sort(samples.begin(), samples.end(), by_alpha); + while (samples.size() < max_points && !done()) { + size_t const N = samples.size(); + size_t ib = 0; + for (size_t k = 1; k < N; ++k) { + if (samples[k].second < samples[ib].second) { + ib = k; + } + } + Real const scale = (std::max)(abs(samples.front().first), abs(samples.back().first)); + Real a; + if (N == 1) { + a = std::clamp(step, amin, amax); + if (a == samples[0].first) { + a = std::clamp(-step, amin, amax); + } + } + else if (ib == 0 && samples[0].first > amin) { + a = (std::max)(amin, samples[0].first - 2 * (samples[1].first - samples[0].first)); + } + else if (ib == N - 1 && samples[N - 1].first < amax) { + a = (std::min)(amax, samples[N - 1].first + 2 * (samples[N - 1].first - samples[N - 2].first)); + } + else { + // Bracketed, or the best point is on the bound: refine around it. + size_t const l = ib == 0 ? 0 : (ib == N - 1 ? N - 2 : ib - 1); + size_t const r = ib == 0 ? 1 : (ib == N - 1 ? N - 1 : ib + 1); + Real const al = samples[l].first; + Real const ar = samples[r].first; + Real const ab = samples[ib].first; + if (ar - al <= tol * (1 + scale)) { + break; + } + a = std::numeric_limits::quiet_NaN(); + if (ib != 0 && ib != N - 1) { + Real c[3]; + quadratic_through(al, samples[l].second, ab, samples[ib].second, ar, samples[r].second, c); + if (c[2] > 0) { + a = -c[1] / (2 * c[2]); + } + } + Real const min_gap = tol * (1 + scale); + if (!(al + min_gap < a && a < ar - min_gap && abs(a - ab) > min_gap)) { + // Golden section step into the larger side of the best point: + a = ab - al > ar - ab ? ab - (1 - golden()) * (ab - al) : ab + (1 - golden()) * (ar - ab); + } + } + bool known = false; + for (auto const &sample : samples) { + if (abs(sample.first - a) <= tol * (1 + scale) / 4) { + known = true; + } + } + if (known) { + break; + } + samples.emplace_back(a, evaluate(along(x, p, a))); + std::sort(samples.begin(), samples.end(), by_alpha); + } + } + + // The local quadratic model of Section 5.1 around xb: f(xb + h) ~ fb + g^T h + h^T G h/2. + struct model { + std::vector xb; + Real fb; + std::vector g; + std::vector G; + }; + + // Triple search (Section 5.1) over the coordinates `coords`, in order. For each coordinate i, `triple(i, M)` gives + // the two other values t1, t2 of the ith coordinate (xb_i being the third) and f at xb with x_i changed to them. + // It returns `skip` for a coordinate with no usable triple, and `stop` to end the search. + // Without `full`, the off-diagonal entries of G are those of the previous model. + enum class triple_status { ok, skip, stop }; + template + void triple_search(model &M, std::vector const &coords, bool full, Triple triple) { + size_t const n = n_; + // The triple of each coordinate: its value at the time and the two others. + std::vector t0(n), t1(n), t2(n); + std::vector done_coords; + for (size_t i : coords) { + Real f1, f2; + triple_status const status = triple(i, M, t1[i], f1, t2[i], f2); + if (status == triple_status::stop) { + return; + } + if (status == triple_status::skip) { + continue; + } + t0[i] = M.xb[i]; + std::vector xnew = M.xb; + Real fnew = M.fb; + Real const h1 = t1[i] - M.xb[i]; + Real const h2 = t2[i] - M.xb[i]; + Real const d1 = (f1 - M.fb) / h1; + Real const d2 = (f2 - M.fb) / h2; + M.G[i * n + i] = 2 * (d1 - d2) / (h1 - h2); + M.g[i] = d1 - M.G[i * n + i] * h1 / 2; + if (f1 < fnew || f2 < fnew) { + fnew = f1 < f2 ? f1 : f2; + xnew[i] = f1 < f2 ? t1[i] : t2[i]; + } + if (full) { + // x^{ik}: the ith coordinate at the better of t1, t2, the kth at the one with the smaller model value. + Real const ti = f1 < f2 ? t1[i] : t2[i]; + Real const hi = ti - M.xb[i]; + for (size_t k : done_coords) { + // The two values of coordinate k's triple other than xb_k, which may have moved within the triple: + Real ck[3] = {t0[k], t1[k], t2[k]}; + Real candidates[2]; + size_t nc = 0; + for (Real c : ck) { + if (c != M.xb[k] && nc < 2 && (nc == 0 || c != candidates[0])) { + candidates[nc++] = c; + } + } + if (nc == 0) { + continue; + } + auto qk = [&](Real c) { + Real const h = c - M.xb[k]; + return M.g[k] * h + M.G[k * n + k] * h * h / 2; + }; + Real const tk = (nc == 1 || qk(candidates[0]) <= qk(candidates[1])) ? candidates[0] : candidates[1]; + Real const hk = tk - M.xb[k]; + std::vector x = M.xb; + x[i] = ti; + x[k] = tk; + Real const fik = evaluate(x); + if (done()) { + return; + } + Real const Gik = (fik - M.fb - M.g[k] * hk - M.g[i] * hi - M.G[k * n + k] * hk * hk / 2 - + M.G[i * n + i] * hi * hi / 2) / (hk * hi); + M.G[i * n + k] = M.G[k * n + i] = Gik; + if (fik < fnew) { + fnew = fik; + xnew = x; + } + } + } + done_coords.push_back(i); + if (fnew < M.fb) { + // Re-expand the model around the new best point: + std::vector dx(n); + for (size_t l = 0; l < n; ++l) { + dx[l] = xnew[l] - M.xb[l]; + } + for (size_t k = 0; k < n; ++k) { + for (size_t l = 0; l < n; ++l) { + M.g[k] += M.G[k * n + l] * dx[l]; + } + } + M.xb = xnew; + M.fb = fnew; + } + } + } + + // The local search of Section 5.3 from x (value f). Returns the best point found and its value. + std::pair, Real> local_search(std::vector x, Real f) { + using std::abs; + size_t const n = n_; + std::vector x0(n); + for (size_t i = 0; i < n; ++i) { + x0[i] = std::clamp(Real(0), u_[i], v_[i]); + } + std::vector d(n); + for (size_t i = 0; i < n; ++i) { + d[i] = (std::min)((std::min)(v_[i] - x[i], x[i] - u_[i]), Real(0.25) * (1 + abs(x[i] - x0[i]))); + } + model M{x, f, std::vector(n, Real(0)), std::vector(n * n, Real(0))}; + std::vector all(n); + for (size_t i = 0; i < n; ++i) { + all[i] = i; + } + // Step 1: a full triple search on the triples found by line searches along the coordinates (Section 5.2). + bool first = true; + bool failed = false; + triple_search(M, all, true, [&](size_t i, model &Mi, Real &t1, Real &f1, Real &t2, Real &f2) { + std::vector e(n, Real(0)); + e[i] = 1; + Real const step = d[i] > 0 ? d[i] : Real(0.25) * (1 + abs(Mi.xb[i] - x0[i])); + std::vector> samples{{Real(0), Mi.fb}}; + line_search(Mi.xb, e, u_[i] - Mi.xb[i], v_[i] - Mi.xb[i], samples, 6, step); + if (samples.size() < 3) { + failed = true; + return triple_status::stop; + } + size_t ib = 0, i0 = 0; + for (size_t k = 0; k < samples.size(); ++k) { + if (samples[k].second < samples[ib].second) { + ib = k; + } + if (samples[k].first == 0) { + i0 = k; + } + } + size_t const N = samples.size(); + size_t pick[2]; + if (first || ib == i0) { + // The best point and its nearest neighbors on both sides, or the two nearest on one side: + size_t const l = ib == 0 ? 0 : (ib == N - 1 ? N - 3 : ib - 1); + size_t trio[3] = {l, l + 1, l + 2}; + size_t k = 0; + for (size_t c : trio) { + if (c != ib) { + pick[k++] = c; + } + } + if (first) { + // Nothing has been fitted yet, so the model can move to the best point. + Mi.xb[i] += samples[ib].first; + Mi.fb = samples[ib].second; + Real const shift = samples[ib].first; + t1 = Mi.xb[i] + samples[pick[0]].first - shift; + t2 = Mi.xb[i] + samples[pick[1]].first - shift; + f1 = samples[pick[0]].second; + f2 = samples[pick[1]].second; + first = false; + return triple_status::ok; + } + } + else { + // The best point, the start of the line search, and the best point's neighbor on the other side: + pick[0] = ib; + if (ib > i0) { + pick[1] = ib + 1 < N ? ib + 1 : (i0 > 0 ? i0 - 1 : i0 + 1); + } + else { + pick[1] = ib > 0 ? ib - 1 : (i0 + 1 < N ? i0 + 1 : i0 - 1); + } + if (pick[1] == ib || pick[1] == i0) { + for (size_t k = 0; k < N; ++k) { + if (k != ib && k != i0) { + pick[1] = k; + break; + } + } + } + } + t1 = Mi.xb[i] + samples[pick[0]].first; + t2 = Mi.xb[i] + samples[pick[1]].first; + f1 = samples[pick[0]].second; + f2 = samples[pick[1]].second; + return triple_status::ok; + }); + if (failed || done()) { + return {M.xb, M.fb}; + } + // Step 2 starts from the point Step 1 found: + for (size_t i = 0; i < n; ++i) { + d[i] = (std::min)((std::min)(v_[i] - M.xb[i], M.xb[i] - u_[i]), Real(0.25) * (1 + abs(M.xb[i] - x0[i]))); + } + using std::cbrt; + Real const delta = cbrt(std::numeric_limits::epsilon()); + // Steps 2 and 6: minimize the model over the box [-d, d] intersected with [u - x, v - x], then search along the + // minimizer. Returns r, the ratio of the actual to the predicted gain. + auto quadratic_step = [&](std::vector const &free) { + std::vector lo(n), hi(n); + for (size_t i = 0; i < n; ++i) { + lo[i] = free[i] ? (std::max)(-d[i], u_[i] - M.xb[i]) : Real(0); + hi[i] = free[i] ? (std::min)(d[i], v_[i] - M.xb[i]) : Real(0); + } + std::vector p = minimize_quadratic_on_box(M.g, M.G, lo, hi); + Real predicted = 0; + bool moves = false; + for (size_t i = 0; i < n; ++i) { + Real Gp = 0; + for (size_t j = 0; j < n; ++j) { + Gp += M.G[i * n + j] * p[j]; + } + predicted += p[i] * (M.g[i] + Gp / 2); + moves = moves || p[i] != 0; + } + if (!moves || !(predicted < 0)) { + return Real(0); + } + Real amin = -std::numeric_limits::infinity(); + Real amax = std::numeric_limits::infinity(); + for (size_t i = 0; i < n; ++i) { + if (p[i] > 0) { + amax = (std::min)(amax, (v_[i] - M.xb[i]) / p[i]); + amin = (std::max)(amin, (u_[i] - M.xb[i]) / p[i]); + } + else if (p[i] < 0) { + amax = (std::min)(amax, (u_[i] - M.xb[i]) / p[i]); + amin = (std::max)(amin, (v_[i] - M.xb[i]) / p[i]); + } + } + amax = (std::max)(amax, Real(1)); + amin = (std::min)(amin, Real(0)); + Real const fold = M.fb; + std::vector> samples{{Real(0), fold}, {Real(1), evaluate(along(M.xb, p, 1))}}; + line_search(M.xb, p, amin, amax, samples, params_.line_search_points, Real(1)); + size_t ib = 0; + for (size_t k = 1; k < samples.size(); ++k) { + if (samples[k].second < samples[ib].second) { + ib = k; + } + } + if (samples[ib].second < fold) { + // Re-expand the gradient of the model around the new point: + auto const xnew = along(M.xb, p, samples[ib].first); + for (size_t k = 0; k < n; ++k) { + for (size_t l = 0; l < n; ++l) { + M.g[k] += M.G[k * n + l] * (xnew[l] - M.xb[l]); + } + } + M.xb = xnew; + M.fb = samples[ib].second; + } + return (fold - M.fb) / -predicted; + }; + std::vector free(n, true); + Real r = quadratic_step(free); + Real f_step_start = f; + std::vector xold = x; + bool last_full = true; + for (size_t step = 0; step < params_.max_local_search_steps && !done(); ++step) { + // Step 3: + bool improved = M.fb < f_step_start; + Real gx = 0; + for (size_t i = 0; i < n; ++i) { + gx += abs(M.g[i]) * (std::max)(abs(M.xb[i]), abs(xold[i])); + } + bool const criterion = !improved || gx < Real(1e-18) * (f0_ - M.fb); + std::vector inner; + std::vector at_bound; + for (size_t i = 0; i < n; ++i) { + free[i] = u_[i] < M.xb[i] && M.xb[i] < v_[i]; + (free[i] ? inner : at_bound).push_back(i); + } + if (at_bound.empty() && last_full && criterion) { + break; + } + f_step_start = M.fb; + xold = M.xb; + // Step 4: try to get away from the boundary. + if (!at_bound.empty() && criterion) { + for (size_t i : at_bound) { + std::vector e(n, Real(0)); + e[i] = 1; + std::vector> samples{{Real(0), M.fb}}; + Real const step0 = Real(0.25) * (1 + abs(M.xb[i] - x0[i])); + line_search(M.xb, e, u_[i] - M.xb[i], v_[i] - M.xb[i], samples, params_.line_search_points, step0); + for (auto const &sample : samples) { + if (sample.second < M.fb) { + M.fb = sample.second; + M.xb[i] = std::clamp(xold[i] + sample.first, u_[i], v_[i]); + } + } + } + if (!(M.fb < f_step_start)) { + break; + } + inner.clear(); + for (size_t i = 0; i < n; ++i) { + free[i] = u_[i] < M.xb[i] && M.xb[i] < v_[i]; + if (free[i]) { + inner.push_back(i); + } + } + } + // Step 5: a triple search at distance delta in the free coordinates. + last_full = abs(r - 1) > Real(0.25) || criterion; + triple_search(M, inner, last_full, [&](size_t i, model &Mi, Real &t1, Real &f1, Real &t2, Real &f2) { + Real const h = delta * (std::max)(abs(Mi.xb[i]), Real(1)); + if (Mi.xb[i] - h < u_[i]) { + t1 = Mi.xb[i] + h; + t2 = Mi.xb[i] + 2 * h; + } + else if (Mi.xb[i] + h > v_[i]) { + t1 = Mi.xb[i] - h; + t2 = Mi.xb[i] - 2 * h; + } + else { + t1 = Mi.xb[i] - h; + t2 = Mi.xb[i] + h; + } + t1 = std::clamp(t1, u_[i], v_[i]); + t2 = std::clamp(t2, u_[i], v_[i]); + if (t1 == Mi.xb[i] || t2 == Mi.xb[i] || t1 == t2) { + return triple_status::skip; + } + std::vector y = Mi.xb; + y[i] = t1; + f1 = evaluate(y); + y[i] = t2; + f2 = evaluate(y); + return done() ? triple_status::stop : triple_status::ok; + }); + if (done()) { + break; + } + // Step 6: + if (r < Real(0.25)) { + for (auto &di : d) { + di /= 2; + } + } + else if (r > Real(0.75)) { + for (auto &di : d) { + di *= 2; + } + } + r = quadratic_step(free); + } + return {M.xb, M.fb}; + } + + Real distance2(std::vector const &a, std::vector const &b) const { + Real s = 0; + for (size_t i = 0; i < n_; ++i) { + s += (a[i] - b[i]) * (a[i] - b[i]); + } + return s; + } + + std::vector between(std::vector const &x, std::vector const &w, Real t) const { + std::vector y(n_); + for (size_t i = 0; i < n_; ++i) { + y[i] = x[i] + t * (w[i] - x[i]); + } + return y; + } + + // The basket points sorted by distance to x. + std::vector basket_by_distance(std::vector const &x) const { + std::vector order(basket_.size()); + for (size_t k = 0; k < order.size(); ++k) { + order[k] = k; + } + std::stable_sort(order.begin(), order.end(), [&](size_t a, size_t b) { + return distance2(basket_[a].first, x) < distance2(basket_[b].first, x); + }); + return order; + } + + // Section 5.4: the candidates of this sweep, by ascending function value, start local searches unless they seem to + // lie in the basin of attraction of a point already in the shopping basket. + void shopping_basket() { + std::vector, Real>> candidates; + for (size_t b : candidates_) { + candidates.emplace_back(boxes_[b].x, boxes_[b].f); + } + std::stable_sort(candidates.begin(), candidates.end(), + [](auto const &a, auto const &b) { return a.second < b.second; }); + for (auto &candidate : candidates) { + auto const base_point = candidate.first; + auto &x = candidate.first; + Real &fx = candidate.second; + // Step 1: + if (std::find(started_.begin(), started_.end(), base_point) != started_.end()) { + continue; + } + bool discard = false; + for (size_t k : basket_by_distance(x)) { + auto const &w = basket_[k].first; + Real const fw = basket_[k].second; + if (done()) { + return; + } + if (fw > fx || w == x) { + continue; + } + // Step 2: + auto const x1 = between(x, w, Real(1) / 3); + Real const f1 = evaluate(x1); + if (f1 > fx) { + continue; + } + // Step 3: + auto const x2 = between(x, w, Real(2) / 3); + Real const f2 = evaluate(x2); + if (f2 > (std::max)(f1, fw)) { + if (f1 < fx) { + x = x1; + fx = f1; + } + } + else if ((std::min)(f1, f2) < fw) { + x = f1 <= f2 ? x1 : x2; + fx = (std::min)(f1, f2); + } + else { + discard = true; + break; + } + } + if (discard || done()) { + continue; + } + // Step 4: + started_.push_back(base_point); + auto result = local_search(x, fx); + bool is_new = true; + for (size_t k : basket_by_distance(result.first)) { + auto &w = basket_[k]; + if (done()) { + return; + } + bool barrier = false; + if (w.first != result.first) { + Real const top = (std::max)(result.second, w.second); + barrier = evaluate(between(result.first, w.first, Real(1) / 3)) > top || + evaluate(between(result.first, w.first, Real(2) / 3)) > top; + } + if (!barrier) { + // The same minimizer as w: keep the better of the two. + is_new = false; + if (result.second < w.second) { + w = result; + } + break; + } + } + if (is_new) { + basket_.push_back(std::move(result)); + } + } + } + + Func const &cost_function_; + multilevel_coordinate_search_parameters const ¶ms_; + ResultType target_value_; + std::atomic *cancellation_; + std::atomic *current_minimum_cost_; + std::vector> *queries_; + size_t n_; + size_t smax_; + size_t max_function_calls_; + size_t function_calls_ = 0; + std::vector u_, v_; + std::vector xbest_; + Real fbest_ = std::numeric_limits::infinity(); + // The smallest function value found in the initialization procedure: + Real f0_ = std::numeric_limits::infinity(); + std::vector> init_list_, init_f_; + std::vector l0_; + // rank_[i] = 0 for the coordinate with the highest variability (pi_i - 1 in the paper): + std::vector rank_; + std::vector boxes_; + // The non-split boxes of each level below s_max, as min-heaps on f: + std::vector> heaps_; + // Boxes that reached level s_max in this sweep: + std::vector candidates_; + std::vector> started_; + std::vector, Real>> basket_; +}; + +} // namespace detail + +BOOST_MATH_EXPORT template +ArgumentContainer multilevel_coordinate_search( + const Func cost_function, multilevel_coordinate_search_parameters const ¶ms, + std::invoke_result_t target_value = + std::numeric_limits>::quiet_NaN(), + std::atomic *cancellation = nullptr, + std::atomic> *current_minimum_cost = nullptr, + std::vector>> *queries = nullptr) { + detail::multilevel_coordinate_search_engine engine(cost_function, params, target_value, cancellation, + current_minimum_cost, queries); + engine.initialize(); + while (engine.sweep()) { + } + return engine.best_argument(); +} + +} BOOST_MATH_NAMESPACE_END // namespace boost::math::optimization +#endif diff --git a/nanobind/src/optimization.cpp b/nanobind/src/optimization.cpp index f6014c1b62..9e9e3b5791 100644 --- a/nanobind/src/optimization.cpp +++ b/nanobind/src/optimization.cpp @@ -2,8 +2,8 @@ // Distributed under the Boost Software License, Version 1.0. // https://www.boost.org/LICENSE_1_0.txt // -// The derivative-free optimizers (differential_evolution, jso, random_search, and cma_es where Eigen is found) on -// std::vector, +// The derivative-free optimizers (differential_evolution, jso, random_search, multilevel_coordinate_search, and cma_es +// where Eigen is found) on std::vector, // minimizing a Python callable f(x) -> float over the box [lower_bounds, upper_bounds]. Each returns the // best argument found. The C++ calls f from its own worker threads, so the GIL is released for the run and // taken for each call; an exception raised by f cancels the run and is re-raised when it returns. @@ -14,6 +14,7 @@ #include #include +#include #include #if __has_include() #include @@ -170,6 +171,27 @@ void bind_optimization(nb::module_& m) arg("threads") = 1u, arg("target_value") = no_target, arg("seed") = 0, "optimization::random_search: the best argument found in the box"); + m.def("multilevel_coordinate_search", [](nb::callable f, vector lower_bounds, vector upper_bounds, + std::optional initial_guess, std::size_t max_function_calls, + std::size_t max_level, bool local_search, std::size_t max_local_search_steps, + std::size_t line_search_points, double target_value) { + check_bounds(lower_bounds, upper_bounds, "multilevel_coordinate_search"); + bo::multilevel_coordinate_search_parameters params; + params.lower_bounds = std::move(lower_bounds); + params.upper_bounds = std::move(upper_bounds); + params.initial_guess = initial_guess ? &*initial_guess : nullptr; + params.max_function_calls = max_function_calls; + params.max_level = max_level; + params.local_search = local_search; + params.max_local_search_steps = max_local_search_steps; + params.line_search_points = line_search_points; + python_cost cost(f); + return run_released(cost, [&] { return bo::multilevel_coordinate_search(std::cref(cost), params, target_value, cost.cancellation()); }); + }, arg("f"), arg("lower_bounds"), arg("upper_bounds"), arg("initial_guess") = nb::none(), arg("max_function_calls") = 0, + arg("max_level") = 0, arg("local_search") = true, arg("max_local_search_steps") = 50, arg("line_search_points") = 15, + arg("target_value") = no_target, + "optimization::multilevel_coordinate_search: the argument in the box that minimizes f (0 means the defaults from the paper)"); + #ifdef BOOST_MATH_NB_HAS_CMA_ES m.def("cma_es", [](nb::callable f, vector lower_bounds, vector upper_bounds, std::optional initial_guess, std::size_t max_generations, std::size_t population_size, double learning_rate, double target_value, diff --git a/nanobind/tests/test_smoke.py b/nanobind/tests/test_smoke.py index c232eb96d1..af8f7c7e85 100644 --- a/nanobind/tests/test_smoke.py +++ b/nanobind/tests/test_smoke.py @@ -278,7 +278,7 @@ def test_array_matches_scalar(case): SIEVE = {"prime_sieve", "prime_range", "prime_count", "polynomial_roots", "polynomial_roots_aberth", "polynomial_roots_companion_matrix"} QUADRATURE = {"tanh_sinh", "exp_sinh", "sinh_sinh", "trapezoidal", "gauss", "gauss_kronrod", "gauss_hermite", "ooura_fourier_sin", "ooura_fourier_cos", "van_den_bos_unit_square"} -OPTIMIZATION = {"differential_evolution", "jso", "random_search", "cma_es"} +OPTIMIZATION = {"differential_evolution", "jso", "random_search", "multilevel_coordinate_search", "cma_es"} def test_every_function_has_a_case(): @@ -292,7 +292,7 @@ def f(x): return (x[0] - 1) ** 2 + (x[1] + 2) ** 2 lo, hi = [-5.0, -5.0], [5.0, 5.0] runs = [bm.differential_evolution(f, lo, hi, population=50, max_generations=200, threads=4), bm.jso(f, lo, hi), - bm.random_search(f, lo, hi, max_function_calls=20000)] + bm.random_search(f, lo, hi, max_function_calls=20000), bm.multilevel_coordinate_search(f, lo, hi)] if hasattr(bm, "cma_es"): # Only where Eigen was found. runs.append(bm.cma_es(f, lo, hi, max_generations=200)) for x in runs: diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index aec73a5d77..27ca2eda55 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -1153,8 +1153,10 @@ test-suite interpolators : [ run test_differential_evolution.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run test_jso.cpp : : : clang:-fno-fast-math gcc:-fno-fast-math msvc:"/fp:precise" [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run test_random_search.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] + [ run test_multilevel_coordinate_search.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run test_cma_es.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] [ compile compile_test/test_random_search_incl.cpp : [ requires cxx17_if_constexpr cxx17_std_apply ] ] + [ compile compile_test/test_multilevel_coordinate_search_incl.cpp : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ compile compile_test/test_differential_evolution_incl.cpp : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ compile compile_test/test_jso_incl.cpp : [ requires cxx17_if_constexpr cxx17_std_apply ] ] ; diff --git a/test/compile_test/test_multilevel_coordinate_search_incl.cpp b/test/compile_test/test_multilevel_coordinate_search_incl.cpp new file mode 100644 index 0000000000..9f305d4487 --- /dev/null +++ b/test/compile_test/test_multilevel_coordinate_search_incl.cpp @@ -0,0 +1,21 @@ +// Copyright Nick Thompson 2026. +// 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) +// +// Basic sanity check that header +// #includes all the files that it needs to. +// +#include +// +// Note this header includes no other headers, this is +// important if this test is to be meaningful: +// +#include "test_compile_result.hpp" + +void compile_and_link_test() +{ + auto func = [](std::vector const & v) { return v[0]*v[0]; }; + boost::math::optimization::multilevel_coordinate_search_parameters> params; + auto v = boost::math::optimization::multilevel_coordinate_search(func, params); +} diff --git a/test/test_multilevel_coordinate_search.cpp b/test/test_multilevel_coordinate_search.cpp new file mode 100644 index 0000000000..2a490dc3bc --- /dev/null +++ b/test/test_multilevel_coordinate_search.cpp @@ -0,0 +1,235 @@ +/* + * Copyright Nick Thompson, 2026 + * 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) + */ + +#ifndef BOOST_MATH_BUILD_MODULE +#include +#else +import boost.math; +#endif + +#include "math_unit_test.hpp" +#include "test_functions_for_optimization.hpp" +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#endif +using boost::math::optimization::multilevel_coordinate_search; +using boost::math::optimization::multilevel_coordinate_search_parameters; +using boost::math::optimization::validate_multilevel_coordinate_search_parameters; + + +// The bounds are asymmetric throughout, so that the midpoint of the box (the default base point) is not the minimizer. +template void test_ackley() { + std::cout << "Testing multilevel coordinate search on Ackley function . . .\n"; + using ArgType = std::array; + auto mcs_params = multilevel_coordinate_search_parameters(); + mcs_params.lower_bounds = {-5, -4}; + mcs_params.upper_bounds = {3, 5}; + auto local_minima = multilevel_coordinate_search(ackley, mcs_params); + CHECK_LE(std::abs(local_minima[0]), std::sqrt(std::numeric_limits::epsilon())); + CHECK_LE(std::abs(local_minima[1]), std::sqrt(std::numeric_limits::epsilon())); + + // Does it work with a lambda? + auto ack = [](std::array const &x) { return ackley(x); }; + local_minima = multilevel_coordinate_search(ack, mcs_params); + CHECK_LE(std::abs(local_minima[0]), std::sqrt(std::numeric_limits::epsilon())); + CHECK_LE(std::abs(local_minima[1]), std::sqrt(std::numeric_limits::epsilon())); + + // Test that if an intial guess is the exact solution, the returned solution is the exact solution: + std::array initial_guess{0, 0}; + mcs_params.initial_guess = &initial_guess; + local_minima = multilevel_coordinate_search(ack, mcs_params); + CHECK_EQUAL(local_minima[0], Real(0)); + CHECK_EQUAL(local_minima[1], Real(0)); + + std::atomic cancel = false; + Real target_value = ackley(initial_guess); + std::atomic current_minimum_cost = std::numeric_limits::quiet_NaN(); + // Test query storage: + std::vector> queries; + local_minima = multilevel_coordinate_search(ack, mcs_params, target_value, &cancel, ¤t_minimum_cost, &queries); + CHECK_EQUAL(local_minima[0], Real(0)); + CHECK_EQUAL(local_minima[1], Real(0)); + CHECK_EQUAL(current_minimum_cost.load(), target_value); + // The initial guess reaches the target, so it is the only call: + CHECK_EQUAL(queries.size(), size_t(1)); + for (auto const & q : queries) { + CHECK_EQUAL(ackley(q.first), q.second); + } +} + + +template void test_rosenbrock_saddle() { + std::cout << "Testing multilevel coordinate search on Rosenbrock saddle . . .\n"; + using ArgType = std::array; + auto mcs_params = multilevel_coordinate_search_parameters(); + mcs_params.lower_bounds = {0.5, 0.5}; + mcs_params.upper_bounds = {2.048, 2.048}; + auto local_minima = multilevel_coordinate_search(rosenbrock_saddle, mcs_params); + CHECK_ABSOLUTE_ERROR(Real(1), local_minima[0], Real(1e-7)); + CHECK_ABSOLUTE_ERROR(Real(1), local_minima[1], Real(1e-7)); + + // Does cancellation work? Cancel from inside the cost function after a few calls: + std::atomic cancel = false; + size_t calls = 0; + auto f = [&](ArgType const &x) { + if (++calls == 10) { + cancel = true; + } + return rosenbrock_saddle(x); + }; + std::vector> queries; + local_minima = multilevel_coordinate_search(f, mcs_params, std::numeric_limits::quiet_NaN(), &cancel, + nullptr, &queries); + CHECK_EQUAL(queries.size(), size_t(10)); + CHECK_GE(std::abs(local_minima[0] - Real(1)), std::sqrt(std::numeric_limits::epsilon())); + + // Cancelled before the first call: no queries, and the midpoint of the box comes back. + queries.clear(); + local_minima = multilevel_coordinate_search(f, mcs_params, std::numeric_limits::quiet_NaN(), &cancel, + nullptr, &queries); + CHECK_EQUAL(queries.size(), size_t(0)); + CHECK_EQUAL(local_minima[0], (mcs_params.lower_bounds[0] + mcs_params.upper_bounds[0]) / 2); + + // max_function_calls is respected: + mcs_params.max_function_calls = 50; + queries.clear(); + multilevel_coordinate_search(rosenbrock_saddle, mcs_params, std::numeric_limits::quiet_NaN(), nullptr, + nullptr, &queries); + CHECK_LE(queries.size(), size_t(50)); +} + + +template void test_rastrigin() { + std::cout << "Testing multilevel coordinate search on Rastrigin function (global minimum = (0,0,...,0))\n"; + using ArgType = std::vector; + auto mcs_params = multilevel_coordinate_search_parameters(); + mcs_params.lower_bounds = {-5.12, -4, -3}; + mcs_params.upper_bounds = {4, 5.12, 5}; + std::vector> queries; + auto local_minima = multilevel_coordinate_search(rastrigin, mcs_params, std::numeric_limits::quiet_NaN(), + nullptr, nullptr, &queries); + CHECK_EQUAL(local_minima.size(), size_t(3)); + CHECK_LE(rastrigin(local_minima), Real(1e-8)); + size_t const calls = queries.size(); + + // By definition, the value of the function which a target value is provided must be <= target_value. + Real target_value = 2.0; + queries.clear(); + local_minima = multilevel_coordinate_search(rastrigin, mcs_params, target_value, nullptr, nullptr, &queries); + CHECK_LE(rastrigin(local_minima), target_value); + CHECK_LE(queries.size() + 1, calls); +} + + +// Tests NaN return types and return type != input type: +void test_sphere() { + std::cout << "Testing multilevel coordinate search on sphere . . .\n"; + using ArgType = std::vector; + auto mcs_params = multilevel_coordinate_search_parameters(); + // The sphere function is NaN outside the unit ball, which includes the corners and most of the faces of this box: + mcs_params.lower_bounds.resize(4, -1); + mcs_params.upper_bounds.resize(4, 0.6f); + auto local_minima = multilevel_coordinate_search(sphere, mcs_params); + for (auto x : local_minima) { + CHECK_ABSOLUTE_ERROR(0.0f, x, 0.001f); + } + // An easy function converges without the local search, too: + mcs_params.local_search = false; + local_minima = multilevel_coordinate_search(sphere, mcs_params); + for (auto x : local_minima) { + CHECK_ABSOLUTE_ERROR(0.0f, x, 0.02f); + } +} + + +template +void test_beale() { + std::cout << "Testing multilevel coordinate search on the Beale function . . .\n"; + using ArgType = std::array; + auto mcs_params = multilevel_coordinate_search_parameters(); + mcs_params.lower_bounds = {-5, -5}; + mcs_params.upper_bounds = {5, 5}; + auto local_minima = multilevel_coordinate_search(beale, mcs_params); + CHECK_ABSOLUTE_ERROR(Real(3), local_minima[0], Real(1e-8)); + CHECK_ABSOLUTE_ERROR(Real(1)/Real(2), local_minima[1], Real(1e-8)); +} + + +// Two of the Dixon-Szego test functions of Section 7 of the paper, on the boxes used there. +template +void test_six_hump_camel() { + std::cout << "Testing multilevel coordinate search on the six-hump camel function . . .\n"; + using ArgType = std::array; + auto camel = [](ArgType const &v) { + Real x = v[0]; + Real y = v[1]; + return (4 - Real(21) / 10 * x * x + x * x * x * x / 3) * x * x + x * y + (-4 + 4 * y * y) * y * y; + }; + auto mcs_params = multilevel_coordinate_search_parameters(); + mcs_params.lower_bounds = {-3, -2}; + mcs_params.upper_bounds = {3, 2}; + auto local_minima = multilevel_coordinate_search(camel, mcs_params); + CHECK_ABSOLUTE_ERROR(Real(-1.0316284534898774), camel(local_minima), Real(1e-12)); +} + + +template +void test_branin() { + std::cout << "Testing multilevel coordinate search on the Branin function . . .\n"; + using ArgType = std::array; + auto branin = [](ArgType const &v) { + using boost::math::constants::pi; + using std::cos; + Real x = v[0]; + Real y = v[1]; + Real t = y - Real(5.1) / (4 * pi() * pi()) * x * x + 5 * x / pi() - 6; + return t * t + 10 * (1 - 1 / (8 * pi())) * cos(x) + 10; + }; + auto mcs_params = multilevel_coordinate_search_parameters(); + mcs_params.lower_bounds = {-5, 0}; + mcs_params.upper_bounds = {10, 15}; + auto local_minima = multilevel_coordinate_search(branin, mcs_params); + CHECK_ABSOLUTE_ERROR(Real(0.39788735772973816), branin(local_minima), Real(1e-12)); +} + + +void test_validation() { + std::cout << "Testing multilevel coordinate search parameter validation . . .\n"; + using ArgType = std::array; + auto mcs_params = multilevel_coordinate_search_parameters(); + mcs_params.lower_bounds = {-1, -1}; + mcs_params.upper_bounds = {1, 1}; + validate_multilevel_coordinate_search_parameters(mcs_params); + + auto bad_params = mcs_params; + bad_params.upper_bounds[1] = -1; + CHECK_THROW(validate_multilevel_coordinate_search_parameters(bad_params), std::domain_error); + auto run = [&]() { return multilevel_coordinate_search(rosenbrock_saddle, bad_params); }; + CHECK_THROW(run(), std::domain_error); + + bad_params = mcs_params; + bad_params.line_search_points = 2; + CHECK_THROW(validate_multilevel_coordinate_search_parameters(bad_params), std::invalid_argument); + + bad_params = mcs_params; + bad_params.max_level = 1; + CHECK_THROW(validate_multilevel_coordinate_search_parameters(bad_params), std::invalid_argument); +} + +int main() { + test_ackley(); + test_ackley(); + test_rosenbrock_saddle(); + test_rastrigin(); + test_sphere(); + test_beale(); + test_six_hump_camel(); + test_branin(); + test_validation(); + return boost::math::test::report_errors(); +}