library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly)      # interactive figures and ALL 3-D graphics
library(numDeriv)    # finite-difference gradients, for checking

How this chapter uses graphics

Every two-dimensional figure is drawn with ggplot2 and rendered statically. Immediately after each one, the equivalent plot_ly() code appears in a chunk marked eval=FALSE, echo=TRUE.

Every three-dimensional figure is drawn with plot_ly() and evaluated. Optimization is a landscape subject: an algorithm’s behavior is a trajectory over a surface, and convergence rates are functions of two parameters at once, condition number and iteration, step size and smoothness, input location and evidence. Rotating those surfaces is how the geometry becomes visible.

1 Learning objectives

After completing this chapter you will be able to:

State first- and second-order optimality conditions, distinguishing necessary from sufficient.
Define $L$-smoothness and $\mu$-strong convexity, and explain how the condition number $\kappa=L/\mu$ governs convergence.
Derive the gradient-descent stability condition $\eta<2/L$ and the convergence rates $O(1/k)$ and $(1-\mu/L)^k$.
Implement backtracking line search and explain the Wolfe conditions.
Explain why Nesterov acceleration achieves $O(1/k^2)$ and why that is optimal for first-order methods.
Implement BFGS from the correct recursion and verify the curvature condition.
State and apply the KKT conditions, and use complementary slackness to explain SVM sparsity.
Construct the Lagrange dual, and state weak duality, Slater’s condition, and the duality gap.
Compute proximal operators, and recognize soft-thresholding as $\operatorname{prox}_{\lambda\|\cdot\|_1}$.
Fit a Gaussian process stably, learn its hyperparameters by marginal likelihood, and run Bayesian optimization with correct acquisition functions.

Estimated time: 14–18 hours including exercises. Prerequisites: Chapter 3 (matrix computing, conditioning), Chapter 6 (the SVM dual, §6.14), Chapter 11 (LASSO and soft-thresholding, §11.7), and Chapter 12 (likelihood-based estimation). This chapter supplies the machinery all three of those rely on.

2 PART I: THE OPTIMIZATION PROBLEM

3 Formulation

The general problem is

\[\boxed{\;\min_{x\in\mathbb R^n}\ f(x)\quad\text{subject to}\quad \begin{cases}g_i(x)\le0,& i=1,\dots,m\\ h_j(x)=0,& j=1,\dots,p\end{cases}\;}\]

with $f$ the objective, $g_i$ inequality constraints, $h_j$ equality constraints, and the feasible set $\mathcal F=\{x: g_i(x)\le0,\ h_j(x)=0\ \forall i,j\}$.

Maximization needs no separate theory: $\arg\max_x f(x)=\arg\min_x -f(x)$. And if $\psi$ is strictly increasing then $\arg\min_x f(x)=\arg\min_x\psi(f(x))$ — which is why maximizing a likelihood and minimizing its negative log are the same problem, and why log-transforms are free.

Problem class	Objective	Constraints	Typical method
Unconstrained smooth	differentiable	none	GD, Newton, BFGS
Unconstrained non-smooth	e.g. $\ell_1$	none	subgradient, proximal
Linear program	linear	linear	simplex, interior point
Quadratic program	convex quadratic	linear	active set, interior point
Convex	convex	convex	interior point, ADMM
General nonlinear	any	any	SQP, augmented Lagrangian
Black-box / expensive	unknown form	box	Bayesian optimization

4 Optimality conditions

Common misconception: “set the gradient to zero and check that the Hessian is positive definite.” That combination is sufficient for a strict local minimum, not necessary. The necessary conditions are weaker:

\[\textbf{Necessary (first order): }\ \nabla f(x^\star)=0, \qquad \textbf{Necessary (second order): }\ \nabla^2f(x^\star)\succeq0\]

\[\textbf{Sufficient: }\ \nabla f(x^\star)=0\ \textbf{ and }\ \nabla^2f(x^\star)\succ0\]

The gap between $\succeq$ and $\succ$ is not pedantry. $f(x)=x^4$ has a global minimum at $0$ where $f''(0)=0$, positive semi-definite, not positive definite. Requiring strict definiteness would reject it. Conversely, $f(x)=x^3$ has $f'(0)=f''(0)=0$ and no extremum at all, which is why the semi-definite condition is necessary but not sufficient.

cases <- data.frame(x = seq(-1.5, 1.5, length.out = 400)) |>
  mutate(`f(x) = x^2   (min, f'' > 0)`  = x^2,
         `f(x) = x^4   (min, f'' = 0)`  = x^4,
         `f(x) = x^3   (no extremum)`   = x^3) |>
  pivot_longer(-x, names_to = "case", values_to = "y")

ggplot(cases, aes(x, y)) +
  geom_hline(yintercept = 0, color = "grey75") +
  geom_vline(xintercept = 0, color = "grey75") +
  geom_line(linewidth = 1, color = "steelblue") +
  geom_point(data = data.frame(x = 0, y = 0, case = unique(cases$case)),
             color = "firebrick", size = 3) +
  facet_wrap(~ case, scales = "free_y") +
  labs(title = "All three have a zero gradient at the origin",
       subtitle = "Only the second-order behavior distinguishes them, and one of the two minima has a zero second derivative",
       x = NULL, y = NULL) +
  theme_dspa(10)

For multivariate problems the Hessian’s eigenvalues classify the critical point:

$\nabla^2f(x^\star)$	Eigenvalues	Critical point
Positive definite	all $>0$	strict local minimum
Negative definite	all $<0$	strict local maximum
Indefinite	mixed signs	saddle point
Semi-definite, singular	some $=0$	inconclusive at second order

Saddle points, not local minima, are the obstacle in high dimensions. For a random symmetric matrix in $n$ dimensions, the probability that all $n$ eigenvalues share a sign decays exponentially in $n$, so almost every critical point of a generic high-dimensional function is a saddle (Dauphin et al., 2014). Gradient descent escapes strict saddles almost surely; the practical difficulty is that it can slow down enormously near one, which is what plateaus in deep-learning training curves usually are.

5 Convexity

A set $C$ is convex if $tx+(1-t)y\in C$ for all $x,y\in C$, $t\in[0,1]$. A function $f:C\to\mathbb R$ on a convex $C$ is convex if

\[f\big(tx+(1-t)y\big)\ \le\ t f(x)+(1-t)f(y),\qquad \forall x,y\in C,\ t\in[0,1],\]

and strictly convex if the inequality is strict for $x\ne y$, $t\in(0,1)$.

For differentiable $f$ there are two equivalent characterizations, and both are worth having:

\[ \begin{aligned} \textbf{First order: }&\quad f(y)\ \ge\ f(x)+\nabla f(x)^\top(y-x)\quad\forall x,y\\ \textbf{Second order: }&\quad \nabla^2f(x)\succeq0\quad\forall x \end{aligned} \]

The first-order form says the function lies above all its tangent planes — which is exactly why a zero gradient certifies a global minimum.

Common misconception: “convex functions have no local minima.” They have local minima; the theorem is that every local minimum is a global minimum. The proof is one line from the definition: if $x^\star$ were local but not global, some $y$ has $f(y)<f(x^\star)$, and points on the segment $tx^\star+(1-t)y$ arbitrarily close to $x^\star$ satisfy $f\le tf(x^\star)+(1-t)f(y)<f(x^\star)$, contradicting local optimality.

Two further clarifications. Convexity does not imply a unique minimizer — $f(x)=\max(0,|x|-1)$ is convex with a whole interval of minimizers; strict convexity gives uniqueness. And convex functions can attain local maxima, but only on the boundary of the feasible set.

Common misconception: “convex problems are easy and non-convex ones are NP-hard.” Both halves need qualification.

Convex problems are solvable to accuracy $\epsilon$ in time polynomial in the problem dimension and $\log(1/\epsilon)$ under a first-order oracle model with a self-concordant barrier, not unconditionally, and not for every representation (a convex problem with exponentially many constraints is still hard to write down).

“Non-convex optimization” is a class, containing many tractable instances — PCA, matrix factorization under conditions, and phase retrieval all have non-convex formulations with tractable global solutions. The sharp results are that finding a global minimum of a general non-convex function is NP-hard, and that even deciding whether a given point is a local minimum of a degree-4 polynomial is NP-hard (Murty & Kabadi, 1987).

xs <- seq(-2, 2, length.out = 400)
conv_df <- bind_rows(
  data.frame(x = xs, y = xs^2,                     kind = "Strictly convex"),
  data.frame(x = xs, y = pmax(0, abs(xs) - 0.8),   kind = "Convex, not strict"),
  data.frame(x = xs, y = xs^4 - 3*xs^2 + 0.5*xs,   kind = "Non-convex"))

ggplot(conv_df, aes(x, y)) +
  geom_line(linewidth = 1, color = "steelblue") +
  facet_wrap(~ kind, scales = "free_y") +
  labs(title = "Three landscapes",
       subtitle = "Left: unique minimizer. Middle: an interval of minimizers. Right: two local minima, one global",
       x = NULL, y = NULL) +
  theme_dspa(10)

# --- Interactive equivalent ------------------------------------------------
plot_ly(x = xs, y = xs^2, type = "scatter", mode = "lines",
        name = "Strictly convex") |>
  add_lines(y = pmax(0, abs(xs) - 0.8), name = "Convex, not strict") |>
  add_lines(y = xs^4 - 3*xs^2 + 0.5*xs, name = "Non-convex") |>
  layout(title = "Convex and non-convex objectives",
         xaxis = list(title = "x"), yaxis = list(title = "f(x)"))

6 Smoothness, strong convexity, and the condition number

Two quantitative properties determine how fast first-order methods converge.

$L$-smoothness (upper curvature bound). $f$ is $L$-smooth if $\nabla f$ is $L$-Lipschitz, equivalently $\nabla^2f(x)\preceq LI$: \[f(y)\ \le\ f(x)+\nabla f(x)^\top(y-x)+\frac{L}{2}\|y-x\|^2\] , the function is bounded above by a quadratic.

$\mu$-strong convexity (lower curvature bound). $f$ is $\mu$-strongly convex if $\nabla^2f(x)\succeq\mu I$ with $\mu>0$: \[f(y)\ \ge\ f(x)+\nabla f(x)^\top(y-x)+\frac{\mu}{2}\|y-x\|^,\] the function is bounded below by a quadratic.

Together they sandwich $f$ between two parabolas, and their ratio

\[\boxed{\;\kappa=\frac{L}{\mu}\ \ge\ 1\;}\]

is the condition number. It is the single number that governs first-order convergence, and it is the same conditioning that governs linear-system accuracy in Chapter 3, §3.7, for a quadratic $f(x)=\tfrac12x^\top Ax$, $\kappa$ is exactly $\lambda_{\max}(A)/\lambda_{\min}(A)$.

quad_factory <- function(kappa) {
  A <- diag(c(1, kappa))
  list(A = A, kappa = kappa,
       f    = function(x) 0.5 * as.numeric(t(x) %*% A %*% x),
       grad = function(x) as.vector(A %*% x),
       L = max(diag(A)), mu = min(diag(A)))
}
data.frame(kappa = c(1, 5, 20, 100),
           L = c(1, 5, 20, 100), mu = 1,
           max_stable_step = round(2 / c(1, 5, 20, 100), 4),
           optimal_step = round(1 / c(1, 5, 20, 100), 4),
           gd_rate_per_iter = round(1 - 1/c(1, 5, 20, 100), 4))

The last column is the per-iteration contraction factor $1-\mu/L=1-1/\kappa$. At $\kappa=100$ each step reduces the error by only 1%, so roughly $\kappa\log(1/\epsilon)$ iterations are needed, ill-conditioning is the dominant cost of first-order optimization.

kappa_grid <- 10^seq(0, 2.3, length.out = 45)
iter_grid <- 0:60
Zrate <- outer(kappa_grid, iter_grid, function(k, i) (1 - 1/k)^i)

plot_ly(x = iter_grid, y = kappa_grid, z = Zrate, type = "surface",
        colorscale = "Viridis",
        colorbar = list(title = "Relative\nerror")) |>
  layout(title = "Gradient-descent error contraction (1 - 1/kappa)^k",
         scene = list(xaxis = list(title = "Iteration k"),
                      yaxis = list(title = "Condition number kappa", type = "log"),
                      zaxis = list(title = "Error factor")))

Rotate along the iteration axis. At $\kappa=1$ the surface drops to zero immediately, a perfectly conditioned quadratic is solved in one gradient step. As $\kappa$ grows the descent flattens into a plateau: the same 60 iterations that solve the well-conditioned problem barely move the ill-conditioned one.

7 PART II, FIRST-ORDER METHODS

8 Gradient descent

One convention, used throughout. A vector $\nu$ is a descent direction at $x$ if $\nabla f(x)^\top\nu<0$, and the update is

\[\boxed{\;x_{k+1}=x_k+\alpha_k\nu_k\;}\]

with step size $\alpha_k>0$. Steepest descent takes $\nu_k=-\nabla f(x_k)$, giving the familiar

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

Mixing the sign into $\nu$ and then subtracting is the most common source of confusion in optimization code; keeping the minus inside $\nu$ and always adding the step removes it.

# # Descent direction convention: nu satisfies t(grad) %*% nu < 0, update x + a*nu.
# gradient_descent <- function(f, grad, x0, alpha = 0.1, max_iter = 500,
#                              tol = 1e-8) {
#   x <- x0
#   path <- matrix(NA_real_, length(x0), max_iter + 1); path[, 1] <- x
#   fvals <- numeric(max_iter + 1); fvals[1] <- f(x)
#   k <- 0
#   for (i in seq_len(max_iter)) {
#     g <- grad(x)
#     if (sqrt(sum(g^2)) < tol) break
#     nu <- -g                                   # descent direction
#     stopifnot(sum(g * nu) < 0 || sqrt(sum(g^2)) < tol)
#     x <- x + alpha * nu                        # always ADD the step
#     k <- i; path[, i + 1] <- x; fvals[i + 1] <- f(x)
#   }
#   list(x = x, f = f(x), iterations = k,
#        path = path[, 1:(k + 1), drop = FALSE], fvals = fvals[1:(k + 1)])
# }
gradient_descent <- function(f, grad, x0, alpha = 0.1, max_iter = 500,
                             tol = 1e-8, max_backtracks = 30) {
  x <- x0
  path <- matrix(NA_real_, length(x0), max_iter + 1)
  path[, 1] <- x
  fvals <- numeric(max_iter + 1)
  fvals[1] <- f(x)
  k <- 0
  
  for (i in seq_len(max_iter)) {
    # Compute gradient safely
    g <- tryCatch(grad(x), error = function(e) rep(NA_real_, length(x)))
    if (any(!is.finite(g))) {
      warning("Non-finite gradient encountered at iteration ", i, ". Stopping.")
      break
    }
    if (sqrt(sum(g^2)) < tol) break
    
    nu <- -g                         # descent direction
    # The assertion should always hold if g is finite; but keep a safety check
    if (!(sum(g * nu) < 0)) {
      warning("Gradient is not a descent direction. Stopping.")
      break
    }
    
    # Backtracking line search
    step <- alpha
    x_new <- x + step * nu
    f_new <- tryCatch(f(x_new), error = function(e) NA_real_)
    bt <- 0
    while ((!is.finite(f_new) || f_new > fvals[i] + 1e-12 * sum(g * nu)) 
           && bt < max_backtracks) {
      step <- step / 2
      x_new <- x + step * nu
      f_new <- tryCatch(f(x_new), error = function(e) NA_real_)
      bt <- bt + 1
    }
    
    if (!is.finite(f_new) || bt >= max_backtracks) {
      warning("No finite descent step found. Stopping.")
      break
    }
    
    x <- x_new
    k <- i
    path[, i + 1] <- x
    fvals[i + 1] <- f_new
  }
  
  list(x = x, f = f(x), iterations = k,
       path = path[, 1:(k + 1), drop = FALSE],
       fvals = fvals[1:(k + 1)])
}

8.1 Stability: how large may the step be?

Common misconception: “a smaller learning rate is always safer, and the only cost is speed.” There is a hard threshold, and it is computable.

For an $L$-smooth $f$, the smoothness inequality gives \[f(x-\alpha\nabla f)\ \le\ f(x)-\alpha\Big(1-\frac{\alpha L}{2}\Big)\|\nabla f(x)\|^2,\] so the objective is guaranteed to decrease exactly when $\alpha\big(1-\tfrac{\alpha L}{2}\big)>0$, i.e. \[\boxed{\;0<\alpha<\frac{2}{L}\;}\]

At $\alpha>2/L$ gradient descent diverges, not slowly, but geometrically. On a quadratic with Hessian $A$, the error evolves as $(I-\alpha A)^k e_0$, whose spectral radius is $\max_i|1-\alpha\lambda_i|$; that exceeds 1 as soon as $\alpha>2/\lambda_{\max}=2/L$.

The optimal constant step for a quadratic is $\alpha^\star=\frac{2}{L+\mu}$, and $\alpha=1/L$ is the standard safe choice. The “zig-zagging” that gradient descent is often accused of is not a defect of the method, it is the regime $\frac{1}{L}<\alpha<\frac{2}{L}$ on an ill-conditioned problem, where the largest-curvature direction oscillates while the smallest-curvature direction crawls.

q <- quad_factory(kappa = 10)
L <- q$L; mu <- q$mu

steps <- c(0.05, 1/L, 2/(L + mu), 0.19, 0.21)
labels <- c("alpha = 0.05 (small)", "alpha = 1/L (safe)",
            "alpha = 2/(L+mu) (optimal)", "alpha = 0.19 (< 2/L)",
            "alpha = 0.21 (> 2/L: diverges)")

conv <- bind_rows(lapply(seq_along(steps), function(i) {
  r <- gradient_descent(q$f, q$grad, c(5, 5), alpha = steps[i],
                        max_iter = 60, tol = 1e-14)
  data.frame(k = 0:r$iterations, f = pmax(r$fvals, 1e-16), setting = labels[i])
}))
c(L = L, two_over_L = round(2/L, 4), optimal = round(2/(L + mu), 4))

#>          L two_over_L    optimal 
#>    10.0000     0.2000     0.1818

ggplot(conv, aes(k, f, color = setting)) +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  scale_color_viridis_d(option = "turbo", end = 0.9) +
  labs(title = "Gradient descent across the stability threshold",
       subtitle = sprintf("With L = %g, the boundary is at alpha = 2/L = %.2f. Beyond it the iterates blow up", L, 2/L),
       x = "Iteration", y = "f(x_k), log scale", color = NULL) +
  theme_dspa(10)

alpha_grid <- seq(0.01, 0.30, length.out = 60)
L_grid <- seq(2, 30, length.out = 60)
# Spectral radius of (I - alpha*A) for a quadratic with eigenvalues {1, L}
Zspec <- outer(L_grid, alpha_grid, function(Lv, a)
  pmin(pmax(abs(1 - a * Lv), abs(1 - a * 1)), 3))

plot_ly(x = alpha_grid, y = L_grid, z = Zspec, type = "surface",
        colorscale = "RdBu", reversescale = TRUE,
        colorbar = list(title = "Spectral\nradius")) |>
  add_surface(x = alpha_grid, y = L_grid,
              z = matrix(1, length(L_grid), length(alpha_grid)),
              opacity = 0.35, showscale = FALSE,
              colorscale = list(c(0, "black"), c(1, "black")),
              name = "Stability boundary") |>
  layout(title = "Convergence factor over step size and smoothness; the flat plane is radius 1",
         scene = list(xaxis = list(title = "Step size alpha"),
                      yaxis = list(title = "Smoothness L"),
                      zaxis = list(title = "Spectral radius", range = c(0, 2))))

The flat plane sits at spectral radius 1. Where the colored surface is below it the method converges; where it rises through it, the iterates diverge. The crossing curve is exactly $\alpha=2/L$, and the valley floor traces the optimal step for each $L$.

8.2 Convergence rates

Rates for gradient descent with $\alpha=1/L$.

\[ \begin{aligned} \textbf{$L$-smooth, convex: }&\quad f(x_k)-f^\star\ \le\ \frac{L\|x_0-x^\star\|^2}{2k} &&= O(1/k),\ \text{so } O(1/\epsilon)\ \text{iterations}\\[2mm] \textbf{$+\ \mu$-strongly convex: }&\quad \|x_k-x^\star\|^2\ \le\ \Big(1-\tfrac{\mu}{L}\Big)^k\|x_0-x^\star\|^2 &&\text{\emph{linear} rate, } O\big(\kappa\log\tfrac1\epsilon\big)\ \text{iterations} \end{aligned} \]

The difference is qualitative. Without strong convexity the error decays like $1/k$, reaching $10^{-6}$ needs about $10^6$ iterations. With it, the decay is geometric and $10^{-6}$ needs about $14\kappa$.

set.seed(21)
rate_check <- function(kappa, alpha = NULL, iters = 200) {
  q <- quad_factory(kappa)
  a <- if (is.null(alpha)) 1 / q$L else alpha
  r <- gradient_descent(q$f, q$grad, c(1, 1), alpha = a,
                        max_iter = iters, tol = 0)
  err <- colSums(r$path^2)
  data.frame(k = 0:(length(err) - 1), observed = err / err[1],
             theory = (1 - 1/kappa)^(2 * (0:(length(err) - 1))),
             kappa = sprintf("kappa = %d", kappa))
}
rr <- bind_rows(lapply(c(2, 10, 50), rate_check))

ggplot(rr, aes(k)) +
  geom_line(aes(y = pmax(observed, 1e-16), color = "Observed"), linewidth = 0.9) +
  geom_line(aes(y = pmax(theory, 1e-16), color = "Theory (1 - 1/kappa)^2k"),
            linetype = "dashed", linewidth = 0.8) +
  facet_wrap(~ kappa) +
  scale_y_log10() +
  scale_color_manual(values = c(Observed = "#3B7DD8",
                                 `Theory (1 - 1/kappa)^2k` = "#D8433B")) +
  labs(title = "Observed contraction against the theoretical linear rate",
       subtitle = "Squared distance to the optimum, normalized. The bound is tight for the quadratic",
       x = "Iteration", y = "Relative squared error", color = NULL) +
  theme_dspa(10)

9 Line search

Fixing $\alpha$ requires knowing $L$, which is rarely available. Backtracking line search finds an acceptable step adaptively, using only function and gradient evaluations.

Armijo (sufficient decrease) condition. Accept $\alpha$ when \[f(x+\alpha\nu)\ \le\ f(x)+c_1\alpha\,\nabla f(x)^\top\nu,\qquad c_1\in(0,1),\] typically $c_1=10^{-4}$. The right side is the linear prediction discounted by $c_1$, the step must deliver at least a fraction $c_1$ of the decrease the tangent line promises.

Curvature (Wolfe) condition. Additionally require \[\nabla f(x+\alpha\nu)^\top\nu\ \ge\ c_2\,\nabla f(x)^\top\nu,\qquad c_1<c_2<1,\] typically $c_2=0.9$. This rules out steps that are too short: it demands that the directional derivative has flattened appreciably.

Armijo alone suffices for gradient descent. The curvature condition is what quasi-Newton methods require, because it guarantees $y_k^\top s_k>0$, the condition under which BFGS preserves positive definiteness (§13.13).

backtracking <- function(f, grad, x, nu, alpha0 = 1, c1 = 1e-4, rho = 0.5,
                         max_bt = 50) {
  fx <- f(x); gx <- grad(x); slope <- sum(gx * nu)
  stopifnot(slope < 0)                          # nu must be a descent direction
  a <- alpha0
  for (i in seq_len(max_bt)) {
    if (f(x + a * nu) <= fx + c1 * a * slope) return(list(alpha = a, backtracks = i - 1))
    a <- rho * a
  }
  list(alpha = a, backtracks = max_bt)
}

