Gradient (first-order)

For a scalar-valued function \(f : \mathbb{R}^n \to \mathbb{R}\), the gradient is the vector of partial derivatives,

\[ \nabla f(x) = \begin{bmatrix} \frac{\partial f}{\partial x_1} \\ \vdots \\ \frac{\partial f}{\partial x_n} \end{bmatrix} \]

• Points in the direction of steepest ascent at the current location \(x\)

• Its magnitude \(\|\nabla f(x)\|\) gives the rate of steepest ascent

• At a (local) minimum, \(\nabla f(x^*) = 0\) (also true for maximum / saddle points)

Gradient - examples

\[ f(x, y) = x^2 + 3xy + y^2 \]

Hessian (second-order)

The Hessian is the \(n \times n\) matrix of second partial derivatives,

\[ \nabla^2 f(x_1,\ldots,x_n) = \begin{bmatrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_1 \partial x_n} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} & \cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_n^2} \end{bmatrix} \]

• Symmetric when \(f\) has continuous second partial derivatives

• Describes the local curvature of \(f\) - how the gradient changes

Hessian - examples

\[ f(x, y) = x^2 + 3xy + y^2 \]

Naive Gradient Descent

The basic idea behind this approach is that the gradient of a function tells us the direction of steepest ascent (or descent). Therefore, to find the minimum we should take our next step in the direction of the negative gradient to most quickly approach the nearest minima.

Given an \(n\)-dimensional function \(f(x_1, \ldots, x_n)\), and an current position \(x_k\) then our update rule is,

\[ x_{k+1} = x_{k} - \alpha \nabla f(x_k) \]

here \(\alpha\) refers to the step length or the learning rate which determines how big a step we will take.

Implementation

grad_desc_1d = function(x, f, grad, step, max_step = 100, tol = 1e-6) {
  res = list(x = x, f = f(x))

  for (i in seq_len(max_step)) {
    x_new = x - grad(x) * step

    if (!is.finite(x_new) || !is.finite(f(x_new))) {
      warning("Non-finite value encountered! Try a smaller step size.")
      break
    }
    if (abs(x_new - x) < tol)
      break

    x = x_new
    res$x = c(res$x, x)
    res$f = c(res$f, f(x))
  }

  if (i == max_step)
    warning("Failed to converge!")

  res
}

A basic example

f = \(x) x^2
grad = \(x) 2 * x

opt = grad_desc_1d(-2, f, grad, step = 0.25)
plot_1d_traj(c(-2, 2), f, opt)

opt = grad_desc_1d(-2, f, grad, step = 0.5)
plot_1d_traj(c(-2, 2), f, opt)

Where can it go wrong?

If you pick a bad step size …

opt = grad_desc_1d(-2, f, grad, step = 0.9)
plot_1d_traj(c(-2, 2), f, opt)

opt = grad_desc_1d(-2, f, grad, step = 0.01)
plot_1d_traj(c(-2, 2), f, opt)

Local minima

The function below has multiple minima, both starting point and step size affect the solution we obtain,

f = \(x) x^4 + x^3 - x^2 - x
grad = \(x) 4*x^3 + 3*x^2 - 2*x - 1

opt = grad_desc_1d(-1.5, f, grad, step = 0.2)
plot_1d_traj(c(-1.5, 1.5), f, opt)

opt = grad_desc_1d(-1.5, f, grad, step = 0.25)
plot_1d_traj(c(-1.5, 1.5), f, opt)

Alternative starting points

opt = grad_desc_1d(1.5, f, grad, step = 0.2)
plot_1d_traj(c(-1.75, 1.5), f, opt)

opt = grad_desc_1d(1.25, f, grad, step = 0.2)
plot_1d_traj(c(-1.75, 1.5), f, opt)

Problematic step sizes

If the step size is too large it is possible for the algorithm to overflow,

opt = grad_desc_1d(-1.5, f, grad, step = 0.75)

Warning in grad_desc_1d(-1.5, f, grad, step = 0.75): Non-finite value
encountered! Try a smaller step size.

plot_1d_traj(c(-1.5, 1.5), f, opt)

opt$x

[1] -1.500000e+00  2.062500e+00 -2.998608e+01  7.879000e+04
[5] -1.467367e+15  9.478445e+45

opt$f

[1]  9.375000e-01  2.055299e+01  7.806665e+05  3.853806e+19
[5]  4.636118e+60 8.071392e+183

Gradient Descent w/ backtracking

Having too large of a step can be problematic, one solution is to allow the step size to adapt.

Backtracking involves checking if the proposed move is advantageous, i.e.  \[ f(x_k-\alpha \nabla f(x_k)) < f(x_k) \]

• If it is downhill then accept \(x_{k+1} = x_k-\alpha \nabla f(x_k)\).

