Unconstrained Optimization

🎯 One Master Equation, Five Ways To Fill It In

There are dozens of optimization algorithms with intimidating names — conjugate gradient, BFGS, Adam, trust region. Underneath, every one of them is doing the same thing: standing somewhere, deciding which way to go, deciding how far to go, and stepping. Everything else is bookkeeping.

The universal skeleton

Written as code, every method on this page is the same loop with three blanks to fill in. Learn this shape once and each algorithm below becomes a variation you can read at a glance:

x = x0;  init HISTORY
for k = 0,1,2,...:
    g = grad(x)                       # the information you have
    if |g| < tol: break               # flat ground — you're done

    d = DIRECTION(g, HISTORY)         # ── SLOT A : which way
    a = LENGTH(f, x, d, g)            # ── SLOT B : how far
    x = x + a*d                       # ── THE MASTER EQUATION

    HISTORY = UPDATE(...)             # ── SLOT C : what to remember

A note on scope

Everything here is unconstrained: minimise \(f(x)\) over all of \(\mathbb{R}^n\), with no side conditions. The target is characterised by \(\nabla f = 0\) — find flat ground. Add constraints and the optimum typically sits on a boundary where \(\nabla f \neq 0\), the condition becomes KKT rather than a vanishing gradient, and you need a different toolkit. That's a separate page — though most constrained methods end up calling the machinery below as their inner solver.

📏 Slot B Is Where Methods Actually Differ

Direction rules vary a lot in appearance and surprisingly little in consequence. The step length \(\alpha_k\) is the real fork in the road — and for any given method you are in exactly one of three situations. You are either forced to compute the step, forced to fix it, or free to choose.

The two outer columns carry no ambiguity — there is no decision to make. But they are determined for different reasons, and the difference matters:

Free does not mean safe

When you fix \(\eta\) yourself, there is still a hard limit you must respect. Linearising \(x_{k+1}=x_k-\eta\nabla f(x_k)\) near a minimiser gives \(e_{k+1}=(I-\eta H)e_k\), so the iteration converges only if the spectral radius stays below one:

Above that bound it diverges — guaranteed, every time. The catch is that \(L\) is a global property of a function you only ever see one point at a time, so a fixed \(\eta\) is fundamentally a guess. And nothing warns you: no error, no exception, just silently worse numbers. On a landscape whose curvature varies by two orders of magnitude, no single \(\eta\) is simultaneously stable in the steep parts and useful in the flat ones.

ARMIJO_BACKTRACK(f, x, d, g, alpha0=1.0, c1=1e-4, rho=0.5):
    alpha = alpha0
    while f(x + alpha*d) > f(x) + c1*alpha*(g @ d):   # ◄── INNER LOOP
        alpha = rho * alpha                              #     halve and retry
    return alpha

# Wolfe adds a second test — |grad(x+a*d) @ d| <= c2*|g @ d| —
# which is the version CG and BFGS specifically require.

🧬 The Five Families

Master equation, then pseudocode, in the same three slots every time. Watch slot B change — and watch the loop count change with it.

A:  d = -g                                    # batch GD
    # v = beta*v - eta*g ;  d = v/eta           # momentum
    # s = rho*s + (1-rho)*g**2                  # RMSprop / Adam
    #   d = -g/(sqrt(s)+eps)
B:  a = eta                                   # FIXED — or line-search, your call
C:  keep v and/or s                           # O(n) state, or nothing at all

1 LOOP One gradient per iteration, \(O(n)\) memory, and exactly one number you must supply. Two sub-variants worth naming, because they behave very differently:

• Momentum / Nesterov carry velocity, so the recurrence becomes second-order: \(r^2-(1+\beta-\eta\lambda)r+\beta=0\). When the gradient is persistent the velocity saturates at \(v\to-\eta g/(1-\beta)\) — a tenfold amplification of displacement at \(\beta=0.9\).
• Adam / RMSprop divide by accumulated gradient magnitude, so in the persistent-gradient limit \(\Delta x \to -\eta\,\mathrm{sign}(g)\). The step size becomes \(\approx\eta\) independent of \(\|g\|\) — a cheap diagonal stand-in for a line search, estimated from data already collected rather than solved for.