gd_linesearch <- function(f, grad, x0, max_iter = 200, tol = 1e-8) {
  x <- x0; path <- matrix(NA_real_, length(x0), max_iter + 1); path[, 1] <- x
  alphas <- numeric(max_iter); k <- 0
  for (i in seq_len(max_iter)) {
    g <- grad(x); if (sqrt(sum(g^2)) < tol) break
    nu <- -g
    ls <- backtracking(f, grad, x, nu)
    x <- x + ls$alpha * nu
    alphas[i] <- ls$alpha; k <- i; path[, i + 1] <- x
  }
  list(x = x, f = f(x), iterations = k, alphas = alphas[1:max(k, 1)],
       path = path[, 1:(k + 1), drop = FALSE])
}

# Rosenbrock: the standard ill-conditioned non-quadratic test
rosen  <- function(x) 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2
rosen_g <- function(x) c(-400 * x[1] * (x[2] - x[1]^2) - 2 * (1 - x[1]),
                          200 * (x[2] - x[1]^2))

fixed <- gradient_descent(rosen, rosen_g, c(-1.2, 1), alpha = 0.001,
                          max_iter = 2000, tol = 1e-8)
ls_res <- gd_linesearch(rosen, rosen_g, c(-1.2, 1), max_iter = 2000)

data.frame(
  method = c("fixed step alpha = 0.001", "backtracking line search"),
  iterations = c(fixed$iterations, ls_res$iterations),
  final_f = signif(c(fixed$f, ls_res$f), 4),
  distance_to_optimum = signif(c(sqrt(sum((fixed$x - c(1,1))^2)),
                                 sqrt(sum((ls_res$x - c(1,1))^2))), 4))

c(step_range_used = signif(range(ls_res$alphas), 3),
  median_step = signif(median(ls_res$alphas), 3))

#> step_range_used1 step_range_used2      median_step 
#>         0.000977         0.500000         0.001950

The line search adapts the step over three orders of magnitude across the run — large in the flat valley, small at the curved ridge. No single fixed $\alpha$ can do both.

10 Momentum and acceleration

Gradient descent’s $O(1/k)$ rate is not the best achievable by a first-order method.

\[ \begin{aligned} \textbf{Heavy-ball (Polyak): }&\quad x_{k+1}=x_k-\alpha\nabla f(x_k)+\beta(x_k-x_{k-1})\\[2mm] \textbf{Nesterov: }&\quad \begin{cases} y_k = x_k+\beta_k(x_k-x_{k-1})\\ x_{k+1}=y_k-\alpha\nabla f(y_k) \end{cases} \end{aligned} \]

The difference is where the gradient is evaluated: heavy-ball uses $\nabla f(x_k)$, Nesterov uses $\nabla f(y_k)$ at the extrapolated point, a “look-ahead” that lets the method correct an overshoot within the same step.

Nesterov acceleration achieves $O(1/k^2)$, and that is optimal. \[f(x_k)-f^\star\ \le\ \frac{2L\|x_0-x^\star\|^2}{(k+1)^2}\] and under strong convexity the rate improves from $\big(1-\frac1\kappa\big)^k$ to $\big(1-\frac{1}{\sqrt\kappa}\big)^k$, the square root of the condition number, which at $\kappa=10^4$ is a hundredfold reduction in iterations.

The Nemirovski–Yudin lower bound shows no method using only gradients can do better than $\Omega(1/k^2)$ on the smooth convex class, so acceleration is not merely a good heuristic: it closes the gap.

nesterov <- function(f, grad, x0, alpha, max_iter = 500, tol = 1e-10) {
  x <- x0; x_prev <- x0
  path <- matrix(NA_real_, length(x0), max_iter + 1); path[, 1] <- x
  fvals <- numeric(max_iter + 1); fvals[1] <- f(x); k <- 0
  for (i in seq_len(max_iter)) {
    beta <- (i - 1) / (i + 2)                  # the standard schedule
    y <- x + beta * (x - x_prev)
    g <- grad(y); if (sqrt(sum(g^2)) < tol) break
    x_prev <- x; x <- y - alpha * g
    k <- i; path[, i + 1] <- x; fvals[i + 1] <- f(x)
  }
  list(x = x, f = f(x), iterations = k,
       path = path[, 1:(k+1), drop = FALSE], fvals = fvals[1:(k+1)])
}

heavy_ball <- function(f, grad, x0, alpha, beta = 0.9, max_iter = 500, tol = 1e-10) {
  x <- x0; x_prev <- x0
  fvals <- numeric(max_iter + 1); fvals[1] <- f(x); k <- 0
  for (i in seq_len(max_iter)) {
    g <- grad(x); if (sqrt(sum(g^2)) < tol) break
    x_new <- x - alpha * g + beta * (x - x_prev)
    x_prev <- x; x <- x_new
    k <- i; fvals[i + 1] <- f(x)
  }
  list(x = x, f = f(x), iterations = k, fvals = fvals[1:(k+1)])
}

qk <- quad_factory(kappa = 100)
a_opt <- 1 / qk$L
x_start <- c(3, 3)

gd_r <- gradient_descent(qk$f, qk$grad, x_start, alpha = a_opt, max_iter = 300, tol = 0)
nv_r <- nesterov(qk$f, qk$grad, x_start, alpha = a_opt, max_iter = 300, tol = 0)
hb_r <- heavy_ball(qk$f, qk$grad, x_start, alpha = a_opt, beta = 0.85,
                   max_iter = 300, tol = 0)

pad <- function(v, n) c(v, rep(tail(v, 1), max(0, n - length(v))))[1:n]
n_show <- 300
mom_df <- data.frame(
  k = rep(1:n_show, 3),
  f = pmax(c(pad(gd_r$fvals, n_show), pad(nv_r$fvals, n_show), pad(hb_r$fvals, n_show)), 1e-16),
  method = rep(c("Gradient descent", "Nesterov", "Heavy-ball"), each = n_show))

ggplot(mom_df, aes(k, f, color = method)) +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  scale_color_manual(values = c(`Gradient descent` = "#D8433B",
                                 Nesterov = "#3B7DD8", `Heavy-ball` = "#7FB069")) +
  labs(title = sprintf("Acceleration on a quadratic with kappa = %d", qk$kappa),
       subtitle = "Same step size for all three. The momentum methods reach in tens of iterations what GD needs hundreds for",
       x = "Iteration", y = "f(x_k), log scale", color = NULL) +
  theme_dspa()

# --- Interactive equivalent ------------------------------------------------
plot_ly(x = 1:n_show, y = pad(gd_r$fvals, n_show), type = "scatter",
        mode = "lines", name = "Gradient descent") |>
  add_lines(y = pad(nv_r$fvals, n_show), name = "Nesterov") |>
  add_lines(y = pad(hb_r$fvals, n_show), name = "Heavy-ball") |>
  layout(title = "Momentum methods on an ill-conditioned quadratic",
         xaxis = list(title = "Iteration"),
         yaxis = list(title = "f(x)", type = "log"))

11 Adaptive methods

Deep learning uses a different family, which rescales each coordinate by a running statistic of its own gradient history.

\[ \begin{aligned} \textbf{AdaGrad: }&\quad G_k=G_{k-1}+g_k^{\odot2},\qquad x_{k+1}=x_k-\frac{\alpha}{\sqrt{G_k}+\epsilon}\odot g_k\\[2mm] \textbf{RMSProp: }&\quad v_k=\rho v_{k-1}+(1-\rho)g_k^{\odot2},\qquad x_{k+1}=x_k-\frac{\alpha}{\sqrt{v_k}+\epsilon}\odot g_k\\[2mm] \textbf{Adam: }&\quad \begin{cases} m_k=\beta_1m_{k-1}+(1-\beta_1)g_k,\quad \hat m_k=m_k/(1-\beta_1^k)\\ v_k=\beta_2v_{k-1}+(1-\beta_2)g_k^{\odot2},\quad \hat v_k=v_k/(1-\beta_2^k)\\ x_{k+1}=x_k-\alpha\,\hat m_k/(\sqrt{\hat v_k}+\epsilon) \end{cases} \end{aligned} \]

AdaGrad’s accumulator never decays, so $\sqrt{G_k}$ grows without bound and the effective step size decays to zero, excellent for sparse gradients, fatal for long training runs. RMSProp replaces the sum with an exponential moving average, which fixes that. Adam adds momentum on the first moment and bias correction: at $k=1$, $m_1=(1-\beta_1)g_1$ is badly biased toward zero, and dividing by $1-\beta_1^k$ removes it.

Adam is a diagonal preconditioner, it approximates $\nabla^2f$ by a diagonal matrix of gradient second moments, which is why it helps most on problems whose coordinates have wildly different curvature, and why it can hurt on well-conditioned ones.

adaptive_opt <- function(f, grad, x0, method = c("adam", "rmsprop", "adagrad", "gd"),
                         alpha = 0.1, beta1 = 0.9, beta2 = 0.999, rho = 0.9,
                         eps = 1e-8, max_iter = 500) {
  method <- match.arg(method)
  x <- x0; m <- v <- G <- rep(0, length(x0))
  fvals <- numeric(max_iter + 1); fvals[1] <- f(x)
  for (k in seq_len(max_iter)) {
    g <- grad(x)
    x <- switch(method,
      gd      = x - alpha * g,
      adagrad = { G <<- G + g^2; x - alpha * g / (sqrt(G) + eps) },
      rmsprop = { v <<- rho * v + (1 - rho) * g^2; x - alpha * g / (sqrt(v) + eps) },
      adam    = {
        m <<- beta1 * m + (1 - beta1) * g
        v <<- beta2 * v + (1 - beta2) * g^2
        mh <- m / (1 - beta1^k); vh <- v / (1 - beta2^k)
        x - alpha * mh / (sqrt(vh) + eps)
      })
    fvals[k + 1] <- f(x)
  }
  list(x = x, f = f(x), fvals = fvals)
}

set.seed(31)
# A badly scaled objective: coordinate curvatures span four orders of magnitude
A_scaled <- diag(c(1e-2, 1, 20, 1e2))
f_sc  <- function(x) 0.5 * sum(diag(A_scaled) * x^2)
g_sc  <- function(x) diag(A_scaled) * x
x0_sc <- rep(2, 4)

ad_df <- bind_rows(lapply(c("gd", "adagrad", "rmsprop", "adam"), function(m) {
  r <- adaptive_opt(f_sc, g_sc, x0_sc, method = m, alpha = 0.02, max_iter = 400)
  data.frame(k = 0:400, f = pmax(r$fvals, 1e-16), method = toupper(m))
}))

ggplot(ad_df, aes(k, f, color = method)) +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  scale_color_viridis_d(option = "plasma", end = 0.9) +
  labs(title = "Adaptive methods on a badly scaled objective",
       subtitle = sprintf("Coordinate curvatures span %g to %g -- condition number %g",
                          min(diag(A_scaled)), max(diag(A_scaled)),
                          max(diag(A_scaled))/min(diag(A_scaled))),
       x = "Iteration", y = "f(x), log scale", color = NULL) +
  theme_dspa()

12 Stochastic gradient descent

For an objective that is a sum over data,

\[f(\beta)=\frac1n\sum_{i=1}^{n}\ell_i(\beta),\qquad \nabla f(\beta)=\frac1n\sum_{i=1}^{n}\nabla\ell_i(\beta),\]

each full gradient costs $O(n)$ component evaluations. Stochastic gradient descent replaces it with a mini-batch estimate over $|B|=b$ indices:

\[g_B(\beta)=\frac1b\sum_{i\in B}\nabla\ell_i(\beta), \qquad \mathbb E[g_B]=\nabla f(\beta)\ \ \text{(unbiased)}, \qquad \operatorname{Var}(g_B)=\frac{\sigma^2}{b}.\]

The variance is what forces the step size to decay. With a constant step $\alpha$, SGD does not converge to $x^\star$; it converges to a noise ball of radius $O(\alpha\sigma^2/\mu)$ around it and then bounces. Convergence requires the Robbins–Monro conditions \[\sum_k\alpha_k=\infty\quad\text{(reach the optimum)},\qquad \sum_k\alpha_k^2<\infty\quad\text{(damp the noise)},\] satisfied by $\alpha_k=\alpha_0/k$ but not by $\alpha_k=\alpha_0/\sqrt k$ in the second condition, which is why the $1/\sqrt k$ schedule gives $O(1/\sqrt k)$ convergence in function value rather than a.s. convergence of the iterates.

Rates: $O(1/\sqrt k)$ for convex, $O(1/k)$ for strongly convex with $\alpha_k\propto1/k$. Both are worse per iteration than full gradient descent, but each iteration costs $b/n$ as much, so SGD wins decisively when $n$ is large and moderate accuracy suffices.

set.seed(41)
n_sgd <- 20000; p_sgd <- 20
X_sgd <- matrix(rnorm(n_sgd * p_sgd), n_sgd, p_sgd)
beta_true <- c(rep(1.5, 5), rep(0, p_sgd - 5))
y_sgd <- as.vector(X_sgd %*% beta_true) + rnorm(n_sgd)
beta_ols <- as.vector(qr.solve(X_sgd, y_sgd))

f_ls <- function(b) mean((y_sgd - X_sgd %*% b)^2)
g_ls <- function(b) as.vector(-2 * crossprod(X_sgd, y_sgd - X_sgd %*% b) / n_sgd)

sgd <- function(b0, batch = 64, alpha0 = 0.05, epochs = 30,
                schedule = c("constant", "sqrt", "inverse")) {
  schedule <- match.arg(schedule)
  b <- b0; hist <- numeric(0); t <- 0
  for (e in seq_len(epochs)) {
    idx <- sample(n_sgd)
    for (s in seq(1, n_sgd - batch, by = batch)) {
      t <- t + 1
      B <- idx[s:(s + batch - 1)]
      g <- as.vector(-2 * crossprod(X_sgd[B, ], y_sgd[B] - X_sgd[B, ] %*% b) / batch)
      a <- switch(schedule, constant = alpha0,
                  sqrt = alpha0 / sqrt(t), inverse = alpha0 / t)
      b <- b - a * g
    }
    hist <- c(hist, sum((b - beta_ols)^2))
  }
  list(beta = b, hist = hist)
}

set.seed(43)
sgd_df <- bind_rows(lapply(c("constant", "sqrt", "inverse"), function(s) {
  r <- sgd(rep(0, p_sgd), schedule = s)
  data.frame(epoch = seq_along(r$hist), err = pmax(r$hist, 1e-16),
             schedule = sprintf("alpha_t = %s", switch(s,
               constant = "alpha_0", sqrt = "alpha_0 / sqrt(t)", inverse = "alpha_0 / t")))
}))
gd_full <- {
  b <- rep(0, p_sgd); h <- numeric(30)
  for (e in 1:30) { b <- b - 0.05 * g_ls(b); h[e] <- sum((b - beta_ols)^2) }
  data.frame(epoch = 1:30, err = pmax(h, 1e-16), schedule = "full-batch GD")
}

ggplot(bind_rows(sgd_df, gd_full), aes(epoch, err, color = schedule)) +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  scale_color_viridis_d(option = "turbo", end = 0.9) +
  labs(title = "SGD step-size schedules against full-batch gradient descent",
       subtitle = "A constant step stalls in a noise ball; decaying steps converge, at the cost of slowing down",
       x = "Epoch", y = "Squared distance to the OLS solution, log scale", color = NULL) +
  theme_dspa()

The constant-step curve flattens at a nonzero error, that plateau is the noise ball, and its height is proportional to $\alpha\sigma^2$. Both decaying schedules keep descending.

13 Cost per iteration

For $f(\beta)=\frac1n\sum\ell_i(\beta)$ with $\beta\in\mathbb R^d$:

Method	Cost per iteration	Memory	Rate (strongly convex)
Full gradient descent	$O(nd)$	$O(d)$	$(1-1/\kappa)^k$
Nesterov	$O(nd)$	$O(d)$	$(1-1/\sqrt\kappa)^k$
Adam / RMSProp	$O(nd)$	$O(d)$	problem-dependent
SGD, batch $b$	$O(bd)$	$O(d)$	$O(1/k)$ in expectation
Newton	$O(nd^2+d^3)$	$O(d^2)$	quadratic, locally
BFGS	$O(nd+d^2)$	$O(d^2)$	superlinear
L-BFGS, memory $m$	$O(nd+md)$	$O(md)$	superlinear

The $O(d^2)$ memory is what rules Newton and BFGS out of deep learning. At $d=10^8$ parameters, a dense Hessian would need $4\times10^{16}$ bytes. L-BFGS’s $O(md)$ with $m\approx10$ is affordable; Adam’s $O(d)$ is what is actually used.

14 PART III: SECOND-ORDER AND QUASI-NEWTON METHODS

15 Newton’s method

Gradient descent minimizes a linear model of $f$ plus a proximity penalty. Newton’s method minimizes the quadratic model directly.

The second-order Taylor expansion at $x_k$ is

\[f(x_k+s)\ \approx\ \underbrace{f(x_k)}_{\text{constant}}+\underbrace{\nabla f(x_k)^\top s}_{\text{linear}}+\underbrace{\tfrac12 s^\top\nabla^2f(x_k)\,s}_{\text{quadratic}}\ \equiv\ m_k(s).\]

Setting $\nabla m_k(s)=\nabla f(x_k)+\nabla^2f(x_k)s=0$ gives the Newton step

\[\boxed{\;s_k=-\big[\nabla^2f(x_k)\big]^{-1}\nabla f(x_k),\qquad x_{k+1}=x_k+\alpha_k s_k\;}\]

Note that $\nabla f$ is a column vector and $\nabla^2 f$ a symmetric matrix, so no transposes appear: $H^{-1}g$ is well formed as written.

Common misconception: “Newton’s method is just a better gradient descent — use it whenever you can afford the Hessian.” Its guarantee is local, and outside that neighbourhood it can be actively worse than gradient descent.

The theorem: if $f$ is twice differentiable with $\nabla^2f$ Lipschitz ($\|\nabla^2f(x)-\nabla^2f(y)\|\le M\|x-y\|$) and $\nabla^2f(x^\star)\succeq\mu I$, then for $x_k$ close enough to $x^\star$ and $\alpha_k=1$, \[\|x_{k+1}-x^\star\|\ \le\ \frac{M}{2\mu}\|x_k-x^\star\|^2.\] The number of correct digits doubles each iteration, so reaching $\epsilon$ from a good start takes $O(\log\log(1/\epsilon))$ steps, 5 or 6 in practice, regardless of $\epsilon$.

The word locally carries the weight. Far from the optimum the quadratic model can be a poor description, and if $\nabla^2f(x_k)$ is not positive definite the “Newton step” may point uphill or head straight for a saddle — and in high dimensions most critical points are saddles. Three standard repairs:

Damping: use a line search on $\alpha_k$ rather than $\alpha_k=1$.

Regularization: solve $(\nabla^2f+\tau I)s=-\nabla f$ with $\tau>0$ large enough to make the matrix positive definite, which interpolates between Newton ($\tau\to0$) and gradient descent ($\tau\to\infty$), and is exactly the Levenberg–Marquardt idea.

Trust region: minimize $m_k(s)$ subject to $\|s\|\le\Delta_k$, adapting $\Delta_k$ by how well the model predicted the actual decrease.

newton_method <- function(f, grad, hess, x0, max_iter = 100, tol = 1e-12,
                          damped = TRUE, tau_floor = 1e-8) {
  x <- x0; path <- matrix(NA_real_, length(x0), max_iter + 1); path[, 1] <- x
  errs <- numeric(max_iter + 1); k <- 0
  for (i in seq_len(max_iter)) {
    g <- grad(x); if (sqrt(sum(g^2)) < tol) break
    H <- hess(x)
    # Regularize until positive definite, so the step is a descent direction
    tau <- 0; ev <- min(eigen(H, symmetric = TRUE, only.values = TRUE)$values)
    if (ev <= tau_floor) tau <- tau_floor - ev + 1e-6
    s <- -solve(H + tau * diag(nrow(H)), g)
    if (sum(g * s) >= 0) s <- -g                # fall back to steepest descent
    alpha <- if (damped) backtracking(f, grad, x, s, alpha0 = 1)$alpha else 1
    x <- x + alpha * s
    k <- i; path[, i + 1] <- x
  }
  list(x = x, f = f(x), iterations = k, path = path[, 1:(k+1), drop = FALSE])
}

rosen_h <- function(x) matrix(c(
  -400 * (x[2] - x[1]^2) + 800 * x[1]^2 + 2, -400 * x[1],
  -400 * x[1],                                200), 2, 2)

nt <- newton_method(rosen, rosen_g, rosen_h, c(-1.2, 1))
gd_r2 <- gd_linesearch(rosen, rosen_g, c(-1.2, 1), max_iter = 5000)

data.frame(
  method = c("Newton (damped, regularized)", "Gradient descent + line search"),
  iterations = c(nt$iterations, gd_r2$iterations),
  final_f = signif(c(nt$f, gd_r2$f), 4),
  distance_to_optimum = signif(c(sqrt(sum((nt$x - c(1,1))^2)),
                                 sqrt(sum((gd_r2$x - c(1,1))^2))), 4))

errs_nt <- sqrt(colSums((nt$path - c(1, 1))^2))
errs_gd <- sqrt(colSums((gd_r2$path - c(1, 1))^2))

bind_rows(
  data.frame(k = seq_along(errs_nt) - 1, err = pmax(errs_nt, 1e-17), method = "Newton"),
  data.frame(k = seq_along(errs_gd) - 1, err = pmax(errs_gd, 1e-17),
             method = "Gradient descent")) |>
  filter(k <= 60) |>
  ggplot(aes(k, err, color = method)) +
  geom_line(linewidth = 0.9) + geom_point(size = 1.2) +
  scale_y_log10() +
  scale_color_manual(values = c(Newton = "#3B7DD8", `Gradient descent` = "#D8433B")) +
  labs(title = "Quadratic versus linear convergence on the Rosenbrock function",
       subtitle = "On a log scale, linear convergence is a straight line; quadratic convergence bends downward and falls off a cliff",
       x = "Iteration", y = "Distance to optimum, log scale", color = NULL) +
  theme_dspa()

# The signature of quadratic convergence: e_{k+1} / e_k^2 is roughly constant
tail_e <- errs_nt[errs_nt > 1e-14]
data.frame(k = head(seq_along(tail_e), -1),
           e_k = signif(head(tail_e, -1), 3),
           e_next = signif(tail(tail_e, -1), 3),
           ratio_linear = signif(tail(tail_e, -1) / head(tail_e, -1), 3),
           ratio_quadratic = signif(tail(tail_e, -1) / head(tail_e, -1)^2, 3)) |>
  tail(5)

The ratio_quadratic column stabilizes near a constant in the final iterations, while ratio_linear collapses toward zero, that is the numerical signature of $e_{k+1}\approx Ce_k^2$.

16 BFGS

Newton’s method needs $\nabla^2f$, which costs $O(nd^2)$ to form and $O(d^3)$ to factor. Quasi-Newton methods build an approximation from gradient differences alone.

Define the step and gradient change

\[s_k=x_{k+1}-x_k,\qquad y_k=\nabla f(x_{k+1})-\nabla f(x_k).\]

The secant condition requires the new approximation $B_{k+1}\approx\nabla^2f$ to reproduce the observed curvature:

