Minimizing or Maximizing a Function

1

Scope of the Chapter

This chapter provides routines for solving various mathematical optimization problems by solvers based on local stopping criteria. The main classes of problems covered in this chapter are:

**Linear Programming (LP)**– dense and sparse;**Quadratic Programming (QP)**– convex and nonconvex, dense and sparse;**Nonlinear Programming (NLP)**– dense and sparse, based on active-set SQP methods or interior point methods (IPM);**Semidefinite Programming (SDP)**– both linear matrix inequalities (LMI) and bilinear matrix inequalities (BMI);**Derivative-free Optimization (DFO)**;**Least Squares (LSQ)**, data fitting – linear and nonlinear, constrained and unconstrained.

For a full overview of the functionality offered in this chapter, see Section 5 or the Chapter Contents (Chapter E04).

See also other chapters in the Library relevant to optimization:

- Chapter E05 contains routines to solve
**global optimization**problems; - Chapter H addresses problems arising in
**operational research**and focuses on**Mixed Integer Programming (MIP)**; - Chapters F07 and F08 include routines for linear algebra and in particular unconstrained linear least squares;
- Chapter E02 focuses on curve and surface fitting, in which linear data fitting in ${l}_{1}$ or ${l}_{\infty}$ norm might be of interest.

This introduction is only a brief guide to the subject of optimization. It discusses a classification of the optimization problems and presents an overview of the algorithms and their stopping criteria to help with the choice of a correct solver for a particular problem. Anyone with a difficult or protracted problem to solve will find it beneficial to consult a more detailed text, see Gill et al. (1981), Fletcher (1987) or Nocedal and Wright (2006). If you are unfamiliar with the mathematics of the subject you may find Sections 2.1, 2.2, 2.3, 2.6 and 3 a useful starting point.

2

Background to the Problems

2.1

Introduction to Mathematical Optimization

Mathematical Optimization, also known as Mathematical Programming, refers to the problem of finding values of the inputs from a given set so that a function (called the **objective function**) is minimized or maximized. The inputs are called **decision variables**, primal variables or just variables. The given set from which the decision variables are selected is referred to as a **feasible set** and might be defined as a domain where **constraints** expressed as functions of the decision variables hold certain values. Each point of the feasible set is called a **feasible point**.

A general mathematical formulation of such a problem might be written as

where $x$ denotes the decision variables, $f\left(x\right)$ the objective function and $\mathcal{F}$ the feasibility set. In this chapter we assume that $\mathcal{F}\subset {\mathbb{R}}^{n}$. Since **maximization** of the objective function $f\left(x\right)$ is equivalent to minimizing $-f\left(x\right)$, only minimization is considered further in the text. Some routines allow you to specify whether you are solving a minimization or maximization problem, carrying out the required transformation of the objective function in the latter case.

$$\begin{array}{ll}\mathrm{minimize}& f\left(x\right)\\ \text{subject to}& x\in \mathcal{F}\end{array}$$ |

A point ${x}^{*}$ is said to be a **local minimum** of a function $f$ if it is feasible (${x}^{*}\in \mathcal{F}$) and if $f\left(x\right)\ge f\left({x}^{*}\right)$ for all $x\in \mathcal{F}$ near ${x}^{*}$. A point ${x}^{*}$ is a **global minimum** if it is a local minimum and $f\left(x\right)\ge f\left({x}^{*}\right)$ for all feasible $x$. The solvers in this chapter are based on algorithms which seek only a local minimum, however, many problems (such as convex optimization problems) have only one local minimum. This is also the global minimum. In such cases the Chapter E04 solvers find the global minimum. See Chapter E05 for solvers which try to find a global solution even for nonconvex functions.

2.2

Classification of Optimization Problems

There is no single efficient solver for all optimization problems. Therefore it is important to choose a solver which matches the problem and any specific needs as closely as possible. A more generic solver might be applied, however the performance suffers in some cases, depending on the underlying algorithm.

There are various criteria to help to classify optimization problems into particular categories. The main criteria are as follows:

- Type of objective function;
- Type of constraints;
- Size of the problem;
- Smoothness of the data and available derivative information.

Each of the criteria is discussed below to give the necessary information to identify the class of the optimization problem. Section 2.5 presents the basic properties of the algorithms and Section 3 advises on the choice of particular routines in the chapter.

2.2.1

Types of objective functions

In general, if there is a structure in the problem the solver should benefit from it. For example, a solver for problems with the sum of squares objective should work better than when this objective is treated as a general nonlinear objective. Therefore it is important to recognize typical types of the objective functions.

An optimization problem which has **no objective** is equivalent to having a constant objective, i.e., $f\left(x\right)=0$. It is usually called a **feasible point problem**. The task is to then find any point which satisfies the constraints.

A **linear objective function** is a function which is linear in all variables and therefore can be represented as

where $c\in {\mathbb{R}}^{n}$. Scalar ${c}_{0}$ has no influence on the choice of decision variables $x$ and is usually omitted. It will not be used further in this text.

$$f\left(x\right)={c}^{\mathrm{T}}x+{c}_{0}$$ |

A **quadratic objective function** is an extension of a linear function with quadratic terms as follows:

Here $H$ is a real symmetric $n\times n$ matrix. In addition, if $H$ is positive semidefinite (all its eigenvalues are non-negative), the objective is **convex**.

$$f\left(x\right)=\frac{1}{2}{x}^{\mathrm{T}}Hx+{c}^{\mathrm{T}}x\text{.}$$ |

A general **nonlinear objective function** is any $f:{\mathbb{R}}^{n}\to \mathbb{R}$ without a special structure.

Special consideration is given to the objective function in the form of a **sum of squares** of functions, such as

where ${r}_{i}:{\mathbb{R}}^{n}\to \mathbb{R}$; often called residual functions. This form of the objective plays a key role in **data fitting** solved as a **least squares problem** as shown in Section 2.2.3.

$$f\left(x\right)=\sum _{i=1}^{m}{r}_{i}^{2}\left(x\right)$$ |

2.2.2

Types of constraints

Not all optimization problems have to have constraints. If there are no restrictions on the choice of $x$ except that $x\in \mathcal{F}={\mathbb{R}}^{n}$, the problem is called **unconstrained** and thus every point is a feasible point.

$${l}_{{x}_{i}}\le {x}_{i}\le {u}_{{x}_{i}}\text{,\hspace{1em}for}i=1,\dots ,n$$ |

$${l}_{x}\le x\le {u}_{x}$$ |

The same format of bounds is adopted to linear and nonlinear constraints in the whole chapter. Note that for the purpose of passing infinite bounds to the routines, all values above a certain threshold (typically ${10}^{20}$) are treated as $+\infty $.

$${l}_{B}\le Bx\le {u}_{B}$$ |

Although the bounds on ${x}_{i}$ could be included in the definition of linear constraints, we recommend you distinguish between them for reasons of computational efficiency as most of the solvers treat simple bounds explicitly.

A set of ${m}_{g}$ **nonlinear constraints** may be defined in terms of a nonlinear function $g:{\mathbb{R}}^{n}\to {\mathbb{R}}^{{m}_{g}}$ and the bounds ${l}_{g}$ and ${u}_{g}$ which follow the same format as simple bounds and linear constraints:

$${l}_{g}\le g\left(x\right)\le {u}_{g}\text{.}$$ |

Although the linear constraints could be included in the definition of nonlinear constraints, again we prefer to distinguish between them for reasons of computational efficiency.

A **matrix constraint** (or matrix inequality) is a constraint on eigenvalues of a matrix operator. More precisely, let ${\mathbb{S}}^{m}$ denote the space of real symmetric matrices $m$ by $m$ and let $\mathcal{A}$ be a matrix operator $\mathcal{A}:{\mathbb{R}}^{n}\to {\mathbb{S}}^{m}$, i.e., it assigns a symmetric matrix $\mathcal{A}\left(x\right)$ for each $x$. The matrix constraint can be expressed as

where the inequality $S\u2ab00$ for $S\in {\mathbb{S}}^{m}$ is meant in the eigenvalue sense, namely all eigenvalues of the matrix $S$ should be non-negative (the matrix should be positive semidefinite).

$$\mathcal{A}\left(x\right)\u2ab00$$ |

There are two types of matrix constraints allowed in the current mark of the Library. The first is **linear matrix inequality (LMI)** formulated as

and the second one, **bilinear matrix inequality (BMI)**, stated as

Here all matrices ${A}_{i}$, ${Q}_{ij}$ are given real symmetric matrices of the same dimension. Note that the latter type is in fact quadratic in $x$, nevertheless, it is referred to as bilinear for historical reasons.

