Optimisation

3.3.1 Objective function transects are a partial view

Clearly plotting some transects through your objective function is a good idea, but they can only give a limited and partial view when the $\bm \theta$ is multidimensional. Consider this $x$-transect through a function, $f(x,y)$

Apparently the function has many local minima, and optimisation will be difficult. Actually the function is this…

It has one local minimum, its global minimum. Head downhill from any point $x,y$ and you will eventually reach the minimum.

The opposite problem can also occur. Here are $x$ and $y$ transects through a second function $g(x,y)$.

…apparently no problem. From these plots you would be tempted to conclude that $g$ is a well behaved and uni-modal. The problem is that $g$ looks like this…

So, generally speaking, it is a good idea to plot transects through your objective function before optimisation, and to plot transects passing through the apparent optimum after optimisation, but bear in mind that transects only give partial views.

3.4.1 Taylor’s Theorem (and conditions for an optimum)

We need only a limited version of Taylor’s theorem. Suppose that $f$ is a twice continuously differentiable function of $\bm \theta$ and that $\bm \Delta$ is of the same dimension as $\bm \theta$. Then

for some $t \in (0,1)$, where

From Taylor’s theorem it is straightforward to establish that the conditions for $f({\bm \theta }^*)$ to be a minimum of $f$ are that $\nabla f({\bm \theta }^*) = {\bf 0}$ and $\nabla ^2 f({\bm \theta }^*)$ is positive semi-definite.

Here we have given an exact form of Taylor’s theorem for a real-valued many-variable function, for a first order approximation plus a second order remainder term in Langrange’s form. In practice, in applications, we typically work with approximate versions of the above result. In some cases we work directly with a linear (first-order) approximation, by setting the remainder term to zero. In other cases we work with the quadratic (second-order) approximation, obtained by evaluating the remainder term at $t=0$.

3.4.2 Steepest descent (AKA gradient descent)

Given some parameter vector ${\bm \theta }^{[k]}$, it might seem reasonable to simply ask which direction would lead to the most rapid decrease in $f$ for a sufficiently small step, and to use this as the search direction. If $\bm \Delta$ is small enough then Taylor’s theorem implies that

For $\bm \Delta$ of fixed length the greatest decrease will be achieved by minimising the inner product $\nabla f({\bm \theta }^{[k]}) ^{\sf T}{\bm \Delta }$, which is achieved by making $\bm \Delta$ parallel to $-\nabla f({\bm \theta }^{[k]})$, i.e. by choosing

To firm things up, let’s set ${\bm \Delta }= -\nabla f({\bm \theta }^{[k]})/\| \nabla f({\bm \theta }^{[k]}) \|$, and consider a new proposed value of the form ${\bm \theta }^{[k]}+\alpha {\bm \Delta }$ for $\alpha >0$. Then, by plugging back into our linear approximation, we have

and so as we increase $\alpha$ from zero, we decrease $f(\cdot )$ at a rate proportional to the magnitude of the gradient vector. Any other choice of ${\bm \Delta }$ will decrease the function value at a slower rate.

Given a search direction, we need a method for choosing how far to move from ${\bm \theta }^{[k]}$ along this direction. Here there is a trade-off. We need to make reasonable progress, but do not want to spend too much time choosing the exactly optimal distance to move, since the point will only be abandoned at the next iteration anyway. Usually it is best to try and choose the step length $\alpha$ so that $f({\bm \theta }^{[k]}+\alpha {\bm \Delta })$ is sufficiently much lower that $f({\bm \theta }^{[k]})$, and also so that the magnitude of the gradient of $f$ in the ${\bm \Delta }$ direction is sufficiently reduced by the step. See e.g. Section 3.1 of Nocedal and Wright (2006) for further details. For the moment presume that we have a step selection method.

The following illustrates the progress of the steepest descent method on Rosenbrock’s function. The contours are $\log _{10}$ of Rosenbrock’s function, and the thick line shows the path of the algorithm.

It took 14,845 steps to reach the minimum. This excruciatingly slow progress is typical of the steepest descent method. Another disadvantage of steepest descent is scale dependence: imagine linearly re-scaling so that the $z$ axis in the above plot ran from -5 to 15, but the contours were unchanged. This would change all the steepest descent directions. By contrast, the methods considered next are scale invariant.

3.4.3 Newton’s method

As steepest descent approaches the minimum, the first derivative term retained in the Taylor approximation of $f$ becomes negligible relative to the neglected second derivative term. This largely explains the method’s poor performance. Retaining the second derivative term in the Taylor expansion yields the method of choice for a wide range of problems: Newton’s method. Taylor’s theorem implies that

