RFC: Per-variable relative tolerance criteria for x

doc/docs/NLopt_Reference.md

@@ -239,25 +239,25 @@ nlopt_result nlopt_set_xtol_rel(nlopt_opt opt, double tol);
 double nlopt_get_xtol_rel(const nlopt_opt opt);
 ```
 
-
-Set relative tolerance on optimization parameters: stop when an optimization step (or an estimate of the optimum) changes every parameter by less than `tol` multiplied by the absolute value of the parameter. (If there is any chance that an optimal parameter is close to zero, you might want to set an absolute tolerance with `nlopt_set_xtol_abs` as well.) Criterion is disabled if `tol` is non-positive.
+Set relative tolerance on optimization parameters: stop when an optimization step (or an estimate of the optimum) causes a relative change the parameters $x$ by less than `tol`, i.e. $\Vert \Delta x \Vert_w < \mbox{tol}\cdot\Vert x \Vert_w$ as measured by a weighted L₁ norm $\Vert x \Vert_w = \sum_i w_i |x_i|$, where the weights $w_i$ default to 1.
+(If there is any chance that the optimal $\Vert x \Vert$ is close to zero, you might want to set an absolute tolerance with `nlopt_set_xtol_abs` as well.) Criterion is disabled if `tol` is non-positive.
 
 ```
-nlopt_result nlopt_set_xtol_abs(nlopt_opt opt, const double* tol);
-nlopt_result nlopt_get_xtol_abs(const nlopt_opt opt, double *tol);
+nlopt_result nlopt_set_x_weights(nlopt_opt opt, const double *w);
+nlopt_result nlopt_set_x_weights1(nlopt_opt opt, const double w);
+nlopt_result nlopt_get_x_weights(const nlopt_opt opt, double *w);
 ```
 
-
-Set absolute tolerances on optimization parameters. `tol` is a pointer to an array of length `n` (the dimension from `nlopt_create`) giving the tolerances: stop when an optimization step (or an estimate of the optimum) changes every parameter `x[i]` by less than `tol[i]`. (Note that this function makes a copy of the `tol` array, so subsequent changes to the caller's `tol` have no effect on `opt`.) In `nlopt_get_xtol_abs`, `tol` must be an array of length `n`, which upon successful return contains a copy of the current tolerances.
-
-For convenience, the following function may be used to set the absolute tolerances in all `n` optimization parameters to the same value:
+Set/get the weights used when the computing L₁ norm for the `xtol_rel` stopping criterion above, where `*w` must point to an array of length equal to the number of optimization parameters in `opt`.   `nlopt_set_x_weights1` can be used to set all of the weights to the same value `w`.   The weights default to `1`, but non-constant weights can be used to handle situations where the different parameters `x` have different units or importance, for example.
 
 ```
+nlopt_result nlopt_set_xtol_abs(nlopt_opt opt, const double *tol);
 nlopt_result nlopt_set_xtol_abs1(nlopt_opt opt, double tol);
+nlopt_result nlopt_get_xtol_abs(const nlopt_opt opt, double *tol);
 ```
 
 
-Criterion is disabled if `tol` is non-positive.
+Set absolute tolerances on optimization parameters. `tol` is a pointer to an array of length `n` (the dimension from `nlopt_create`) giving the tolerances: stop when an optimization step (or an estimate of the optimum) changes every parameter `x[i]` by less than `tol[i]`. (Note that `nlopt_set_xtol_abs` makes a copy of the `tol` array, so subsequent changes to the caller's `tol` have no effect on `opt`.) In `nlopt_get_xtol_abs`, `tol` must be an array of length `n`, which upon successful return contains a copy of the current tolerances.  For convenience, the `nlopt_set_xtol_abs1` may be used to set the absolute tolerances in all `n` optimization parameters to the same value.  Criterion is disabled if `tol` is non-positive.
 
 ```
 nlopt_result nlopt_set_maxeval(nlopt_opt opt, int maxeval);

src/api/nlopt-internal.h

@@ -57,6 +57,7 @@ extern "C" {
         double stopval;         /* stop when f reaches stopval or better */
         double ftol_rel, ftol_abs;      /* relative/absolute f tolerances */
         double xtol_rel, *xtol_abs;     /* rel/abs x tolerances */
+        double *x_weights;      /* weights for relative x tolerance */
         int maxeval;            /* max # evaluations */
         int numevals;           /* number of evaluations */
         double maxtime;         /* max time (seconds) */

src/api/nlopt.h

@@ -243,6 +243,9 @@ NLOPT_EXTERN(double) nlopt_get_xtol_rel(const nlopt_opt opt);
 NLOPT_EXTERN(nlopt_result) nlopt_set_xtol_abs1(nlopt_opt opt, double tol);
 NLOPT_EXTERN(nlopt_result) nlopt_set_xtol_abs(nlopt_opt opt, const double *tol);
 NLOPT_EXTERN(nlopt_result) nlopt_get_xtol_abs(const nlopt_opt opt, double *tol);
+NLOPT_EXTERN(nlopt_result) nlopt_set_x_weights1(nlopt_opt opt, double tol);
+NLOPT_EXTERN(nlopt_result) nlopt_set_x_weights(nlopt_opt opt, const double *tol);
+NLOPT_EXTERN(nlopt_result) nlopt_get_x_weights(const nlopt_opt opt, double *tol);
 
 NLOPT_EXTERN(nlopt_result) nlopt_set_maxeval(nlopt_opt opt, int maxeval);
 NLOPT_EXTERN(int) nlopt_get_maxeval(const nlopt_opt opt);

src/api/optimize.c

@@ -426,6 +426,7 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf)
     stop.ftol_abs = opt->ftol_abs;
     stop.xtol_rel = opt->xtol_rel;
     stop.xtol_abs = opt->xtol_abs;