$$\mathcal{A}\left(x\right)=\sum _{\mathit{i}=1}^{n}{x}_{\mathit{i}}{A}_{\mathit{i}}-{A}_{0}\u2ab00$$ |

$$\mathcal{A}\left(x\right)=\sum _{\mathit{i,j}=1}^{n}{x}_{\mathit{i}}{x}_{\mathit{j}}{Q}_{\mathit{i}\mathit{j}}+\sum _{\mathit{i}=1}^{n}{x}_{\mathit{i}}{A}_{\mathit{i}}-{A}_{0}\u2ab00\text{.}$$ |

2.2.3

Typical classes of optimization problems

Specific combinations of the types of the objective functions and constraints give rise to various classes of optimization problems. The common ones are presented below. It is always advisable to consider the closest formulation which covers your problem when choosing the solver. For more information see classical texts such as Dantzig (1963), Gill et al. (1981), Fletcher (1987), Nocedal and Wright (2006) or Chvátal (1983).

A **Linear Programming (LP)** problem is a problem with a linear objective function, linear constraints and simple bounds. It can be written as follows:

$$\begin{array}{ll}\underset{x\in {\mathbb{R}}^{n}}{\mathrm{minimize}}\phantom{\rule{0.25em}{0ex}}& {c}^{\mathrm{T}}x\\ \text{subject to}& {l}_{B}\le Bx\le {u}_{B}\\ & {l}_{x}\le x\le {u}_{x}\end{array}$$ |

$$\begin{array}{ll}\underset{x\in {\mathbb{R}}^{n}}{\mathrm{minimize}}\phantom{\rule{0.25em}{0ex}}& \frac{1}{2}{x}^{\mathrm{T}}Hx+{c}^{\mathrm{T}}x\\ \text{subject to}& {l}_{B}\le Bx\le {u}_{B}\\ & {l}_{x}\le x\le {u}_{x}\end{array}$$ |

$$\begin{array}{ll}\underset{x\in {\mathbb{R}}^{n}}{\mathrm{minimize}}\phantom{\rule{0.25em}{0ex}}& f\left(x\right)\\ \text{subject to}& {l}_{g}\le g\left(x\right)\le {u}_{g}\\ & {l}_{B}\le Bx\le {u}_{B}\\ & {l}_{x}\le x\le {u}_{x}\end{array}$$ |

$$\begin{array}{ll}\underset{x\in {\mathbb{R}}^{n}}{\mathrm{minimize}}\phantom{\rule{0.25em}{0ex}}& {c}^{\mathrm{T}}x\\ \text{subject to \hspace{1em}}& \sum _{\mathit{i}=1}^{n}{x}_{i}{A}_{i}^{k}-{A}_{0}^{k}\u2ab00\text{,\hspace{1em}}k=1,\dots ,{m}_{A}\\ & {l}_{B}\le Bx\le {u}_{B}\\ & {l}_{x}\le x\le {u}_{x}\end{array}$$ |

$$\begin{array}{ll}\underset{x\in {\mathbb{R}}^{n}}{\mathrm{minimize}}\phantom{\rule{0.25em}{0ex}}& \frac{1}{2}{x}^{\mathrm{T}}Hx+{c}^{\mathrm{T}}x\\ \text{subject to \hspace{1em}}& \sum _{\mathit{i,j}=1}^{n}{x}_{i}{x}_{j}{Q}_{ij}^{k}+\sum _{\mathit{i}=1}^{n}{x}_{i}{A}_{i}^{k}-{A}_{0}^{k}\u2ab00\text{,\hspace{1em}}k=1,\dots ,{m}_{A}\\ & {l}_{B}\le Bx\le {u}_{B}\\ & {l}_{x}\le x\le {u}_{x}\end{array}$$ |

A **least squares (LSQ)** problem is a problem where the objective function in the form of sum of squares is minimized subject to usual constraints. If the residual functions ${r}_{i}\left(x\right)$ are linear or nonlinear, the problem is known as **linear** or **nonlinear least squares**, respectively. Not all types of the constraints need to be present which brings up special cases of unconstrained, bound-constrained or linearly-constrained least squares problems as in NLP .

$$\begin{array}{ll}\underset{x\in {\mathbb{R}}^{n}}{\mathrm{minimize}}\phantom{\rule{0.25em}{0ex}}& \sum _{i=1}^{m}{r}_{i}^{2}\left(x\right)\\ \text{subject to}& {l}_{g}\le g\left(x\right)\le {u}_{g}\\ & {l}_{B}\le Bx\le {u}_{B}\\ & {l}_{x}\le x\le {u}_{x}\end{array}$$ |

This form of the problem is very common in **data fitting** as demonstrated on the following example. Let us consider a process that is observed at times ${t}_{i}$ and measured with results ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$. Furthermore, the process is assumed to behave according to a model $\varphi \left(t;x\right)$ where $x$ are parameters of the model. Given the fact that the measurements might be inaccurate and the process might not exactly follow the model, it is beneficial to find model parameters $x$ so that the error of the fit of the model to the measurements is minimized. This can be formulated as an optimization problem in which $x$ are decision variables and the objective function is the sum of squared errors of the fit at each individual measurement, thus:

$$\begin{array}{llll}\underset{x\in {\mathbb{R}}^{n}}{\mathrm{minimize}}\phantom{\rule{0.25em}{0ex}}& \sum _{i=1}^{m}{r}_{i}^{2}\left(x\right)& \text{where}& {r}_{i}\left(x\right)=\varphi \left({t}_{i};x\right)-{y}_{i}\end{array}$$ |

2.2.4

Problem size, dense and sparse problems

The size of the optimization problem plays an important role in the choice of the solver. The size is usually understood to be the number of variables $n$ and the number (and the type) of the constraints. Depending on the size of the problem we talk about small-scale, medium-scale or large-scale problems.

It is often more practical to look at the data and its structure rather than just the size of the problem. Typically in a large-scale problem not all variables interact with everything else. It is natural that only a small portion of the constraints (if any) involves all variables and the majority of the constraints depends only on small different subsets of the variables. This creates many explicit zeros in the data representation which it is beneficial to capture and pass to the solver. In such a case the problem is referred to as **sparse**. The data representation usually has the form of a sparse matrix which defines the linear constraint matrix $B$, Jacobian matrix of the nonlinear constraints ${g}_{i}$ or the Hessian of the objective $H$. Common sparse matrix formats are used, such as coordinate storage (CS) and compressed column storage (CCS) (see Section 2.1 in the F11 Chapter Introduction).

The counterpart to a sparse problem is a **dense** problem in which the matrices are stored in general full format and no structure is assumed or exploited. Whereas passing a dense problem to a sparse solver presents typically only a small overhead, calling a dense solver on a large-scale sparse problem is ill-advised; it leads to a significant performance degradation and memory overuse.

2.2.5

Derivative information, smoothness, noise and derivative-free optimization (DFO)

Most of the classical optimization algorithms rely heavily on derivative information. It plays a key role in necessary and sufficient conditions (see Section 2.4) and in the computation of the search direction at each iteration (see Section 2.5). Therefore it is important that accurate derivatives of the nonlinear objective and nonlinear constraints are provided whenever possible.

Unless stated otherwise, it is assumed that the nonlinear functions are sufficiently smooth. The solvers will usually solve optimization problems even if there are isolated discontinuities away from the solution, however you should always consider whether an alternative smooth representation of the problem exists. A typical example is an absolute value $\left|{x}_{i}\right|$ which does not have a first derivative for ${x}_{i}=0$. Nevertheless, if the model allows it can be transformed as

which avoids the discontinuity of the first derivative. If many discontinuities are present, alternative methods need to be applied such as e04cbf or stochastic algorithms in Chapter E05, e05saf or e05sbf.

$$\begin{array}{llll}{x}_{i}={x}_{i}^{+}-{x}_{i}^{-}\text{,}& \left|{x}_{i}\right|={x}_{i}^{+}+{x}_{i}^{-}\text{,}& \text{where}& {x}_{i}^{+}\text{,}{x}_{i}^{-}\ge 0\end{array}$$ |

The vector of first partial derivatives of a function is called the **gradient vector**, i.e.,

the matrix of second partial derivatives is termed the **Hessian matrix**, i.e.,

and the matrix of first partial derivatives of the vector-valued function $f:{\mathbb{R}}^{n}\to {\mathbb{R}}^{m}$ is known as the **Jacobian matrix**:

