Mize

Defaults

By default, mize will grind away on the function you give it for up to 100 iterations. But it will give up early if it senses that no further progress is being made.

The return value provides both the optimized parameter as well as some information about how much work was done to get there.

res <- mize(rb0, rb_fg)
# What were the final parameter values? (should be close to c(1, 1))
res$par
#> [1] 0.9999972 0.9999923

# What was the function value at that point (should be close to 0)
res$f
#> [1] 4.521673e-10

# How many iterations did it take?
res$iter
#> [1] 37

# How many function evaluations?
res$nf
#> [1] 49

# How many gradient evaluations?
res$ng
#> [1] 49

# Why did the optimization terminate?
res$terminate
#> $what
#> [1] "abs_tol"
#> 
#> $val
#> [1] 1.06609e-08

As you can see from looking at res$iter, the L-BFGS didn’t use all 100 iterations it had available. The difference between consecutive function values, res$f got sufficiently low, that the absolute tolerance criterion was passed: res$terminate$what will contain a string code that communicates what happened, while res$terminate$val provides the specific value that causes the termination, although it’s rarely interesting.

There is an entire vignette on the different convergence criteria available. The defaults should do a reasonable, albeit fairly conservative, job.

Logging progress to console

If you want to keep track of the progress of the optimization, set the verbose option to TRUE:

res <- mize(rb0, rb_fg, grad_tol = 1e-3, ginf_tol = 1e-3, max_iter = 10, 
            verbose = TRUE)
#> 03:39:34 iter 0 f = 24.2 g2 = 232.9 ginf = 215.6 nf = 1 ng = 1 step = 0 alpha = 0
#> 03:39:34 iter 1 f = 19.5 g2 = 200.9 ginf = 185.6 nf = 3 ng = 3 step = 0.02169 alpha = 9.312e-05
#> 03:39:34 iter 2 f = 11.57 g2 = 135.4 ginf = 124.6 nf = 4 ng = 4 step = 0.04729 alpha = 0.3469
#> 03:39:34 iter 3 f = 4.281 g2 = 17.28 ginf = 16.27 nf = 5 ng = 5 step = 0.09809 alpha = 1
#> 03:39:34 iter 4 f = 4.144 g2 = 2.588 ginf = 2.469 nf = 6 ng = 6 step = 0.01426 alpha = 1
#> 03:39:34 iter 5 f = 4.139 g2 = 1.773 ginf = 1.626 nf = 7 ng = 7 step = 0.002565 alpha = 1
#> 03:39:34 iter 6 f = 4.132 g2 = 2.44 ginf = 2.327 nf = 8 ng = 8 step = 0.004716 alpha = 1
#> 03:39:34 iter 7 f = 4.12 g2 = 3.727 ginf = 3.033 nf = 9 ng = 9 step = 0.009086 alpha = 0.4616
#> 03:39:34 iter 8 f = 4.098 g2 = 5.427 ginf = 3.863 nf = 10 ng = 10 step = 0.01724 alpha = 0.2631
#> 03:39:34 iter 9 f = 2.569 g2 = 8.248 ginf = 7.431 nf = 16 ng = 16 step = 0.8311 alpha = 4.743
#> 03:39:34 iter 10 f = 2.553 g2 = 11.67 ginf = 9.996 nf = 18 ng = 18 step = 0.04657 alpha = 0.05068

As you can see, information about the state of the optimization is logged to the console every time the convergence criteria are checked. The values logged are:

• f - The current value of the function.
• g2 - The current value of the gradient 2-norm (if you set a non-NULL grad_tol).
• ginf - The current value of the gradient infinity norm (if you set a non-NULL ginf_tol).
• nf - The total number of function evaluations so far.
• ng - The total number of gradient evaluations so far.
• step - The magnitude of the change between the previous version of par and the current value.

Note that values are also provided for iteration 0, which represents the state of par at its starting point. Because this requires a function evaluation that might otherwise not be needed, you may get results where nf and ng are one larger than you’d otherwise expect if verbose was FALSE.

Logging every iteration might be more output than you want. You can also set the log_every parameter to make logging slightly less noisy:

