PRIMA

This package is a Julia interface to the PRIMA library, a Reference Implementation for Powell's Methods with Modernization and Amelioration, by Zaikun Zhang who re-implemented and improved algorithms originally by M.J.D. Powell for minimizing a multi-variate objective function possibly under constraints and without derivatives.

Depending on the problem to solve, other Julia package(s) with similar objectives may be of interest:

• BlackBoxOptim is a global optimization package for single- and multi-variate problems that does not require the objective function to be differentiable.

• NOMAD is an interface to the Mesh Adaptive Direct Search algorithm (MADS) designed for difficult black-box optimization problems.

Formally, the algorithms provided by PRIMA are designed to solve problems of the form:

min f(x)    subject to   x ∈ Ω ⊆ ℝⁿ

where f: Ω → ℝ is the function to minimize, Ω ⊆ ℝⁿ is the set of feasible variables, and n ≥ 1 is the number of variables. The most general feasible set is:

Ω = { x ∈ ℝⁿ | xl ≤ x ≤ xu, Aₑ⋅x = bₑ, Aᵢ⋅x ≤ bᵢ, cₑ(x) = 0, and cᵢ(x) ≤ 0 }

where xl ∈ ℝⁿ and xu ∈ ℝⁿ are lower and upper bounds, Aₑ and bₑ implement linear equality constraints, Aᵢ and bᵢ implement linear inequality constraints, cₑ: ℝⁿ → ℝᵐ implements m non-linear equality constraints, and cᵢ: ℝⁿ → ℝʳ implements r non-linear inequality constraints.

The five Powell's algorithms of the PRIMA library are provided by the PRIMA package:

• uobyqa (Unconstrained Optimization BY Quadratic Approximations) is for unconstrained optimization, that is Ω = ℝⁿ.

• newuoa is also for unconstrained optimization. According to M.J.D. Powell, newuoa is superior to uobyqa.

• bobyqa (Bounded Optimization BY Quadratic Approximations) is for simple bound constrained problems, that is Ω = { x ∈ ℝⁿ | xl ≤ x ≤ xu }.

• lincoa (LINearly Constrained Optimization) is for constrained optimization problems with bound constraints, linear equality constraints, and linear inequality constraints.

• cobyla (Constrained Optimization BY Linear Approximations) is for general constrained problems with bound constraints, non-linear constraints, linear equality constraints, and linear inequality constraints.

All these algorithms are trust region methods where the variables are updated according to an affine or a quadratic local approximation interpolating the objective function at a given number of points (set by keyword npt by some of the algorithms). No derivatives of the objective function are needed. These algorithms are well suited to problems with a non-analytic objective function that takes time to be evaluated.

The table below summarizes the characteristics of the different Powell's methods ("linear" constraints includes equality and inequality linear constraints).

Method	Model	Constraints
newuoa	quadratic	none
uobyqa	quadratic	none
bobyqa	quadratic	bounds
lincoa	quadratic	bounds, linear
cobyla	affine	bounds, linear, non-linear

To use the most suitable Powell's algorithm depending on the constraints, simply do:

x, info = prima(f, x0; kwds...)

To explicitly use one the specific Powell's algorithms, call one of:

x, info = uobyqa(f, x0; kwds...)
x, info = newuoa(f, x0; kwds...)
x, info = bobyqa(f, x0; kwds...)
x, info = cobyla(f, x0; kwds...)
x, info = lincoa(f, x0; kwds...)

In any case, f is the objective function and x0 specifies the initial values of the variables (and is left unchanged). Constraints and options may be specified by keywords kwds... (see below).

The objective function is called as f(x) with x the variables, it must implement the following signature:

f(x::Vector{Cdouble})::Real

All the algorithms return a 2-tuple (x, info) with x the variables and info a structured object collecting all other information. If issuccess(info) is true, then the algorithm was successful and x is an approximate solution of the problem.

The output info has the following properties:

info.fx       # value of the objective function f(x) on return
info.nf       # number of calls to the objective function
info.status   # final status code
info.cstrv    # amount of constraints violation, 0.0 if unconstrained
info.nl_eq    # non-linear equality constraints, empty vector if none
info.nl_ineq  # non-linear inequality constraints, empty vector if none

Calling one of:

issuccess(info)
issuccess(info.status)

yield whether the algorithm has converged. If this is the case, info.status can be one of:

• PRIMA.SMALL_TR_RADIUS if the radius of the trust region becomes smaller or equal the value of keyword rhobeg, in other words, the algorithm has converged in terms of variable precision;

• PRIMA.FTARGET_ACHIEVED if the objective function is smaller of equal the value of keyword ftarget, in other words, the algorithm has converged in terms of function value.

There are other possibilities which all indicate a failure. Calling one of:

PRIMA.reason(info)
PRIMA.reason(info.status)