• If not, adjust \(\alpha\) by a factor \(\tau\) (e.g. 0.5) and check again.

Pick larger \(\alpha\) to start (but not so large so as to overflow) and then let the backtracking tune things.

Implementation

grad_desc_1d_bt = function(x, f, grad, step, tau = 0.5, max_step = 100, max_back = 10, tol = 1e-6) {
  res = list(x = x, f = f(x))

  for (i in seq_len(max_step)) {
    grad_f = grad(x)
    for (j in seq_len(max_back)) {
      x_new = x - step * grad_f
      f_new = f(x_new)
      if (f_new < tail(res$f, 1))
        break
      step = step * tau
    }

    if (abs(x_new - x) < tol)
      break

    x = x_new
    res$x = c(res$x, x)
    res$f = c(res$f, f_new)
  }

  if (i == max_step)
    warning("Failed to converge!")

  res
}

A 2d function

We will be using mk_quad() to create quadratic functions with varying conditioning (as specified by the epsilon parameter).

\[ \begin{align} f(x,y) &= 0.33(x^2 + \epsilon^2 y^2 ) \\ \nabla f(x,y) &= \left[ \begin{matrix} 0.66 \, x \\ 0.66 \, \epsilon^2 \, y \end{matrix} \right] \\ \nabla^2 f(x,y) &= \left[\begin{array}{cc} 0.66 & 0 \\ 0 & 0.66 \, \epsilon^2 \end{array}\right] \end{align} \]

Examples

q = mk_quad(0.7)
plot_2d_traj(
  c(-1, 2), c(-1, 2), q$f, 
  title = "ε = 0.7"
)

q = mk_quad(0.05)
plot_2d_traj(
  c(-1, 2), c(-2, 2), q$f, 
  title = "ε = 0.05"
)

\(n\)-d gradient descent w/ backtracking

grad_desc = function(x, f, grad, step, tau = 0.5, max_step = 100, max_back = 10, tol = 1e-8) {
  res = list(x = list(x), f = f(x))

  for (i in seq_len(max_step)) {
    grad_f = grad(x)

    for (j in seq_len(max_back)) {
      x_new = x - grad_f * step
      f_new = f(x_new)
      if (f_new < tail(res$f, 1))
        break
      step = step * tau
    }

    if (sqrt(sum((x_new - x)^2)) < tol)
      break

    x = x_new
    res$x = c(res$x, list(x))
    res$f = c(res$f, f_new)
  }

  if (i == max_step)
    warning("Failed to converge!")

  res
}

Well conditioned cost function

q = mk_quad(0.7)
opt = grad_desc(c(1.6, 1.1), q$f, q$grad, step = 1)
plot_2d_traj(
  c(-1, 2), c(-1, 2), q$f, traj = opt, title = "ε = 0.7"
)

q = mk_quad(0.7)
opt = grad_desc(c(1.6, 1.1), q$f, q$grad, step = 2)
plot_2d_traj(
  c(-1, 2), c(-1, 2), q$f, traj = opt, title = "ε = 0.7"
)

Ill-conditioned cost function

q = mk_quad(0.05)
opt = grad_desc(c(1.6, 1.1), q$f, q$grad, step = 1)
plot_2d_traj(
  c(-1, 2), c(-2, 2), q$f, traj = opt, title = "ε = 0.05"
)

q = mk_quad(0.05)
opt = grad_desc(c(1.6, 1.1), q$f, q$grad, step = 2)
plot_2d_traj(
  c(-1, 2), c(-2, 2), q$f, traj = opt, title = "ε = 0.05"
)

Rosenbrock’s function

This is another classic ill-conditioned function commonly used to benchmark optimization algorithms. It has a narrow, curved (banana-shaped) valley - easy to find the valley, hard to follow it to the minimum

Classic definition is, \[ f(u, v) = \tfrac{1}{2}(1 - u)^2 + (v - u^2)^2 \]

Our version is scaled anduses \(u = 4x+1\) and \(v = 4y+3\),

\[ f(x, y) = 8x^2 + \big(4y + 3 - (4x+1)^2\big)^2 \]

The minimum is at \((u^*, v^*) = (1, 1)\), \((x^*, y^*) = (0, -\tfrac{1}{2})\) where \(f = 0\)

Rosenbrock function (very ill conditioned)

rb = mk_rosenbrock()

opt = grad_desc(c(1.6, 1.1), rb$f, rb$grad, step = 0.25)
plot_2d_traj(c(-2, 2), c(-2, 2), rb$f, traj = opt)

opt = grad_desc(c(-0.5, 0), rb$f, rb$grad, step = 0.25)
plot_2d_traj(c(-2, 2), c(-2, 2), rb$f, traj = opt)

Limitations of gradient descent