res <- mize(rb0, rb_fg, grad_tol = 1e-3, verbose = TRUE, log_every = 10)
#> 03:39:34 iter 0 f = 24.2 g2 = 232.9 nf = 1 ng = 1 step = 0 alpha = 0
#> 03:39:34 iter 10 f = 2.553 g2 = 11.67 nf = 18 ng = 18 step = 0.04657 alpha = 0.05068
#> 03:39:34 iter 20 f = 0.5142 g2 = 3.001 nf = 31 ng = 31 step = 0.1319 alpha = 1
#> 03:39:34 iter 30 f = 0.009862 g2 = 3.333 nf = 42 ng = 42 step = 0.03554 alpha = 0.1706
#> 03:39:34 iter 37 f = 4.522e-10 g2 = 0.0009379 nf = 49 ng = 49 step = 0.0002131 alpha = 1

On the assumption that you probably want to see information about the last iteration, this is also reported, even if it wouldn’t normally be reported according to the setting for log_every.

Storing progress

You may wish to do more than stare at the progress numbers logged to the console during optimization when verbose = TRUE. If you would like access to these values after optimization, set the store_progress parameter to true, and a data frame called progress will be part of the return value:

res <- mize(rb0, rb_fg, store_progress = TRUE, log_every = 10)
res$progress
#>               f nf ng         step      alpha
#> 0  2.420000e+01  1  0 0.0000000000 0.00000000
#> 10 2.552598e+00 18 18 0.0465699833 0.05067737
#> 20 5.142086e-01 31 31 0.1319347626 1.00000000
#> 30 9.862134e-03 42 42 0.0355364981 0.17057755
#> 37 4.521673e-10 49 49 0.0002130984 1.00000000

As can be seen, store_progress will also obey the log_every parameter and not store every iteration of convergence information if requested.

Steepest Descent ("SD")

The good news is that you always get some improvement with steepest descent. The bad news is everything else: it often takes a long time to make any progress, unless you are very clever with the step size. It is however, the basis for most other gradient descent techniques, and many of them will revert back to this direction, at least temporarily, if things start going wrong.

res <- mize(rb0, rb_fg, max_iter = 10, method = "SD")

Broyden-Fletcher-Goldfarb-Shanno ("BFGS")

The BFGS method is a quasi-Newton method: it constructs an approximation to the inverse of the Hessian. This saves on the annoyance of working out what the second derivatives are for your problem, as well as on the time it takes to calculate them on every iteration, but not on the storage cost, which is O(N^2). It’s therefore not appropriate for large problems.

res <- mize(rb0, rb_fg, max_iter = 10, method = "BFGS")

By default, after the first iteration, the inverse Hessian approximation is scaled to make it more similar to the real inverse Hessian (in terms of eigenvalues). If you don’t want to do this, set the scale_hess parameter to FALSE:

res <- mize(rb0, rb_fg, max_iter = 10, method = "BFGS", scale_hess = FALSE)

Limited-Memory BFGS ("L-BFGS")

The L-BFGS method builds the inverse Hessian approximation implicitly, yielding something close to the BFGS direction without ever directly creating a matrix. This makes it suitable for larger problems. The closeness to the BFGS result is determined by the size of memory parameter, which determines the number of previous updates that are used to create the direction for the next iteration. The higher this value, the closer to the BFGS result will be obtained, but the more storage is required. In practice, the recommended memory size is between 3 and 7. The default is 5.

res <- mize(rb0, rb_fg, max_iter = 10, method = "L-BFGS", memory = 7)

The L-BFGS method also applies a similar Hessian scaling method as with the BFGS method, but does so at every iteration. If you don’t want it for some reason, you can also set scale_hess to FALSE:

res <- mize(rb0, rb_fg, max_iter = 10, method = "L-BFGS", scale_hess = FALSE)

Conjugate Gradient ("CG")

When someone is pulling a decent default optimization off the shelf, if they’re not reaching for L-BFGS, they’re probably going for a conjugate gradient method. In their favor, they require less storage than the L-BFGS. However, they tend to be a lot more sensitive to poor scaling of the function (i.e. when changing different parameters causes very different changes to the objective function) and so need tighter line search criteria to make progress.

There are in fact multiple different methods out there under the name ‘conjugate gradient’, which differ in how the direction is determined from the current and previous gradient and the previous direction. The default method in mize is the “PR+” method (sometimes referred to as “PRP+”), which refers to its inventors, Polak and Ribiere (the third “P” is Polyak who independently discovered it). The “+” indicates that a restart to steepest descent is included under certain conditions, which prevents the possibility of convergence failing.

res <- mize(rb0, rb_fg, max_iter = 10, method = "CG")