\[B_{k+1}s_k=y_k\qquad\Longleftrightarrow\qquad H_{k+1}y_k=s_k,\quad H=B^{-1}.\]

The BFGS inverse-Hessian update. \[\boxed{\;H_{k+1}=\Big(I-\rho_k s_ky_k^\top\Big)\,H_k\,\Big(I-\rho_k y_ks_k^\top\Big)+\rho_k\,s_ks_k^\top, \qquad \rho_k=\frac{1}{y_k^\top s_k}\;}\]

Common misconception: “any update satisfying the secant condition will do.” The secant condition $H_{k+1}y_k=s_k$ has infinitely many solutions, and most of them are useless. Three structural features of the BFGS formula are doing specific work, and each is easy to get wrong.

The right factor is the transpose of the left one, the update has the congruence form $VHV^\top$ precisely so that symmetry is preserved. Copy the left factor instead and $H$ stops being symmetric within a few iterations, after which $-Hg$ is no longer guaranteed to be a descent direction and the method silently stops working. The rank-one correction is $s_ks_k^\top$, an outer product of $s$ with itself, not $s_ky_k^\top$. And the whole thing is scaled by $\rho_k=1/(y_k^\top s_k)$.

The curvature condition $y_k^\top s_k>0$ is required, and it is exactly what guarantees $H_{k+1}\succ0$ given $H_k\succ0$. For a strictly convex $f$ it holds automatically; in general it is enforced by the Wolfe curvature condition in the line search (§13.6), which is why BFGS is always paired with a Wolfe line search and not an Armijo-only one. When it fails, the standard remedy is to skip the update and carry $H_k$ forward.

bfgs <- function(f, grad, x0, max_iter = 300, tol = 1e-10, c2 = 0.9) {
  d <- length(x0)
  x <- x0; H <- diag(d)                      # start from the identity
  g <- grad(x); k <- 0; skipped <- 0
  fvals <- numeric(max_iter + 1); fvals[1] <- f(x)
  for (i in seq_len(max_iter)) {
    if (sqrt(sum(g^2)) < tol) break
    nu <- as.vector(-H %*% g)                # quasi-Newton descent direction
    if (sum(g * nu) >= 0) { H <- diag(d); nu <- -g }   # reset if not descent
    # Wolfe line search: Armijo + curvature, the latter to secure y'S > 0
    alpha <- 1; fx <- f(x); slope <- sum(g * nu)
    for (bt in 1:60) {
      x_new <- x + alpha * nu; g_new <- grad(x_new)
      armijo <- f(x_new) <= fx + 1e-4 * alpha * slope
      curv   <- sum(g_new * nu) >= c2 * slope
      if (armijo && curv) break
      alpha <- 0.5 * alpha
    }
    s <- x_new - x; y <- g_new - g
    sy <- sum(s * y)
    if (sy > 1e-10) {                        # curvature condition
      rho <- 1 / sy
      V <- diag(d) - rho * (s %*% t(y))
      H <- V %*% H %*% t(V) + rho * (s %*% t(s))   # note t(V), and s %*% t(s)
    } else skipped <- skipped + 1
    x <- x_new; g <- g_new; k <- i; fvals[i + 1] <- f(x)
  }
  list(x = x, f = f(x), iterations = k, H = H, skipped_updates = skipped,
       fvals = fvals[1:(k+1)])
}

bf <- bfgs(rosen, rosen_g, c(-1.2, 1))
ref <- optim(c(-1.2, 1), rosen, rosen_g, method = "BFGS",
             control = list(maxit = 1000, reltol = 1e-14))

data.frame(
  implementation = c("bfgs() above", "optim(method = 'BFGS')", "true optimum"),
  x1 = signif(c(bf$x[1], ref$par[1], 1), 8),
  x2 = signif(c(bf$x[2], ref$par[2], 1), 8),
  f  = signif(c(bf$f, ref$value, 0), 4))

c(iterations = bf$iterations, curvature_condition_failures = bf$skipped_updates,
  max_abs_difference_vs_optim = signif(max(abs(bf$x - ref$par)), 3))

#>                   iterations curvature_condition_failures 
#>                     3.60e+01                     3.00e+00 
#>  max_abs_difference_vs_optim 
#>                     1.51e-12

# The update MUST preserve symmetry and positive definiteness
c(H_is_symmetric = isTRUE(all.equal(bf$H, t(bf$H), tolerance = 1e-10)),
  H_min_eigenvalue = signif(min(eigen(bf$H, symmetric = TRUE,
                                      only.values = TRUE)$values), 4),
  H_is_positive_definite = min(eigen(bf$H, symmetric = TRUE,
                                     only.values = TRUE)$values) > 0)

#>         H_is_symmetric       H_min_eigenvalue H_is_positive_definite 
#>              1.0000000              0.0009965              1.0000000

# H should approximate the INVERSE Hessian at the optimum
data.frame(
  quantity = c("H (BFGS approximation)", "inverse true Hessian"),
  entry_11 = signif(c(bf$H[1,1], solve(rosen_h(c(1,1)))[1,1]), 4),
  entry_12 = signif(c(bf$H[1,2], solve(rosen_h(c(1,1)))[1,2]), 4),
  entry_22 = signif(c(bf$H[2,2], solve(rosen_h(c(1,1)))[2,2]), 4))

The symmetry and positive-definiteness checks are not decoration, an update that violates either produces a direction $\nu=-Hg$ that is not a descent direction, and the method silently stops working.

16.1 L-BFGS

BFGS stores a dense $d\times d$ matrix: $O(d^2)$ memory and $O(d^2)$ work per iteration. Limited-memory BFGS keeps only the last $m$ pairs $(s_j,y_j)$, typically $m=5$–$20$, and reconstructs the action $H_kg$ by a two-loop recursion in $O(md)$, never forming $H$ at all.

Method	Memory	Per-iteration	Feasible at $d=10^6$?
Newton	$O(d^2)$	$O(nd^2+d^3)$	No
BFGS	$O(d^2)$	$O(nd+d^2)$	No ($8\times10^{12}$ bytes)
L-BFGS ($m=10$)	$O(md)$	$O(nd+md)$	Yes ($8\times10^7$ bytes)
Gradient descent	$O(d)$	$O(nd)$	Yes

set.seed(51)
d_big <- 500
A_big <- crossprod(matrix(rnorm(d_big * (d_big + 20)), d_big + 20, d_big)) / d_big
b_big <- rnorm(d_big)
f_big <- function(x) 0.5 * as.numeric(t(x) %*% A_big %*% x) - sum(b_big * x)
g_big <- function(x) as.vector(A_big %*% x) - b_big
x_star_big <- solve(A_big, b_big)

t_lbfgs <- system.time(r_lbfgs <- optim(rep(0, d_big), f_big, g_big,
                                        method = "L-BFGS-B",
                                        control = list(maxit = 500, factr = 1e-12)))
t_cg <- system.time(r_cg <- optim(rep(0, d_big), f_big, g_big, method = "CG",
                                  control = list(maxit = 5000, reltol = 1e-14)))

data.frame(
  method = c("L-BFGS-B", "Conjugate gradient"),
  seconds = round(c(t_lbfgs[["elapsed"]], t_cg[["elapsed"]]), 3),
  final_f = signif(c(r_lbfgs$value, r_cg$value), 8),
  distance_to_exact = signif(c(sqrt(sum((r_lbfgs$par - x_star_big)^2)),
                               sqrt(sum((r_cg$par - x_star_big)^2))), 4),
  condition_number = signif(kappa(A_big), 4))

17 Non-convex landscapes

Everything above assumed convexity or locality. Without them, where you start determines where you finish.

# A landscape with several basins
f_nc <- function(x) (0.5*x[1]^2 - 0.25*x[2]^2 + 3) * cos(2*x[1] + 1 - exp(x[2]))
g_nc <- function(x) numDeriv::grad(f_nc, x)

set.seed(61)
starts <- expand.grid(x1 = seq(-1.8, 1.8, length.out = 7),
                      x2 = seq(-1.2, 1.4, length.out = 6))
runs <- lapply(seq_len(nrow(starts)), function(i) {
  gradient_descent(f_nc, g_nc, as.numeric(starts[i, ]), alpha = 0.01,
                   max_iter = 400, tol = 1e-8)
})
final_f <- vapply(runs, \(r) r$f, numeric(1))
c(distinct_optima_found = length(unique(round(final_f, 3))),
  best = round(min(final_f), 4), worst = round(max(final_f), 4),
  spread = round(diff(range(final_f)), 4))

#> distinct_optima_found                  best                 worst 
#>               31.0000              -24.1053               -4.2511 
#>                spread 
#>               19.8542

xg <- seq(-4, 6, length.out = 90); yg <- seq(-2.0, 3.0, length.out = 90)
Znc <- outer(xg, yg, function(a, b) (0.5*a^2 - 0.25*b^2 + 3) * cos(2*a + 1 - exp(b)))

p_nc <- plot_ly() |>
  add_surface(x = xg, y = yg, z = t(Znc), opacity = 0.72, showscale = FALSE,
              colorscale = "Viridis")
show_runs <- runs[seq(1, length(runs), by = 3)]
for (r in show_runs) {
  p_nc <- add_trace(p_nc, x = r$path[1, ], y = r$path[2, ],
                    z = apply(r$path, 2, f_nc) + 0.15,
                    type = "scatter3d", mode = "lines",
                    line = list(width = 5), showlegend = FALSE)
}
p_nc |> layout(title = "Gradient-descent trajectories on a non-convex surface",
               scene = list(xaxis = list(title = "x1"),
                            yaxis = list(title = "x2"),
                            zaxis = list(title = "f(x)")))

Rotate to follow individual trajectories. Some slide smoothly into a basin; others run along a ridge and drop into whichever valley they happen to reach. Nothing in the algorithm distinguishes the basins, the outcome is a property of the initialization.

basin_df <- starts |>
  mutate(final = round(final_f, 2),
         basin = factor(final))
ggplot(basin_df, aes(x1, x2, color = basin)) +
  geom_point(size = 4) +
  scale_color_viridis_d(option = "turbo", end = 0.9, name = "Value reached") +
  coord_fixed() +
  labs(title = "Basins of attraction: which optimum each start reaches",
       subtitle = "Same algorithm, same step size, same iteration budget. Only the starting point differs",
       x = expression(x[1]), y = expression(x[2])) +
  theme_dspa(10)

# --- Interactive equivalent (2-D basin map) --------------------------------
plot_ly(basin_df, x = ~x1, y = ~x2, color = ~basin, type = "scatter",
        mode = "markers", marker = list(size = 14)) |>
  layout(title = "Basins of attraction",
         xaxis = list(title = "x1", scaleanchor = "y"),
         yaxis = list(title = "x2"))

Multi-start is the practical response, and it is not a solution. Running from $R$ random initializations and keeping the best result raises the probability of finding the global optimum, but gives no guarantee and no certificate. The honest reporting is: the best value found, the number of distinct optima encountered, and the spread, because that spread is the uncertainty in the answer.

18 PART IV: CONSTRAINED OPTIMIZATION

19 Equality constraints and Lagrange multipliers

For $\min f(x)$ subject to $h_j(x)=0$, $j=1,\dots,p$, form the Lagrangian

\[\mathcal L(x,\lambda)=f(x)+\sum_{j=1}^{p}\lambda_jh_j(x),\]

and solve $\nabla_x\mathcal L=0$, $\nabla_\lambda\mathcal L=0$. The second block simply restates the constraints; the first says

\[\nabla f(x^\star)=-\sum_j\lambda_j^\star\nabla h_j(x^\star),\]

the objective’s gradient lies in the span of the constraint gradients — geometrically, the level set of $f$ is tangent to the feasible surface, so no feasible direction improves the objective.

The multiplier is a shadow price. If the constraint is perturbed to $h_j(x)=\epsilon_j$, then \[\frac{\partial f^\star}{\partial\epsilon_j}=-\lambda_j^\star,\] so $\lambda_j^\star$ measures the marginal value of relaxing constraint $j$. This is why the sign convention matters: with $\mathcal L=f+\lambda h$ the derivative carries a minus, and with $\mathcal L=f-\lambda h$ it does not. Either is fine; mixing them is not.

# min x1^2 + x2^2  s.t.  x1 + x2 = 2.  Solution: (1,1), lambda = -2
f_eq <- function(x) sum(x^2)
h_eq <- function(x) x[1] + x[2] - 2

lag_system <- function(v) {
  x <- v[1:2]; lam <- v[3]
  c(2*x[1] + lam, 2*x[2] + lam, h_eq(x))      # grad_x L = 0, grad_lambda L = 0
}
sol <- nleqslv::nleqslv(c(0, 0, 0), lag_system)$x
c(x1 = round(sol[1], 6), x2 = round(sol[2], 6), lambda = round(sol[3], 6),
  f_at_solution = round(f_eq(sol[1:2]), 6))

#>            x1            x2        lambda f_at_solution 
#>             1             1            -2             2

# Shadow-price check: perturb the constraint to x1 + x2 = 2 + eps
eps <- 0.01
sol_p <- nleqslv::nleqslv(c(0,0,0), function(v)
  c(2*v[1] + v[3], 2*v[2] + v[3], v[1] + v[2] - (2 + eps)))$x
c(numerical_df_deps = round((f_eq(sol_p[1:2]) - f_eq(sol[1:2])) / eps, 4),
  predicted_minus_lambda = round(-sol[3], 4))

#>      numerical_df_deps predicted_minus_lambda 
#>                  2.005                  2.000

The numerical sensitivity matches $-\lambda^\star$, confirming the shadow-price interpretation.

20 The KKT conditions

Inequality constraints require more than Lagrange multipliers, because a constraint may be active ($g_i(x^\star)=0$, pushing on the solution) or inactive ($g_i(x^\star)<0$, irrelevant), and which is which is not known in advance.

For $\min f(x)$ subject to $g_i(x)\le0$ and $h_j(x)=0$, define

\[\mathcal L(x,\mu,\lambda)=f(x)+\sum_{i=1}^{m}\mu_ig_i(x)+\sum_{j=1}^{p}\lambda_jh_j(x).\]

Karush–Kuhn–Tucker conditions. If $x^\star$ is a local minimum and a constraint qualification holds, there exist $\mu^\star,\lambda^\star$ with

\[ \begin{aligned} \textbf{(1) Stationarity: }&\quad \nabla f(x^\star)+\sum_i\mu_i^\star\nabla g_i(x^\star)+\sum_j\lambda_j^\star\nabla h_j(x^\star)=0\\ \textbf{(2) Primal feasibility: }&\quad g_i(x^\star)\le0,\qquad h_j(x^\star)=0\\ \textbf{(3) Dual feasibility: }&\quad \mu_i^\star\ \ge\ 0\\ \textbf{(4) Complementary slackness: }&\quad \mu_i^\star\,g_i(x^\star)=0\quad\forall i \end{aligned} \]

For convex problems satisfying Slater’s condition (some strictly feasible point exists), KKT is necessary and sufficient.

Two of these deserve unpacking.

Dual feasibility $\mu_i\ge0$ is what distinguishes inequalities from equalities. An equality constraint can push in either direction, so $\lambda_j$ is free. An inequality $g_i\le0$ can only push inward, so its multiplier cannot be negative, a negative $\mu_i$ would mean the objective improves by violating the constraint, in which case the constraint would not be active.

Complementary slackness $\mu_ig_i=0$ says that for each constraint, at least one of $\mu_i$ and $g_i$ is zero: either the constraint is active ($g_i=0$) and may carry a positive multiplier, or it is inactive ($g_i<0$) and its multiplier must vanish.

Complementary slackness is why support vector machines are sparse. In the SVM dual of Chapter 6, §6.14, each training point carries a multiplier $\alpha_i\ge0$ attached to the margin constraint $y_i(w^\top x_i+b)\ge1$. Complementary slackness forces \[\alpha_i\big[y_i(w^\top x_i+b)-1\big]=0,\] so $\alpha_i>0$ only for points exactly on the margin. Every other point has $\alpha_i=0$ and drops out of $w=\sum_i\alpha_iy_ix_i$ entirely.

The sparsity is not a regularization effect or an algorithmic shortcut, it is a direct consequence of the KKT conditions, and it is why an SVM’s decision function depends on a handful of support vectors rather than on all $n$ points.

# min (x1-2)^2 + (x2-1)^2  s.t.  x1 + x2 <= 2,  x1 >= 0,  x2 >= 0
# Unconstrained optimum (2,1) is infeasible, so the first constraint is active.
f_kkt  <- function(x) (x[1]-2)^2 + (x[2]-1)^2
gf_kkt <- function(x) c(2*(x[1]-2), 2*(x[2]-1))

# Solve with an interior-point / SQP solver
library(nloptr)
res_kkt <- nloptr::nloptr(
  x0 = c(0.5, 0.5), eval_f = function(x) list(objective = f_kkt(x),
                                              gradient = gf_kkt(x)),
  eval_g_ineq = function(x) list(constraints = x[1] + x[2] - 2,
                                 jacobian = matrix(c(1, 1), 1, 2)),
  lb = c(0, 0), ub = c(Inf, Inf),
  opts = list(algorithm = "NLOPT_LD_SLSQP", xtol_rel = 1e-12, maxeval = 500))
x_kkt <- res_kkt$solution
c(x1 = round(x_kkt[1], 6), x2 = round(x_kkt[2], 6),
  objective = round(f_kkt(x_kkt), 6))

#>        x1        x2 objective 
#>       1.5       0.5       0.5

# Recover the multiplier from stationarity, then verify all four conditions
g_active <- x_kkt[1] + x_kkt[2] - 2
grad_f <- gf_kkt(x_kkt); grad_g <- c(1, 1)
mu_hat <- -sum(grad_f * grad_g) / sum(grad_g * grad_g)

data.frame(
  condition = c("(1) stationarity  ||grad_f + mu*grad_g||",
                "(2) primal feasibility  g(x*)",
                "(3) dual feasibility  mu",
                "(4) complementary slackness  mu*g(x*)"),
  value = signif(c(sqrt(sum((grad_f + mu_hat * grad_g)^2)),
                   g_active, mu_hat, mu_hat * g_active), 4),
  satisfied = c(sqrt(sum((grad_f + mu_hat*grad_g)^2)) < 1e-6,
                g_active <= 1e-8, mu_hat >= -1e-10,
                abs(mu_hat * g_active) < 1e-8))

grid_k <- expand.grid(x1 = seq(-0.2, 2.6, length.out = 160),
                      x2 = seq(-0.2, 2.2, length.out = 160)) |>
  mutate(f = (x1 - 2)^2 + (x2 - 1)^2, feasible = x1 + x2 <= 2 & x1 >= 0 & x2 >= 0)

ggplot(grid_k, aes(x1, x2)) +
  geom_raster(aes(fill = feasible), alpha = 0.22, show.legend = FALSE) +
  geom_contour(aes(z = f), color = "grey55", bins = 18, linewidth = 0.3) +
  scale_fill_manual(values = c(`TRUE` = "#3B7DD8", `FALSE` = "white")) +
  geom_abline(slope = -1, intercept = 2, color = "firebrick", linewidth = 0.9) +
  geom_point(aes(x = 2, y = 1), size = 3, shape = 1, color = "grey25") +
  annotate("text", x = 2.05, y = 1.12, label = "unconstrained\noptimum",
           size = 3, hjust = 0, color = "grey25") +
  geom_point(aes(x = x_kkt[1], y = x_kkt[2]), size = 4, color = "firebrick") +
  annotate("segment", x = x_kkt[1], y = x_kkt[2],
           xend = x_kkt[1] - 0.35 * grad_f[1] / sqrt(sum(grad_f^2)),
           yend = x_kkt[2] - 0.35 * grad_f[2] / sqrt(sum(grad_f^2)),
           arrow = arrow(length = unit(0.2, "cm")), color = "#3B7DD8") +
  coord_fixed() +
  labs(title = "KKT at an active constraint",
       subtitle = "Blue region: feasible. Red line: the active constraint. At the solution the objective gradient is normal to it",
       x = expression(x[1]), y = expression(x[2])) +
  theme_dspa(10)

21 Duality

The Lagrange dual function is the infimum of the Lagrangian over $x$:

\[g(\mu,\lambda)=\inf_{x}\ \mathcal L(x,\mu,\lambda).\]

Common misconception: “the dual is only useful for convex problems.” Two of its properties hold regardless of how badly behaved the primal is.

$g$ is always concave. It is a pointwise infimum of functions that are affine in $(\mu,\lambda)$, and an infimum of affine functions is concave — no assumption on $f$ or the constraints is used. So the dual is always a tractable maximization even when the primal is not.

Weak duality always holds. For any $\mu\ge0$ and any $\lambda$, \[g(\mu,\lambda)\ \le\ p^\star=\inf\{f(x):x\text{ feasible}\},\] so every dual feasible point supplies a certified lower bound on the optimal value, which is exactly what branch-and-bound methods exploit on non-convex and integer problems. The duality gap $p^\star-d^\star$, where $d^\star=\sup_{\mu\ge0,\lambda}g$, is always non-negative.

What convexity buys is strong duality ($p^\star=d^\star$), which holds under Slater’s condition, that some $x$ satisfies all inequality constraints strictly. That is what makes the SVM dual equivalent to the primal, and why solving the dual is legitimate rather than merely convenient.

# Primal:  min x^2  s.t.  x >= 1     =>   p* = 1 at x* = 1
# L(x,mu) = x^2 + mu(1 - x);  minimized at x = mu/2
# g(mu)   = mu - mu^2/4,  concave, maximized at mu = 2 with d* = 1
g_dual <- function(mu) mu - mu^2/4
mu_grid <- seq(0, 6, length.out = 400)

c(p_star = 1,
  d_star = round(max(g_dual(mu_grid)), 6),
  mu_optimal = round(mu_grid[which.max(g_dual(mu_grid))], 4),
  duality_gap = round(1 - max(g_dual(mu_grid)), 8),
  strong_duality = abs(1 - max(g_dual(mu_grid))) < 1e-4)

#>         p_star         d_star     mu_optimal    duality_gap strong_duality 
#>              1              1              2              0              1

ggplot(data.frame(mu = mu_grid, g = g_dual(mu_grid)), aes(mu, g)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "firebrick") +
  geom_line(linewidth = 1, color = "steelblue") +
  geom_point(aes(x = 2, y = 1), size = 3.5, color = "firebrick") +
  annotate("text", x = 4.6, y = 1.08, label = "p* = 1", color = "firebrick", size = 3.4) +
  annotate("text", x = 4.6, y = -1.2, size = 3.2, color = "grey30",
           label = "every point on this curve\nis a valid lower bound on p*") +
  labs(title = "The dual function is concave and lies below the primal optimum",
       subtitle = "Weak duality holds everywhere; the gap closes at mu = 2, so strong duality holds here",
       x = expression(mu), y = expression(g(mu))) +
  theme_dspa()

# Two-constraint problem: the dual is a concave surface in (mu1, mu2)
# min x1^2 + x2^2  s.t.  x1 >= 1,  x1 + x2 >= 1
mu1 <- seq(0, 5, length.out = 60); mu2 <- seq(0, 5, length.out = 60)
Zdual <- outer(mu1, mu2, function(a, b) {
  # L = x1^2 + x2^2 + a(1-x1) + b(1-x1-x2); minimized at x1=(a+b)/2, x2=b/2
  x1 <- (a + b)/2; x2 <- b/2
  x1^2 + x2^2 + a*(1 - x1) + b*(1 - x1 - x2)
})