+    stop.x_weights = opt->x_weights;
     opt->numevals = 0;
     stop.nevals_p = &(opt->numevals);
     stop.maxeval = opt->maxeval;

src/api/options.c

@@ -52,6 +52,7 @@ void NLOPT_STDCALL nlopt_destroy(nlopt_opt opt)
         free(opt->lb);
         free(opt->ub);
         free(opt->xtol_abs);
+        free(opt->x_weights);
         free(opt->fc);
         free(opt->h);
         nlopt_destroy(opt->local_opt);
@@ -88,6 +89,7 @@ nlopt_opt NLOPT_STDCALL nlopt_create(nlopt_algorithm algorithm, unsigned n)
         opt->stopval = -HUGE_VAL;
         opt->ftol_rel = opt->ftol_abs = 0;
         opt->xtol_rel = 0;
+        opt->x_weights = NULL;
         opt->xtol_abs = NULL;
         opt->maxeval = 0;
         opt->numevals = 0;
@@ -133,7 +135,7 @@ nlopt_opt NLOPT_STDCALL nlopt_copy(const nlopt_opt opt)
         nlopt_munge munge;
         nopt = (nlopt_opt) malloc(sizeof(struct nlopt_opt_s));
         *nopt = *opt;
-        nopt->lb = nopt->ub = nopt->xtol_abs = NULL;
+        nopt->lb = nopt->ub = nopt->xtol_abs = nopt->x_weights = NULL;
         nopt->fc = nopt->h = NULL;
         nopt->m_alloc = nopt->p_alloc = 0;
         nopt->local_opt = NULL;
@@ -157,6 +159,12 @@ nlopt_opt NLOPT_STDCALL nlopt_copy(const nlopt_opt opt)
             nopt->xtol_abs = (double *) malloc(sizeof(double) * (opt->n));
             if (!opt->xtol_abs)
                 goto oom;
+            if (opt->x_weights) {
+                nopt->x_weights = (double *) malloc(sizeof(double) * (opt->n));
+                if (!opt->x_weights)
+                    goto oom;
+                memcpy(nopt->x_weights, opt->x_weights, sizeof(double) * (opt->n));
+            }
 
             memcpy(nopt->lb, opt->lb, sizeof(double) * (opt->n));
             memcpy(nopt->ub, opt->ub, sizeof(double) * (opt->n));
@@ -628,6 +636,56 @@ nlopt_result NLOPT_STDCALL nlopt_get_xtol_abs(const nlopt_opt opt, double *xtol_
     return NLOPT_SUCCESS;
 }
 
+nlopt_result NLOPT_STDCALL nlopt_set_x_weights(nlopt_opt opt, const double *x_weights)
+{
+    if (opt) {
+        nlopt_unset_errmsg(opt);
+        if (!opt->x_weights) {
+            opt->x_weights = (double *) calloc(opt->n, sizeof(double));
+            if (!opt->x_weights) {
+                return NLOPT_OUT_OF_MEMORY;
+            }
+        }
+        memcpy(opt->x_weights, x_weights, opt->n * sizeof(double));
+        return NLOPT_SUCCESS;
+    }
+    return NLOPT_INVALID_ARGS;
+}
+
+nlopt_result NLOPT_STDCALL nlopt_set_x_weights1(nlopt_opt opt, double x_weights)
+{
+    if (opt) {
+        unsigned i;
+        nlopt_unset_errmsg(opt);
+        if (!opt->x_weights) {
+            opt->x_weights = (double *) calloc(opt->n, sizeof(double));
+            if (!opt->x_weights) {
+                return NLOPT_OUT_OF_MEMORY;
+            }
+        }
+        for (i = 0; i < opt->n; ++i)
+            opt->x_weights[i] = x_weights;
+        return NLOPT_SUCCESS;
+    }
+    return NLOPT_INVALID_ARGS;
+}
+
+nlopt_result NLOPT_STDCALL nlopt_get_x_weights(const nlopt_opt opt, double *x_weights)
+{
+    nlopt_unset_errmsg(opt);
+    if (opt && (opt->n == 0 || x_weights)) {
+        if (opt->x_weights) {
+            memcpy(x_weights, opt->x_weights, sizeof(double) * (opt->n));
+        } else {
+            unsigned i;
+            for (i = 0; i < opt->n; ++i)
+                x_weights[i] = 1;
+        }
+        return NLOPT_SUCCESS;
+    }
+    return NLOPT_INVALID_ARGS;
+}
+
 GETSET(maxeval, int, maxeval)
 
     GET(numevals, int, numevals)