Other CG update methods are available via the cg_update parameter. Methods which seem competitive or even superior to the "PR+" update is the "HS+" method (named after Hestenes and Stiefel) and the "HZ+" update (named after Hager and Zhang, and part of the CG_DESCENT optimizer):

res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", cg_update = "HZ+")

Nesterov Accelerated Gradient ("NAG")

The Nesterov Accelerated Gradient method consists of doing a steepest descent step immediately followed by a step that resembles momentum (see below), but with the previous gradient term replaced with the current gradient.

res <- mize(rb0, rb_fg, max_iter = 10, method = "NAG")

Be aware that the convergence of NAG is not guaranteed to be monotone, that is you may (and probably will) see the function value increase between iterations. Getting the optimum convergence rates guaranteed by theory requires you to be able to estimate how strong the convexity of the function you are minimizing is.

To see this in action, let’s increase the iteration count to 100 and plot the progress of the function:

res <- mize(rb0, rb_fg, max_iter = 100, method = "NAG", store_progress = TRUE)
plot(res$progress$nf, log(res$progress$f), type = "l")

res$f
#> [1] 0.04420106

Yikes. If you let this optimization continue, it does eventually settle down, but it’s a real test of patience. Clearly, in this case, there’s probably a better estimate of the strong convexity parameter.

You can control the strong convexity estimate with the nest_q parameter, which is inversely related to the amount of momentum that is applied. It should take a value between 0 (you think your function is strongly convex) and 1 (method becomes steepest descent). The default is zero, and this is in fact the assumed value in most applications and papers.

Let’s try turning down the momentum a little bit and repeat:

resq <- mize(rb0, rb_fg, max_iter = 100, method = "NAG", nest_q = 0.001, 
            store_progress = TRUE)
plot(res$progress$nf, log(res$progress$f), type = "l",
     ylim = range(log(res$progress$f), log(resq$progress$f)))
lines(resq$progress$nf, log(resq$progress$f), col = "red")

resq$f
#> [1] 1.147505

The red line is the result with nest_q = 0.001. That’s better. Unfortunately, I have no advice on how you should go about estimating the right value for nest_q. However, see the section on Adaptive Restart for ways to ameliorate the problem.

To get the theoretical guarantees about convergence, you actually need to follow a fairly strict formula about what the step size should be (it’s related to the reciprocal of the Lipschitz constant of the function), and Sutskever notes that for non-convex functions (with sparse stochastic gradient) the theoretically optimal momentum and step size values are too aggressive. Alas, the best advice (at least for deep learning applications) given was “manual tuning”.

There are connections between NAG and classical momentum, and there is some evidence (at least in deep learning) that replacing classical momentum with the NAG method, but using the same momentum schedule (see below) can improve results. I have written more (much more) on the subject here.

Classical Momentum ("MOM")

Momentum methods use steepest descent and then add a fraction of the previous iteration’s direction. Often the effectiveness of this approach is likened to inertia and has the effect of damping the oscillation in direction that slows the progress of steepest descent with poorly-scaled functions.

Originally popular for training neural networks, the momentum value is often a constant and set to a fairly high value, like 0.9. You can set a constant momentum by setting the mom_schedule parameter:

res <- mize(rb0, rb_fg, max_iter = 10, method = "MOM", mom_schedule = 0.9)

Momentum-based methods share the same issues as Nesterov Accelerated Gradient: convergence need not be monotone. The plot below demonstrates:

res <- mize(rb0, rb_fg, max_iter = 100, method = "MOM", mom_schedule = 0.9,
            store_progress = TRUE)
plot(res$progress$nf, log(res$progress$f), type = "l")

res$f
#> [1] 2.414765

To be fair, there is little reason to believe that the sort of settings used for neural network training would be effective for minimizing the Rosenbrock function.

Non-constant momentum schedules are also available, by passing a string to mom_schedule. These include:

• A switch function, that steps the momentum from one value (mom_init), to another (mom_final) at a specified iteration (mom_switch_iter):

# Switch from a momentum of 0.4 to 0.8 at iteration 5
res <- mize(rb0, rb_fg, max_iter = 10, method = "MOM", mom_schedule = "switch",
            mom_init = 0.4, mom_final = 0.8, mom_switch_iter = 5)

• A ramp function, the linearly increases the momentum from mom_init to mom_final over max_iter iterations:

res <- mize(rb0, rb_fg, max_iter = 10, method = "MOM", mom_schedule = "ramp",
            mom_init = 0.4, mom_final = 0.8)