$$\nabla f\left(x\right)={\left[\frac{\partial f\left(x\right)}{\partial {x}_{1}},\frac{\partial f\left(x\right)}{\partial {x}_{2}},\dots ,\frac{\partial f\left(x\right)}{\partial {x}_{n}}\right]}^{\mathrm{T}}\text{,}$$ |

$${\nabla}^{2}f\left(x\right)={\left[\frac{{\partial}^{2}f\left(x\right)}{\partial {x}_{i}\partial {x}_{j}}\right]}_{i,j=1,\dots ,n}$$ |

$$J\left(x\right)={\left[\frac{\partial {f}_{i}\left(x\right)}{\partial {x}_{j}}\right]}_{i=1,\dots ,m,j=1,\dots ,n}\text{.}$$ |

If the function is smooth and the derivative is unavailable, it is possible to approximate it by **finite differences**, a change in function values in response to small perturbations of the variables. Many routines in the Library estimate missing elements of the gradients automatically this way. The choice of the size of the perturbations strongly affects the quality of the approximation. Too small perturbations might spoil the approximation due to the cancellation errors in floating-point arithmetic and too big reduce the match of the finite differences and the derivative (see e04xaf/e04xaa for optimal balance of the factors). In addition, finite differences are very sensitive to the accuracy of $f\left(x\right)$. They might be unreliable or fail completely if the function evaluation is inaccurate or noisy such as when $f\left(x\right)$ is a result of a stochastic simulation or an approximate solution of a PDE.

2.2.6

Minimization subject to bounds on the objective function

In all of the above problem categories it is assumed that

where $a=-\infty $ and $b=+\infty $. Problems in which $a$ and/or $b$ are finite can be solved by adding an extra constraint of the appropriate type (i.e., linear or nonlinear) depending on the form of $f\left(x\right)$. Further advice is given in Section 3.7.

$$a\le f\left(x\right)\le b$$ |

2.2.7

Multi-objective optimization

Sometimes a problem may have two or more objective functions which are to be optimized at the same time. Such problems are called **multi-objective**, multi-criteria or multi-attribute optimization. If the constraints are linear and the objectives are all linear then the terminology **goal programming** is also used.

Although there is no routine dealing with this type of problems explicitly in this mark of the Library, techniques used in this chapter and in Chapter E05 may be employed to address such problems, see Section 2.5.5.

2.3

Geometric Representation

To illustrate the nature of optimization problems it is useful to consider the following example:

(This function is used as the example function in the documentation for the unconstrained routines.)
**Figure 1**

$$f\left(x\right)={e}^{{x}_{1}}\left(4{x}_{1}^{2}+2{x}_{2}^{2}+4{x}_{1}{x}_{2}+2{x}_{2}+1\right)\text{.}$$ |

Figure 1 is a contour diagram of $f\left(x\right)$. The contours labelled ${F}_{0},{F}_{1},\dots ,{F}_{4}$ are isovalue contours, or lines along which the function $f\left(x\right)$ takes specific constant values. The point ${x}^{*}={\left(\frac{1}{2},-1\right)}^{\mathrm{T}}$ is a **local unconstrained minimum**, that is, the value of $f\left({x}^{*}\right)$ ($\text{}=0$) is less than at all the neighbouring points. A function may have several such minima. The point ${x}_{s}$ is said to be a **saddle point** because it is a minimum along the line AB, but a maximum along CD.

If we add the constraint ${x}_{1}\ge 0$ (a simple bound) to the problem of minimizing $f\left(x\right)$, the solution remains unaltered. In Figure 1 this constraint is represented by the straight line passing through ${x}_{1}=0$, and the shading on the line indicates the unacceptable region (i.e., ${x}_{1}<0$).

If we add the nonlinear constraint ${g}_{1}\left(x\right):{x}_{1}+{x}_{2}-{x}_{1}{x}_{2}-\frac{3}{2}\ge 0$, represented by the curved shaded line in Figure 1, then ${x}^{*}$ is not a feasible point because ${g}_{1}\left({x}^{*}\right)<0$. The solution of the new constrained problem is ${x}_{b}\simeq {\left(1.1825,-1.7397\right)}^{\mathrm{T}}$, the feasible point with the smallest function value (where $f\left({x}_{b}\right)\simeq 3.0607$).

2.4

Sufficient Conditions for a Solution

All nonlinear functions will be assumed to have continuous second derivatives in the neighbourhood of the solution.

2.4.1

Unconstrained minimization

The following conditions are sufficient for the point ${x}^{*}$ to be an unconstrained local minimum of $f\left(x\right)$:

where $\Vert \xb7\Vert $ denotes the Euclidean norm.

(i) | $\Vert \nabla f\left({x}^{*}\right)\Vert =0$ and |

(ii) | ${\nabla}^{2}f\left({x}^{*}\right)$ is positive definite, |

2.4.2

Minimization subject to bounds on the variables

At the solution of a bounds-constrained problem, variables which are not on their bounds are termed **free variables**. If it is known in advance which variables are on their bounds at the solution, the problem can be solved as an unconstrained problem in just the free variables; thus, the sufficient conditions for a solution are similar to those for the unconstrained case, applied only to the free variables.

Sufficient conditions for a feasible point ${x}^{*}$ to be the solution of a bounds-constrained problem are as follows:

where $\stackrel{-}{g}\left(x\right)$ is the gradient of $f\left(x\right)$ with respect to the free variables, and $\stackrel{-}{G}\left(x\right)$ is the Hessian matrix of $f\left(x\right)$ with respect to the free variables. The extra condition (iii) ensures that $f\left(x\right)$ cannot be reduced by moving off one or more of the bounds.

(i) | $\Vert \stackrel{-}{g}\left({x}^{*}\right)\Vert =0$; and |

(ii) | $\stackrel{-}{G}\left({x}^{*}\right)$ is positive definite; and |

(iii) | $\frac{\partial}{\partial {x}_{j}}f\left({x}^{*}\right)<0,{x}_{j}={u}_{j}$; $\frac{\partial}{\partial {x}_{j}}f\left({x}^{*}\right)>0,{x}_{j}={l}_{j}$, |

2.4.3

Linearly-constrained minimization

For the sake of simplicity, the following description does not include a specific treatment of bounds or range constraints, since the results for general linear inequality constraints can be applied directly to these cases.

At a solution ${x}^{*}$, of a linearly-constrained problem, the constraints which hold as equalities are called the **active** or **binding** constraints. Assume that there are $t$ active constraints at the solution ${x}^{*}$, and let $\hat{A}$ denote the matrix whose columns are the columns of $A$ corresponding to the active constraints, with $\hat{b}$ the vector similarly obtained from $b$; then

The matrix $Z$ is defined as an $n\times \left(n-t\right)$ matrix satisfying:

The columns of $Z$ form an orthogonal basis for the set of vectors orthogonal to the columns of $\hat{A}$.

$${\hat{A}}^{\mathrm{T}}{x}^{*}=\hat{b}\text{.}$$ |

$$\begin{array}{l}{\hat{A}}^{\mathrm{T}}Z=0\text{;}\\ {Z}^{\mathrm{T}}Z=I\text{.}\end{array}$$ |

Define

- ${g}_{Z}\left(x\right)={Z}^{\mathrm{T}}\nabla f\left(x\right)$, the
**projected gradient vector**of $f\left(x\right)$; - ${G}_{Z}\left(x\right)={Z}^{\mathrm{T}}{\nabla}^{2}f\left(x\right)Z$, the
**projected Hessian matrix**of $f\left(x\right)$.

At the solution of a linearly-constrained problem, the projected gradient vector must be zero, which implies that the gradient vector $\nabla f\left({x}^{*}\right)$ can be written as a linear combination of the columns of $\hat{A}$, i.e., $\nabla f\left({x}^{*}\right)={\displaystyle \sum _{i=1}^{t}}{\lambda}_{i}^{*}{\hat{a}}_{i}=\hat{A}{\lambda}^{*}$. The scalar ${\lambda}_{i}^{*}$ is defined as the **Lagrange multiplier** corresponding to the $i$th active constraint. A simple interpretation of the $i$th Lagrange multiplier is that it gives the gradient of $f\left(x\right)$ along the $i$th active constraint normal; a convenient definition of the Lagrange multiplier vector (although not a recommended method for computation) is:

Sufficient conditions for ${x}^{*}$ to be the solution of a linearly-constrained problem are:

$${\lambda}^{*}={\left({\hat{A}}^{\mathrm{T}}\hat{A}\right)}^{-1}{\hat{A}}^{\mathrm{T}}\nabla f\left({x}^{*}\right)\text{.}$$ |

(i) | ${x}^{*}$ is feasible, and ${\hat{A}}^{\mathrm{T}}{x}^{*}=\hat{b}$; and |

