Gradient and Hessian approximations in Derivative Free Optimization

This work considers unconstrained optimization problems of the form

$\displaystyle\min_{x\in\mathbf{R}^{n}}f(x),$

(1)

where $f:\mathbf{R}^{n}\to\mathbf{R}$ is smooth, but the derivatives of $f$ are unavailable. This is a typical derivative free optimization setting. Problems of this form arise, for example when the function is smooth but it is computationally impractical to obtain derivative information, or when the function itself is not known explicitly although function values are available from a black-box or oracle, (for example via physical measurements).

A key ingredient of many algorithms for solving (1) in a derivative free setting is an estimate of the gradient of $f$. Sometimes gradients are used to determine search directions (for example, in implicit filtering algorithms [5, 18, 19]), and sometimes they are used in the definition of an algorithm’s stopping criteria. Several techniques exist for determining estimates to derivatives, including automatic differentiation [21, 3, 25], finite difference approximations, interpolation models [8, 11], Gaussian smoothing [24, 29], and smoothing on the unit sphere [17]. This work investigates finite differences and interpolation models, so henceforth, we focus on these approaches.

A well known, simple, and widely applicable approach to estimating derivatives is finite difference approximations (or finite differences). Finite differences are so called because they provide estimates to derivatives by taking the difference between function values at certain known points. For example, for a function $f:\mathbf{R}^{n}\to\mathbf{R}$, the forward differences approximation $g(x)$ to the (first) derivative is

$$[g(x)]_{i}=\frac{f(x+he_{i})-f(x)}{h},\qquad\text{for}\;i=1,\dots,n,$$

(2)

which is an $\mathcal{O}(h)$ accurate approximation. Many other (well known) formulations exist, including approximations that provide a higher accuracy estimates of the gradient, as well as formulations for higher order derivatives. The standard finite difference approximation formulations can be used whenever function values evaluated along coordinate directions are available. The formulations involve minimal computations ($\mathcal{O}(n)$ operations), and they have low memory requirements because one does not need to store the coordinate directions, but finite differences can suffer from cancellation errors.

Another approach to approximating derivatives is to use interpolation models [9, 33, 20]. Interpolation models can be viewed as a generalization of finite differences, where, instead of being restricted to the coordinate directions, function evaluations along a general set of directions are utilized, and derivative approximations are found by solving an associated system of equations (i.e., minimizing the interpolation model). Linear models lead to approximations to the gradient [4], while quadratic models provide approximations to the gradient and Hessian [8]. Much research has focused on the ‘quality’ of the set of directions, because this impacts the accuracy of the derivative estimates [2, 6, 7, 16].

The purpose of this work is to show that, if one uses a (diagonal) quadratic interpolation model with carefully chosen interpolation directions and $2n+2$ function values, then approximations to the gradient and diagonal of the Hessian can be obtained in $\mathcal{O}(n)$ computations and with only $\mathcal{O}(1)$ vectors of storage. The solution to this specific instance of an interpolation model can be thought of, in some sense, as a general finite difference type approximation to the gradient and diagonal of the Hessian, because the solution to the interpolation model is a formula that involves sums/differences of function values at the interpolation points. Although the focus of this work is on derivative estimates and not algorithmic development, the estimates here are widely applicable, and could be employed by many algorithms that require gradient approximations. Moreover, the diagonal Hessian approximation is readily available, and could be used as preconditioner within an algorithm. (See the numerical experiments in Section 6, which include a derivative free preconditioned conjugate gradients algorithm.)

It remains to note that the approaches mentioned above all require access to function values, and the cost of evaluating a function in a derivative free optimization setting can vary greatly depending upon the particular application being considered. It is up to the user to determine their function evaluation budget, and the approaches in this work require between $n+1$ and $2n+2$ function evaluations to generate derivative approximations.

This section describes how to approximate derivatives of a function via interpolation models. The Taylor series for a function $f$ about a point $x\in\mathbf{R}^{n}$ is

$$f(y)=f(x)+(y-x)^{T}\nabla f(x)+\tfrac{1}{2}(y-x)^{T}\nabla^{2}f(x)(y-x)+\mathcal{O}(\|y-x\|_{2}^{3}).$$

(7)

Let $g(x)$ denote an approximation to the gradient $\nabla f(x)$ (i.e., $g(x)\approx\nabla f(x)$), and let $H(x)$ denote an approximation to the Hessian $\nabla^{2}f(x)$ (i.e., $H(x)\approx\nabla^{2}f(x)$). Then (7) becomes

