Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_opt_bounds_bobyqa_func (e04jc)

## Purpose

nag_opt_bounds_bobyqa_func (e04jc) is an easy-to-use algorithm that uses methods of quadratic approximation to find a minimum of an objective function $F$ over $\mathbf{x}\in {R}^{n}$, subject to fixed lower and upper bounds on the independent variables ${x}_{1},{x}_{2},\dots ,{x}_{n}$. Derivatives of $F$ are not required.
The function is intended for functions that are continuous and that have continuous first and second derivatives (although it will usually work even if the derivatives have occasional discontinuities). Efficiency is maintained for large $n$.

## Syntax

[x, f, nf, user, ifail] = e04jc(objfun, npt, x, bl, bu, rhobeg, rhoend, monfun, maxcal, 'n', n, 'user', user)
[x, f, nf, user, ifail] = nag_opt_bounds_bobyqa_func(objfun, npt, x, bl, bu, rhobeg, rhoend, monfun, maxcal, 'n', n, 'user', user)

## Description

nag_opt_bounds_bobyqa_func (e04jc) is applicable to problems of the form:
 $minimize x∈Rn Fx subject to ℓ ≤ x ≤ u and ℓ≤u ,$
where $F$ is a nonlinear scalar function whose derivatives may be unavailable, and where the bound vectors are elements of ${R}^{n}$. Relational operators between vectors are interpreted elementwise.
Fixing variables (that is, setting ${\ell }_{i}={u}_{i}$ for some $i$) is allowed in nag_opt_bounds_bobyqa_func (e04jc).
You must supply a function to calculate the value of $F$ at any given point $\mathbf{x}$.
The method used by nag_opt_bounds_bobyqa_func (e04jc) is based on BOBYQA, the method of Bound Optimization BY Quadratic Approximation described in Powell (2009). In particular, each iteration of nag_opt_bounds_bobyqa_func (e04jc) generates a quadratic approximation $Q$ to $F$ that agrees with $F$ at $m$ automatically chosen interpolation points. The value of $m$ is a constant prescribed by you. Updates to the independent variables mostly occur from approximate solutions to trust-region subproblems, using the current quadratic model.

## References

Powell M J D (2009) The BOBYQA algorithm for bound constrained optimization without derivatives Report DAMTP 2009/NA06 University of Cambridge http://www.damtp.cam.ac.uk/user/na/NA_papers/NA2009_06.pdf

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{objfun}$ – function handle or string containing name of m-file
objfun must evaluate the objective function $F$ at a specified vector $\mathbf{x}$.
[f, user, inform] = objfun(n, x, user)

Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
$n$, the number of independent variables.
2:     $\mathrm{x}\left({\mathbf{n}}\right)$ – double array
$\mathbf{x}$, the vector at which the objective function is to be evaluated.
3:     $\mathrm{user}$ – Any MATLAB object
objfun is called from nag_opt_bounds_bobyqa_func (e04jc) with the object supplied to nag_opt_bounds_bobyqa_func (e04jc).

Output Parameters