(ii) | $\Vert {g}_{Z}\left({x}^{*}\right)\Vert =0$, or equivalently, $\nabla f\left({x}^{*}\right)=\hat{A}{\lambda}^{*}$; and |

(iii) | ${G}_{Z}\left({x}^{*}\right)$ is positive definite; and |

(iv) | ${\lambda}_{i}^{*}>0$ if ${\lambda}_{i}^{*}$ corresponds to a constraint ${\hat{a}}_{i}^{\mathrm{T}}{x}^{*}\ge {\hat{b}}_{i}$;
${\lambda}_{i}^{*}<0$ if ${\lambda}_{i}^{*}$ corresponds to a constraint ${\hat{a}}_{i}^{\mathrm{T}}{x}^{*}\le {\hat{b}}_{i}$.
The sign of ${\lambda}_{i}^{*}$ is immaterial for equality constraints, which by definition are always active. |

2.4.4

Nonlinearly-constrained minimization

For nonlinearly-constrained problems, much of the terminology is defined exactly as in the linearly-constrained case. To simplify the notation, let us assume that all nonlinear constraints are in the form $c\left(x\right)\ge 0$. The set of active constraints at $x$ again means the set of constraints that hold as equalities at $x$, with corresponding definitions of $\hat{c}$ and $\hat{A}$: the vector $\hat{c}\left(x\right)$ contains the active constraint functions, and the columns of $\hat{A}\left(x\right)$ are the gradient vectors of the active constraints. As before, $Z$ is defined in terms of $\hat{A}\left(x\right)$ as a matrix such that:

where the dependence on $x$ has been suppressed for compactness.

$$\begin{array}{l}{\hat{A}}^{\mathrm{T}}Z=0\text{;}\\ {Z}^{\mathrm{T}}Z=I\end{array}$$ |

The projected gradient vector ${g}_{Z}\left(x\right)$ is the vector ${Z}^{\mathrm{T}}\nabla f\left(x\right)$. At the solution ${x}^{*}$ of a nonlinearly-constrained problem, the projected gradient must be zero, which implies the existence of Lagrange multipliers corresponding to the active constraints, i.e., $\nabla f\left({x}^{*}\right)=\hat{A}\left({x}^{*}\right){\lambda}^{*}$.

The **Lagrangian function** is given by:

We define ${g}_{L}\left(x\right)$ as the gradient of the Lagrangian function; ${G}_{L}\left(x\right)$ as its Hessian matrix, and ${\hat{G}}_{L}\left(x\right)$ as its projected Hessian matrix, i.e., ${\hat{G}}_{L}={Z}^{\mathrm{T}}{G}_{L}Z$.

$$L\left(x,\lambda \right)=f\left(x\right)-{\lambda}^{\mathrm{T}}\hat{c}\left(x\right)\text{.}$$ |

Sufficient conditions for ${x}^{*}$ to be the solution of a nonlinearly-constrained problem are:

(i) | ${x}^{*}$ is feasible, and $\hat{c}\left({x}^{*}\right)=0$; and |

(ii) | $\Vert {g}_{Z}\left({x}^{*}\right)\Vert =0$, or, equivalently, $\nabla f\left({x}^{*}\right)=\hat{A}\left({x}^{*}\right){\lambda}^{*}$; and |

(iii) | ${\hat{G}}_{L}\left({x}^{*}\right)$ is positive definite; and |

(iv) | ${\lambda}_{i}^{*}>0$ if ${\lambda}_{i}^{*}$ corresponds to a constraint of the form ${\hat{c}}_{i}\ge 0$.
The sign of ${\lambda}_{i}^{*}$ is immaterial for equality constraints, which by definition are always active. |

Note that condition (ii) implies that the projected gradient of the Lagrangian function must also be zero at ${x}^{*}$, since the application of ${Z}^{\mathrm{T}}$ annihilates the matrix $\hat{A}\left({x}^{*}\right)$.

2.5

Background to Optimization Methods

All the algorithms contained in this chapter generate an iterative sequence $\left\{{x}^{\left(k\right)}\right\}$ that converges to the solution ${x}^{*}$ in the limit, except for some special problem categories (i.e., linear and quadratic programming). To terminate computation of the sequence, a convergence test is performed to determine whether the current estimate of the solution is an adequate approximation. The convergence tests are discussed in Section 2.7.

Most of the methods construct a sequence $\left\{{x}^{\left(k\right)}\right\}$ satisfying:

where the vector ${p}^{\left(k\right)}$ is termed the **direction of search**, and ${\alpha}^{\left(k\right)}$ is the **steplength**. The steplength ${\alpha}^{\left(k\right)}$ is chosen so that $f\left({x}^{\left(k+1\right)}\right)<f\left({x}^{\left(k\right)}\right)$ and is computed using one of the techniques for one-dimensional optimization referred to in Section 2.5.1.

$${x}^{\left(k+1\right)}={x}^{\left(k\right)}+{\alpha}^{\left(k\right)}{p}^{\left(k\right)}\text{,}$$ |

2.5.1

One-dimensional optimization

The Library contains two special routines for minimizing a function of a single variable. Both routines are based on safeguarded polynomial approximation. One routine requires function evaluations only and fits a quadratic polynomial whilst the other requires function and gradient evaluations and fits a cubic polynomial. See Section 4.1 of Gill et al. (1981).

2.5.2

Methods for unconstrained optimization

The distinctions between methods arise primarily from the need to use varying levels of information about derivatives of $f\left(x\right)$ in defining the search direction. We describe three basic approaches to unconstrained problems, which may be extended to other problem categories. Since a full description of the methods would fill several volumes, the discussion here can do little more than allude to the processes involved, and direct you to other sources for a full explanation.

(a) | Newton-type Methods (Modified Newton Methods)
Newton-type methods use the Hessian matrix ${\nabla}^{2}f\left({x}^{\left(k\right)}\right)$, or its finite difference approximation , to define the search direction. The routines in the Library either require a subroutine that computes the elements of the Hessian directly, or they approximate them by finite differences.
Newton-type methods are the most powerful methods available for general problems and will find the minimum of a quadratic function in one iteration. See Sections 4.4 and 4.5.1 of Gill et al. (1981). |

(b) | Quasi-Newton Methods
Quasi-Newton methods approximate the Hessian ${\nabla}^{2}f\left({x}^{\left(k\right)}\right)$ by a matrix ${B}^{\left(k\right)}$ which is modified at each iteration to include information obtained about the curvature of $f$ along the current search direction ${p}^{\left(k\right)}$. Although not as robust as Newton-type methods, quasi-Newton methods can be more efficient because the Hessian is not computed directly, or approximated by finite differences. Quasi-Newton methods minimize a quadratic function in $n$ iterations, where $n$ is the number of variables. See Section 4.5.2 of Gill et al. (1981). |

(c) | Conjugate-gradient Methods
Unlike Newton-type and quasi-Newton methods, conjugate-gradient methods do not require the storage of an $n$ by $n$ matrix and so are ideally suited to solve large problems. Conjugate-gradient type methods are not usually as reliable or efficient as Newton-type, or quasi-Newton methods. See Section 4.8.3 of Gill et al. (1981). |

2.5.3

Methods for nonlinear least squares problems

These methods are similar to those for general nonlinear optimization, but exploit the special structure of the Hessian matrix to give improved computational efficiency.

Since

the Hessian matrix is of the form

where $J\left(x\right)$ is the Jacobian matrix of $r\left(x\right)$.

$$f\left(x\right)=\sum _{i=1}^{m}{r}_{i}^{2}\left(x\right)$$ |

$${\nabla}^{2}f\left(x\right)=2\left(J{\left(x\right)}^{\mathrm{T}}J\left(x\right)+\sum _{i=1}^{m}{r}_{i}\left(x\right){\nabla}^{2}{r}_{i}\left(x\right)\right)\text{,}$$ |

In the neighbourhood of the solution, $\Vert r\left(x\right)\Vert $ is often small compared to $\Vert J{\left(x\right)}^{\mathrm{T}}J\left(x\right)\Vert $ (for example, when $r\left(x\right)$ represents the goodness-of-fit of a nonlinear model to observed data). In such cases, $2J{\left(x\right)}^{\mathrm{T}}J\left(x\right)$ may be an adequate approximation to ${\nabla}^{2}f\left(x\right)$, thereby avoiding the need to compute or approximate second derivatives of $\left\{{r}_{i}\left(x\right)\right\}$. See Section 4.7 of Gill et al. (1981).

2.5.4

Methods for handling constraints

