e05jb is designed to find the global minimum or maximum of an arbitrary function, subject to simple bound-constraints using a multi-level coordinate search method. Derivatives are not required, but convergence is only guaranteed if the objective function is continuous in a neighbourhood of a global optimum. It is not intended for large problems.

The initialization method (E05JAF not in this release) must have been called before calling e05jb.

Syntax

C#
public static void e05jb( int n, E05..::..E05JB_OBJFUN objfun, int ibound, int iinit, double[] bl, double[] bu, double[,] list, int[] numpts, int[] initpt, E05..::..E05JB_MONIT monit, double[] x, out double obj, E05..::..e05jbOptions options, out int ifail )

public static void e05jb(
	int n,
	E05..::..E05JB_OBJFUN objfun,
	int ibound,
	int iinit,
	double[] bl,
	double[] bu,
	double[,] list,
	int[] numpts,
	int[] initpt,
	E05..::..E05JB_MONIT monit,
	double[] x,
	out double obj,
	E05..::..e05jbOptions options,
	out int ifail
)

Visual Basic
Public Shared Sub e05jb ( _ n As Integer, _ objfun As E05..::..E05JB_OBJFUN, _ ibound As Integer, _ iinit As Integer, _ bl As Double(), _ bu As Double(), _ list As Double(,), _ numpts As Integer(), _ initpt As Integer(), _ monit As E05..::..E05JB_MONIT, _ x As Double(), _ <OutAttribute> ByRef obj As Double, _ options As E05..::..e05jbOptions, _ <OutAttribute> ByRef ifail As Integer _ )

Visual Basic

Public Shared Sub e05jb ( _
	n As Integer, _
	objfun As E05..::..E05JB_OBJFUN, _
	ibound As Integer, _
	iinit As Integer, _
	bl As Double(), _
	bu As Double(), _
	list As Double(,), _
	numpts As Integer(), _
	initpt As Integer(), _
	monit As E05..::..E05JB_MONIT, _
	x As Double(), _
	<OutAttribute> ByRef obj As Double, _
	options As E05..::..e05jbOptions, _
	<OutAttribute> ByRef ifail As Integer _
)

Visual C++
public: static void e05jb( int n, E05..::..E05JB_OBJFUN^ objfun, int ibound, int iinit, array<double>^ bl, array<double>^ bu, array<double,2>^ list, array<int>^ numpts, array<int>^ initpt, E05..::..E05JB_MONIT^ monit, array<double>^ x, [OutAttribute] double% obj, E05..::..e05jbOptions^ options, [OutAttribute] int% ifail )

Visual C++

public:
static void e05jb(
	int n, 
	E05..::..E05JB_OBJFUN^ objfun, 
	int ibound, 
	int iinit, 
	array<double>^ bl, 
	array<double>^ bu, 
	array<double,2>^ list, 
	array<int>^ numpts, 
	array<int>^ initpt, 
	E05..::..E05JB_MONIT^ monit, 
	array<double>^ x, 
	[OutAttribute] double% obj, 
	E05..::..e05jbOptions^ options, 
	[OutAttribute] int% ifail
)

F#
static member e05jb : n : int * objfun : E05..::..E05JB_OBJFUN * ibound : int * iinit : int * bl : float[] * bu : float[] * list : float[,] * numpts : int[] * initpt : int[] * monit : E05..::..E05JB_MONIT * x : float[] * obj : float byref * options : E05..::..e05jbOptions * ifail : int byref -> unit

static member e05jb : 
        n : int * 
        objfun : E05..::..E05JB_OBJFUN * 
        ibound : int * 
        iinit : int * 
        bl : float[] * 
        bu : float[] * 
        list : float[,] * 
        numpts : int[] * 
        initpt : int[] * 
        monit : E05..::..E05JB_MONIT * 
        x : float[] * 
        obj : float byref * 
        options : E05..::..e05jbOptions * 
        ifail : int byref -> unit

Parameters

n: Type: System..::..Int32
On entry: $n$ , the number of variables.

Constraint: $n > 0$ .

objfun: Type: NagLibrary..::..E05..::..E05JB_OBJFUN
objfun must evaluate the objective function $F (x)$ for a specified $n$ -vector $x$ .
A delegate of type E05JB_OBJFUN.

ibound

Type: System..::..Int32

On entry: indicates whether the facility for dealing with bounds of special forms is to be used. ibound must be set to one of the following values.

$ibound = 0$: You will supply $ℓ$ and $u$ individually.
$ibound = 1$: There are no bounds on $x$ .
$ibound = 2$: There are semi-infinite bounds $0 \leq x$ .
$ibound = 3$: There are constant bounds $ℓ = ℓ_{1}$ and $u = u_{1}$ .

Note that it only makes sense to fix any components of

x

when

ibound = 0

Constraint:

ibound = 0

1

2

3

iinit

Type: System..::..Int32

On entry: selects which initialization method to use.

$iinit = 0$: Simple initialization (boundary and midpoint), with
$numpts [i - 1] = 3$ , $initpt [i - 1] = 2$ and
$list [i - 1, j - 1] = (bl [i - 1], (bl [i - 1] + bu [i - 1]) / 2, bu [i - 1])$ ,
for $i = 1, 2, \dots, n$ and $j = 1, 2, 3$ .
$iinit = 1$: Simple initialization (off-boundary and midpoint), with
$numpts [i - 1] = 3$ , $initpt [i - 1] = 2$ and
$list [i - 1, j - 1] = ((5 bl [i - 1] + bu [i - 1]) / 6, (bl [i - 1] + bu [i - 1]) / 2, (bl [i - 1] + 5 bu [i - 1]) / 6)$ ,
for $i = 1, 2, \dots, n$ and $j = 1, 2, 3$ .
$iinit = 2$: Initialization using linesearches.
$iinit = 3$: You are providing your own initialization list.
$iinit = 4$: Generate a random initialization list.