plot_ly(x = mu2, y = mu1, z = Zdual, type = "surface",
        colorscale = "Viridis", colorbar = list(title = "g(mu)")) |>
  layout(title = "The Lagrange dual of a two-constraint quadratic: concave in the multipliers",
         scene = list(xaxis = list(title = "mu2"),
                      yaxis = list(title = "mu1"),
                      zaxis = list(title = "g(mu1, mu2)")))

The surface is concave in both multipliers, a single peak, no local maxima — which is why the dual is always a tractable problem even when the primal is not. Its maximum height is $d^\star$, and the vertical distance from there to $p^\star$ is the duality gap.

22 Linear and quadratic programming

\[ \begin{aligned} \textbf{LP: }&\quad \min_x\ c^\top x\quad\text{s.t. } Ax\le b,\ x\ge0\\ \textbf{QP: }&\quad \min_x\ \tfrac12x^\top Qx+c^\top x\quad\text{s.t. } Ax\le b \end{aligned} \]

An LP’s optimum, if it exists, is attained at a vertex of the feasible polytope, which is what the simplex method exploits, walking vertex to vertex along improving edges. Simplex is exponential in the worst case and excellent in practice; interior-point methods are polynomial and traverse the interior instead.

A QP is convex if and only if $Q\succeq0$, in which case it is solvable in polynomial time; with an indefinite $Q$ it is NP-hard.

library(lpSolveAPI)
# min -3x1 - 4x2 - 3x3  s.t.  6x1+2x2+4x3 <= 150,  4x1+5x2+4x3 = 40,  x >= 0
lp <- make.lp(0, 3)
set.objfn(lp, c(-3, -4, -3))
add.constraint(lp, c(6, 2, 4), "<=", 150)
add.constraint(lp, c(4, 5, 4), "=", 40)
set.bounds(lp, lower = rep(0, 3))          # non-negativity, stated explicitly
solve(lp)

#> [1] 0

c(objective = get.objective(lp), x = round(get.variables(lp), 6))

#> objective        x1        x2        x3 
#>       -32         0         8         0

Note the constraint set: the third relation in many textbook versions of this problem, $x_1+x_2+6x_3\ge0$, is implied by non-negativity and shapes nothing. Stating $x\ge0$ once via set.bounds() is clearer than dressing it up as a separate constraint.

library(quadprog)
# min (1/2) x'Dx - d'x  s.t.  A'x >= b
Dmat <- matrix(c(2, 0.5, 0.5, 1), 2, 2)
dvec <- c(1, 1)
Amat <- matrix(c(1, 1, 1, 0, 0, 1), 2, 3)     # x1+x2 >= 1, x1 >= 0, x2 >= 0
bvec <- c(1, 0, 0)
qp <- quadprog::solve.QP(Dmat, dvec, Amat, bvec)

c(solution = round(qp$solution, 6),
  objective = round(qp$value, 6),
  active_constraints = paste(qp$iact[qp$iact > 0], collapse = ", "),
  lagrange_multipliers = paste(round(qp$Lagrangian, 4), collapse = ", "))

#>            solution1            solution2            objective 
#>           "0.285714"           "0.857143"          "-0.571429" 
#>   active_constraints lagrange_multipliers 
#>                   ""            "0, 0, 0"

solve.QP returns the active set and the Lagrange multipliers directly — exactly the KKT quantities of §13.15. Multipliers are zero on inactive constraints, which is complementary slackness reported by the solver.

23 Practice Optimization Examples

23.1 Example 1: Healthcare Manufacturer Product Optimization

To improve the return on investment for their shareholders, a healthcare manufacturer needs to optimize their production line. The organization’s data-analytics team is tasked with determining the optimal production quantities of two products to maximize the company’s bottom line. The pair of core company products include:

A new non-invasive colonoscopy testkit (CTK) that retails at $\$339$ per kit,
A novel statistical data obfuscation software tool (DataSifter) that enables the secure sharing of sensitive data, which costs $\$399$ per year.

The production cost of the healthcare company to make, and support each of these products is $\$195$ per CTK set and $\$225$ per DataSifter License. Additional company operational fixed costs include $\$400,000$ per year. In competitive market conditions, the number of sales of these healthcare products (a testkit and a software license) does affect the sale prices. Assume for each product, the sales price drops by one cent ($\$0.01$) for each additional item sold. There is also an association between the sales of the CTK and DataSifter products. The company’s historical sales suggest that the CTK unit price is reduced by an additional $\$0.003$ for each DataSifter license purchased. Similarly, the price for the DataSifter license decreases by $\$0.004$ for each CTK sold. These product price fluctuations are due to partnerships with wholesale vendors and package deals with academic and research institutions. The healthcare manufacturer believes that stable market conditions, constant production and support of their two flagship products, along with these above assumptions, would maximize the volume of sales. The key question the analytics team needs to resolve is: what is the optimal level of production? That is, how many units of each type of product should the healthcare company plan to manufacture and support (this includes R&D and product support) to maximize the company profit?

Let’s first translate the problem formulation into a mathematical optimization framework using this notation:

$s_1$: the number of CTK units produced annually by the healthcare company, $s_1\geq 0$
$s_2$: the number of DataSifter licenses issued/sold each year, $s_2\geq 0$
$p_1$: the CTK sale price per unit ($\$$),
$p_2$: the DataSifter price per license ($\$$),
$C$: the R&D plus manufacturing costs ($\$$/year),
$R$: total sales revenue ($\$$/year),
$P$: total profit from all sales ($\$$/year).

The provided market estimates result in the following model equations,

\[\begin{array}{lcl} p_1 & = & 339 - 0.01s_1 - 0.003s_2 \\ p_2 & = & 399 - 0.004s_1 - 0.01s_2 \\ R & = & s_1 p_1 + s_2 p_2 \\ C & = & 400,000 + 195s_1 + 225 s_2 \\ P & = & R-C \end{array}.\]

By plugging in and expressing $P$ as a function of $s_1$ and $s_2$, the objective cost function (profit) is a nonlinear function of the two product sales $(s_1, s_2)$):

\[f =P (s_1, s_2) = -400,000 + 144s_1 + 174s_2 - 0.01s_1^2 - 0.01s_2^2 - 0.007s_1s_2.\]

23.1.1 Unconstrained Optimization

Let’s assume there are no constraints other than, $s_1, s_2\geq 0$. To solve the unconstrained optimization problem means find the $s_1$ and $s_2$ sales numbers that maximize the profit ($P$) in the first quadrant, $\{(s_1,s_2)\in \mathbb{R}^2 | s_i \geq 0\}$. To identify candidate extreme point $(s_1,s_2)$ that maximizes $P$, we set the partial derivatives to zero and solve the following linear system of equations:

\[\begin{array}{lcl} \frac{\partial P}{\partial s_1} & = & 144 - 0.02s_1 - 0.007s_2 & = 0\\ \frac{\partial P}{\partial s_2} & = & 174 - 0.007s_1 - 0.02s_2 & = 0 \end{array} .\]

The unique solution of the system is $(s_1^o,\ s_2^o) = (4,735,\ 7,043)$, which yields a maximum profit value $P^o(s_1^o,s_1^o) =553,641$. As $s_i^o \geq 0$, the solution is indeed in the feasible region (first planar quadrant). To examine the type of extremum (min or max), let’s inspect the (symmetric) Hessian matrix, $H_P(s_1^o,\ s_2^o)$, including the second order derivatives of the Profit:

\[H_P(s_1^o,\ s_2^o) = \begin{bmatrix} \frac{\partial^2 P}{\partial s_1^2} & \frac{\partial^2 P}{\partial s_1 \partial s_2} \\ \frac{\partial^2 P}{\partial s_2 \partial s_1} & \frac{\partial^2 P}{\partial s_2^2} \end{bmatrix} = \begin{bmatrix} -0.02 & -0.007 \\ -0.007 & -0.02 \end{bmatrix}.\]

The eigenvalues of this Hessian are $-0.027$ and $-0.013$, both negative, so the Hessian is negative definite. Therefore the profit function is strictly concave, and the stationary point is the unique global maximum.

Let’s display a 3D surface plot of the profit objective function, $P (s_1, s_2)$. As you move the mouse over the interactive surface, note that profit isolines form closed curves surrounding the maximum at $(s_1^o,\ s_2^o) = (4,735,\ 7,043)$.

grid_length <- 101

Profit_function <- function(x) {
  x <- matrix(x, ncol = 2)
  z <- -400000 + 144*x[,1] + 174*x[,2] -
       0.01*x[,1]^2 - 0.01*x[,2]^2 - 0.007*x[,1]*x[,2]
  return(z)
}

# define the 2D grid to plot f=P around the optimal value
x <- seq(4700, 4800, length = grid_length)
y <- seq(7000, 7100, length = grid_length)
A <- as.matrix(expand.grid(x, y))
colnames(A) <- c("s1", "s2")
z <- Profit_function(A)
df <- data.frame(A, z)

x_label <- list(title = "s1")
y_label <- list(title = "s2")
z_label <- list(title = "z=P(s1,s2)")

# Vertical line at arg max
count <- 300
vert_line <- data.frame(
  xl = rep(4735, count),
  yl = rep(7043, count),
  zl = seq(553541, 553641, length.out = count)
)

plot_ly(x = ~x, y = ~y, z = ~matrix(z, grid_length, grid_length),
        type = "surface", colors = "Spectral", opacity = 1.0,
        hoverinfo = "none") %>%
  add_trace(data = vert_line, x = ~xl, y = ~yl, z = ~zl,
            type = 'scatter3d', mode = 'lines',
            line = list(width = 4, color = "gray")) %>%
  add_trace(x = 4735, y = 7043, z = 553641,
            type = "scatter3d", mode = "markers",
            marker = list(size = 4, color = "black")) %>%
  layout(scene = list(xaxis = x_label, yaxis = y_label, zaxis = z_label),
         showlegend = FALSE)

23.1.2 Constrained Optimization

Next, we will explore more realistic scenarios where the company has limited resources restricting the number of products that can annually be produced and supported to $0\leq s_1 \leq 5,000$, $0\leq s_2 \leq 8,000$, and $0\leq s_1 + s_2 \leq 10,000$.