There are two main approaches for handling constraints in optimization algorithms – the **active-set sequential quadratic programming** method (or just SQP) and the **interior point method** (IPM). It is important to understand their very distinct features as both algorithms complement each other. The easiest method of comparison is to look at how the inequality constraints are treated and how the solver approaches the optimal solution (the progress of the KKT optimality measures: optimality, feasibility, complementarity).

Inequality constraints are the hard part of the optimization because of their ‘twofold nature’. If the optimal solution strictly satisfies the inequality, i.e., the optimal point is in the interior of the constraint, the inequality constraint does not influence the result and could be removed from the model. On the other hand, if the inequality is satisfied as an equality (is active at the solution), the constraint must be present and could be treated as an equality from the very beginning. This is expressed by the complementarity in KKT conditions.

Solvers, based on the active-set method, solve at each iteration a quadratic approximation of the original problem; they try to estimate which constraints need to be kept (are active) and which can be ignored. A practical consequence is that the algorithm partly ‘walks along the boundary’ of the feasible region given by the constraints. The iterates are thus feasible early on with regard to all linear constraints (and a local linearization of the nonlinear constraints) which is preserved through the iterations. The complementarity is satisfied by default, and once the active set is determined correctly and optimality is within the tolerance, the solver finishes. The number of iterations might be high but each is relatively cheap. See Chapter 6 of Gill et al. (1981) for further details.

In contrast, an interior point method generates iterations that avoid the boundary defined by the inequality constraints. As the solver progresses the iterates are allowed to get closer and closer to the boundary and converge to the optimal solution which might lie on the boundary. From the practical point of view, IPM typically requires only tens of iterations. Each iteration consists of solving a large linear system of equations taking into account all variables and constraints, so each iteration is fairly expensive. All three optimality measures are reduced simultaneously.

In many cases it is difficult to predict which of the algorithms will behave better on a particular problem, however, some initial guidance can be given in the following table:

IPM advantages |
SQP advantages |

Can exploit second derivatives and its structure Efficient on unconstrained or loosely constrained problems Efficient also for (both convex and nonconvex) quadratic problems (QP) Better use of multi-core architecture (SMP library only) New interface, easier to use |
Stay feasible with regard to linear constraints through most of the iterations Very efficient for highly constrained problems Better results on pathological problems in our experience Generally requires less function evaluations (efficient for problems with expensive function evaluations) Requires first derivatives but can work only with function values Can capitalize on a good initial point Allows warm starts (good for a sequence of similar problems) Infeasibility detection |

Unless some of the specific features are required which are offered only by one algorithm, the initial decision should be based on the availability of the derivatives of the problem and the number of constraints (for example, expressed as a ratio between the numbers of variables and the sum of the number of linear and nonlinear constraints). Readiness of exact second derivatives is a clear advantage for IPM so unless the number of constraints is close to the number of variables, IPM will probably work better. Similarly, if a large-scale problem has relatively few constraints (e.g., less than 40%) IPM might be more successful, especially as the problem gets bigger. On the other hand, if no derivatives are available, either the SQP or a specialized algorithm from the Library (see Derivative Free Optimization, Section 2.2.5) needs to be used. With more and more constraints SQP might be faster. For problems which do not fall in either of the categories, it is not easy to anticipate which solver will work better and some experimentation might be required.

2.5.5

Methods for handling multi-objective optimization

Suppose we have objective functions ${f}_{i}\left(x\right)$, $i>1$, all of which we need to minimize at the same time. There are two main approaches to this problem:

(a) | Combine the individual objectives into one composite objective. Typically this might be a weighted sum of the objectives, e.g.,
Here you choose the weights to express the relative importance of the corresponding objective. Ideally each of the ${f}_{i}\left(x\right)$ should be of comparable size at a solution. |
||

(b) | Order the objectives in order of importance. Suppose ${f}_{\mathit{i}}$ are ordered such that ${f}_{\mathit{i}}\left(x\right)$ is more important than ${f}_{\mathit{i}+1}\left(x\right)$, for $\mathit{i}=1,2,\dots ,n-1$. Then in the lexicographical approach to multi-objective optimization a sequence of subproblems are solved. Firstly solve the problem for objective function ${f}_{1}\left(x\right)$ and denote by ${r}_{1}$ the value of this minimum. If $\left(\mathit{i}-1\right)$ subproblems have been solved with results ${r}_{\mathit{i}-1}$ then subproblem $\mathit{i}$ becomes $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({f}_{\mathit{i}}\left(x\right)\right)$ subject to ${r}_{k}\le {f}_{k}\left(x\right)\le {r}_{k}$, for $\mathit{k}=1,2,\dots ,i-1$ plus the other constraints. |

Clearly the bounds on ${f}_{k}$ might be relaxed at your discretion.

In general, if NAG routines from the Chapter E04 are used then only local minima are found. This means that a better solution to an individual objective might be found without worsening the optimal solutions to the other objectives. Ideally you seek a Pareto solution; one in which an improvement in one objective can only be achieved by a worsening of another objective.

To obtain a Pareto solution routines from Chapter E05 might be used or, alternatively, a pragmatic attempt to derive a global minimum might be tried (see e05ucf). In this approach a variety of different minima are computed for each subproblem by starting from a range of different starting points. The best solution achieved is taken to be the global minimum. The more starting points chosen the greater confidence you might have in the computed global minimum.

2.6

Scaling

Scaling (in a broadly defined sense) often has a significant influence on the performance of optimization methods.

Since convergence tolerances and other criteria are necessarily based on an implicit definition of ‘small’ and ‘large’, problems with unusual or unbalanced scaling may cause difficulties for some algorithms.

Although there are currently no user-callable scaling routines in the Library, scaling can be performed automatically in routines which solve sparse LP, QP or NLP problems and in some dense solver routines. Such routines have an optional parameter ‘Scale Option’ which can be set by the user; see individual routine documents for details.

The following sections present some general comments on problem scaling.

2.6.1

Transformation of variables

One method of scaling is to transform the variables from their original representation, which may reflect the physical nature of the problem, to variables that have certain desirable properties in terms of optimization. It is generally helpful for the following conditions to be satisfied:

(i) | the variables are all of similar magnitude in the region of interest; |

(ii) | a fixed change in any of the variables results in similar changes in $f\left(x\right)$. Ideally, a unit change in any variable produces a unit change in $f\left(x\right)$; |

(iii) | the variables are transformed so as to avoid cancellation error in the evaluation of $f\left(x\right)$. |

Normally, you should restrict yourself to linear transformations of variables, although occasionally nonlinear transformations are possible. The most common such transformation (and often the most appropriate) is of the form

where $D$ is a diagonal matrix with constant coefficients. Our experience suggests that more use should be made of the transformation

where $v$ is a constant vector.

$${x}_{\mathrm{new}}=D{x}_{\mathrm{old}}\text{,}$$ |

$${x}_{\mathrm{new}}=D{x}_{\mathrm{old}}+v\text{,}$$ |

Consider, for example, a problem in which the variable ${x}_{3}$ represents the position of the peak of a Gaussian curve to be fitted to data for which the extreme values are $150$ and $170$; therefore ${x}_{3}$ is known to lie in the range $150$–$170$. One possible scaling would be to define a new variable ${\stackrel{-}{x}}_{3}$, given by

A better transformation, however, is given by defining ${\stackrel{-}{x}}_{3}$ as

Frequently, an improvement in the accuracy of evaluation of $f\left(x\right)$ can result if the variables are scaled before the routines to evaluate $f\left(x\right)$ are coded. For instance, in the above problem just mentioned of Gaussian curve-fitting, ${x}_{3}$ may always occur in terms of the form $\left({x}_{3}-{x}_{m}\right)$, where ${x}_{m}$ is a constant representing the mean peak position.

$${\stackrel{-}{x}}_{3}=\frac{{x}_{3}}{170}\text{.}$$ |

$${\stackrel{-}{x}}_{3}=\frac{{x}_{3}-160}{10}\text{.}$$ |

2.6.2

Scaling the objective function

The objective function has already been mentioned in the discussion of scaling the variables. The solution of a given problem is unaltered if $f\left(x\right)$ is multiplied by a positive constant, or if a constant value is added to $f\left(x\right)$. It is generally preferable for the objective function to be of the order of unity in the region of interest; thus, if in the original formulation $f\left(x\right)$ is always of the order of ${10}^{+5}$ (say), then the value of $f\left(x\right)$ should be multiplied by ${10}^{-5}$ when evaluating the function within an optimization routine. If a constant is added or subtracted in the computation of $f\left(x\right)$, usually it should be omitted, i.e., it is better to formulate $f\left(x\right)$ as ${x}_{1}^{2}+{x}_{2}^{2}$ rather than as ${x}_{1}^{2}+{x}_{2}^{2}+1000$ or even ${x}_{1}^{2}+{x}_{2}^{2}+1$. The inclusion of such a constant in the calculation of $f\left(x\right)$ can result in a loss of significant figures.