yield a textual explanation about the reason that leads the algorithm to stop.

The keywords allowed by the different algorithms are summarized by the following table.

Keyword	Description	Algorithms
rhobeg	Initial trust region radius	all
rhoend	Final trust region radius	all
ftarget	Target objective function value	all
maxfun	Maximum number of function evaluations	all
iprint	Verbosity level	all
npt	Number of points in local model	bobyqa, lincoa, newuoa
xl	Lower bound	bobyqa, cobyla, lincoa
xu	Upper bound	bobyqa, cobyla, lincoa
nonlinear_eq	Non-linear equality constraints	cobyla
nonlinear_ineq	Non-linear inequality constraints	cobyla
linear_eq	Linear equality constraints	cobyla, lincoa
linear_ineq	Linear inequality constraints	cobyla, lincoa

Assuming n = length(x) is the number of variables, then:

• scale (default value nothing) may be set with a vector of n positive scaling factors. If specified, the problem is solved in the scaled variables u ∈ ℝⁿ such that u[i] = x[i]/scale[i]. If unspecified, it is assumed that scale[i] = 1 for all variables. Note that the objective function, the constraints (linear and non-linear), and the bounds remain specified in the variables. Scaling the variables is useful to improve the conditioning of the problem, to make the scaled variables u having approximately the same magnitude, and to adapt to heterogeneous variables or with different units.

• rhobeg (default value 1.0) is the initial radius of the trust region. The radius of the trust region is given by the Euclidean norm of the scaled variables (see keyword scale above).

• rhoend (default value 1e-6*rhobeg) is the final radius of the trust region. The algorithm stops when the trust region radius becomes smaller or equal rhoend and the status PRIMA.SMALL_TR_RADIUS is returned.

• ftarget (default value -Inf) is another convergence setting. The algorithm stops as soon as f(x) ≤ ftarget and the status PRIMA.FTARGET_ACHIEVED is returned.

• iprint (default value PRIMA.MSG_NONE) sets the level of verbosity of the algorithm. Possible values are PRIMA.MSG_EXIT, PRIMA.MSG_RHO, or PRIMA.MSG_FEVL. Note that the values that are printed by the software are those of the scaled variables (see keyword scale above).

• maxfun (default 500*n) is the maximum number of function evaluations allowed for the algorithm. If the number of calls to f(x) exceeds this value, the algorithm is stopped and the status PRIMA.MAXFUN_REACHED is returned.

• npt (default value 2n + 1) is the number of points used to approximate the local behavior of the objective function and such that n + 2 ≤ npt ≤ (n + 1)*(n + 2)/2. The default value corresponds to the one recommended by M.J.D. Powell.

• xl and xu (default fill(+Inf, n) and fill(-Inf, n)) are the element-wise lower and upper bounds for the variables. Feasible variables are such that xl ≤ x ≤ xu.

• nonlinear_eq (default nothing) may be specified with a function, say c_eq, implementing non-linear equality constraints defined by c_eq(x) = 0. On return, the values of the non-linear equality constraints are given by info.nl_eq to save calling c_eq(x).

• nonlinear_ineq (default nothing) may be specified with a function, say c_ineq, implementing non-linear inequality constraints defined by c_ineq(x) ≤ 0. On return, the values of the non-linear inequality constraints are given by info.nl_ineq to save calling c_ineq(x).

• linear_eq (default nothing) may be specified as a tuple (Aₑ,bₑ) to impose the linear equality constraints Aₑ⋅x = bₑ.

• linear_ineq (default nothing) may be specified as a tuple (Aᵢ,bᵢ) to impose the linear inequality constraints Aᵢ⋅x ≤ bᵢ.

References

The following 5 first references respectively describe Powell's algorithms cobyla, uobyqa, newuoa, bobyqa, and lincoa while the last reference provides a good and comprehensive introduction to these algorithms.

1. 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 Mathematics and Its Applications, vol. 275, pp. 51-67 (1994).

2. M.J.D. Powell, "UOBYQA: unconstrained optimization by quadratic approximation", in Mathematical Programming, vol. 92, pp. 555-582 (2002) DOI.

3. M.J.D. Powell, "The NEWUOA software for unconstrained optimization without derivatives", in Nonconvex Optimization and Its Applications, pp. 255-297 (2006).

4. M.J.D. Powell, "The BOBYQA algorithm for bound constrained optimization without derivatives", Department of Applied Mathematics and Theoretical Physics, Cambridge, England, Technical Report NA2009/06 (2009).

5. M.J.D. Powell, "On fast trust region methods for quadratic models with linear constraints", Technical Report of the Department of Applied Mathematics and Theoretical Physics, Cambridge University (2014).

6. T.M. Ragonneau & Z. Zhang, "PDFO: A Cross-Platform Package for Powell's Derivative-Free Optimization Solvers", arXiv (2023).