Numerical Optimization and Gradient Descent

Almost every optimization in this course has been discrete: pick the best subset, the shortest path, the cheapest spanning tree, the longest common subsequence. The search space is finite (if astronomically large), and the tools are combinatorial — exchange arguments, dynamic programming, network flow. This lesson changes the setting. The variable $x$ ranges over $R^{n}$ , a continuum, and the objective $f : R^{n} \to R$ is a smooth, differentiable function we want to minimize. There is no finite set to enumerate and no subproblem table to fill. Instead we use calculus: the derivative tells us which way is downhill, and we keep stepping that way. This method underlies curve fitting, logistic regression, and the training of every neural network, and it rests on the gradient.

The continuous problem

We are given a differentiable $f : R^{n} \to R$ and want a point $x^{⋆}$ where $f$ is smallest:

x^{⋆} = ar g x \in R^{n} min f (x) .

In one dimension, calculus says a minimum of a smooth $f$ occurs where the slope vanishes, $f^{'} (x) = 0$ , and the curve bends upward, $f^{''} (x) \geq 0$ . In $n$ dimensions the slope becomes a vector, the gradient

\nabla f (x) = (\frac{\partial f}{\partial x _{1}}, \dots, \frac{\partial f}{\partial x _{n}}),

and a minimum requires $\nabla f (x) = 0$ . Solving $\nabla f (x) = 0$ in closed form is possible only for the simplest $f$ (a quadratic gives a linear system). For everything else we need an iterative method that starts somewhere and improves, step by step. Every such method starts from the geometry of the gradient.

The gradient points uphill

Fix a point $x$ and ask: among all unit directions $u$ , which one increases $f$ the fastest? The first-order Taylor expansion answers it. Moving a small distance $ε$ along $u$ ,

f (x + ε u) \approx f (x) + ε ⟨ \nabla f (x), u ⟩,

so the rate of change in direction $u$ is the inner product $⟨ \nabla f (x), u ⟩$ . By Cauchy–Schwarz this is largest when $u$ aligns with $\nabla f (x)$ and most negative when $u$ points the opposite way. So the gradient is the direction of steepest ascent, and its negation $- \nabla f (x)$ is the direction of steepest descent. The gradient is also perpendicular to the level curve ${x : f (x) = c}$ through the point: along the level curve $f$ does not change, so the steepest-change direction must be orthogonal to it.

The gradient

\nabla f

at a point is perpendicular to the level curve through it and points toward higher values; the descent direction is

- \nabla f

Gradient descent

The algorithm follows directly. Start at any $x_{0}$ . At each step, move a little in the steepest-descent direction:

x_{t + 1} = x_{t} - η \nabla f (x_{t}),

where $η > 0$ is the learning rate (or step size). Each step strictly decreases $f$ for small enough $η$ , because the first-order change is $f (x_{t + 1}) - f (x_{t}) \approx - η ∥ \nabla f (x_{t}) ∥^{2} \leq 0$ . We stop when the gradient is nearly zero — the slope has flattened, and we are at (or near) a stationary point.

Algorithm 1:

\textsc{Gradient-Descent}(f, x_0, \eta)

— step against the gradient until flat

1
$x \gets x_0$
2
while $\lVert \nabla f(x) \rVert > \varepsilon$ do
not yet flat
3
$g \gets \nabla f(x)$
steepest-ascent direction
4
$x \gets x - \eta \cdot g$
step downhill
5
return $x$

gradient_descent.pypython

from dataclasses import dataclass, field
from typing import Callable, Sequence

# A point and a gradient are both real vectors; the gradient maps one to the
# other. f is the scalar objective evaluated at a point.
Vector = list[float]
Gradient = Callable[[Sequence[float]], Sequence[float]]
Objective = Callable[[Sequence[float]], float]

def _euclidean_norm(vector: Sequence[float]) -> float:
  """
    Length of `vector` under the usual Euclidean (L2) norm.\n
  """
  return sum(component * component for component in vector) ** 0.5

@dataclass
class DescentResult:
  """
    The outcome of a descent run: the final point, the objective value\n
    there, how many steps were taken, and whether the gradient flattened\n
    below tolerance (converged) or the step budget ran out.\n
  """
  point: Vector
  value: float
  iterations: int
  converged: bool
  trajectory: list[Vector] = field(default_factory=lambda: [])

def gradient_descent(
  objective: Objective,
  gradient: Gradient,
  start: Sequence[float],
  learning_rate: float = 0.1,
  tolerance: float = 1e-8,
  max_iterations: int = 10_000,
  record_trajectory: bool = False,
) -> DescentResult:
  """
    Minimize `objective` by stepping against `gradient` from `start`.\n
    Each iteration moves x <- x - learning_rate * gradient(x); the loop\n
    stops once the gradient norm falls below `tolerance` (converged) or the\n
    iteration budget is exhausted. Set `record_trajectory` to keep every\n
    visited point for inspection or plotting.\n
  """
  current: Vector = list(start)
  trajectory: list[Vector] = [list(current)] if record_trajectory else []

  for iteration in range(max_iterations):

    # stop early once the gradient flattens below tolerance.
    slope: Sequence[float] = gradient(current)
    if _euclidean_norm(slope) <= tolerance:
      return DescentResult(
        point=current,
        value=objective(current),
        iterations=iteration,
        converged=True,
        trajectory=trajectory,
      )

    # step downhill: x <- x - eta * grad(x), component by component.
    current = [
      position - learning_rate * partial
      for position, partial in zip(current, slope)
    ]
    if record_trajectory:
      trajectory.append(list(current))

  return DescentResult(
    point=current,
    value=objective(current),
    iterations=max_iterations,
    converged=_euclidean_norm(gradient(current)) <= tolerance,
    trajectory=trajectory,
  )

The method has one free parameter, and it matters.

If $η$ is too small, the iterates creep toward the minimum, needing huge numbers of steps. If $η$ is too large, a step can jump clear over the valley floor and land higher than it started; successive overshoots then bounce between the walls, oscillating, or spiral outward and diverge. Much of first-order optimization comes down to choosing (or adapting) $η$ .

Worked example (three learning rates on a parabola). Minimize $f (x) = x^{2}$ , for which $f^{'} (x) = 2 x$ and the update is $x_{t + 1} = x_{t} - η \cdot 2 x_{t} = (1 - 2 η) x_{t}$ — each step multiplies the iterate by the fixed contraction factor $r = 1 - 2 η$ . Starting from $x_{0} = 1$ , the behaviour is decided entirely by $∣ r ∣$ :

$η = 0.1$ gives $r = 0.8$ : the iterates $1, 0.8, 0.64, 0.512, \dots$ shrink steadily to $0$ — healthy convergence.
$η = 0.9$ gives $r = - 0.8$ : $1, - 0.8, 0.64, - 0.512, \dots$ still converges, but oscillates in sign, overshooting the minimum each step before creeping in.
$η = 1.1$ gives $r = - 1.2$ : $1, - 1.2, 1.44, - 1.728, \dots$ grows without bound — the step is so large it lands farther from the minimum each time, and the method diverges.

The threshold is exactly $η < 1$ here (so that $∣1 - 2 η ∣ < 1$ ); more generally, for an $L$ -smooth objective the safe range is $η < 2/ L$ , and $η = 1/ L$ is the textbook choice. The condition $r = 1 - 2 η$ is a one-line preview of why the strong-convexity rate below is geometric: a fixed multiplicative shrink per step.

A well-chosen

η

(left) descends steadily to the minimum; too large an

η

(right) overshoots and oscillates across the valley.

Convexity, and why local is global

Gradient descent finds a point where the gradient vanishes — a stationary point. For a general nonconvex $f$ that could be a local minimum, a saddle, or a poor local optimum far from the best value. Convexity eliminates this problem.

A convex function is a single bowl with no false bottoms, so local information determines the global optimum.

So on a convex objective, gradient descent's walk downhill until flat cannot get trapped: the only flat spot is the bottom.

A convex bowl: the graph lies above every tangent and below every chord, so the single stationary point is the global minimum.

A convergence guarantee

How fast does gradient descent reach the bottom? The answer depends on two regularity constants. We say $f$ is $L$ -smooth if its gradient does not change too abruptly, $∥ \nabla f (x) - \nabla f (y)∥ \leq L ∥ x - y ∥$ — equivalently $f$ is bounded above by a quadratic of curvature $L$ . We say $f$ is $μ$ -strongly convex if it is bounded below by a quadratic of curvature $μ > 0$ (a bowl that curves at least as much as $\frac{μ}{2} ∥ x ∥^{2}$ ). With these two constants the rate is classical.

To reach error $ϵ$ , the convex case needs $O (1/ ϵ)$ steps, while the strongly convex case needs only $O (κ log (1/ ϵ))$ — exponentially fewer. The condition number $κ = L / μ$ is the elongation of the bowl: a round bowl ( $κ \approx 1$ ) converges in a few steps; a long thin valley ( $κ$ large) makes gradient descent zig-zag slowly down the narrow axis. Contrast this with the discrete optimization of the rest of the course, where rate is measured in the input size $n$ (see asymptotic analysis); here the iteration count depends on the conditioning of $f$ and the target accuracy $ϵ$ , not on a combinatorial size.

Stochastic gradient descent

In machine learning the objective is an average over data, $f (x) = \frac{1}{m} \sum_{i = 1}^{m} f_{i} (x)$ , where $f_{i}$ is the loss on the $i$ -th of $m$ training examples. The exact gradient $\nabla f$ sums $m$ per-example gradients — one full sweep over the dataset per step, which is ruinous when $m$ is in the millions. Stochastic gradient descent (SGD) replaces the full gradient with a cheap, noisy estimate from a small random mini-batch $B$ :

x_{t + 1} = x_{t} - η \nabla f_{B} (x_{t}), \nabla f_{B} (x) = \frac{1}{∣ B ∣} i \in B \sum \nabla f_{i} (x) .

Each step costs $∣ B ∣ ≪ m$ gradients, so SGD takes vastly more, vastly cheaper steps. The mini-batch gradient is an unbiased estimate of the true gradient, so on average it points downhill; the variance is noise that jitters the path.

stochastic_gradient_descent.pypython

import random
from dataclasses import dataclass
from typing import Callable, Optional, Sequence

# Each training example contributes its own gradient at the current point.
Vector = list[float]
PerExampleGradient = Callable[[Sequence[float], int], Sequence[float]]
StepSchedule = Callable[[int], float]

@dataclass
class SGDResult:
  """
    The outcome of an SGD run: the final parameter vector, the number of\n
    mini-batch updates applied, and the number of full passes (epochs).\n
  """
  point: Vector
  updates: int
  epochs: int

def constant_schedule(rate: float) -> StepSchedule:
  """
    A step schedule that returns `rate` regardless of the update count.\n
  """
  return lambda _update: rate

def decaying_schedule(initial: float, decay: float = 1.0) -> StepSchedule:
  """
    A 1/t-style schedule: step `initial / (1 + decay * update)` shrinks the\n
    learning rate as updates accumulate, letting the iterate settle.\n
  """
  return lambda update: initial / (1.0 + decay * update)

def stochastic_gradient_descent(
  example_gradient: PerExampleGradient,
  start: Sequence[float],
  num_examples: int,
  step_schedule: StepSchedule,
  batch_size: int = 1,
  epochs: int = 100,
  generator: Optional[random.Random] = None,
) -> SGDResult:
  """
    Minimize a mean-over-data objective from `start` by mini-batch SGD.\n
    `example_gradient(point, index)` returns the gradient of the loss on the\n
    example at `index`; each update averages it over a random batch of\n
    `batch_size` examples and steps against it with the current scheduled\n
    rate. One epoch shuffles and sweeps all `num_examples` examples once.\n
    Pass a seeded `generator` for reproducible shuffles.\n
  """
  # fall back to a fresh, unseeded generator when none is supplied.
  randomness: random.Random = (
    generator if generator is not None else random.Random()
  )
  current: Vector = list(start)
  dimension: int = len(current)
  update_count: int = 0

  for _epoch in range(epochs):

    # reshuffle the example order at the start of every epoch.
    order: list[int] = list(range(num_examples))
    randomness.shuffle(order)

    # walk the shuffled examples in contiguous mini-batches.
    for batch_start in range(0, num_examples, batch_size):
      batch: list[int] = order[batch_start : batch_start + batch_size]

      # sum the per-example gradients over the batch.
      averaged: Vector = [0.0 for _ in range(dimension)]
      for example_index in batch:
        partials: Sequence[float] = example_gradient(current, example_index)
        for component in range(dimension):
          averaged[component] += partials[component]

      # rescale the sum to a mean (unbiased estimate of the full gradient).
      scale: float = 1.0 / len(batch)
      averaged = [partial * scale for partial in averaged]

      # step against the averaged gradient at the current scheduled rate.
      learning_rate: float = step_schedule(update_count)
      current = [
        position - learning_rate * partial
        for position, partial in zip(current, averaged)
      ]
      update_count += 1

  return SGDResult(point=current, updates=update_count, epochs=epochs)

Second order: Newton's method

Gradient descent uses only the slope. Newton's method also uses the curvature — the matrix of second derivatives, the Hessian $\nabla^{2} f (x)$ — to take a better-scaled step. Approximate $f$ near $x_{t}$ by its second-order Taylor expansion (a quadratic) and jump straight to that quadratic's minimum. Setting the gradient of the approximation to zero gives

x_{t + 1} = x_{t} - [\nabla^{2} f (x_{t})]^{- 1} \nabla f (x_{t}) .

Where gradient descent crawls down a long thin valley, Newton's method rescales by the curvature and heads straight for the bottom — and on a quadratic it reaches the exact minimum in one step. Near a minimum its convergence is quadratic: the number of correct digits roughly doubles each iteration.

Algorithm 2:

\textsc{Newton}(f, x_0)

— second-order steps; quadratic local convergence

1
$x \gets x_0$
2
while $\lVert \nabla f(x) \rVert > \varepsilon$ do
3
$g \gets \nabla f(x)$
gradient
4
$H \gets \nabla^2 f(x)$
Hessian (curvature)
5
$x \gets x - H^{-1} g$
solve $H\,\Delta = g$ , then step
6
return $x$

newtons_method.pypython

from dataclasses import dataclass
from typing import Callable, Sequence

Vector = list[float]
Matrix = list[list[float]]
Gradient = Callable[[Sequence[float]], Sequence[float]]
Hessian = Callable[[Sequence[float]], Sequence[Sequence[float]]]
Objective = Callable[[Sequence[float]], float]

def _euclidean_norm(vector: Sequence[float]) -> float:
  """
    Length of `vector` under the Euclidean (L2) norm.\n
  """
  return sum(component * component for component in vector) ** 0.5

def solve_linear_system(matrix: Sequence[Sequence[float]],
                        rhs: Sequence[float]) -> Vector:
  """
    Solve `matrix * x = rhs` by Gaussian elimination with partial pivoting.\n
    Used to take a Newton step without inverting the Hessian. Raises\n
    ValueError if the matrix is singular (no unique solution).\n
  """
  size: int = len(rhs)

  # work on an augmented copy [matrix | rhs] so the inputs stay untouched.
  augmented: Matrix = [
    [float(matrix[row][column]) for column in range(size)] + [float(rhs[row])]
    for row in range(size)
  ]

  for pivot in range(size):

    # partial pivoting: swap in the largest-magnitude entry in this column.
    best_row: int = max(
      range(pivot, size), key=lambda row: abs(augmented[row][pivot])
    )
    if abs(augmented[best_row][pivot]) < 1e-15:
      raise ValueError("singular matrix: Newton step is undefined")
    augmented[pivot], augmented[best_row] = augmented[best_row], augmented[pivot]

    # eliminate this column from every row below the pivot.
    for row in range(pivot + 1, size):
      factor: float = augmented[row][pivot] / augmented[pivot][pivot]
      for column in range(pivot, size + 1):
        augmented[row][column] -= factor * augmented[pivot][column]

  # back-substitute from the last variable up to the first.
  solution: Vector = [0.0 for _ in range(size)]
  for row in range(size - 1, -1, -1):

    # subtract the already-solved variables, then divide by the pivot.
    accumulated: float = augmented[row][size]
    for column in range(row + 1, size):
      accumulated -= augmented[row][column] * solution[column]
    solution[row] = accumulated / augmented[row][row]
  return solution

@dataclass
class NewtonResult:
  """
    The outcome of a Newton run: the final point, the objective value there,\n
    the number of steps taken, and whether the gradient flattened below\n
    tolerance.\n
  """
  point: Vector
  value: float
  iterations: int
  converged: bool

def newtons_method(
  objective: Objective,
  gradient: Gradient,
  hessian: Hessian,
  start: Sequence[float],
  tolerance: float = 1e-10,
  max_iterations: int = 100,
) -> NewtonResult:
  """
    Minimize `objective` by Newton's method from `start`.\n
    Each step solves hessian(x) * delta = gradient(x) and moves\n
    x <- x - delta, stopping once the gradient norm drops below `tolerance`.\n
    Best started near a minimizer where the Hessian is positive definite.\n
  """
  current: Vector = list(start)

  for iteration in range(max_iterations):

    # stop once the gradient flattens below tolerance.
    slope: Sequence[float] = gradient(current)
    if _euclidean_norm(slope) <= tolerance:
      return NewtonResult(
        point=current,
        value=objective(current),
        iterations=iteration,
        converged=True,
      )

    # take a curvature-corrected step: solve H * delta = grad, then x -= delta.
    curvature: Sequence[Sequence[float]] = hessian(current)
    step: Vector = solve_linear_system(curvature, slope)
    current = [position - delta for position, delta in zip(current, step)]

  return NewtonResult(
    point=current,
    value=objective(current),
    iterations=max_iterations,
    converged=_euclidean_norm(gradient(current)) <= tolerance,
  )

The cost per step is high. Each step forms the $n \times n$ Hessian and solves a linear system with it, costing $Θ (n^{3})$ (or $Θ (n^{2})$ memory just to store $H$ ) — prohibitive for the high-dimensional problems of machine learning. And the quadratic convergence is only local: started far away, Newton's step can overshoot wildly. Practical solvers therefore use quasi-Newton methods (BFGS), which approximate $H^{- 1}$ from successive gradients, and add a line search or trust region for global safety.

One dimension: Newton–Raphson and bisection

The cleanest case of Newton's method is root-finding in one variable: solve $g (x) = 0$ . Replacing the Hessian with the scalar derivative, the update becomes the classic Newton–Raphson iteration

x_{t + 1} = x_{t} - \frac{g ( x _{t} )}{g ^{'} ( x _{t} )} .

Geometrically, follow the tangent at $x_{t}$ down to where it crosses the axis, and take that crossing as the next guess. (Minimizing $f$ is finding a root of $g = f^{'}$ ; this is the same iteration applied to the derivative.) It converges quadratically when it converges — computing $a$ as a root of $g (x) = x^{2} - a$ gives the famous $x_{t + 1} = \frac{1}{2} (x_{t} + a / x_{t})$ , doubling correct digits each step. Sqrt(x) asks for this iteration verbatim.

Newton–Raphson: the tangent at

x_{t}

meets the axis at

x_{t + 1} = x_{t} - g (x_{t}) / g^{'} (x_{t})

, leaping toward the root.

newton_raphson.pypython

from dataclasses import dataclass
from typing import Callable

ScalarFunction = Callable[[float], float]

@dataclass
class RootResult:
  """
    The outcome of a root search: the located root, the residual |g(root)|,\n
    the number of iterations, and whether the residual met the tolerance.\n
  """
  root: float
  residual: float
  iterations: int
  converged: bool

def newton_raphson(
  function: ScalarFunction,
  derivative: ScalarFunction,
  start: float,
  tolerance: float = 1e-12,
  max_iterations: int = 100,
) -> RootResult:
  """
    Find a root of `function` from `start` using its `derivative`.\n
    Iterates x <- x - function(x) / derivative(x) until |function(x)| falls\n
    below `tolerance`. Raises ValueError if the derivative vanishes, which\n
    would send the tangent horizontal and the next guess to infinity.\n
  """
  current: float = start

  for iteration in range(max_iterations):

    # stop once the residual |g(x)| is within tolerance.
    value: float = function(current)
    if abs(value) <= tolerance:
      return RootResult(
        root=current,
        residual=abs(value),
        iterations=iteration,
        converged=True,
      )

    # a horizontal tangent has no axis crossing to step toward.
    slope: float = derivative(current)
    if abs(slope) < 1e-15:
      raise ValueError("derivative vanished: Newton-Raphson step is undefined")

    # tangent crossing: x_{t+1} = x_t - g(x_t) / g'(x_t).
    current = current - value / slope

  return RootResult(
    root=current,
    residual=abs(function(current)),
    iterations=max_iterations,
    converged=abs(function(current)) <= tolerance,
  )

def square_root(value: float, tolerance: float = 1e-12) -> float:
  """
    Newton-Raphson square root of a non-negative `value`, the root of\n
    x^2 - value with the classic update x <- (x + value/x) / 2.\n
  """
  if value < 0:
    raise ValueError("square root of a negative number is undefined")
  if value == 0:
    return 0.0

  # stop on a relative residual so small values converge as tightly as large
  # ones (an absolute residual would leave coarse error when value << 1).
  guess: float = value if value >= 1 else 1.0
  while abs(guess * guess - value) > tolerance * value:
    guess = 0.5 * (guess + value / guess)
  return guess

Newton–Raphson is fast but fragile: a near-zero derivative sends the tangent nearly horizontal and flings the next guess far away, and a bad start may not converge at all. When robustness matters more than speed, use bisection. If $g$ is continuous and changes sign across an interval $[a, b]$ — that is, $g (a)$ and $g (b)$ have opposite signs — the intermediate value theorem guarantees a root inside. Test the midpoint, keep the half-interval that still brackets the sign change, and repeat. The bracket halves every step.

Bisection maintains a sign-changing bracket

[a, b]

; testing the midpoint

m

keeps the half that still brackets the root, halving the interval each step.

bisection.pypython

from dataclasses import dataclass
from typing import Callable

ScalarFunction = Callable[[float], float]
MonotonePredicate = Callable[[int], bool]

@dataclass
class BisectionResult:
  """
    The outcome of a bisection search: the located root, the residual\n
    |g(root)|, the final bracket width, and the number of midpoint tests.\n
  """
  root: float
  residual: float
  width: float
  iterations: int

def bisection(
  function: ScalarFunction,
  low: float,
  high: float,
  tolerance: float = 1e-12,
  max_iterations: int = 1000,
) -> BisectionResult:
  """
    Find a root of continuous `function` bracketed by `low` and `high`.\n
    Requires function(low) and function(high) to have opposite signs (a\n
    bracket); raises ValueError otherwise. Halves the bracket until its\n
    width drops below `tolerance` or the iteration budget is spent.\n
  """
  if low > high:
    low, high = high, low

  low_value: float = function(low)
  high_value: float = function(high)

  # an exact hit at an endpoint is already a root.
  if low_value == 0.0:
    return BisectionResult(low, 0.0, high - low, 0)
  if high_value == 0.0:
    return BisectionResult(high, 0.0, high - low, 0)

  # same-sign endpoints give no guaranteed crossing to bracket.
  if low_value * high_value > 0.0:
    raise ValueError("endpoints do not bracket a sign change")

  iteration: int = 0
  midpoint: float = 0.5 * (low + high)
  while iteration < max_iterations and (high - low) > tolerance:

    # test the midpoint; an exact zero there ends the search.
    midpoint = 0.5 * (low + high)
    midpoint_value: float = function(midpoint)
    if midpoint_value == 0.0:
      return BisectionResult(midpoint, 0.0, high - low, iteration + 1)

    # keep the half whose endpoints still straddle zero.
    if low_value * midpoint_value < 0.0:
      high = midpoint
    else:
      low, low_value = midpoint, midpoint_value
    iteration += 1

  return BisectionResult(
    root=midpoint,
    residual=abs(function(midpoint)),
    width=high - low,
    iterations=iteration,
  )

def binary_search_answer(
  predicate: MonotonePredicate,
  low: int,
  high: int,
) -> int:
  """
    Smallest integer in [low, high] satisfying a monotone `predicate`.\n
    `predicate` must be False on a (possibly empty) prefix and True\n
    thereafter; this brackets the threshold the way bisection brackets a\n
    root. Returns `high + 1` if the predicate is never satisfied. This is\n
    the discrete answer-search behind problems like Koko Eating Bananas.\n
  """
  result: int = high + 1

  # binary-search the threshold: on a True hit, record it and look lower.
  while low <= high:
    middle: int = low + (high - low) // 2
    if predicate(middle):
      result = middle
      high = middle - 1
    else:
      low = middle + 1
  return result

Bisection converges linearly — one bit of accuracy per step, $O (log \frac{b - a}{ϵ})$ steps to width $ϵ$ — slower than Newton's quadratic doubling, but it cannot fail on a continuous sign-changing function. It is, in disguise, the monotone binary search behind the discrete answer-search problems Koko Eating Bananas and Minimize Max Distance to Gas Station: a feasibility predicate is monotone in the answer, so bracket the threshold and halve. The robust strategy in practice is hybrid — Newton steps when they stay inside the bracket, a bisection step as a safe fallback when they do not.¹

Modern optimizers and neural-network training

The optimizers that actually train modern models are refinements of plain gradient descent.

Momentum and adaptive rates. The zig-zag down a long thin valley (large condition number $κ$ ) is the central weakness of vanilla gradient descent. Momentum damps it by accumulating a velocity, $v_{t + 1} = β v_{t} + \nabla f (x_{t})$ , $x_{t + 1} = x_{t} - η v_{t + 1}$ , so consistent directions build speed while oscillating ones cancel; Nesterov's accelerated gradient sharpens this to the optimal $O (1/ t^{2})$ rate for smooth convex $f$ , a provable improvement over the $O (1/ t)$ derived above.² Adaptive methods — AdaGrad, RMSProp, and especially Adam — give each coordinate its own effective step size from a running estimate of gradient magnitudes, which is why Adam is the default optimizer for deep networks.³ All of them are still $x_{t + 1} = x_{t} - (step)$ ; only the step changes.

Backpropagation is the chain rule. Training a neural network minimizes a loss $f (x)$ over millions of parameters $x$ , and gradient descent needs $\nabla f$ . The gradient is computed by backpropagation — the multivariate chain rule applied in reverse across the network's computation graph — which yields all partial derivatives in a single backward pass costing the same order as one forward evaluation. This is reverse-mode automatic differentiation, and it is what makes gradient descent tractable in millions of dimensions where forming the Hessian for Newton's method (the $Θ (n^{3})$ cost above) is infeasible.⁴

Non-convex reality. The clean local is global guarantee holds only for convex $f$ ; neural-network losses are highly non-convex. Several empirical facts explain why training still works: in very high dimension most critical points are saddles, not bad local minima, and the stochastic noise of SGD helps the iterate slip off them; and the many minima that do exist tend to have similar loss values, so finding a good one suffices. Understanding why first-order methods succeed on non-convex deep networks remains an active research question, but the working answer — SGD plus momentum plus a well-tuned schedule — trains essentially every model in use today.

Takeaways

Continuous optimization minimizes a differentiable $f$ over $R^{n}$ with calculus, not enumeration. The gradient $\nabla f$ is the steepest-ascent direction and is orthogonal to level sets; its negation points downhill.
Gradient descent iterates $x_{t + 1} = x_{t} - η \nabla f (x_{t})$ . The learning rate $η$ is the one critical parameter: too small is slow, too large overshoots and diverges.
Convexity makes every local minimum global, so descend until flat finds the true optimum. For convex $L$ -smooth $f$ the rate is $O (1/ t)$ ; under $μ$ -strong convexity it is geometric, $O (κ log \frac{1}{ϵ})$ steps with $κ = L / μ$ .
Stochastic gradient descent estimates $\nabla f$ from a random mini-batch — far cheaper, noisier steps; the noise/speed trade-off and step decay make it the standard method for neural-network training.
Newton's method uses the Hessian for quadratic local convergence but costs $Θ (n^{3})$ per step; in 1-D it is Newton–Raphson root-finding, with bisection the robust linear-time bracketing fallback.¹

Skiena, § — Numerical Problems / Optimization, and CLRS, Ch. — Newton's method for root-finding: gradient descent, Newton–Raphson, and bisection, and the trade-off between fast-but-fragile and slow-but-robust convergence. ↩ ↩²
Y. Nesterov, A method of solving a convex programming problem with convergence rate $O (1/ k^{2})$ , Soviet Mathematics Doklady 27, 1983; and Nesterov, Introductory Lectures on Convex Optimization, Springer, 2004. ↩
D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, ICLR 2015; J. Duchi, E. Hazan, Y. Singer, Adaptive subgradient methods, JMLR 12, 2011 (AdaGrad). ↩
D. E. Rumelhart, G. E. Hinton, R. J. Williams, Learning representations by back-propagating errors, Nature 323, 1986 — backpropagation as reverse-mode differentiation for training neural networks. ↩

The continuous problem

The gradient points uphill

Gradient descent

Convexity, and why local is global

A convergence guarantee

Stochastic gradient descent

Second order: Newton's method

One dimension: Newton–Raphson and bisection

Modern optimizers and neural-network training

Takeaways

Footnotes