2.6.3

Scaling the constraints

A ‘well scaled’ set of constraints has two main properties. Firstly, each constraint should be well-conditioned with respect to perturbations of the variables. Secondly, the constraints should be balanced with respect to each other, i.e., all the constraints should have ‘equal weight’ in the solution process.

The solution of a linearly- or nonlinearly-constrained problem is unaltered if the $i$th constraint is multiplied by a positive weight ${w}_{i}$. At the approximation of the solution determined by an active-set solver, any active linear constraints will (in general) be satisfied ‘exactly’ (i.e., to within the tolerance defined by machine precision) if they have been properly scaled. This is in contrast to any active nonlinear constraints, which will not (in general) be satisfied ‘exactly’ but will have ‘small’ values (for example, ${\hat{g}}_{1}\left({x}^{*}\right)={10}^{-8}$, ${\hat{g}}_{2}\left({x}^{*}\right)={-10}^{-6}$, and so on). In general, this discrepancy will be minimized if the constraints are weighted so that a unit change in $x$ produces a similar change in each constraint.

A second reason for introducing weights is related to the effect of the size of the constraints on the Lagrange multiplier estimates and, consequently, on the active-set strategy. This means that different sets of weights may cause an algorithm to produce different sequences of iterates. Additional discussion is given in Gill et al. (1981).

2.7

Analysis of Computed Results

2.7.1

Convergence criteria

The convergence criteria inevitably vary from routine to routine, since in some cases more information is available to be checked (for example, is the Hessian matrix positive definite?), and different checks need to be made for different problem categories (for example, in constrained minimization it is necessary to verify whether a trial solution is feasible). Nonetheless, the underlying principles of the various criteria are the same; in non-mathematical terms, they are:

(i) | is the sequence $\left\{{x}^{\left(k\right)}\right\}$ converging? |

(ii) | is the sequence $\left\{{f}^{\left(k\right)}\right\}$ converging? |

(iii) | are the necessary and sufficient conditions for the solution satisfied? |

The decision as to whether a sequence is converging is necessarily speculative. The criterion used in the present routines is to assume convergence if the relative change occurring between two successive iterations is less than some prescribed quantity. Criterion (iii) is the most reliable but often the conditions cannot be checked fully because not all the required information may be available.

2.7.2

Checking results

Little a priori guidance can be given as to the quality of the solution found by a nonlinear optimization algorithm, since no guarantees can be given that the methods will not fail. Therefore, you should always check the computed solution even if the routine reports success. Frequently a ‘solution’ may have been found even when the routine does not report a success. The reason for this apparent contradiction is that the routine needs to assess the accuracy of the solution. This assessment is not an exact process and consequently may be unduly pessimistic. Any ‘solution’ is in general only an approximation to the exact solution, and it is possible that the accuracy you have specified is too stringent.

Further confirmation can be sought by trying to check whether or not convergence tests are almost satisfied, or whether or not some of the sufficient conditions are nearly satisfied. When it is thought that a routine has returned a
nonzero value of **ifail**
only because the requirements for ‘success’ were too stringent it may be worth restarting with increased convergence tolerances.

For constrained problems, check whether the solution returned is feasible, or nearly feasible; if not, the solution returned is not an adequate solution.

Confidence in a solution may be increased by restarting the solver with a different initial approximation to the solution. See Section 8.3 of Gill et al. (1981) for further information.

2.7.3

Monitoring progress

Many of the routines in the chapter have facilities to allow you to monitor the progress of the minimization process, and you are encouraged to make use of these facilities. Monitoring information can be a great aid in assessing whether or not a satisfactory solution has been obtained, and in indicating difficulties in the minimization problem or in the ability of the routine to cope with the problem.

The behaviour of the function, the estimated solution and first derivatives can help in deciding whether a solution is acceptable and what to do in the event of a return with a
nonzero value of **ifail**.

2.7.4

Confidence intervals for least squares solutions

When estimates of the parameters in a nonlinear least squares problem have been found, it may be necessary to estimate the variances of the parameters and the fitted function. These can be calculated from the Hessian of the objective $f\left(x\right)$ at the solution.

In many least squares problems, the Hessian is adequately approximated at the solution by $G=2{J}^{\mathrm{T}}J$ (see Section 2.5.3). The Jacobian, $J$, or a factorization of $J$ is returned by all the comprehensive least squares routines and, in addition, a routine is available in the Library to estimate variances of the parameters following the use of most of the nonlinear least squares routines, in the case that $G=2{J}^{\mathrm{T}}J$ is an adequate approximation.

Let $H$ be the inverse of $G$, and $S$ be the sum of squares, both calculated at the solution $\stackrel{-}{x}$; an unbiased estimate of the **variance** of the $i$th parameter ${x}_{i}$ is

and an unbiased estimate of the covariance of ${\stackrel{-}{x}}_{i}$ and ${\stackrel{-}{x}}_{j}$ is

If ${x}^{*}$ is the true solution, then the $100\left(1-\beta \right)\%\text{}$ **confidence interval** on $\stackrel{-}{x}$ is

where ${t}_{\left(1-\beta /2,m-n\right)}$ is the $100\left(1-\beta \right)/2$ percentage point of the $t$-distribution with $m-n$ degrees of freedom.

$$\mathrm{var}{\stackrel{-}{x}}_{i}=\frac{2S}{m-n}{H}_{ii}$$ |

$$\mathrm{covar}\left({\stackrel{-}{x}}_{i},{\stackrel{-}{x}}_{j}\right)=\frac{2S}{m-n}{H}_{ij}\text{.}$$ |

$${\stackrel{-}{x}}_{i}-\sqrt{\mathrm{var}{\stackrel{-}{x}}_{i}}.{t}_{\left(1-\beta /2,m-n\right)}<{x}_{i}^{*}<{\stackrel{-}{x}}_{i}+\sqrt{\mathrm{var}{\stackrel{-}{x}}_{i}}.{t}_{\left(1-\beta /2,m-n\right)}\text{, \hspace{1em}}i=1,2,\dots ,n$$ |

In the majority of problems, the residuals ${r}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$, contain the difference between the values of a model function $\varphi \left(z,x\right)$ calculated for $m$ different values of the independent variable $z$, and the corresponding observed values at these points. The minimization process determines the parameters, or constants $x$, of the fitted function $\varphi \left(z,x\right)$. For any value, $\stackrel{-}{z}$, of the independent variable $z$, an unbiased estimate of the **variance** of $\varphi $ is

$$\mathrm{var}\varphi =\frac{2S}{m-n}\sum _{i=1}^{n}\sum _{j=1}^{n}{\left[\frac{\partial \varphi}{\partial {x}_{i}}\right]}_{\stackrel{-}{z}}{\left[\frac{\partial \varphi}{\partial {x}_{j}}\right]}_{\stackrel{-}{z}}{H}_{ij}\text{.}$$ |

The $100\left(1-\beta \right)\%$ **confidence interval** on $f$ at the point $\stackrel{-}{z}$ is

For further details on the analysis of least squares solutions see Bard (1974) and Wolberg (1967).

$$\varphi \left(\stackrel{-}{z},\stackrel{-}{x}\right)-\sqrt{\mathrm{var}\varphi}.{t}_{\left(\beta /2,m-n\right)}<\varphi \left(\stackrel{-}{z},{x}^{*}\right)<\varphi \left(\stackrel{-}{z},\stackrel{-}{x}\right)+\sqrt{\mathrm{var}\varphi}.{t}_{\left(\beta /2,m-n\right)}\text{.}$$ |

3

Recommendations on Choice and Use of Available Routines

The choice of routine depends on several factors: the type of problem (unconstrained, etc.); the level of derivative information available (function values only, etc.); your experience (there are easy-to-use versions of some routines); whether or not a problem is sparse; whether or not the routine is to be used in a multithreaded environment; and whether computational time has a high priority. Not all choices are catered for in the current version of the Library.

3.1

Easy-to-use and Comprehensive Routines

Many routines appear in the Library in two forms: a comprehensive form and an easy-to-use form. The purpose of the easy-to-use forms is to make the routine simple to use by including in the calling sequence only those arguments absolutely essential to the definition of the problem, as opposed to arguments relevant to the solution method. If you are an experienced user the comprehensive routines have additional arguments which enable you to improve their efficiency by ‘tuning’ the method to a particular problem. If you are a casual or inexperienced user, this feature is of little value and may in some cases cause a failure because of a poor choice of some arguments.

