NLopt🔗ℹ

I consider this package to be in a somewhat beta state. I don’t yet promise to keep the API from changing. It needs some feedback yet. Feel free to comment on it.

This package provides a wrapper for the NLopt nonlinear optimization package[NLopt], which is a common interface for a number of different optimization routines.

The nlopt module currently reexports the contents of nlopt/highlevel. The interface provided by nlopt/highlevel will be more stable than the interface provided by nlopt. The goal of nlopt will be to facilitate interactive programming, scripting, and generally getting things done quickly with a minimal amount of typing. If nlopt/highlevel is no longer the best way to do that, nlopt will provide different stuff which better satisfies those goals, while nlopt/highlevel will be left alone. (Perhaps nlopt would reexport a hypothetical new "nlopt/highlevel2".)

1 Installation🔗ℹ

This Racket wrapper was developed and tested against NLopt 2.4.2; I’d expect anything later in the 2.4.x series to suffice as well, but there are no guarantees.

You’ll need to get the appropriate NLopt shared library for your Racket installation. Precompiled Windows binaries are available from the NLopt website; make sure you download the appropriate bitness.

Many Linux distributions provide NLopt. Look for libnlopt0 or similar.

In FreeBSD, NLopt is available in the ports collection as nlopt.

And on Mac OS X, I believe you’re stuck compiling it for yourself.

If you have to compile it yourself, and on Windows, where you’re handed a DLL, you’ll need to ensure that the shared library ends up somewhere that Racket will be able to find it with minimal effort. Placing it in the same directory as your code might work, or you can modify your PATH or LD_LIBRARY_PATH environmental variables to point to it.

2 High Level Interface🔗ℹ

This is the most unstable part of the package. Not only will things here change, they might not even work right now.

The highlevel package provides interfaces to NLopt functionality that are more convenient and self-contained than the low-level C-like interfaces. Variants are provided to cover a wide variety of styles that one might use to represent n-ary functions and their corresponding gradient functions.

• /flvector: for mathematical functions that take and produce one n-element flvector;

• /f64vector: corresponding versions for n-element f64vectors;

• /vector: corresponding versions for n-element vectors;

• /list: corresponding versions for n-element lists;

• /arg: operate seamlessly with n-argument Racket functions, where the result is either a single number or a n-element list of numbers.

2.1 A Quick-start Example🔗ℹ

