--- title: "Mize" author: "James Melville" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Mize} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE, echo = FALSE, message = FALSE} knitr::opts_chunk$set(echo = TRUE, collapse = TRUE, comment = "#>") library(mize) ``` `mize` is a package for doing unconstrained numerical optimization. This vignette describes the use of the `mize` function, which takes as input a function to optimize, its gradient and a starting point and returns the optimized parameters. First, I'll describe the required structure of the function and gradient to optimize, then the various convergence options that apply no matter what optimization methods are being used. In the last section, I'll describe the methods which are available and any options that apply to them specifically. ## The `fg` list If you look at `stats::optim`, you'll see you pass the function to be optimized and its gradient separately. By contrast, `mize` asks you to bundle them up in a list, which I lazily refer to as an `fg` list. Here, for example, is the venerable Rosenbrock function: ```{r Definining a function and gradient to optimize} rb_fg <- list( fn = function(x) { 100 * (x[2] - x[1] * x[1]) ^ 2 + (1 - x[1]) ^ 2 }, gr = function(x) { c( -400 * x[1] * (x[2] - x[1] * x[1]) - 2 * (1 - x[1]), 200 * (x[2] - x[1] * x[1])) }) ``` For minimizing a simple function, there's not a lot of advantage to this. But for more complex optimizations, there are some advantages: * If `fn` and `gr` require extra parameterization, you can store the extra state in a function that creates and returns the `fg` list, and then `fn` and `gr` can access that state via the created closure. `stats::optim` takes an alternative approach where extra data is passed as `...` to the `optim` function and then internally this is passed to `fn` and `gr` when needed. * For some optimization methods (e.g. those using a strong Wolfe line search), both `fn` and `gr` need to be evaluated for the same value of `par`. If there's a lot of shared work between `fn` and `gr`, this can be quite inefficient. If you provide an optional routine called `fg`, which calculates the function and gradient simultaneously, this will be used by some parts of `mize` in preference to calling `fn` and `gr` separately, which can save a bit of time. A version of `rb_fg` with an optional `fg` function could look like: ```{r A function list with an optional fg item} rb_fg <- list( fn = function(x) { 100 * (x[2] - x[1] * x[1]) ^ 2 + (1 - x[1]) ^ 2 }, gr = function(x) { c( -400 * x[1] * (x[2] - x[1] * x[1]) - 2 * (1 - x[1]), 200 * (x[2] - x[1] * x[1])) }, fg = function(x) { a <- x[2] - x[1] * x[1] b <- 1 - x[1] list( fn = 100 * a ^ 2 + b ^ 2, gr = c( -400 * x[1] * a - 2 * b, 200 * a) ) } ) ``` In the case of the Rosenbrock function, I'm not suggesting you should actually do this. See the vignette on [Metric Multi-Dimensional Scaling](mmds.html) for the sort of scenario where it could make a difference. For the rest of this vignette, we'll just assume you're only providing `fn` and `gr`. ## The parameters to optimize The parameters to optimize as passed as a numeric vector, `par`. As well as using the `rb_fg` list we defined above to describe the Rosenbrock function, we'll also start from the following starting point: ```{r Defining a starting point} rb0 <- c(-1.2, 1) ``` ## Basic options If you don't particularly care about how your function gets optimized, the default setting is to use the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) optimizer, which is a good choice for pretty much anything. You may still need to tweak some options that are independent of the method, like the maximum number of iterations to use, or other convergence criteria. ### 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. ```{r Defaults} res <- mize(rb0, rb_fg) # What were the final parameter values? (should be close to c(1, 1)) res$par # What was the function value at that point (should be close to 0) res$f # How many iterations did it take? res$iter # How many function evaluations? res$nf # How many gradient evaluations? res$ng # Why did the optimization terminate? res$terminate ``` 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](convergence.html) 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`: ```{r Verbose mode} res <- mize(rb0, rb_fg, grad_tol = 1e-3, ginf_tol = 1e-3, max_iter = 10, verbose = TRUE) ``` 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: ```{r Log every 10 iterations} res <- mize(rb0, rb_fg, grad_tol = 1e-3, verbose = TRUE, log_every = 10) ``` 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: ```{r Returning stored progress} res <- mize(rb0, rb_fg, store_progress = TRUE, log_every = 10) res$progress ``` As can be seen, `store_progress` will also obey the `log_every` parameter and not store every iteration of convergence information if requested. ## Optimization Methods The available optimization methods are listed below. I won't be describing the algorithms behind them in any depth. For that, you should look in the book by [Nocedal and Wright](https://doi.org/10.1007/978-0-387-40065-5). Many of these methods can be thought of as coming up with a direction to move in, but with the distance to travel in that direction as being someone else's problem, i.e. a line search method needs to also be used to determine the step size. We'll discuss the available step size methods separately. ### 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](https://doi.org/10.1093/imanum/8.1.141). 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. ```{r Steepest descent} 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. ```{r BFGS} 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`: ```{r BFGS without scaled Hessian} 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. ```{r LBFGS} 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`: ```{r LBFGS without scaled Hessian} 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. ```{r CG with PR+} 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](https://users.clas.ufl.edu/hager/papers/Software/) optimizer): ```{r CG with HZ+} 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. ```{r NAG} 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: ```{r NAG with 100 steps} 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 ``` 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: ```{r NAG with 100 steps and less aggressive momentum} 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 ``` 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](https://hdl.handle.net/1807/36012) 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](https://jlmelville.github.io/mize/nesterov.html). ### 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: ```{r Momentum} 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: ```{r Momentum plot} 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 ``` 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`): ```{r momentum with a switch function} # 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: ```{r momentum with a ramp function} 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. ```{r momentum with nesterov 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: ```{r momentum with nesterov schedule and non-zero q} 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: ```{r momentum with random momentum} 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](https://proceedings.mlr.press/v28/sutskever13.html) and [Bengio and co-workers](https://arxiv.org/abs/1212.0901) 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: ```{r Simplified Nesterov momentum} 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: ```{r Nesterov versus classical momentum} 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 # Best f found for classical momentum resc$f 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"`: ```{r Nesterov momentum with convex approximation} res <- mize(rb0, rb_fg, max_iter = 10, method = "MOM", mom_schedule = "nsconvex", nest_convex_approx = TRUE, mom_type = "nesterov") ``` ## Line Searches All the above methods need a line search to find out how far to move in the gradient descent (i.e. non-momentum) part. ### 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](https://doi.org/10.1145/192115.192132), converted from [Dianne O'Leary's Matlab code](https://www.cs.umd.edu/users/oleary/software/). 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](https://gaussianprocess.org/gpml/code/matlab/doc/), Mark Schmidt's [minimization routines](https://www.cs.ubc.ca/~schmidtm/Software/minFunc.html), and an implementation of [Hager-Zhang](https://doi.org/10.1137/030601880) line search are also available (although they are all quite similar): ```{r other Wolfe line search} 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.: ```{r Line search parameters} res <- mize(rb0, rb_fg, max_iter = 10, method = "CG", cg_update = "HZ+", c2 = 0.5, c1 = 0.1) ``` #### Initial Step Size Guess Another important choice is what step size to start each iteration from. Nocedal and Wright suggest two methods based on the result achieved for the previous iteration, one involving the ratio of the slopes at consecutive iterations, and one involving a quadratic interpolation. By default, for Wolfe line searches other than `"Hager-Zhang"`, the quadratic interpolation method is tried. To try the slope ratio method, set `step_next_init = "slope"`. ```{r Line search with slope ratio} res <- mize(rb0, rb_fg, max_iter = 10, step_next_init = "slope") ``` In the [CG_DESCENT](https://doi.org/10.1145/1132973.1132979) software, Hager and Zhang suggest using a "QuadStep" method. In this case, a small trial step is attempted followed by a quadratic interpolation. If the resulting quadratic is strongly convex (for a 1D quadratic $Ax^2 + Bx + c$, this just means $A > 0$), then the minimizer of the quadratic is used as the initial guess. Otherwise, a multiple of the previous step length is used (in `mize`, the old step length is doubled). This is the default if `line_search = "hz"` is used. To explicitly set it for other methods, use `step_next_init = "hz"` (or `"hager-zhang"`): ```{r Line search with Hager-Zhang QuadStep} res <- mize(rb0, rb_fg, max_iter = 10, step_next_init = "hz", line_search = "mt") ``` For the budget-conscious, be aware that using this method always incurs an extra function evaluation per iteration, even if the quadratic minimizer isn't used as the guess. #### Initial Step Size Guess on First Iteration This only leaves what to do on the very first iteration, where there is no previous iteration to look at. You don't really have much information at that point apart from the gradient, so all the implementations I've seen use a step size that is effectively inversely proportional to a norm of the gradient: * Carl Edward Rasmussen's [conjugate gradient routine](https://gaussianprocess.org/gpml/code/matlab/doc/) uses $\frac{1}{1+\left\| g \right\|_2^2}$ * [SciPy](https://scipy.org/)'s `optimize.py` uses $\frac{1}{\left\| g \right\|_2}$ * Mark Schmidt's [minimization routines](https://www.cs.ubc.ca/~schmidtm/Software/minFunc.html) uses $\frac{1}{1+\left\| g \right\|_1}$ * Hager and Zhang use a slightly more complicated procedure for their [CG_DESCENT](https://doi.org/10.1145/1132973.1132979) software: use $\psi_0\frac{\left\|x\right\|_\infty}{\left\|g\right\|_\infty}$, where $x$ is the initial parameter vector and $\psi_0$ is a small positive value (e.g. 0.01), or if that doesn't produce a finite positive value, use $\psi_0\frac{\left|f\right|}{\left\|g\right\|_{2}^{2}}$ where $\left|f\right|$ is the absolute value of the function evaluated at $x$. And if that doesn't work, just use 1. By default, the Hager-Zhang line search uses the CG_DESCENT method. For other Wolfe line searches I have arbitrarily plumped for the Rasmussen method. If you think it will make a huge difference you can explicitly choose by passing one of `"rasmussen"`, `"scipy"`, `"schmidt"` or `"hz"` to the `step0` argument, e.g.: ```{r Line search with scipy initialization} res <- mize(rb0, rb_fg, max_iter = 10, step0 = "scipy") ``` You may also pass a number to `step0` and it will be used as-is: ```{r Line search with initial step length of 1} # 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") ``` The L-BFGS and BFGS methods will also attempt the "natural" Newton step length of 1, in the way suggested in Nocedal and Wright. If you don't want to do this for some reason, it can be explicitly turned off or on via the `try_newton_step` parameter: ```{r BFGS with no Newton step} res <- mize(rb0, rb_fg, max_iter = 10, method = "BFGS", try_newton_step = FALSE) ``` For every other method this is off by default anyway with no particular reason to turn it on. All Wolfe line searches except `Hager-Zhang` use the strong curvature condition to determine when to terminate the line search. The Hager-Zhang line search uses the standard curvature conditions by default. It also uses an approximation to the Armijo sufficient descent conditions which may prevent premature termination of a line search when the minimizer lies close to the initial step size under some circumstances. Use of these variations on the Wolfe conditions can be applied to any of the Wolfe line searches by supplying the `strong_curvature` and `approx_armijo` options: ```{r alternative Wolfe conditions} # 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) ``` It's not obvious to me that any proof of convergence that exists for a given line search method will hold up when the curvature and Armijo conditions are altered, so be careful when using specifying these options. ### 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: ```{r constant step size} 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](https://www.cs.ubc.ca/~schmidtm/Software/minFunc.html) is used to find the step size. ```{r backtracking with cubic interpolation} 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: ```{r backtracking with halved step size} 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`. ```{r bold driver} # 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: ```{r max line search functions} # No more than 10 gradient evaluations allowed per line search res <- mize(rb0, rb_fg, max_iter = 10, ls_max_gr = 10) ``` ## Adaptive Learning Rate Methods So far, the `method` and `line_search` arguments can be freely mixed and matched. However, some methods choose both the direction and the step size directly. The [Delta-Bar-Delta method](https://doi.org/10.1016/0893-6080(88)90003-2) has its roots in neural network research but was also used in the currently popular [t-Distributed Stochastic Neighbor Embedding](https://lvdmaaten.github.io/tsne/) data visualization method. The way DBD works is to consider each component of the gradient and a weighted average of previous iterates gradients. If the sign has changed for that component, that suggests the minimum has been skipped, so the step size (where each component is initialized with value `step0`) is reduced (by a factor `step_down`) in that direction. Otherwise the step size is increased (by a factor `step_up`). The weighting of the average of previous gradients is controlled by `dbd_theta`, which varies between 0 and 1. The same line search parameters that control the bold driver method can be used in DBD: ```{r} res <- mize(rb0, rb_fg, max_iter = 10, method = "DBD", step0 = "rasmussen", step_down = 0.5, step_up = 1.1, dbd_weight = 0.5) ``` One extra parameter is available: `step_up_fn`. In the original DBD paper, when the new step size is to be increased for a component, it is done by adding the value of `step_up` to the current value of the step size. [Janet and co-workers](https://doi.org/10.1109/IJCNN.1998.687205) noted that when the step sizes are small compared to the value of `step_up` this can lead to far too large a step size increase. Instead, they suggest multiplying by `step_up` so that a percentage increase in the step size is achieved. This is the default behavior in `mize`, but the t-SNE optimization method uses the addition method. To get this behavior, set the `step_up_fun` value to `"+"`: ```{r t-SNE style DBD parameters} res <- mize(rb0, rb_fg, max_iter = 10, method = "DBD", step0 = "rasmussen", step_down = 0.8, step_up = 0.2, step_up_fun = "+") ``` The DBD method can be used with momentum. If it is, then the `dbd_weight` parameter is ignored and instead of storing its own history of weighted gradients, the update vector stored by the momentum routines is used instead. This is not mentioned in the original paper, but is used in the t-SNE implementation. The DBD method is interesting in that it uses no function evaluations at all in its optimization. It doesn't even avail itself of the values of any of the gradient components, only their signs. However, as a result, it's really, really easy for DBD to diverge horribly and quickly. You will therefore want to ensure that you are checking the convergence. This is one of the methods where using the `grad_tol` or `ginf_tol` is a better idea than `abs_tol` and `rel_tol`, because it provides some check on divergence without calculating the function ever iteration only for convergence checking purposes. ```{r} # DBD with rel_tol and abs_tol is explicitly set res <- mize(rb0, rb_fg, max_iter = 10, method = "DBD", step0 = "rasmussen", step_down = 0.8, step_up = 0.2, step_up_fun = "+", rel_tol = 1e-8, abs_tol = 1e-8) # 10 gradient calculations as expected res$ng # But 10 function calculations too, only used in the tolerance check res$nf # Turn off the rel_tol and abs_tol and let max_iter handle termination res <- mize(rb0, rb_fg, max_iter = 10, method = "DBD", step0 = "rasmussen", step_down = 0.8, step_up = 0.2, step_up_fun = "+", rel_tol = NULL, abs_tol = NULL, grad_tol = 1e-5) # 11 gradient calculations res$ng # Only one function evalation needed (to calculate res$f) res$nf ``` Using `grad_tol` results in one extra gradient calculation (because it using gradients for the convergence check), but only one function evaluation. The more iterations this runs for, the better trade off this is, although in the event of a diverging optimization, the "best" `par` that is returned may not be as good as would have been returned if we were also using function evaluations. ## Adaptive Restart [O'Donoghue and Candes](https://arxiv.org/abs/1204.3982) described a restart scheme for NAG as a cure for the non-monotonic convergence. They suggest using either a function-based check (did the function value after the gradient descent stage decrease between iterations?) or a gradient-based check (is the momentum direction a descent direction?). If the check fails, then the momentum is restarted: the previous update vector is discarded and the momentum set back to zero. `mize` allows these checks to be used to restart any momentum scheme, not just NAG. The only change is that the function-based check is carried out at the beginning and end of each iteration, rather than at the end of the gradient descent stage. Restarting is applied using the `restart` argument, supplying either `fn` or `gr` for function-based and gradient-based restarting, respectively: ```{r momentum with restart} resc <- mize(rb0, rb_fg, max_iter = 100, method = "MOM", mom_schedule = 0.9, store_progress = TRUE) resf <- mize(rb0, rb_fg, max_iter = 100, method = "MOM", mom_schedule = 0.9, store_progress = TRUE, restart = "fn") resg <- mize(rb0, rb_fg, max_iter = 100, method = "MOM", mom_schedule = 0.9, store_progress = TRUE, restart = "gr") plot(resc$progress$nf, log(resc$progress$f), type = "l", ylim = range(log(resc$progress$f), log(resf$progress$f), log(resg$progress$f))) lines(resf$progress$nf, log(resf$progress$f), col = "red") lines(resg$progress$nf, log(resg$progress$f), col = "blue") ``` Adaptive restart (red and blue lines) certainly seems to help here, although with either gradient or function-based restart seeming equally effective. The only downside to adaptive restart is that you may want to choose the type of restart to not require more overhead of function and gradient calls than necessary. In practice, the gradient-based restart is economical because in most cases, the gradient calculated at the end of the iteration can be re-used for the gradient descent at the start of the next iteration. The one exception here is when using nesterov momentum, where the gradients aren't calculated in the right location to be of use. In that case, you are better off using function-based restarting, although if you aren't checking convergence with a function calculation or using a line search method that makes use of the function at the starting point, you will obviously incur the overhead of the function call. This is one reason to favour using the NAG method directly rather than emulating it with Nesterov momentum. There is also a `restart_wait` parameter. This determines how many iterations to wait between attempted restarts. By default, using the value from the paper of Su and co-workers, we wait 10 iterations before allowing another restart. Making the wait time too short may lead to premature convergence: ```{r momentum with restart and wait time} resfw <- mize(rb0, rb_fg, max_iter = 100, method = "MOM", mom_schedule = 0.9, store_progress = TRUE, restart = "fn", restart_wait = 1) resgw <- mize(rb0, rb_fg, max_iter = 100, method = "MOM", mom_schedule = 0.9, store_progress = TRUE, restart = "gr", restart_wait = 1) plot(resc$progress$nf, log(resc$progress$f), type = "l", ylim = range(log(resc$progress$f), log(resf$progress$f), log(resg$progress$f), log(resfw$progress$f), log(resgw$progress$f))) lines(resf$progress$nf, log(resf$progress$f), col = "red") lines(resfw$progress$nf, log(resfw$progress$f), col = "blue") lines(resgw$progress$nf, log(resgw$progress$f), col = "orange") ``` Here, reducing the wait time to 1 had a great effect on gradient restarting (the orange line), it substantially outperforms the other methods. On the other hand, function restarting more often (the blue line) led to worse performance. There is also a `"speed"` restart type, related to a restart strategy introduced by [Su and co-workers](https://papers.nips.cc/paper/5322-a-differential-equation-for-modeling-nesterovs-accelerated-gradient-method-theory-and-insights), which restarts if the update vector doesn't increase in length (as measured by Euclidean 2-norm) between iterations. This seems less performant in their experiments than gradient-based restarting, but might provide more stability under some circumstances.