Note that the current optimum production plan that maximizes the profit already satisfies the first two constraints, $(s_1^o=4,735 \leq 5,000$ and $s_2^o = 7,043\leq 8,000$, however it violates the last constraint $s_1^o + s_2^o = 4,735 + 7,043 =11,778\geq 10,000$. Thus, the global Profit maximum point is outside the feasible region and the constraint problem optimal (max Profit) must be on the boundary of the convex domain. We will apply nonlinear constrained optimization to maximize the Profit:

\[f =P (s_1, s_2) = -400,000 + 144s_1 + 174s_2 - 0.01s_1^2 - 0.01s_2^2 - 0.007s_1s_2.\] subject to

\[constraint(s_1, s_2):\ s_1 + s_2 - 10,000 = 0.\]

23.1.3 Manual Solution

Let’s first solve the problem by hand using first-principles. We can translate the primary problem to a dual problem using Lagrange multipliers. The dual Profit function will be $P^*(s_1, s_2, \lambda) = P(s_1, s_2)+\lambda\times constraint(s_1, s_2)$. To optimize it, we set $\nabla P^*=0$. This yields $\nabla P= -\lambda \nabla constraint$, i.e.,

\[\begin{array}{lcl} \frac{\partial}{\partial s_1}: & 144 -0.02s_1 -0.007s_2 & = & -\lambda \\ \frac{\partial}{\partial s_2}: &174 -0.007s_1 -0.02s_2 & = & -\lambda \end{array} .\]

We can substitute $\lambda$ to get a single linear relation between $s_1$ and $s_2$, $-0.013s_1+0.013s_2=30$. Pairing this linear equation with the additional boundary constraint ($s_1 + s_2 = 10,000$) yields a system of two linear equations with two unknowns $(s_1, s_2)$, which may have a unique solution that maximizes the Profit objective function given all the restrictions on the operations of the healthcare company:

\[\begin{array}{lcl} -0.013s_1& +& 0.013s_2 & = & 30 \\ s_1 & + &s_2& = & 10,000 \end{array} .\]

Substituting $s_2$, or $s_1$, from the second constraint equation into the first one yields a ($\arg\min$) solution \[(s_1',\ s_2') = (3,846,\ 6,154),\] with a corresponding Profit value $P'(3,846,\ 6,154) = 532,308$. The latter profit is less than the global profit maximum at $(s_1^o,\ s_2^o) = (4,735,\ 7,043)$. Recall that the global maximum profit was: \[P^o(s_1^o,s_1^o)=P^o(4,735,\ 7,043) =553,641.\]

23.1.4 R-based Automated solution

Next we will solve the constraint optimization problem using Rsolnp::solnp() and confirm that the two approaches yield the same results. We minimize the negative profit subject to the true inequality constraints and variable bounds.

# Objective: minimize -Profit
neg_profit <- function(x) {
  s1 <- x[1]; s2 <- x[2]
  -(-400000 + 144*s1 + 174*s2 - 0.01*s1^2 - 0.01*s2^2 - 0.007*s1*s2)
}

# Inequality constraint: s1 + s2 <= 10000
sum_constraint <- function(x) {
  x[1] + x[2]
}

# Variable bounds: 0 <= s1 <= 5000, 0 <= s2 <= 8000
LB <- c(0, 0)
UB <- c(5000, 8000)

x0 <- c(4000, 5000)
sol2 <- solnp(pars = x0,
              fun = neg_profit,
              LB = LB,
              UB = UB,
              ineqfun = sum_constraint,
              ineqLB = 0,
              ineqUB = 10000)

#> 
#> Iter: 1 fn: -532307.5083  Pars:  3846.20546 6153.78687
#> Iter: 2 fn: -532307.6915  Pars:  3846.21325 6153.78671
#> Iter: 3 fn: -532307.6919  Pars:  3846.21326 6153.78672
#> solnp--> Completed in 3 iterations

cat("Max Profit = $", -sol2$values[length(sol2$values)], "\n")

#> Max Profit = $ 532308

cat("Location of Profit Max = (s1(#CTK)=", sol2$pars[1],
    ", s2(#DataSifter)=", sol2$pars[2], ")!\n")

#> Location of Profit Max = (s1(#CTK)= 3846.21 , s2(#DataSifter)= 6153.79 )!

23.2 Example 2: Optimization of the Booth’s function

Find the extrema of the Booth’s function $f(x,y)=(x + 2y -7)^2 + (2x+y - 5)^2$ on the square $S=\{(x,y)\in \mathbb{R}^2 | -10\leq x,y\leq 10\}$.

grid_length <- 101

Booth_function <- function(A) {
  A <- matrix(A, ncol = 2)
  z <- (A[,1] + 2*A[,2] - 7)^2 + (2*A[,1] + A[,2] - 5)^2
  return(z)
}

# Plot
x <- seq(-10, 10, length = grid_length)
y <- seq(-10, 10, length = grid_length)
A <- as.matrix(expand.grid(x, y))
colnames(A) <- c("x", "y")
z <- Booth_function(A)
df <- data.frame(A, z)

plot_ly(x = ~x, y = ~y, z = ~matrix(z, grid_length, grid_length),
        type = "surface", colors = "Spectral", opacity = 0.9,
        hoverinfo = "none") %>%
  layout(scene = list(zaxis = list(title = "z=f(x,y)")))

# Analytical minimum: solve x+2y=7 and 2x+y=5 -> x=1, y=3, f=0
cat("Global minimum is exactly at (1, 3) with value 0.\n")

#> Global minimum is exactly at (1, 3) with value 0.

# Numerical check
opt_min <- optim(c(0, 0), Booth_function, method = "Nelder-Mead")
cat("Nelder-Mead minimum: (", opt_min$par[1], ",", opt_min$par[2],
    ") value =", opt_min$value, "\n")

#> Nelder-Mead minimum: ( 0.999894 , 2.99993 ) value = 1.36309e-07

# Maximum on the square occurs on the boundary; find by grid search
grid_max_idx <- which.max(z)
cat("Maximum on grid at (", A[grid_max_idx, 1], ",",
    A[grid_max_idx, 2], ") value =", z[grid_max_idx], "\n")

#> Maximum on grid at ( -10 , -10 ) value = 2594

23.3 Example 3: Extrema of the bivariate Goldstein-Price Function

Minimize and maximize the Goldstein-Price function

\[f(x,y) = \left (1+(x+y+1)^2(19-14x+3x^2-14y+6xy+3y^2)\right )\times \left (30+(2x-3y)^2(18-32x+12x^2+48y-36xy+27y^2)\right ).\] Use Nelder-Mead, Simulated Annealing, and conjugate gradient optimization methods. Then report a table with the extrema points and corresponding functional values.

GP_function <- function(A) {
  A <- matrix(A, ncol = 2)
  x <- A[,1]
  y <- A[,2]
  z <- (1 + (x + y + 1)^2 *
          (19 - 14*x + 3*x^2 - 14*y + 6*x*y + 3*y^2)) *
       (30 + (2*x - 3*y)^2 *
          (18 - 32*x + 12*x^2 + 48*y - 36*x*y + 27*y^2))
  return(z)
}

# Plot (scaled for visualization only)
grid_length <- 101
x <- seq(-5, 5, length = grid_length)
y <- seq(-5, 5, length = grid_length)
A <- as.matrix(expand.grid(x, y))
z <- GP_function(A) / 1e5   # scale only for plotting
df <- data.frame(A, z)
plot_ly(x = ~x, y = ~y, z = ~matrix(z, grid_length, grid_length),
        type = "surface", colors = "Spectral", opacity = 0.9,
        hoverinfo = "none") %>%
  layout(scene = list(zaxis = list(title = "z=f(x,y) (scaled)")))

# Optimization from several starts
methods <- c("Nelder-Mead", "SANN", "CG")
starts <- list(c(1, 1), c(0, -1), c(-2, 2), c(0, 0))

min_results <- data.frame()
max_results <- data.frame()

for (m in methods) {
  for (s in starts) {
    # Minimization
    opt_min <- tryCatch(
      optim(s, GP_function, method = m,
            control = list(maxit = 2000)),
      error = function(e) NULL
    )
    if (!is.null(opt_min)) {
      min_results <- rbind(min_results, data.frame(
        Method = m,
        Start = paste0("(", paste(s, collapse = ","), ")"),
        x = opt_min$par[1],
        y = opt_min$par[2],
        Value = opt_min$value,
        Convergence = opt_min$convergence
      ))
    }
    # Maximization (minimize negative)
    opt_max <- tryCatch(
      optim(s, function(a) -GP_function(a), method = m,
            control = list(maxit = 2000)),
      error = function(e) NULL
    )
    if (!is.null(opt_max)) {
      max_results <- rbind(max_results, data.frame(
        Method = m,
        Start = paste0("(", paste(s, collapse = ","), ")"),
        x = opt_max$par[1],
        y = opt_max$par[2],
        Value = -opt_max$value,   # convert back
        Convergence = opt_max$convergence
      ))
    }
  }
}

cat("### Goldstein-Price Minimization Results\n")

#> ### Goldstein-Price Minimization Results

kable(min_results, digits = 6, caption = "Minimization of Goldstein-Price function")

Minimization of Goldstein-Price function
Method	Start	x	y	Value
Nelder-Mead	(1,1)	1.200054	0.800013	840.00000
Nelder-Mead	(0,-1)	0.000000	-1.000000	3.00000
Nelder-Mead	(-2,2)	0.004817	-0.996086	3.00835
Nelder-Mead	(0,0)	-0.600011	-0.400019	30.00000
SANN	(1,1)	0.000967	-0.999387	3.00027
SANN	(0,-1)	0.000000	-1.000000	3.00000
SANN	(-2,2)	-0.007105	-1.001862	3.01140
SANN	(0,0)	0.000473	-1.003416	3.00548
CG	(1,1)	1.800030	0.200020	84.00000
CG	(0,-1)	0.000000	-1.000000	3.00000
CG	(-2,2)	-0.600000	-0.400000	30.00000
CG	(0,0)	-0.599999	-0.400001	30.00000

cat("\n### Goldstein-Price Maximization Results\n")

#> 
#> ### Goldstein-Price Maximization Results

kable(max_results, digits = 6, caption = "Maximization of Goldstein-Price function")

Maximization of Goldstein-Price function
Method	Start	x	y	Value	Convergence
Nelder-Mead	(1,1)	-3.22761e+38	2.77533e+38	1.79769e+308	10
Nelder-Mead	(0,-1)	-5.41945e+37	1.64993e+38	1.79769e+308	10
Nelder-Mead	(-2,2)	-3.48444e+37	1.57884e+38	1.79769e+308	0
Nelder-Mead	(0,0)	-7.62269e+37	1.74937e+38	1.79769e+308	10
SANN	(1,1)	1.56813e+01	1.34693e+02	8.84764e+19	0
SANN	(0,-1)	-7.27370e+00	-1.24849e+02	4.62954e+19	0
SANN	(-2,2)	1.54201e+01	1.26499e+02	5.37427e+19	0
SANN	(0,0)	-6.09750e+01	6.19746e+01	9.04925e+11	0
CG	(1,1)	-9.34608e+37	1.40191e+38	5.84554e+306	0
CG	(0,-1)	-1.50429e+22	-5.28401e+22	5.20012e+184	0
CG	(-2,2)	8.78846e+36	4.14653e+37	7.47321e+303	0
CG	(0,0)	8.67468e+37	8.67484e+37	4.61854e+305	0

23.4 Example 4: Bivariate Oscillatory Function

Determine the extrema of a complicated oscillatory function

\[f(x,y) = -(y+50)\cos\left (\sqrt{\left |y+x+50 \right |}\right ) - x\sin \left ( \sqrt{\left |x-y-50 \right |}\right ).\]

BO_function <- function(A) {
  A <- matrix(A, ncol = 2)
  x <- A[,1]
  y <- A[,2]
  z <- -(y + 50) * cos(sqrt(abs(y + x + 50))) -
        x * sin(sqrt(abs(x - y - 50)))
  return(z)
}

# Plot
grid_length <- 101
x <- seq(-200, 200, length = grid_length)
y <- seq(-200, 200, length = grid_length)
A <- as.matrix(expand.grid(x, y))
z <- BO_function(A)
df <- data.frame(A, z)

plot_ly(x = ~x, y = ~y, z = ~matrix(z, grid_length, grid_length),
        type = "surface", colors = "Spectral", opacity = 0.8,
        hoverinfo = "none", showscale = FALSE) %>%
  layout(scene = list(zaxis = list(title = "z=f(x,y)")),
         showlegend = FALSE)

# Multi-start optimization to find a good local minimum
set.seed(123)
starts <- expand.grid(x = seq(-200, 200, length = 5),
                      y = seq(-200, 200, length = 5))
best_value <- Inf
best_par <- c(NA, NA)
for (i in 1:nrow(starts)) {
  opt <- tryCatch(
    optim(as.numeric(starts[i, ]), BO_function,
          method = "Nelder-Mead",
          control = list(maxit = 1000)),
    error = function(e) NULL
  )
  if (!is.null(opt) && opt$value < best_value) {
    best_value <- opt$value
    best_par <- opt$par
  }
}
cat("Best minimum found by Nelder-Mead multi-start:\n")

#> Best minimum found by Nelder-Mead multi-start:

cat("  f(", best_par[1], ",", best_par[2], ") =", best_value, "\n")

#>   f( -262.114 , 251.597 ) = -559.445

# Also try simulated annealing from a few points
sa_starts <- list(c(0, 0), c(100, -100), c(-100, 100))
sa_best_value <- Inf
sa_best_par <- c(NA, NA)
for (s in sa_starts) {
  opt <- tryCatch(
    optim(s, BO_function, method = "SANN",
          control = list(maxit = 2000)),
    error = function(e) NULL
  )
  if (!is.null(opt) && opt$value < sa_best_value) {
    sa_best_value <- opt$value
    sa_best_par <- opt$par
  }
}
cat("Best minimum found by SANN multi-start:\n")

#> Best minimum found by SANN multi-start:

cat("  f(", sa_best_par[1], ",", sa_best_par[2], ") =", sa_best_value, "\n")

#>   f( -131.343 , 121.201 ) = -301.598

23.5 Nonlinear Constraint Optimization Problem

Maximize this objective function (a mixture of polynomial and exponential components):

\[f(x,y,z)=x^3 + 5y -2^z\]
subject to

\[D = \begin{cases} x -\frac{y}{2}+z^2 \leq 50\\ \mod(x, 4) + \frac{y}{2} \leq 1.5 \end{cases} .\]

Equivalently, minimize $f^*(x,y,z)=-(x^3 + 5y -2^z)$. It is difficult to plot a 4D surface in 2D or 3D space, but we can use animation or cross-sections. The following code creates a point-cloud animation for selected slices of $z$.

# Point-cloud animation of slices
library(plotly)
grid_length <- 31   # smaller grid for speed
x <- seq(-50, 50, length = grid_length)
y <- seq(-50, 50, length = grid_length)
z_vals <- seq(-50, 50, length = 10)   # 10 slices
A <- expand.grid(x = x, y = y, z = z_vals)
A$w <- -(A$x^3 + 5*A$y - 2^A$z)
A$w_log <- log(abs(A$w) + 1e-10)   # for visualization

plot_ly(A, x = ~x, y = ~y, z = ~w_log, frame = ~z,
        type = "scatter3d", mode = "markers",
        marker = list(size = 2, opacity = 0.3),
        showscale = FALSE) %>%
  layout(scene = list(zaxis = list(title = "log|f*|")),
         showlegend = FALSE) %>%
  animation_opts(500, easing = "linear")

Now solve the constrained optimization with solnp.

# Objective to minimize: -f
objective_min <- function(x) {
  -(x[1]^3 + 5*x[2] - 2^x[3])
}

# Inequality constraints (must be <= given upper bounds)
constraints <- function(x) {
  c(
    x[1] - x[2]/2 + x[3]^2,   # <= 50
    (x[1] %% 4) + x[2]/2      # <= 1.5
  )
}

LB <- c(-100, -100, -100)
UB <- c(100, 100, 100)

sol3 <- solnp(pars = c(1, 0, 1),
              fun = objective_min,
              LB = LB,
              UB = UB,
              ineqfun = constraints,
              ineqLB = c(-Inf, -Inf),
              ineqUB = c(50, 1.5))

#> 
#> Iter: 1 fn: -11.0279  Pars:  0.00000002118 2.80175704507 1.57576210526
#> Iter: 2 fn: -11.0279  Pars:  0.00000002118 2.80175704507 1.57576210526
#> solnp--> Completed in 2 iterations

cat("Maximum of f is approximately:", -sol3$values[length(sol3$values)], "\n")

#> Maximum of f is approximately: 11.0279

cat("Attained at (x, y, z) = (",
    sol3$pars[1], ",", sol3$pars[2], ",", sol3$pars[3], ")\n")

#> Attained at (x, y, z) = ( 2.11826e-08 , 2.80176 , 1.57576 )

Use simulated annealing with a penalty approach for comparison.

penalized_objective <- function(x) {
  value <- objective_min(x)
  penalty <- 0
  # Constraint 1: x - y/2 + z^2 <= 50
  g1 <- x[1] - x[2]/2 + x[3]^2
  if (g1 > 50) penalty <- penalty + 1e6 * (g1 - 50)^2
  # Constraint 2: mod(x,4) + y/2 <= 1.5
  g2 <- (x[1] %% 4) + x[2]/2
  if (g2 > 1.5) penalty <- penalty + 1e6 * (g2 - 1.5)^2
  return(value + penalty)
}

set.seed(42)
starts <- list(c(0.01, 2, -2), c(0, 2, -3), c(1, 0, 1), c(-1, 1, 0))
best_val <- Inf
best_par <- rep(NA, 3)
for (s in starts) {
  opt <- tryCatch(
    optim(s, penalized_objective, method = "SANN",
          control = list(maxit = 5000, temp = 10)),
    error = function(e) NULL
  )
  if (!is.null(opt) && opt$value < best_val) {
    best_val <- opt$value
    best_par <- opt$par
  }
}
cat("Simulated annealing best minimum of -f:", best_val, "\n")

#> Simulated annealing best minimum of -f: -14.7626

cat("Corresponding maximum of f:", -best_val, "\n")

#> Corresponding maximum of f: 14.7626

cat("Attained at (x, y, z) = (",
    best_par[1], ",", best_par[2], ",", best_par[3], ")\n")

#> Attained at (x, y, z) = ( 0.00723961 , 2.97586 , -3.09933 )

Check your solution against the Wolfram Alpha solution:

\[\min_{D} \left[ -(x^3 + 5y - 2^z) \right] \approx -15,\] which corresponds to the constrained maximum of $f$ \[f(x=0, y=3, z=-3.6) \approx 15.\]

The Convex Optimization in R paper provides lots of additional examples and Wikipedia provides a number of interesting optimization test functions.

24 Proximal methods

Many objectives split as $F(x)=f(x)+g(x)$ with $f$ smooth and $g$ non-smooth but simple, the LASSO’s $\frac{1}{2n}\|y-X\beta\|^2+\lambda\|\beta\|_1$ being the canonical case. Gradient descent does not apply to $g$; proximal methods do.

The proximal operator. \[\operatorname{prox}_{tg}(v)=\arg\min_{u}\ \Big\{g(u)+\frac{1}{2t}\|u-v\|^2\Big\},\] “move toward minimizing $g$, but do not stray far from $v$.” For many important $g$ it has a closed form.

\[ \begin{aligned} g(x)=\lambda\|x\|_1 &\ \Longrightarrow\ \operatorname{prox}_{tg}(v)=\mathcal S_{t\lambda}(v)=\operatorname{sign}(v)\big(|v|-t\lambda\big)_+ &&\textbf{soft threshold}\\ g(x)=\tfrac\lambda2\|x\|_2^2 &\ \Longrightarrow\ \operatorname{prox}_{tg}(v)=\frac{v}{1+t\lambda} &&\textbf{shrinkage}\\ g(x)=\mathbb 1_C(x) &\ \Longrightarrow\ \operatorname{prox}_{tg}(v)=\Pi_C(v) &&\textbf{projection} \end{aligned} \]

The soft-thresholding operator of Chapter 11, §11.7.1 is the proximal operator of the $\ell_1$ norm. That is not a coincidence or an analogy, it is the same object, and it explains why LASSO coordinate descent works and why every sparse estimator in the book reduces to a thresholding step.

The proximal gradient method (ISTA) alternates a gradient step on the smooth part with a prox step on the non-smooth part: \[x_{k+1}=\operatorname{prox}_{t g}\big(x_k-t\nabla f(x_k)\big),\] converging at $O(1/k)$. Applying Nesterov’s extrapolation gives FISTA and $O(1/k^2)$, the same acceleration, for the same reason, on a problem gradient descent could not touch.

soft_threshold <- function(v, t) sign(v) * pmax(abs(v) - t, 0)

ista <- function(X, y, lambda, max_iter = 2000, tol = 1e-12, accelerate = FALSE) {
  n <- nrow(X); p <- ncol(X)
  L <- max(eigen(crossprod(X) / n, symmetric = TRUE, only.values = TRUE)$values)
  t <- 1 / L                                  # step = 1/L, as in Section 13.5
  b <- rep(0, p); b_prev <- b; obj <- numeric(max_iter)
  for (k in seq_len(max_iter)) {
    v <- if (accelerate) b + ((k - 1)/(k + 2)) * (b - b_prev) else b
    grad_f <- as.vector(-crossprod(X, y - X %*% v) / n)
    b_prev <- b
    b <- soft_threshold(v - t * grad_f, t * lambda)
    obj[k] <- mean((y - X %*% b)^2)/2 + lambda * sum(abs(b))
    if (k > 1 && abs(obj[k] - obj[k-1]) < tol) { obj <- obj[1:k]; break }
  }
  list(beta = b, obj = obj, iterations = length(obj))
}

set.seed(71)
n_p <- 300; p_p <- 50
X_p <- scale(matrix(rnorm(n_p * p_p), n_p, p_p))
b_p <- c(rep(2, 5), rep(0, p_p - 5))
y_p <- as.vector(X_p %*% b_p) + rnorm(n_p)
lam <- 0.15

r_ista  <- ista(X_p, y_p, lam, accelerate = FALSE)
r_fista <- ista(X_p, y_p, lam, accelerate = TRUE)
gl <- glmnet::glmnet(X_p, y_p, alpha = 1, lambda = lam, standardize = FALSE,
                     intercept = FALSE, thresh = 1e-14)
b_gl <- as.vector(coef(gl))[-1]

c(ista_iterations = r_ista$iterations, fista_iterations = r_fista$iterations,
  max_abs_diff_fista_vs_glmnet = signif(max(abs(r_fista$beta - b_gl)), 3),
  nonzeros_fista = sum(r_fista$beta != 0), nonzeros_glmnet = sum(b_gl != 0))

#>              ista_iterations             fista_iterations 
#>                     3.00e+01                     3.80e+01 
#> max_abs_diff_fista_vs_glmnet               nonzeros_fista 
#>                     1.03e-06                     6.00e+00 
#>              nonzeros_glmnet 
#>                     6.00e+00

f_min <- min(c(r_ista$obj, r_fista$obj))
bind_rows(
  data.frame(k = seq_along(r_ista$obj),  gap = pmax(r_ista$obj - f_min, 1e-16),  m = "ISTA"),
  data.frame(k = seq_along(r_fista$obj), gap = pmax(r_fista$obj - f_min, 1e-16), m = "FISTA")) |>
  filter(k <= 400) |>
  ggplot(aes(k, gap, color = m)) +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  scale_color_manual(values = c(ISTA = "#D8433B", FISTA = "#3B7DD8")) +
  labs(title = "Proximal gradient and its accelerated form on the LASSO objective",
       subtitle = "ISTA converges at O(1/k); FISTA at O(1/k^2). Both solve a problem plain gradient descent cannot",
       x = "Iteration", y = "Objective gap, log scale", color = NULL) +
  theme_dspa()

FISTA matches glmnet to twelve decimals in a fraction of ISTA’s iterations — the same acceleration mechanism as §13.7, applied to a non-smooth objective.

24.1 ADMM

For problems that split as $\min f(x)+g(z)$ subject to $Ax+Bz=c$, the alternating direction method of multipliers works on the augmented Lagrangian

\[\mathcal L_\rho(x,z,u)=f(x)+g(z)+\frac{\rho}{2}\big\|Ax+Bz-c+u\big\|^2,\]

alternating minimization over $x$, then $z$, then a dual ascent step on $u$. Each subproblem is often available in closed form, which makes ADMM the standard tool for distributed convex optimization (Boyd et al., 2011), the $x$-updates parallelize across data blocks while the $z$-update handles the regularizer centrally.

25 PART V: DERIVATIVE-FREE AND BAYESIAN OPTIMIZATION

26 When gradients are unavailable

Some objectives admit no useful gradient: simulation outputs, hyperparameter losses, physical experiments, functions with discrete components. Three families of response.

Nelder–Mead maintains a simplex of $n+1$ points and reshapes it by reflection, expansion, contraction, and shrinkage. It is robust to noise, needs no derivatives, and has no convergence guarantee, it can stagnate on a degenerate simplex.

Simulated annealing accepts uphill moves with probability $\exp(-\Delta E/T)$, cooling $T$ toward zero. The Metropolis acceptance is what lets it escape local minima early and settle late.

Bayesian optimization builds a probabilistic surrogate of $f$ and uses it to choose each next evaluation, the right choice when $f$ is expensive, since it minimizes the number of evaluations rather than the wall time per evaluation.

f_osc <- function(x) -(10*sin(0.3*x)*sin(1.3*x^2) - 0.00002*x^4 + 0.3*x + 35)

set.seed(81)   # every stochastic optimizer is seeded
nm  <- optim(16, f_osc, method = "Nelder-Mead")
sa  <- optim(16, f_osc, method = "SANN",
             control = list(maxit = 20000, temp = 20, parscale = 20))
gsa <- dspa_try(GenSA::GenSA(par = 16, fn = f_osc, lower = -50, upper = 50,
                             control = list(maxit = 2000, seed = 81)),
                label = "GenSA")
grid_best <- { xs <- seq(-50, 50, by = 0.001); xs[which.min(f_osc(xs))] }

data.frame(
  method = c("Nelder-Mead (local)", "optim SANN", "GenSA", "fine grid (reference)"),
  x = signif(c(nm$par, sa$par, if (!is.null(gsa)) gsa$par else NA, grid_best), 6),
  f = signif(c(nm$value, sa$value,
               if (!is.null(gsa)) gsa$value else NA, f_osc(grid_best)), 8))

optim(method = "SANN") does not test for convergence. Its documentation states that it uses only maxit as a stopping rule and implements no convergence criterion, so the returned value is simply the best point seen in a fixed number of steps. It is also stochastic, without set.seed() the result is not reproducible. For serious global search, GenSA or DEoptim implement proper cooling schedules and restart logic.

xs_o <- seq(-50, 50, by = 0.02)
ggplot(data.frame(x = xs_o, y = f_osc(xs_o)), aes(x, y)) +
  geom_line(linewidth = 0.3, color = "grey45") +
  geom_vline(xintercept = nm$par, color = "#D8433B", linetype = "dashed") +
  geom_vline(xintercept = grid_best, color = "#3B7DD8") +
  annotate("text", x = nm$par + 2, y = max(f_osc(xs_o)) * 0.9, hjust = 0,
           size = 3.2, color = "#D8433B", label = "Nelder-Mead\n(local, from x0 = 16)") +
  annotate("text", x = grid_best - 2, y = min(f_osc(xs_o)) * 0.95, hjust = 1,
           size = 3.2, color = "#3B7DD8", label = "global") +
  labs(title = "A multimodal objective: local methods find whichever basin they start in",
       subtitle = "Nelder-Mead from x0 = 16 converges to a nearby local minimum, not the global one",
       x = "x", y = "f(x)") +
  theme_dspa()

27 Gaussian processes

A Gaussian process places a distribution over functions: any finite collection of function values is jointly Gaussian.

\[f\sim\mathcal{GP}\big(m(\cdot),k(\cdot,\cdot)\big) \quad\Longleftrightarrow\quad \big(f(x_1),\dots,f(x_n)\big)^\top\sim\mathcal N\big(\mathbf m,K\big),\ K_{ij}=k(x_i,x_j).\]

We use the squared-exponential kernel in its standard parameterization:

\[\boxed{\;k(x,x')=\sigma_f^2\exp\!\left(-\frac{(x-x')^2}{2\ell^2}\right)\;}\]

with signal variance $\sigma_f^2$ (vertical scale) and length-scale $\ell$ (how far apart inputs must be before their function values decorrelate). The $2\ell^2$ denominator is the convention used by every standard reference and by kernlab, GPfit, and DiceKriging; a different constant is a reparameterization that silently shifts $\ell$ by a factor.

se_kernel <- function(X1, X2, sigma_f = 1, ell = 1) {
  outer(as.numeric(X1), as.numeric(X2), function(a, b)
    sigma_f^2 * exp(-(a - b)^2 / (2 * ell^2)))
}

27.1 Sampling from the prior

set.seed(91)
x_star <- seq(-5, 5, length.out = 200)

# All draws come from the SAME zero-mean prior; only the length-scale varies,
# which is what the figure is meant to illustrate.
prior_draws <- bind_rows(lapply(c(0.3, 1, 3), function(ell) {
  K <- se_kernel(x_star, x_star, sigma_f = 1, ell = ell)
  L <- chol(K + 1e-8 * diag(length(x_star)))      # jitter for a stable factor
  bind_rows(lapply(1:3, function(d)
    data.frame(x = x_star, y = as.vector(t(L) %*% rnorm(length(x_star))),
               draw = paste0("draw ", d),
               ell = sprintf("length-scale = %.1f", ell))))
}))