$$f(y)\approx f(x)+(y-x)^{T}g(x)+\tfrac{1}{2}(y-x)^{T}H(x)(y-x).$$

(8)

In this work, only a diagonal approximation to the Hessian is considered, so define

$$D(x):=\begin{bmatrix}d_{1}&&\\ &\ddots&\\ &&d_{n}\end{bmatrix}\in\mathbf{R}^{n\times n},\qquad d(x):=\begin{bmatrix}d_{1}\\ \vdots\\ d_{n}\end{bmatrix}\in\mathbf{R}^{n},$$

(9)

where $D(x)\approx{\rm diag}(\nabla^{2}f(x))$. Combining (8) and (9) leads to the (diagonal) quadratic model

$$m(y)=f(x)+(y-x)^{T}g(x)+\tfrac{1}{2}(y-x)^{T}D(x)(y-x).$$

(10)

The motivation for this diagonal quadratic model is as follows. It is well known that, if one uses a linear model with $n$ (or $n+1$) interpolation points (i.e., $n$ or $n+1$ function evaluations), then an estimate of the gradient is available that is $\mathcal{O}(h)$ accurate (see also Section 5 which summarises results for linear models). For some situations, an $\mathcal{O}(h)$ accurate gradient is not of sufficient accuracy to make algorithmic progress, so an $\mathcal{O}(h^{2})$ accurate gradient is used instead. In this case, one can still use a linear model, but $2n$ (or $2n+2$) interpolation points (i.e., $2n$ (or $2n+2$) function evaluations) are needed. In derivative free optimization, function evaluations are often expensive to obtain, so it is prudent to ‘squeeze out’ as much information from them as possible. Hence, suppose one had available $2n$ (or $2n+2$) function evaluations because an accurate gradient was desired. The model (10) contains $2n$ unknowns, so $2n$ interpolation points and function values is enough to uniquely determine an $\mathcal{O}(h^{2})$ accurate gradient $g(x)$, as well as an approximation to the pure second derivatives on the diagonal of $D(x)$, i.e., no additional function evaluations are required to compute the pure second derivatives. Although $D(x)$ can be determined ‘for free’ in terms of function evaluations (assuming $2n$ are available for the gradient estimate), we will also show that the diagonal matrix $D(x)$ can be computed in $\mathcal{O}(n)$ flops, so the linear algebra costs are low. This is similar to the case for finite differences (i.e., one requires $2n(+1)$ function evaluations to compute $g(x)$ via the central differences formula, and those same function values can also be used to compute the pure second derivatives), and so this diagonal quadratic model (10) is useful because one can obtain analogous results for a certain set of interpolation points.

The goal is to determine $g(x)$ and $D(x)$ from the interpolation model (10) using a set of known points and the corresponding function values. Hence, one must build a sample set (a set of interpolation points). In general, one can use any number of interpolation points, although this affects the model. For a linear model (omitting the third term in (10)), one can determine a unique solution $g(x)$ if there are $n$ interpolation/sample points with known function values at those points (see Section 5). For a general quadratic model with Hessian approximation $H(x)$, one requires $(n+1)(n+2)/2$ function values to determine unique solutions $g(x)$ and $H(x)$. If one has access to more or fewer sample points/function values, then one can use, for example, a least squares approach to find approximations to $g(x)$ and $H(x)$. In the rest of this section, $2n$ interpolation points will be used, and it will become clear that this allows unique solutions $g(x)$ and $D(x)$ to the model (10).

Given $x$, a set of directions $\{u_{1},\dots,u_{n}\}\subset\mathbf{R}^{n}$, and a scalar $h$ (sometimes referred to as the sampling radius), a set of points $\{x_{1},\dots,x_{n}\}$ is defined via the equations

$\displaystyle x_{j}=x+hu_{j},\quad\|u_{j}\|_{2}=\tfrac{1}{|h|}\|x_{j}-x\|_{2},\qquad j=1,\dots,n.$

(11)

The model is constructed to satisfy the interpolation conditions

$$f(x+hu_{j})=m(x+hu_{j}),\qquad j=1,\dots,n.$$

(12)

Substituting $y=x+hu_{j}=x_{j}$ $\forall j$ into the model (10) gives the equations

$$m(x_{j})=f(x)+hu_{j}^{T}g(x)+\tfrac{1}{2}h^{2}u_{j}^{T}D(x)u_{j}.$$

(13)

