SNewton: safeguarded Newton methods for function minimization

In snewton: Safeguarded Newton Approach for Function Minimization

Safeguarded Newton algorithms

So-called Newton methods are among the most commonly mentioned in the solution of nonlinear equations or function minimization. However, as discussed in https://en.wikipedia.org/wiki/Newton%27s_method#History, the Newton or Newton-Raphson method as we know it today was not what either of its supposed originators knew.

This vignette discusses the development of simple safeguarded variants of the Newton method for function minimization in R. Note that there are some resources in R for solving nonlinear equations by Newton-like methods in the packages nleqslv and pracma.

The basic approach

If we have a function $f(x)$, with gradient $g(x)$ and second derivative (Hessian) $H(x)$ the first order condition for an extremum (min or max) is

$g(x) = 0$

To ensure a minimum, we want

$ H(x) > 0 $

The first order condition leads to a root-finding problem.

It turns out that $x$ need not be a scalar. We can consider it to be a vector of parameters to be determined. This renders $g(x)$ a vector also, and $H(x)$ a matrix. The conditions of optimality then require a zero gradient and positive-definite Hessian.

The Newton approach to such equations is to provide a guess to the root $x_try$ and to then solve the equation

$ H(x_t) * s = - g(x_t)$

for the search vector $s$. We update $x_t$ to $x_t + s$ and repeat until we have a very small gradient $g(x_t)$. If $H(x)$ is positive definite, we have a reasonable approximation to a (local) minimum.

Motivations

A particular interest in Newton-like methods its theoretical quadratic convergence. See https://en.wikipedia.org/wiki/Newton%27s_method. That is, the method will converge in one step for a quadratic function $f(x)$, and for "reasonable" functions will converge very rapidly. There are, however, a number of conditions, and practical programs need to include safequards against mis-steps in the iterations.

The principal issues concern the possiblity that $H(x)$ may not be positive definite, at least in some parts of the domain, and that the curvature may be such that a unit step $x_t + s$ does not reduce the function $f$. We therefore get a number of possible variants of the method when different possible safeguards are applied.

Algorithm possibilities

There are many choices we can make in building a practical code to implement the ideas above. In tandem with the two main issues expressed above, we will consider

• the modification of the solution of the main equation $ H(x_t) * s = - g(x_t)$ so that a reasonable search vector $s$ is always generated by avoiding Hessian matrices that are not positive definite.

• the selection of a new set of parameters $x_new = x_t + step * s$ so that the function value $f(x_new)$ is less than $f(x_t)$.

The second choice above could be made slightly more stringent so that the Armijo (??ref) condition of sufficient-decrease is met. Adding a curvature requirement gives the Wolfe condisions. See https://en.wikipedia.org/wiki/Wolfe_conditions. The Armijo requirement is generally written

$f(x_t + steps) < f(x_t) + c * step * g(x_t)^Ts$

where c is some number less than 1. Typically $ c = 1e-4 = 0.0001 $. Note that the product of gradient times search vector is negative for any reasonable situation, since we are trying to go "downhill".

As a result of the ideas in this section, the code snewton() uses a solution of the Newton equations with the Hessian provided (if this is possible, else we stop), along with a backtracking line search. The code snewtonm uses a Marquardt stabilization of the Hessian to create

$ Haug = H + 1_n * lambda$

That is, we add lambda times the unit matrix to $H$. Then we try the set of parameters found by adding the solution of the Newton equations with $Haug$ in place of $H$ to the current "best" set of parameters. If this new set of parameters has a higher function value than the "best" so far, we increase lambda and try again. Note that we do not need to re-evaluate the gradient or Hessian to do this. Moreover, for some value of lambda, the step is clearly down the gradient (i.e., steepest descents) or we have converged and no progress is possible. This leads to a very compact and elegant code. It is reliable, but may be less efficient than using the un-modified Hessian.

A choice to compute the search vector

The primary concern in solving for $s$ is that the Hessian may not be positive definite. This means that we cannot apply fast and stable methods like the Cholesky decomposition to the matrix. At the time of writing, we use the following approach:

• We attempt to solve $ H(x_t) * s = - g(x_t)$ with R directly, and rely on internal checks to catch any cases where the solution fails. We then use try() to stop the program in this case.

Choosing the step size

The traditional Newton approach is that the stepsize is taken to be 1. In practice, this can sometimes mean that the function value is not reduced. As an alternative, we can use a simple backtrack search. We start with $step = 1$ (actually the program allows for the element defstep of the control list to be set to a value other than 1). If the Armijo condition is not met, we replace $step$ with $ r * step $ where $r$ is less than 1. Here we suggest control$stepdec = 0.2. We repeat until $x_t$ satisfies the Armijo condition or $x_t$ is essentially unchanged by the step.

Here "essentially unchanged" is determined by a test using an offset value, that is, the test

$ (x_t + offset) == (x_t + step * d + offset) $

where $d$ is the search direction. control$offset = 100 is used. We could also, and almost equivalently, use the R identical function.

Examples

These examples are coded as test to the package snewton.

A simple example

The Rosenbrock function

The Wood function

Note that we have NOT found the minimum for the Wood function.

snewton documentation built on May 2, 2019, 6:47 p.m.