For more information on methods

iinit = 2

3

4

see [Initialization and Sweeps].

If ‘infinite’ values (as determined by the value of the optional parameter Infinite Bound Size) are detected by e05jb when you are using a simple initialization method (

iinit = 0

1

), a safeguarded initialization procedure will be attempted, to avoid overflow.

Suggested value:

iinit = 0

Constraint:

iinit = 0

1

2

3

4

bl

Type: array<System..::..Double>[]()[][]

An array of size [n]

On entry:

bl

ℓ

, the array of lower bounds.

bu

u

, the array of upper bounds.

ibound = 0

, you must set

bl [i - 1]

ℓ_{i}

and

bu [i - 1]

u_{i}

, for

i = 1, 2, \dots, n

. If a particular

x_{i}

is to be unbounded below, the corresponding

bl [i - 1]

should be set to

- infbnd

, where

infbnd

is the value of the optional parameter Infinite Bound Size. Similarly, if a particular

x_{i}

is to be unbounded above, the corresponding

bu [i - 1]

should be set to

infbnd

ibound = 1

2

, arrays bl and bu need not be set on input.

ibound = 3

, you must set

bl [0]

ℓ_{1}

and

bu [0]

u_{1}

. The remaining elements of bl and bu will then be populated by these initial values.

On exit: unless

ifail = 1

2

on exit, bl and bu are the actual arrays of bounds used by e05jb.

Constraints:

if $ibound = 0$ , $bl [i - 1] \leq bu [i - 1]$ , for $i = 1, 2, \dots, n$ ;
if $ibound = 3$ , $bl [0] < bu [0]$ .

bu

Type: array<System..::..Double>[]()[][]

An array of size [n]

On entry:

bl

ℓ

, the array of lower bounds.

bu

u

, the array of upper bounds.

ibound = 0

, you must set

bl [i - 1]

ℓ_{i}

and

bu [i - 1]

u_{i}

, for

i = 1, 2, \dots, n

. If a particular

x_{i}

is to be unbounded below, the corresponding

bl [i - 1]

should be set to

- infbnd

, where

infbnd

is the value of the optional parameter Infinite Bound Size. Similarly, if a particular

x_{i}

is to be unbounded above, the corresponding

bu [i - 1]

should be set to

infbnd

ibound = 1

2

, arrays bl and bu need not be set on input.

ibound = 3

, you must set

bl [0]

ℓ_{1}

and

bu [0]

u_{1}

. The remaining elements of bl and bu will then be populated by these initial values.

On exit: unless

ifail = 1

2

on exit, bl and bu are the actual arrays of bounds used by e05jb.

Constraints:

if $ibound = 0$ , $bl [i - 1] \leq bu [i - 1]$ , for $i = 1, 2, \dots, n$ ;
if $ibound = 3$ , $bl [0] < bu [0]$ .

list: Type: array<System..::..Double,2>[,](,)[,][,]
An array of size [n, sdlist]
On entry: this parameter need not be set on entry if you wish to use one of the preset initialization methods ( $iinit \neq 3$ ).
list is the ‘initialization list’: whenever a sub-box in the algorithm is split for the first time (either during the initialization procedure or later), for each non-fixed coordinate $i$ the split is done at the values $list [i - 1, 0 : numpts [i - 1] - 1]$ , as well as at some adaptively chosen intermediate points. The array sections $list [i - 1, 0 : numpts [i - 1] - 1]$ , for $i = 1, 2, \dots, n$ , must be in ascending order with each entry being distinct. In this context, ‘distinct’ should be taken to mean relative to the safe-range parameter (see x02am).

On exit: unless $ifail = 1$ , $2$ or $- 999$ on exit, the actual initialization data used by e05jb. If you wish to monitor the contents of list you are advised to do so solely through monit, not through the output value here.

Constraint: if $x [i - 1]$ is not fixed, $list [i - 1, 0 : numpts [i - 1] - 1]$ is in ascending order with each entry being distinct, for $i = 1, 2, \dots, n$ $bl [i - 1] \leq list [i - 1, j - 1] \leq bu [i - 1]$ , for $i = 1, 2, \dots, n$ and $j = 1, 2, \dots, numpts [i - 1]$ .

numpts

Type: array<System..::..Int32>[]()[][]

An array of size [n]