In the easy-to-use routines, these extra arguments are determined either by fixing them at a known safe and reasonably efficient value, or by an auxiliary routine which generates a ‘good’ value automatically.

For routines introduced since Mark 12 of the Library a different approach has been adopted between the choice of easy-to-use and comprehensive routines. The optimization routine has an easy-to-use argument list, but additional arguments may be changed from their default values by calling an ‘option’ setting routine before the call to the main optimization routine. This approach has the advantages of allowing the options to be given in the form of keywords and requiring only those options that are to be different from their default values to be set.

3.2

Thread Safe Routines

Many of the routines in this chapter come in pairs, with each routine in the pair having exactly the same functionality, except that one of them has additional arguments in order to make it safe for use in multithreaded applications. The routine that is safe for use in multithreaded applications has an ‘a’ as the last character in the name, in place of the usual ‘f’.

3.3

Reverse Communication Routines

Most of the routines in this chapter are called just once in order to compute the minimum of a given objective function subject to a set of constraints on the variables. The objective function and nonlinear constraints (if any) are specified by you and written as subroutines to a very rigid format described in the relevant routine document.

This chapter also contains
a pair of **reverse communication** routines,
e04uff/e04ufa, which solve dense NLP problems using a sequential quadratic programming method.
These
may be convenient to use when the minimization routine is being called from a computer language which does not fully support procedure arguments in a way that is compatible with the Library.
These routines are
also useful if a large amount of data needs to be transmitted into the routine. See Section 3.3.3 in How to Use the NAG Library and its Documentation for more information about reverse communication routines.

3.4

Choosing Between Variant Routines for Some Problems

As evidenced by the wide variety of routines available in Chapter E04, it is clear that no single algorithm can solve all optimization problems. It is important to try to match the problem to the most suitable routine, and that is what the decision trees in Section 4 help to do.

Sometimes in Chapter E04 more than one routine is available to solve precisely the same optimization problem. If their differences lay in the underlying method, refer to the sections above. Section 2.5.4 discusses key features of interior point methods (represented by e04stf) and active-set SQP methods (for example, e04ugf/e04uga or e04vhf). Alternatively, there are routines implementing slightly different variants of the same method (such as e04ucf/e04uca and e04wdf). Experience shows that in this case although both routines can usually solve the same problem and get similar results, sometimes one routine will be faster, sometimes one might find a different local minimum to the other, or, in difficult cases, one routine may obtain a solution when the other one fails.

After using one of these routines, if the results obtained are unacceptable for some reason, it may be worthwhile trying the other routine instead. In the absence of any other information, in the first instance you are recommended to try using e04ucf/e04uca, and if that proves unsatisfactory, try using e04wdf. Although the algorithms used are very similar, the two routines each have slightly different optional parameters which may allow the course of the computation to be altered in different ways.

Other pairs of routines which solve the same kind of problem are e04nqf (recommended first choice) or e04nkf/e04nka, for sparse quadratic or linear programming problems, and e04vhf (recommended) or e04ugf/e04uga, for sparse nonlinear programming. In these cases the argument lists are not so similar as e04ucf/e04uca or e04wdf, but the same considerations apply.

3.5

NAG Optimization Modelling Suite

Mark 26 of the Library introduced the NAG optimization modelling suite, a suite of routines which allows you to define and solve various optimization problems in a uniform manner. The first key feature of the suite is that the definition of the optimization problem and the call to the solver have been separated so it is possible to set up a problem in the same way for different solvers. The second feature is that the problem representation is built up from basic components (for example, a QP problem is composed of a quadratic objective, simple bounds and linear constraints), therefore different types of problems reuse the same routines for their common parts.

A connecting element to all routines in the suite is a handle, a pointer to an internal data structure, which is passed among the routines. It holds all information about the problem, the solution and the solver. Each handle should go through four stages in its life: initialization, problem formulation, problem solution and deallocation.

The initialization is performed by e04raf which creates an empty problem with $n$ decision variables. A call to e04rzf marks the end of the life of the handle as it deallocates all the allocated memory and data within the handle and destroys the handle itself. After the initialization, the objective may be defined as one of the following:
The routines for constraint definition are

These routines may be called in an arbitrary order, however, a call to e04rnf must precede a call to e04rpf for the matrix inequalities with bilinear terms and the nonlinear objective or constraints (e04rgf or e04rkf) must precede the definition of the second derivatives by e04rlf. For further details please refer to the documentation of the individual routines.

The suite also includes the following service routines:

When the problem is fully formulated, the handle can be passed to a solver which is compatible with the defined problem. At the current mark of the Library the NAG optimization modelling suite comprises of e04fff, e04mtf, e04stf and e04svf. The solver indicates by an error flag if it cannot deal with the given formulation. A diagram of the life cycle of the handle is depicted in Figure 2.

3.6

Service Routines

One of the most common errors in the use of optimization routines is that user-supplied subroutines do not evaluate the relevant partial derivatives correctly. Because exact gradient information normally enhances efficiency in all areas of optimization, you are encouraged to provide analytical derivatives whenever possible. However, mistakes in the computation of derivatives can result in serious and obscure run-time errors. Consequently, **service routines** are provided to perform an elementary check on the gradients you supplied. These routines are inexpensive to use in terms of the number of calls they require to user-supplied subroutines.

The appropriate checking routines are as follows:

It should be noted that routines
e04stf, e04ucf/e04uca, e04uff/e04ufa, e04ugf/e04uga, e04usf/e04usa, e04vhf and e04wdf
each incorporate a check on the derivatives being supplied. This involves verifying the gradients at the first point that satisfies the linear constraints and bounds. There is also an option to perform a more reliable (but more expensive) check on the individual gradient elements being supplied. Note that the checks are not infallible.

A second type of service routine computes a set of finite differences to be used when approximating first derivatives. Such differences are required as input arguments by some routines that use only function evaluations.

e04ycf estimates selected elements of the variance-covariance matrix for the computed regression parameters following the use of a nonlinear least squares routine.

e04xaf/e04xaa estimates the gradient and Hessian of a function at a point, given a routine to calculate function values only, or estimates the Hessian of a function at a point, given a routine to calculate function and gradient values.

3.7

Function Evaluations at Infeasible Points

All the solvers for constrained problems based on an active-set method will ensure that any evaluations of the objective function occur at points which **approximately** (up to the given tolerance) satisfy any **simple bounds** or **linear constraints**.

There is no attempt to ensure that the current iteration satisfies any nonlinear constraints. If you wish to prevent your objective function being evaluated outside some known region (where it may be undefined or not practically computable), you may try to confine the iteration within this region by imposing suitable simple bounds or linear constraints (but beware as this may create new local minima where these constraints are active).

Note also that some routines allow you to return the argument
(**iflag**, **inform**, **mode** or **status**)
with a negative value to indicate when the objective function (or nonlinear constraints where appropriate) cannot be evaluated. In case the routine cannot recover (e.g., cannot find a different trial point), it forces an immediate clean exit from the routine.

3.8

Related Problems

Apart from the standard types of optimization problem, there are other related problems which can be solved by routines in this or other chapters of the Library.

h02bbf solves **dense integer LP** problems, h02cbf solves **dense integer QP** problems, h02cef solves **sparse integer QP** problems, h02daf solves **dense mixed integer NLP** problems and h03abf solves a special type of such problem known as a **‘transportation’ problem**.

Several routines in Chapters F04 and F08 solve **linear least squares problems**, i.e., $\mathrm{minimize}{\displaystyle \sum _{i=1}^{m}}{r}_{i}{\left(x\right)}^{2}$ where ${r}_{i}\left(x\right)={b}_{i}-{\displaystyle \sum _{j=1}^{n}}{a}_{ij}{x}_{j}$.

e02gaf solves an overdetermined system of linear equations in the ${l}_{1}$ norm, i.e., minimizes $\sum _{i=1}^{m}}\left|{r}_{i}\left(x\right)\right|$, with ${r}_{i}$ as above, and e02gbf solves the same problem subject to linear inequality constraints.

e02gcf solves an overdetermined system of linear equations in the ${l}_{\infty}$ norm, i.e., minimizes $\underset{i}{\mathrm{max}}}\phantom{\rule{0.25em}{0ex}}\left|{r}_{i}\left(x\right)\right|$, with ${r}_{i}$ as above.

Chapter E05 contains routines for global minimization.

Section 2.5.5 describes how a multi-objective optimization problem might be addressed using routines from this chapter and from Chapter E05.