A:  beta = (g @ g) / (g_old @ g_old)          # Fletcher-Reeves
    # beta = (g @ (g-g_old)) / (g_old @ g_old)  # Polak-Ribiere
    beta = max(beta, 0)                       # PR+ safeguard
    if k % n == 0: beta = 0                   # restart every n iterations
    d = -g + beta*d_old
B:  a = WOLFE(f, x, d, g)                     # COMPUTED — required, not optional
C:  g_old = g ;  d_old = d

2 LOOPS \(O(n)\) memory, and — this is the trade — zero hyperparameters. Nothing to tune, because nothing was left to you.

A:  d = -B @ g
B:  a = WOLFE(f, x, d, g)                     # COMPUTED — required
C:  s = a*d ;  y = grad(x_new) - g
    if y @ s > 0:                             # curvature condition
        r = 1/(y @ s)
        B = (I - r*outer(s,y)) @ B @ (I - r*outer(y,s)) + r*outer(s,s)

2 LOOPS Memory is \(O(n^2)\) — L-BFGS drops that to \(O(mn)\) by storing only the last \(m\) pairs \((s,y)\) instead of the whole matrix. That is why L-BFGS, not anything in family 1, is the default workhorse in scipy.optimize.minimize for deterministic smooth problems.

# ── damped Newton ──
A:  H = hess(x)
    d = solve(H + tau*I, -g)                  # tau > 0 if H not positive-definite
B:  a = ARMIJO(f, x, d, g)                    # COMPUTED  (a → 1 near the solution)
C:  —                                         # H rebuilt from scratch each step

# ── trust region: the skeleton BREAKS here ──
A+B: p = SOLVE_SUBPROBLEM(g, hess(x), Delta)   # direction AND length, jointly
     rho = (f(x) - f(x+p)) / -(g@p + 0.5*p@H@p)
     x = x + p  if rho > 0.1  else  x        # the step can be REJECTED
C:   Delta = 0.25*Delta          if rho < 0.25
     Delta = min(2*Delta, Dmax)  if rho > 0.75 and |p| == Delta

2 LOOPS Uses second derivatives, so \(O(n^2)\) memory and an \(O(n^3)\) linear solve per iteration. Pure Newton with \(\alpha=1\) is technically single-loop — but it only converges locally, and from a poor start it diverges. Anything you would actually ship is damped.

m = x0;  sigma = sigma0;  C = I;  p_sigma = 0;  p_c = 0
for k = 1..N:
    for i = 1..lam:                           # ◄── INNER LOOP (population, not step size)
        x[i]   = m + sigma * chol(C) @ randn(n)
        fit[i] = f(x[i])                      # no gradient anywhere
    sort x by fit
    m_old = m
    m     = sum(w[i]*x[i] for i in 1..mu)     # weighted mean of the best mu
    sigma = sigma * exp((c_s/d_s)*(|p_sigma|/E_norm - 1))   # step size adapts itself
    C     = (1-c1-cmu)*C + c1*outer(p_c,p_c) + cmu*rank_mu_update

2 LOOPS But the inner loop is over a population, not a step size. There is no \(d_k\), no \(\alpha_k\), no single iterate at all — the state is \((m,\sigma,C)\) and what moves is a distribution, not a point. That is the honest signal that global search is a different kind of object, not another row in the same table.

Side by side

🗺️ The Two Axes Underneath

"Gradient descent versus conjugate gradient" is a common mental model, and it is too small. Nearly every method above is really

and the families are just choices of the matrix \(B_k\), crossed with choices of \(\alpha_k\). Two axes, not one list:

⌨️ Coding It Yourself

Start with gradient descent, and start with it fixed. It is the only single-loop method, so the entire algorithm is the master equation and nothing else:

for k in range(N):
    x -= eta * grad(x)

Write conjugate gradient instead and you need a line-search subroutine, which is more code than the method itself — plus Wolfe conditions, plus bracketing and zoom. You end up debugging the scaffolding rather than seeing the idea.