Because $D(x)$ is diagonal, $u_{j}^{T}D(x)u_{j}=\sum_{i=1}^{n}d_{i}(u_{j})_{i}^{2}=(u_{j}\odot u_{j})^{T}d(x),$ $\forall j$. Defining

	$\displaystyle U$	$\displaystyle:=$	$\displaystyle\begin{bmatrix}u_{1}&\dots&u_{n}\end{bmatrix}$		(14)
	$\displaystyle W$	$\displaystyle:=$	$\displaystyle\begin{bmatrix}u_{1}\odot u_{1}&u_{2}\odot u_{2}&\dots&u_{n}\odot u_{n}\end{bmatrix},$		(15)

where $U,W\in\mathbf{R}^{n\times n}$, allows (13) to be written in matrix form as

$\displaystyle\delta\mathbf{f}\overset{\eqref{eq:df}}{=}hU^{T}g(x)+\tfrac{1}{2}h^{2}W^{T}d(x).$

(16)

The system (16) is underdetermined, (there are only $n$ equations in $2n$ unknowns), and as the goal is to find unique solutions $g(x)$ and $D(x)={\rm diag}(d(x))$, an additional $n$ equations are required. Define

$\displaystyle h^{\prime}:=\eta h,\qquad\text{where }\eta\neq 1.$

(17)

Now, keep the point $x$, and the set of directions $\{u_{1},\dots,u_{n}\}\subset\mathbf{R}^{n}$ fixed and choose $h^{\prime}$ satisfying (17). A new set of points $\{x_{1}^{\prime},\dots,x_{n}^{\prime}\}$ is defined via the equations

$\displaystyle x_{j}^{\prime}=x+h^{\prime}u_{j},\quad\|u_{j}\|_{2}=\tfrac{1}{|h^{\prime}|}\|x_{j}^{\prime}-x\|_{2},\qquad j=1,\dots,n.$

(18)

The definition of $h^{\prime}$ in (17) ensures that $x_{j}^{\prime}\neq x_{j}$ $\forall j$. The model is constructed to satisfy the interpolation conditions

$$f(x+\eta hu_{j})=m(x+\eta hu_{j})\qquad j=1,\dots,n.$$

(19)

(That is, $f(x_{j})=m(x_{j})$ and $f(x_{j}^{\prime})=m(x_{j}^{\prime})$ for $j=1,\dots,n$.) Following similar arguments to those previously established, one arrives at the system

	$\displaystyle\delta\mathbf{f}^{\prime}$	$\displaystyle\overset{\eqref{eq:df}}{=}$	$\displaystyle h^{\prime}U^{T}g(x)+\tfrac{1}{2}(h^{\prime})^{2}W^{T}d(x)$		(20)
		$\displaystyle=$	$\displaystyle\eta hU^{T}g(x)+\tfrac{1}{2}\eta^{2}h^{2}W^{T}d(x).$

The systems (16) and (20) can be combined into block matrix form as follows:

$\displaystyle\begin{bmatrix}\delta\mathbf{f}\\ \delta\mathbf{f}^{\prime}\end{bmatrix}=\begin{bmatrix}hU^{T}&\tfrac{1}{2}h^{2}W^{T}\\ \eta hU^{T}&\tfrac{1}{2}\eta^{2}h^{2}W^{T}\end{bmatrix}\begin{bmatrix}g(x)\\ d(x)\end{bmatrix}.$

(21)

It is clear that (22) represents a system of $2n$ equations in 2n unknowns. Performing block row operations on the system (22) (i.e., multiplying the first row by $\eta$ and subtracting the second row, and also multiplying the first row by $\eta^{2}$ and subtracting the second row) gives

$\displaystyle\begin{bmatrix}\eta\delta\mathbf{f}-\delta\mathbf{f}^{\prime}\\ \eta^{2}\delta\mathbf{f}-\delta\mathbf{f}^{\prime}\end{bmatrix}=\begin{bmatrix}0&\tfrac{1}{2}h^{2}(\eta-\eta^{2})W^{T}\\ (\eta^{2}-\eta)hU^{T}&0\end{bmatrix}\begin{bmatrix}g(x)\\ d(x)\end{bmatrix}.$

(22)

Hence, estimates of the gradient and diagonal of the Hessian can be found by solving

$\displaystyle y=hU^{T}g(x),\qquad\text{and}\qquad z=\tfrac{1}{2}h^{2}W^{T}d(x),$

(23)

where

	$\displaystyle y$	$\displaystyle:=$	$\displaystyle\tfrac{1}{\eta(\eta-1)}(\eta^{2}\delta\mathbf{f}-\delta\mathbf{f}^{\prime})$		(24)
	$\displaystyle z$	$\displaystyle:=$	$\displaystyle\tfrac{1}{\eta(1-\eta)}(\eta\delta\mathbf{f}-\delta\mathbf{f}^{\prime}).$		(25)