• The schedule used in NAG can be used here by asking for a "nsconvex" schedule.

res <- mize(rb0, rb_fg, max_iter = 10, method = "MOM", mom_schedule = "nsconvex")

• The nest_q parameter can also be used here as it can be with NAG:

res <- mize(rb0, rb_fg, max_iter = 10, method = "MOM", mom_schedule = "nsconvex",
            nest_q = 0.001)

If none of this meets your needs, you can pass a function that takes the current iteration number and the maximum number of iterations as input and returns a scalar, which will be used as the momentum value, although the value will be clamped between 0 and 1.

For example, you could create a momentum schedule that randomly picks a value between 0 and 1:

mom_fn <- function(iter, max_iter) {
  runif(n = 1, min = 0, max = 1)
}
res <- mize(rb0, rb_fg, max_iter = 10, method = "MOM", mom_schedule = mom_fn)

I’m not saying it’s a good idea.

Simplified Nesterov Momentum

Interest in NAG in the neural network community (specifically deep learning), was sparked by papers from Sutskever and co-workers and Bengio and co-workers who showed that it could be related to classical momentum in a way that made it easy to implement if you already had a working classical momentum implementation, and more importantly, that it seemed to give improvements over classical momentum.

As mize already offers NAG by setting method = "NAG", being able to use it in a different form related to classical momentum may seem superfluous. But in fact, the use of what’s often called “NAG” in deep learning is slightly more general (and as far as I’m aware, has no convergence theory to go with it, so improvements are not guaranteed). It’s probably better to refer to it as (simplified) Nesterov momentum.

At any rate, the key observation of Sutskever was that NAG could be considered a momentum method where the momentum stage occurs before the gradient descent stage. In this formulation, you can have whatever momentum schedule you like, and the hope is that the gradient descent stage can more easily correct if the momentum schedule is too aggressive because it gets to “see” the result of the momentum update and hence the descent occurs from a different location.

To get Nesterov momentum, you need only add a mom_type = "nesterov" parameter to a momentum optimizer:

res <- mize(rb0, rb_fg, max_iter = 10, method = "MOM", mom_schedule = 0.9, 
            mom_type = "nesterov")

I know you’re itching to see whether Nesterov momentum can calm down the bad case of non-monotone convergence we got with a large classical momentum value on the Rosenbrock function:

resc <- mize(rb0, rb_fg, max_iter = 100, method = "MOM", mom_schedule = 0.9, 
             store_progress = TRUE)
resn <- mize(rb0, rb_fg, max_iter = 100, method = "MOM", mom_schedule = 0.9, 
             mom_type = "nesterov", 
             store_progress = TRUE)
# Best f found for Nesterov momentum
resn$f
#> [1] 1.069263
# Best f found for classical momentum
resc$f
#> [1] 2.414765
plot(resc$progress$nf, log(resc$progress$f), type = "l",
     ylim = range(log(resc$progress$f), log(resn$progress$f)))
lines(resn$progress$nf, log(resn$progress$f), col = "red")

The Nesterov momentum is in red. And yes, it’s a lot smoother. On the other hand, you can see that the cost function increases and then levels off. Clearly, some tweaking of the correct momentum schedule to use is needed.

Sutskever also came up with a simplified approximation to the momentum schedule used in NAG. It doesn’t have a big effect on the run time, only applies for the default nest_q of 0, and at early iterations the momentum is much larger. To use it, set the nest_convex_approx = TRUE when using the mom_schedule = "nesterov":

res <- mize(rb0, rb_fg, max_iter = 10, method = "MOM", 
            mom_schedule = "nsconvex", nest_convex_approx = TRUE, 
            mom_type = "nesterov")

Wolfe Line Search

The default for all methods mentioned so far is a Wolfe line search that satisfies the strong Wolfe conditions. Specifically, More-Thuente line search, converted from Dianne O’Leary’s Matlab code. This uses cubic interpolation to find a point close to a minimize along the search direction for each iteration. This requires a function and gradient evaluation at each candidate point.

If the More-Thuente line search fails to work for some reason, alternative Wolfe line search method from Carl Edward Rasmussen’s conjugate gradient routine, Mark Schmidt’s minimization routines, and an implementation of Hager-Zhang line search are also available (although they are all quite similar):