Differentiating the r.h.s. w.r.t. $\bm \Delta$ and setting the result to $\bf 0$ we have

as the equation to solve for the step $\bm \Delta$ which will take us to the turning point of the (quadratic) Taylor approximation of $f$. This turning point will be a minimum (of the approximation) if $\nabla ^2 f({\bm \theta }^{[k]})$ is positive definite, and can also be derived directly without vector calculus by completing the square.

Notice how Newton’s method automatically suggests a search direction and a step length. Usually we simply accept the Newton step, $\bm \Delta$, unless it fails to lead to sufficient decrease in $f$. In the latter case the step length is repeatedly halved until sufficient decrease is achieved.

The second practical issue with Newton’s method is the requirement that $\nabla ^2 f({\bm \theta }^{[k]})$ is positive definite. There is no guarantee that it will be, especially if ${\bm \theta }^{[k]}$ is far from the optimum. This would be a major problem, were it not for the fact that the solution to

is a descent direction for any positive definite ${\bf H}$. To see this, put ${\bm \Delta } = -{\bf H}^{-1}\nabla f({\bm \theta }^{[k]})$ as a search direction in our first-order approximation,

for small $\alpha >0$ due to the positive-definiteness of ${\bf H}^{-1}$. This gives us the freedom to modify $\nabla ^2 f({\bm \theta }^{[k]})$ to make it positive definite, and still be guaranteed a descent direction. One approach is to take the symmetric eigen-decomposition of $\nabla ^2 f({\bm \theta }^{[k]})$, reverse the sign of any negative eigenvalues, replace any zero eigenvalues with a small positive number and then replace $\nabla ^2 f({\bm \theta }^{[k]})$ by the matrix with this modified eigen-decomposition. Computationally cheaper alternatives simply add a multiple of the identity matrix to $\nabla ^2 f({\bm \theta }^{[k]})$, sufficiently large to make the result positive definite. Note that steepest descent corresponds to the choice ${\bm H} = {\bm I}$, and so inflating the diagonal of the Hessian corresponds to nudging the search direction from that of the Newton method towards that of steepest descent, and this is usually a fairly safe and conservative thing to do.

Combining these steps gives a practical Newton algorithm. Starting with $k=0$ and a guesstimate ${\bm \theta }^{[0]}$…

1. 1.

   Evaluate $f({{\bm \theta }^{[k]}})$, $\nabla f({\bm \theta }^{[k]})$ and $\nabla ^2 f({\bm \theta }^{[k]})$, implying a quadratic approximation to $f$, via Taylor’s theorem.

2. 2.

   Test whether ${\bm \theta }^{[k]}$ is a minimum, and terminate if it is.

3. 3.

   If $\nabla ^2 f({\bm \theta }^{[k]})$ is not positive definite, perturb it so that it is (thereby modifying the quadratic model of $f$).

4. 4.

   The search direction $\bm \Delta$ is now defined by the solution to

   $$\nabla ^2 f({\bm \theta }^{[k]}) {\bm \Delta } = - \nabla f({\bm \theta }^{[k]}),$$

   i.e. the $\bm \Delta$ minimising the current quadratic model.

5. 5.

   If $f({\bm \theta }^{[k]}+{\bm \Delta })$ is not sufficiently lower than $f({\bm \theta }^{[k]})$, repeatedly halve $\bm \Delta$ until it is.

6. 6.

   Set ${\bm \theta }^{[k+1]} = {\bm \theta }^{[k]}+{\bm \Delta }$, increment $k$ by one and return to step 1.

The following figure illustrates the progress of the method on Rosenbrock’s function, working across rows from top left. Grey contours show the function, dotted show the quadratic approximation at the current best point, and, if present, black contours are the quadratic approximation before modification of the Hessian. The thick black line tracks the steps accepted, while the thin black line shows any suggested steps that were reduced. Note the fact that the method converged in 20 steps, compared to the 14,845 steps taken by steepest descent.

3.4.4 Quasi-Newton methods

Newton’s method usually performs very well, especially near the optimum, but to use it we require the second derivative (Hessian) matrix of the objective. This can be tedious to obtain. Furthermore, if the dimension of $\bm \theta$ is large, then solving for the Newton direction can be a substantial part of the computational burden of optimisation. Is it possible to produce methods that only require first derivative information, but with efficiency rivalling Newton’s method, rather than steepest descent?

The answer is yes. One approach is to replace the exact Hessian, of Newton’s method, with an approximate Hessian based on finite differences of gradients. The finite difference intervals must be chosen with care (see section 4.2), but the approach is reliable (recall that any positive definite matrix in place of the Hessian will yield a descent direction, so an approximate Hessian should be unproblematic). However, finite differencing for the Hessian is often more costly than exact evaluation, and there is no saving in solving for the search direction.