As a brief example of the high-level interface in action, we recreate the example from the NLopt Tutorial (https://nlopt.readthedocs.io/en/latest/NLopt_Tutorial/) here.

Our goal is to solves the following nonlinearly-constrained minimization problem:

consider the function:

f(x0,x1) = (sqrt x1)

find the pair x0,x1 that (approximately) minimizes f over the region

x1 >= (-x0 + 1)^3;

x1 >= (2*x0 + 0)^3;

Next, render the inequality constraints. NLopt represents each constraint as a function C(x), which is interpreted as imposing the inequality C(x) <= 0. Since both constraints have the parametric shape:

(a*x0 + b)^3

Finally we’re ready: run the search!

(define-values (fn-x x)

(minimize/args fn

initial-x

#:bounds bounds

#:ineq-constraints ineq-constraints))

The minimize/args function returns two values: the minimal value fn-x and the point x known to yield that value.

(define (digits n) (real->decimal-string n 3))

(printf "Result: x = ~a; f(x) = ~a.\n"

(map digits x)

(digits fn-x))

2.1.1 Some extra options🔗ℹ

The example above is pretty minimal, but the high-level interface provides more options to tailor the search. For instance, the above code lets the interface choose the optimization algorithm. Furthermore, the interface synthesizes its own approximate gradient (a.k.a. Jacobian) functions for the target function and the constraint functions. We can supply our own more accurate gradient functions, as well as select the algorithm and the optimization direction. The following additions execute that plan.

First, introduce gradient functions for both the target and constraint functions. The constraint gradient functions are parameterized, just like the originals.

(define (grad-fn x0 x1)

(list 0.0 (/ 0.5 (sqrt x1))))

(define ((grad-cn a b) x0 x1)

(list (* 3 a (expt (+ (* a x0) b) 2))

0.0))

Now supply constraint gradient functions alongside the constraint functions

(define ineq-constraint-grads

(list (cons (cn 2.0 0.0) (grad-cn 2.0 0.0))

(cons (cn -1.0 1.0) (grad-cn -1.0 1.0))))

Finally, supply the #:minimize #t option to the general optimize/args function to choose minimization, and provide the target gradient function and constraint-gradient function pairs using keyword options.

(define-values (fn-x^ x^)

(optimize/args fn

initial-x

#:minimize #t

#:jac grad-fn

#:method 'LD_MMA

#:bounds bounds

#:ineq-constraints ineq-constraint-grads))

(printf "Result: x = ~a; f(x) = ~a.\n"

(map digits x^)

(digits fn-x^))

2.2 The Programming Interface🔗ℹ

These /flvector variants require flonum? values. (Which is largely enforced by the flvectors themselves.) /flvector objective functions should be of the form (fun x) where x is an flvector? and the Jacobians should be of the form (jac y x grad) where y is (fun x), and grad is an flvector? to be populated with the gradient.

fun is the procedure to be optimized. It shouldn’t be invoked significantly more than maxeval times, over maxtime seconds. x0 is your initial guess for the optimization; some algorithms are more sensitive to the quality of your initial guess than others. data is passed to every invocation of fun or jac. You may use force-stop inside of fun.

#:maximize and #:minimize determine whether a maximization, or minimzation will be performed. Exactly one of them must be #t. Anything else will result in an exception being raised.

method is a symbol indicating which optimization algorithm to run. See the Algorithms section for your options. If you omit it, or set it to #f, an algorithm will be chosen automatically based on jac, bounds, ineq-constraints and eq-constraints. It should run without error; performance is not guaranteed.

jac is the Jacobian of fun. If you omit it, or supply #f then a very simple approximation will be constructed, by determining how much fun varies when the current x is varied by epsilon. If you provide a Jacobian, or it is not used by the algorithm you select, then epsilon is unused.

bounds may be #f in which case no upper or lower bounds are applied. If it isn’t #f then it should be a sequence of the same length as x0, each element in the sequence should be a pair. The car of each pair will be the lower bound, and the cdr the upper bound. You may supply +max.0 and -max.0 if you don’t wish to bound a dimension above or below, respectively.

ineq-constraints and eq-constraints are sequences of constraint functions (or #f). They must have the same interface as an objective function. An inequality constraint f will constrain the optimization so that (<= (f x _) 0.0) is #t. An equality constraint requires that (= (f x _) 0.0) remain #t. You may provide just the constraint function itself, or a pair, containing the constraint function in car, and the Jacobian of the constraint function in cdr.

tolerance is not currently used. Sorry!

This procedure’s interface was based on scipy’s optimize function.

The /f64vector variants should perform about as well as the /flvector variants. They accept any real? values. /flvector objective functions should be of the form (fun x) where x is an f64vector? and the Jacobians should be of the form (jac y x grad) where y is (fun x), and grad is an f64vector? to be populated with the gradient.

The /vector variants will be less efficient than the /flvector variants. They accept any real? values. /vector objective functions should be of the form (fun x) where x is a vector? and the Jacobians should be of the form (jac y x grad) where y is (fun x), and grad is a vector? to be populated with the gradient.

The /list variants will be the slowest. Like /vector variants, will handle any real? values. /list objective functions should be of the form (fun x) where x is a list? and the Jacobians should be of the form (jac y x) where y is (fun x). They should return the gradient as a list?.

(maximize/args sin

'(0.0)

#:jac cos

#:bounds '((-inf.0 . 0.0)))

3 Safe Interface🔗ℹ

This module is the safe, fully-contracted version of the interface to the C library. For the unsafe, contractless version, see nlopt/unsafe.

3.1 A FullExample🔗ℹ

The following code reconstructs the NLOpt tutorial https://nlopt.readthedocs.io/en/latest/NLopt_Tutorial/ in terms of the safe interface. The same example is presented in terms of the high-level and unsafe Racket interfaces. Consult those sections of this documentation for more detailed explanation.

(define DIMENSIONS 2)

(define (myfunc x grad _)

(define x0 (flvector-ref x 0))

(define x1 (flvector-ref x 1))

(when grad

(flvector-set! grad 0 0.0)

(flvector-set! grad 1 (/ 0.5 (sqrt x1))))

(sqrt x1))

(define ((myconstraint a b) x grad _)

(define x0 (flvector-ref x 0))

(define x1 (flvector-ref x 1))

(when grad

(flvector-set! grad

0

(* 3 a (expt (+ (* a x0) b) 2)))

(flvector-set! grad

1

0.0))

(- (expt (+ (* a x0) b) 3) x1))

(define opt (create 'LD_MMA DIMENSIONS))

(set-lower-bounds opt (flvector -inf.0 0.0))

(set-min-objective opt myfunc #f)

(add-inequality-constraint opt (myconstraint 2.0 0.0) #f 1e-8)

(add-inequality-constraint opt (myconstraint -1.0 0.0) #f 1e-8)

(set-xtol-rel opt 0.0001)

(define x (flvector 1.234 5.678))

(define-values (result minf) (optimize opt x))

(define HARD-FAILURE '(FAILURE INVALID_ARGS OUT_OF_MEMORY))

(when (member result HARD-FAILURE)

(error "nlopt failed: ~a\n" result))

(when (equal? result 'ROUNDOFF_LIMITED)

(printf "warning: roundoff limited!\n"))

(printf "found minimum at f(~a,~a) = ~a\n"

(real->decimal-string (flvector-ref x 0) 3)

(real->decimal-string (flvector-ref x 1) 3)

(real->decimal-string minf 3))

3.2 Basics🔗ℹ

The general pattern for using this library is to create an options structure, apply the various setup options to it (making sure to include a stopping condition!), and then run optimize.

x must be at least as large as the dimension of opt.

3.3 Constraints🔗ℹ

procedure (set-lower-bounds opt bounds) → nlopt-result/c

opt : nlopt-opt?

bounds : (flvector/length? (get-dimension opt))

procedure (set-upper-bounds opt bounds) → nlopt-result/c

opt : nlopt-opt?

bounds : (flvector/length? (get-dimension opt))

procedure (get-lower-bounds opt)

→   [res nlopt-result/c] [bounds (flvector/length? (get-dimension opt))]

opt : nlopt-opt?

procedure (get-upper-bounds opt)

→   [res nlopt-result/c] [bounds (flvector/length? (get-dimension opt))]

opt : nlopt-opt?

3.4 Stopping Criteria🔗ℹ

procedure (set-stopval opt stopval) → nlopt-result/c

opt : nlopt-opt?

stopval : real?

procedure (get-stopval opt) → flonum?

opt : nlopt-opt?

procedure (set-ftol-rel opt ftol-rel) → nlopt-result/c

opt : nlopt-opt?

ftol-rel : real?

procedure (get-ftol-rel opt) → flonum?

opt : nlopt-opt?

procedure (set-ftol-abs opt ftol-abs) → nlopt-result/c

opt : nlopt-opt?

ftol-abs : real?

procedure (get-ftol-abs opt) → flonum?

opt : nlopt-opt?

procedure (set-xtol-rel opt xtol-rel) → nlopt-result/c

opt : nlopt-opt?

xtol-rel : real?

procedure (get-xtol-rel opt) → flonum?

opt : nlopt-opt?

procedure (set-xtol-abs1 opt xtol-abs) → nlopt-result/c

opt : nlopt-opt?

xtol-abs : real?

procedure (set-xtol-abs opt xtols) → nlopt-result/c

opt : nlopt-opt?

xtols : (flvector/length? (get-dimension opt))

procedure (get-xtol-abs opt)   →   nlopt-result/c (flvector/length? (get-dimension opt))

opt : nlopt-opt?

procedure (set-maxeval opt maxeval) → nlopt-result/c

opt : nlopt-opt?

maxeval : natural-number/c

procedure (get-maxeval opt) → natural-number/c

opt : nlopt-opt?

procedure (set-maxtime opt maxtime) → nlopt-result/c

opt : nlopt-opt?

maxtime : real?

procedure (get-maxtime opt) → flonum?

opt : nlopt-opt?

3.5 Algorithm-Specific Parameters🔗ℹ

procedure (set-local-optimizer opt local-opt) → nlopt-result/c

opt : nlopt-opt?

local-opt : nlopt-opt?

4 Unsafe Interface🔗ℹ

This module is the unsafe, contractless version of the interface to the C library. For the safe, fully-contracted version, see nlopt/safe.

The one bit of safety provided by this module is that nlopt_opt structures will be cleaned up properly, and Racket values passed to NLopt procedures will be held onto until NLopt no longer refers to them.

The API of this module is probably the most stable thing in the whole package, as it is a direct mapping of the C API.

Failure to obey any of the requirements mentioned below will probably cause Racket to crash.

4.1 Quickstart Example🔗ℹ

To demonstrate the unsafe interface, the following code reconstructs the NLOpt tutorial https://nlopt.readthedocs.io/en/latest/NLopt_Tutorial/ in terms of it. To understand the following, it is recommended that you familiarize yourself with the problem description and analogous code for the high-level nlopt interface.

In order to properly work with the unsafe interface, we require features of Racket’s unsafe FFI library. At times we also exploit flonum? and flvectors, which are vectors of flonum? objects that have the same memory layout as arrays of double from C.

(require nlopt/unsafe)

(require racket/flonum)

(require ffi/unsafe)

The unsafe interface must occasionally indicate the dimensionality of the target function’s search space.

(define DIMENSIONS 2)

To emphasize the low-level nature of the unsafe interface, the unsafe implementation mirrors the implementation that can be found in the NLOpt tutorial. If you are familiar with C, then you might wish to compare this Racket code to that C code.

(define (myfunc n x grad data)

(define x0 (ptr-ref x _double 0))

(define x1 (ptr-ref x _double 1))

(when grad

(ptr-set! grad _double 0 0.0)

(ptr-set! grad _double 1 (/ 0.5 (sqrt x1))))

(sqrt x1))

This function uses Racket’s unsafe memory operations ptr-ref and ptr-set! to manipulate arrays of double values. Racket’s FFI translates null C pointers to #f, so a simple Boolean check on grad is used to determine whether NLopt needs the function to compute its gradient.

Following the NLOpt tutorial, we use the NLOpt’s client data mechanism to register different values for a and b to be used in two separate constraints. Note that this approach is not recommended: client data is used in C to simulate closures, but Racket’s FFI supports passing closures to C as callbacks (see the highlevel API’s quick-start example). We use client data here simply to illustrate what is possible.

Since NLOpt retains pointers to client data, we use Racket’s custom memory allocation to construct objects to create two-element arrays that the garbage collector guarantees not to move (i.e., interior memory), and are known not to contain pointers (i.e., atomic memory). See the malloc documentation for details. Later, we show how to pass a Racket struct through this interface, using a specially-allocated pointer to the struct.

(define (make-constraint-data a b)

(define ptr (malloc _double 2 'atomic-interior))

(ptr-set! ptr _double 0 a)

(ptr-set! ptr _double 1 b)

ptr)

(define (constraint-data-a cd)

(ptr-ref cd _double 0))

(define (constraint-data-b cd)

(ptr-ref cd _double 1))

The myconstraint function implements a family of inequality constraints, parameterized on values a and b, which are retrieved from client-data.

(define (myconstraint n x grad data)

(define x0 (ptr-ref x _double 0))

(define x1 (ptr-ref x _double 1))

(define a (constraint-data-a data))

(define b (constraint-data-b data))

(when grad

(ptr-set! grad _double 0

(* 3 a (expt (+ (* a x0) b) 2)))

(ptr-set! grad _double 1 0.0))

(- (expt (+ (* a x0) b) 3) x1))

To prepare for optimization, we create an optimization options object that commits to its algorithm—a Local, graDient-exploiting variant of MMA (see NLopt documentation for details and other algorithms).

(define opt (create 'LD_MMA DIMENSIONS))

Next, we set lower bounds on the search space.

(define lower-bounds (flvector -inf.0 0.0))

(set-lower-bounds opt (flvector->cpointer lower-bounds))

Next, we configure the optimization to minimize the myfunc function. Client data is initialized to #f (equivalent to a null pointer in C), since myfunc uses no client data.

(set-min-objective opt myfunc #f)

We configure two constraints using a combination of the myconstraint function and distinct client data for each constraint.

(define data1 (make-constraint-data 2.0 0.0))

(define data2 (make-constraint-data -1.0 1.0))

(add-inequality-constraint opt myconstraint  data1 1e-8)

(add-inequality-constraint opt myconstraint  data2 1e-8)

NLOpt needs at least one stopping criterium to ensure that optimization terminates. Here we set a lower-bound on the the relative change of x during search. See the documentation for alternative termination criteria.

(set-xtol-rel opt 0.0001)

Finally, we allocate an interior array that is set to an initial search point, and that NLOpt will modify to hold the final search value.

(define x (malloc _double DIMENSIONS 'atomic-interior))

(ptr-set! x _double 0 1.234)

(ptr-set! x _double 1 5.678)

With all settings in place, we initiate the optimization.

(define-values (result minf) (optimize opt x))

optimize returns three values. result is a status flag, minf is the minimum value found during search, and x is modified to hold the point at which the function produced minf.

We can interrogate the status flag and possibly report errors.

(define HARD-FAILURE '(FAILURE INVALID_ARGS OUT_OF_MEMORY))

(when (member result HARD-FAILURE)

(error "nlopt failed: ~a\n" result))

According to the NLOpt documentation, the result of optimization may still be useful even if the optimization reports a 'ROUNDOFF_LIMITED result. Nonetheless it is good to report that this happened.

(when (equal? result 'ROUNDOFF_LIMITED)

(printf "warning: roundoff limited!\n"))

Finally, if optimization was successful (which it should be in this example), we can inspect and manipulate the results.

(printf "found minimum at f(~a,~a) = ~a\n"

(real->decimal-string (ptr-ref x _double 0) 3)

(real->decimal-string (ptr-ref x _double 1) 3)

(real->decimal-string minf 3))

4.1.1 Passing Racket Structured Data Directly🔗ℹ

Though closures are the preferred technique for associating client data to an NLOpt callback, one may wish to directly pass structured client data. This can be done using an extra indirection. The following code, without commentary, presents changes to the above code that create a safe memory block that points to a traditional Racket struct, from which a callback can retrieve it.

Note that the 'interior option to malloc only works for versions of Racket after Version 7.1 due to a bug in the memory allocator.

(define (cbox s)

(define ptr (malloc racket 'interior))

(ptr-set! ptr racket s)

ptr)

(define (cunbox cb)

(ptr-ref cb racket))

(define-struct constraint-data (a b))

(define (myconstraint n x grad data)

(define x0 (ptr-ref x _double 0))

(define x1 (ptr-ref x _double 1))

(define cd (cunbox data))

(define a (constraint-data-a cd))

(define b (constraint-data-b cd))

(when grad

(ptr-set! grad _double 0

(* 3 a (expt (+ (* a x0) b) 2)))

(ptr-set! grad _double 1 0.0))

(- (expt (+ (* a x0) b) 3) x1))

(define cbdata1 (cbox (make-constraint-data 2.0 0.0)))

(define cbdata2 (cbox (make-constraint-data -1.0 1.0)))

(add-inequality-constraint opt myconstraint cbdata1 1e-8)

(add-inequality-constraint opt myconstraint cbdata2 1e-8)

4.2 Differences in the Basics🔗ℹ

You should ensure that the x pointer will not be moved over the course of execution. For short runs, or one-off hacky bits of code it probably doesn’t matter. But if you start running long optimizations, sooner or later the garbage collector will move anything that can be moved. malloc with a mode of 'atomic-interior is suggested.

4.3 Differences in the Constraints🔗ℹ