res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", line_search = "Rasmussen")
# Use Mark Schmidt's minFunc line search 
res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", line_search = "Schmidt")
# Hager-Zhang line search  
res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", line_search = "Hager-Zhang")
# Hager-Zhang can be abbreviated to "HZ"
res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", line_search = "HZ")
# You can explicitly set More-Thuente too
res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", line_search = "More-Thuente")
# More-Thuente can be abbreviated to "MT"
res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", line_search = "MT")

A variety of interpolation methods are available in the Matlab minfunc code, but mize only uses the default: use cubic interpolation and extrapolation methods, falling back to back-tracking Armijo search (also using cubic interpolation) if a non-finite value is encountered.

Default values for the line search tolerance use the advice given in the Nocedal and Wright book: the c1 parameter is set to 1e-4, while c2 is set to 0.9 for quasi-Newton methods (BFGS and L-BFGS) and 0.1 for everything else. For a line search involving the strong Wolfe conditions, the c2 parameter measures how strict the line search is: the smaller it is (it cannot be set smaller than c1 or greater than 1), the tighter the line search. This is an important property for some conjugate gradient methods to maintain their convergence properties.

You may be able to get away with a looser line search (smaller c2) for some methods, e.g.:

res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", cg_update = "HZ+", 
            c2 = 0.5, c1 = 0.1)

res <- mize(rb0, rb_fg, max_iter = 10, step_next_init = "slope")

res <- mize(rb0, rb_fg, max_iter = 10, step_next_init = "hz", 
            line_search = "mt")

res <- mize(rb0, rb_fg, max_iter = 10, step0 = "scipy")

# An initial guess of 1 for the step length isn't bad for L-BFGS
res <- mize(rb0, rb_fg, max_iter = 10, step0 = 1, method = "L-BFGS")

res <- mize(rb0, rb_fg, max_iter = 10, method = "BFGS", try_newton_step = FALSE)

# Rasmussen line search with standard Wolfe conditions
res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", line_search = "Rasmussen",
            strong_curvature = FALSE)
# Hager-Zhang with strong Wolfe conditions
res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", line_search = "HZ",
            strong_curvature = TRUE, approx_armijo = FALSE)
# More-Thuente with approx Armijo conditions
res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", line_search = "MT",
            approx_armijo = TRUE)

Other Line Searches

The Wolfe line search is your best bet, but there are some alternatives (of somewhat dubious value).

You can set a constant value for the step size. But note that the total step size is actually determined by the product of the step size and the length of the gradient vector, so you may want to set norm_direction too:

res <- mize(rb0, rb_fg, max_iter = 10, method = "SD", line_search = "constant",
            norm_direction = TRUE, step0 = 0.01)

There is also a backtracking line search, which repeatedly reduces the step size starting from step0, until the sufficient decrease condition (also known as the Armijo condition) of the Wolfe line search conditions is fulfilled. This is controlled by c1, the larger the value, the larger the decrease in the function value is needed for a step size to be accepted. By default, the cubic interpolation method from minFunc is used to find the step size.

res <- mize(rb0, rb_fg, max_iter = 10, line_search = "backtracking", step0 = 1, 
            c1 = 0.1)

You will often see a simple fixed step size reduction used in combination with this line search, mainly due to ease of implementation. Cubic interpolation is probably always preferable, but if you want to use a constant step size reduction, provide a non-NULL value to step_down. The step size will be multiplied by this value so it should be a value between zero and one. Set it to 0.5, for example, to halve the step size while backtracking:

res <- mize(rb0, rb_fg, max_iter = 10, line_search = "backtracking",
            step0 = 1, c1 = 0.1, step_down = 0.5)

A similar approach is the "bold driver", which also backtracks, but starts from the step size found for the last iteration, increased by a factor step_up.

# increase step size by 10%, but reduce by 50%
res <- mize(rb0, rb_fg, max_iter = 10, line_search = "bold",
            step0 = 1, step_down = 0.5, step_up = 1.1)

Unlike backtracking line search, the bold driver accepts any function decrease, and does not use c1.

Maximum function or gradient evaluations per line search

It’s a good idea to specify a maximum number of function evaluations that can be spent in any one line search. This guards against the possibility of poorly chosen parameters, badly behaved functions or other numerical issues causing the line search to fail to converge. By default, no more than 20 function evaluations are allowed per line search. The ls_max_fn, ls_max_gr and ls_max_fg can be used to control this:

# No more than 10 gradient evaluations allowed per line search
res <- mize(rb0, rb_fg, max_iter = 10, ls_max_gr = 10)