• Converges to a local minima - sensitive to starting location

  • Global convergence guarantees only hold for convex functions

• Requires gradient computation at every step - expensive when \(n\) is large or the gradient has no closed form

• Sensitive to the choice of learning rate - too large diverges, too small converges slowly

• Uses a scalar step size applied uniformly in all directions

  • Performs poorly on ill-conditioned problems where the optimal step varies by direction

Newton’s Method in 1d

Let’s assume we have a 1d function \(f(x)\) we are trying to optimize, our current guess is \(x\) and we want to know how to generate our next step \(\Delta x\).

We start by constructing the 2nd order Taylor approximation of our function at \(x+\Delta x\),

\[ f(x + \Delta x) \approx \widehat{f}(x + \Delta x) = f(x) + \Delta x \, f'(x) + \frac{1}{2} (\Delta x)^2 f''(x) \]

Finding the Newton step

Our optimal step then becomes the value of \(\Delta x\) that minimizes the quadratic given by our Taylor approximation.

Generalizing to \(n\)d

Based on the same argument we can see the following result for a function in \(\mathbb{R}^n\),

\[ f(x + \Delta x) \approx \widehat{f}(x) = f(x) + \Delta x^T \nabla f(x) + \frac{1}{2} \Delta x^T \, \nabla^2 f(x) \,\Delta x \]

Implementation

newtons_method = function(x, f, grad, hess, max_iter = 100, tol = 1e-8) {
  res = list(x = list(x), f = f(x))

  for (i in seq_len(max_iter)) {
    g = grad(x)
    H = hess(x)
    x_new = x - solve(H, g)

    if (sqrt(sum((x_new - x)^2)) < tol)
      break

    x = x_new
    res$x = c(res$x, list(x))
    res$f = c(res$f, f(x))
  }

  res
}

A basic example

f = \(x) x^2
grad = \(x) 2 * x
hess = \(x) matrix(2)

opt = newtons_method(-2, f, \(x) matrix(grad(x)), \(x) hess(x))
plot_1d_traj(c(-2, 2), f, opt)

opt = newtons_method(1, f, \(x) matrix(grad(x)), \(x) hess(x))
plot_1d_traj(c(-2, 2), f, opt)

1d Quartic

f = \(x) x^4 + x^3 - x^2 - x
grad = \(x) 4*x^3 + 3*x^2 - 2*x - 1
hess = \(x) 12*x^2 + 6*x - 2

opt = newtons_method(-1.5, f, \(x) matrix(grad(x)), \(x) matrix(hess(x)))
plot_1d_traj(c(-1.5, 1.5), f, opt)

opt = newtons_method(1.5, f, \(x) matrix(grad(x)), \(x) matrix(hess(x)))
plot_1d_traj(c(-1.5, 1.5), f, opt)

2d quadratic cost functions

q = mk_quad(0.7)
opt = newtons_method(c(1.6, 1.1), q$f, q$grad, q$hess)
plot_2d_traj(c(-1, 2), c(-1, 2), q$f, traj = opt)

q = mk_quad(0.05)
opt = newtons_method(c(1.6, 1.1), q$f, q$grad, q$hess)
plot_2d_traj(c(-1, 2), c(-1, 2), q$f, traj = opt)

Rosenbrock function

rb = mk_rosenbrock()
opt = newtons_method(c(1.6, 1.1), rb$f, rb$grad, rb$hess)
plot_2d_traj(c(-2, 2), c(-2, 2), rb$f, traj = opt)

rb = mk_rosenbrock()
opt = newtons_method(c(-0.5, 0), rb$f, rb$grad, rb$hess)
plot_2d_traj(c(-2, 2), c(-2, 2), rb$f, traj = opt)

Damped / backtracking implementation

newtons_method_damped = function(
  x, f, grad, hess, max_iter = 100, max_back = 10, tol = 1e-8,
  alpha = 0.5, beta = 0.75
) {
  res = list(x = list(x), f = f(x))

  for (i in seq_len(max_iter)) {
    grad_f = grad(x)
    step = -solve(hess(x), grad_f)
    t = 1
    for (j in seq_len(max_back)) {
      if (f(x + t * step) < f(x) + alpha * t * sum(grad_f * step))
        break
      t = t * beta
    }

    x_new = x + t * step

    if (sqrt(sum((x_new - x)^2)) < tol)
      break

    x = x_new
    res$x = c(res$x, list(x))
    res$f = c(res$f, f(x))
  }

  res
}

2d quadratic cost function

q = mk_quad(0.7)
opt = newtons_method_damped(c(1.6, 1.1), q$f, q$grad, q$hess)
plot_2d_traj(c(-1, 2), c(-1, 2), q$f, traj = opt)

q = mk_quad(0.05)
opt = newtons_method_damped(c(1.6, 1.1), q$f, q$grad, q$hess)
plot_2d_traj(c(-1, 2), c(-1, 2), q$f, traj = opt)

Rosenbrock function

rb = mk_rosenbrock()
opt = newtons_method_damped(c(1.6, 1.1), rb$f, rb$grad, rb$hess)
plot_2d_traj(c(-2, 2), c(-2, 2), rb$f, traj = opt)

rb = mk_rosenbrock()
opt = newtons_method_damped(c(-0.5, 0), rb$f, rb$grad, rb$hess)
plot_2d_traj(c(-2, 2), c(-2, 2), rb$f, traj = opt)

Limitations of Newton’s Method

• Requires both the gradient and Hessian - second derivatives may be unavailable or expensive to compute

• Storing the \(n \times n\) Hessian costs \(O(n^2)\) memory and inverting it costs \(O(n^3)\) per step - impractical for large \(n\)

• Can diverge or behave poorly if the Hessian is indefinite (not positive definite) away from the minimum

• Like GD, finds local minima and is sensitive to starting location

• Damping (backtracking) helps with robustness but adds cost and complexity

Method Comparison

R’s optim()

R’s optim() function (from the stats package) provides a set of general-purpose optimization methods:

Nelder-Mead

optim() interface

optim(par, fn, gr = NULL, method = "Nelder-Mead",
      lower = -Inf, upper = Inf, control = list(), hessian = FALSE)

Key arguments:

• par - initial parameter values (numeric vector)
• fn - function to minimize (must take a vector and return a scalar)
• gr - gradient function (optional; if NULL, finite differences are used)
• method - optimization method
• lower, upper - bounds (only for "L-BFGS-B")
• control - list of control parameters (e.g. maxit, reltol, fnscale)
• hessian - if TRUE, return a numerically approximated Hessian

Example - Well-conditioned quadratic

q = mk_quad(0.7)

optim(
  c(1.6, 1.1), q$f, q$grad,
  method = "BFGS"
)

$par
[1]  5.862023e-11 -4.031139e-11

$value
[1] 1.396753e-21

$counts
function gradient 
      18       15 

$convergence
[1] 0

$message
NULL

optim(
  c(1.6, 1.1), q$f,
  method = "BFGS"
)

$par
[1]  5.862023e-11 -4.031139e-11

$value
[1] 1.396753e-21

$counts
function gradient 
      18       15 

$convergence
[1] 0

$message
NULL

Comparing methods - fixed start

cost_funcs = list(
  "Well-cond quad" = mk_quad(0.7),
  "Ill-cond quad" = mk_quad(0.02),
  "Rosenbrock" = mk_rosenbrock()
)

results_fixed = bench::press(
  cost_func = names(cost_funcs) |> forcats::as_factor(),
  {
    fns = cost_funcs[[cost_func]]
    bench::mark(
      "Nelder-Mead"   = optim(c(1.6, 1.1), fns$f, fns$grad, method = "Nelder-Mead"),
      "BFGS"          = optim(c(1.6, 1.1), fns$f, fns$grad, method = "BFGS"),
      "BFGS w/o Grad" = optim(c(1.6, 1.1), fns$f,           method = "BFGS"),
      "CG"            = optim(c(1.6, 1.1), fns$f, fns$grad, method = "CG"),
      "L-BFGS-B"      = optim(c(1.6, 1.1), fns$f, fns$grad, method = "L-BFGS-B"),
      check = FALSE
    )
  }
)

Comparing methods - random start

results_random = bench::press(
  cost_func = names(cost_funcs) |> forcats::as_factor(),
  x0 = runif(200, -2, 2) |> matrix(ncol=2) |> apply(1,c, simplify = FALSE),
  {
    fns = cost_funcs[[cost_func]] 
    x0 = unlist(x0)
    bench::mark(
      "Nelder-Mead"   = optim(x0, fns$f, fns$grad, method = "Nelder-Mead"),
      "BFGS"          = optim(x0, fns$f, fns$grad, method = "BFGS"),
      "BFGS w/o Grad" = optim(x0, fns$f,           method = "BFGS"),
      "CG"            = optim(x0, fns$f, fns$grad, method = "CG"),
      "L-BFGS-B"      = optim(x0, fns$f, fns$grad, method = "L-BFGS-B"),
      check = FALSE,
      iterations = 1
    )
  }
) |>
  select(-x0)

Some general advice

• Having access to the gradient is almost always helpful / necessary

• Having access to the Hessian can be helpful, but usually does not significantly improve things

• The curse of dimensionality is real - everything gets much harder as \(n\) grows

  • Be careful with tolerance settings - what is “close enough” in 2d may not be close enough in 1000d

• In general, BFGS or L-BFGS are a good place to start for most problems (either well- or ill-conditioned)