src/util/nlopt-util.h

@@ -81,6 +81,7 @@ extern "C" {
         double ftol_abs;
         double xtol_rel;
         const double *xtol_abs;
+        const double *x_weights;
         int *nevals_p, maxeval;
         double maxtime, start;
         int *force_stop;

src/util/stop.c

@@ -29,6 +29,55 @@
 
 /* utility routines to implement the various stopping criteria */
 
+static double sc(double x, double smin, double smax)
+{
+    return smin + x * (smax - smin);
+}
+
+static double vector_norm(unsigned n, const double *vec, const double *w, const double *scale_min, const double *scale_max)
+{
+    unsigned i;
+    double ret = 0;
+    if (scale_min && scale_max) {
+        if (w)
+            for (i = 0; i < n; i++)
+                ret += w[i] * fabs(sc(vec[i], scale_min[i], scale_max[i]));
+        else
+            for (i = 0; i < n; i++)
+                ret += fabs(sc(vec[i], scale_min[i], scale_max[i]));
+    } else {
+        if (w)
+            for (i = 0; i < n; i++)
+                ret += w[i] * fabs(vec[i]);
+        else
+            for (i = 0; i < n; i++)
+                ret += fabs(vec[i]);
+    }
+    return ret;
+}
+
+static double diff_norm(unsigned n, const double *x, const double *oldx, const double *w, const double *scale_min, const double *scale_max)
+{
+    unsigned i;
+    double ret = 0;
+    if (scale_min && scale_max) {
+        if (w)
+            for (i = 0; i < n; i++)
+                ret += w[i] * fabs(sc(x[i], scale_min[i], scale_max[i]) - sc(oldx[i], scale_min[i], scale_max[i]));
+        else
+            for (i = 0; i < n; i++)
+                ret += fabs(sc(x[i], scale_min[i], scale_max[i]) - sc(oldx[i], scale_min[i], scale_max[i]));
+    } else {
+        if (w)
+            for (i = 0; i < n; i++)
+                ret += w[i] * fabs(x[i] - oldx[i]);
+        else
+            for (i = 0; i < n; i++)
+                ret += fabs(x[i] - oldx[i]);
+    }
+    return ret;
+}
+
 static int relstop(double vold, double vnew, double reltol, double abstol)
 {
     if (nlopt_isinf(vold))
@@ -49,33 +98,34 @@ int nlopt_stop_f(const nlopt_stopping * s, double f, double oldf)
 int nlopt_stop_x(const nlopt_stopping * s, const double *x, const double *oldx)
 {
     unsigned i;
+    if (diff_norm(s->n, x, oldx, s->x_weights, NULL, NULL) <= s->xtol_rel * vector_norm(s->n, x, s->x_weights, NULL, NULL))
+        return 1;
     for (i = 0; i < s->n; ++i)
-        if (!relstop(oldx[i], x[i], s->xtol_rel, s->xtol_abs[i]))
+        if (fabs(x[i] - oldx[i]) > s->xtol_abs[i])
             return 0;
     return 1;
 }
 
 int nlopt_stop_dx(const nlopt_stopping * s, const double *x, const double *dx)
 {
     unsigned i;
+    if (vector_norm(s->n, dx, s->x_weights, NULL, NULL) <= s->xtol_rel * vector_norm(s->n, x, s->x_weights, NULL, NULL))
+        return 1;
     for (i = 0; i < s->n; ++i)
-        if (!relstop(x[i] - dx[i], x[i], s->xtol_rel, s->xtol_abs[i]))
+        if (fabs(dx[i]) > s->xtol_abs[i])
             return 0;
     return 1;
 }
 
-static double sc(double x, double smin, double smax)
-{
-    return smin + x * (smax - smin);
-}
-
 /* some of the algorithms rescale x to a unit hypercube, so we need to
    scale back before we can compare to the tolerances */
 int nlopt_stop_xs(const nlopt_stopping * s, const double *xs, const double *oldxs, const double *scale_min, const double *scale_max)
 {
     unsigned i;
+    if (diff_norm(s->n, xs, oldxs, s->x_weights, scale_min, scale_max) <= s->xtol_rel * vector_norm(s->n, xs, s->x_weights, scale_min, scale_max))
+        return 1;
     for (i = 0; i < s->n; ++i)
-        if (!relstop(sc(oldxs[i], scale_min[i], scale_max[i]), sc(xs[i], scale_min[i], scale_max[i]), s->xtol_rel, s->xtol_abs[i]))
+        if (fabs(sc(xs[i], scale_min[i], scale_max[i]) - sc(oldxs[i], scale_min[i], scale_max[i])) > s->xtol_abs[i])
             return 0;
     return 1;
 }