SOCR ≫ DSPA ≫ DSPA3 Topics ≫ 

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:

  1. State first- and second-order optimality conditions, distinguishing necessary from sufficient.
  2. Define \(L\)-smoothness and \(\mu\)-strong convexity, and explain how the condition number \(\kappa=L/\mu\) governs convergence.
  3. Derive the gradient-descent stability condition \(\eta<2/L\) and the convergence rates \(O(1/k)\) and \((1-\mu/L)^k\).
  4. Implement backtracking line search and explain the Wolfe conditions.
  5. Explain why Nesterov acceleration achieves \(O(1/k^2)\) and why that is optimal for first-order methods.
  6. Implement BFGS from the correct recursion and verify the curvature condition.
  7. State and apply the KKT conditions, and use complementary slackness to explain SVM sparsity.
  8. Construct the Lagrange dual, and state weak duality, Slater’s condition, and the duality gap.
  9. Compute proximal operators, and recognize soft-thresholding as \(\operatorname{prox}_{\lambda\|\cdot\|_1}\).
  10. 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)

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 =