RFC: Per-variable relative tolerance criteria for x (#183)

* xtol_rel: L1 norm for whole-vector stopping criterion

* vector_norm: fix uninitialized ret, move branch outside for loop, make it const-correct

* diff_norm: like vector_norm, but for x-oldx

* xtol_abs: test directly when diff_norm tests for weighted xtol_rel

* add double * x_weights to nlopt and nlopt_stopping structures

 - use the weights when checking relative x stopping criterion
 - allow setting the weights like it's implemented for xtol_abs
 - weights default to 1

* don't forget to free the weights; actually declare the weight-related functions

* provide scale support for vector_norm and diff_norm

* add argument checks to nlopt_get_x_weights

* do not allocate x_weights unless requested by

Unallocated (NULL) x_weights behave as if all weigths are set to 1.

* start documenting the new functions

TODO: add *_x_weights functions to C++, Fortran, Python interface,
document those

* fixes, add equations

* Update NLopt_Reference.md

* param names

* check for w < 0 in set_x_weights

* add x_weights to C++ interface and document them

* add x_weights to F77 interface and document them

* add & document x_weights to SWIG-based interfaces

This includes Guile and Python.

* add and document x_weights to Octave&Matlab

* set_x_weights: ISO C90 compatibility

* set_x_weights: provide informative error on w <= 0

* tweaks

doc/docs/NLopt_C-plus-plus_Reference.md

@@ -203,6 +203,16 @@ double nlopt::opt::get_xtol_rel() const;
 
 Set relative tolerance on optimization parameters.
 
+```
+void nlopt::opt::set_x_weights(const std::vector`<double>` &w);
+void nlopt::opt::set_x_weights(double w);
+void nlopt::opt::get_x_weights(std::vector`<double>` &w) const;
+std::vector`<double>` nlopt::opt::get_x_weights() const;
+```
+
+
+Set/get the weights used when the computing L₁ norm for the `xtol_rel` stopping criterion above.
+
 ```
 void nlopt::opt::set_xtol_abs(const std::vector`<double>` &tol);
 void nlopt::opt::set_xtol_abs(double tol);

doc/docs/NLopt_Fortran_Reference.md

@@ -241,6 +241,15 @@ call nlo_get_xtol_rel(tol, opt)
 
 Set relative tolerance on optimization parameters.
 
+```
+call nlo_set_x_weights(ires, opt, w)
+call nlo_set_x_weights1(ires, opt, w1)
+call nlo_get_x_weights(ires, opt, w)
+```
+
+
+Set/get the weights used when the computing L₁ norm for the `xtol_rel` stopping criterion above.
+
 ```
 call nlo_set_xtol_abs(ires, opt, tol)
 call nlo_set_xtol_abs1(ires, opt, tol1)

doc/docs/NLopt_Guile_Reference.md

@@ -178,6 +178,14 @@ Set relative tolerance on optimization parameters.
 
 Set absolute tolerances on optimization parameters. The `tol` input must be a vector or list of length `n` (the dimension specified in the `nlopt.opt` constructor); alternatively, you can pass a single number in order to set the same tolerance for all optimization parameters. `get-xtol-abs()` returns the tolerances as a vector.
 
+```
+(nlopt-opt-set-x-weights opt x)
+(nlopt-opt-get-x-weights opt x)
+```
+
+
+Set the weights used when the computing L₁ norm for the `xtol_rel` stopping criterion above.
+
 ```
 (nlopt-opt-set-maxeval opt maxeval)
 (nlopt-opt-get-maxeval opt)

doc/docs/NLopt_Matlab_Reference.md

@@ -101,6 +101,13 @@ opt.xtol_rel
 
 Set relative tolerance on optimization parameters.
 
+```
+opt.x_weights
+```
+
+
+Set the weights used when the computing L₁ norm for the `xtol_rel` stopping criterion above. The `opt.x_weights` value must be a vector of length `n` (the same length as the initial guess passed to `nlopt_optimize`).
+
 ```
 opt.xtol_abs
 ```

doc/docs/NLopt_Python_Reference.md

@@ -223,6 +223,14 @@ opt.get_xtol_abs()
 
 Set absolute tolerances on optimization parameters. The `tol` input must be an array (NumPy array or Python list) of length `n` (the dimension specified in the `nlopt.opt` constructor); alternatively, you can pass a single number in order to set the same tolerance for all optimization parameters. `get_xtol_abs()` returns the tolerances as a NumPy array.
 
+```
+opt.set_x_weights(w)
+opt.get_x_weights()
+```
+
+
+Set the weights used when the computing L₁ norm for the `xtol_rel` stopping criterion above.
+
 ```
 opt.set_maxeval(maxeval)
 opt.get_maxeval()

doc/docs/NLopt_Reference.md

@@ -247,25 +247,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/f77funcs_.h

@@ -152,6 +152,7 @@ F77_GETSET(stopval, STOPVAL, double)
 F77_GETSET(ftol_rel, FTOL_REL, double)
 F77_GETSET(ftol_abs, FTOL_ABS, double)
 F77_GETSET(xtol_rel, XTOL_REL, double) F77_GETSETA(xtol_abs, XTOL_ABS, double) F77_GETSET(maxeval, MAXEVAL, int) F77_GET(numevals, NUMEVALS, int) F77_GETSET(maxtime, MAXTIME, double)
+F77_GETSETA(x_weights, X_WEIGHTS, double)
  F77_GETSET(force_stop, FORCE_STOP, int)
 NLOPT_EXTERN(void) F77_(nlo_force_stop, NLO_FORCE_STOP) (int *ret, nlopt_opt * opt)
 {

src/api/nlopt-in.hpp

@@ -470,6 +470,7 @@ namespace nlopt {
     NLOPT_GETSET(double, ftol_abs)
     NLOPT_GETSET(double, xtol_rel)
     NLOPT_GETSET_VEC(xtol_abs)
+    NLOPT_GETSET_VEC(x_weights)
     NLOPT_GETSET(int, maxeval)
 
     int get_numevals() const {

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

@@ -253,6 +253,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 w);
+NLOPT_EXTERN(nlopt_result) nlopt_set_x_weights(nlopt_opt opt, const double *w);
+NLOPT_EXTERN(nlopt_result) nlopt_get_x_weights(const nlopt_opt opt, double *w);
 
 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));
@@ -632,6 +640,58 @@ nlopt_result NLOPT_STDCALL nlopt_get_xtol_abs(const nlopt_opt opt, double *xtol_
     return NLOPT_INVALID_ARGS;
 }
 
+nlopt_result NLOPT_STDCALL nlopt_set_x_weights(nlopt_opt opt, const double *x_weights)
+{
+    if (opt) {
+        unsigned i;
+        nlopt_unset_errmsg(opt);
+        for (i = 0; i < opt->n; i++)
+            if (x_weights[i] < 0)
+                return ERR(NLOPT_INVALID_ARGS, opt, "invalid negative weight");
+        if (!opt->x_weights && opt->n > 0) {
+            opt->x_weights = (double *) calloc(opt->n, sizeof(double));
+            if (!opt->x_weights) return NLOPT_OUT_OF_MEMORY;
+        }
+        if (opt->n > 0) 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_weight)
+{
+    if (opt) {
+        unsigned i;
+        if (x_weight < 0) return ERR(NLOPT_INVALID_ARGS, opt, "invalid negative weight");
+        nlopt_unset_errmsg(opt);
+        if (!opt->x_weights && opt->n > 0) {
+            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_weight;
+        return NLOPT_SUCCESS;
+    }
+    return NLOPT_INVALID_ARGS;
+}
+
+nlopt_result NLOPT_STDCALL nlopt_get_x_weights(const nlopt_opt opt, double *x_weights)
+{
+    if (opt) {
+	if (opt->n > 0 && !x_weights) return ERR(NLOPT_INVALID_ARGS, opt, "invalid NULL weights");
+        nlopt_unset_errmsg(opt);
+        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/octave/nlopt_optimize-mex.c

@@ -171,6 +171,8 @@ nlopt_opt make_opt(const mxArray *opts, unsigned n)
      nlopt_set_xtol_rel(opt, struct_val_default(opts, "xtol_rel", 0.0));
      nlopt_set_xtol_abs(opt, struct_arrval(opts, "xtol_abs", n,
 					   fill(tmp, n, 0.0)));
+     nlopt_set_x_weights(opt, struct_arrval(opts, "x_weights", n,
+					   fill(tmp, n, 1.0)));
      nlopt_set_maxeval(opt, struct_val_default(opts, "maxeval", 0.0) < 0 ?
 		       0 : struct_val_default(opts, "maxeval", 0.0));
      nlopt_set_maxtime(opt, struct_val_default(opts, "maxtime", 0.0));

src/octave/nlopt_optimize-oct.cc

@@ -200,6 +200,12 @@ nlopt_opt make_opt(octave_map &opts, int n)
     CHECK1(n == xtol_abs.numel(), "stop.xtol_abs must have same length as x");
     CHECK1(nlopt_set_xtol_abs(opt, xtol_abs.data())>0, "nlopt: out of memory");
   }
+  {
+    Matrix ones(1, n, 1.0);
+    Matrix x_weights = struct_val_default(opts, "x_weights", ones);
+    CHECK1(n == x_weights.numel(), "stop.x_weights must have same length as x");
+    CHECK1(nlopt_set_x_weights(opt, x_weights.data())>0, "nlopt: invalid x_weights or out of memory");
+  }
 
   nlopt_set_maxeval(opt, struct_val_default(opts, "maxeval", 0) < 0 ?
 		    0 : struct_val_default(opts, "maxeval", 0));

src/swig/nlopt-exceptions.i

@@ -53,6 +53,7 @@ GETSET_EXCEPT(ftol_rel, double)
 GETSET_EXCEPT(ftol_abs, double)
 GETSET_EXCEPT(xtol_rel, double)
 GETSETVEC_EXCEPT(xtol_abs)
+GETSETVEC_EXCEPT(x_weights)
 GETSET_EXCEPT(maxeval, int)
 GETSET_EXCEPT(maxtime, double)
 GETSET_EXCEPT(force_stop, int)

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;
 }