The second approach is to use a quasi-Newton method. The basic idea is to use first derivative information gathered from past steps of the algorithm to accumulate an approximation to the Hessian at the current best point. In principle, this approximation can be used exactly like the exact Hessian in Newton’s method, but if structured carefully, such an algorithm can also work directly on an approximation to the inverse of the Hessian, thereby reducing the cost of calculating the actual step. It is also possible to ensure that the approximate Hessian is always positive definite.

Quasi-Newton methods were invented in the mid 1950’s by W.C. Davidon (a physicist). In the mathematical equivalent of not signing the Beatles, his paper on the method was rejected (it was eventually published in 1991). There are now many varieties of Quasi-Newton method, but the most popular is the BFGS method ² , which will be briefly covered here.

We suppose that at step $k$ we have already computed ${\bm \theta }^{[k]}$ and an approximate Hessian ${\bf H}^{[k]}$. So, just as with a regular Newton method, we solve

to obtain a new search direction, leading either to a new ${\bm \theta }^{[k+1]}$ directly, or more likely via line search. We now have the problem of how to appropriately update the approximate Hessian, ${\bf H}^{[k+1]}$. We want the approximate Hessian to give us a good quadratic approximation of our function about ${\bm \theta }^{[k+1]}$. If ${\bf H}^{[k+1]}$ is our new approximate positive definite Hessian at the $(k+1)^{\rm th}$ step, we have the function approximation

The basic requirement of a quasi Newton method is that the new Hessian should be chosen so that the approximation should exactly match $\nabla f({\bm \theta }^{[k]})$ — i.e. the approximation should get the gradient vector at the previous best point, ${\bm \theta }^{[k]}$, exactly right. Taking the gradient of the above approximation gives

Evaluating this at ${\bm \theta }^{[k]}$ leads to our requirement

and hence

where ${\bf s}_ k = {\bm \theta }^{[k+1]} - {\bm \theta }^{[k]}$ and ${\bf y}_ k = \nabla f({\bm \theta }^{[k+1]}) - \nabla f({\bm \theta }^{[k]})$. The so-called secant equation (3.2) will only be feasible for positive definite ${\bf H}^{[k+1]}$ under certain conditions on ${\bf s}_ k$ and ${\bf y}_ k$, but we can always arrange for these to be met by appropriate step length selection (see Wolfe conditions, Nocedal & Wright, 2006, §3.1). Note, in particular, that pre-multiplying (3.2) by ${\bf s}_ k^{\sf T}$ gives ${\bf s}_ k^{\sf T}{\bf y}_ k = {\bf s}_ k^{\sf T}{\bf H}^{[k+1]}{\bf s}_ k>0$ for positive definite Hessian, so ${\bf s}_ k^{\sf T}{\bf y}_ k$ must be positive. Still, the secant equation is under-determined, so we need to constrain the problem further in order to identify an appropriate unique solution. For quasi Newton methods, this is typically done by proposing a low-rank update to the previous Hessian. Methods exist which use a rank 1 update, but BFGS uses a rank 2 update of the form

For (good) reasons that we don’t have time to explore here, the choice ${\bf u} = {\bf y}_ k$, ${\bf v} = {\bf H}^{[k]} {\bf s}_{k}$ is made, and substituting this update into the secant equation allows solving for $\alpha$ and $\beta$, giving the BFGS update,

where $\rho _ k = 1/{\bf s}_ k^{\sf T}{\bf y}_ k$. By substitution, we can directly verify that this satisfies the secant equation. Now let’s work in terms of the inverse approximate Hessian, ${\bf B}^{[k]} \equiv \left({\bf H}^{[k]}\right)^{-1}$, since this will allow us to avoid solving a linear system at each iteration. A naive re-write of our update would be

but this is a low-rank update of an inverse, and hence a good candidate for application of Woodbury’s formula. It is actually a bit tedious to apply Woodbury directly, but if we can guess the solution,

we can verify that it is correct by direct multiplication. We can also double-check that this update satisfies the modified secant equation, ${\bf B}^{[k+1]} {\bf y}_ k = {\bf s}_ k$ (it does). The BFGS update is usually given in the form above, but for actual implementation on a computer, it can be expanded as

and then some careful thinking about the order and re-use of computational steps within the update can lead to a very efficient algorithm.

The BFGS method then works exactly like Newton’s method, but with ${\bf B}^{[k]}$ in place of the inverse of $\nabla ^2 f({\bm \theta }^{[k]})$, and without the need for second derivative evaluation or for perturbing the Hessian to achieve positive definiteness. It also allows direct computation of the search direction via