On entry: this parameter need not be set on entry if you wish to use one of the preset initialization methods (

iinit \neq 3

numpts encodes the number of splitting points in each non-fixed dimension.

On exit: unless

ifail = 1

2

- 999

on exit, the actual initialization data used by e05jb.

Constraints:

if $x [i - 1]$ is not fixed, $numpts [i - 1] \leq sdlist$ ;
$numpts [i - 1] \geq 3$ , for $i = 1, 2, \dots, n$ .

initpt: Type: array<System..::..Int32>[]()[][]
An array of size [n]
On entry: this parameter need not be set on entry if you wish to use one of the preset initialization methods ( $iinit \neq 3$ ).
You must designate a point stored in list that you wish e05jb to consider as an ‘initial point’ for the purposes of the splitting procedure. Call this initial point $x^{*}$ . The coordinates of $x^{*}$ correspond to a set of indices $J_{i}$ , for $i = 1, 2, \dots, n$ , such that $x_{i}^{*}$ is stored in $list [i - 1, J_{i} - 1]$ , for $i = 1, 2, \dots, n$ . You must set $initpt [i - 1] = J_{i}$ , for $i = 1, 2, \dots, n$ .

On exit: unless $ifail = 1$ , $2$ or $- 999$ on exit, the actual initialization data used by e05jb.

Constraint: if $x [i - 1]$ is not fixed, $1 \leq initpt [i - 1] \leq sdlist$ , for $i = 1, 2, \dots, n$ .

monit: Type: NagLibrary..::..E05..::..E05JB_MONIT
monit may be used to monitor the optimization process. It is invoked upon every successful completion of the procedure in which a sub-box is considered for splitting. It will also be called just before e05jb exits if that splitting procedure was not successful.
If no monitoring is required, monit may be the dummy monitoring method E05JBK supplied by the NAG Library.
A delegate of type E05JB_MONIT.

x: Type: array<System..::..Double>[]()[][]
An array of size [n]
On exit: if $ifail = 0$ , contains an estimate of the global optimum (see also [Accuracy]).

obj: Type: System..::..Double%
On exit: if $ifail = 0$ , contains the function value at x.
If you request early termination of e05jb using inform in objfun or the analogous inform in monit, there is no guarantee that the function value at x equals obj.

options: Type: NagLibrary..::..E05..::..e05jbOptions
An Object of type E05.e05jbOptions. Used to configure optional parameters to this method.

ifail: Type: System..::..Int32%
On exit: $ifail = 0$ unless the method detects an error or a warning has been flagged (see [Error Indicators and Warnings]).

Description

e05jb is designed to solve modestly sized global optimization problems having simple bound-constraints only; it finds the global optimum of a nonlinear function subject to a set of bound constraints on the variables. Without loss of generality, the problem is assumed to be stated in the following form:

\underset{x \in R^{n}}{minimize} F (x) subject to ℓ \leq x \leq u and ℓ \leq u,

where

F (x)

(the objective function) is a nonlinear scalar function (assumed to be continuous in a neighbourhood of a global minimum), and the bound vectors are elements of

{\bar{R}}^{n}

, where

\bar{R}

denotes the extended reals

R \cup \{- \infty, \infty\}

. Relational operators between vectors are interpreted elementwise.

The optional parameter Maximize should be set if you wish to solve maximization, rather than minimization, problems.

If certain bounds are not present, the associated elements of

ℓ

u

can be set to special values that will be treated as

- \infty

+ \infty

. See the description of the optional parameter Infinite Bound Size. Phrases in this document containing terms like ‘unbounded values’ should be understood to be taken relative to this optional parameter.

Fixing variables (that is, setting

l_{i} = u_{i}

for some

i

) is allowed in e05jb.

A typical excerpt from a method calling e05jb is: where (E05JDF not in this release) sets the optional parameter and value specified in optstr.

The initialization method (E05JAF not in this release) does not need to be called before each invocation of e05jb. You should be aware that a call to the initialization method will reset each optional parameter to its default value, and, if you are using repeatable randomized initialization lists (see the description of the parameter iinit), the random state stored in the array comm will be destroyed.

You must supply a method that evaluates

F (x)

; derivatives are not required.

The method used by e05jb is based on MCS, the Multi-level Coordinate Search method described in Huyer and Neumaier (1999), and the algorithm it uses is described in detail in [Algorithmic Details].

References

Huyer W and Neumaier A (1999) Global optimization by multi-level coordinate search Journal of Global Optimization 14 331–355

Error Indicators and Warnings

Errors or warnings detected by the method:

$ifail = 1$: Initialization method (E05JAF not in this release) has not been called.

On entry, $lcomm = 〈value〉$ .
Constraint: $lcomm \geq 100$ .

$ifail = 2$: A value of Splits Limit ( $smax$ ) smaller than $n_{r} + 3$ was set: $smax = 〈value〉$ , $n_{r} = 〈value〉$ .

On entry, $ibound = 〈value〉$ .
Constraint: $ibound = 0$ , $1$ , $2$ or $3$ .

On entry, $ibound = 0$ or $3$ and $bl [i - 1] = 〈value〉$ , $bu [i - 1] = 〈value〉$ and $i = 〈value〉$ .
Constraint: if $ibound = 0$ then $bl [i - 1] \leq bu [i - 1]$ , for $i = 1, 2, \dots, n$ ; if $ibound = 3$ then $bl [0] < bu [0]$ .

On entry, $ibound = 3$ and $bl [0] = bu [0] = 〈value〉$ .
Constraint: if $ibound = 3$ then $bl [0] < bu [0]$ .

On entry, $iinit = 〈value〉$ .
Constraint: $iinit = 0$ , $1$ , $2$ , $3$ or $4$ .

On entry, $iinit = 2$ and $sdlist = 〈value〉$ .
Constraint: if $iinit = 2$ then $sdlist \geq 192$ .

On entry, $iinit = 〈value〉$ and $sdlist = 〈value〉$ .
Constraint: if $iinit \neq 2$ then $sdlist \geq 3$ .

On entry, $n = 〈value〉$ .
Constraint: $n > 0$ .

On entry, user-supplied $initpt [i - 1] = 〈value〉$ , $i = 〈value〉$ .
Constraint: if $x [i - 1]$ is not fixed then $initpt [i - 1] \geq 1$ , for $i = 1, 2, \dots, n$ .

On entry, user-supplied $initpt [i - 1] = 〈value〉$ , $i = 〈value〉$ and $sdlist = 〈value〉$ .
Constraint: if $x [i - 1]$ is not fixed then $initpt [i - 1] \leq sdlist$ , for $i = 1, 2, \dots, n$ .

On entry, user-supplied $list [i - 1, j - 1] = 〈value〉$ , $i = 〈value〉$ , $j = 〈value〉$ , and $bl [i - 1] = 〈value〉$ .
Constraint: if $x [i - 1]$ is not fixed then $list [i - 1, j - 1] \geq bl [i - 1]$ , for $i = 1, 2, \dots, n$ and $j = 1, 2, \dots, numpts [i - 1]$ .

On entry, user-supplied $list [i - 1, j - 1] = 〈value〉$ , $i = 〈value〉$ , $j = 〈value〉$ , and $bu [i - 1] = 〈value〉$ .
Constraint: if $x [i - 1]$ is not fixed then $list [i - 1, j - 1] \leq bu [i - 1]$ , for $i = 1, 2, \dots, n$ and $j = 1, 2, \dots, numpts [i - 1]$ .

On entry, user-supplied $numpts [i - 1] = 〈value〉$ , $i = 〈value〉$ .
Constraint: if $x [i - 1]$ is not fixed then $numpts [i - 1] \geq 3$ , for $i = 1, 2, \dots, n$ .

On entry, user-supplied $numpts [i - 1] = 〈value〉$ , $i = 〈value〉$ and $sdlist = 〈value〉$ .
Constraint: if $x [i - 1]$ is not fixed then $numpts [i - 1] \leq sdlist$ , for $i = 1, 2, \dots, n$ .

On entry, user-supplied section $list [i, 1 : numpts [i - 1]]$ contained $ndist$ distinct elements, and $ndist < numpts [i - 1]$ : $ndist = 〈value〉$ , $numpts [i - 1] = 〈value〉$ , $i = 〈value〉$ .

On entry, user-supplied section $list [i, 1 : numpts [i - 1]]$ was not in ascending order: $numpts [i - 1] = 〈value〉$ , $i = 〈value〉$ .

The number of non-fixed variables $n_{r} = 0$ .
Constraint: $n_{r} > 0$ .

$ifail = 3$: A finite initialization list could not be computed internally. Consider reformulating the bounds on the problem, try providing your own initialization list, use the randomization option ( $iinit = 4$ ) or vary the value of Infinite Bound Size.

The user-supplied initialization list contained infinite values, as determined by the optional parameter Infinite Bound Size.

$ifail = 4$: The division procedure completed but your target value could not be reached.
Despite every sub-box being processed Splits Limit times, the target value you provided in Target Objective Value could not be found to the tolerances given in Target Objective Error and Target Objective Safeguard. You could try reducing Splits Limit or the objective tolerances.

$ifail = 5$: The function evaluations limit was exceeded.
Approximately Function Evaluations Limit function calls have been made without your chosen termination criterion being satisfied.

$ifail = 6$: User-supplied monitoring method requested termination.

User-supplied objective function requested termination.

$ifail = 7$: An error occurred during initialization. It is likely that points from the initialization list are very close together. Try relaxing the bounds on the variables or use a different initialization method.

An error occurred during linesearching. It is likely that your objective function is badly scaled: try rescaling it. Also, try relaxing the bounds or use a different initialization method. If the problem persists, please contact NAG quoting error code $〈value〉$ .

$ifail = - 999$: Dynamic memory allocation failed.

$ifail = -9000$: An error occured, see message report.
$ifail = -6000$: Invalid Parameters $〈value〉$
$ifail = -4000$: Invalid dimension for array $〈value〉$
$ifail = -8000$: Negative dimension for array $〈value〉$
$ifail = -6000$: Invalid Parameters $〈value〉$

Accuracy

ifail = 0

on exit, then the vector returned in the array x is an estimate of the solution

x

whose function value satisfies your termination criterion: the function value was static for Static Limit sweeps through levels, or

F (x) - objval \leq \max (objerr \times |objval|, objsfg),

where

objval

is the value of the optional parameter Target Objective Value,

objerr

is the value of the optional parameter Target Objective Error, and

objsfg

is the value of the optional parameter Target Objective Safeguard.

Parallelism and Performance

None.

Further Comments

For each invocation of e05jb, local workspace arrays of fixed length are allocated internally. The total size of these arrays amounts to

13 n_{r} + smax - 1

integer elements, where

smax

is the value of the optional parameter Splits Limit and

n_{r}

is the number of non-fixed variables, and

(2 + n_{r}) sdlist + 2 n + 21 n_{r} + 3 n_{r}^{2} + 1

real elements. In addition, if you are using randomized initialization lists (see the description of the parameter iinit), a further

21

integer elements are allocated internally.

In order to keep track of the regions of the search space that have been visited while looking for a global optimum, e05jb internally allocates arrays of increasing sizes depending on the difficulty of the problem. Two of the main factors that govern the amount allocated are the number of sub-boxes (call this quantity

nboxes

) and the number of points in the ‘shopping basket’ (the parameter nbaskt on entry to monit). Safe, pessimistic upper bounds on these two quantities are so large as to be impractical. In fact, the worst-case number of sub-boxes for even the most simple initialization list (when

ninit = 3

on entry to monit) grows like

{n_{r}}^{n_{r}}

. Thus e05jb does not attempt to estimate in advance the final values of

nboxes

or nbaskt for a given problem. There are a total of

5

integer arrays and

4 + n_{r} + ninit

real arrays whose lengths depend on

nboxes

, and there are a total of

2

integer arrays and

3 + n + n_{r}

real arrays whose lengths depend on nbaskt. e05jb makes a fixed initial guess that the maximum number of sub-boxes required will be

10000

and that the maximum number of points in the ‘shopping basket’ will be

1000

. If ever a greater amount of sub-boxes or more room in the ‘shopping basket’ is required, e05jb performs reallocation, usually doubling the size of the inadequately-sized arrays. Clearly this process requires periods where the original array and its extension exist in memory simultaneously, so that the data within can be copied, which compounds the complexity of e05jb's memory usage. It is possible (although not likely) that if your problem is particularly difficult to solve, or of a large size (hundreds of variables), you may run out of memory.

One array that could be dynamically resized by e05jb is the ‘shopping basket’ (xbaskt on entry to monit). If the initial attempt to allocate

1000 n_{r}

reals for this array fails, monit will not be called on exit from e05jb.

e05jb performs better if your problem is well-scaled. It is worth trying (by guesswork perhaps) to rescale the problem if necessary, as sensible scaling will reduce the difficulty of the optimization problem, so that e05jb will take less computer time.

Example

This example finds the global minimum of the ‘peaks’ function in two dimensions

F (x, y) = 3 {(1 - x)}^{2} \exp (- x^{2} - {(y + 1)}^{2}) - 10 (\frac{x}{5} - x^{3} - y^{5}) \exp (- x^{2} - y^{2}) - \frac{1}{3} \exp (- {(x + 1)}^{2} - y^{2})

on the box

[- 3, 3] \times [- 3, 3]

The function

F

has several local minima and one global minimum in the given box. The global minimum is approximately located at

(0.23, - 1.63)

, where the function value is approximately

- 6.55

We use default values for all the optional parameters, and we instruct e05jb to use the simple initialization list corresponding to

iinit = 0

. In particular, this will set for us the initial point

(0, 0)

(see []).

Example program (C#): e05jbe.cs

Example program data: e05jbe.d

Example program results: e05jbe.r

Algorithmic Details

Here we summarise the main features of the MCS algorithm used in e05jb, and we introduce some terminology used in the description of the method and its arguments. We assume throughout that we will only do any work in coordinates

i

in which

x_{i}

is free to vary. The MCS algorithm is fully described in Huyer and Neumaier (1999).

Initialization and Sweeps

Each sub-box is determined by a basepoint

x

and an opposite point

y

. We denote such a sub-box by

B [x, y]

. The basepoint is allowed to belong to more than one sub-box, is usually a boundary point, and is often a vertex.

An initialization procedure produces an initial set of sub-boxes. Whenever a sub-box is split along a coordinate

i

for the first time (in the initialization procedure or later), the splitting is done at three or more user-defined values

{\{x_{i}^{j}\}}_{j}

at which the objective function is sampled, and at some adaptively chosen intermediate points. At least four children are generated. More precisely, we assume that we are given

ℓ_{i} \leq x_{i}^{1} < x_{i}^{2} < \dots < x_{i}^{L_{i}} \leq u_{i}, L_{i} \geq 3,  for ​ i = 1, 2, \dots, n

and a vector

p

that, for each

i

, locates within

{\{x_{i}^{j}\}}_{j}

the

i

th coordinate of an initial point

x^{0}

; that is, if

x_{i}^{0} = x_{i}^{j}

for some

j = 1, 2, \dots, L_{i}

, then

p_{i} = j

. A good guess for the global optimum can be used as

x^{0}

The initialization points and the vectors

ℓ

and

p

are collectively called the initialization list (and sometimes we will refer to just the initialization points as ‘the initialization list’, whenever this causes no confusion). The initialization data may be input by you, or they can be set to sensible default values by e05jb: if you provide them yourself,

list [i - 1, j - 1]

should contain

x_{i}^{j}

numpts [i - 1]

should contain

L_{i}

, and

initpt [i - 1]

should contain

p_{i}

, for

i = 1, 2, \dots, n

and

j = 1, 2, \dots, L_{i}

; if you wish e05jb to use one of its preset initialization methods, you could choose one of two simple, three-point methods (see Figure 1). If the list generated by one of these methods contains infinite values, attempts are made to generate a safeguarded list using the function

subint (x, y)

(which is also used during the splitting procedure, and is described in [Splitting]). If infinite values persist, e05jb exits with

ifail = 3

. There is also the option to generate an initialization list with the aid of linesearches (by setting

iinit = 2

). Starting with the absolutely smallest point in the root box, linesearches are made along each coordinate. For each coordinate, the local minimizers found by the linesearches are put into the initialization list. If there were fewer than three minimizers, they are augmented by close-by values. The final preset initialization option (

iinit = 4

) generates a randomized list, so that independent multiple runs may be made if you suspect a global optimum has not been found. Each call to the initialization method (E05JAF not in this release) resets the initial-state vector for the Wichmann–Hill base-generator that is used. Depending on whether you set the optional parameter Repeatability to

'ON'

'OFF'

, the random state is initialized to give a repeatable or non-repeatable sequence. Then, a random integer between

3

and sdlist is selected, which is then used to determine the number of points to be generated in each coordinate; that is, numpts becomes a constant vector, set to this value. The components of list are then generated, from a uniform distribution on the root box if the box is finite, or else in a safeguarded fashion if any bound is infinite. The array

initpt

is set to point to the best point in list.

Given an initialization list (preset or otherwise), e05jb evaluates

F

x^{0}

, and sets the initial estimate of the global minimum,

x^{*}

, to

x^{0}

. Then, for

i = 1, 2, \dots, n

, the objective function

F

is evaluated at

L_{i} - 1

points that agree with

x^{*}

in all but the

i

th coordinate. We obtain pairs

({\hat{x}}^{j}, f_{i}^{j})

, for

j = 1, 2, \dots, L_{i}

, with:

x^{*} = {\hat{x}}^{j_{1}}

, say; with, for

j \neq j_{1}

{\hat{x}}_{k}^{j} = \{\begin{cases} x_{k}^{*} & if ​ k \neq i; \\ x_{k}^{j} & otherwise; \end{cases}

and with

f_{i}^{j} = F ({\hat{x}}^{j}) .

The point having the smallest function value is renamed

x^{*}

and the procedure is repeated with the next coordinate.

Once e05jb has a full set of initialization points and function values, it can generate an initial set of sub-boxes. Recall that the root box is

B [x, y] = [ℓ, u]

, having basepoint

x = x^{0}

. The opposite point

y

is a corner of

[ℓ, u]

farthest away from

x

, in some sense. The point

x

need not be a vertex of

[ℓ, u]

, and

y

is entitled to have infinite coordinates. We loop over each coordinate

i

, splitting the current box along coordinate

i

into

2 L_{i} - 2

2 L_{i} - 1

2 L_{i}

sub-intervals with exactly one of the

{\hat{x}}_{i}^{j}

as endpoints, depending on whether two, one or none of the

{\hat{x}}_{i}^{j}

are on the boundary. Thus, as well as splitting at

{\hat{x}}_{i}^{j}

, for

j = 1, 2, \dots, L_{i}

, we split at additional points

z_{i}^{j}

, for

j = 2, 3, \dots, L_{i}

. These additional

z_{i}^{j}

are such that

z_{i}^{j} = {\hat{x}}_{i}^{j - 1} + q^{m} ({\hat{x}}_{i}^{j} - {\hat{x}}_{i}^{j - 1}), j = 2, \dots, L_{i},

where

q

is the golden-section ratio

(\sqrt{5} - 1) / 2

, and the exponent

m

takes the value

1

2

, chosen so that the sub-box with the smaller function value gets the larger fraction of the interval. Each child sub-box gets as basepoint the point obtained from

x^{*}

by changing

x_{i}^{*}

to the

x_{i}^{j}

that is a boundary point of the corresponding

i

th coordinate interval; this new basepoint therefore has function value

f_{i}^{j}

. The opposite point is derived from

y

by changing

y_{i}

to the other end of that interval.

e05jb can now rank the coordinates based on an estimated variability of

F

. For each

i

we compute the union of the ranges of the quadratic interpolant through any three consecutive

{\hat{x}}_{i}^{j}

, taking the difference between the upper and lower bounds obtained as a measure of the variability of

F

in coordinate

i

. A vector

π

is populated in such a way that coordinate

i

has the

π_{i}

th highest estimated variability. For tiebreaks, when the

x^{*}

obtained after splitting coordinate

i

belongs to two sub-boxes, the one that contains the minimizer of the quadratic models is designated the current sub-box for coordinate

i + 1

Boxes are assigned levels in the following manner. The root box is given level

1

. When a sub-box of level

s

is split, the child with the smaller fraction of the golden-section split receives level

s + 2

; all other children receive level

s + 1

. The box with the better function value is given the larger fraction of the splitting interval and the smaller level because then it is more likely to be split again more quickly. We see that after the initialization procedure the first level is empty and the non-split boxes have levels

2, \dots, n_{r} + 2

, so it is meaningful to choose

s_{\max}

much larger than

n_{r}

. Note that the internal structure of e05jb demands that

s_{\max}

be at least

n_{r} + 3

Examples of initializations in two dimensions are given in Figure 1. In both cases the initial point is

x^{0} = (ℓ + u) / 2

; on the left the initialization points are

x^{1} = ℓ, x^{2} = (ℓ + u) / 2, x^{3} = u,

while on the right the points are

x^{1} = (5 ℓ + u) / 6, x^{2} = (ℓ + u) / 2, x^{3} = (ℓ + 5 u) / 6 .

In Figure 1, basepoints and levels after initialization are displayed. Note that these initialization lists correspond to

iinit = 0

and

iinit = 1

, respectively.

Figure 1: Examples of the initialization procedure

After initialization, a series of sweeps through levels is begun. A sweep is defined by three steps:

(i)	scan the list of non-split sub-boxes. Fill a record list $b$ according to $b_{s} = 0$ if there is no box at level $s$ , and with $b_{s}$ pointing to a sub-box with the lowest function value among all sub-boxes with level $s$ otherwise, for $0 < s < s_{\max}$ ;
(ii)	the sub-box with label $b_{s}$ is a candidate for splitting. If the sub-box is not to be split, according to the rules described in [Splitting], increase its level by $1$ and update $b_{s + 1}$ if necessary. If the sub-box is split, mark it so, insert its children into the list of sub-boxes, and update $b$ if any child with level $s^{'}$ yields a strict improvement of $F$ over those sub-boxes at level $s^{'}$ ;
(iii)	increment $s$ by $1$ . If $s = s_{\max}$ then displaying monitoring information and start a new sweep; else if $b_{s} = 0$ then repeat this step; else display monitoring information and go to the previous step.

Clearly, each sweep ends after at most

s_{\max} - 1

visits of the third step.

Splitting

Each sub-box is stored by e05jb as a set of information about the history of the sub-box: the label of its parent, a label identifying which child of the parent it is, etc. Whenever a sub-box

B [x, y]

of level

s < s_{\max}

is a candidate for splitting, as described in [Initialization and Sweeps], we recover

x

y

, and the number,

n_{j}

, of times coordinate

j

has been split in the history of

B

. Sub-box

B

could be split in one of two ways.

(i)	Splitting by rank If $s > 2 n_{r} (\min n_{j} + 1)$ , the box is always split. The splitting index is set to a coordinate $i$ such that $n_{i} = \min n_{j}$ .
(ii)	Splitting by expected gain If $s \leq 2 n_{r} (\min n_{j} + 1)$ , the sub-box could be split along a coordinate where a maximal gain in function value is expected. This gain is estimated according to a local separable quadratic model obtained by fitting to $2 n_{r} + 1$ function values. If the expected gain is too small the sub-box is not split at all, and its level is increased by $1$ .

Eventually, a sub-box that is not eligible for splitting by expected gain will reach level

2 n_{r} (\min n_{j} + 1) + 1

and then be split by rank, as long as

s_{\max}

is large enough. As

s_{\max} \to \infty

, the rule for splitting by rank ensures that each coordinate is split arbitrarily often.

Before describing the details of each splitting method, we introduce the procedure for correctly handling splitting at adaptive points and for dealing with unbounded intervals. Suppose we want to split the

i

th coordinate interval

▯ \{x_{i}, y_{i}\}

, where we define

▯ \{x_{i}, y_{i}\} = [\min (x_{i}, y_{i}), \max (x_{i}, y_{i})]

, for

x_{i} \in R

and

y_{i} \in \bar{R}

, and where

x

is the basepoint of the sub-box being considered. The descendants of the sub-box should shrink sufficiently fast, so we should not split too close to

x_{i}

. Moreover, if

y_{i}

is large we want the new splitting value to not be too large, so we force it to belong to some smaller interval

▯ \{ξ^{'}, ξ^{''}\}

, determined by

ξ^{''} = subint (x_{i}, y_{i}), ξ^{'} = x_{i} + (ξ^{''} - x_{i}) / 10,

where the function

subint

is defined by

subint (x, y) = \{\begin{cases} sign (y) & if ​ 1000 |x| < 1 ​ and ​ |y| > 1000; \\ 10 sign (y) |x| & if ​ 1000 |x| \geq 1 ​ and ​ |y| > 1000 |x|; \\ y & otherwise. \end{cases}

Splitting by rank

Consider a sub-box

B

with level

s > 2 n_{r} (\min n_{j} + 1)

. Although the sub-box has reached a high level, there is at least one coordinate along which it has not been split very often. Among the

i

such that

n_{i} = \min n_{j}

for

B

, select the splitting index to be the coordinate with the lowest

π_{i}

(and hence highest variability rank). ‘Splitting by rank’ refers to the ranking of the coordinates by

n_{i}

and

π_{i}

n_{i} = 0

, so that

B

has never been split along coordinate

i

, the splitting is done according to the initialization list and the adaptively chosen golden-section split points, as described in [Initialization and Sweeps]. Also as covered there, new basepoints and opposite points are generated. The children having the smaller fraction of the golden-section split (that is, those with larger function values) are given level

\min \{s + 2, s_{\max}\}

. All other children are given level

s + 1

Otherwise,

B

ranges between

x_{i}

and

y_{i}

in the

i

th coordinate direction. The splitting value is selected to be

z_{i} = x_{i} + 2 (subint (x_{i}, y_{i}) - x_{i}) / 3

; we are not attempting to split based on a large reduction in function value, merely in order to reduce the size of a large interval, so

z_{i}

may not be optimal. Sub-box

B

is split at

z_{i}

and the golden-section split point, producing three parts and requiring only one additional function evaluation, at the point

x^{'}

obtained from

x

by changing the

i

th coordinate to

z_{i}

. The child with the smaller fraction of the golden-section split is given level

\min \{s + 2, s_{\max}\}

, while the other two parts are given level

s + 1

. Basepoints are assigned as follows: the basepoint of the first child is taken to be

x

, and the basepoint of the second and third children is the point

x^{'}

. Opposite points are obtained by changing

y_{i}

to the other end of the

i

th coordinate-interval of the corresponding child.

Splitting by expected gain

When a sub-box

B

has level

s \leq 2 n_{r} (\min n_{j} + 1)

, we compute the optimal splitting index and splitting value from a local separable quadratic used as a simple local approximation of the objective function. To fit this curve, for each coordinate we need two additional points and their function values. Such data may be recoverable from the history of

B

: whenever the

i

th coordinate was split in the history of

B

, we obtained values that can be used for the current quadratic interpolation in coordinate

i

We loop over

i

; for each coordinate we pursue the history of

B

back to the root box, and we take the first two points and function values we find, since these are expected to be closest to the current basepoint

x

. If the current coordinate has not yet been split we use the initialization list. Then we generate a local separable model

e (ξ)

for

F (ξ)

by interpolation at

x

and the

2 n_{r}

additional points just collected:

e (ξ) = F (x) + \sum_{i = 1}^{n} e_{i} (ξ_{i}) .

We define the expected gain

{\hat{e}}_{i}

in function value when we evaluate at a new point obtained by changing coordinate

i

in the basepoint, for each

i

, based on two cases:

(i)

n_{i} = 0

. We compute the expected gain as

{\hat{e}}_{i} = \min_{1 \leq j \leq L_{i}} \{f_{i}^{j}\} - f_{i}^{p_{i}} .

Again, we split according to the initialization list, with the new basepoints and opposite points being as before.

(ii)

n_{i} > 0

. Now, the

i

th component of our sub-box ranges from

x_{i}

y_{i}

. Using the quadratic partial correction function

e_{i} (ξ_{i}) = α_{i} (ξ_{i} - x_{i}) + β_{i} {(ξ_{i} - x_{i})}^{2}

we can approximate the maximal gain expected when changing

x_{i}

only. We will choose the splitting value from

▯ \{ξ^{'}, ξ^{''}\}

. We compute

{\hat{e}}_{i} = \min_{ξ_{i} \in ▯ \{ξ^{'}, ξ^{''}\}} e_{i} (ξ_{i})

and call

z_{i}

the minimizer in

▯ \{ξ^{'}, ξ^{''}\}

If the expected best function value

f_{\exp}

satisfies

f_{\exp} = F (x) + \min_{1 \leq i \leq n} {\hat{e}}_{i} < f_{best},

(1)

where

f_{best}

is the current best function value (including those function values obtained by local optimization), we expect the sub-box to contain a better point and so we split it, using as splitting index the component with minimal

{\hat{e}}_{i}

. Equation (1) prevents wasting function calls by avoiding splitting sub-boxes whose basepoints have bad function values. These sub-boxes will eventually be split by rank anyway.

We now have a splitting index and a splitting value

z_{i}

. The sub-box is split at

z_{i}

as long as

z_{i} \neq y_{i}

, and at the golden-section split point; two or three children are produced. The larger fraction of the golden-section split receives level

s + 1

, while the smaller fraction receives level

\min \{s + 2, s_{\max}\}

. If it is the case that

z_{i} \neq y_{i}

and the third child is larger than the smaller of the two children from the golden-section split, the third child receives level

s + 1

. Otherwise it is given the level

\min \{s + 2, s_{\max}\}

. The basepoint of the first child is set to

x

, and the basepoint of the second (and third if it exists) is obtained by changing the

i

th coordinate of

x

z_{i}

. The opposite points are again derived by changing

y_{i}

to the other end of the

i

th coordinate interval of

B

If equation (1) does not hold, we expect no improvement. We do not split, and we increase the level of

B

1

Local Search

The local optimization algorithm used by e05jb uses linesearches along directions that are determined by minimizing quadratic models, all subject to bound constraints. Triples of vectors are computed using coordinate searches based on linesearches. These triples are used in triple search procedures to build local quadratic models for

F

. A trust-region-type approach to minimize these models is then carried out, and more information about the coordinate search and the triple search can be found in Huyer and Neumaier (1999).

The local search starts by looking for better points without being too local, by making a triple search using points found by a coordinate search. This yields a new point and function value, an approximation of the gradient of the objective, and an approximation of the Hessian of the objective. Then the quadratic model for

F

is minimized over a small box, with the solution to that minimization problem then being used as a linesearch direction to minimize the objective. A measure

r

is computed to quantify the predictive quality of the quadratic model.

The third stage is the checking of termination criteria. The local search will stop if more than

loclim

visits to this part of the local search have occurred, where

loclim

is the value of the optional parameter Local Searches Limit. If that is not the case, it will stop if the limit on function calls has been exceeded (see the description of the optional parameter Function Evaluations Limit). The final criterion checks if no improvement can be made to the function value, or whether the approximated gradient

g

is small, in the sense that

{|g|}^{T} \max (|x|, |x_{old}|) < loctol (f_{0} - f) .

The vector

x_{old}

is the best point at the start of the current loop in this iterative local-search procedure, the constant

loctol

is the value of the optional parameter Local Searches Tolerance,

f

is the objective value at

x

, and

f_{0}

is the smallest function value found by the initialization procedure.

Next, e05jb attempts to move away from the boundary, if any components of the current point lie there, using linesearches along the offending coordinates. Local searches are terminated if no improvement could be made.

The fifth stage carries out another triple search, but this time it does not use points from a coordinate search, rather points lying within the trust-region box are taken.

The final stage modifies the trust-region box to be bigger or smaller, depending on the quality of the quadratic model, minimizes the new quadratic model on that box, and does a linesearch in the direction of the minimizer. The value of

r

is updated using the new data, and then we go back to the third stage (checking of termination criteria).

The Hessians of the quadratic models generated by the local search may not be positive definite, so e05jb uses the general nonlinear optimizer e04vh to minimize the models.

Syntax

Parameters

Description

References

Error Indicators and Warnings

Accuracy

Parallelism and Performance

Further Comments

Example

Algorithmic Details

Initialization and Sweeps

Splitting

Splitting by rank

Splitting by expected gain

Local Search

See Also