1:     $\mathrm{f}$ – double scalar
Must be set to the value of the objective function at $\mathbf{x}$, unless you have specified termination of the current problem using inform.
2:     $\mathrm{user}$ – Any MATLAB object
3:     $\mathrm{inform}$int64int32nag_int scalar
Must be set to a value describing the action to be taken by the solver on return from objfun. Specifically, if the value is negative the solution of the current problem will terminate immediately; otherwise, computations will continue.
2:     $\mathrm{npt}$int64int32nag_int scalar
Suggested value: ${\mathbf{npt}}=2×{n}_{r}+1$, where ${n}_{r}$ denotes the number of non-fixed variables.
$m$, the number of interpolation conditions imposed on the quadratic approximation at each iteration.
Constraint: ${n}_{r}+2\le {\mathbf{npt}}\le \frac{\left({n}_{r}+1\right)×\left({n}_{r}+2\right)}{2}$, where ${n}_{r}$ denotes the number of non-fixed variables.
3:     $\mathrm{x}\left({\mathbf{n}}\right)$ – double array
An estimate of the position of the minimum. If any component is out-of-bounds it is replaced internally by the bound it violates.
4:     $\mathrm{bl}\left({\mathbf{n}}\right)$ – double array
5:     $\mathrm{bu}\left({\mathbf{n}}\right)$ – double array
The fixed vectors of bounds: the lower bounds $\mathbf{\ell }$ and the upper bounds $\mathbf{u}$, respectively. To signify that a variable is unbounded you should choose a large scalar $r$ appropriate to your problem, then set the lower bound on that variable to $-r$ and the upper bound to $r$. For well-scaled problems $r={r}_{\mathrm{max}}^{\frac{1}{4}}$ may be suitable, where ${r}_{\mathrm{max}}$ denotes the largest positive model number (see nag_machine_real_largest (x02al)).
Constraints:
• if ${\mathbf{x}}\left(\mathit{i}\right)$ is to be fixed at ${\mathbf{bl}}\left(\mathit{i}\right)$, then ${\mathbf{bl}}\left(\mathit{i}\right)={\mathbf{bu}}\left(\mathit{i}\right)$;
• otherwise ${\mathbf{bu}}\left(\mathit{i}\right)-{\mathbf{bl}}\left(\mathit{i}\right)\ge 2.0×{\mathbf{rhobeg}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{n}}$.
6:     $\mathrm{rhobeg}$ – double scalar
Suggested value: rhobeg should be about one tenth of the greatest expected overall change to a variable: the initial quadratic model will be constructed by taking steps from the initial x of length rhobeg along each coordinate direction.
An initial lower bound on the value of the trust-region radius.
Constraints:
• ${\mathbf{rhobeg}}>0.0$;
• ${\mathbf{rhobeg}}\ge {\mathbf{rhoend}}$.
7:     $\mathrm{rhoend}$ – double scalar
Suggested value: rhoend should indicate the absolute accuracy that is required in the final values of the variables.
A final lower bound on the value of the trust-region radius.
Constraint: ${\mathbf{rhoend}}>0.0$.
8:     $\mathrm{monfun}$ – function handle or string containing name of m-file
monfun may be used to monitor the optimization process. It is invoked every time a new trust-region radius is chosen.
If no monitoring is required, monfun may be string nag_opt_bounds_bobyqa_func_dummy_monfun (e04jcp)
[user, inform] = monfun(n, nf, x, f, rho, user)

Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
$n$, the number of independent variables.
2:     $\mathrm{nf}$int64int32nag_int scalar
The cumulative number of calls made to objfun.
3:     $\mathrm{x}\left({\mathbf{n}}\right)$ – double array
The current best point.
4:     $\mathrm{f}$ – double scalar
The value of objfun at x.
5:     $\mathrm{rho}$ – double scalar
A lower bound on the current trust-region radius.
6:     $\mathrm{user}$ – Any MATLAB object
monfun is called from nag_opt_bounds_bobyqa_func (e04jc) with the object supplied to nag_opt_bounds_bobyqa_func (e04jc).

Output Parameters

1:     $\mathrm{user}$ – Any MATLAB object
2:     $\mathrm{inform}$int64int32nag_int scalar
Must be set to a value describing the action to be taken by the solver on return from monfun. Specifically, if the value is negative the solution of the current problem will terminate immediately; otherwise, computations will continue.
9:     $\mathrm{maxcal}$int64int32nag_int scalar
The maximum permitted number of calls to objfun.
Constraint: ${\mathbf{maxcal}}\ge 1$.

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the dimension of the arrays x, bl, bu. (An error is raised if these dimensions are not equal.)
$n$, the number of independent variables.
Constraint: ${\mathbf{n}}\ge 2$ and ${n}_{r}\ge 2$, where ${n}_{r}$ denotes the number of non-fixed variables.
2:     $\mathrm{user}$ – Any MATLAB object
user is not used by nag_opt_bounds_bobyqa_func (e04jc), but is passed to objfun and monfun. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

### Output Parameters

1:     $\mathrm{x}\left({\mathbf{n}}\right)$ – double array
The lowest point found during the calculations. Thus, if ${\mathbf{ifail}}={\mathbf{0}}$ on exit, x is the position of the minimum.
2:     $\mathrm{f}$ – double scalar
The function value at the lowest point found (x).
3:     $\mathrm{nf}$int64int32nag_int scalar
Unless ${\mathbf{ifail}}={\mathbf{1}}$ or $-{\mathbf{999}}$ on exit, the total number of calls made to objfun.
4:     $\mathrm{user}$ – Any MATLAB object
5:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).
nag_opt_bounds_bobyqa_func (e04jc) returns with ${\mathbf{ifail}}={\mathbf{0}}$ if the final trust-region radius has reached its lower bound rhoend.

## Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\mathbf{ifail}}=1$
Constraint: if ${\mathbf{bl}}\left(i\right)\ne {\mathbf{bu}}\left(i\right)$ in coordinate $i$, then ${\mathbf{bu}}\left(i\right)-{\mathbf{bl}}\left(i\right)\ge 2×{\mathbf{rhobeg}}$.
Constraint: ${\mathbf{maxcal}}\ge 1$.
Constraint: ${\mathbf{rhobeg}}>0.0$.
Constraint: ${\mathbf{rhoend}}>0.0$.
Constraint: ${\mathbf{rhoend}}\le {\mathbf{rhobeg}}$.
There were ${n}_{r}=_$ unequal bounds and ${\mathbf{npt}}=_$ on entry.
Constraint: ${n}_{r}+2\le {\mathbf{npt}}\le \frac{\left({n}_{r}+1\right)×\left({n}_{r}+2\right)}{2}$.
There were unequal bounds.
Constraint: ${n}_{r}\ge 2$.
${\mathbf{ifail}}=2$
The function evaluations limit was reached: objfun has been called maxcal times.
${\mathbf{ifail}}=3$
The predicted reduction in a trust-region step was non-positive. Check your specification of objfun and whether the function needs rescaling. Try a different initial x.
${\mathbf{ifail}}=4$
A rescue procedure has been called in order to correct damage from rounding errors when computing an update to a quadratic approximation of $F$, but no further progess could be made. Check your specification of objfun and whether the function needs rescaling. Try a different initial x.
W  ${\mathbf{ifail}}=5$
User-supplied monitoring function requested termination.
User-supplied objective function requested termination.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