A finite difference approximation to the Hessian is often used to start the method. The following figure illustrates its application (showing every second step of 38).

Notice how this quasi-Newton method has taken nearly twice as many steps as the Newton method, for the same problem, and shows a greater tendency to ‘zig-zag’. However, it is still computationally cheaper than finite differencing for the Hessian, which would have required an extra 20 gradient evaluations, compared to quasi-Newton. Compared to steepest descent it’s miraculous.

3.4.5 The Nelder–Mead polytope method

What if even gradient evaluation is too taxing, or if our objective is not smooth enough for Taylor approximations to be valid (or perhaps even possible)? What can be done with function values alone? The Nelder–Mead polytope ³ method provides a rather successful answer, and as with the other methods we have met, it is beautifully simple.

Let $n$ be the dimension of $\bm \theta$. At each stage of the method we maintain $n+1$ distinct $\bm \theta$ vectors, defining a polytope in $\Theta$ (e.g. for a 2-dimensional $\bm \theta$, the polytope is a triangle). The following is iterated until a minimum is reached/ the polytope collapses to a point.

1. 1.

   The search direction is defined as the vector from the worst point (vertex of the polytope with highest objective value) through the average of the remaining $n$ points.

2. 2.

   The initial step length is set to twice the distance from the worst point to the centroid of the others. If it succeeds (meaning that the new point is no longer the worst point) then a step length of 1.5 times that is tried, and the better of the two accepted.

3. 3.

   If the previous step did not find a successful new point then step lengths of half and one and a half times the distance from the worst point to the centroid are tried.

4. 4.

   If the last two steps failed to locate a successful point then the polytope is reduced in size, by linear rescaling towards the current best point (which remains fixed.)

5. 5.

   When a point is accepted, the current worst point is discarded, to maintain $n+1$ points, and the algorithm returns to step 1.

Variations are possible, in particular with regard to the step lengths and shrinkage factors.

The following figure illustrates the polytope method in action. Each polytope is plotted, with the line style cycling through, black, grey and dashed black. The worst point in each polytope is highlighted with a circle.

In this case it took 79 steps to reach the minimum of Rosenbrock’s function. This is somewhat more than Newton or BFGS, but given that we need no derivatives in this case, the amount of computation is actually less.

On the basis of this example you might be tempted to suppose that Nelder-Mead is all you ever need, but this is generally not the case. If you need to know the optimum very accurately then Nelder-Mead will often take a long time to get an answer that Newton based methods would give very quickly. Also, the polytope can get ‘stuck’, so that it is usually a good idea to restart the optimisation from any apparent minimum (with a new polytope having the apparent optimum as one vertex), to check that further progress is really not possible. Essentially, Nelder-Mead is good if the answer does not need to be too accurate, and derivatives are hard to come by, but you wouldn’t use it as the underlying optimiser for general purpose modelling software, for example.

3.4.6 Other methods

These notes can only provide a brief overview of some important methods. There are many other optimisation methods, and quite a few are actually useful. For example:

• In likelihood maximisation contexts the method of scoring is often used. This replaces the Hessian in Newton’s method with the expected Hessian.

• The Gauss-Newton method is another way of getting a Newton type method when only first derivatives are available. It is applicable when the objective is somehow related to a sum of squares of differences between some data and the model predictions of those data, and has a nice interpretation in terms of approximating linear models. It is closely related to the algorithm for fitting GLMs.

• Conjugate gradient methods are another way of getting a good method using only first derivative information. The problem with steepest descent is that if you minimise along one steepest descent direction then the next search direction is always at right angles to it. This can mean that it takes a huge number of tiny zig-zagging steps to minimise even a perfectly quadratic function. If $n$ is the dimension of $\bm \theta$, then Conjugate Gradient methods aim to construct a set of $n$ consecutive search directions such that successively minimising along each will lead you exactly to the minimum of a quadratic function.

• Simulated annealing is a gradient-free method for optimising difficult objectives, such as those arising in discrete optimisation problems, or with very spiky objective functions. The basic idea is to propose random changes to ${\bm \theta }$, always accepting a change which reduces $f({\bm \theta })$, but also accepting moves which increase the objective, with a probability which decreases with the size of increase, and also decreases as the optimisation progresses.

• Stochastic gradient descent (SGD) is an increasingly popular approach to optimisation of likelihoods arising from very large data sets. Each step of SGD is a gradient descent based on a new objective derived from a random subset of the data. The stochastic aspect can speed up convergence and the use of a relatively small subset of the data can speed up each iteration: https://en.wikipedia.org/wiki/Stochastic_gradient_descent

• The EM algorithm is a useful approach to likelihood maximisation when there are missing data or random effects that are difficult to integrate out.