From d5c7e2dc7c6947b3debbe0144f1fc8500e994098 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Thu, 22 Oct 2026 22:54:27 -0700 Subject: [PATCH] Add TuRBO, trust-region Bayesian optimization (Eriksson et al. 2019) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit TuRBO-1 with a Matérn-5/2 Gaussian process, Thompson sampling and hand-rolled linear algebra; no dependencies. Co-Authored-By: Claude Opus 5.5 --- doc/math.qbk | 1 + doc/optimization/turbo.qbk | 195 +++++ include/boost/math/optimization/turbo.hpp | 931 ++++++++++++++++++++++ module/CMakeLists.txt | 1 + module/math.cppm | 1 + test/Jamfile.v2 | 2 + test/compile_test/test_turbo_incl.cpp | 22 + test/test_turbo.cpp | 290 +++++++ 8 files changed, 1443 insertions(+) create mode 100644 doc/optimization/turbo.qbk create mode 100644 include/boost/math/optimization/turbo.hpp create mode 100644 test/compile_test/test_turbo_incl.cpp create mode 100644 test/test_turbo.cpp diff --git a/doc/math.qbk b/doc/math.qbk index fd79cdd91b..5539bfb53a 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -740,6 +740,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include optimization/jso.qbk] [include optimization/random_search.qbk] [include optimization/multilevel_coordinate_search.qbk] +[include optimization/turbo.qbk] [include optimization/cma_es.qbk] [include optimization/gradient_optimizers.qbk] [endmathpart] [/mathpart optimization Optimization] diff --git a/doc/optimization/turbo.qbk b/doc/optimization/turbo.qbk new file mode 100644 index 0000000000..a9a10899e2 --- /dev/null +++ b/doc/optimization/turbo.qbk @@ -0,0 +1,195 @@ +[/ +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:turbo Trust Region Bayesian Optimization] + +[heading Synopsis] + +`` + #include + + namespace boost::math::optimization { + + template struct turbo_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 initial_points = 0; + size_t batch_size = 1; + size_t candidates = 0; + size_t max_training_points = 200; + size_t hyperparameter_steps = 50; + }; + + template + void validate_turbo_parameters(turbo_parameters const ¶ms); + + template + ArgumentContainer turbo(const Func cost_function, + turbo_parameters const ¶ms, + URBG &gen, + 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 `turbo` function searches for a global minimum of a function on a box with TuRBO-1, the trust region Bayesian optimization algorithm of Eriksson, Pearce, Gardner, Turner and Poloczek. +It is meant for expensive functions: it spends milliseconds of its own time per function call, so that it can spend as few calls as possible. + +TuRBO fits a Gaussian process (an ARD Matérn-5/2 kernel with a constant mean and a noise term) to the points it has evaluated, +but uses it only in a /trust region/: a box centred at the best point found so far, stretched along each coordinate in proportion to the lengthscale the Gaussian process learned for it. +Each new point is the minimizer, over a few hundred candidate points in the trust region, of a random function drawn from the posterior (Thompson sampling). +The trust region doubles in size after three batches in a row that improve on its best point, and halves after ⌈max(4, /n/)\/ /q/⌉ batches in a row that do not, where /q/ is the batch size; +when it becomes too small, TuRBO forgets it and starts again from a new random design, so that the search as a whole is global. + +TuRBO is a good choice when each function call takes much longer than a few milliseconds and the budget is a few hundred to a few thousand calls. +It finds good values quickly, but it does not converge to high accuracy: +the trust region never gets smaller than 2[super -7] of the box, and the Gaussian process treats differences in function values smaller than its noise level as noise. +For cheap functions, or when the minimum is needed to many digits, multilevel coordinate search, differential evolution, jSO or CMA-ES are better choices, +or the result of TuRBO can be polished by a local method. + +[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. It is the first function call, in place of the first point of the first random design. +* `max_function_calls`: Defaults to 100 times the dimension. +* `initial_points`: The number of points of the random (Latin hypercube) design that starts each trust region. Defaults to twice the dimension. +* `batch_size`: The number of points /q/ proposed at a time, from independent Thompson samples of the same Gaussian process. Defaults to 1. + The points of a batch could be evaluated in parallel, but `turbo` itself evaluates them one after the other. +* `candidates`: The number of candidate points in the trust region at which each Thompson sample is drawn. + Defaults to min(100/n/, 1000). It must be at least `batch_size`. +* `max_training_points`: The Gaussian process is fitted to at most this many points of the current trust region, those nearest its centre. Defaults to 200. +* `hyperparameter_steps`: The number of steps of the Adam optimizer on the log marginal likelihood each time the Gaussian process is fitted. Defaults to 50. + +`validate_turbo_parameters` throws a `std::domain_error` if the bounds or initial guess are invalid, +and a `std::invalid_argument` if `batch_size`, `candidates` or `max_training_points` is out of range. +`turbo` calls it before the first function evaluation. + +[heading The function] + +`` +template +ArgumentContainer turbo(const Func cost_function, + turbo_parameters const ¶ms, + URBG &gen, + 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. +* `gen`: A uniform random bit generator, like `std::mt19937_64`. +* `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. +If the computation is cancelled before the first function call, it is the midpoint of the box (or the initial guess). + +[heading Termination] + +As in the paper, there is no convergence criterion: TuRBO restarts until it is stopped. +The computation stops when the `target_value` is reached, when `max_function_calls` calls have been made, or when `cancellation` is set. +These are checked before every call, even in the middle of a batch, so `max_function_calls` is never exceeded. + +[heading Cost] + +Each iteration fits the Gaussian process, with `hyperparameter_steps` Cholesky factorizations and inverses of an /N/ x /N/ matrix, +where /N/ is the number of training points, and draws the Thompson samples, with a Cholesky factorization of the /M/ x /M/ posterior covariance of the /M/ candidates. +That is O(`hyperparameter_steps` /N/[super 3] + /M/[super 3] + /M/[super 2]/N/) operations per iteration. +A trust region lasts about 25/n/ function calls, so /N/ reaches the default cap of 200 only for /n/ of about 8 or more. +On an Apple M1 Pro core, an iteration with the default parameters takes about 1 ms for /n/ = 2, 4 ms for /n/ = 4, 8 ms for /n/ = 5, 12 ms for /n/ = 6, and 85 ms for /n/ = 10. +Setting `max_training_points` to 100 and `candidates` to 500 brings the last down to about 15 ms. +All of the linear algebra is hand-rolled, and done in double precision whatever the `Real` type of the arguments. + +[heading Implementation notes] + +The trust region, its update rules, the Thompson sampling, the Gaussian process and the restarts follow Section 2 of the paper and its supplement, +and the parameters follow the authors' code ([@https://github.com/uber-research/TuRBO github.com/uber-research/TuRBO], `turbo_1.py` and `gp.py`) where the two differ. +Each Gaussian process fit starts from the same hyperparameters (lengthscales 0.5, output scale 1, noise variance 0.005, mean 0), +and runs 50 steps of Adam with learning rate 0.1 on the negative log marginal likelihood, with analytic gradients, +over hyperparameters mapped to the real line by scaled sigmoids onto the intervals \[0.005, 2\] (lengthscales), \[0.05, 20\] (output scale) and \[0.0005, 0.2\] (noise variance), as GPyTorch does. +The inputs are scaled to the unit cube, and the function values standardized by their median and standard deviation. +The departures from the paper and the authors' code are: + +* The failure tolerance is ⌈max(4, /n/)\/ /q/⌉, the noise variance bound 0.2, and a success an improvement by more than 10[super -3] times the magnitude of the best value in the trust region, + as in the authors' code. The paper has ⌈/n/\/ /q/⌉, 0.1, and any improvement. +* The candidates are a random Latin hypercube in the trust region rather than a scrambled Sobol sequence, with each coordinate of each candidate kept at the centre of the trust region except with probability min(20\/n, 1), as in the paper. + A candidate with no perturbed coordinate gets one, chosen from all /n/ coordinates (the authors' code never chooses the last). +* The default number of candidates is min(100/n/, 1000), not min(100/n/, 5000). + The joint posterior sample is exact, by a Cholesky factorization of the candidates' covariance, which costs /M/[super 3]/3 operations; + the authors' code also uses a Cholesky factorization up to 2000 candidates, and an approximate (Lanczos) factorization above that. +* The Gaussian process is fitted to at most `max_training_points` points of the trust region, those nearest its centre; the authors' code uses all of them. + The cap rarely matters below 8 dimensions, where trust regions do not last that long, and it bounds the cost of a fit, which grows as the cube of the number of points. +* Values that are +[infin] (or NaN) are replaced, in the training set, by the largest finite value in it. If there is none, the batch is taken at random from the candidates. +* The initial guess, if any, replaces the first point of the first design. +* If the run stops in the middle of a batch, the trust region is not updated. + +Only TuRBO-1, with a single trust region, is implemented; TuRBO-/m/ maintains /m/ trust regions and allocates each batch among them. +The paper's choice of centre for noisy functions (the point with the smallest posterior mean) is not implemented: the centre is always the best point evaluated, as in the authors' code. + +[h4:examples Examples] + +The Branin function has three global minima with value 0.397887 on the box \[-5, 10\]×\[0, 15\]: + +`` +#include +#include +#include +#include +#include + +int main() { + using boost::math::optimization::turbo; + using boost::math::optimization::turbo_parameters; + using ArgType = std::array; + auto branin = [](ArgType const & v) { + double const pi = 3.14159265358979323846; + double x = v[0]; + double y = v[1]; + double t = y - 5.1/(4*pi*pi)*x*x + 5*x/pi - 6; + return t*t + 10*(1 - 1/(8*pi))*std::cos(x) + 10; + }; + turbo_parameters params; + params.lower_bounds = {-5, 0}; + params.upper_bounds = {10, 15}; + params.max_function_calls = 200; + std::mt19937_64 gen(42); + auto x = turbo(branin, params, gen); + std::cout << "f(" << x[0] << ", " << x[1] << ") = " << branin(x) << "\n"; +} +`` + +After 200 calls the value is typically within a few times 10[super -4] of the minimum: over 30 seeds, the median distance was 3.7×10[super -4] and the largest 1.7×10[super -3]. + +[h4:references References] + +* D. Eriksson, M. Pearce, J. Gardner, R. D. Turner and M. Poloczek, ['Scalable global optimization via local Bayesian optimization.] Advances in Neural Information Processing Systems 32 (NeurIPS 2019). [@https://arxiv.org/abs/1910.01739 arXiv:1910.01739]. + +[endsect] [/section:turbo] diff --git a/include/boost/math/optimization/turbo.hpp b/include/boost/math/optimization/turbo.hpp new file mode 100644 index 0000000000..727cf9dba6 --- /dev/null +++ b/include/boost/math/optimization/turbo.hpp @@ -0,0 +1,931 @@ +/* + * 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_TURBO_HPP +#define BOOST_MATH_OPTIMIZATION_TURBO_HPP + +#include +#ifndef BOOST_MATH_BUILD_MODULE +#include +#endif +#include +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +BOOST_MATH_NAMESPACE_BEGIN namespace optimization { + +// TuRBO-1 of: Eriksson, D., Pearce, M., Gardner, J., Turner, R. D., Poloczek, M. (2019). +// Scalable global optimization via local Bayesian optimization. NeurIPS 32. arXiv:1910.01739. +// The defaults follow the authors' code (github.com/uber-research/TuRBO, turbo_1.py and gp.py), "the reference". +// Section numbers in the comments below refer to the paper and its supplement. + +BOOST_MATH_EXPORT template struct turbo_parameters { + using Real = typename ArgumentContainer::value_type; + ArgumentContainer lower_bounds; + ArgumentContainer upper_bounds; + // Replaces the first point of the first initial design. + ArgumentContainer const *initial_guess = nullptr; + // 0 means 100 times the dimension. + size_t max_function_calls = 0; + // The points of the Latin hypercube design that starts each trust region; 0 means twice the dimension. + size_t initial_points = 0; + // q, the number of points proposed (by independent Thompson samples) per Gaussian process fit. + size_t batch_size = 1; + // The candidate points of each Thompson sample; 0 means min(100n, 1000). The paper (Section D) has min(100n, 5000), + // but the joint sample costs a Cholesky factorization of the candidates' covariance, O(candidates^3). + size_t candidates = 0; + // The Gaussian process is fitted to at most this many of the trust region's points, those nearest its centre. + size_t max_training_points = 200; + // Adam steps on the log marginal likelihood per fit, as in the reference. + size_t hyperparameter_steps = 50; +}; + +BOOST_MATH_EXPORT template +void validate_turbo_parameters(turbo_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 << ": TuRBO 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.batch_size == 0) { + oss << __FILE__ << ":" << __LINE__ << ":" << __func__; + oss << ": The batch size must be at least 1."; + throw std::invalid_argument(oss.str()); + } + if (params.candidates != 0 && params.candidates < params.batch_size) { + oss << __FILE__ << ":" << __LINE__ << ":" << __func__; + oss << ": There must be at least as many candidates as the batch size, but candidates = " << params.candidates + << " and batch_size = " << params.batch_size << "."; + throw std::invalid_argument(oss.str()); + } + if (params.max_training_points == 0) { + oss << __FILE__ << ":" << __LINE__ << ":" << __func__; + oss << ": The Gaussian process needs at least one training point."; + throw std::invalid_argument(oss.str()); + } +} + +namespace detail { + +// Hand-rolled dense linear algebra on row-major double arrays. The kernels with long inner loops are written so that +// the compiler can vectorize them without reassociating floating point sums: either as axpy updates of contiguous +// rows, or as several rows' dot products with one shared operand, which gives independent accumulators. + +// s[r] = sum_{k < m} a[r][k] b[k] for four rows a[0..3]. +inline void turbo_dot4(double const *const *a, double const *b, size_t m, double *s) { + double s0 = 0, s1 = 0, s2 = 0, s3 = 0, t0 = 0, t1 = 0, t2 = 0, t3 = 0; + size_t k = 0; + for (; k + 1 < m; k += 2) { + double const x = b[k]; + double const y = b[k + 1]; + s0 += a[0][k] * x; + t0 += a[0][k + 1] * y; + s1 += a[1][k] * x; + t1 += a[1][k + 1] * y; + s2 += a[2][k] * x; + t2 += a[2][k + 1] * y; + s3 += a[3][k] * x; + t3 += a[3][k + 1] * y; + } + if (k < m) { + double const x = b[k]; + s0 += a[0][k] * x; + s1 += a[1][k] * x; + s2 += a[2][k] * x; + s3 += a[3][k] * x; + } + s[0] = s0 + t0; + s[1] = s1 + t1; + s[2] = s2 + t2; + s[3] = s3 + t3; +} + +inline double turbo_dot(double const *a, double const *b, size_t m) { + double s0 = 0, s1 = 0, s2 = 0, s3 = 0; + size_t k = 0; + for (; k + 3 < m; k += 4) { + s0 += a[k] * b[k]; + s1 += a[k + 1] * b[k + 1]; + s2 += a[k + 2] * b[k + 2]; + s3 += a[k + 3] * b[k + 3]; + } + for (; k < m; ++k) { + s0 += a[k] * b[k]; + } + return (s0 + s1) + (s2 + s3); +} + +// C[i][j] -= sum_{k < len} X[i][k] X[j][k] for lo <= j <= i < hi, where row i of X starts at X + i ldx and C has +// leading dimension ldc. Tiles of 64 rows j stay in cache while every block of four rows i below them goes by. +inline void turbo_subtract_gram(double *C, size_t ldc, double const *X, size_t ldx, size_t len, size_t lo, size_t hi) { + constexpr size_t tile = 64; + double s[4]; + double const *rows[4]; + for (size_t j0 = lo; j0 < hi; j0 += tile) { + size_t const j1 = (std::min)(hi, j0 + tile); + for (size_t i = j0; i < hi; i += 4) { + size_t const nb = (std::min)(size_t(4), hi - i); + for (size_t b = 0; b < 4; ++b) { + rows[b] = X + (i + (b < nb ? b : 0)) * ldx; + } + size_t const jend = (std::min)(j1, i + nb); + for (size_t j = j0; j < jend; ++j) { + double const *xj = X + j * ldx; + if (j < i) { + turbo_dot4(rows, xj, len, s); + for (size_t b = 0; b < nb; ++b) { + C[(i + b) * ldc + j] -= s[b]; + } + } + else { + for (size_t b = j - i; b < nb; ++b) { + C[(i + b) * ldc + j] -= turbo_dot(rows[b], xj, len); + } + } + } + } + } +} + +// The Cholesky factorization A = L L^T in place: on input the lower triangle of the n x n matrix A (leading dimension +// n), on output that of L. The upper triangle is not referenced. Returns false if A is not numerically positive +// definite. Right-looking by blocks of 64 columns; within a block, left-looking four rows at a time. +inline bool turbo_cholesky(double *A, size_t n) { + using std::sqrt; + constexpr size_t block = 64; + double s[4]; + double *rows[4]; + double const *parts[4]; + for (size_t k0 = 0; k0 < n; k0 += block) { + size_t const k1 = (std::min)(n, k0 + block); + // Columns k0 to k1 of L, from the diagonal block down: + for (size_t i = k0; i < n; i += 4) { + size_t const nb = (std::min)(size_t(4), n - i); + for (size_t b = 0; b < 4; ++b) { + rows[b] = A + (i + (b < nb ? b : 0)) * n; + parts[b] = rows[b] + k0; + } + size_t const jend = (std::min)(k1, i + nb); + for (size_t j = k0; j < jend; ++j) { + double const *lj = A + j * n; + if (j < i) { + turbo_dot4(parts, lj + k0, j - k0, s); + for (size_t b = 0; b < nb; ++b) { + rows[b][j] = (rows[b][j] - s[b]) / lj[j]; + } + continue; + } + // Column j meets these rows on the diagonal: row j first, then the rows below it. + for (size_t b = j - i; b < nb; ++b) { + double const d = rows[b][j] - turbo_dot(parts[b], lj + k0, j - k0); + if (b == j - i) { + if (!(d > 0)) { + return false; + } + rows[b][j] = sqrt(d); + } + else { + rows[b][j] = d / lj[j]; + } + } + } + } + // The trailing submatrix loses the contribution of these columns: + turbo_subtract_gram(A, n, A + k0, n, k1 - k0, k1, n); + } + return true; +} + +// Solves L X = B in place for the n x m matrix B (leading dimension m), where L is the lower triangle of an n x n +// matrix (leading dimension n). +inline void turbo_forward_solve(double const *L, size_t n, double *B, size_t m) { + for (size_t i = 0; i < n; ++i) { + double *bi = B + i * m; + for (size_t k = 0; k < i; ++k) { + double const lik = L[i * n + k]; + double const *bk = B + k * m; + for (size_t c = 0; c < m; ++c) { + bi[c] -= lik * bk[c]; + } + } + double const inv = 1 / L[i * n + i]; + for (size_t c = 0; c < m; ++c) { + bi[c] *= inv; + } + } +} + +// Solves L L^T x = b in place. +inline void turbo_cholesky_solve(double const *L, size_t n, double *b) { + for (size_t i = 0; i < n; ++i) { + b[i] = (b[i] - turbo_dot(L + i * n, b, i)) / L[i * n + i]; + } + // L^T x = b, by rows of L: once x_i is known, remove its contribution L_ik x_i from b_k for k < i. + for (size_t i = n; i-- > 0;) { + double const *li = L + i * n; + double const xi = b[i] / li[i]; + b[i] = xi; + for (size_t k = 0; k < i; ++k) { + b[k] -= li[k] * xi; + } + } +} + +// The lower triangle of (L L^T)^{-1} = Z^T Z, where Z = L^{-1}. Z is n x n scratch. +inline void turbo_cholesky_inverse(double const *L, size_t n, double *Z, double *Kinv) { + for (size_t i = 0; i < n; ++i) { + double *zi = Z + i * n; + for (size_t j = 0; j < i; ++j) { + zi[j] = 0; + } + for (size_t k = 0; k < i; ++k) { + double const lik = L[i * n + k]; + double const *zk = Z + k * n; + for (size_t j = 0; j <= k; ++j) { + zi[j] -= lik * zk[j]; + } + } + double const inv = 1 / L[i * n + i]; + for (size_t j = 0; j < i; ++j) { + zi[j] *= inv; + } + zi[i] = inv; + } + for (size_t a = 0; a < n; ++a) { + for (size_t b = 0; b <= a; ++b) { + Kinv[a * n + b] = 0; + } + } + for (size_t k = 0; k < n; ++k) { + double const *zk = Z + k * n; + for (size_t a = 0; a <= k; ++a) { + double const zka = zk[a]; + double *row = Kinv + a * n; + for (size_t b = 0; b <= a; ++b) { + row[b] += zka * zk[b]; + } + } + } +} + +// The ARD Matern-5/2 kernel (Section C) on the unit cube: s (1 + sqrt(5) r + 5 r^2/3) exp(-sqrt(5) r), where +// r^2 = sum_i (x_i - y_i)^2 / l_i^2. inv_l2[i] = 1/l_i^2. +inline double turbo_matern52(double const *x, double const *y, double const *inv_l2, size_t d) { + using std::exp; + using std::sqrt; + double r2 = 0; + for (size_t i = 0; i < d; ++i) { + double const t = x[i] - y[i]; + r2 += t * t * inv_l2[i]; + } + double const r = sqrt(5 * r2); + return (1 + r + r * r / 3) * exp(-r); +} + +// The state of a TuRBO-1 run; turbo() below drives it. +template class turbo_engine { +public: + using Real = typename ArgumentContainer::value_type; + using ResultType = std::invoke_result_t; + + turbo_engine(Func const &cost_function, turbo_parameters const ¶ms, URBG &gen, + ResultType target_value, std::atomic *cancellation, std::atomic *current_minimum_cost, + std::vector> *queries) + : cost_function_(cost_function), params_(params), gen_(gen), target_value_(target_value), + cancellation_(cancellation), current_minimum_cost_(current_minimum_cost), queries_(queries) { + validate_turbo_parameters(params); + n_ = params.lower_bounds.size(); + q_ = params.batch_size; + max_function_calls_ = params.max_function_calls ? params.max_function_calls : 100 * n_; + initial_points_ = params.initial_points ? params.initial_points : 2 * n_; + candidates_ = params.candidates ? params.candidates : (std::min)(100 * n_, size_t(1000)); + candidates_ = (std::max)(candidates_, q_); + // tau_fail = ceil(max(4/q, n/q)), as in the reference; the paper (Section D) has ceil(n/q). + failure_tolerance_ = ((std::max)(size_t(4), n_) + q_ - 1) / q_; + if constexpr (detail::has_resize_v) { + best_argument_.resize(n_); + } + for (size_t i = 0; i < n_; ++i) { + best_argument_[i] = params.initial_guess ? (*params.initial_guess)[i] + : params.lower_bounds[i] + (params.upper_bounds[i] - params.lower_bounds[i]) / 2; + } + } + + ArgumentContainer optimize() { + bool first = true; + while (!done()) { + restart(first); + first = false; + while (!done() && length_ >= length_min_) { + step(); + } + } + return best_argument_; + } + +private: + // Section D and the reference: + static constexpr double length_init_ = 0.8; + static constexpr double length_min_ = 0.0078125; // 2^-7 + static constexpr double length_max_ = 1.6; + static constexpr size_t success_tolerance_ = 3; + + // The hyperparameter bounds of the reference (gp.py); the paper (Section C) has 0.1 for the upper noise bound. + static constexpr double lengthscale_bounds_[2] = {0.005, 2.0}; + static constexpr double outputscale_bounds_[2] = {0.05, 20.0}; + static constexpr double noise_bounds_[2] = {0.0005, 0.2}; + + bool done() const { + using std::isnan; + if (cancellation_ && *cancellation_) { + return true; + } + if (calls_ > 0 && !isnan(target_value_) && best_value_ <= target_value_) { + return true; + } + if (calls_ > 0 && best_value_ == -std::numeric_limits::infinity()) { + return true; + } + return calls_ >= max_function_calls_; + } + + ArgumentContainer to_argument(double const *u) const { + ArgumentContainer x; + if constexpr (detail::has_resize_v) { + x.resize(n_); + } + for (size_t i = 0; i < n_; ++i) { + Real const lo = params_.lower_bounds[i]; + Real const hi = params_.upper_bounds[i]; + Real const xi = lo + static_cast(u[i]) * (hi - lo); + x[i] = xi < lo ? lo : (xi > hi ? hi : xi); + } + return x; + } + + // Calls the cost function at x (with unit-cube coordinates u), and appends the point to the trust region's data. + // NaN counts as +infinity. + void evaluate(ArgumentContainer const &x, double const *u) { + using std::isnan; + ResultType const value = cost_function_(x); + ++calls_; + if (queries_) { + queries_->push_back(std::make_pair(x, value)); + } + ResultType const f = isnan(value) ? std::numeric_limits::infinity() : value; + if (calls_ == 1 || f < best_value_) { + best_value_ = f; + best_argument_ = x; + if (current_minimum_cost_) { + *current_minimum_cost_ = f; + } + } + X_.insert(X_.end(), u, u + n_); + y_.push_back(static_cast(f)); + } + + double uniform01() { return std::uniform_real_distribution(0, 1)(gen_); } + + // A random Latin hypercube design of m points in the box [lo, hi] (row-major m x n), as latin_hypercube in the + // reference: one point in each of m equal slices of every coordinate, uniformly distributed within its slice. + std::vector latin_hypercube(size_t m, double const *lo, double const *hi) { + std::vector design(m * n_); + std::vector slices(m); + for (size_t i = 0; i < n_; ++i) { + for (size_t k = 0; k < m; ++k) { + slices[k] = k; + } + std::shuffle(slices.begin(), slices.end(), gen_); + for (size_t k = 0; k < m; ++k) { + double const t = (static_cast(slices[k]) + uniform01()) / static_cast(m); + design[k * n_ + i] = (std::min)(lo[i] + t * (hi[i] - lo[i]), hi[i]); + } + } + return design; + } + + // Starts a new trust region from a new initial design. + void restart(bool first) { + length_ = length_init_; + successes_ = 0; + failures_ = 0; + X_.clear(); + y_.clear(); + std::vector lo(n_, 0.0), hi(n_, 1.0); + std::vector design = latin_hypercube(initial_points_, lo.data(), hi.data()); + for (size_t k = 0; k < initial_points_ && !done(); ++k) { + double *u = design.data() + k * n_; + if (first && k == 0 && params_.initial_guess) { + auto const &g = *params_.initial_guess; + for (size_t i = 0; i < n_; ++i) { + Real const lo_i = params_.lower_bounds[i]; + Real const hi_i = params_.upper_bounds[i]; + u[i] = static_cast((g[i] - lo_i) / (hi_i - lo_i)); + } + evaluate(g, u); + } + else { + evaluate(to_argument(u), u); + } + } + } + + // One iteration of TuRBO-1: fit the Gaussian process, Thompson sample a batch in the trust region, evaluate it, and + // resize the trust region. + void step() { + using std::abs; + using std::isfinite; + size_t const N_all = y_.size(); + size_t center = 0; + for (size_t k = 1; k < N_all; ++k) { + if (y_[k] < y_[center]) { + center = k; + } + } + double const *xc = X_.data() + center * n_; + select_training_set(xc); + bool const model = fit(); + + // The trust region (Section 2): side lengths L l_i / (prod_j l_j)^(1/n), centred at the best point, clipped to the + // unit cube. The reference first divides by the mean lengthscale, for stability. + std::vector weights(n_, 1.0); + if (model) { + double mean = 0; + for (size_t i = 0; i < n_; ++i) { + mean += lengthscale_[i]; + } + mean /= static_cast(n_); + double log_prod = 0; + for (size_t i = 0; i < n_; ++i) { + weights[i] = lengthscale_[i] / mean; + log_prod += std::log(weights[i]); + } + double const scale = std::exp(-log_prod / static_cast(n_)); + for (auto &w : weights) { + w *= scale; + } + } + std::vector lo(n_), hi(n_); + for (size_t i = 0; i < n_; ++i) { + lo[i] = (std::max)(xc[i] - weights[i] * length_ / 2, 0.0); + hi[i] = (std::min)(xc[i] + weights[i] * length_ / 2, 1.0); + } + + // The candidates (Section D): a Latin hypercube in the trust region, then each coordinate is kept at the centre + // unless it is perturbed, with probability min(20/n, 1); a candidate with no perturbed coordinate gets one. + size_t const M = candidates_; + std::vector C = latin_hypercube(M, lo.data(), hi.data()); + double const p_perturb = (std::min)(20.0 / static_cast(n_), 1.0); + if (p_perturb < 1) { + std::uniform_int_distribution pick(0, n_ - 1); + for (size_t k = 0; k < M; ++k) { + double *c = C.data() + k * n_; + bool any = false; + for (size_t i = 0; i < n_; ++i) { + if (uniform01() <= p_perturb) { + any = true; + } + else { + c[i] = xc[i]; + } + } + if (!any) { + size_t const i = pick(gen_); + c[i] = lo[i] + uniform01() * (hi[i] - lo[i]); + } + } + } + + // Thompson sampling: q joint samples of the posterior (with the noise, as the reference samples the likelihood) + // over the candidates, and the argmin of each, never the same candidate twice. + std::vector chosen; + if (model) { + thompson_sample(C, M, chosen); + } + else { + // No finite function value yet: no model, so take random candidates. + std::vector order(M); + for (size_t k = 0; k < M; ++k) { + order[k] = k; + } + std::shuffle(order.begin(), order.end(), gen_); + chosen.assign(order.begin(), order.begin() + static_cast(q_)); + } + + double const local_best = y_[center]; + double batch_best = std::numeric_limits::infinity(); + for (size_t k : chosen) { + if (done()) { + return; + } + double const *u = C.data() + k * n_; + evaluate(to_argument(u), u); + batch_best = (std::min)(batch_best, y_.back()); + } + // A success is an improvement on the trust region's best value by more than 1e-3 of its magnitude (reference). + bool const success = isfinite(local_best) ? batch_best < local_best - 1e-3 * abs(local_best) : batch_best < local_best; + if (success) { + ++successes_; + failures_ = 0; + } + else { + successes_ = 0; + ++failures_; + } + if (successes_ == success_tolerance_) { + length_ = (std::min)(2 * length_, length_max_); + successes_ = 0; + } + else if (failures_ == failure_tolerance_) { + length_ /= 2; + failures_ = 0; + } + } + + // The training set: the max_training_points points of the trust region nearest its centre xc (all of them if + // there are fewer), in the unit cube, with the function values standardized as in the reference (median and + // standard deviation). Non-finite values are replaced by the largest finite one. + void select_training_set(double const *xc) { + size_t const N_all = y_.size(); + std::vector idx(N_all); + for (size_t k = 0; k < N_all; ++k) { + idx[k] = k; + } + size_t const N = (std::min)(N_all, params_.max_training_points); + if (N < N_all) { + std::vector dist(N_all); + for (size_t k = 0; k < N_all; ++k) { + double s = 0; + for (size_t i = 0; i < n_; ++i) { + double const t = X_[k * n_ + i] - xc[i]; + s += t * t; + } + dist[k] = s; + } + std::nth_element(idx.begin(), idx.begin() + static_cast(N), idx.end(), + [&](size_t a, size_t b) { return dist[a] < dist[b] || (dist[a] == dist[b] && a < b); }); + idx.resize(N); + std::sort(idx.begin(), idx.end()); + } + N_ = N; + Xt_.resize(N * n_); + yt_.resize(N); + using std::isfinite; + double worst = -std::numeric_limits::infinity(); + for (size_t k = 0; k < N; ++k) { + std::copy(X_.begin() + static_cast(idx[k] * n_), + X_.begin() + static_cast((idx[k] + 1) * n_), Xt_.begin() + static_cast(k * n_)); + yt_[k] = y_[idx[k]]; + if (isfinite(yt_[k])) { + worst = (std::max)(worst, yt_[k]); + } + } + finite_ = isfinite(worst); + if (!finite_) { + return; + } + for (auto &v : yt_) { + if (!isfinite(v)) { + v = worst; + } + } + std::vector sorted = yt_; + std::nth_element(sorted.begin(), sorted.begin() + static_cast(N / 2), sorted.end()); + double median = sorted[N / 2]; + if (N % 2 == 0) { + median = (median + *std::max_element(sorted.begin(), sorted.begin() + static_cast(N / 2))) / 2; + } + double mean = 0; + for (double v : yt_) { + mean += v; + } + mean /= static_cast(N); + double var = 0; + for (double v : yt_) { + var += (v - mean) * (v - mean); + } + using std::sqrt; + double sigma = sqrt(var / static_cast(N)); + if (!(sigma >= 1e-6)) { + sigma = 1; + } + for (auto &v : yt_) { + v = (v - median) / sigma; + } + } + + static double sigmoid(double t) { return 1 / (1 + std::exp(-t)); } + static double logit(double p) { return std::log(p / (1 - p)); } + static double from_raw(double raw, double const *bounds) { return bounds[0] + (bounds[1] - bounds[0]) * sigmoid(raw); } + static double to_raw(double value, double const *bounds) { return logit((value - bounds[0]) / (bounds[1] - bounds[0])); } + + // Sets the hyperparameters from theta = (raw lengthscales, raw outputscale, raw noise, mean), builds the kernel + // matrix K + noise I in L_, and factors it. pair_m_ and pair_dm_ receive the kernel's radial factor and the factor + // of its lengthscale derivatives at each pair of training points. Returns false if the factorization fails. + bool factor(std::vector const &theta) { + using std::exp; + using std::sqrt; + size_t const N = N_; + for (size_t i = 0; i < n_; ++i) { + lengthscale_[i] = from_raw(theta[i], lengthscale_bounds_); + inv_l2_[i] = 1 / (lengthscale_[i] * lengthscale_[i]); + } + outputscale_ = from_raw(theta[n_], outputscale_bounds_); + noise_ = from_raw(theta[n_ + 1], noise_bounds_); + mean_ = theta[n_ + 2]; + L_.resize(N * N); + pair_m_.resize(N * (N - 1) / 2); + pair_dm_.resize(pair_m_.size()); + double const root5 = sqrt(5.0); + size_t p = 0; + for (size_t a = 0; a < N; ++a) { + double *row = L_.data() + a * N; + for (size_t b = 0; b < a; ++b, ++p) { + double const *D = sqdiff_.data() + p * n_; + double r2 = 0; + for (size_t i = 0; i < n_; ++i) { + r2 += D[i] * inv_l2_[i]; + } + double const r = root5 * sqrt(r2); + double const e = exp(-r); + pair_m_[p] = (1 + r + r * r / 3) * e; + // d k / d l_i = s (5/3) (1 + sqrt(5) rho) exp(-sqrt(5) rho) (x_i - y_i)^2 / l_i^3, with r = sqrt(5) rho: + pair_dm_[p] = outputscale_ * (5.0 / 3) * (1 + r) * e; + row[b] = outputscale_ * pair_m_[p]; + } + row[a] = outputscale_ + noise_; + } + return turbo_cholesky(L_.data(), N); + } + + // The loss, -log p(y | theta)/N (as GPyTorch's ExactMarginalLogLikelihood, up to the constant log(2 pi)/2), and + // its gradient with respect to theta. Leaves L_ and alpha_ = (K + noise I)^{-1} (y - mean). Returns false if the + // kernel matrix cannot be factored. + bool loss_and_gradient(std::vector const &theta, double &loss, std::vector &grad) { + using std::log; + size_t const N = N_; + if (!factor(theta)) { + return false; + } + alpha_.resize(N); + for (size_t a = 0; a < N; ++a) { + alpha_[a] = yt_[a] - mean_; + } + turbo_cholesky_solve(L_.data(), N, alpha_.data()); + loss = 0; + for (size_t a = 0; a < N; ++a) { + loss += (yt_[a] - mean_) * alpha_[a] / 2 + log(L_[a * N + a]); + } + loss /= static_cast(N); + Z_.resize(N * N); + Kinv_.resize(N * N); + turbo_cholesky_inverse(L_.data(), N, Z_.data(), Kinv_.data()); + // d loss / d theta = -tr(W dK/d theta)/(2N), where W = alpha alpha^T - K^{-1}: + double trace_noise = 0; + double trace_scale = 0; + std::vector trace_l(n_, 0.0); + size_t p = 0; + for (size_t a = 0; a < N; ++a) { + double const *kinv = Kinv_.data() + a * N; + for (size_t b = 0; b < a; ++b, ++p) { + double const w = 2 * (alpha_[a] * alpha_[b] - kinv[b]); + trace_scale += w * pair_m_[p]; + double const wd = w * pair_dm_[p]; + double const *D = sqdiff_.data() + p * n_; + for (size_t i = 0; i < n_; ++i) { + trace_l[i] += wd * D[i]; + } + } + double const w = alpha_[a] * alpha_[a] - kinv[a]; + trace_noise += w; + trace_scale += w; + } + double sum_alpha = 0; + for (size_t a = 0; a < N; ++a) { + sum_alpha += alpha_[a]; + } + double const scale = -1 / (2 * static_cast(N)); + // The chain rule through value = lo + (hi - lo) sigmoid(raw): + auto dvalue = [](double raw, double const *bounds) { + double const sg = sigmoid(raw); + return (bounds[1] - bounds[0]) * sg * (1 - sg); + }; + for (size_t i = 0; i < n_; ++i) { + double const l = lengthscale_[i]; + grad[i] = scale * trace_l[i] / (l * l * l) * dvalue(theta[i], lengthscale_bounds_); + } + grad[n_] = scale * trace_scale * dvalue(theta[n_], outputscale_bounds_); + grad[n_ + 1] = scale * trace_noise * dvalue(theta[n_ + 1], noise_bounds_); + grad[n_ + 2] = -sum_alpha / static_cast(N); + return true; + } + + // Fits the Gaussian process (Section C) to the training set by Adam on the loss above, over the hyperparameters + // transformed to the real line as in GPyTorch: value = lo + (hi - lo) sigmoid(raw). Starts from the reference's + // values every time (lengthscales 0.5, outputscale 1, noise 0.005, mean 0). Leaves L_ and alpha_ for the + // posterior. Returns false if there is no model. + bool fit() { + using std::isfinite; + using std::sqrt; + if (!finite_) { + return false; + } + size_t const N = N_; + sqdiff_.resize(N * (N - 1) / 2 * n_); + size_t p = 0; + for (size_t a = 0; a < N; ++a) { + for (size_t b = 0; b < a; ++b, ++p) { + for (size_t i = 0; i < n_; ++i) { + double const t = Xt_[a * n_ + i] - Xt_[b * n_ + i]; + sqdiff_[p * n_ + i] = t * t; + } + } + } + lengthscale_.resize(n_); + inv_l2_.resize(n_); + size_t const dim = n_ + 3; + std::vector theta(dim); + for (size_t i = 0; i < n_; ++i) { + theta[i] = to_raw(0.5, lengthscale_bounds_); + } + theta[n_] = to_raw(1.0, outputscale_bounds_); + theta[n_ + 1] = to_raw(0.005, noise_bounds_); + theta[n_ + 2] = 0; + std::vector m1(dim, 0.0), m2(dim, 0.0), grad(dim), last_good = theta; + // Adam with the PyTorch defaults and the reference's learning rate 0.1: + double const lr = 0.1, beta1 = 0.9, beta2 = 0.999, eps = 1e-8; + double beta1_t = 1, beta2_t = 1; + double loss; + // Set if Adam stops on hyperparameters whose kernel matrix cannot be factored, or whose loss or gradient is not + // finite; the model then uses the last hyperparameters that were fine. + bool bad = false; + for (size_t step = 0; step < params_.hyperparameter_steps; ++step) { + bool finite = loss_and_gradient(theta, loss, grad) && isfinite(loss); + for (size_t k = 0; k < dim && finite; ++k) { + finite = isfinite(grad[k]); + } + if (!finite) { + bad = true; + break; + } + last_good = theta; + beta1_t *= beta1; + beta2_t *= beta2; + for (size_t k = 0; k < dim; ++k) { + m1[k] = beta1 * m1[k] + (1 - beta1) * grad[k]; + m2[k] = beta2 * m2[k] + (1 - beta2) * grad[k] * grad[k]; + theta[k] -= lr * (m1[k] / (1 - beta1_t)) / (sqrt(m2[k] / (1 - beta2_t)) + eps); + } + } + // The model at the final (or last good) hyperparameters: + if (bad || !factor(theta)) { + if (!factor(last_good)) { + return false; + } + } + alpha_.resize(N); + for (size_t a = 0; a < N; ++a) { + alpha_[a] = yt_[a] - mean_; + } + turbo_cholesky_solve(L_.data(), N, alpha_.data()); + return true; + } + + // q joint posterior samples over the M candidates C (row-major M x n), and the argmin of each. + void thompson_sample(std::vector const &C, size_t M, std::vector &chosen) { + using std::sqrt; + size_t const N = N_; + // Kxc (N x M), the posterior mean, then V = L^{-1} Kxc in place: + V_.resize(N * M); + std::vector mu(M, mean_); + for (size_t a = 0; a < N; ++a) { + double const *xa = Xt_.data() + a * n_; + double *row = V_.data() + a * M; + for (size_t c = 0; c < M; ++c) { + row[c] = outputscale_ * turbo_matern52(xa, C.data() + c * n_, inv_l2_.data(), n_); + mu[c] += row[c] * alpha_[a]; + } + } + turbo_forward_solve(L_.data(), N, V_.data(), M); + // The transpose, M x N, so that the rows of the Gram matrix V^T V are contiguous dot products: + Vt_.resize(M * N); + for (size_t a = 0; a < N; ++a) { + for (size_t c = 0; c < M; ++c) { + Vt_[c * N + a] = V_[a * M + c]; + } + } + // The posterior covariance Kcc + noise I - V^T V, then its Cholesky factor, with jitter if it is needed: + S_.resize(M * M); + for (size_t a = 0; a < M; ++a) { + double const *xa = C.data() + a * n_; + double *row = S_.data() + a * M; + for (size_t b = 0; b < a; ++b) { + row[b] = outputscale_ * turbo_matern52(xa, C.data() + b * n_, inv_l2_.data(), n_); + } + row[a] = outputscale_ + noise_; + } + turbo_subtract_gram(S_.data(), M, Vt_.data(), N, N, 0, M); + // The noise bounds the covariance below, so jitter should never be needed; if it is, retry from a copy. + covariance_ = S_; + bool factored = turbo_cholesky(S_.data(), M); + for (double jitter = 1e-8; !factored && jitter < 1; jitter *= 10) { + S_ = covariance_; + for (size_t a = 0; a < M; ++a) { + S_[a * M + a] += jitter; + } + factored = turbo_cholesky(S_.data(), M); + } + std::normal_distribution normal(0, 1); + std::vector z(M); + std::vector taken(M, false); + for (size_t j = 0; j < q_; ++j) { + for (auto &zi : z) { + zi = normal(gen_); + } + size_t best = M; + double best_value = std::numeric_limits::infinity(); + for (size_t c = 0; c < M; ++c) { + // Without a factor (it should not happen, as the noise bounds the covariance below), the marginals: + double const f = factored ? mu[c] + turbo_dot(S_.data() + c * M, z.data(), c + 1) + : mu[c] + sqrt(outputscale_ + noise_) * z[c]; + if (!taken[c] && (best == M || f < best_value)) { + best = c; + best_value = f; + } + } + taken[best] = true; + chosen.push_back(best); + } + } + + Func const &cost_function_; + turbo_parameters const ¶ms_; + URBG &gen_; + ResultType target_value_; + std::atomic *cancellation_; + std::atomic *current_minimum_cost_; + std::vector> *queries_; + size_t n_; + size_t q_; + size_t max_function_calls_; + size_t initial_points_; + size_t candidates_; + size_t failure_tolerance_; + size_t calls_ = 0; + ArgumentContainer best_argument_; + ResultType best_value_ = std::numeric_limits::infinity(); + // The trust region: its base side length, counters, and the points evaluated since it started (unit cube): + double length_ = length_init_; + size_t successes_ = 0; + size_t failures_ = 0; + std::vector X_; + std::vector y_; + // The Gaussian process: training set, hyperparameters, and scratch. + size_t N_ = 0; + bool finite_ = false; + std::vector Xt_, yt_; + std::vector lengthscale_, inv_l2_; + double outputscale_ = 1, noise_ = 0.005, mean_ = 0; + std::vector sqdiff_, pair_m_, pair_dm_; + std::vector L_, Z_, Kinv_, alpha_; + std::vector V_, Vt_, S_, covariance_; +}; + +} // namespace detail + +BOOST_MATH_EXPORT template +ArgumentContainer turbo( + const Func cost_function, turbo_parameters const ¶ms, URBG &gen, + 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::turbo_engine engine(cost_function, params, gen, target_value, cancellation, + current_minimum_cost, queries); + return engine.optimize(); +} + +} BOOST_MATH_NAMESPACE_END // namespace boost::math::optimization +#endif diff --git a/module/CMakeLists.txt b/module/CMakeLists.txt index cc340d4ccd..6b3db3a649 100644 --- a/module/CMakeLists.txt +++ b/module/CMakeLists.txt @@ -84,6 +84,7 @@ set(BOOST_MATH_MODULE_TESTS test_random_search test_cma_es test_multilevel_coordinate_search + test_turbo test_quaternion_constexpr octonion_test_simple test_lambert_w_module diff --git a/module/math.cppm b/module/math.cppm index 337f2e5c78..37de4ac6ab 100644 --- a/module/math.cppm +++ b/module/math.cppm @@ -214,6 +214,7 @@ extern "C++" { #include #include #include +#include // Differentiation. finite_difference and lanczos_smoothing must precede // autodiff: their unqualified detail:: references become ambiguous once diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index b60ff7a8cf..249fcdfde1 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -1154,9 +1154,11 @@ test-suite interpolators : [ 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_turbo.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ run test_cma_es.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ] [ 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_turbo_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_turbo_incl.cpp b/test/compile_test/test_turbo_incl.cpp new file mode 100644 index 0000000000..c7ae4c9aef --- /dev/null +++ b/test/compile_test/test_turbo_incl.cpp @@ -0,0 +1,22 @@ +// 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::turbo_parameters> params; + std::mt19937_64 gen(0); + auto v = boost::math::optimization::turbo(func, params, gen); +} diff --git a/test/test_turbo.cpp b/test/test_turbo.cpp new file mode 100644 index 0000000000..47c412b8c7 --- /dev/null +++ b/test/test_turbo.cpp @@ -0,0 +1,290 @@ +/* + * 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 +#include +#include +#include +#include +#include +#endif +using boost::math::optimization::turbo; +using boost::math::optimization::turbo_parameters; +using boost::math::optimization::validate_turbo_parameters; + +// The random number generators of the standard library are fixed, but its distributions are not, so the same seed +// gives different runs on different standard libraries. The tolerances below hold for any seed we tried. + +template +void check_in_bounds(std::vector> const &queries, ArgType const &lo, ArgType const &hi) { + bool inside = true; + for (auto const &q : queries) { + for (size_t i = 0; i < lo.size(); ++i) { + inside = inside && lo[i] <= q.first[i] && q.first[i] <= hi[i]; + } + } + CHECK_TRUE(inside); +} + +template Real branin(std::array 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; +} + +template void test_branin() { + std::cout << "Testing TuRBO on the Branin function . . .\n"; + using ArgType = std::array; + turbo_parameters params; + params.lower_bounds = {-5, 0}; + params.upper_bounds = {10, 15}; + params.max_function_calls = 200; + std::mt19937_64 gen(12345); + std::vector> queries; + std::atomic current_minimum_cost = std::numeric_limits::quiet_NaN(); + auto x = turbo(branin, params, gen, std::numeric_limits::quiet_NaN(), nullptr, ¤t_minimum_cost, + &queries); + // The global minimum value, attained at (-pi, 12.275), (pi, 2.275) and (9.42478, 2.475): + CHECK_ABSOLUTE_ERROR(Real(0.39788735772973816), branin(x), Real(1e-2)); + CHECK_EQUAL(queries.size(), size_t(200)); + check_in_bounds(queries, params.lower_bounds, params.upper_bounds); + // The returned point and the current minimum are those of the best query: + Real best = std::numeric_limits::infinity(); + for (auto const &q : queries) { + CHECK_EQUAL(branin(q.first), q.second); + best = (std::min)(best, q.second); + } + CHECK_EQUAL(branin(x), best); + CHECK_EQUAL(current_minimum_cost.load(), best); +} + +template void test_six_hump_camel() { + std::cout << "Testing TuRBO on the six-hump camel function . . .\n"; + using ArgType = std::vector; + 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; + }; + turbo_parameters params; + params.lower_bounds = {-3, -2}; + params.upper_bounds = {3, 2}; + params.max_function_calls = 200; + std::mt19937_64 gen(2); + auto x = turbo(camel, params, gen); + CHECK_EQUAL(x.size(), size_t(2)); + CHECK_ABSOLUTE_ERROR(Real(-1.0316284534898774), camel(x), Real(1e-2)); +} + +void test_hartmann6() { + std::cout << "Testing TuRBO on the Hartmann 6 function . . .\n"; + using ArgType = std::vector; + auto hartmann = [](ArgType const &x) { + static const double A[4][6] = {{10, 3, 17, 3.5, 1.7, 8}, {0.05, 10, 17, 0.1, 8, 14}, {3, 3.5, 1.7, 10, 17, 8}, + {17, 8, 0.05, 10, 0.1, 14}}; + static const double P[4][6] = {{0.1312, 0.1696, 0.5569, 0.0124, 0.8283, 0.5886}, + {0.2329, 0.4135, 0.8307, 0.3736, 0.1004, 0.9991}, + {0.2348, 0.1451, 0.3522, 0.2883, 0.3047, 0.6650}, + {0.4047, 0.8828, 0.8732, 0.5743, 0.1091, 0.0381}}; + static const double c[4] = {1, 1.2, 3, 3.2}; + double s = 0; + for (int i = 0; i < 4; ++i) { + double e = 0; + for (int j = 0; j < 6; ++j) { + e += A[i][j] * (x[j] - P[i][j]) * (x[j] - P[i][j]); + } + s -= c[i] * std::exp(-e); + } + return s; + }; + turbo_parameters params; + params.lower_bounds.resize(6, 0); + params.upper_bounds.resize(6, 1); + params.max_function_calls = 300; + // Fewer candidates than the default 600, to keep the test fast: + params.candidates = 300; + std::mt19937_64 gen(3); + auto x = turbo(hartmann, params, gen); + // The global minimum is -3.32237, and the next best local minimum -3.2032. TuRBO-1 is a local method with restarts, + // and about one run in five ends in the second basin after 300 calls (one in twenty after 500). + CHECK_LE(hartmann(x), -3.19); +} + +template void test_rosenbrock() { + std::cout << "Testing TuRBO on the Rosenbrock function . . .\n"; + using ArgType = std::array; + turbo_parameters params; + params.lower_bounds = {-1.5, -0.5}; + params.upper_bounds = {2.5, 3}; + params.max_function_calls = 400; + std::mt19937_64 gen(4); + auto x = turbo(rosenbrock_saddle, params, gen); + // The function varies over five orders of magnitude on this box, which is hard for the Gaussian process. Over 20 + // seeds, the median is 0.006 and the worst 0.12. + CHECK_LE(rosenbrock_saddle(x), Real(0.25)); +} + +// The sphere of test_functions_for_optimization.hpp is NaN outside the unit ball, which covers most of these boxes. +// It also tests a return type different from the argument type. +void test_sphere() { + std::cout << "Testing TuRBO on the sphere, on asymmetric boxes . . .\n"; + using ArgType = std::vector; + turbo_parameters params; + params.lower_bounds = {-1, -0.4f, -0.7f}; + params.upper_bounds = {0.6f, 1, 0.3f}; + params.max_function_calls = 150; + std::mt19937_64 gen(5); + std::vector> queries; + auto x = turbo(sphere, params, gen, std::numeric_limits::quiet_NaN(), nullptr, nullptr, &queries); + CHECK_LE(sphere(x), 1e-3); + check_in_bounds(queries, params.lower_bounds, params.upper_bounds); + + // A box that excludes the minimizer, which is then at the corner (0.2, 0.1): + params.lower_bounds = {0.2f, 0.1f}; + params.upper_bounds = {0.9f, 0.5f}; + queries.clear(); + x = turbo(sphere, params, gen, std::numeric_limits::quiet_NaN(), nullptr, nullptr, &queries); + CHECK_LE(sphere(x) - 0.05, 2e-3); + check_in_bounds(queries, params.lower_bounds, params.upper_bounds); +} + +void test_termination() { + std::cout << "Testing TuRBO termination: budget, cancellation, target value, batches . . .\n"; + using ArgType = std::array; + turbo_parameters params; + params.lower_bounds = {0.5, 0.5}; + params.upper_bounds = {2.048, 2.048}; + params.max_function_calls = 37; + std::mt19937_64 gen(6); + std::vector> queries; + auto const nan = std::numeric_limits::quiet_NaN(); + turbo(rosenbrock_saddle, params, gen, nan, nullptr, nullptr, &queries); + CHECK_EQUAL(queries.size(), size_t(37)); + + // The budget is exact even when it ends a batch early: + params.batch_size = 3; + queries.clear(); + turbo(rosenbrock_saddle, params, gen, nan, nullptr, nullptr, &queries); + CHECK_EQUAL(queries.size(), size_t(37)); + check_in_bounds(queries, params.lower_bounds, params.upper_bounds); + params.batch_size = 1; + + // The default budget is 100 n: + params.max_function_calls = 0; + queries.clear(); + turbo(rosenbrock_saddle, params, gen, nan, nullptr, nullptr, &queries); + CHECK_EQUAL(queries.size(), size_t(200)); + + // Cancellation from inside the cost function stops it after exactly that call: + std::atomic cancel = false; + size_t calls = 0; + auto f = [&](ArgType const &x) { + if (++calls == 10) { + cancel = true; + } + return rosenbrock_saddle(x); + }; + queries.clear(); + turbo(f, params, gen, nan, &cancel, nullptr, &queries); + CHECK_EQUAL(queries.size(), size_t(10)); + // Cancelled before the first call: no queries, and the midpoint of the box comes back. + queries.clear(); + auto x = turbo(f, params, gen, nan, &cancel, nullptr, &queries); + CHECK_EQUAL(queries.size(), size_t(0)); + CHECK_EQUAL(x[0], (params.lower_bounds[0] + params.upper_bounds[0]) / 2); + CHECK_EQUAL(x[1], (params.lower_bounds[1] + params.upper_bounds[1]) / 2); + + // The run stops at the first call that reaches the target value: + double const target = 0.05; + queries.clear(); + x = turbo(rosenbrock_saddle, params, gen, target, nullptr, nullptr, &queries); + CHECK_LE(rosenbrock_saddle(x), target); + CHECK_LE(queries.back().second, target); + for (size_t k = 0; k + 1 < queries.size(); ++k) { + CHECK_GE(queries[k].second, target); + } + CHECK_LE(queries.size(), size_t(199)); + + // The initial guess is the first call; if it reaches the target, it is the only one: + ArgType guess{1, 1}; + params.initial_guess = &guess; + queries.clear(); + x = turbo(rosenbrock_saddle, params, gen, 0.0, nullptr, nullptr, &queries); + CHECK_EQUAL(queries.size(), size_t(1)); + CHECK_EQUAL(x[0], 1.0); + CHECK_EQUAL(x[1], 1.0); +} + +void test_validation() { + std::cout << "Testing TuRBO parameter validation . . .\n"; + using ArgType = std::array; + turbo_parameters params; + params.lower_bounds = {-1, -1}; + params.upper_bounds = {1, 1}; + validate_turbo_parameters(params); + std::mt19937_64 gen(7); + + auto bad = params; + bad.upper_bounds[1] = -1; + CHECK_THROW(validate_turbo_parameters(bad), std::domain_error); + auto run = [&]() { return turbo(rosenbrock_saddle, bad, gen); }; + CHECK_THROW(run(), std::domain_error); + + // TuRBO scales the box to the unit cube, so it must have positive width: + bad = params; + bad.upper_bounds[0] = bad.lower_bounds[0]; + CHECK_THROW(validate_turbo_parameters(bad), std::domain_error); + + bad = params; + bad.upper_bounds[0] = std::numeric_limits::infinity(); + CHECK_THROW(validate_turbo_parameters(bad), std::domain_error); + + ArgType outside{2, 0}; + bad = params; + bad.initial_guess = &outside; + CHECK_THROW(validate_turbo_parameters(bad), std::domain_error); + + bad = params; + bad.batch_size = 0; + CHECK_THROW(validate_turbo_parameters(bad), std::invalid_argument); + + bad = params; + bad.batch_size = 4; + bad.candidates = 3; + CHECK_THROW(validate_turbo_parameters(bad), std::invalid_argument); + + bad = params; + bad.max_training_points = 0; + CHECK_THROW(validate_turbo_parameters(bad), std::invalid_argument); +} + +int main() { + test_branin(); + test_branin(); + test_six_hump_camel(); + test_six_hump_camel(); + test_hartmann6(); + test_rosenbrock(); + test_rosenbrock(); + test_sphere(); + test_termination(); + test_validation(); + return boost::math::test::report_errors(); +}