Questions naturally arise, such as ‘when do unique solutions to (23) exist?’ and ‘what is the cost of computing the solutions?’. Clearly, if the matrices $U$ and $W$ have full rank, then unique solutions to (23) exist. The rank of $U$ and $W$, of course, depends upon how the interpolation points are chosen. If the interpolation directions $\{u_{1},\dots,u_{n}\}$ form a basis for $\mathbf{R}^{n}$, then a unique solution $g(x)$ exists. Similarly, if the vectors $w_{j}=u_{j}\odot u_{j}$, for $j=1,\dots,n$ form a basis for $\mathbf{R}^{n}$, then a unique solution $d(x)$ exists.

Computing the solution to a general unstructured system has a cost of $\mathcal{O}(n^{3})$ flops. The focus of the remainder of this paper is to show that, if one selects the interpolation points/interpolation directions in a specific way, then unique solutions $g(x)$ and $d(x)$ to (23) exist, and moreover, the computational cost of finding the solutions remains low, at $\mathcal{O}(n)$. This is the same cost as finite difference approximations to the gradient and pure second order derivatives. Importantly, it will also be shown that this specific set of directions does not increase the memory footprint either, requiring only $\mathcal{O}(1)$ vectors of storage.

In what follows, it will be useful to notice that, in the special case $\eta=-1$, (24) and (25) become

	$\displaystyle y$	$\displaystyle=$	$\displaystyle\tfrac{1}{2}(\delta\mathbf{f}-\delta\mathbf{f}^{\prime}),\qquad[y]_{i}=\tfrac{1}{2}(f(x+hu_{j})-f(x-hu_{j}))$		(26)
	$\displaystyle z$	$\displaystyle=$	$\displaystyle\tfrac{1}{2}(\delta\mathbf{f}+\delta\mathbf{f}^{\prime}),\qquad[z]_{i}=\tfrac{1}{2}(f(x+hu_{j})+f(x-hu_{j})-2f(x)).$		(27)

The Sherman-Morrison formula for a nonsingular matrix $A\in\mathbf{R}^{n\times n}$, and vectors $u,v\in\mathbf{R}^{n}$ is used many times in this work, and is stated now for convenience:

$$(A+uv^{T})^{-1}=A^{-1}-\frac{A^{-1}uv^{T}A^{-1}}{1+v^{T}A^{-1}u}.$$

(28)

This section describes how to compute solutions to the diagonal quadratic interpolation model when the coordinate basis is used to construct the interpolation points, and also when the regular basis is employed. In both cases the computational cost of determining the solution is $\mathcal{O}(n)$, and when the coordinate basis is used with $\eta=-1$ in (17), the usual finite differences approximations to the first and second (pure) derivatives are recovered.

Spendley et al. [32] and Nelder and Mead [28] carried out some of the pioneering work on derivative-free simplex-based algorithms. This work continues today with studies on properties of simplicies [1, 22] as well as work on simplex algorithm development, [31]. These, and many other algorithms in derivative free optimization, employ positive bases and simplices, and involve geometric operations such as expansions, rotations and contractions. In such cases, the original simplex/positive basis involves $n+1$ (or possibly more) points (with their corresponding function values), and if a geometric operation is performed, then this often doubles the number of interpolation points and corresponding function values that are available. As previously mentioned, function values are expensive, so it is prudent to make full use of them.

This section is motivated in part, by the ongoing research on simplex based methods, and the following types of ideas. Consider a simplex based algorithm, where one arrives at an iteration that involves an expansion step. In such a case, one would have access to $2n+2$ function values (measured at the vertices of the 2 resulting original and expanded simplices) as a by-product. If an approximation to the gradient or diagonal of the Hessian was readily available, then this might be useful curvature information that could be built into a simplex based algorithm to help locate the solution more quickly. Moreover, depending on the type of geometric operation involved, the vertices of the 2 simplices are generated from the same set of base directions, one scaled by $h$, and one scaled by some $h^{\prime}=\eta h$. Therefore, in this work we have not restricted to $\eta=-1$, to ensure the results are applicable to expansions, contractions, and not just rotations.

In this section we investigate the use of minimal positive bases (and their associated simplices) to generate interpolation points that can be used in the model (13). We will show that the results from the previous section carry over, in a natural way, when there are $2n+2$ (rather than $2n$) interpolation points.