diff --git a/doc/optimization/cma_es.qbk b/doc/optimization/cma_es.qbk index 144c1bc597..fa75902a7a 100644 --- a/doc/optimization/cma_es.qbk +++ b/doc/optimization/cma_es.qbk @@ -40,6 +40,8 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) The `cma_es` optimizer searches for a global minimum of a function. Our implementation attempts to follow "The CMA evolution strategy: A tutorial" exactly. +Each sample is projected onto the bounds before the cost function is called, so the cost function is only evaluated inside the box, +and the mean of the search distribution is kept inside the box as well. [heading Parameters] diff --git a/doc/overview/modules.qbk b/doc/overview/modules.qbk index 8fe6cbcf82..d41d0c5a73 100644 --- a/doc/overview/modules.qbk +++ b/doc/overview/modules.qbk @@ -42,7 +42,7 @@ The `cmake-module-test` job in the library's CI runs this configuration on every * The `user_error` policy is not supported through the module, because a user supplied `boost::math::policies::user_..._error` definition cannot attach to the module. * A program should consume Boost.Math either through the module or through textual includes, not both. The only headers that are safe to include textually alongside the import are `boost/math/tools/config.hpp` and `boost/math/tools/assert.hpp`, which the test support headers rely on. -* Excluded from the module: the TR1/C99 interfaces and `common_factor` headers, `boost/math/bindings`, `boost/math/concepts`, `boost/math/cstdfloat`, all deprecated headers, and `optimization/cma_es.hpp`, which requires Eigen. +* Excluded from the module: the TR1/C99 interfaces and `common_factor` headers, `boost/math/bindings`, `boost/math/concepts`, `boost/math/cstdfloat`, and all deprecated headers. These remain available as ordinary includes. * `__float128` support (`BOOST_MATH_USE_FLOAT128`) is not enabled in the module build. * On libc++ the parallel execution policy overloads in statistics compile out (no `` support), matching the ordinary header behavior on that standard library. diff --git a/include/boost/math/optimization/cma_es.hpp b/include/boost/math/optimization/cma_es.hpp index fac3950d28..e4f4d75def 100644 --- a/include/boost/math/optimization/cma_es.hpp +++ b/include/boost/math/optimization/cma_es.hpp @@ -19,11 +19,6 @@ #endif #include #include -#if __has_include() -#include -#else -#error "CMA-ES requires Eigen." -#endif // Follows the notation in: // https://arxiv.org/pdf/1604.00772.pdf @@ -40,6 +35,224 @@ BOOST_MATH_NAMESPACE_BEGIN namespace optimization { +namespace detail { + +// Computes A = V diag(eigenvalues) V^T for a symmetric n x n matrix A stored row-major, +// by Householder tridiagonalization followed by the implicit QL algorithm; +// see Wilkinson & Reinsch, Handbook for Automatic Computation, Vol. II (tred2 and tql2), +// and Golub & Van Loan, Matrix Computations, 4th ed., Section 8.3. +// V is returned row-major, with the eigenvectors as its columns. +// Returns false if some eigenvalue did not converge within max_iterations QL iterations. +template +bool symmetric_eigendecomposition(std::vector A, size_t n, std::vector &V, std::vector &eigenvalues, size_t max_iterations = 30) { + using std::abs; + using std::hypot; + using std::sqrt; + using std::max; + V = std::move(A); + std::vector &d = eigenvalues; + d.resize(n); + std::vector e(n, Real(0)); + // Householder reduction to tridiagonal form; d is the diagonal and e the subdiagonal: + for (size_t j = 0; j < n; ++j) { + d[j] = V[(n - 1)*n + j]; + } + for (size_t i = n - 1; i > 0; --i) { + Real scale = 0; + Real h = 0; + for (size_t k = 0; k < i; ++k) { + scale += abs(d[k]); + } + if (scale == 0) { + e[i] = d[i - 1]; + for (size_t j = 0; j < i; ++j) { + d[j] = V[(i - 1)*n + j]; + V[i*n + j] = 0; + V[j*n + i] = 0; + } + } + else { + for (size_t k = 0; k < i; ++k) { + d[k] /= scale; + h += d[k]*d[k]; + } + Real f = d[i - 1]; + Real g = sqrt(h); + if (f > 0) { + g = -g; + } + e[i] = scale*g; + h -= f*g; + d[i - 1] = f - g; + for (size_t j = 0; j < i; ++j) { + e[j] = 0; + } + for (size_t j = 0; j < i; ++j) { + f = d[j]; + V[j*n + i] = f; + g = e[j] + V[j*n + j]*f; + for (size_t k = j + 1; k < i; ++k) { + g += V[k*n + j]*d[k]; + e[k] += V[k*n + j]*f; + } + e[j] = g; + } + f = 0; + for (size_t j = 0; j < i; ++j) { + e[j] /= h; + f += e[j]*d[j]; + } + const Real hh = f/(h + h); + for (size_t j = 0; j < i; ++j) { + e[j] -= hh*d[j]; + } + for (size_t j = 0; j < i; ++j) { + f = d[j]; + g = e[j]; + for (size_t k = j; k < i; ++k) { + V[k*n + j] -= f*e[k] + g*d[k]; + } + d[j] = V[(i - 1)*n + j]; + V[i*n + j] = 0; + } + } + d[i] = h; + } + // Accumulate the transformations: + for (size_t i = 0; i + 1 < n; ++i) { + V[(n - 1)*n + i] = V[i*n + i]; + V[i*n + i] = 1; + const Real h = d[i + 1]; + if (h != 0) { + for (size_t k = 0; k <= i; ++k) { + d[k] = V[k*n + i + 1]/h; + } + for (size_t j = 0; j <= i; ++j) { + Real g = 0; + for (size_t k = 0; k <= i; ++k) { + g += V[k*n + i + 1]*V[k*n + j]; + } + for (size_t k = 0; k <= i; ++k) { + V[k*n + j] -= g*d[k]; + } + } + } + for (size_t k = 0; k <= i; ++k) { + V[k*n + i + 1] = 0; + } + } + for (size_t j = 0; j < n; ++j) { + d[j] = V[(n - 1)*n + j]; + V[(n - 1)*n + j] = 0; + } + V[(n - 1)*n + n - 1] = 1; + // The QL rotations act on columns of V, so apply them to the rows of V^T, which are contiguous: + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < i; ++j) { + std::swap(V[i*n + j], V[j*n + i]); + } + } + // Implicit QL iterations on the tridiagonal matrix: + for (size_t i = 1; i < n; ++i) { + e[i - 1] = e[i]; + } + e[n - 1] = 0; + Real f = 0; + Real tst1 = 0; + const Real eps = std::numeric_limits::epsilon(); + for (size_t l = 0; l < n; ++l) { + // Find a small subdiagonal element: + tst1 = (max)(tst1, abs(d[l]) + abs(e[l])); + size_t m = l; + while (m < n - 1 && !(abs(e[m]) <= eps*tst1)) { + ++m; + } + // If m == l, d[l] is already an eigenvalue; otherwise, iterate: + size_t iterations = 0; + while (m > l && abs(e[l]) > eps*tst1) { + if (iterations++ == max_iterations) { + break; + } + // Compute the implicit shift: + Real g = d[l]; + Real p = (d[l + 1] - g)/(2*e[l]); + Real r = hypot(p, Real(1)); + if (p < 0) { + r = -r; + } + d[l] = e[l]/(p + r); + d[l + 1] = e[l]*(p + r); + const Real dl1 = d[l + 1]; + Real h = g - d[l]; + for (size_t i = l + 2; i < n; ++i) { + d[i] -= h; + } + f += h; + p = d[m]; + Real c = 1; + Real c2 = c; + Real c3 = c; + const Real el1 = e[l + 1]; + Real s = 0; + Real s2 = 0; + for (size_t i = m; i-- > l;) { + c3 = c2; + c2 = c; + s2 = s; + g = c*e[i]; + h = c*p; + r = hypot(p, e[i]); + e[i + 1] = s*r; + s = e[i]/r; + c = p/r; + p = c*d[i] - s*g; + d[i + 1] = h + s*(c*g + s*d[i]); + for (size_t k = 0; k < n; ++k) { + h = V[(i + 1)*n + k]; + V[(i + 1)*n + k] = s*V[i*n + k] + c*h; + V[i*n + k] = c*V[i*n + k] - s*h; + } + } + p = -s*s2*c3*el1*e[l]/dl1; + e[l] = s*p; + d[l] = c*p; + } + if (iterations > max_iterations) { + return false; + } + d[l] += f; + e[l] = 0; + } + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < i; ++j) { + std::swap(V[i*n + j], V[j*n + i]); + } + } + return true; +} + +// y = Mx, or y = M^{T}x if transpose is true, for the n x n matrix M stored row-major. +// Both accumulate along the rows of M, so the inner loop is over contiguous memory: +template +void matrix_vector_product(std::vector const &M, std::vector const &x, std::vector &y, bool transpose = false) { + const size_t n = x.size(); + for (size_t i = 0; i < n; ++i) { + y[i] = 0; + } + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < n; ++j) { + if (transpose) { + y[j] += M[i*n + j]*x[i]; + } + else { + y[i] += M[i*n + j]*x[j]; + } + } + } +} + +} // namespace detail + BOOST_MATH_EXPORT template struct cma_es_parameters { using Real = typename ArgumentContainer::value_type; using DimensionlessReal = decltype(Real()/Real()); @@ -117,9 +330,9 @@ ArgumentContainer cma_es( std::atomic lowest_cost = std::numeric_limits::infinity(); ArgumentContainer best_vector; // p_{c} := evolution path, equation (24) of the arxiv review: - Eigen::Vector p_c(n); + std::vector p_c(n); // p_{\sigma} := conjugate evolution path, equation (31) of the arxiv review: - Eigen::Vector p_sigma(n); + std::vector p_sigma(n); if constexpr (detail::has_resize_v) { best_vector.resize(n, std::numeric_limits::quiet_NaN()); } @@ -183,7 +396,11 @@ ArgumentContainer cma_es( weights[i] = min_alpha*w_prime[i]/abs(negative_weight_sum); } // mu:= number of parents, lambda := number of offspring. - Eigen::Matrix C = Eigen::Matrix::Identity(n, n); + // C, and B below, are stored row-major: + std::vector C(n*n, DimensionlessReal(0)); + for (size_t i = 0; i < n; ++i) { + C[i*n + i] = 1; + } ArgumentContainer mean_vector; // See the footnote in Figure 6 of the arxiv review: // We should consider the more robust initialization described there. . . @@ -223,11 +440,13 @@ ArgumentContainer cma_es( #endif size_t generation = 0; - std::vector> ys(params.population_size); + std::vector> ys(params.population_size); std::vector xs(params.population_size); std::vector costs(params.population_size, std::numeric_limits::quiet_NaN()); - Eigen::Vector weighted_avg_y(n); - Eigen::Vector z(n); + std::vector weighted_avg_y(n); + std::vector z(n); + std::vector B; + std::vector D; if constexpr (detail::has_resize_v) { for (auto & x : xs) { x.resize(n, std::numeric_limits::quiet_NaN()); @@ -247,24 +466,29 @@ ArgumentContainer cma_es( // max(1, floor(1/(10n(c_1+c_mu)))) generations." // Note that sigma can be dimensionless, in which case C carries the units. // This is a weird decision-we're not gonna do that! - Eigen::SelfAdjointEigenSolver> eigensolver(C); - if (eigensolver.info() != Eigen::Success) { + if (!detail::symmetric_eigendecomposition(C, n, B, D)) { std::ostringstream oss; oss << __FILE__ << ":" << __LINE__ << ":" << __func__; oss << ": Could not decompose the covariance matrix as BDB^{T}."; throw std::logic_error(oss.str()); } - Eigen::Matrix B = eigensolver.eigenvectors(); - // Eigen returns D^2, in the notation of the survey: - auto D = eigensolver.eigenvalues(); - // So make it better: + // The eigenvalues are D^2, in the notation of the survey: for (auto & d : D) { - if (d <= 0 || isnan(d)) { + if (d <= 0 || !isfinite(d)) { std::ostringstream oss; oss << __FILE__ << ":" << __LINE__ << ":" << __func__; oss << ": The covariance matrix is not positive definite. This breaks the evolution path computation downstream.\n"; - oss << "C=\n" << C << "\n"; - oss << "Eigenvalues: " << D; + oss << "C=\n"; + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < n; ++j) { + oss << C[i*n + j] << " "; + } + oss << "\n"; + } + oss << "Eigenvalues: "; + for (auto & e : D) { + oss << e << " "; + } throw std::domain_error(oss.str()); } d = sqrt(d); @@ -282,15 +506,20 @@ ArgumentContainer cma_es( for (size_t i = 0; i < n; ++i) { z[i] = dis(gen); } - Eigen::Vector Dz(n); for (size_t i = 0; i < n; ++i) { - Dz[i] = D[i]*z[i]; + z[i] *= D[i]; } - y = B*Dz; + detail::matrix_vector_product(B, z, y); for (size_t i = 0; i < n; ++i) { BOOST_MATH_ASSERT(!isnan(mean_vector[i])); BOOST_MATH_ASSERT(!isnan(y[i])); x[i] = mean_vector[i] + sigma*y[i]; // equation (40) of Figure 6. + // Repair by projecting onto the bounds, and update with the repaired step; + // see "Boundaries and constraints" in Appendix B of the arxiv review: + if (x[i] < params.lower_bounds[i] || x[i] > params.upper_bounds[i]) { + x[i] = (max)(params.lower_bounds[i], (min)(params.upper_bounds[i], x[i])); + y[i] = (x[i] - mean_vector[i])/sigma; + } } costs[k] = cost_function(x); if (resample_counter++ == 50) { @@ -335,14 +564,24 @@ ArgumentContainer cma_es( } // Equation (42), Figure 6: for (size_t j = 0; j < n; ++j) { + const Real previous_mean = mean_vector[j]; mean_vector[j] = mean_vector[j] + c_m*sigma*weighted_avg_y[j]; + // For c_m <= 1, the new mean is a convex combination of the old mean and points in the bounds, + // so it can only leave the bounds by roundoff or if c_m > 1. Then project it, + // and give the evolution paths the step that the mean actually took: + if (mean_vector[j] < params.lower_bounds[j] || mean_vector[j] > params.upper_bounds[j]) { + mean_vector[j] = (max)(params.lower_bounds[j], (min)(params.upper_bounds[j], mean_vector[j])); + weighted_avg_y[j] = (mean_vector[j] - previous_mean)/(c_m*sigma); + } } // Equation (43), Figure 6: Start with C^{-1/2}_{w} - Eigen::Vector inv_D_B_transpose_y = B.transpose()*weighted_avg_y; - for (long j = 0; j < inv_D_B_transpose_y.size(); ++j) { + std::vector inv_D_B_transpose_y(n); + detail::matrix_vector_product(B, weighted_avg_y, inv_D_B_transpose_y, true); + for (size_t j = 0; j < n; ++j) { inv_D_B_transpose_y[j] /= D[j]; } - Eigen::Vector C_inv_sqrt_y_avg = B*inv_D_B_transpose_y; + std::vector C_inv_sqrt_y_avg(n); + detail::matrix_vector_product(B, inv_D_B_transpose_y, C_inv_sqrt_y_avg); // Equation (43), Figure 6: DimensionlessReal p_sigma_norm = 0; for (size_t j = 0; j < n; ++j) { @@ -361,7 +600,9 @@ ArgumentContainer cma_es( h_sigma = 1; } // Equation (45), Figure 6: - p_c = (1-c_c)*p_c + h_sigma*sqrt(c_c*(2-c_c)*mu_eff)*weighted_avg_y; + for (size_t j = 0; j < n; ++j) { + p_c[j] = (1-c_c)*p_c[j] + h_sigma*sqrt(c_c*(2-c_c)*mu_eff)*weighted_avg_y[j]; + } DimensionlessReal delta_h_sigma = (1-h_sigma)*c_c*(2-c_c); DimensionlessReal weight_sum = 0; for (auto & w : weights) { @@ -369,21 +610,42 @@ ArgumentContainer cma_es( } // Equation (47), Figure 6: DimensionlessReal K = (1 + c_1*delta_h_sigma - c_1 - c_mu*weight_sum); - // Can these operations be sped up using `.selfadjointView`? - // Maybe: A.selfadjointView().rankUpdate(p_c, c_1);? - C = K*C + c_1*p_c*p_c.transpose(); + // C is symmetric, so only its upper triangle is updated, and then copied to the lower triangle: + auto add_outer_product = [&](DimensionlessReal alpha, std::vector const & v) { + for (size_t r = 0; r < n; ++r) { + for (size_t c = r; c < n; ++c) { + C[r*n + c] += alpha*v[r]*v[c]; + } + } + }; + for (auto & c : C) { + c *= K; + } + add_outer_product(c_1, p_c); // Incorporate positive weights of Equation (46): for (size_t i = 0; i < params.population_size/2; ++i) { - C += c_mu*weights[i]*ys[indices[i]]*ys[indices[i]].transpose(); + add_outer_product(c_mu*weights[i], ys[indices[i]]); } + std::vector D_inv_BTy(n); for (size_t i = params.population_size/2; i < params.population_size; ++i) { - Eigen::Vector D_inv_BTy = B.transpose()*ys[indices[i]]; + auto const & y = ys[indices[i]]; + detail::matrix_vector_product(B, y, D_inv_BTy, true); + DimensionlessReal squared_norm = 0; for (size_t j = 0; j < n; ++j) { D_inv_BTy[j] /= D[j]; + squared_norm += D_inv_BTy[j]*D_inv_BTy[j]; + } + // A repaired step can be exactly zero: + if (squared_norm == 0) { + continue; } - DimensionlessReal squared_norm = D_inv_BTy.squaredNorm(); DimensionlessReal K2 = c_mu*weights[i]/squared_norm; - C += K2*ys[indices[i]]*ys[indices[i]].transpose(); + add_outer_product(K2, y); + } + for (size_t r = 0; r < n; ++r) { + for (size_t c = 0; c < r; ++c) { + C[r*n + c] = C[c*n + r]; + } } } while (generation++ < params.max_generations); diff --git a/module/CMakeLists.txt b/module/CMakeLists.txt index e0167a5f88..cc340d4ccd 100644 --- a/module/CMakeLists.txt +++ b/module/CMakeLists.txt @@ -82,6 +82,8 @@ set(BOOST_MATH_MODULE_TESTS test_ljung_box test_differential_evolution test_random_search + test_cma_es + test_multilevel_coordinate_search test_quaternion_constexpr octonion_test_simple test_lambert_w_module diff --git a/module/math.cppm b/module/math.cppm index 47fa0d423f..337f2e5c78 100644 --- a/module/math.cppm +++ b/module/math.cppm @@ -203,13 +203,15 @@ extern "C++" { #include #include -// Optimization (cma_es is excluded: it requires Eigen) +// Optimization +#include #include #include #include #include #include #include +#include #include #include diff --git a/nanobind/CMakeLists.txt b/nanobind/CMakeLists.txt index 69c096565f..63ebc61340 100644 --- a/nanobind/CMakeLists.txt +++ b/nanobind/CMakeLists.txt @@ -51,11 +51,6 @@ set(BOOST_MATH_NB_SOURCES nanobind_add_module(_ext NOMINSIZE ${BOOST_MATH_NB_SOURCES}) target_include_directories(_ext PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/../include" "${CMAKE_CURRENT_SOURCE_DIR}/src") target_compile_definitions(_ext PRIVATE BOOST_MATH_STANDALONE) -# Optional: cma_es needs Eigen. -find_package(Eigen3 CONFIG QUIET) -if(Eigen3_FOUND) - target_link_libraries(_ext PRIVATE Eigen3::Eigen) -endif() # Optional: boost.math.ulps needs Boost.Multiprecision (standalone, which needs Boost.Config) and zlib. Only the # sibling checkouts (or the cache variables) are searched, never a system Boost of another release. find_path(BOOST_MATH_NB_MULTIPRECISION_INCLUDE boost/multiprecision/cpp_bin_float.hpp diff --git a/nanobind/src/optimization.cpp b/nanobind/src/optimization.cpp index 9e9e3b5791..c491804d61 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, multilevel_coordinate_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) +// 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. @@ -12,14 +12,11 @@ #include "bindings.hpp" +#include #include #include #include #include -#if __has_include() -#include -#define BOOST_MATH_NB_HAS_CMA_ES -#endif #include #include @@ -192,7 +189,6 @@ void bind_optimization(nb::module_& m) 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, std::uint64_t seed) { @@ -210,7 +206,6 @@ void bind_optimization(nb::module_& m) }, arg("f"), arg("lower_bounds"), arg("upper_bounds"), arg("initial_guess") = nb::none(), arg("max_generations") = 1000, arg("population_size") = 0, arg("learning_rate") = 1.0, arg("target_value") = no_target, arg("seed") = 0, "optimization::cma_es: the argument in the box that minimizes f"); -#endif } } // namespace bmpy diff --git a/nanobind/tests/test_smoke.py b/nanobind/tests/test_smoke.py index af8f7c7e85..c100a70157 100644 --- a/nanobind/tests/test_smoke.py +++ b/nanobind/tests/test_smoke.py @@ -292,9 +292,8 @@ 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.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)) + bm.random_search(f, lo, hi, max_function_calls=20000), bm.multilevel_coordinate_search(f, lo, hi), + bm.cma_es(f, lo, hi, max_generations=200)] for x in runs: assert f(x) < 1e-2 def bad(x): diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 27ca2eda55..b60ff7a8cf 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -1154,7 +1154,7 @@ 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_cma_es.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] + [ 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_differential_evolution_incl.cpp : [ requires cxx17_if_constexpr cxx17_std_apply ] ] diff --git a/test/cmake_namespace_test/CMakeLists.txt b/test/cmake_namespace_test/CMakeLists.txt index 105bc3a511..e9b6cd6950 100644 --- a/test/cmake_namespace_test/CMakeLists.txt +++ b/test/cmake_namespace_test/CMakeLists.txt @@ -51,7 +51,6 @@ if(BOOST_MATH_NAMESPACE_TEST_HEADERS) boost/math/constants/calculate_constants.hpp boost/math/cstdfloat/cstdfloat_complex_std.hpp boost/math/differentiation/autodiff_cpp11.hpp - boost/math/optimization/cma_es.hpp boost/math/tools/signal_statistics.hpp ) diff --git a/test/test_cma_es.cpp b/test/test_cma_es.cpp index f221094223..bdc3bde124 100644 --- a/test/test_cma_es.cpp +++ b/test/test_cma_es.cpp @@ -5,14 +5,22 @@ * 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" -#include +#ifndef BOOST_MATH_BUILD_MODULE #include #include #include #include #include +#include +#endif using std::abs; using boost::math::optimization::cma_es; @@ -53,8 +61,9 @@ template void test_ackley() { CHECK_EQUAL(local_minima[1], Real(0)); CHECK_LE(size_t(1), queries.size()); for (auto const & q : queries) { + // Not exact: with -ffp-contract=fast, ackley can be fused differently where it is inlined. auto expected = ackley(q.first); - CHECK_EQUAL(expected, q.second); + CHECK_ULP_CLOSE(expected, q.second, 4); } } @@ -146,6 +155,111 @@ void test_beale() { CHECK_ABSOLUTE_ERROR(Real(1)/Real(2), local_minima[1], Real(0.1)); } +// The unconstrained minimum of the quadratic is outside the bounds, +// and the constrained minimum on [1, 2]^4 is at (1, 1, 1, 1): +template +void test_bounds() { + std::cout << "Testing that CMA-ES respects the bounds . . .\n"; + using ArgType = std::vector; + auto params = cma_es_parameters(); + params.lower_bounds.resize(4, Real(1)); + params.upper_bounds.resize(4, Real(2)); + params.max_generations = 200; + std::mt19937_64 gen(12345); + size_t calls = 0; + size_t out_of_bounds_calls = 0; + auto f = [&](ArgType const & x) { + ++calls; + for (size_t i = 0; i < x.size(); ++i) { + if (!(params.lower_bounds[i] <= x[i] && x[i] <= params.upper_bounds[i])) { + ++out_of_bounds_calls; + break; + } + } + auto y = x; + return quadratic(y); + }; + auto local_minima = cma_es(f, params, gen); + CHECK_LE(size_t(1), calls); + CHECK_EQUAL(out_of_bounds_calls, size_t(0)); + for (auto x : local_minima) { + CHECK_ABSOLUTE_ERROR(Real(1), x, Real(0.01)); + } +} + +#ifndef BOOST_MATH_BUILD_MODULE +// detail::symmetric_eigendecomposition is not exported from the module. +template +void check_symmetric_eigendecomposition(std::vector const & C, size_t n) { + using boost::math::optimization::detail::symmetric_eigendecomposition; + std::vector B; + std::vector D; + CHECK_EQUAL(symmetric_eigendecomposition(C, n, B, D), true); + // The residuals are computed in long double, so that they measure the error of the decomposition: + long double C_norm_sq = 0; + long double residual_sq = 0; + long double orthogonality_sq = 0; + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < n; ++j) { + C_norm_sq += static_cast(C[i*n + j])*C[i*n + j]; + long double CB = -static_cast(B[i*n + j])*D[j]; + long double BTB = (i == j) ? -1 : 0; + for (size_t k = 0; k < n; ++k) { + CB += static_cast(C[i*n + k])*B[k*n + j]; + BTB += static_cast(B[k*n + i])*B[k*n + j]; + } + residual_sq += CB*CB; + orthogonality_sq += BTB*BTB; + } + } + const long double tol = (2*n + 8)*static_cast(std::numeric_limits::epsilon()); + CHECK_LE(std::sqrt(residual_sq/C_norm_sq), tol); + CHECK_LE(std::sqrt(orthogonality_sq), tol); +} + +template +void test_symmetric_eigendecomposition() { + std::cout << "Testing the symmetric eigendecomposition used by CMA-ES . . .\n"; + std::mt19937_64 gen(12345); + std::normal_distribution dis(0, 1); + for (size_t n : {1, 2, 3, 10, 50}) { + // G G^T + nI is symmetric positive definite: + std::vector G(n*n); + for (auto & g : G) { + g = dis(gen); + } + std::vector C(n*n, Real(0)); + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < n; ++j) { + for (size_t k = 0; k < n; ++k) { + C[i*n + j] += G[i*n + k]*G[j*n + k]; + } + } + C[i*n + i] += static_cast(n); + } + check_symmetric_eigendecomposition(C, n); + // A diagonal matrix: + std::vector diagonal(n*n, Real(0)); + for (size_t i = 0; i < n; ++i) { + diagonal[i*n + i] = static_cast(i + 1); + } + check_symmetric_eigendecomposition(diagonal, n); + // I + uu^T has the eigenvalue 1 repeated n - 1 times: + std::vector u(n); + for (auto & ui : u) { + ui = dis(gen); + } + std::vector repeated(n*n); + for (size_t i = 0; i < n; ++i) { + for (size_t j = 0; j < n; ++j) { + repeated[i*n + j] = (i == j ? Real(1) : Real(0)) + u[i]*u[j]; + } + } + check_symmetric_eigendecomposition(repeated, n); + } +} +#endif + int main() { #if (defined(__clang__) || defined(_MSC_VER)) test_ackley(); @@ -156,5 +270,11 @@ int main() { test_beale(); #endif test_sphere(); + test_bounds(); + test_bounds(); +#ifndef BOOST_MATH_BUILD_MODULE + test_symmetric_eigendecomposition(); + test_symmetric_eigendecomposition(); +#endif return boost::math::test::report_errors(); }