optimize: transform box constraints

Closes stevengj#528

NEWS.md

@@ -1,5 +1,9 @@
 # NLopt Release Notes
 
+in-progress
+
+* Fixed PRAXIS box constraints ([#528])
+
 ## NLopt 2.9
 
 10 November 2024
@@ -517,6 +521,7 @@
 [#504]: https://github.com/stevengj/nlopt/issues/504
 [#509]: https://github.com/stevengj/nlopt/issues/509
 [#510]: https://github.com/stevengj/nlopt/issues/510
+[#528]: https://github.com/stevengj/nlopt/issues/528
 [#533]: https://github.com/stevengj/nlopt/issues/533
 [#534]: https://github.com/stevengj/nlopt/issues/534
 [#535]: https://github.com/stevengj/nlopt/issues/535

doc/docs/NLopt_Algorithms.md

@@ -258,7 +258,7 @@ The original Fortran code was written by Richard Brent and made available by the
 
 Specified in NLopt as `NLOPT_LN_PRAXIS`
 
-This algorithm was originally designed for unconstrained optimization. In NLopt, bound constraints are "implemented" in PRAXIS by the simple expedient of returning infinity (Inf) when the constraints are violated (this is done automatically—you don't have to do this in your own function). This seems to work, more-or-less, but appears to slow convergence significantly. If you have bound constraints, you are probably better off using COBYLA or BOBYQA.
+This algorithm was originally designed for unconstrained optimization. In NLopt, bound constraints are emulated by variable transformation unlike gradient-free algorithms that natively support them like COBYLA or BOBYQA.
 
 ### Nelder-Mead Simplex
 

src/api/optimize.c

@@ -57,23 +57,68 @@
 #include "slsqp.h"
 
 /*********************************************************************/
+/* to emulate box constraint support with variable transform */
 
-static double f_bound(int n, const double *x, void *data_)
+static void x_bound_inv(int n, double *x, const double *lb, const double *ub)
+{
+    int i;
+    double mid, width, th;
+    if (!lb || !ub)
+        return;
+    for (i = 0; i < n; ++ i)
+    {
+        if (!nlopt_isinf(lb[i]) && !nlopt_isinf(ub[i]))
+        {
+            mid = (lb[i] + ub[i]) * 0.5;
+            width = (ub[i] - lb[i]) * 0.5;
+            th = tanh(x[i]);
+            x[i] = mid + th * width;
+        }
+        else if (!nlopt_isinf(lb[i]))
+            x[i] = sqrt(x[i] - lb[i]);
+        else if (!nlopt_isinf(ub[i]))
+            x[i] = sqrt(ub[i] - x[i]);
+    }
+}
+
+static void x_bound(int n, double *x, const double *lb, const double *ub)
 {
     int i;
+    double mid, width, th;
+    if (!lb || !ub)
+        return;
+    for (i = 0; i < n; ++ i)
+    {
+        if (!nlopt_isinf(lb[i]) && !nlopt_isinf(ub[i]))
+        {
+            mid = (lb[i] + ub[i]) * 0.5;
+            width = (ub[i] - lb[i]) * 0.5;
+            th = tanh(x[i]);
+            x[i] = mid + th * width;
+        }
+        else if (!nlopt_isinf(lb[i]))
+            x[i] = lb[i] + x[i] * x[i];
+        else if (!nlopt_isinf(ub[i]))
+            x[i] = ub[i] - x[i] * x[i];
+    }
+}
+
+static double f_bound(int n, const double *x, void *data_)
+{
     nlopt_opt data = (nlopt_opt) data_;
     double f;
+    double *x_tr = NULL;
 
-    /* some methods do not support bound constraints, but support
-       discontinuous objectives so we can just return Inf for invalid x */
-    for (i = 0; i < n; ++i)
-        if (x[i] < data->lb[i] || x[i] > data->ub[i])
-            return HUGE_VAL;
-
-    f = data->f((unsigned) n, x, NULL, data->f_data);
-    return (nlopt_isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f);
+    x_tr = (double *) malloc(n * sizeof(double));
+    memcpy(x_tr, x, n * sizeof(double));
+    x_bound(n, x_tr, data->lb, data->ub);
+    f = data->f((unsigned) n, x_tr, NULL, data->f_data);
+    free(x_tr);
+    return f;
 }
 
+/*********************************************************************/
+
 static double f_noderiv(int n, const double *x, void *data_)
 {
     nlopt_opt data = (nlopt_opt) data_;
@@ -649,9 +694,19 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf)
     case NLOPT_LN_PRAXIS:
         {
             double step;
+            nlopt_result result;
+
             if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
                 return NLOPT_OUT_OF_MEMORY;
-            return praxis_(nlopt_get_param(opt, "t0_tol", 0.0), DBL_EPSILON, step, ni, x, f_bound, opt, &stop, minf);
+
+            /* transform starting point into unbound space */
+            x_bound_inv(ni, x, opt->lb, opt->ub);
+
+            result = praxis_(nlopt_get_param(opt, "t0_tol", 0.0), DBL_EPSILON, step, ni, x, f_bound, opt, &stop, minf);
+
+            /* transform back optimal point to bounded space */
+            x_bound(ni, x, opt->lb, opt->ub);
+            return result;
         }
 
     case NLOPT_LD_LBFGS:

test/CMakeLists.txt

@@ -15,7 +15,8 @@ macro(NLOPT_add_cpp_test test_name)
 endmacro()
 
 NLOPT_add_cpp_test(t_tutorial 23 24 30 39)
-NLOPT_add_cpp_test(cpp_functor)
+NLOPT_add_cpp_test(cpp_functor 0)
+NLOPT_add_cpp_test(t_fbound 0)
 
 NLOPT_add_cpp_test(t_bounded 0 1 2 3 4 5 6 7 8 18 34 41 42)
 if (NOT NLOPT_CXX)

test/t_fbound.cxx

@@ -0,0 +1,46 @@
+#include <iostream>
+#include <vector>
+#include <string>
+#include <cmath>
+#include <iomanip>
+#include <nlopt.hpp>
+
+double myvfunc(const std::vector<double> &x, std::vector<double> &grad, void *data)
+{
+  (void)data;
+  if (!grad.empty()) {
+    grad[0] = 2*x[0];
+    grad[1] = 2*x[1];
+  }
+  if (std::isnan(x[0]))
+    throw std::invalid_argument("nan");
+  return -x[0]*x[0]-x[1]*x[1];
+}
+
+
+int main() {
+  nlopt::srand(0);
+  for (int repeat = 0; repeat < 100; ++repeat)
+  {
+    std::cout << "repeat="<<repeat<<std::endl;
+    nlopt::opt opt(nlopt::LN_PRAXIS, 2);
+    std::vector<double> lb = {-1.0, -10.0};
+    std::vector<double> ub = {10.0, 1.0};
+    opt.set_lower_bounds(lb);
+    opt.set_upper_bounds(ub);
+    opt.set_min_objective(myvfunc, NULL);
+    opt.set_ftol_rel(1e-6);
+    std::vector<double> x = {0.5, 0.5};
+    double minf;
+    try{
+      opt.optimize(x, minf);
+      std::cout << "found optimum at f(" << x[0] << "," << x[1] << ") = "
+                << std::setprecision(10) << minf <<std::endl;
+    }
+    catch(std::exception &e) {
+      std::cerr << "nlopt failed: " << e.what() << std::endl;
+      return EXIT_FAILURE;
+    }
+  }
+  return EXIT_SUCCESS;
+}