Experience shows that, in many cases, on successful termination the $\infty$-norm distance from the best point $\mathbf{x}$ to a local minimum of $F$ is less than $10×{\mathbf{rhoend}}$, unless rhoend is so small that such accuracy is unattainable.

For each invocation of nag_opt_bounds_bobyqa_func (e04jc), local workspace arrays of fixed length are allocated internally. The total size of these arrays amounts to $\left({\mathbf{npt}}+6\right)×\left({\mathbf{npt}}+{n}_{r}\right)+\frac{{n}_{r}×\left(3{n}_{r}+21\right)}{2}$ double elements and ${n}_{r}$ integer elements, where ${n}_{r}$ denotes the number of non-fixed variables; that is, the total size is $\mathcal{O}\left({n}_{r}^{4}\right)$. If you follow the recommendation for the choice of npt on entry, this total size reduces to $\mathcal{O}\left({n}_{r}^{2}\right)$.
Usually the total number of function evaluations (nf) is substantially less than $\mathcal{O}\left({n}_{r}^{2}\right)$, and often, if ${\mathbf{npt}}=2×{n}_{r}+1$ on entry, nf is only of magnitude ${n}_{r}$ or less.

## Example

This example involves the minimization of
 $F = x1+10x2 2 +5⁢ x3 - x4 2 + x2-2 x3 4 +10⁢ x1 - x4 4$
subject to
 $-1≤x1≤ 3, -2≤x2≤ 0, -1≤x4≤ 3,$
starting from the initial guess $\left(3,-1,0,1\right)$.
```function e04jc_example

fprintf('e04jc example results\n\n');

maxcal = int64(500);
rhobeg = 0.1;
rhoend = 1.0e-6;
n = 4;
npt = int64(2*n + 1);

% x(3) is unconstrained, so we're going to set bl(3) to a large
% negative number and bu(3) to a large positive number.
infbnd = x02al^0.25;
bl = [1, -2, -infbnd, 1];
bu = [3, 0, infbnd, 3];
x  = [3, -1, 0, 1];

[x, f, nf, user, ifail] = ...
e04jc( ...
@objfun, npt, x, bl, bu, rhobeg, rhoend, @monfun, maxcal);

fprintf('\nFunction value at lowest point found =  %12.5e\n',f);
fprintf(' The corresponding X is:\n');
fprintf('%14.5e',x);
fprintf('\n');

function [f, user, inform] = objfun(n, x, user)
inform = int64(0);

f = (x(1)+10*x(2))^2 + 5*(x(3)-x(4))^2 + (x(2)-2*x(3))^4 + 10*(x(1)-x(4))^4;

function [user, inform] = monfun(n, nf, x, f, rho, user)
inform = int64(0);

fprintf('\nNew rho = %13.5e, number of function evaluations = %d\n', rho, nf);
fprintf('Current function value = %13.5e\n', f);
fprintf('The corresponding X is:');
fprintf(' %13.5e', x);
fprintf('\n');
```
```e04jc example results

New rho =   1.00000e-02, number of function evaluations = 25
Current function value =   4.09399e+00
The corresponding X is:   1.60106e+00  -1.03604e-01   4.51135e-01   1.02335e+00

New rho =   1.00000e-03, number of function evaluations = 67
Current function value =   2.43397e+00
The corresponding X is:   1.00000e+00  -8.59741e-02   4.06744e-01   1.00000e+00

New rho =   1.00000e-04, number of function evaluations = 77
Current function value =   2.43379e+00
The corresponding X is:   1.00000e+00  -8.52328e-02   4.09342e-01   1.00000e+00

New rho =   1.00000e-05, number of function evaluations = 81
Current function value =   2.43379e+00
The corresponding X is:   1.00000e+00  -8.52328e-02   4.09342e-01   1.00000e+00

New rho =   1.00000e-06, number of function evaluations = 93
Current function value =   2.43379e+00
The corresponding X is:   1.00000e+00  -8.52329e-02   4.09303e-01   1.00000e+00

Function value at lowest point found =   2.43379e+00
The corresponding X is:
1.00000e+00  -8.52326e-02   4.09303e-01   1.00000e+00
```