47.1 Introduction to cobyla

fmin_cobyla is a Common Lisp translation (via f2cl) of the Fortran constrained optimization routine COBYLA by Powell[1][2][3].

COBYLA minimizes an objective function F(X) subject to M inequality constraints of the form \(g(X) \ge 0\) on X, where X is a vector of variables that has N components.

Equality constraints g(X) = 0 can often be implemented by a pair of inequality constraints \(g(X) \ge 0\) and \(-g(X) \ge 0.\) Maxima’s interface to COBYLA allows equality constraints and internally converts the equality constraints to a pair of inequality constraints.

The algorithm employs linear approximations to the objective and constraint functions, the approximations being formed by linear interpolation at N+1 points in the space of the variables. The interpolation points are regarded as vertices of a simplex. The parameter RHO controls the size of the simplex and it is reduced automatically from RHOBEG to RHOEND. For each RHO the subroutine tries to achieve a good vector of variables for the current size, and then RHO is reduced until the value RHOEND is reached. Therefore, RHOBEG and RHOEND should be set to reasonable initial changes to and the required accuracy in the variables respectively, but this accuracy should be viewed as a subject for experimentation because it is not guaranteed. The routine treats each constraint individually when calculating a change to the variables, rather than lumping the constraints together into a single penalty function. The name of the subroutine is derived from the phrase Constrained Optimization BY Linear Approximations.

References:

[1] Fortran Code is from http://plato.asu.edu/sub/nlores.html#general

[2] M. J. D. Powell, "A direct search optimization method that models the objective and constraint functions by linear interpolation," in Advances in Optimization and Numerical Analysis, eds. S. Gomez and J.-P. Hennart (Kluwer Academic: Dordrecht, 1994), p. 51-67.

[3] M. J. D. Powell, "Direct search algorithms for optimization calculations," Acta Numerica 7, 287-336 (1998). Also available as University of Cambridge, Department of Applied Mathematics and Theoretical Physics, Numerical Analysis Group, Report NA1998/04 from https://web.archive.org/web/20160607190705/http://www.damtp.cam.ac.uk:80/user/na/reports.html

47.2 Functions and Variables for cobyla

Function: fmin_cobyla
    fmin_cobyla (F, X, Y)
    fmin_cobyla (F, X, Y, optional_args) ¶

Returns an approximate minimum of the expression F with respect to the variables X, subject to an optional set of constraints. Y is a list of initial guesses for X.

F must be ordinary expressions, not names of functions or lambda expressions.

optional_args represents additional arguments, specified as symbol = value. The optional arguments recognized are:

constraints

List of inequality and equality constraints that must be satisfied by X. The inequality constraints must be actual inequalities of the form \(g(X) \ge h(X)\) or \(g(X) \le h(X).\) The equality constraints must be of the form \(g(X) = h(X).\)

rhobeg

Initial value of the internal RHO variable which controls the size of simplex. (Defaults to 1.0)

rhoend

The desired final value rho parameter. It is approximately the accuracy in the variables. (Defaults to 1d-6.)

iprint

Verbose output level. (Defaults to 0)

• 0 - No output
• 1 - Summary at the end of the calculation
• 2 - Each new value of RHO and SIGMA is printed, including the vector of variables, some function information when RHO is reduced.
• 3 - Like 2, but information is printed when F(X) is computed.

maxfun

The maximum number of function evaluations. (Defaults to 1000).

On return, a vector is given:

1. The value of the variables giving the minimum. This is a list of elements of the form var = value for each of the variables listed in X.
2. The minimized function value
3. The number of function evaluations.
4. Return code with the following meanings
   • 0 - No errors.
   • 1 - Limit on maximum number of function evaluations reached.
   • 2 - Rounding errors inhibiting progress.
   • -1 - MAXCV value exceeds RHOEND. This indicates that the constraints were probably not satisfied. User should investigate the value of the constraints.

MAXCV stands for “MAXimum Constraint Violation” and is the value of max(0.0, -c1(x), -c2(x),...-cm(x)) where ck(x) denotes the k’th constraint function. (Note that maxima allows constraints of the form f(x) = g(x), which are internally converted to f(x)-g(x) >= 0 and g(x)-f(x) >= 0 which is required by COBYLA).

load("fmin_cobyla") loads this function.

Function: bf_fmin_cobyla
    bf_fmin_cobyla (F, X, Y)
    bf_fmin_cobyla (F, X, Y, optional_args) ¶

This function is identical to fmin_cobyla, except that bigfloat operations are used, and the default value for rhoend is 10^(fpprec/2).

See fmin_cobyla for more information.

load("bf_fmin_cobyla") loads this function.

47.3 Examples for cobyla

Minimize \(x_1 x_2\) with \(1-x_1^2-x_2^2 \ge 0.\) The theoretical solution is $$ \eqalign{ x_1 &= {1\over \sqrt{2}} \cr x_2 &= -{1\over \sqrt{2}} } $$

(%i1) load("fmin_cobyla")$

(%i2) fmin_cobyla(x1*x2, [x1, x2], [1,1], 
                  constraints = [x1^2+x2^2<=1], iprint=1);

   Normal return from subroutine COBYLA

   NFVALS =   66   F =-5.000000E-01    MAXCV = 1.999845E-12
   X = 7.071058E-01  -7.071077E-01
(%o2) [[x1 = 0.70710584934848, x2 = - 0.7071077130248], 
       - 0.49999999999926, [[-1.999955756559757e-12],[]], 66]

Here is the same example but the constraint is \(x_1^2+x_2^2 \le -1\) which is impossible over the reals.

(%i1) fmin_cobyla(x1*x2, [x1, x2], [1,1],
         constraints = [x1^2+x2^2 <= -1], iprint=1);

   Normal return from subroutine COBYLA

   NFVALS =   65   F = 3.016417E-13    MAXCV = 1.000000E+00
   X =-3.375179E-07  -8.937057E-07
(%o1) [[x1 = - 3.375178983064622e-7, x2 = - 8.937056510780022e-7], 
                                                3.016416530564557e-13, 65, - 1]
(%i2) subst(%o1[2], [x1^2+x2^2 <= -1]);
(%o2)                 [- 6.847914590915444e-13 <= - 1]

We see the return code (%o1[4]) is -1 indicating that the constraints may not be satisfied. Substituting the solution into the constraint equation as shown in %o2 shows that the constraint is, of course, violated.

There are additional examples in the share/cobyla/ex directory and in share/cobyla/rtest_cobyla.mac.