4

Decision Trees

no objective | linear | quadratic | nonlinear | sum of squares | |

unconstrained | QP See Tree 2 |
NLP See Tree 3 |
LSQ See Tree 4 |
||

simple bounds | LP See Tree 1 |
LP See Tree 1 |
QP See Tree 2 |
NLP See Tree 3 |
LSQ See Tree 4 |

linear | LP See Tree 1 |
LP See Tree 1 |
QP See Tree 2 |
NLP See Tree 3 |
LSQ See Tree 4 |

nonlinear | NLP See Tree 3 |
NLP See Tree 3 |
NLP See Tree 3 |
NLP See Tree 3 |
LSQ See Tree 4 |

matrix inequalities | e04svf | e04svf | e04svf |

Is the problem sparse/large-scale? | e04mtf, e04nqf, e04nkf | |||

yes | ||||

no | ||||

e04mff, e04ncf | ||||

Is the problem sparse/large-scale? | Is it convex? | e04nqf, e04stf, e04nkf | |||||

yes | yes | ||||||

no | no | ||||||

e04stf, e04vhf, e04ugf | |||||||

Is it convex? | e04ncf | ||||||

yes | |||||||

no | |||||||

e04nff | |||||||

Is the problem sparse/large-scale? | Is it unconstrained? | Are first derivatives available? | e04stf, e04dga, e04vhf, e04uga | |||||||

yes | yes | yes | ||||||||

no | no | no | ||||||||

e04vhf, e04uga | ||||||||||

Are first derivatives available? | Are second derivatives available? | e04stf | ||||||||

yes | yes | |||||||||

no | no | |||||||||

e04vhf, e04stf, e04uga | ||||||||||

e04vhf, e04uga | ||||||||||

Are there linear or nonlinear constraints? | e04uca, e04ufa, e04wdf | |||||||||

yes | ||||||||||

no | ||||||||||

Is there only one variable? | Are first derivatives available? | e04bba | ||||||||

yes | yes | |||||||||

no | no | |||||||||

e04aba | ||||||||||

Is it unconstrained with the objective with many discontinuities? | e04cbf or e05saf | |||||||||

yes | ||||||||||

no | ||||||||||

Are first derivatives available? | Are second derivatives available? | Are you an experienced user? | e04lbf | |||||||

yes | yes | yes | ||||||||

no | no | no | ||||||||

e04lyf | ||||||||||

Are many function evaluations problematic? | Are you an experienced user? | e04uca, e04ufa, e04wdf | ||||||||

yes | yes | |||||||||

no | no | |||||||||

e04kyf | ||||||||||

Are you an experienced user? | e04kdf | |||||||||

yes | ||||||||||

no | ||||||||||

e04kzf | ||||||||||

Is the objective expensive to evaluate or noisy? | e04jcf | |||||||||

yes | ||||||||||

no | ||||||||||

Are you an experienced user? | e04uca, e04ufa, e04wdf | |||||||||

yes | ||||||||||

no | ||||||||||

e04jyf | ||||||||||

Is the objective sum of squared linear functions and no nonlinear constraints? | Are there linear constraints? | e04ncf | ||||||||

yes | yes | |||||||||

no | no | |||||||||

Are there simple bounds? | e04pcf, e04ncf | |||||||||

yes | ||||||||||

no | ||||||||||

Chapters F04, F07 or F08 or e04pcf, e04ncf | ||||||||||

Are there linear or nonlinear constraints? | e04usf | |||||||||

yes | ||||||||||

no | ||||||||||

Are there bound constraints? | Are first derivatives available? | e04usf | ||||||||

yes | yes | |||||||||

no | no | |||||||||

e04fff | ||||||||||

Are first derivatives available? | Are second derivatives available? | Are you an experienced user? | e04hef | |||||||

yes | yes | yes | ||||||||

no | no | no | ||||||||

e04hyf | ||||||||||

Are many function evaluations problematic? | Are you an experienced user? | e04gbf | ||||||||

yes | yes | |||||||||

no | no | |||||||||

e04gyf | ||||||||||

Are you an experienced user? | e04gdf | |||||||||

yes | ||||||||||

no | ||||||||||

e04gzf | ||||||||||

e04fff, e04fcf | ||||||||||

5

Functionality Index

6

Auxiliary Routines Associated with Library Routine Arguments

e04cbk | nagf_opt_uncon_simplex_dummy_monit See the description of the argument monit in e04cbf. |

e04fcv | nagf_opt_lsq_uncon_quasi_deriv_comp_lsqlin_fun See the description of the argument lsqlin in e04gbf. |

e04fdz | nagf_opt_lsq_dummy_lsqmon See the description of the argument lsqmon in e04fcf, e04gdf and e04hef. |

e04ffu | nagf_opt_bobyqa_ls_dummy_monit See the description of the argument mon in e04fff. |

e04hev | nagf_opt_lsq_uncon_quasi_deriv_comp_lsqlin_deriv See the description of the argument lsqlin in e04gbf. |

e04jcp | nagf_opt_bounds_bobyqa_func_dummy_monfun See the description of the argument monfun in e04jcf. |

e04mtu | nagf_opt_lp_imp_dummy_monit See the description of the argument monit in e04mtf. |

e54nfu | nagf_opt_qp_dense_sample_qphess See the description of the argument qphess in e04nff/e04nfa and h02cbf. |

e04nfu | nagf_opt_qp_dense_sample_qphess_old See the description of the argument qphess in e04nff/e04nfa and h02cbf. |

e54nku | nagf_opt_qpconvex1_sparse_dummy_qphx See the description of the argument qphx in e04nkf/e04nka and h02cef. |

e04nku | nagf_opt_qpconvex1_sparse_dummy_qphx_old See the description of the argument qphx in e04nkf/e04nka and h02cef. |

e04nsh | nagf_opt_qpconvex2_sparse_dummy_qphx See the description of the argument qphx in e04nqf. |

e04stu | nagf_opt_ipopt_dummy_mon See the description of the argument mon in e04stf. |

e04stv | nagf_opt_ipopt_dummy_objfun See the description of the argument objfun in e04stf. |

e04stw | nagf_opt_ipopt_dummy_objgrd See the description of the argument objgrd in e04stf. |

e04stx | nagf_opt_ipopt_dummy_confun See the description of the argument confun in e04stf. |

e04sty | nagf_opt_ipopt_dummy_congrd See the description of the argument congrd in e04stf. |

e04stz | nagf_opt_ipopt_dummy_hess See the description of the argument hess in e04stf. |

e04udm | nagf_opt_nlp1_dummy_confun See the description of the argument confun in e04ucf/e04uca and e04usf/e04usa. |

e04ugm | nagf_opt_nlp1_sparse_dummy_confun See the description of the argument confun in e04ugf/e04uga. |

e04ugn | nagf_opt_nlp1_sparse_dummy_objfun See the description of the argument objfun in e04ugf/e04uga. |

e04wdp | nagf_opt_nlp2_dummy_confun See the description of the argument confun in e04wdf. |

7

Routines Withdrawn or Scheduled for Withdrawal

The following lists all those routines that have been withdrawn since Mark 19 of the Library or are scheduled for withdrawal at one of the next two marks.

WithdrawnRoutine | Mark ofWithdrawal | Replacement Routine(s) |

e04ccf/e04cca | 24 | e04cbf |

e04fdf | 19 | e04fyf |

e04gcf | 19 | e04gyf |

e04gef | 19 | e04gzf |

e04hff | 19 | e04hyf |

e04jaf | 19 | e04jyf |

e04kaf | 19 | e04kyf |

e04kcf | 19 | e04kzf |

e04laf | 19 | e04lyf |

e04unf | 22 | e04usf/e04usa |

e04upf | 19 | e04usf/e04usa |

e04zcf/e04zca | 24 | No longer required |

8

References

Bard Y (1974) *Nonlinear Parameter Estimation* Academic Press

Chvátal V (1983) *Linear Programming* W.H. Freeman

Dantzig G B (1963) *Linear Programming and Extensions* Princeton University Press

Fletcher R (1987) *Practical Methods of Optimization* (2nd Edition) Wiley

Gill P E and Murray W (ed.) (1974) *Numerical Methods for Constrained Optimization* Academic Press

Gill P E, Murray W and Wright M H (1981) *Practical Optimization* Academic Press

Murray W (ed.) (1972) *Numerical Methods for Unconstrained Optimization* Academic Press

Nocedal J and Wright S J (2006) *Numerical Optimization* (2nd Edition) Springer Series in Operations Research, Springer, New York

Wolberg J R (1967) *Prediction Analysis* Van Nostrand