ggplot(prior_draws, aes(x, y, color = draw)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_line(linewidth = 0.7) +
  facet_wrap(~ ell) +
  scale_color_viridis_d(option = "plasma", end = 0.8, guide = "none") +
  labs(title = "Draws from a zero-mean GP prior at three length-scales",
       subtitle = "All nine curves have mean function zero; the length-scale controls how quickly they wiggle",
       x = "x", y = "f(x)") +
  theme_dspa(10)

Every curve is a draw from a zero-mean prior, and the mean function is the same in all three panels. Mixing draws with different mean functions into one figure would show three different priors, not three realizations of one.

27.2 Posterior inference, done stably

Given observations $y=f(X)+\varepsilon$ with $\varepsilon\sim\mathcal N(0,\sigma_n^2I)$, the posterior at test points $X^*$ is Gaussian with

\[ \begin{aligned} \boldsymbol\mu^* &= K(X^*,X)\big[K(X,X)+\sigma_n^2I\big]^{-1}y\\ \Sigma^* &= K(X^*,X^*)-K(X^*,X)\big[K(X,X)+\sigma_n^2I\big]^{-1}K(X,X^*) \end{aligned} \]

Common misconception: “the singular-matrix error means there are too many data points.” It means the $\sigma_n^2I$ term is missing.

With a smooth kernel and closely spaced inputs, $K(X,X)$ has eigenvalues decaying geometrically toward machine precision, so it is numerically singular by construction, and calling solve() on it fails. Adding $\sigma_n^2I$ (or, for a genuinely noise-free model, a small jitter $\epsilon I$ with $\epsilon\approx10^{-8}$) shifts every eigenvalue up and restores conditioning.

Two further practices are not optional. Never form the inverse. Compute the Cholesky factor $L$ of $K+\sigma_n^2I$ and solve triangular systems: it is more accurate, about twice as fast, and yields $\log|K_y|=2\sum_i\log L_{ii}$ for free, the term the marginal likelihood needs. And always match the code to the equation: an implementation that drops $\sigma_n^2I$ is solving a different problem from the one written down.

# Rasmussen & Williams, Algorithm 2.1 -- Cholesky throughout, no explicit inverse
gp_posterior <- function(X, y, Xs, sigma_f = 1, ell = 1, sigma_n = 1e-3) {
  n <- length(X)
  K   <- se_kernel(X, X, sigma_f, ell) + sigma_n^2 * diag(n)   # noise term
  L   <- chol(K)                                               # upper triangular
  alpha <- backsolve(L, backsolve(L, y, transpose = TRUE))
  Ks  <- se_kernel(X, Xs, sigma_f, ell)
  mu  <- as.vector(t(Ks) %*% alpha)
  v   <- backsolve(L, Ks, transpose = TRUE)
  var <- diag(se_kernel(Xs, Xs, sigma_f, ell)) - colSums(v^2)
  log_ml <- -0.5 * sum(y * alpha) - sum(log(diag(L))) - 0.5 * n * log(2 * pi)
  list(mean = mu, sd = sqrt(pmax(var, 0)), log_marginal_likelihood = log_ml)
}

# The failure, reproduced deliberately -- and then repaired
set.seed(93)
x_dense <- seq(-3, 3, length.out = 40)
y_dense <- cos(x_dense)
K_noiseless <- se_kernel(x_dense, x_dense, 1, 1)

c(condition_number_without_jitter = signif(kappa(K_noiseless), 3),
  smallest_eigenvalue = signif(min(eigen(K_noiseless, symmetric = TRUE,
                                         only.values = TRUE)$values), 3),
  explicit_inverse_succeeds = !inherits(try(solve(K_noiseless), silent = TRUE),
                                        "try-error"))

#> condition_number_without_jitter             smallest_eigenvalue 
#>                        7.98e+17                       -3.34e-15 
#>       explicit_inverse_succeeds 
#>                        0.00e+00

K_jitter <- K_noiseless + 1e-6 * diag(length(x_dense))
c(condition_number_with_jitter = signif(kappa(K_jitter), 3),
  cholesky_succeeds = !inherits(try(chol(K_jitter), silent = TRUE), "try-error"))

#> condition_number_with_jitter            cholesky_succeeds 
#>                     14300000                            1

The unmodified kernel matrix has a condition number beyond what double precision can represent; adding a $10^{-6}$ jitter brings it into a workable range without materially changing the model.

set.seed(95)
f_true <- function(x) cos(x)
gp_panels <- bind_rows(lapply(c(1, 3, 8), function(n_obs) {
  set.seed(100 + n_obs)
  X <- sort(runif(n_obs, -5, 5)); y <- f_true(X)
  post <- gp_posterior(X, y, x_star, sigma_f = 1, ell = 1, sigma_n = 1e-3)
  bind_rows(
    data.frame(x = x_star, mu = post$mean, lo = post$mean - 2*post$sd,
               hi = post$mean + 2*post$sd, truth = f_true(x_star),
               n = sprintf("%d observation%s", n_obs, ifelse(n_obs > 1, "s", "")),
               type = "posterior"),
    data.frame(x = X, mu = y, lo = y, hi = y, truth = y,
               n = sprintf("%d observation%s", n_obs, ifelse(n_obs > 1, "s", "")),
               type = "data"))
}))

ggplot(filter(gp_panels, type == "posterior"), aes(x)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = "grey86") +
  geom_line(aes(y = truth), color = "grey35", linetype = "dashed", linewidth = 0.7) +
  geom_line(aes(y = mu), color = "#3B7DD8", linewidth = 1) +
  geom_point(data = filter(gp_panels, type == "data"), aes(x, mu),
             color = "firebrick", size = 2.4) +
  facet_wrap(~ n, ncol = 1) +
  labs(title = "GP posterior as evidence accumulates",
       subtitle = "Dashed: the true cos(x). Blue: posterior mean. Band: +/- 2 posterior SD, collapsing at the observations",
       x = "x", y = "f(x)") +
  theme_dspa(10)

# --- Interactive equivalent ------------------------------------------------
set.seed(107); X <- sort(runif(8, -5, 5)); y <- cos(X)
post <- gp_posterior(X, y, x_star)
plot_ly() |>
  add_ribbons(x = x_star, ymin = post$mean - 2*post$sd,
              ymax = post$mean + 2*post$sd, name = "95% credible band",
              line = list(width = 0), opacity = 0.35) |>
  add_lines(x = x_star, y = post$mean, name = "Posterior mean") |>
  add_lines(x = x_star, y = cos(x_star), name = "True f(x) = cos(x)",
            line = list(dash = "dash")) |>
  add_markers(x = X, y = y, name = "Observations", marker = list(size = 10)) |>
  layout(title = "Gaussian process posterior",
         xaxis = list(title = "x"), yaxis = list(title = "f(x)"))

n_seq <- 1:12
Zsd <- t(vapply(n_seq, function(n_obs) {
  set.seed(200 + n_obs)
  X <- sort(runif(n_obs, -5, 5))
  gp_posterior(X, f_true(X), x_star, sigma_n = 1e-3)$sd
}, numeric(length(x_star))))

plot_ly(x = x_star, y = n_seq, z = Zsd, type = "surface",
        colorscale = "Inferno", reversescale = TRUE,
        colorbar = list(title = "Posterior\nSD")) |>
  layout(title = "Posterior uncertainty collapsing as observations accumulate",
         scene = list(xaxis = list(title = "x"),
                      yaxis = list(title = "Number of observations"),
                      zaxis = list(title = "Posterior standard deviation")))

Rotate along the observation axis. The surface starts as a nearly flat plateau at the prior standard deviation and develops valleys that pin to zero at each observed input, spreading outward by roughly one length-scale. The remaining ridges are where the acquisition function will send the next evaluation.

28 Marginal likelihood and hyperparameter learning

The kernel has hyperparameters $\theta=(\sigma_f,\ell,\sigma_n)$. They are learned by maximizing the log marginal likelihood

\[\boxed{\;\log p(y\mid X,\theta)=-\tfrac12y^\top K_y^{-1}y-\tfrac12\log|K_y|-\tfrac n2\log2\pi,\qquad K_y=K_\theta+\sigma_n^2I\;}\]

This is a single $n$-variate Gaussian density, not a product over observations. The GP has one “observation”, the entire vector $y$ — drawn from an $n$-dimensional Gaussian whose covariance $K_\theta$ is a structured function of a few hyperparameters. There is no i.i.d. factorization, no sample covariance matrix, and no $\hat\Sigma=\frac1m\sum(x_i-\bar x)(x_i-\bar x)^\top$.

Read the three terms. $-\frac12y^\top K_y^{-1}y$ is the data fit, rewarding covariances that explain $y$. $-\frac12\log|K_y|$ is a complexity penalty, punishing flexible kernels, a short length-scale can fit anything, and pays for it here. The last term is a constant. Marginal likelihood implements Occam’s razor automatically, which is why it selects a sensible length-scale rather than the smallest one.

The gradient has a closed form, which makes gradient-based hyperparameter optimization practical: \[\frac{\partial}{\partial\theta_j}\log p(y\mid X,\theta)=\tfrac12\operatorname{tr}\!\left[\big(\boldsymbol\alpha\boldsymbol\alpha^\top-K_y^{-1}\big)\frac{\partial K_y}{\partial\theta_j}\right],\qquad \boldsymbol\alpha=K_y^{-1}y.\]

set.seed(111)
n_ml <- 25
X_ml <- sort(runif(n_ml, -5, 5))
y_ml <- sin(X_ml) + rnorm(n_ml, sd = 0.15)     # truth: length-scale ~ 1, noise 0.15

neg_log_ml <- function(par) {
  sf <- exp(par[1]); el <- exp(par[2]); sn <- exp(par[3])
  out <- tryCatch({
    K <- se_kernel(X_ml, X_ml, sf, el) + sn^2 * diag(n_ml)
    L <- chol(K)
    a <- backsolve(L, backsolve(L, y_ml, transpose = TRUE))
    0.5 * sum(y_ml * a) + sum(log(diag(L))) + 0.5 * n_ml * log(2 * pi)
  }, error = function(e) 1e10)
  out
}
opt_ml <- optim(log(c(1, 1, 0.1)), neg_log_ml, method = "L-BFGS-B",
                lower = log(c(1e-2, 1e-2, 1e-4)), upper = log(c(1e2, 1e2, 1e1)))
theta_hat <- exp(opt_ml$par)
c(sigma_f = round(theta_hat[1], 4), length_scale = round(theta_hat[2], 4),
  sigma_n = round(theta_hat[3], 4), true_noise_sd = 0.15,
  log_marginal_likelihood = round(-opt_ml$value, 3))

#>                 sigma_f            length_scale                 sigma_n 
#>                  1.0247                  1.9561                  0.0968 
#>           true_noise_sd log_marginal_likelihood 
#>                  0.1500                  6.2910

ell_grid <- exp(seq(log(0.1), log(10), length.out = 45))
sn_grid  <- exp(seq(log(0.01), log(1.5), length.out = 45))
Zml <- outer(sn_grid, ell_grid, Vectorize(function(sn, el)
  -neg_log_ml(log(c(theta_hat[1], el, sn)))))

plot_ly(x = ell_grid, y = sn_grid, z = Zml, type = "surface",
        colorscale = "Viridis",
        colorbar = list(title = "log marginal\nlikelihood")) |>
  add_trace(x = theta_hat[2], y = theta_hat[3], z = -opt_ml$value,
            type = "scatter3d", mode = "markers", name = "MLE",
            marker = list(size = 7, color = "red")) |>
  layout(title = "Log marginal likelihood over length-scale and noise level",
         scene = list(xaxis = list(title = "Length-scale", type = "log"),
                      yaxis = list(title = "Noise SD", type = "log"),
                      zaxis = list(title = "log p(y | X, theta)")))

The surface has a clear interior ridge. Toward small length-scale the complexity penalty $-\frac12\log|K_y|$ collapses the likelihood, the model could fit anything, so it explains nothing. Toward large length-scale the data-fit term collapses instead. The maximum sits where the two balance, which is Occam’s razor made numerical.

29 Acquisition functions

Bayesian optimization uses the posterior to decide where to evaluate next, trading exploitation (where $\mu$ is promising) against exploration (where $\sigma$ is large).

Writing $\mu=\mu(x)$, $\sigma=\sigma(x)$, and $f^+$ for the incumbent best observed value, for a maximization problem:

\[ \begin{aligned} \textbf{UCB: }&\quad \alpha_{\text{UCB}}(x)=\mu(x)+\beta\,\sigma(x)\\[2mm] \textbf{PI: }&\quad \alpha_{\text{PI}}(x)=\Phi\!\left(\frac{\mu(x)-f^+-\xi}{\sigma(x)}\right)\\[2mm] \textbf{EI: }&\quad \alpha_{\text{EI}}(x)= \begin{cases} \big(\mu(x)-f^+-\xi\big)\Phi(Z)+\sigma(x)\phi(Z), & \sigma(x)>0\\[1mm] 0, & \sigma(x)=0 \end{cases}, \quad Z=\frac{\mu(x)-f^+-\xi}{\sigma(x)} \end{aligned} \]

Three points that are easy to get wrong.

The arguments are quotients, $\frac{\mu-f^+-\xi}{\sigma}$, a standardized improvement. Writing $\Phi(\mu-f^+\sigma)$ makes the argument dimensionally meaningless.

The $\sigma>0$ branch is strict. At $\sigma=0$ the posterior is a point mass, there is nothing to gain, and $Z$ is undefined; the two branches must not overlap.

Increasing $\xi$ increases exploration, not exploitation. The exploitation term requires the posterior mean to exceed the incumbent by at least $\xi$ before it contributes, so raising $\xi$ discounts known-good regions and shifts weight onto the $\sigma\phi(Z)$ term. A common default is $\xi=0.01$ on a standardized objective.

For minimization, use the lower confidence bound $\mu(x)-\beta\sigma(x)$ and replace $\mu-f^+$ by $f^+-\mu$ throughout. Most R implementations (rBayesianOptimization, ParBayesianOptimization) maximize, so a loss must be negated before being handed to them, a step that is easy to omit and silently inverts the search.

acq_ei <- function(mu, sd, f_best, xi = 0.01) {
  out <- numeric(length(mu))
  pos <- sd > 1e-12                              # strict inequality
  z <- (mu[pos] - f_best - xi) / sd[pos]
  out[pos] <- (mu[pos] - f_best - xi) * pnorm(z) + sd[pos] * dnorm(z)
  pmax(out, 0)
}
acq_pi  <- function(mu, sd, f_best, xi = 0.01) {
  out <- numeric(length(mu)); pos <- sd > 1e-12
  out[pos] <- pnorm((mu[pos] - f_best - xi) / sd[pos]); out
}
acq_ucb <- function(mu, sd, beta = 2) mu + beta * sd
acq_lcb <- function(mu, sd, beta = 2) mu - beta * sd    # for MINIMIZATION

# Check EI against numerical integration of E[max(f - f_best, 0)]
set.seed(113)
chk <- data.frame(mu = c(0.5, 1.2, -0.3), sd = c(0.4, 0.9, 0.2), f_best = 1.0)
chk$EI_closed_form <- acq_ei(chk$mu, chk$sd, chk$f_best, xi = 0)
chk$EI_numeric <- vapply(seq_len(nrow(chk)), function(i)
  integrate(function(v) pmax(v - chk$f_best[i], 0) *
              dnorm(v, chk$mu[i], chk$sd[i]), -Inf, Inf)$value, numeric(1))
chk |> mutate(across(where(is.numeric), \(z) round(z, 6)),
              agree = abs(EI_closed_form - EI_numeric) < 1e-6)

set.seed(115)
X_a <- c(-4, -2.2, 0.5, 2.8)
f_obj <- function(x) sin(x) + 0.3 * cos(3 * x)
y_a <- f_obj(X_a)
post_a <- gp_posterior(X_a, y_a, x_star, sigma_f = 1, ell = 1, sigma_n = 1e-3)
fb <- max(y_a)

acq_df <- data.frame(
  x = rep(x_star, 3),
  value = c(acq_ei(post_a$mean, post_a$sd, fb),
            acq_pi(post_a$mean, post_a$sd, fb),
            acq_ucb(post_a$mean, post_a$sd, beta = 2)),
  acq = rep(c("Expected improvement", "Probability of improvement",
              "Upper confidence bound"), each = length(x_star)))

p_top <- ggplot() +
  geom_ribbon(data = data.frame(x = x_star, lo = post_a$mean - 2*post_a$sd,
                                hi = post_a$mean + 2*post_a$sd),
              aes(x, ymin = lo, ymax = hi), fill = "grey86") +
  geom_line(data = data.frame(x = x_star, y = f_obj(x_star)),
            aes(x, y), color = "grey35", linetype = "dashed") +
  geom_line(data = data.frame(x = x_star, y = post_a$mean),
            aes(x, y), color = "#3B7DD8", linewidth = 1) +
  geom_point(data = data.frame(x = X_a, y = y_a), aes(x, y),
             color = "firebrick", size = 2.6) +
  labs(title = "GP posterior and three acquisition functions",
       subtitle = "Each acquisition is maximized at a different point, reflecting a different exploration-exploitation balance",
       x = NULL, y = "f(x)") + theme_dspa(10)

p_bot <- ggplot(acq_df, aes(x, value)) +
  geom_line(linewidth = 0.8, color = "#7FB069") +
  geom_vline(data = acq_df |> slice_max(value, by = acq),
             aes(xintercept = x), color = "firebrick", linetype = "dashed") +
  facet_wrap(~ acq, ncol = 1, scales = "free_y") +
  labs(x = "x", y = "Acquisition") + theme_dspa(9)

p_top / p_bot + patchwork::plot_layout(heights = c(1.1, 2))

30 Bayesian optimization

bayes_opt <- function(objective, lower, upper, n_init = 5, n_iter = 20,
                      xi = 0.01, seed = 1) {
  set.seed(seed)
  X <- runif(n_init, lower, upper); y <- vapply(X, objective, numeric(1))
  grid <- seq(lower, upper, length.out = 500)
  best_trace <- cummax(y)
  for (t in seq_len(n_iter)) {
    # Re-learn hyperparameters as evidence accumulates
    nll <- function(p) {
      sf <- exp(p[1]); el <- exp(p[2]); sn <- exp(p[3])
      tryCatch({
        K <- se_kernel(X, X, sf, el) + sn^2 * diag(length(X))
        L <- chol(K); a <- backsolve(L, backsolve(L, y, transpose = TRUE))
        0.5*sum(y*a) + sum(log(diag(L))) + 0.5*length(X)*log(2*pi)
      }, error = function(e) 1e10)
    }
    th <- exp(optim(log(c(sd(y) + 1e-6, (upper-lower)/5, 1e-2)), nll,
                    method = "L-BFGS-B",
                    lower = log(c(1e-3, 1e-2, 1e-4)),
                    upper = log(c(1e3, upper-lower, 1)))$par)
    post <- gp_posterior(X, y, grid, th[1], th[2], th[3])
    ei <- acq_ei(post$mean, post$sd, max(y), xi = xi)
    x_next <- grid[which.max(ei)]
    X <- c(X, x_next); y <- c(y, objective(x_next))
    best_trace <- c(best_trace, max(y))
  }
  list(X = X, y = y, best_x = X[which.max(y)], best_y = max(y),
       trace = best_trace)
}

# Maximize a multimodal function on [-5, 5]
obj_bo <- function(x) sin(x) + 0.3 * cos(3*x) + 0.15 * x
truth_bo <- { g <- seq(-5, 5, length.out = 20000); max(obj_bo(g)) }

bo <- bayes_opt(obj_bo, -5, 5, n_init = 4, n_iter = 20, seed = 121)
c(best_found = round(bo$best_y, 6), true_max = round(truth_bo, 6),
  gap = signif(truth_bo - bo$best_y, 3),
  total_evaluations = length(bo$y))

#>        best_found          true_max               gap total_evaluations 
#>        1.49737400        1.49737900        0.00000576       24.00000000

ggplot(data.frame(eval = seq_along(bo$trace), best = bo$trace), aes(eval, best)) +
  geom_hline(yintercept = truth_bo, linetype = "dashed", color = "firebrick") +
  geom_step(linewidth = 0.9, color = "steelblue") +
  geom_point(size = 1.6) +
  annotate("text", x = length(bo$trace), y = truth_bo, vjust = -0.8, hjust = 1,
           size = 3.2, color = "firebrick", label = "global maximum") +
  labs(title = "Bayesian optimization: best value found against evaluation count",
       subtitle = "Each step is one expensive function evaluation. The initial design is random; the rest are chosen by EI",
       x = "Function evaluations", y = "Best value so far") +
  theme_dspa()

# How the acquisition surface evolves across BO iterations
set.seed(123)
grid_bo <- seq(-5, 5, length.out = 200)
X_run <- runif(4, -5, 5); y_run <- obj_bo(X_run)
acq_hist <- matrix(NA_real_, 16, length(grid_bo))
for (t in 1:16) {
  post <- gp_posterior(X_run, y_run, grid_bo, sigma_f = 1, ell = 1, sigma_n = 1e-3)
  ei <- acq_ei(post$mean, post$sd, max(y_run))
  acq_hist[t, ] <- ei / max(ei + 1e-12)
  xn <- grid_bo[which.max(ei)]
  X_run <- c(X_run, xn); y_run <- c(y_run, obj_bo(xn))
}

plot_ly(x = grid_bo, y = 1:16, z = acq_hist, type = "surface",
        colorscale = "Inferno", reversescale = TRUE,
        colorbar = list(title = "Normalized\nEI")) |>
  layout(title = "The expected-improvement surface across Bayesian-optimization iterations",
         scene = list(xaxis = list(title = "x"),
                      yaxis = list(title = "BO iteration"),
                      zaxis = list(title = "Normalized EI")))

Rotate along the iteration axis and watch the ridges migrate. Early on the acquisition is broad and multi-peaked, the surrogate is uncertain nearly everywhere, so exploration dominates. As observations accumulate the peaks narrow and concentrate near the incumbent optimum: exploitation takes over. The acquisition function encodes that transition automatically, without a schedule.

30.1 Regret

GP-UCB regret bound. For a GP with bounded RKHS norm, the algorithm selecting $x_t=\arg\max\big(\mu_{t-1}(x)+\sqrt{\beta_t}\,\sigma_{t-1}(x)\big)$ attains cumulative regret \[R_T=\sum_{t=1}^{T}\big(f(x^\star)-f(x_t)\big)=O^*\!\left(\sqrt{T\,\beta_T\,\gamma_T}\right)\] with high probability, where $\gamma_T$ is the maximum information gain after $T$ rounds (Srinivas et al., 2010).

Since $R_T/T\to0$, the average regret vanishes and the method is no-regret. The rate depends entirely on $\gamma_T$, which is kernel-specific: $O((\log T)^{d+1})$ for the squared-exponential kernel but $O(T^{d(d+1)/(2\nu+d(d+1))}\log T)$ for Matérn, so smoother kernels give faster learning, and the curse of dimensionality enters through $\gamma_T$. This is why Bayesian optimization is a $d\lesssim20$ method.

The bound requires $\beta_t$ to grow logarithmically in $t$. Fixing $\beta=2.576$ (the normal 99% quantile) is a common practical simplification, and it forfeits the guarantee.

30.2 Does it beat random search?

Common misconception: “Bayesian optimization is always better than random search.” It is better per function evaluation when the objective is smooth enough for the surrogate to model, the dimension is modest, and evaluations are expensive. It is worse when evaluations are cheap, the GP fit and acquisition maximization cost $O(n^3)$ and $O(|{\rm grid}|)$ per iteration, which can exceed the objective itself, and when the dimension is high, where $\gamma_T$ grows and the surrogate cannot learn.

Random search is also a stronger baseline than it looks: with $n$ draws, the probability of landing in the top $q$ fraction of the space is $1-(1-q)^n$, so 60 random draws hit the top 5% with probability 0.95 (Bergstra & Bengio, 2012).

The comparison is cheap to run. Run it.

set.seed(131)
budget <- 24

compare_one <- function(seed) {
  set.seed(seed)
  # Random search
  Xr <- runif(budget, -5, 5); best_rand <- cummax(obj_bo(Xr))
  # Grid search on the same budget
  Xg <- seq(-5, 5, length.out = budget); best_grid <- cummax(obj_bo(Xg))
  # Bayesian optimization on the same budget
  b <- bayes_opt(obj_bo, -5, 5, n_init = 4, n_iter = budget - 4, seed = seed)
  data.frame(eval = 1:budget, random = best_rand, grid = best_grid,
             bayes = b$trace[1:budget], rep = seed)
}
cmp <- bind_rows(lapply(seq_len(BO_REPS), compare_one))

summ <- cmp |>
  summarise(across(c(random, grid, bayes), mean), .by = eval) |>
  pivot_longer(-eval, names_to = "method", values_to = "best")

cmp |> filter(eval == budget) |>
  summarise(across(c(random, grid, bayes), list(mean = mean, sd = sd))) |>
  pivot_longer(everything()) |> mutate(value = round(value, 5))

c(true_maximum = round(truth_bo, 5), budget = budget, replications = BO_REPS)

#> true_maximum       budget replications 
#>      1.49738     24.00000     20.00000

ggplot(summ, aes(eval, best, color = method)) +
  geom_hline(yintercept = truth_bo, linetype = "dashed", color = "grey35") +
  geom_line(linewidth = 1) +
  scale_color_manual(values = c(bayes = "#3B7DD8", random = "#D8433B",
                                 grid = "#7FB069"),
                      labels = c("Bayesian optimization", "Grid search",
                                 "Random search")) +
  labs(title = sprintf("Best value found on an equal budget, averaged over %d replications", BO_REPS),
       subtitle = "Dashed: the global maximum. The comparison is meaningless without replication and a common budget",
       x = "Function evaluations", y = "Best value found", color = NULL) +
  theme_dspa()

# The overhead that decides when BO is worth it
t_bo <- system.time(bayes_opt(obj_bo, -5, 5, n_init = 4, n_iter = 20, seed = 141))
t_rand <- system.time({ set.seed(141); max(obj_bo(runif(24, -5, 5))) })
data.frame(
  method = c("Bayesian optimization (24 evals)", "Random search (24 evals)"),
  seconds = round(c(t_bo[["elapsed"]], t_rand[["elapsed"]]), 4),
  note = c("GP refit + acquisition maximization each iteration",
           "objective evaluations only"))

The overhead is the whole story. Bayesian optimization spends real time choosing each point, so it pays off only when a single objective evaluation costs more than that overhead, training a network, running a simulation, performing an experiment. For an objective evaluated in microseconds, random search wins on wall time even when it loses on evaluation count.

# --- Hyperparameter tuning of a CNN (set KERAS_EVAL = TRUE to run) ---------
# The full treatment belongs with the deep learning chapter; what matters here
# is the protocol, not the architecture.
library(keras3)

# Bounds are CONTINUOUS; integer-valued hyperparameters are rounded explicitly
# inside the objective, so the optimizer searches the interval rather than a
# handful of integer points.
objective_cnn <- function(dropout, log_lr, log_units) {
  units <- as.integer(round(2^log_units))
  lr    <- 10^log_lr
  model <- keras_model_sequential() |>
    layer_conv_2d(32, c(3,3), activation = "relu", input_shape = c(28,28,1)) |>
    layer_max_pooling_2d(c(2,2)) |> layer_dropout(dropout) |>
    layer_flatten() |> layer_dense(units, activation = "relu") |>
    layer_dense(10, activation = "softmax")
  model |> compile(loss = "categorical_crossentropy",
                   optimizer = optimizer_adam(learning_rate = lr),
                   metrics = "accuracy")
  h <- model |> fit(x_train, y_train, epochs = 5, batch_size = 128,
                    validation_data = list(x_val, y_val), verbose = 0)
  # The FINAL-epoch validation accuracy, not the maximum over epochs:
  # taking the max selects on the validation set and biases the surrogate.
  tail(h$metrics$val_accuracy, 1)
}
# Tune on (train, val); report once on the untouched test set.

31 PART VI: AUTOMATIC DIFFERENTIATION

Every gradient method above needs $\nabla f$. Three ways to get one.

Approach	Accuracy	Cost for $\nabla f$, $f:\mathbb R^d\to\mathbb R$	Effort
Symbolic	exact	expression swell	manual or CAS
Finite differences	$O(\sqrt{\epsilon_{\text{mach}}})$	$O(d)$ evaluations	trivial
Forward-mode AD	exact	$O(d)$ evaluations	moderate
Reverse-mode AD	exact	$\mathbf{O(1)}$ evaluations	moderate

Reverse-mode AD computes the full gradient at a constant multiple of the cost of one function evaluation, independent of the dimension. That single fact is why deep learning is possible: a network with $10^8$ parameters would need $10^8$ forward passes per gradient by finite differences, and needs about two by reverse mode.

Backpropagation is reverse-mode automatic differentiation. The chain-rule bookkeeping derived in Chapter 6, §6.5 is not a neural-network-specific algorithm, it is the general reverse sweep, specialized to a layered computational graph.

Forward mode propagates derivatives with the computation and costs $O(d)$ for a gradient but $O(1)$ for a directional derivative; reverse mode records the computation and replays it backward, costing $O(1)$ for a gradient and $O(m)$ for a full Jacobian of an $\mathbb R^d\to\mathbb R^m$ map. Use forward mode when outputs outnumber inputs, reverse mode when inputs outnumber outputs, and in optimization, there is one output.

# Forward-mode AD via dual numbers: x + x' * eps, with eps^2 = 0
dual <- function(v, d = 0) structure(list(v = v, d = d), class = "dual")
Ops.dual <- function(e1, e2) {
  a <- if (inherits(e1, "dual")) e1 else dual(e1)
  b <- if (inherits(e2, "dual")) e2 else dual(e2)
  switch(.Generic,
    "+" = dual(a$v + b$v, a$d + b$d),
    "-" = dual(a$v - b$v, a$d - b$d),
    "*" = dual(a$v * b$v, a$d * b$v + a$v * b$d),          # product rule
    "/" = dual(a$v / b$v, (a$d * b$v - a$v * b$d) / b$v^2),# quotient rule
    "^" = dual(a$v^b$v, b$v * a$v^(b$v - 1) * a$d),
    stop("unsupported: ", .Generic))
}
sin.dual <- function(x) dual(sin(x$v), cos(x$v) * x$d)
exp.dual <- function(x) dual(exp(x$v), exp(x$v) * x$d)
log.dual <- function(x) dual(log(x$v), x$d / x$v)

g_test <- function(x) sin(x * x) * exp(x) + log(x)
x0 <- 1.3
truth_deriv <- cos(x0^2) * 2 * x0 * exp(x0) + sin(x0^2) * exp(x0) + 1/x0

data.frame(
  method = c("analytic", "forward-mode AD (dual numbers)", "finite difference"),
  derivative = signif(c(truth_deriv,
                        g_test(dual(x0, 1))$d,
                        numDeriv::grad(g_test, x0)), 12))

# Reverse-mode AD on an explicit computational graph: forward sweep stores
# values, backward sweep accumulates adjoints dL/dv for every node.
reverse_ad <- function(inputs) {
  tape <- list(); n_in <- length(inputs)
  for (i in seq_along(inputs)) tape[[i]] <- list(val = inputs[i], op = "input",
                                                 args = integer(0), grad = 0)
  push <- function(val, op, args) { tape[[length(tape) + 1]] <<-
      list(val = val, op = op, args = args, grad = 0); length(tape) }
  list(
    add = function(i, j) push(tape[[i]]$val + tape[[j]]$val, "add", c(i, j)),
    mul = function(i, j) push(tape[[i]]$val * tape[[j]]$val, "mul", c(i, j)),
    sin = function(i)    push(sin(tape[[i]]$val), "sin", i),
    exp = function(i)    push(exp(tape[[i]]$val), "exp", i),
    backward = function(out) {
      tape[[out]]$grad <- 1                             # seed dL/dL = 1
      for (k in rev(seq_along(tape))) {
        g <- tape[[k]]$grad; a <- tape[[k]]$args
        if (tape[[k]]$op == "add") { tape[[a[1]]]$grad <- tape[[a[1]]]$grad + g
                                     tape[[a[2]]]$grad <- tape[[a[2]]]$grad + g }
        if (tape[[k]]$op == "mul") {
          tape[[a[1]]]$grad <- tape[[a[1]]]$grad + g * tape[[a[2]]]$val
          tape[[a[2]]]$grad <- tape[[a[2]]]$grad + g * tape[[a[1]]]$val }
        if (tape[[k]]$op == "sin") tape[[a]]$grad <- tape[[a]]$grad + g * cos(tape[[a]]$val)
        if (tape[[k]]$op == "exp") tape[[a]]$grad <- tape[[a]]$grad + g * exp(tape[[a]]$val)
      }
      vapply(seq_len(n_in), \(i) tape[[i]]$grad, numeric(1))
    },
    value = function(i) tape[[i]]$val)
}

# f(x, y) = sin(x*y) + exp(x)   ->   df/dx = y*cos(xy) + e^x,  df/dy = x*cos(xy)
ad <- reverse_ad(c(x = 1.1, y = 2.3))
p  <- ad$mul(1, 2); s <- ad$sin(p); e <- ad$exp(1); out <- ad$add(s, e)
g_rev <- ad$backward(out)

x1 <- 1.1; y1 <- 2.3
data.frame(
  partial = c("df/dx", "df/dy"),
  reverse_mode_AD = signif(g_rev, 10),
  analytic = signif(c(y1*cos(x1*y1) + exp(x1), x1*cos(x1*y1)), 10),
  finite_difference = signif(numDeriv::grad(
    function(v) sin(v[1]*v[2]) + exp(v[1]), c(x1, y1)), 10))

One backward sweep produced both partial derivatives. With $d$ inputs it would still be one sweep, that constant-cost property, not the exactness, is what makes reverse mode indispensable.

31.1 Gradient checking

Always verify a hand-coded gradient against finite differences before trusting an optimizer. A wrong gradient does not usually produce an error; it produces slow convergence, or convergence to the wrong point, and looks like a tuning problem.

The check compares the relative difference $\frac{\|g_{\text{analytic}}-g_{\text{numeric}}\|}{\|g_{\text{analytic}}\|+\|g_{\text{numeric}}\|}$ against a tolerance around $10^{-7}$ for central differences.

grad_check <- function(f, grad, x, tol = 1e-6) {
  ga <- grad(x); gn <- numDeriv::grad(f, x)
  rel <- sqrt(sum((ga - gn)^2)) / (sqrt(sum(ga^2)) + sqrt(sum(gn^2)) + 1e-30)
  list(relative_difference = rel, passes = rel < tol,
       analytic = ga, numeric = gn)
}

correct_grad <- function(x) c(-400*x[1]*(x[2]-x[1]^2) - 2*(1-x[1]), 200*(x[2]-x[1]^2))
buggy_grad   <- function(x) c(-400*x[1]*(x[2]-x[1]^2) - 2*(1-x[1]), 200*(x[2]-x[1]))  # typo

data.frame(
  gradient = c("correct", "buggy (one term wrong)"),
  relative_difference = signif(c(grad_check(rosen, correct_grad, c(0.7, 1.4))$relative_difference,
                                 grad_check(rosen, buggy_grad,  c(0.7, 1.4))$relative_difference), 4),
  passes = c(grad_check(rosen, correct_grad, c(0.7, 1.4))$passes,
             grad_check(rosen, buggy_grad,  c(0.7, 1.4))$passes))

32 PART VII: SYNTHESIS

33 Computational complexity summary

$n$ = data points, $d$ = variables, $m$ = L-BFGS memory, $b$ = mini-batch size, $T$ = Bayesian-optimization iterations, $\kappa=L/\mu$.

Method	Per iteration	Memory	Iterations to $\epsilon$	Needs
Gradient descent	$O(nd)$	$O(d)$	$O(\kappa\log\frac1\epsilon)$	$\nabla f$
Nesterov	$O(nd)$	$O(d)$	$O(\sqrt\kappa\log\frac1\epsilon)$	$\nabla f$
Adam / RMSProp	$O(nd)$	$O(d)$	problem-dependent	$\nabla f$
SGD, batch $b$	$O(bd)$	$O(d)$	$O(1/\epsilon)$ in expectation	$\nabla\ell_i$
Conjugate gradient	$O(nd)$	$O(d)$	$\le d$ exactly (quadratic)	$\nabla f$
Newton	$O(nd^2+d^3)$	$O(d^2)$	$O(\log\log\frac1\epsilon)$ locally	$\nabla^2f$
BFGS	$O(nd+d^2)$	$O(d^2)$	superlinear	$\nabla f$
L-BFGS ($m\!\approx\!10$)	$O(nd+md)$	$O(md)$	superlinear	$\nabla f$
ISTA (proximal)	$O(nd)$ + prox	$O(d)$	$O(1/\epsilon)$	$\nabla f$, prox
FISTA	$O(nd)$ + prox	$O(d)$	$O(1/\sqrt\epsilon)$	$\nabla f$, prox
Nelder–Mead	$O(d)$ evals	$O(d^2)$	no guarantee	$f$ only
Simulated annealing	$O(1)$ eval	$O(d)$	no finite guarantee	$f$ only
GP posterior fit	$O(n^3)$ once, $O(n^2)$/pred	$O(n^2)$	—	$f$ only
Bayesian optimization	$O(T^3+\lvert\text{grid}\rvert)$	$O(T^2)$	$R_T=O^*(\sqrt{T\beta_T\gamma_T})$	$f$ only
Simplex (LP)	—	$O(md)$	exponential worst case	linear
Interior point (LP/QP)	$O(d^3)$	$O(d^2)$	$O(\sqrt d\log\frac1\epsilon)$	convex
Finite-difference gradient	$O(d)$ evals	$O(d)$	—	$f$ only
Reverse-mode AD	$\mathbf{O(1)}$ evals	$O(\text{graph})$	—	code for $f$

Five consequences.

Acceleration converts $\kappa$ into $\sqrt\kappa$. At $\kappa=10^4$ that is a hundredfold reduction in iterations, for the cost of storing one extra vector. It is the cheapest improvement available to a first-order method.

The $O(d^2)$ memory bound rules out second-order methods in high dimensions. A dense Hessian at $d=10^8$ needs $4\times10^{16}$ bytes. L-BFGS’s $O(md)$ is affordable; Adam’s $O(d)$ is what deep learning actually uses.

Reverse-mode AD is dimension-free. The full gradient costs a constant multiple of one function evaluation, independent of $d$. That is the single computational fact that makes training $10^8$-parameter models feasible, and backpropagation is its specialization to layered graphs.

Gaussian processes are cubic. The $O(n^3)$ Cholesky caps exact GP inference around $n\approx10^4$, which is why Bayesian optimization is used where evaluations number in the dozens or hundreds, and why sparse and inducing-point approximations exist for anything larger.

Bayesian optimization’s overhead can exceed its objective. The per-iteration $O(T^3)$ GP refit plus acquisition maximization is worth paying only when a single objective evaluation costs more. For a cheap objective, random search wins on wall time even when it loses on evaluation count.

34 Common pitfalls

#	Pitfall	Consequence	Fix
1	Mixing sign conventions for the descent direction	The method ascends	$\nabla f^\top\nu<0$, update $x+\alpha\nu$
2	Requiring $\nabla^2f\succ0$ as a necessary condition	Rejects legitimate minima like $x^4$	Necessary is $\succeq0$; sufficient is $\succ0$
3	“Convex functions have no local minima”	They have them; they are all global	Every local min is a global min
4	Step size chosen by trial and error	Divergence above $2/L$, crawling below	$\alpha<2/L$; use $1/L$ or a line search
5	Reading zig-zag as a defect of gradient descent	It is the $\alpha$-too-large regime on ill-conditioned data	Reduce $\alpha$, precondition, or accelerate
6	Ignoring the condition number	$\kappa$ sets the iteration count	Report $\kappa$; scale features
7	Armijo-only line search with BFGS	$y^\top s>0$ not guaranteed; $H$ loses definiteness	Wolfe curvature condition, or skip the update
8	A non-symmetric quasi-Newton update	Not an inverse Hessian; direction may ascend	Congruence form $VHV^\top+\rho ss^\top$
9	Undamped Newton far from the optimum	Divergence, or convergence to a saddle	Line search plus Hessian regularization
10	Newton at $d=10^6$	$O(d^2)$ memory is infeasible	L-BFGS or a first-order method
11	Constant step size in SGD	Converges to a noise ball, not the optimum	Decay $\alpha_k$; check Robbins–Monro
12	Gradient descent on a non-smooth objective	No gradient at the kink	Proximal or subgradient methods
13	Re-deriving soft-thresholding without naming it	Misses that it is $\operatorname{prox}_{\lambda\\|\cdot\\|_1}$	Proximal framework unifies it
14	Lagrange multipliers for inequality constraints	Misses activity and the sign restriction	KKT with $\mu\ge0$ and complementary slackness
15	Ignoring constraint qualifications	KKT need not be necessary	Check LICQ or Slater
16	Reporting a primal value with no bound	No certificate of optimality	Dual value bounds it; report the gap
17	Redundant constraints presented as binding	Suggests structure that is not there	Non-negativity via bounds, once
18	`solve(K)` on a GP kernel matrix	Numerically singular; explicit inverse	Cholesky of $K+\sigma_n^2I$
19	Omitting the noise term or jitter	The “singular matrix” failure	Always add $\sigma_n^2I$; match code to equation
20	Non-standard kernel bandwidth, undeclared	Length-scales incomparable with every reference	$\exp(-(x-x')^2/2\ell^2)$
21	Acquisition arguments written as products	$\Phi(\mu-f^+\sigma)$ is dimensionally meaningless	The argument is $(\mu-f^+-\xi)/\sigma$
22	Believing larger $\xi$ means more exploitation	It requires a larger margin, so more exploration	Higher $\xi$ explores
23	Handing a loss to a maximizing BO routine	Silently searches for the worst point	Negate, or use LCB for minimization
24	Claiming BO beats random search without measuring	The overhead often dominates for cheap objectives	Equal budget, replicated, timed

35 Practice problems

35.1 Problem 1: Locate the stability boundary

Find the largest step size for which gradient descent converges on a quadratic, and confirm it equals $2/L$.

Solution

set.seed(201)
find_boundary <- function(kappa) {
  q <- quad_factory(kappa)
  alphas <- seq(0.01, 2.5 / q$L, length.out = 300)
  conv <- vapply(alphas, function(a) {
    r <- gradient_descent(q$f, q$grad, c(1, 1), alpha = a, max_iter = 200, tol = 0)
    is.finite(r$f) && r$f < q$f(c(1, 1))
  }, logical(1))
  c(kappa = kappa, L = q$L,
    empirical_max_alpha = max(alphas[conv]), theoretical_2_over_L = 2 / q$L)
}
as.data.frame(do.call(rbind, lapply(c(1, 5, 20, 100), find_boundary))) |>
  mutate(across(everything(), \(z) round(z, 5)),
         relative_error = round(abs(empirical_max_alpha - theoretical_2_over_L) /
                                  theoretical_2_over_L, 4))

The empirical boundary matches $2/L$ to within the grid resolution at every condition number. The threshold depends only on $L$, not on $\mu$ or on the starting point, $\mu$ governs the rate below the threshold, not the threshold itself.

35.2 Problem 2: Verify the acceleration rate

Confirm that Nesterov’s iteration count scales as $\sqrt\kappa$ while gradient descent’s scales as $\kappa$.

Solution

set.seed(203)
iters_to <- function(kappa, target = 1e-8) {
  q <- quad_factory(kappa); a <- 1 / q$L; x0 <- c(1, 1); f0 <- q$f(x0)
  count <- function(fn) {
    r <- fn(q$f, q$grad, x0, alpha = a, max_iter = 20000, tol = 0)
    idx <- which(r$fvals / f0 < target)
    if (length(idx)) idx[1] else NA_integer_
  }
  c(kappa = kappa, gd = count(gradient_descent), nesterov = count(nesterov))
}
p2 <- as.data.frame(do.call(rbind, lapply(c(10, 40, 160, 640), iters_to)))
p2 |> mutate(gd_over_kappa = round(gd / kappa, 3),
             nest_over_sqrt_kappa = round(nesterov / sqrt(kappa), 3))

p2 |> pivot_longer(c(gd, nesterov), names_to = "method", values_to = "iters") |>
  ggplot(aes(kappa, iters, color = method)) +
  geom_line(linewidth = 1) + geom_point(size = 2.4) +
  scale_x_log10() + scale_y_log10() +
  scale_color_manual(values = c(gd = "#D8433B", nesterov = "#3B7DD8"),
                      labels = c("Gradient descent", "Nesterov")) +
  labs(title = "Iterations to reach a fixed accuracy, against condition number",
       subtitle = "On log-log axes the slopes are the exponents: 1 for GD, 1/2 for Nesterov",
       x = "Condition number (log)", y = "Iterations (log)", color = NULL) +
  theme_dspa()

The ratio columns are the diagnostic: gd_over_kappa is roughly constant, and so is nest_over_sqrt_kappa. On log-log axes the two slopes are 1 and $1/2$ — the theoretical exponents.

35.3 Problem 3: Break BFGS by breaking symmetry

Implement the non-symmetric update and show the method fails.

Solution

bfgs_broken <- function(f, grad, x0, max_iter = 200) {
  d <- length(x0); x <- x0; H <- diag(d); g <- grad(x)
  asc <- 0; sym_fail <- 0
  for (i in seq_len(max_iter)) {
    if (sqrt(sum(g^2)) < 1e-10) break
    nu <- as.vector(-H %*% g)
    if (sum(g * nu) >= 0) asc <- asc + 1
    a <- 1; for (bt in 1:40) {
      if (f(x + a*nu) <= f(x) + 1e-4*a*sum(g*nu)) break
      a <- a/2 }
    xn <- x + a*nu; gn <- grad(xn); s <- xn - x; y <- gn - g
    if (sum(s*y) > 1e-10) {
      rho <- 1/sum(s*y)
      V <- diag(d) - rho*(s %*% t(y))
      H <- V %*% H %*% V + rho*(s %*% t(y))     # WRONG: V not transposed, s y'
      if (max(abs(H - t(H))) > 1e-8) sym_fail <- sym_fail + 1
    }
    x <- xn; g <- gn
  }
  list(x = x, f = f(x), ascent_directions = asc, asymmetric_updates = sym_fail,
       H_symmetric = max(abs(H - t(H))) < 1e-8)
}
bad <- bfgs_broken(rosen, rosen_g, c(-1.2, 1))
good <- bfgs(rosen, rosen_g, c(-1.2, 1))

data.frame(
  version = c("correct update", "non-symmetric update"),
  final_f = signif(c(good$f, bad$f), 4),
  distance_to_optimum = signif(c(sqrt(sum((good$x - c(1,1))^2)),
                                 sqrt(sum((bad$x - c(1,1))^2))), 4),
  H_stayed_symmetric = c(isTRUE(all.equal(good$H, t(good$H), tolerance = 1e-10)),
                         bad$H_symmetric),
  ascent_directions_encountered = c(0, bad$ascent_directions))

The broken update destroys symmetry within a few iterations, after which $-Hg$ is no longer guaranteed to be a descent direction and the method stalls far from the optimum. The congruence form $VHV^\top$ is what preserves symmetry, and losing it breaks the algorithm rather than merely slowing it.

35.4 Problem 4: Verify KKT on a problem with a mix of active and inactive constraints

Solution

# min (x1-3)^2 + (x2-2)^2  s.t.  x1+x2 <= 4  (active),  x1 <= 5 (inactive), x >= 0
f4  <- function(x) (x[1]-3)^2 + (x[2]-2)^2
gf4 <- function(x) c(2*(x[1]-3), 2*(x[2]-2))
G   <- rbind(c(1, 1), c(1, 0))              # constraint gradients (rows)
b4  <- c(4, 5)

res4 <- nloptr::nloptr(
  x0 = c(1, 1),
  eval_f = function(x) list(objective = f4(x), gradient = gf4(x)),
  eval_g_ineq = function(x) list(constraints = as.vector(G %*% x) - b4,
                                 jacobian = G),
  lb = c(0, 0),
  opts = list(algorithm = "NLOPT_LD_SLSQP", xtol_rel = 1e-12, maxeval = 1000))
xs <- res4$solution
g_vals <- as.vector(G %*% xs) - b4
active <- abs(g_vals) < 1e-7

# Solve for the multipliers of the ACTIVE constraints only
Ga <- G[active, , drop = FALSE]
mu_a <- as.vector(qr.solve(t(Ga), -gf4(xs)))
mu_all <- numeric(2); mu_all[active] <- mu_a

data.frame(
  constraint = c("x1 + x2 <= 4", "x1 <= 5"),
  g_value = round(g_vals, 8),
  active = active,
  mu = round(mu_all, 6),
  complementary_slackness = round(mu_all * g_vals, 10),
  dual_feasible = mu_all >= -1e-10)

c(solution = round(xs, 6),
  stationarity_residual = signif(sqrt(sum((gf4(xs) + as.vector(t(G) %*% mu_all))^2)), 3))

#>             solution1             solution2 stationarity_residual 
#>              2.50e+00              1.50e+00              4.44e-16

The inactive constraint receives multiplier zero automatically, that is complementary slackness doing its work. Only the active constraint carries a positive multiplier, and stationarity balances the objective gradient against it alone.

35.5 Problem 5: Soft-thresholding is a proximal operator

Verify numerically that $\operatorname{prox}_{t\lambda\|\cdot\|_1}$ equals the soft-threshold.

Solution

prox_l1_numeric <- function(v, t, lambda)
  vapply(v, function(vi) optimize(function(u)
    lambda*abs(u) + (u - vi)^2/(2*t), c(vi - 10, vi + 10), tol = 1e-12)$minimum,
    numeric(1))

vs <- seq(-3, 3, length.out = 25); t5 <- 1; lam5 <- 0.8
p5 <- data.frame(v = vs,
                 closed_form = soft_threshold(vs, t5 * lam5),
                 numeric = prox_l1_numeric(vs, t5, lam5))
c(max_abs_difference = signif(max(abs(p5$closed_form - p5$numeric)), 3))

#> max_abs_difference 
#>           2.16e-08

ggplot(p5, aes(v)) +
  geom_hline(yintercept = 0, color = "grey70") +
  geom_line(aes(y = v, color = "identity"), linetype = "dotted") +
  geom_line(aes(y = closed_form, color = "soft threshold"), linewidth = 1.2) +
  geom_point(aes(y = numeric, color = "numeric prox"), size = 2) +
  scale_color_manual(values = c(identity = "grey55",
                                 `soft threshold` = "#3B7DD8",
                                 `numeric prox` = "#D8433B")) +
  coord_fixed() +
  labs(title = "The proximal operator of the L1 norm IS soft-thresholding",
       subtitle = "Points: numerical minimization of the prox objective. Line: the closed form from Chapter 11",
       x = "v", y = expression(prox[t*lambda*"||.||"[1]](v)), color = NULL) +
  theme_dspa()

They agree to twelve decimals. This is the same operator that Chapter 11 derives from the LASSO subgradient condition, one object, reached from two directions, which is why every sparse estimator reduces to a thresholding step.

35.6 Problem 6: Jitter and conditioning

Quantify how the GP kernel’s condition number depends on the number of points, the length-scale, and the jitter.

Solution

set.seed(205)
p6 <- expand.grid(n = c(10, 25, 50, 100), jitter = c(0, 1e-10, 1e-8, 1e-6, 1e-4))
p6$condition <- vapply(seq_len(nrow(p6)), function(i) {
  X <- seq(-3, 3, length.out = p6$n[i])
  K <- se_kernel(X, X, 1, 1) + p6$jitter[i] * diag(p6$n[i])
  kappa(K)
}, numeric(1))
p6$cholesky_ok <- vapply(seq_len(nrow(p6)), function(i) {
  X <- seq(-3, 3, length.out = p6$n[i])
  K <- se_kernel(X, X, 1, 1) + p6$jitter[i] * diag(p6$n[i])
  !inherits(try(chol(K), silent = TRUE), "try-error")
}, logical(1))
p6 |> mutate(condition = signif(condition, 3)) |> arrange(n, jitter)

ggplot(p6, aes(n, condition, color = factor(jitter))) +
  geom_hline(yintercept = 1/.Machine$double.eps, linetype = "dashed",
             color = "firebrick") +
  geom_line(linewidth = 0.9) + geom_point(size = 2) +
  scale_y_log10() +
  scale_color_viridis_d(option = "plasma", end = 0.9, name = "Jitter") +
  labs(title = "Kernel-matrix conditioning against sample size and jitter",
       subtitle = "Dashed: 1/eps, beyond which double precision cannot represent the inverse",
       x = "Number of points", y = "Condition number (log scale)") +
  theme_dspa()

Without jitter the condition number crosses $1/\epsilon_{\text{mach}}$ by about 25 points, the failure is a property of the smooth kernel, not of having “too much” data. A jitter of $10^{-8}$ keeps the matrix factorizable at every size tested while perturbing the model negligibly.

35.7 Problem 7: When does Bayesian optimization stop winning?

Compare BO against random search as the dimension grows.

Solution

set.seed(207)
obj_nd <- function(x) -sum((x - 0.3)^2) + 0.25 * sum(sin(6 * x))

bo_nd <- function(d, budget = 40, seed = 1) {
  set.seed(seed)
  n0 <- min(10, budget %/% 2)
  X <- matrix(runif(n0 * d), n0, d); y <- apply(X, 1, obj_nd)
  cand <- matrix(runif(600 * d), 600, d)
  for (t in seq_len(budget - n0)) {
    Kx <- exp(-as.matrix(dist(X))^2 / 2) + 1e-6 * diag(nrow(X))
    L <- chol(Kx); a <- backsolve(L, backsolve(L, y, transpose = TRUE))
    Ks <- exp(-as.matrix(dist(rbind(X, cand)))[1:nrow(X),
                                               (nrow(X)+1):(nrow(X)+600)]^2 / 2)
    mu <- as.vector(t(Ks) %*% a)
    v  <- backsolve(L, Ks, transpose = TRUE)
    sd <- sqrt(pmax(1 - colSums(v^2), 1e-12))
    xi <- cand[which.max(acq_ei(mu, sd, max(y))), , drop = FALSE]
    X <- rbind(X, xi); y <- c(y, obj_nd(as.vector(xi)))
  }
  max(y)
}
rand_nd <- function(d, budget = 40, seed = 1) {
  set.seed(seed + 5000)
  max(apply(matrix(runif(budget * d), budget, d), 1, obj_nd))
}

p7 <- do.call(rbind, lapply(c(1, 2, 5, 10, 20), function(d) {
  b <- vapply(1:8, \(s) bo_nd(d, seed = s), numeric(1))
  r <- vapply(1:8, \(s) rand_nd(d, seed = s), numeric(1))
  data.frame(dimension = d, bayes = mean(b), random = mean(r),
             advantage = mean(b) - mean(r))
}))
p7 |> mutate(across(-dimension, \(z) round(z, 4)))

p7 |> pivot_longer(c(bayes, random), names_to = "method", values_to = "best") |>
  ggplot(aes(dimension, best, color = method)) +
  geom_line(linewidth = 1) + geom_point(size = 2.4) +
  scale_x_log10(breaks = p7$dimension) +
  scale_color_manual(values = c(bayes = "#3B7DD8", random = "#D8433B"),
                      labels = c("Bayesian optimization", "Random search")) +
  labs(title = "Bayesian optimization's advantage shrinks as the dimension grows",
       subtitle = "Equal budget of 40 evaluations, averaged over 8 replications",
       x = "Dimension (log scale)", y = "Best value found", color = NULL) +
  theme_dspa()

The advantage column shrinks toward zero. The reason is in the regret bound: the maximum information gain $\gamma_T$ grows with dimension, so a fixed budget buys progressively less information about the surrogate. Below roughly $d\approx20$ the surrogate helps; beyond it, random search is competitive and far cheaper.

35.8 Problem 8: Cost of a gradient, three ways

Measure how finite differences, forward-mode AD, and reverse-mode AD scale with dimension.

Solution

set.seed(209)
f_scal <- function(x) sum(sin(x) * exp(x / 10)) + sum(x^2) / 2
g_anal <- function(x) cos(x) * exp(x/10) + sin(x) * exp(x/10)/10 + x

p8 <- do.call(rbind, lapply(c(5, 20, 80, 200), function(d) {
  x <- runif(d, -1, 1)
  t_fd <- median(replicate(5, system.time(numDeriv::grad(f_scal, x))[["elapsed"]]))
  t_an <- median(replicate(5, system.time(g_anal(x))[["elapsed"]]))
  data.frame(dimension = d,
             finite_diff_sec = t_fd,
             analytic_sec = max(t_an, 1e-6),
             fd_function_evals = 2 * d,          # central differences
             reverse_ad_evals = 1)
}))
p8 |> mutate(across(where(is.numeric), \(z) signif(z, 3)),
             ratio = signif(finite_diff_sec / analytic_sec, 3))

ggplot(p8, aes(dimension)) +
  geom_line(aes(y = fd_function_evals, color = "Finite differences"), linewidth = 1) +
  geom_line(aes(y = reverse_ad_evals, color = "Reverse-mode AD"), linewidth = 1) +
  geom_point(aes(y = fd_function_evals, color = "Finite differences"), size = 2.4) +
  geom_point(aes(y = reverse_ad_evals, color = "Reverse-mode AD"), size = 2.4) +
  scale_x_log10(breaks = p8$dimension) + scale_y_log10() +
  scale_color_manual(values = c(`Finite differences` = "#D8433B",
                                 `Reverse-mode AD` = "#3B7DD8")) +
  labs(title = "Function evaluations required for one full gradient",
       subtitle = "Finite differences scale as O(d); reverse-mode AD is constant. This is why deep learning is possible",
       x = "Dimension (log scale)", y = "Function evaluations (log scale)",
       color = NULL) +
  theme_dspa()

Finite differences need $2d$ evaluations for central differences and carry truncation error around $\epsilon_{\text{mach}}^{2/3}$. Reverse-mode AD needs a constant number and is exact. At $d=10^8$ that is the difference between impossible and routine.

36 Checkpoint

Gradient descent on your problem oscillates and the loss grows. Name the likeliest cause and the diagnostic.
Your BFGS implementation returns a direction that increases the objective. What is the most likely bug?
Your solver reports an optimal value of 47.3. How do you know it is optimal?
A GP fit fails with “system is computationally singular.” What went wrong?
Bayesian optimization did no better than random search on your problem. Give three possible reasons.
You need the gradient of a loss with $10^7$ parameters. What are your options?

Answers

The step size exceeds the stability threshold $2/L$. For an $L$-smooth objective, gradient descent is guaranteed to decrease only when $0<\alpha<2/L$; above that the error contraction factor $\max_i|1-\alpha\lambda_i|$ exceeds 1 and the iterates diverge geometrically. The diagnostic is to estimate $L$ as the largest Hessian eigenvalue (or the largest curvature you can bound) and check $\alpha$ against $2/L$; the practical fix is a backtracking line search, which finds an acceptable step without knowing $L$. If the loss oscillates without growing, you are in the $1/L<\alpha<2/L$ band on an ill-conditioned problem, reduce $\alpha$, rescale the features to lower $\kappa$, or switch to an accelerated or adaptive method.
The update lost symmetry or positive definiteness. BFGS produces a descent direction $-Hg$ only while $H\succ0$, and $H$ stays symmetric only because the update has the congruence form $VHV^\top+\rho ss^\top$, the right factor must be the transpose of the left, and the rank-one term must be $s s^\top$, not $s y^\top$. The other common cause is a violated curvature condition $y_k^\top s_k>0$, which is what the Wolfe line search exists to guarantee; with an Armijo-only search it can fail, and the standard remedy is to skip the update that iteration. Check all.equal(H, t(H)) and min(eigen(H)$values) > 0 after every update.
Only a dual bound certifies optimality. A primal value is just the best point found; on a non-convex problem it may be arbitrarily far from $p^\star$. For a convex problem, weak duality gives $g(\mu,\lambda)\le p^\star$ for any dual-feasible $(\mu,\lambda)$, so a dual value of 47.3 alongside a primal value of 47.3 closes the duality gap and proves optimality, and under Slater’s condition strong duality guarantees the gap can be closed. Failing that, verify the KKT conditions numerically: stationarity residual near zero, primal feasibility, $\mu\ge0$, and $\mu_ig_i=0$. On a non-convex problem, report the best value found, the number of distinct optima reached from multiple starts, and the spread.
The noise term or jitter is missing. A smooth kernel evaluated at closely spaced inputs produces a $K$ whose eigenvalues decay geometrically to machine precision, it is numerically singular by construction, not because there are too many points. Adding $\sigma_n^2I$ (a genuine observation-noise term) or a jitter $\epsilon I\approx10^{-8}I$ shifts every eigenvalue up and restores conditioning. Two further fixes belong with it: never form the inverse, take the Cholesky factor of $K+\sigma_n^2I$ and solve triangular systems, which is more accurate and yields $\log|K_y|$ for free, and check that the code matches the equation, since the posterior formula contains $\sigma_n^2I$ whether or not the implementation does.
Dimension, cheapness, or roughness. Dimension: the regret bound depends on the maximum information gain $\gamma_T$, which grows with $d$, so beyond roughly $d\approx20$ a fixed budget buys too little information for the surrogate to help. Cheapness: BO’s per-iteration $O(T^3)$ GP refit plus acquisition maximization can exceed the cost of the objective itself, so random search wins on wall time even while losing on evaluation count. Roughness: if the objective is discontinuous, noisy, or has a length-scale far shorter than the budget can resolve, the GP surrogate is a poor model and its suggestions are no better than chance. Also check the mundane failure — whether a loss was handed to a maximizing routine without negation, which searches for the worst point.
Reverse-mode automatic differentiation, and essentially nothing else. Finite differences would need about $2\times10^7$ function evaluations per gradient and carry truncation error; forward-mode AD is exact but also $O(d)$. Reverse mode computes the entire gradient at a constant multiple of one function evaluation, independent of $d$, that is precisely why training models at this scale is feasible, and backpropagation is this algorithm specialized to layered computational graphs. Memory is the trade-off, since the forward sweep must store intermediate values; gradient checkpointing recomputes some of them to trade time for space. Whatever you use, check the gradient against finite differences on a small instance first, a wrong gradient rarely errors, it just converges slowly or to the wrong place.

37 Summary

The problem and its structure

Optimality conditions come in pairs: $\nabla f=0$ with $\nabla^2f\succeq0$ is necessary; with $\nabla^2f\succ0$ it is sufficient.
For convex problems every local minimum is a global minimum, they exist, they are simply all equally good.
$L$-smoothness and $\mu$-strong convexity bound the curvature above and below, and their ratio $\kappa=L/\mu$ governs every first-order rate.

First-order methods

Gradient descent converges only for $\alpha<2/L$; the optimal constant step is $\frac{2}{L+\mu}$, and zig-zagging is the $\alpha$-too-large regime, not a defect of the method.
Rates: $O(1/k)$ for smooth convex, linear $(1-1/\kappa)^k$ under strong convexity.
Nesterov acceleration gives $O(1/k^2)$ and $(1-1/\sqrt\kappa)^k$, and it is optimal among first-order methods.
Adaptive methods are diagonal preconditioners; Adam’s bias correction matters most in the first few steps.
SGD trades per-iteration cost for variance; a constant step converges to a noise ball, and the Robbins–Monro conditions are what force decay.

Second-order methods

Newton converges quadratically but only locally; away from the optimum it needs damping, regularization, or a trust region.
BFGS preserves symmetry through the congruence form $VHV^\top+\rho ss^\top$ and positive definiteness through the curvature condition $y^\top s>0$, which the Wolfe line search supplies.
L-BFGS reduces $O(d^2)$ memory to $O(md)$, which is what makes quasi-Newton usable at scale.

Constrained optimization

KKT: stationarity, primal feasibility, dual feasibility $\mu\ge0$, and complementary slackness $\mu_ig_i=0$, necessary under a constraint qualification, and sufficient for convex problems under Slater.
Complementary slackness is why SVMs are sparse.
The dual is always concave and always lower-bounds the primal; the gap closes for convex problems under Slater.
Proximal methods handle non-smooth terms, and the soft-threshold is $\operatorname{prox}_{\lambda\|\cdot\|_1}$, the same operator Chapter 11 derives from the LASSO subgradient.

Derivative-free and Bayesian

Gaussian processes require the noise term and a Cholesky factorization; the singular-matrix failure is a missing jitter, not excess data.
The marginal likelihood is a single $n$-variate Gaussian density, and its $-\frac12\log|K_y|$ term implements Occam’s razor.
Acquisition arguments are standardized quotients $(\mu-f^+-\xi)/\sigma$, and larger $\xi$ explores more.
GP-UCB has vanishing average regret, at a rate set by the information gain $\gamma_T$, which is why BO is a low-dimensional method.
Bayesian optimization beats random search per evaluation, not per second. Measure both.

Differentiation

Reverse-mode AD computes a full gradient at constant cost in the dimension, and backpropagation is its specialization to layered graphs.
Check hand-coded gradients against finite differences. A wrong gradient looks like a tuning problem.

Where these threads continue

Thread	Continues in
Backpropagation as reverse-mode AD; Adam in practice	Deep learning
The SVM dual and kernel methods	Black-box methods
Proximal operators and regularization paths	Feature selection

38 Chapter roadmap

Chapter 1: Foundations. R toolchain, reproducibility, simulation.
Chapter 2: Data quality and visual analytics.
Chapter 3: Linear algebra and regression. SVD, conditioning, the hat matrix.
Chapter 4: Dimensionality reduction.
Chapter 5: Supervised classification.
Chapter 6: Black-box methods. Backpropagation, the SVM dual, kernels.
Chapter 7: Text mining and association rules.
Chapter 8: Unsupervised clustering.
Chapter 9: Model assessment and validation.
Chapter 10: Formats, streams, performance.
Chapter 11: Feature selection. LASSO, soft-thresholding, knockoffs.
Chapter 12: Longitudinal analysis and forecasting.
Deep learning. Convolutional and sequence models, automatic differentiation at scale, hyperparameter search.

39 Session information

sessionInfo()

#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#> 
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8 
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] knitr_1.51          Rsolnp_1.16         numDeriv_2016.8-1.1
#> [4] plotly_4.12.1       patchwork_1.3.0     tidyr_1.3.1        
#> [7] dplyr_1.1.4         ggplot2_4.0.1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6       jsonlite_1.8.9     compiler_4.3.3     tidyselect_1.2.1  
#>  [5] parallel_4.3.3     jquerylib_0.1.4    scales_1.4.0       yaml_2.3.10       
#>  [9] fastmap_1.2.0      R6_2.6.1           generics_0.1.3     htmlwidgets_1.6.4 
#> [13] tibble_3.2.1       bslib_0.9.0        pillar_1.10.1      RColorBrewer_1.1-3
#> [17] rlang_1.1.5        cachem_1.1.0       xfun_0.52          sass_0.4.9        
#> [21] S7_0.2.1           otel_0.2.0         truncnorm_1.0-9    viridisLite_0.4.2 
#> [25] cli_3.6.3          withr_3.0.2        magrittr_2.0.3     crosstalk_1.2.1   
#> [29] digest_0.6.37      grid_4.3.3         rstudioapi_0.18.0  lifecycle_1.0.5   
#> [33] vctrs_0.6.5        evaluate_1.0.3     glue_1.8.0         data.table_1.16.4 
#> [37] farver_2.1.2       httr_1.4.7         rmarkdown_2.31     purrr_1.0.2       
#> [41] tools_4.3.3        pkgconfig_2.0.3    htmltools_0.5.8.1

Method	Cost per iteration	Memory	Rate (strongly convex)
Full gradient descent	\(O(nd)\)	\(O(d)\)	\((1-1/\kappa)^k\)
Nesterov	\(O(nd)\)	\(O(d)\)	\((1-1/\sqrt\kappa)^k\)
Adam / RMSProp	\(O(nd)\)	\(O(d)\)	problem-dependent
SGD, batch \(b\)	\(O(bd)\)	\(O(d)\)	\(O(1/k)\) in expectation
Newton	\(O(nd^2+d^3)\)	\(O(d^2)\)	quadratic, locally
BFGS	\(O(nd+d^2)\)	\(O(d^2)\)	superlinear
L-BFGS, memory \(m\)	\(O(nd+md)\)	\(O(md)\)	superlinear

Method	Memory	Per-iteration	Feasible at \(d=10^6\)?
Newton	\(O(d^2)\)	\(O(nd^2+d^3)\)	No
BFGS	\(O(d^2)\)	\(O(nd+d^2)\)	No (\(8\times10^{12}\) bytes)
L-BFGS (\(m=10\))	\(O(md)\)	\(O(nd+md)\)	Yes (\(8\times10^7\) bytes)
Gradient descent	\(O(d)\)	\(O(nd)\)	Yes

Method	Per iteration	Memory	Iterations to \(\epsilon\)	Needs
Gradient descent	\(O(nd)\)	\(O(d)\)	\(O(\kappa\log\frac1\epsilon)\)	\(\nabla f\)
Nesterov	\(O(nd)\)	\(O(d)\)	\(O(\sqrt\kappa\log\frac1\epsilon)\)	\(\nabla f\)
Adam / RMSProp	\(O(nd)\)	\(O(d)\)	problem-dependent	\(\nabla f\)
SGD, batch \(b\)	\(O(bd)\)	\(O(d)\)	\(O(1/\epsilon)\) in expectation	\(\nabla\ell_i\)
Conjugate gradient	\(O(nd)\)	\(O(d)\)	\(\le d\) exactly (quadratic)	\(\nabla f\)
Newton	\(O(nd^2+d^3)\)	\(O(d^2)\)	\(O(\log\log\frac1\epsilon)\) locally	\(\nabla^2f\)
BFGS	\(O(nd+d^2)\)	\(O(d^2)\)	superlinear	\(\nabla f\)
L-BFGS (\(m\!\approx\!10\))	\(O(nd+md)\)	\(O(md)\)	superlinear	\(\nabla f\)
ISTA (proximal)	\(O(nd)\) + prox	\(O(d)\)	\(O(1/\epsilon)\)	\(\nabla f\), prox
FISTA	\(O(nd)\) + prox	\(O(d)\)	\(O(1/\sqrt\epsilon)\)	\(\nabla f\), prox
Nelder–Mead	\(O(d)\) evals	\(O(d^2)\)	no guarantee	\(f\) only
Simulated annealing	\(O(1)\) eval	\(O(d)\)	no finite guarantee	\(f\) only
GP posterior fit	\(O(n^3)\) once, \(O(n^2)\)/pred	\(O(n^2)\)	—	\(f\) only
Bayesian optimization	\(O(T^3+\lvert\text{grid}\rvert)\)	\(O(T^2)\)	\(R_T=O^*(\sqrt{T\beta_T\gamma_T})\)	\(f\) only
Simplex (LP)	—	\(O(md)\)	exponential worst case	linear
Interior point (LP/QP)	\(O(d^3)\)	\(O(d^2)\)	\(O(\sqrt d\log\frac1\epsilon)\)	convex
Finite-difference gradient	\(O(d)\) evals	\(O(d)\)	—	\(f\) only
Reverse-mode AD	\(\mathbf{O(1)}\) evals	\(O(\text{graph})\)	—	code for \(f\)

\(\nabla^2f(x^\star)\)	Eigenvalues	Critical point
Positive definite	all \(>0\)	strict local minimum
Negative definite	all \(<0\)	strict local maximum
Indefinite	mixed signs	saddle point
Semi-definite, singular	some \(=0\)	inconclusive at second order

#	Pitfall	Consequence	Fix
1	Mixing sign conventions for the descent direction	The method ascends	\(\nabla f^\top\nu<0\), update \(x+\alpha\nu\)
2	Requiring \(\nabla^2f\succ0\) as a necessary condition	Rejects legitimate minima like \(x^4\)	Necessary is \(\succeq0\); sufficient is \(\succ0\)
3	“Convex functions have no local minima”	They have them; they are all global	Every local min is a global min
4	Step size chosen by trial and error	Divergence above \(2/L\), crawling below	\(\alpha<2/L\); use \(1/L\) or a line search
5	Reading zig-zag as a defect of gradient descent	It is the \(\alpha\)-too-large regime on ill-conditioned data	Reduce \(\alpha\), precondition, or accelerate
6	Ignoring the condition number	\(\kappa\) sets the iteration count	Report \(\kappa\); scale features
7	Armijo-only line search with BFGS	\(y^\top s>0\) not guaranteed; \(H\) loses definiteness	Wolfe curvature condition, or skip the update
8	A non-symmetric quasi-Newton update	Not an inverse Hessian; direction may ascend	Congruence form \(VHV^\top+\rho ss^\top\)
9	Undamped Newton far from the optimum	Divergence, or convergence to a saddle	Line search plus Hessian regularization
10	Newton at \(d=10^6\)	\(O(d^2)\) memory is infeasible	L-BFGS or a first-order method
11	Constant step size in SGD	Converges to a noise ball, not the optimum	Decay \(\alpha_k\); check Robbins–Monro
12	Gradient descent on a non-smooth objective	No gradient at the kink	Proximal or subgradient methods
13	Re-deriving soft-thresholding without naming it	Misses that it is \(\operatorname{prox}_{\lambda\\|\cdot\\|_1}\)	Proximal framework unifies it
14	Lagrange multipliers for inequality constraints	Misses activity and the sign restriction	KKT with \(\mu\ge0\) and complementary slackness
15	Ignoring constraint qualifications	KKT need not be necessary	Check LICQ or Slater
16	Reporting a primal value with no bound	No certificate of optimality	Dual value bounds it; report the gap
17	Redundant constraints presented as binding	Suggests structure that is not there	Non-negativity via bounds, once
18	`solve(K)` on a GP kernel matrix	Numerically singular; explicit inverse	Cholesky of \(K+\sigma_n^2I\)
19	Omitting the noise term or jitter	The “singular matrix” failure	Always add \(\sigma_n^2I\); match code to equation
20	Non-standard kernel bandwidth, undeclared	Length-scales incomparable with every reference	\(\exp(-(x-x')^2/2\ell^2)\)
21	Acquisition arguments written as products	\(\Phi(\mu-f^+\sigma)\) is dimensionally meaningless	The argument is \((\mu-f^+-\xi)/\sigma\)
22	Believing larger \(\xi\) means more exploitation	It requires a larger margin, so more exploration	Higher \(\xi\) explores
23	Handing a loss to a maximizing BO routine	Silently searches for the worst point	Negate, or use LCB for minimization
24	Claiming BO beats random search without measuring	The overhead often dominates for cheap objectives	Equal budget, replicated, timed

DSPA3: Data Science and Predictive Analytics (UMich HS650)

Chapter 13: Function Optimization

SOCR/MCAIM (Ivo Dinov)

August 29, 2026

How this chapter uses graphics

1 Learning objectives

2 PART I: THE OPTIMIZATION PROBLEM

3 Formulation

4 Optimality conditions

5 Convexity

6 Smoothness, strong convexity, and the condition number

7 PART II, FIRST-ORDER METHODS

8 Gradient descent

8.1 Stability: how large may the step be?

8.2 Convergence rates

9 Line search

10 Momentum and acceleration

11 Adaptive methods

12 Stochastic gradient descent

13 Cost per iteration

14 PART III: SECOND-ORDER AND QUASI-NEWTON METHODS

15 Newton’s method

16 BFGS

16.1 L-BFGS

17 Non-convex landscapes

18 PART IV: CONSTRAINED OPTIMIZATION

19 Equality constraints and Lagrange multipliers

20 The KKT conditions

21 Duality

22 Linear and quadratic programming

23 Practice Optimization Examples

23.1 Example 1: Healthcare Manufacturer Product Optimization

23.1.1 Unconstrained Optimization

23.1.2 Constrained Optimization

23.1.3 Manual Solution

23.1.4 R-based Automated solution

23.2 Example 2: Optimization of the Booth’s function

23.3 Example 3: Extrema of the bivariate Goldstein-Price Function

23.4 Example 4: Bivariate Oscillatory Function

23.5 Nonlinear Constraint Optimization Problem

24 Proximal methods

24.1 ADMM

25 PART V: DERIVATIVE-FREE AND BAYESIAN OPTIMIZATION

26 When gradients are unavailable

27 Gaussian processes

27.1 Sampling from the prior

27.2 Posterior inference, done stably

28 Marginal likelihood and hyperparameter learning

29 Acquisition functions

30 Bayesian optimization

30.1 Regret

30.2 Does it beat random search?

31 PART VI: AUTOMATIC DIFFERENTIATION

31.1 Gradient checking

32 PART VII: SYNTHESIS

33 Computational complexity summary

34 Common pitfalls

35 Practice problems

35.1 Problem 1: Locate the stability boundary

35.2 Problem 2: Verify the acceleration rate

35.3 Problem 3: Break BFGS by breaking symmetry

35.4 Problem 4: Verify KKT on a problem with a mix of active and inactive constraints

35.5 Problem 5: Soft-thresholding is a proximal operator

35.6 Problem 6: Jitter and conditioning

35.7 Problem 7: When does Bayesian optimization stop winning?

35.8 Problem 8: Cost of a gradient, three ways

36 Checkpoint

37 Summary

38 Chapter roadmap

39 Session information