5. Nonlinear Optimization
Linear models assume constant returns: twice the production costs twice as much, and every constraint is a straight wall. Real systems bend: holding and ordering costs trade off along a curve, risk grows with the square of exposure, and transport cost grows with Euclidean distance. This chapter covers how to recognize an optimum when the objective and constraints are curved (gradients, Hessians, KKT conditions), the algorithms that find one (gradient, Newton and quasi-Newton methods, penalty, augmented Lagrangian and SQP methods), and the convex problem classes (QP, SOCP, SDP) that can be solved as reliably as LPs and that supply bounds and relaxations for the hard problems of Chapter 4.
<a id="51-unconstrained-optimization"></a>
5.1 Unconstrained Optimization
This subchapter studies $\min_{x \in \mathbb{R}^n} f(x)$ for a smooth function $f$, with no constraints. Unconstrained problems appear directly (EOQ, regression for forecasting, single-facility location), and every constrained method in Section 5.2 solves a sequence of them, so the conditions and algorithms here are the building blocks for the rest of the chapter.
Two test functions recur throughout:
- Quadratics $f(x) = \tfrac{1}{2} x^\top Q x - b^\top x$ with $Q$ symmetric positive definite. Every smooth function looks like one near a minimum, so they predict how algorithms behave near the end of a run.
- The Rosenbrock function $f(x, y) = (1 - x)^2 + 100,(y - x^2)^2$, with minimum $f = 0$ at $(1, 1)$. It has a long, curved, narrow valley, and it is the standard stress test for optimization codes. The classic starting point is $(-1.2, 1)$.
Optimality Conditions
Why it matters. An algorithm needs a test that says "stop, this is a minimum", and an analyst needs a way to verify a solver's answer. Optimality conditions provide both. They also turn optimization into equation solving: instead of searching over all of $\mathbb{R}^n$, we look for points where the gradient is zero, and then check curvature.
Theory. Picture a ball on a hilly landscape. It can rest only where the ground is flat, so a minimum must be a point with zero slope. Flat ground is not enough, though: the ball also rests (precariously) on a hilltop or a mountain pass. Curvature tells these cases apart.
Let $f$ be twice continuously differentiable, with gradient $\nabla f$ and Hessian $\nabla^2 f$ (next topic).
| Condition | Statement | Role |
|---|---|---|
| First-order necessary | $x^$ local minimum $\Rightarrow \nabla f(x^) = 0$ | candidates are the stationary points |
| Second-order necessary | $x^$ local minimum $\Rightarrow \nabla^2 f(x^)$ positive semidefinite | rules out maxima and saddles |
| Second-order sufficient | $\nabla f(x^) = 0$ and $\nabla^2 f(x^)$ positive definite $\Rightarrow$ strict local minimum | certifies a local minimum |
| Convex case | $f$ convex: $\nabla f(x^) = 0 \Leftrightarrow x^$ global minimum | no second-order check needed |
Why the first-order condition holds: by Taylor's theorem, $f(x - t\nabla f(x)) \approx f(x) - t\lVert \nabla f(x) \rVert^2$. If the gradient were nonzero, a small step against it would lower $f$, so $x$ could not be a minimum.
At a stationary point, the eigenvalues of the Hessian classify the point:
Hessian PD (all eig > 0) Hessian ND (all eig < 0) Hessian indefinite (mixed signs)
\ / .-'''-. \ along y: up /
\ / / \ '.____ x* ____.'
'._.' / \ along x: down (saddle)
local min local max neither: a mountain pass
When the Hessian is only semidefinite (some eigenvalue is 0), the test is inconclusive: at $x = 0$, the functions $x^4$ (minimum), $-x^4$ (maximum) and $x^3$ (neither) all have $f' = f'' = 0$. Higher-order terms decide.
Example: classify stationary points. For $f(x, y) = x^3 - 3x + y^2$, the gradient $(3x^2 - 3,\ 2y)$ vanishes at $(1, 0)$ and $(-1, 0)$. The Hessian is $\operatorname{diag}(6x, 2)$. The code also applies the 1D conditions to the EOQ cost $C(Q) = KD/Q + hQ/2$ (ordering cost $K$ per order, annual demand $D$, holding cost $h$ per unit per year, Section 10.1), whose stationary point is the familiar $Q^* = \sqrt{2KD/h}$.
import numpy as np
# f(x, y) = x^3 - 3x + y^2
def grad(p):
x, y = p
return np.array([3 * x**2 - 3, 2 * y])
def hess(p):
x, y = p
return np.array([[6 * x, 0.0], [0.0, 2.0]])
def classify(p):
eig = np.linalg.eigvalsh(hess(p))
if np.all(eig > 0): return eig, "strict local minimum (Hessian PD)"
if np.all(eig < 0): return eig, "strict local maximum (Hessian ND)"
if eig.min() < 0 < eig.max(): return eig, "saddle point (Hessian indefinite)"
return eig, "inconclusive (Hessian singular)"
for p in [np.array([1.0, 0.0]), np.array([-1.0, 0.0])]:
eig, verdict = classify(p)
print(f"point {p}: grad = {grad(p)}, eigenvalues = {eig}, -> {verdict}")
# EOQ: total cost C(Q) = K D / Q + h Q / 2, Q > 0
K, D, h = 100.0, 1200.0, 6.0
Q = np.sqrt(2 * K * D / h)
dC = -K * D / Q**2 + h / 2
d2C = 2 * K * D / Q**3
print(f"EOQ Q* = {Q:.1f}, C'(Q*) = {dC:.1e}, C''(Q*) = {d2C:.4f} > 0, cost = {K*D/Q + h*Q/2:.1f}")
# point [1. 0.]: grad = [0. 0.], eigenvalues = [2. 6.], -> strict local minimum (Hessian PD)
# point [-1. 0.]: grad = [0. 0.], eigenvalues = [-6. 2.], -> saddle point (Hessian indefinite)
# EOQ Q* = 200.0, C'(Q*) = 0.0e+00, C''(Q*) = 0.0300 > 0, cost = 1200.0
Note that $(1, 0)$ is only a local minimum: $f \to -\infty$ as $x \to -\infty$, so no global minimum exists. The EOQ cost is convex for $Q > 0$ ($C'' = 2KD/Q^3 > 0$), so its stationary point is the global minimum.
Key takeaways
- Local minima of smooth functions are stationary points: $\nabla f(x^*) = 0$.
- At a stationary point, a positive definite Hessian certifies a strict local minimum, and an indefinite one means a saddle.
- A semidefinite Hessian leaves the question open. Higher-order terms or problem structure must decide.
- For convex $f$, the first-order condition alone certifies a global minimum.
Practice
- (Beginner) Find and classify all stationary points of $f(x, y) = x^2 + xy + y^2 - 3x$. Is the minimum global? Justify with convexity.
- (Intermediate) For $f(x, y) = x^4 + y^4 - 4xy$, find all stationary points and classify them with the Hessian. Verify your answer numerically by evaluating $f$ on a small grid around each point.
- (Interview) "Your optimizer stopped with gradient norm $10^{-9}$. Can you report the result as the optimum?" Hint: a zero gradient is also satisfied at saddles and maxima. Think about second-order checks, convexity and multi-start runs.
Gradient and Hessian
Why it matters. Every method in this chapter is driven by derivatives. The gradient says which way is downhill. The Hessian says how the slope changes, which decides how far to step and how hard the problem is. Getting derivatives right (and cheap) is often the difference between a solver that converges and one that wanders.
Theory. For $f: \mathbb{R}^n \to \mathbb{R}$:
- The gradient $\nabla f(x) = \left(\tfrac{\partial f}{\partial x_1}, \dots, \tfrac{\partial f}{\partial x_n}\right)$ points in the direction of steepest ascent, and its length is the slope in that direction. The rate of change along a direction $d$ is the directional derivative $\nabla f(x)^\top d$. A direction with $\nabla f(x)^\top d < 0$ is a descent direction.
- The gradient is perpendicular to the level curves ${x : f(x) = c}$. Along a contour, $f$ does not change, so the gradient has no component along it.
- The Hessian $\nabla^2 f(x)$ is the $n \times n$ matrix of second partial derivatives $\tfrac{\partial^2 f}{\partial x_i \partial x_j}$. It is symmetric for twice continuously differentiable $f$. The curvature along a direction $d$ is $d^\top \nabla^2 f(x), d$.
Together they define the second-order Taylor model, the local approximation every method relies on:
$$ f(x + s) \approx f(x) + \nabla f(x)^\top s + \tfrac{1}{2}, s^\top \nabla^2 f(x), s. $$
Gradient methods use the first two terms (a tilted plane). Newton's method uses all three (a bowl).
level curves of an elongated bowl gradient at P is perpendicular to the curve
(condition number = ratio of axes^2) and points uphill; -grad points downhill
.-------------------------.
.' .-----------------. '. ^ grad f(P)
| | .---------. | P | \
| | | * | | * | * P
| | '---------' | | / <- contour through P
'. '-----------------' .'
'-------------------------'
The eigenvalues of the Hessian are the curvatures along its principal axes. Their ratio, the condition number $\kappa = \lambda_{\max} / \lambda_{\min}$, measures how elongated the level curves are. A large $\kappa$ means a narrow valley: steep across, nearly flat along. That is what makes problems hard for gradient methods.
Computing derivatives in practice:
| Method | Cost of $\nabla f$ | Accuracy | Notes |
|---|---|---|---|
| Analytic (by hand) | cheap | exact | error-prone for large models. Always check against finite differences |
| Forward differences $\tfrac{f(x + he_i) - f(x)}{h}$ | $n$ extra evaluations | error $O(h)$, about $10^{-8}$ relative with $h \approx 10^{-8}$ | simplest fallback |
| Central differences $\tfrac{f(x + he_i) - f(x - he_i)}{2h}$ | $2n$ evaluations | error $O(h^2)$, about $10^{-10}$ to $10^{-11}$ | good for verification |
| Automatic differentiation (reverse mode) | a small constant times the cost of $f$, independent of $n$ | exact to rounding | used by AMPL, Pyomo/CasADi, JAX, PyTorch |
Example: Rosenbrock at the classic start. The code compares the analytic gradient with central differences, looks at the Hessian's eigenvalues, and measures how well the linear and quadratic Taylor models predict $f$ for shrinking steps.
import numpy as np
def f(p): # Rosenbrock "banana" function
x, y = p
return (1 - x)**2 + 100 * (y - x**2)**2
def grad(p):
x, y = p
return np.array([-2 * (1 - x) - 400 * x * (y - x**2),
200 * (y - x**2)])
def hess(p):
x, y = p
return np.array([[2 - 400 * (y - x**2) + 800 * x**2, -400 * x],
[-400 * x, 200.0]])
def fd_grad(f, p, h=1e-6): # central differences: error O(h^2)
e = np.eye(len(p))
return np.array([(f(p + h * e[i]) - f(p - h * e[i])) / (2 * h) for i in range(len(p))])
p = np.array([-1.2, 1.0]) # the classic starting point
print("analytic grad :", grad(p))
print("finite diff :", fd_grad(f, p).round(4))
H = hess(p)
print("Hessian :", H.tolist())
print("eigenvalues :", np.linalg.eigvalsh(H).round(2), " condition number:", round(np.linalg.cond(H), 1))
# Quality of the first- and second-order Taylor models along d
d = np.array([0.1, -0.05])
for t in [1.0, 0.1, 0.01]:
s = t * d
lin = f(p) + grad(p) @ s
quad = lin + 0.5 * s @ H @ s
print(f"step {t:5.2f}: true {f(p + s):9.5f} linear {lin:9.5f} quadratic {quad:9.5f}")
print("Hessian at optimum (1,1):", hess(np.array([1.0, 1.0])).tolist(),
"cond =", round(np.linalg.cond(hess(np.array([1.0, 1.0]))), 0))
# analytic grad : [-215.6 -88. ]
# finite diff : [-215.6 -88. ]
# Hessian : [[1330.0, 480.0], [480.0, 200.0]]
# eigenvalues : [ 23.63 1506.37] condition number: 63.7
# step 1.00: true 11.17000 linear 7.04000 quadratic 11.54000
# step 0.10: true 22.52862 linear 22.48400 quadratic 22.52900
# step 0.01: true 24.02885 linear 24.02840 quadratic 24.02885
# Hessian at optimum (1,1): [[802.0, -400.0], [-400.0, 200.0]] cond = 2508.0
As the step shrinks by 10, the linear model's error shrinks by about 100 (it is $O(\lVert s \rVert^2)$), and the quadratic model's error shrinks by about 1000 ($O(\lVert s \rVert^3)$). For a long step the quadratic model is wrong too, which is why Newton's method needs safeguards. The Hessian at the minimum has $\kappa \approx 2500$: the Rosenbrock valley is narrow.
Key takeaways
- The gradient points uphill, is perpendicular to level curves, and gives the slope along any direction through $\nabla f^\top d$.
- The Hessian holds curvature. Its eigenvalues give the curvature along the principal axes, and their ratio $\kappa$ measures how ill-conditioned the problem is.
- The second-order Taylor model is the local picture behind every method in this chapter.
- Always verify hand-coded derivatives with central differences, or use automatic differentiation.
Practice
- (Beginner) Compute the gradient and Hessian of $f(x, y) = e^{x + y} + x^2 - 2y$ by hand, and check them with the finite-difference function from the example.
- (Intermediate) For $f(x) = \tfrac{1}{2} x^\top Q x - b^\top x$, show that $\nabla f(x) = Qx - b$ and $\nabla^2 f(x) = Q$ (assuming $Q$ symmetric). What changes if $Q$ is not symmetric?
- (Interview) "A model has 50,000 variables and its gradient is computed with forward differences. Each run takes hours. What would you change?" Hint: count function evaluations per gradient, then compare with reverse-mode automatic differentiation and with exploiting sparsity.
Line Search Methods
Why it matters. Knowing the downhill direction is not enough. A step that is too long overshoots into a region where $f$ is larger (or diverges), and a step that is too short wastes thousands of iterations. Line searches choose the step length so that every iteration makes guaranteed progress, which is what turns a local model into a globally convergent algorithm.
Theory. A line search method iterates
$$ x_{k+1} = x_k + \alpha_k d_k, $$
where $d_k$ is a descent direction ($\nabla f(x_k)^\top d_k < 0$) chosen by the method (gradient, Newton, BFGS, ...), and the step length $\alpha_k > 0$ is chosen by the line search. Along the ray, $f$ becomes a 1D function $\phi(\alpha) = f(x_k + \alpha d_k)$ with negative slope $\phi'(0) = \nabla f(x_k)^\top d_k$ at $\alpha = 0$.
- Exact line search minimizes $\phi$. It is only practical when $\phi$ has a closed form. For a quadratic $f$ with Hessian $Q$, $\alpha^* = -\nabla f^\top d / (d^\top Q d)$.
- Inexact line search accepts any step that is "good enough", tested by the Wolfe conditions with constants $0 < c_1 < c_2 < 1$:
$$ \begin{aligned} &\text{(Armijo, sufficient decrease)} && f(x_k + \alpha d_k) \le f(x_k) + c_1, \alpha, \nabla f(x_k)^\top d_k \ &\text{(curvature)} && \nabla f(x_k + \alpha d_k)^\top d_k \ge c_2, \nabla f(x_k)^\top d_k \end{aligned} $$
Armijo forbids steps that are too long: the decrease must be at least a fraction $c_1$ of what the tangent line promises. The curvature condition forbids steps that are too short: the slope must have flattened noticeably. The strong Wolfe variant bounds $\lvert \nabla f(x_k + \alpha d_k)^\top d_k \rvert$, which also stops steps far past the minimum along the ray. Typical values are $c_1 = 10^{-4}$, and $c_2 = 0.9$ for Newton and quasi-Newton methods or $c_2 = 0.1$ for conjugate gradients.
phi(alpha) = f(x + alpha d)
phi |
|* *
| '*. . . . . . . . . . . . . . . . . . . . . * <- Armijo line, slope c1 phi'(0)
| \ '*. *
| \ '*. .*
| \ '*._ _.*'
| \ ''*--.._________..--*''
| tangent line, slope phi'(0)
+-----|========================================|----> alpha
too short: acceptable too long:
slope still steep above the Armijo line
(curvature fails) (Armijo fails)
Backtracking is the simplest practical rule: try $\alpha = 1$ (the natural step for Newton-type methods), and multiply by $\rho \in (0, 1)$, typically $0.5$, until Armijo holds. It never produces a step that is too short by more than a factor $\rho$, so it usually works without the curvature test.
Convergence guarantee (Zoutendijk). If $\nabla f$ is Lipschitz continuous, $f$ is bounded below, the steps satisfy Wolfe, and the directions stay away from orthogonal to the gradient ($\cos \theta_k \ge \delta > 0$), then $\nabla f(x_k) \to 0$. Line searches thus turn locally motivated directions into a globally convergent method: it reaches a stationary point from any start (not necessarily a global minimum).
Trust regions are the main alternative. Choose a maximum step length (the radius) first, minimize the quadratic model within that radius, then expand or shrink the radius depending on how well the model predicted the actual decrease.
| Line search | Trust region | |
|---|---|---|
| Order of choices | direction first, then length | length (radius) first, then direction |
| Indefinite Hessian | needs a modification to get a descent direction | handled naturally by the constrained subproblem |
| Typical codes | BFGS, L-BFGS, nonlinear CG | trust-exact, trust-ncg in SciPy, KNITRO's Interior/CG algorithm |
Example: backtracking on Rosenbrock. At $(-1.2, 1)$ the gradient is huge, so a unit step along $-\nabla f$ lands at a point with $f \approx 2 \times 10^{11}$. Backtracking halves the step until Armijo accepts it:
import numpy as np
def f(p):
x, y = p
return (1 - x)**2 + 100 * (y - x**2)**2
def grad(p):
x, y = p
return np.array([-2 * (1 - x) - 400 * x * (y - x**2), 200 * (y - x**2)])
def backtracking(f, g, x, d, alpha=1.0, rho=0.5, c1=1e-4, verbose=False):
"""Armijo backtracking: shrink alpha until the sufficient-decrease test passes."""
fx, slope = f(x), g(x) @ d # slope < 0 for a descent direction
while f(x + alpha * d) > fx + c1 * alpha * slope:
if verbose:
print(f" alpha = {alpha:.6f}: f = {f(x + alpha * d):12.4f} > Armijo bound {fx + c1 * alpha * slope:9.4f} reject")
alpha *= rho
if verbose:
print(f" alpha = {alpha:.6f}: f = {f(x + alpha * d):12.4f} <= Armijo bound {fx + c1 * alpha * slope:9.4f} accept")
return alpha
def wolfe_ok(f, g, x, d, alpha, c1=1e-4, c2=0.9):
armijo = f(x + alpha * d) <= f(x) + c1 * alpha * g(x) @ d
curvature = g(x + alpha * d) @ d >= c2 * g(x) @ d # slope has flattened enough
return armijo, curvature
x0 = np.array([-1.2, 1.0])
d = -grad(x0) # steepest-descent direction
print(f"f(x0) = {f(x0):.2f}, directional derivative = {grad(x0) @ d:.1f}")
a = backtracking(f, grad, x0, d, verbose=True)
for alpha in [a, 1e-6]:
arm, curv = wolfe_ok(f, grad, x0, d, alpha)
print(f"alpha = {alpha:.2e}: Armijo {arm}, curvature {curv}")
# Exact line search on a quadratic: alpha* = g'g / g'Qg (closed form)
Q = np.array([[10.0, 0.0], [0.0, 1.0]]); b = np.zeros(2)
x = np.array([1.0, 10.0]); g = Q @ x - b
alpha_star = (g @ g) / (g @ Q @ g)
print(f"quadratic: exact step alpha* = {alpha_star:.4f}; new gradient is orthogonal to old:",
round(float((Q @ (x - alpha_star * g) - b) @ g), 12) + 0.0)
# f(x0) = 24.20, directional derivative = -54227.4
# alpha = 1.000000: f = 210482437168.5200 > Armijo bound 18.7773 reject
# alpha = 0.500000: f = 12810991198.7200 > Armijo bound 21.4886 reject
# alpha = 0.250000: f = 758614013.3000 > Armijo bound 22.8443 reject
# alpha = 0.125000: f = 42388850.4531 > Armijo bound 23.5222 reject
# alpha = 0.062500: f = 2078788.2100 > Armijo bound 23.8611 reject
# alpha = 0.031250: f = 72456.4239 > Armijo bound 24.0305 reject
# alpha = 0.015625: f = 543.5463 > Armijo bound 24.1153 reject
# alpha = 0.007812: f = 211.3521 > Armijo bound 24.1576 reject
# alpha = 0.003906: f = 149.6412 > Armijo bound 24.1788 reject
# alpha = 0.001953: f = 35.1074 > Armijo bound 24.1894 reject
# alpha = 0.000977: f = 5.1011 <= Armijo bound 24.1947 accept <- first step passing Armijo (10 halvings)
# alpha = 9.77e-04: Armijo True, curvature True
# alpha = 1.00e-06: Armijo True, curvature False <- too short: rejected by the curvature condition
# quadratic: exact step alpha* = 0.1818; new gradient is orthogonal to old: 0.0
A tiny step like $10^{-6}$ passes Armijo but fails the curvature condition, which is exactly the case the second Wolfe condition exists to reject. On the quadratic, the exact step makes the new gradient orthogonal to the old one. That property explains the zigzag of gradient descent in the next topic.
Key takeaways
- A line search method chooses a descent direction, then a step length along it.
- Armijo rejects steps that are too long. The curvature condition rejects steps that are too short. Together they are the Wolfe conditions.
- Backtracking from $\alpha = 1$ is simple, cheap and adequate for Newton-type methods.
- With Wolfe steps and directions not orthogonal to the gradient, $\nabla f(x_k) \to 0$ from any starting point.
Practice
- (Beginner) For $f(x) = x^2$ at $x_0 = 1$ with $d = -f'(1) = -2$, find all step lengths that satisfy Armijo with $c_1 = 0.5$. Which step does backtracking from $\alpha = 1$ with $\rho = 0.5$ return?
- (Intermediate) Derive the exact step $\alpha^* = -\nabla f^\top d / (d^\top Q d)$ for a quadratic by setting $\phi'(\alpha) = 0$. Then prove that after an exact step, $\nabla f(x_{k+1})^\top d_k = 0$.
- (Interview) "Why don't production solvers do an exact line search?" Hint: compare the cost of many function evaluations per iteration with the benefit, given that the direction itself is only approximately right far from the solution.
Gradient Descent
Why it matters. Gradient descent (steepest descent) is the simplest optimization algorithm and the ancestor of the stochastic methods that train machine learning models. Each iteration is cheap: one gradient, no matrices. Its weakness on ill-conditioned problems is the reason Newton, quasi-Newton and conjugate gradient methods exist, so understanding why it is slow explains the rest of this subchapter.
Theory. Gradient descent takes the direction of steepest descent $d_k = -\nabla f(x_k)$:
$$ x_{k+1} = x_k - \alpha_k \nabla f(x_k). $$
Analogy: a hiker in fog who always walks straight down the steepest local slope. On a round hill this leads straight to the bottom. In a long, narrow ravine, the steepest direction points mostly at the opposite wall, so the hiker crosses the ravine again and again while making slow progress along it.
Step length. With a fixed step, convergence requires $\alpha < 2/L$, where $L$ is the Lipschitz constant of $\nabla f$ (for a quadratic, the largest Hessian eigenvalue). The standard safe choice is $\alpha = 1/L$. With a line search, no knowledge of $L$ is needed.
Rate of convergence. On a strongly convex quadratic with exact line search,
$$ f(x_{k+1}) - f^* \le \left( \frac{\kappa - 1}{\kappa + 1} \right)^2 \big( f(x_k) - f^* \big), $$
so the method converges linearly and needs on the order of $\kappa \log(1/\varepsilon)$ iterations. With $\kappa = 1$ (round level curves) one step suffices. With $\kappa = 1000$ each step removes only about 0.4% of the error. The cause is geometric: an exact line search makes consecutive gradients orthogonal, so the path zigzags.
y f = 5x^2 + 0.5y^2 (kappa = 10), exact line search from (1, 10)
10.0 | 1 the iterates 1, 2, 3, ... alternate sides
8.2 | 2 of the valley floor x = 0. Consecutive
6.7 | 3 steps are perpendicular, so the path
5.5 | 4 zigzags across the valley instead of
4.5 | 5 running down it to the minimum o
3.7 | 6
... |
0.0 +----------o---------- x
-1.0 0.0 1.0
Remedies. Rescale variables so the level curves are rounder (preconditioning), add momentum (heavy-ball or Nesterov acceleration, which reduces the iteration count to the order of $\sqrt{\kappa}$), or use curvature information (Newton, BFGS, CG, the next topics). For objectives that are sums over millions of data points, stochastic gradient descent uses the gradient of a random sample of terms per step.
Example: the cost of ill-conditioning. Gradient descent with exact steps on $f = \tfrac{1}{2} x^\top Q x$ with $Q = \operatorname{diag}(\kappa, 1)$:
import numpy as np
def gd_quadratic(Q, x0, tol=1e-6, max_iter=100_000):
"""Steepest descent with exact line search on f(x) = 1/2 x'Qx (minimum at 0)."""
x, path = x0.astype(float), [x0.astype(float)]
for k in range(max_iter):
g = Q @ x
if np.linalg.norm(g) < tol:
return k, path
alpha = (g @ g) / (g @ Q @ g) # exact step for a quadratic
x = x - alpha * g
path.append(x)
return max_iter, path
for kappa in [1, 10, 100, 1000]:
Q = np.diag([kappa, 1.0]) # condition number = kappa
iters, path = gd_quadratic(Q, np.array([1.0, kappa])) # worst-case start
rate = ((kappa - 1) / (kappa + 1))**2
print(f"kappa = {kappa:5d}: {iters:6d} iterations (theory: f shrinks by {rate:.4f} per step)")
_, path = gd_quadratic(np.diag([10.0, 1.0]), np.array([1.0, 10.0]))
print("first iterates (kappa = 10):", [p.round(3).tolist() for p in path[:5]])
# Gradient descent with a FIXED step: must satisfy alpha < 2 / L (L = largest eigenvalue)
Q = np.diag([10.0, 1.0])
for alpha in [0.05, 0.19, 0.21]:
x = np.array([1.0, 1.0])
for _ in range(100):
x = x - alpha * (Q @ x)
print(f"fixed alpha = {alpha:.2f} (2/L = 0.20): |x_100| = {np.linalg.norm(x):.2e}")
# kappa = 1: 1 iterations (theory: f shrinks by 0.0000 per step)
# kappa = 10: 83 iterations (theory: f shrinks by 0.6694 per step)
# kappa = 100: 939 iterations (theory: f shrinks by 0.9608 per step)
# kappa = 1000: 10535 iterations (theory: f shrinks by 0.9960 per step)
# first iterates (kappa = 10): [[1.0, 10.0], [-0.818, 8.182], [0.669, 6.694], [-0.548, 5.477], [0.448, 4.481]]
# fixed alpha = 0.05 (2/L = 0.20): |x_100| = 5.92e-03
# fixed alpha = 0.19 (2/L = 0.20): |x_100| = 2.66e-05 <- converges, oscillating
# fixed alpha = 0.21 (2/L = 0.20): |x_100| = 1.38e+04 <- diverges
The iteration count grows linearly with $\kappa$, as the theory predicts. The fixed-step runs show the sharp threshold $2/L$: just below it the method converges, and just above it the iterates oscillate with growing amplitude. On Rosenbrock from $(-1.2, 1)$, gradient descent with backtracking needs about 19,000 iterations to reach $\lVert \nabla f \rVert < 10^{-8}$. Newton's method (next topic) needs 22.
Key takeaways
- Gradient descent steps along $-\nabla f$. Each iteration is cheap: one gradient evaluation.
- Fixed steps must satisfy $\alpha < 2/L$. Line searches remove the need to know $L$.
- Convergence is linear with ratio $\big(\tfrac{\kappa - 1}{\kappa + 1}\big)^2$, so iterations grow with the condition number.
- The zigzag comes from orthogonal consecutive gradients. Preconditioning, momentum and curvature-based methods fix it.
Practice
- (Beginner) Do two iterations of gradient descent by hand on $f(x, y) = x^2 + 4y^2$ from $(4, 1)$ with fixed step $\alpha = 0.1$. Is $\alpha = 0.3$ a safe fixed step? What is $2/L$ here?
- (Intermediate) Rescale the variables of $f = 5x^2 + 0.5y^2$ with $x = u/\sqrt{10}$, $y = v$, and run the example code on the scaled function. How many iterations does it take now, and why?
- (Interview) "Gradient descent is still the default for training neural networks with billions of parameters. Why not use Newton's method?" Hint: compare memory and cost per iteration ($n^2$ entries, $O(n^3)$ solves), and consider stochastic gradients over huge datasets.
Newton's Method
Why it matters. Newton's method uses curvature to jump straight to the minimum of the local quadratic model. Near a solution it converges quadratically, roughly doubling the number of correct digits per iteration. It is the engine inside interior point methods (Section 2.6) and sequential quadratic programming (Section 5.2).
Theory. Minimize the second-order Taylor model $f(x_k) + \nabla f(x_k)^\top s + \tfrac{1}{2} s^\top \nabla^2 f(x_k), s$. Setting its gradient to zero gives the Newton step:
$$ \nabla^2 f(x_k), s_k = -\nabla f(x_k), \qquad x_{k+1} = x_k + s_k. $$
In 1D this is $x_{k+1} = x_k - f'(x_k) / f''(x_k)$, which is Newton's root-finding method applied to $f'(x) = 0$. If $f$ is a convex quadratic, one step lands exactly on the minimum, whatever the condition number. Newton's method is affine invariant: rescaling the variables does not change its iterates, so ill-conditioning that cripples gradient descent does not affect it.
f | Newton fits a parabola (the quadratic model)
| \ .' at x_k and jumps to the parabola's minimum.
| \ f .-'
| \ .' Gradient descent only uses the tangent line,
| '. .-' which says "downhill" but not "how far".
| '-.____.-'
| ^ ^
| x_k+1 x_k
+-----------------------> x
Convergence. If $\nabla^2 f(x^)$ is positive definite and the Hessian is Lipschitz continuous, then for $x_0$ close enough to $x^$,
$$ \lVert x_{k+1} - x^* \rVert \le C, \lVert x_k - x^* \rVert^2 . $$
An error of $10^{-2}$ becomes about $10^{-4}$, then $10^{-8}$, then $10^{-16}$.
What can go wrong far from the solution, and the fixes:
| Problem | Symptom | Standard fix |
|---|---|---|
| Hessian not positive definite | Newton step points uphill or toward a saddle or maximum | modify the Hessian: add $\delta I$ until a Cholesky factorization succeeds, or fall back to $-\nabla f$ |
| Quadratic model inaccurate | full step increases $f$ or diverges | damped Newton: line search on $\alpha$ (backtracking from 1), or a trust region |
| Hessian unavailable or too costly | $n(n+1)/2$ second derivatives and an $O(n^3)$ solve per iteration | quasi-Newton (BFGS), or Newton-CG with Hessian-vector products |
Example: damped Newton on Rosenbrock. The code uses a Cholesky factorization as the positive definiteness test and backtracking for the step length:
import numpy as np
def f(p):
x, y = p
return (1 - x)**2 + 100 * (y - x**2)**2
def grad(p):
x, y = p
return np.array([-2 * (1 - x) - 400 * x * (y - x**2), 200 * (y - x**2)])
def hess(p):
x, y = p
return np.array([[2 - 400 * (y - x**2) + 800 * x**2, -400 * x], [-400 * x, 200.0]])
def damped_newton(x, tol=1e-10, max_iter=100):
for k in range(max_iter):
g, H = grad(x), hess(x)
if np.linalg.norm(g) < tol:
return x, k
try:
np.linalg.cholesky(H) # succeeds only if H is positive definite
d = np.linalg.solve(H, -g) # Newton direction
except np.linalg.LinAlgError:
d = -g # fallback: steepest descent
alpha = 1.0
while f(x + alpha * d) > f(x) + 1e-4 * alpha * g @ d: # Armijo backtracking
alpha *= 0.5
x = x + alpha * d
print(f"iter {k + 1:2d}: step {alpha:6.4f} x = ({x[0]: .8f}, {x[1]: .8f}) "
f"error = {np.linalg.norm(x - 1):.2e}")
return x, max_iter
x, k = damped_newton(np.array([-1.2, 1.0]))
print(f"converged in {k} iterations")
# iter 1: step 1.0000 x = (-1.17528090, 1.38067416) error = 2.21e+00
# iter 2: step 0.1250 x = (-0.93298143, 0.81121066) error = 1.94e+00 <- full step rejected, damped
# iter 3: step 1.0000 x = (-0.78254008, 0.58973638) error = 1.83e+00
# ...
# iter 16: step 1.0000 x = ( 0.94207869, 0.88133620) error = 1.32e-01
# iter 17: step 1.0000 x = ( 0.96799182, 0.93633667) error = 7.13e-02
# iter 18: step 1.0000 x = ( 0.99621031, 0.99163870) error = 9.18e-03
# iter 19: step 1.0000 x = ( 0.99947938, 0.99894834) error = 1.17e-03 <- local phase: error roughly squares
# iter 20: step 1.0000 x = ( 0.99999889, 0.99999751) error = 2.73e-06
# iter 21: step 1.0000 x = ( 1.00000000, 1.00000000) error = 1.35e-10
# iter 22: step 1.0000 x = ( 1.00000000, 1.00000000) error = 0.00e+00
# converged in 22 iterations
The run has two phases. In the global phase (iterations 1 to about 17) the method creeps along the curved valley and backtracking cuts some steps. In the local phase the full step is accepted every time and the error is roughly squared at each iteration: $10^{-3} \to 10^{-6} \to 10^{-10} \to 0$.
Key takeaways
- The Newton step solves $\nabla^2 f, s = -\nabla f$: it minimizes the local quadratic model.
- Near a minimum with a positive definite Hessian, convergence is quadratic and independent of conditioning.
- Far away, safeguard it with a line search (damped Newton) or a trust region, and modify an indefinite Hessian.
- The price is second derivatives and an $O(n^3)$ linear solve per iteration, which motivates quasi-Newton methods.
Practice
- (Beginner) Apply Newton's method by hand to $f(x) = x - \ln x$ from $x_0 = 0.5$. Compute three iterates and the error $\lvert x_k - 1 \rvert$ at each. Do the digits double?
- (Intermediate) Run pure Newton (no line search, no Hessian check) on $f(x) = \sqrt{1 + x^2}$ from $x_0 = 0.5$ and from $x_0 = 1.5$. Derive the iteration in closed form ($x_{k+1} = -x_k^3$) and explain the divergence. Then add backtracking and repeat.
- (Interview) "Newton's method converges in a handful of iterations. Why is it not the default everywhere?" Hint: consider the cost of forming and factorizing the Hessian for large $n$, the effort of coding second derivatives, and behavior on non-convex functions.
Quasi-Newton Methods (BFGS)
Why it matters. BFGS gives most of Newton's speed using only gradients. It is the default unconstrained method in SciPy and many other libraries, and its limited-memory variant L-BFGS scales to millions of variables. If you remember one general-purpose method for smooth unconstrained problems, it should be BFGS.
Theory. Instead of computing the Hessian, build an approximation $B_k \approx \nabla^2 f(x_k)$ from the gradients already observed. Analogy: you can estimate how sharply a road bends by noting how much the direction of travel changed between two points, without a map of the curvature.
After a step $s_k = x_{k+1} - x_k$, the gradient change is $y_k = \nabla f(x_{k+1}) - \nabla f(x_k)$. For a quadratic, $y_k = Q s_k$ exactly. So we require the new approximation to satisfy the secant condition:
$$ B_{k+1}, s_k = y_k . $$
This does not determine $B_{k+1}$ uniquely. BFGS (Broyden, Fletcher, Goldfarb, Shanno, 1970) chooses the symmetric matrix closest to the previous one. It is usually written for the inverse approximation $H_k \approx \nabla^2 f(x_k)^{-1}$, so the direction costs only a matrix-vector product:
$$ H_{k+1} = \left(I - \rho_k s_k y_k^\top\right) H_k \left(I - \rho_k y_k s_k^\top\right) + \rho_k s_k s_k^\top, \qquad \rho_k = \frac{1}{y_k^\top s_k}. $$
The algorithm is: $d_k = -H_k \nabla f(x_k)$, a Wolfe line search, then the update. Two facts make it robust:
- If $y_k^\top s_k > 0$ (the curvature condition, which the Wolfe line search guarantees), the update keeps $H_{k+1}$ positive definite, so $d_k$ is always a descent direction.
- Near the solution, $H_k$ approaches the true inverse Hessian along the directions explored, and convergence is superlinear: the ratio of successive errors goes to zero.
Variants:
| Method | Idea | Use |
|---|---|---|
| BFGS | rank-2 update of $H$, stays positive definite | default for small and medium $n$ (up to a few thousand) |
| L-BFGS | store only the last $m$ pairs $(s_k, y_k)$, typically $m = 5$ to $20$, and apply $H_k$ implicitly | large $n$: memory $O(mn)$ instead of $O(n^2)$ |
| DFP | the "dual" update, of $B$ instead of $H$ | historical. Less robust than BFGS |
| SR1 | symmetric rank-1 update, may become indefinite | trust-region methods, where indefinite models are useful |
Example: BFGS on Rosenbrock, and a comparison. The hand-written BFGS below starts from $H_0 = I$. The comparison uses SciPy's implementations with the same tolerance.
import numpy as np
from scipy.optimize import minimize, rosen, rosen_der, rosen_hess
def bfgs(f, grad, x, tol=1e-8, max_iter=500):
n = len(x)
Hinv = np.eye(n) # initial inverse-Hessian guess
g = grad(x)
for k in range(max_iter):
if np.linalg.norm(g) < tol:
return x, k, Hinv
d = -Hinv @ g # quasi-Newton direction
alpha = 1.0
while f(x + alpha * d) > f(x) + 1e-4 * alpha * g @ d:
alpha *= 0.5 # Armijo backtracking
s = alpha * d # step taken
x_new = x + s
g_new = grad(x_new)
y = g_new - g # change in gradient
if y @ s > 1e-12: # curvature condition: keep Hinv PD
rho = 1.0 / (y @ s)
I = np.eye(n)
Hinv = (I - rho * np.outer(s, y)) @ Hinv @ (I - rho * np.outer(y, s)) \
+ rho * np.outer(s, s) # BFGS inverse update
x, g = x_new, g_new
return x, max_iter, Hinv
x0 = np.array([-1.2, 1.0])
x, k, Hinv = bfgs(rosen, rosen_der, x0)
print(f"hand-written BFGS: x = {x.round(6)}, {k} iterations")
print("BFGS inverse-Hessian estimate:", Hinv.round(3).tolist())
print("true inverse Hessian at x* :", np.linalg.inv(rosen_hess(x)).round(3).tolist())
print(f"{'method':12s} {'iters':>6s} {'f evals':>8s} {'grad evals':>11s} {'Hess evals':>11s}")
for method, kw in [("CG", {}), ("BFGS", {}), ("L-BFGS-B", {}),
("trust-exact", {"hess": rosen_hess})]:
r = minimize(rosen, x0, jac=rosen_der, method=method, tol=1e-10, **kw)
print(f"{method:12s} {r.nit:6d} {r.nfev:8d} {r.njev:11d} {getattr(r, 'nhev', 0):11d}"
f" x = {r.x.round(6)}")
# hand-written BFGS: x = [1. 1.], 35 iterations
# BFGS inverse-Hessian estimate: [[0.495, 0.989], [0.989, 1.983]]
# true inverse Hessian at x* : [[0.5, 1.0], [1.0, 2.005]]
# method iters f evals grad evals Hess evals
# CG 37 80 79 0 x = [1. 1.]
# BFGS 34 41 41 0 x = [1. 1.]
# L-BFGS-B 37 45 45 0 x = [1. 1.]
# trust-exact 26 27 24 27 x = [1. 1.]
BFGS learned the inverse Hessian at the solution to within about 1% without computing a single second derivative. Compare the methods on this problem:
| Method | Derivatives used | Cost per iteration | Memory | Local convergence | Iterations (Rosenbrock) |
|---|---|---|---|---|---|
| Gradient descent | $\nabla f$ | $O(n)$ | $O(n)$ | linear, depends on $\kappa$ | about 19,000 |
| Nonlinear CG | $\nabla f$ | $O(n)$ | $O(n)$ | linear | 37 |
| BFGS | $\nabla f$ | $O(n^2)$ | $O(n^2)$ | superlinear | 34 |
| L-BFGS | $\nabla f$ | $O(mn)$ | $O(mn)$ | linear, fast in practice | 37 |
| Newton (trust region) | $\nabla f$, $\nabla^2 f$ | $O(n^3)$ | $O(n^2)$ | quadratic | 26 |
Key takeaways
- Quasi-Newton methods build curvature information from gradient differences through the secant condition $B_{k+1} s_k = y_k$.
- BFGS keeps the approximation positive definite when $y_k^\top s_k > 0$, which the Wolfe line search guarantees.
- Convergence is superlinear with only first derivatives and $O(n^2)$ work per iteration.
- L-BFGS stores a few vector pairs instead of a matrix and is the method of choice for large smooth problems.
Practice
- (Beginner) For $f(x) = x^2$ in 1D, the secant condition gives $B_{k+1} = y_k / s_k$. Compute it for $x_k = 3$, $x_{k+1} = 1$ and compare with $f''$. Why is it exact here?
- (Intermediate) Verify numerically that the BFGS update satisfies $H_{k+1} y_k = s_k$, and that $H_{k+1}$ stays positive definite when $y^\top s > 0$. Then pick a step with $y^\top s < 0$ and show that positive definiteness can be lost.
- (Interview) "You need to minimize a smooth function of 2 million variables, and you only have code for the function and its gradient. What do you use and how do you tune it?" Hint: L-BFGS, the history length $m$, gradient checking, scaling, and the line search parameters.
Conjugate Gradient Method
Why it matters. The conjugate gradient (CG) method minimizes a convex quadratic, or equivalently solves a symmetric positive definite linear system $Ax = b$, using only matrix-vector products and a few vectors of memory. It is the standard solver for huge sparse systems, it computes Newton steps inside Newton-CG and some interior point codes, and its nonlinear extension is a low-memory alternative to BFGS.
Theory. Minimizing $f(x) = \tfrac{1}{2} x^\top A x - b^\top x$ is the same as solving $Ax = b$, since $\nabla f(x) = Ax - b$. Gradient descent is slow because each step partly undoes earlier progress: improving along one direction spoils the minimization along a previous one. CG avoids this by choosing directions that are $A$-conjugate:
$$ d_i^\top A, d_j = 0 \quad \text{for } i \neq j . $$
Conjugate directions are orthogonal in the geometry where the level curves become circles. Minimizing exactly along each of them in turn never disturbs the minimization along the previous ones. So after $k$ steps, $x_k$ is the exact minimizer of $f$ over the subspace spanned by the first $k$ directions (a Krylov subspace). Analogy: tuning an old radio with separate bass and treble knobs that do not interact. Set each knob once and you are done. Gradient descent is the radio where every knob affects every other.
Each new direction is built from the current residual $r_k = b - Ax_k = -\nabla f(x_k)$ and only the previous direction:
$$ \alpha_k = \frac{r_k^\top r_k}{d_k^\top A d_k}, \quad x_{k+1} = x_k + \alpha_k d_k, \quad r_{k+1} = r_k - \alpha_k A d_k, \quad \beta_k = \frac{r_{k+1}^\top r_{k+1}}{r_k^\top r_k}, \quad d_{k+1} = r_{k+1} + \beta_k d_k . $$
Convergence.
- In exact arithmetic, CG terminates in at most $n$ steps, and in at most as many steps as $A$ has distinct eigenvalues.
- In general, the error in the $A$-norm shrinks at least by the factor $\frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}$ per iteration, against $\frac{\kappa - 1}{\kappa + 1}$ for steepest descent. For $\kappa = 100$, that is $0.82$ against $0.98$.
- Preconditioning (solving an equivalent system $M^{-1}Ax = M^{-1}b$ with $M \approx A$ easy to invert) clusters the eigenvalues and makes CG much faster.
Nonlinear CG applies the same recurrence to a general $f$, with $r_k = -\nabla f(x_k)$, a line search for $\alpha_k$, and a $\beta_k$ formula such as Fletcher-Reeves (as above) or Polak-Ribiere $\beta_k = \max!\left{0,\ \frac{\nabla f_{k+1}^\top (\nabla f_{k+1} - \nabla f_k)}{\nabla f_k^\top \nabla f_k}\right}$, which restarts automatically when progress stalls. It needs only $O(n)$ memory, like gradient descent, but it is usually far faster.
Example: CG against steepest descent on a 200-variable system. Two SPD matrices with the same size: one with eigenvalues spread over $[1, 100]$, one with only five distinct eigenvalues.
import numpy as np
def conjugate_gradient(A, b, tol=1e-8):
"""Linear CG: minimizes 1/2 x'Ax - b'x (A symmetric positive definite)."""
x = np.zeros_like(b)
r = b - A @ x # residual = -gradient
d = r.copy()
for k in range(len(b)):
if np.linalg.norm(r) < tol * np.linalg.norm(b):
return x, k
Ad = A @ d # the only matrix-vector product per step
alpha = (r @ r) / (d @ Ad) # exact line search along d
x = x + alpha * d
r_new = r - alpha * Ad
beta = (r_new @ r_new) / (r @ r) # makes the next d A-conjugate to d
d = r_new + beta * d
r = r_new
return x, len(b)
def steepest_descent(A, b, tol=1e-8, max_iter=200_000):
x = np.zeros_like(b)
for k in range(max_iter):
r = b - A @ x
if np.linalg.norm(r) < tol * np.linalg.norm(b):
return x, k
x = x + (r @ r) / (r @ A @ r) * r
return x, max_iter
rng = np.random.default_rng(0)
n = 200
U, _ = np.linalg.qr(rng.standard_normal((n, n))) # random orthogonal basis
b = rng.standard_normal(n)
spectra = {"spread 1..100": np.linspace(1, 100, n),
"5 clusters": np.repeat([1.0, 10.0, 30.0, 60.0, 100.0], n // 5)}
for name, eig in spectra.items():
A = U @ np.diag(eig) @ U.T # SPD matrix with the chosen eigenvalues
x_cg, k_cg = conjugate_gradient(A, b)
x_sd, k_sd = steepest_descent(A, b)
print(f"{name:14s}: CG {k_cg:3d} iterations, steepest descent {k_sd:5d} iterations, "
f"CG error {np.linalg.norm(x_cg - np.linalg.solve(A, b)):.1e}")
# spread 1..100 : CG 64 iterations, steepest descent 788 iterations, CG error 1.1e-08
# 5 clusters : CG 5 iterations, steepest descent 860 iterations, CG error 8.1e-14
With spread eigenvalues, CG needs 64 iterations against 788 for steepest descent. With five eigenvalue clusters, CG finishes in exactly five iterations, as the theory predicts, while steepest descent still needs hundreds because its speed depends only on $\kappa = 100$.
Key takeaways
- CG minimizes a convex quadratic (solves an SPD system) with one matrix-vector product per iteration and $O(n)$ memory.
- Conjugate directions never undo earlier progress, so CG finishes in at most $n$ steps, or as many steps as there are distinct eigenvalues.
- Its rate depends on $\sqrt{\kappa}$ instead of $\kappa$. Preconditioning makes it faster still.
- Nonlinear CG (Fletcher-Reeves, Polak-Ribiere) extends the idea to general smooth functions with gradient-descent-level memory.
Practice
- (Beginner) Run two iterations of CG by hand on $A = \begin{pmatrix} 4 & 1 \ 1 & 3 \end{pmatrix}$, $b = (1, 2)$ from $x_0 = 0$. Check that $d_0^\top A d_1 = 0$ and that $x_2$ solves $Ax = b$.
- (Intermediate) Add a diagonal (Jacobi) preconditioner to the example code and test it on a matrix whose diagonal entries vary from 1 to $10^4$. Compare the iteration counts with and without the preconditioner.
- (Interview) "Your system $Ax = b$ has $10^7$ unknowns and $A$ is sparse and SPD. Why not just call a direct solver?" Hint: consider fill-in during factorization, memory, the fact that CG only needs a routine computing $Av$, and when a direct solver still wins.
<a id="52-constrained-optimization"></a>
5.2 Constrained Optimization
This subchapter adds constraints to the problem:
$$ \min_{x \in \mathbb{R}^n}\ f(x) \quad \text{s.t.} \quad g_i(x) \le 0 \ \ (i = 1, \dots, m), \qquad h_j(x) = 0 \ \ (j = 1, \dots, p). $$
It first derives optimality conditions (Lagrange, KKT) and the regularity assumptions they need, then presents the main algorithm families. Notation used throughout: $\lambda_j$ is the multiplier of the equality constraint $h_j$, $\mu_i \ge 0$ is the multiplier of the inequality $g_i$, $\rho$ is a penalty parameter and $\tau$ a barrier parameter.
Most examples use one running example, small enough to solve by hand:
$$ \min\ (x_1 - 2)^2 + (x_2 - 1)^2 \quad \text{s.t.} \quad g_1(x) = x_1^2 - x_2 \le 0, \qquad g_2(x) = x_1 + x_2 - 2 \le 0 . $$
It asks for the feasible point closest to $(2, 1)$. The region lies above the parabola $x_2 = x_1^2$ and below the line $x_1 + x_2 = 2$:
x2
4.0 | # # = feasible region
3.5 | # # x1^2 <= x2 <= 2 - x1
3.0 | # # # left/lower edge: parabola g1 = 0
2.5 | # # # # # right/upper edge: line g2 = 0
2.0 | # # # # # #
1.5 | # # # # # # # * = x* = (1, 1): both constraints active
1.0 | # # # # # # # # * o o = (2, 1): unconstrained minimizer
0.5 | # # # # # of f, infeasible
0.0 +--------------------#---------------- x1
-2 -1 0 1 2
The solution is $x^* = (1, 1)$ with $f^* = 1$ and multipliers $\mu^* = (\tfrac{2}{3}, \tfrac{2}{3})$. The topics below derive these values and compute them with each algorithm.
Lagrange Multipliers
Why it matters. Many classical models have equality constraints: a production target, a budget that must be spent, a fixed volume. Lagrange multipliers turn such a problem into a system of equations. The multipliers themselves are prices: they say how much the optimal cost changes per unit change in the constraint, just like the LP shadow prices of Section 2.5.
Theory. Consider $\min f(x)$ subject to $h_j(x) = 0$, $j = 1, \dots, p$. Analogy: a hiker must stay on a trail ($h(x) = 0$) and wants its lowest point. At that point the trail must run along a contour line of the terrain. If it crossed contours, walking one way along the trail would go further downhill. "The trail is tangent to the contour" means the terrain gradient is perpendicular to the trail, and so is the gradient of $h$. The two gradients are therefore parallel:
$$ \nabla f(x^) + \sum_{j=1}^{p} \lambda_j \nabla h_j(x^) = 0 . $$
contours of f At x*, the constraint curve h = 0
.----------------------. touches a contour of f without crossing it.
.-' .--------------. '-. grad f and grad h are both perpendicular
.' .' .------. '. '. to the curve there, hence parallel:
/ / ( min ) \ \ grad f = -lambda grad h
\ \ '------' / /
'. '. x* .' .'
'-. '-------------*-' .-'
'-------------- / -----'
h(x) = 0 ------------'
Define the Lagrangian $L(x, \lambda) = f(x) + \sum_j \lambda_j h_j(x)$. The conditions $\nabla_x L = 0$ and $h(x) = 0$ give $n + p$ equations in the $n + p$ unknowns $(x, \lambda)$. They are necessary for a local minimum when the gradients $\nabla h_j(x^*)$ are linearly independent (a constraint qualification, discussed below).
Second-order check. A stationary point of $L$ is a strict local minimum if $d^\top \nabla^2_{xx} L(x^, \lambda^), d > 0$ for every $d \neq 0$ in the tangent space ${d : \nabla h_j(x^*)^\top d = 0 \ \forall j}$. Curvature only matters along directions that stay on the constraint.
Multipliers as prices. Suppose a constraint has a right-hand side $b$ and the optimal value is $f^(b)$. By the envelope theorem, $\tfrac{d f^}{d b} = \tfrac{\partial L}{\partial b}$ at the optimum. If the constraint is written as $b - a(x) = 0$, then $\tfrac{d f^*}{d b} = \lambda$. The multiplier is the marginal cost of one more unit of requirement.
Example: splitting production between two plants. Plant 1 produces $x$ units at cost $x^2$, plant 2 produces $y$ units at cost $2y^2$, and together they must deliver $b = 3$ units. Stationarity says $2x = \lambda$ and $4y = \lambda$: at the optimum, both plants have the same marginal cost, and that common marginal cost is $\lambda$. If one plant were cheaper at the margin, shifting a unit to it would save money.
import numpy as np
# Two plants share a production target b: min x^2 + 2 y^2 s.t. x + y = b
# Write the constraint as h = b - x - y = 0 (unmet demand). L = x^2 + 2y^2 + lam (b - x - y)
# Stationarity: 2x - lam = 0, 4y - lam = 0
def solve(b):
K = np.array([[2.0, 0.0, -1.0], # d/dx: 2x - lam = 0
[0.0, 4.0, -1.0], # d/dy: 4y - lam = 0
[1.0, 1.0, 0.0]]) # constraint: x + y = b
x, y, lam = np.linalg.solve(K, [0.0, 0.0, b])
return x, y, lam, x**2 + 2 * y**2
x, y, lam, cost = solve(3.0)
print(f"x = {x:.3f}, y = {y:.3f}, lambda = {lam:.3f}, cost = {cost:.3f}")
print(f"marginal costs: plant 1 = {2 * x:.3f}, plant 2 = {4 * y:.3f} <- equal to lambda")
eps = 1e-4 # perturb the target and re-solve
dcost = (solve(3.0 + eps)[3] - solve(3.0 - eps)[3]) / (2 * eps)
print(f"d(cost)/db by finite differences = {dcost:.3f} <- lambda is the shadow price")
# x = 2.000, y = 1.000, lambda = 4.000, cost = 6.000
# marginal costs: plant 1 = 4.000, plant 2 = 4.000 <- equal to lambda
# d(cost)/db by finite differences = 4.000 <- lambda is the shadow price
Here the conditions are linear, so one linear solve finds $(x, y, \lambda) = (2, 1, 4)$. The Hessian of $L$ is $\operatorname{diag}(2, 4)$, positive definite, so this is a minimum. Raising the target by one unit raises the cost by about $\lambda = 4$.
Key takeaways
- At a constrained optimum with equality constraints, $\nabla f$ is a linear combination of the constraint gradients: $\nabla_x L = 0$.
- Stationarity plus feasibility give $n + p$ equations in $n + p$ unknowns.
- Second-order conditions only check curvature of $L$ along the tangent space of the constraints.
- A multiplier is the marginal change in the optimal value per unit change in the constraint's right-hand side: an economic price.
Practice
- (Beginner) Solve $\min x^2 + y^2$ subject to $x + 2y = 5$ with a Lagrange multiplier. Interpret the answer geometrically (closest point of a line to the origin).
- (Intermediate) A box with a square base $x \times x$ and height $y$ must hold volume $V = 32$. The base costs 2 per unit area and the four sides cost 1 per unit area (no lid). Minimize the cost with a Lagrange multiplier, and use $\lambda$ to estimate the extra cost if $V$ grows to 33. Check by re-solving.
- (Interview) "Three data centers serve a fixed total load and each has a convex power-cost curve. How should the load be split, and what single number summarizes the solution?" Hint: equal marginal costs across centers that are used, and the multiplier as the system's marginal price. What changes when a center hits its capacity?
Karush-Kuhn-Tucker Conditions
Why it matters. The KKT conditions extend Lagrange multipliers to inequality constraints. They are the central optimality theorem of nonlinear programming: every NLP solver uses them as its stopping test, and LP duality with complementary slackness (Section 2.4) is a special case. Reading a solver's output ("dual infeasibility", "complementarity") is reading KKT residuals.
Theory. Form the Lagrangian $L(x, \mu, \lambda) = f(x) + \sum_i \mu_i g_i(x) + \sum_j \lambda_j h_j(x)$. A point $x^$ is a KKT point if there exist multipliers $\mu^$ and $\lambda^*$ such that:
| Condition | Formula | Meaning |
|---|---|---|
| Stationarity | $\nabla f(x^) + \sum_i \mu_i^ \nabla g_i(x^) + \sum_j \lambda_j^ \nabla h_j(x^*) = 0$ | forces balance |
| Primal feasibility | $g_i(x^) \le 0$, $h_j(x^) = 0$ | $x^*$ is feasible |
| Dual feasibility | $\mu_i^* \ge 0$ | walls can only push, not pull |
| Complementary slackness | $\mu_i^, g_i(x^) = 0$ for every $i$ | only touching walls push |
Analogy: a ball rolls downhill inside a box. When it stops, gravity ($-\nabla f$) is balanced by the walls it touches. Each wall can only push outward along its normal ($\mu_i \ge 0$), and a wall the ball does not touch exerts no force ($g_i < 0 \Rightarrow \mu_i = 0$).
Geometrically, $-\nabla f(x^*)$ must lie in the cone spanned by the gradients of the active constraints (plus any combination of equality gradients). If it pointed outside the cone, some feasible direction would decrease $f$.
KKT at x* = (1, 1) in the running example
grad g2 = (1, 1)
^
/ -grad f(x*) = (2, 0) lies inside the cone spanned
/ by the gradients of the two active constraints:
x* = (1,1) *-------> -grad f
\ (2, 0) = 2/3 * (2, -1) + 2/3 * (1, 1)
\
v both weights are >= 0, so x* is a KKT point
grad g1 = (2, -1)
When are KKT conditions necessary, and when sufficient?
- Necessary: if $x^$ is a local minimum and a constraint qualification holds at $x^$ (next topic), multipliers satisfying KKT exist.
- Sufficient: if the problem is convex ($f$ and all $g_i$ convex, all $h_j$ affine), every KKT point is a global minimum. For non-convex problems, second-order conditions on $\nabla^2_{xx} L$ over the critical directions certify a local minimum.
Solving by hand: enumerate active sets. Complementary slackness splits the problem into cases: each constraint is either active ($g_i = 0$) or has $\mu_i = 0$. For the running example:
| Active set | Candidate | Check |
|---|---|---|
| none | $(2, 1)$ | $g_1 = 3 > 0$: infeasible |
| ${g_1}$ | $2x_1^3 - x_1 - 2 = 0 \Rightarrow (1.165, 1.358)$ | $g_2 = 0.52 > 0$: infeasible |
| ${g_2}$ | $(1.5, 0.5)$, $\mu_2 = 1$ | $g_1 = 1.75 > 0$: infeasible |
| ${g_1, g_2}$ | $(-2, 4)$: $\mu_1 = -\tfrac{14}{3}$ | $\mu_1 < 0$: not KKT |
| ${g_1, g_2}$ | $(1, 1)$: $\mu = (\tfrac{2}{3}, \tfrac{2}{3})$ | all conditions hold: KKT point |
The problem is convex ($f$ and $g_1$ are convex, $g_2$ is linear), so $(1, 1)$ is the global minimum. Enumeration needs $2^m$ cases, so it is a hand technique only. Algorithms find the active set iteratively.
Example: a KKT checker. Given a candidate $x$, the code identifies the active constraints and finds the nonnegative multipliers that best satisfy stationarity (nonnegative least squares). A zero residual means a KKT point.
import numpy as np
from scipy.optimize import nnls
# min (x1-2)^2 + (x2-1)^2 s.t. g1 = x1^2 - x2 <= 0, g2 = x1 + x2 - 2 <= 0
f_grad = lambda x: np.array([2 * (x[0] - 2), 2 * (x[1] - 1)])
g = lambda x: np.array([x[0]**2 - x[1], x[0] + x[1] - 2])
g_jac = lambda x: np.array([[2 * x[0], -1.0], [1.0, 1.0]]) # rows = constraint gradients
def kkt_check(x, tol=1e-8):
gx = g(x)
if np.any(gx > tol):
return f"infeasible (g = {gx.round(3)})"
active = np.where(np.abs(gx) <= tol)[0]
if active.size == 0:
res, mu = np.linalg.norm(f_grad(x)), np.array([])
else:
# stationarity: grad f + sum mu_i grad g_i = 0 with mu >= 0 (inactive mu_i = 0)
mu, res = nnls(g_jac(x)[active].T, -f_grad(x))
status = "KKT point" if res < 1e-6 else "NOT a KKT point"
return f"active = {[f'g{i + 1}' for i in active]}, mu = {np.round(mu, 4)}, residual = {res:.2e} -> {status}"
for x in [np.array([1.0, 1.0]), np.array([0.0, 0.0]), np.array([0.5, 0.5]), np.array([1.5, 0.5])]:
print(f"x = {x}: {kkt_check(x)}")
# x = [1. 1.]: active = ['g1', 'g2'], mu = [0.6667 0.6667], residual = 0.00e+00 -> KKT point
# x = [0. 0.]: active = ['g1'], mu = [0.], residual = 4.47e+00 -> NOT a KKT point
# x = [0.5 0.5]: active = [], mu = [], residual = 3.16e+00 -> NOT a KKT point
# x = [1.5 0.5]: infeasible (g = [1.75 0. ])
Connection to LP. For $\min c^\top x$ s.t. $Ax \ge b$, $x \ge 0$, the KKT conditions are exactly primal feasibility, dual feasibility ($A^\top y \le c$, $y \ge 0$) and complementary slackness from Section 2.4. LP duality is KKT theory for linear functions.
Key takeaways
- KKT = stationarity + primal feasibility + dual feasibility ($\mu \ge 0$) + complementary slackness ($\mu_i g_i = 0$).
- Geometrically, $-\nabla f$ lies in the cone of active constraint gradients.
- KKT conditions are necessary under a constraint qualification, and sufficient for global optimality in convex problems.
- The multiplier $\mu_i$ is the marginal value of relaxing $g_i$: relaxing $g_i \le 0$ to $g_i \le \varepsilon$ lowers the optimal cost by about $\mu_i \varepsilon$.
Practice
- (Beginner) Write the KKT conditions for $\min x^2$ s.t. $x \ge 1$, and solve them. What is the multiplier, and what does it say about moving the bound to $x \ge 1.1$?
- (Intermediate) Change the running example's target to $(0, 2)$, i.e. minimize $x_1^2 + (x_2 - 2)^2$. Enumerate active sets, find the KKT point, and confirm with the checker. Which constraints are active now?
- (Interview) "A solver reports an optimal solution where a capacity constraint has zero slack but a zero multiplier. Is that a bug?" Hint: complementary slackness allows both to be zero. Think about degeneracy and what the zero price says about adding capacity.
Constraint Qualifications
Why it matters. The KKT theorem has a fine-print assumption: multipliers are guaranteed to exist only if the constraints are "regular" at the optimum. When regularity fails, a true minimum may satisfy no KKT system, and solvers that search for KKT points can stall, report huge multipliers or fail outright. The way a constraint is written can matter as much as the set it describes.
Theory. KKT necessity comes from linearizing the constraints: near $x^$, the feasible directions should be the directions $d$ with $\nabla g_i(x^)^\top d \le 0$ for active $i$ and $\nabla h_j(x^*)^\top d = 0$. A constraint qualification (CQ) is any condition guaranteeing that this linearization describes the true local geometry. Analogy: using a road map to plan a route works only if the map matches the terrain. A CQ guarantees that the first-order "map" of the constraints is faithful.
| Constraint qualification | Condition at $x^*$ | Remarks |
|---|---|---|
| LICQ (linear independence) | gradients of active inequalities and all equalities are linearly independent | most common. Implies unique multipliers |
| MFCQ (Mangasarian-Fromovitz) | equality gradients independent, and some $d$ has $\nabla h_j^\top d = 0$ and $\nabla g_i^\top d < 0$ for active $i$ | weaker than LICQ. Multiplier set nonempty and bounded |
| Slater (convex problems) | some point $\bar x$ with $g_i(\bar x) < 0$ for all $i$ and $h(\bar x) = 0$ | global and easy to check. Also gives strong duality (Section 5.3) |
| Linearity | all active constraints are affine | always holds, so LPs and linearly constrained problems never need a CQ |
The hierarchy is LICQ $\Rightarrow$ MFCQ $\Rightarrow$ (Abadie CQ) $\Rightarrow$ KKT multipliers exist at a local minimum. CQs are sufficient, not necessary: KKT can still hold when all of them fail.
Typical failures.
- A feasible set that is a single point or has a cusp, so no direction is strictly feasible.
- Constraints written with vanishing gradients, such as $x^2 \le 0$ instead of $x = 0$, or $(x - a)^2 \le 0$.
- Complementarity constraints $x \ge 0$, $y \ge 0$, $xy = 0$ (equilibrium models, bilevel problems). They violate MFCQ at every feasible point and need specialized reformulations (MPCC methods).
Example: three failures and a fix. The code checks LICQ and the best achievable stationarity residual with $\mu \ge 0$:
import numpy as np
from scipy.optimize import nnls
def analyze(name, grad_f, active_grads):
G = np.array(active_grads, dtype=float) # rows = gradients of active constraints
licq = np.linalg.matrix_rank(G) == len(G) # linearly independent?
mu, res = nnls(G.T, -np.asarray(grad_f, dtype=float))
print(f"{name:34s} LICQ {str(licq):5s} best mu = {mu.round(3)} stationarity residual = {res:.3f}")
# (a) min x s.t. x^2 <= 0. Feasible set {0}, so x* = 0, but grad g(0) = 0.
analyze("(a) x^2 <= 0 at x* = 0", [1.0], [[0.0]])
# (b) Kuhn-Tucker cusp: min -x1 s.t. x2 - (1 - x1)^3 <= 0, -x2 <= 0. x* = (1, 0)
analyze("(b) cusp at x* = (1, 0)", [-1.0, 0.0], [[0.0, 1.0], [0.0, -1.0]])
# (c) same feasible set as (a), described linearly: x <= 0, -x <= 0
analyze("(c) x <= 0, -x <= 0 at x* = 0", [1.0], [[1.0], [-1.0]])
# (d) running example of the KKT topic at x* = (1, 1)
analyze("(d) running example at (1, 1)", [-2.0, 0.0], [[2.0, -1.0], [1.0, 1.0]])
# (a) x^2 <= 0 at x* = 0 LICQ False best mu = [0.] stationarity residual = 1.000
# (b) cusp at x* = (1, 0) LICQ False best mu = [0. 0.] stationarity residual = 1.000
# (c) x <= 0, -x <= 0 at x* = 0 LICQ False best mu = [0. 1.] stationarity residual = 0.000
# (d) running example at (1, 1) LICQ True best mu = [0.667 0.667] stationarity residual = 0.000
- (a) $\min x$ s.t. $x^2 \le 0$: the only feasible point is $x^* = 0$, but $\nabla g(0) = 0$, so $1 + \mu \cdot 0 = 0$ has no solution. The minimum exists and no KKT point does.
- (b) The Kuhn-Tucker cusp. At $(1, 0)$ the two active gradients $(0, 1)$ and $(0, -1)$ point in opposite directions and cannot balance $\nabla f = (-1, 0)$.
- (c) The same feasible set as (a), written with two linear constraints. LICQ still fails (two active constraints in one dimension), but the linearity CQ holds and KKT multipliers exist. Rewriting the model fixed the problem.
- (d) The running example satisfies LICQ, so its multipliers are unique.
Key takeaways
- KKT conditions are necessary only under a constraint qualification: the linearized constraints must describe the true local geometry.
- LICQ gives unique multipliers. MFCQ gives a bounded set. Slater is the standard CQ for convex problems. Linear constraints need none.
- Failures come from cusps, single-point feasible sets, squared constraints and complementarity constraints.
- When a solver struggles, check whether the constraints can be rewritten (often linearly) to restore regularity.
Practice
- (Beginner) Check LICQ at $x^* = (1, 1)$ for the running example by hand: list the active constraints and test their gradients for linear independence.
- (Intermediate) For $\min x_1 + x_2$ s.t. $(x_1 - 1)^2 + x_2^2 \le 1$ and $(x_1 + 1)^2 + x_2^2 \le 1$, find the feasible set and the optimum. Show that no KKT multipliers exist, and explain which CQ fails and why Slater's condition does not hold.
- (Interview) "An NLP solver returns 'converged to a point of local infeasibility' or reports multipliers of order $10^{8}$ on a model you know is feasible. What do you look for in the formulation?" Hint: degenerate or redundant constraints, squared equalities, complementarity-type products, and missing constraint qualifications at the solution.
Penalty and Barrier Methods
Why it matters. Constrained problems are harder than unconstrained ones, and Section 5.1 already provides excellent unconstrained algorithms. Penalty and barrier methods convert a constrained problem into a sequence of unconstrained ones. The barrier idea, followed carefully, is the interior point method that solves LPs (Section 2.6) and most large NLPs today (IPOPT, KNITRO).
Theory. Both methods replace the hard walls of the feasible region by a cost.
Quadratic penalty (exterior). Charge for constraint violation:
$$ P(x; \rho) = f(x) + \frac{\rho}{2} \sum_i \max{0,\ g_i(x)}^2 + \frac{\rho}{2} \sum_j h_j(x)^2 . $$
Minimize $P$ for an increasing sequence $\rho_1 < \rho_2 < \dots$, warm-starting each solve from the previous solution. The minimizers $x(\rho)$ are usually infeasible and approach $x^$ from outside. The violation shrinks like $\mu^ / \rho$, and $\rho \max{0, g_i(x(\rho))}$ estimates the multiplier $\mu_i^*$.
Log barrier (interior). Forbid leaving the interior:
$$ B(x; \tau) = f(x) - \tau \sum_i \ln\big(-g_i(x)\big), \qquad \text{defined only where } g_i(x) < 0 . $$
Minimize $B$ for a decreasing sequence $\tau \to 0$. The barrier is infinite on the boundary, so the iterates stay strictly feasible, and the minimizers $x(\tau)$ trace the central path toward $x^$. Stationarity of $B$ gives $\nabla f + \sum_i \frac{\tau}{-g_i} \nabla g_i = 0$, so $\tau / (-g_i(x(\tau)))$ estimates $\mu_i^$. The barrier solution satisfies the perturbed complementarity $\mu_i g_i = -\tau$ instead of $\mu_i g_i = 0$. For convex problems the suboptimality is at most $m \tau$.
extra cost for the constraint x <= 1
quadratic penalty (rho/2) max(0, x-1)^2 log barrier -tau ln(1 - x)
| / | |
| / larger rho: | | wall at x = 1
| / steeper | / smaller tau:
| / | .' wall hugs x = 1
|______________________/ |__________....---''
+----------------------+----- x +----------------------+--- x
1 1
zero inside, grows outside finite inside, infinite at the wall
| Quadratic penalty | Log barrier | |
|---|---|---|
| Iterates | usually infeasible, approach from outside | strictly feasible, approach from inside |
| Parameter | $\rho \to \infty$ | $\tau \to 0$ |
| Equality constraints | handled directly | cannot be handled by the barrier (combine with a penalty or Newton-KKT steps) |
| Needs a feasible start | no | yes, a strictly feasible one |
| Weakness | Hessian condition number grows like $\rho$ | the same as $\tau \to 0$, and line searches must stay inside |
| Modern descendant | augmented Lagrangian (next topic) | primal-dual interior point methods |
Exact penalties. The nonsmooth penalty $f(x) + \rho \sum_i \max{0, g_i(x)} + \rho \sum_j \lvert h_j(x) \rvert$ recovers $x^*$ exactly for any finite $\rho > \max$ multiplier, with no ill-conditioning, but it has kinks that smooth methods dislike. It is mainly used as a merit function that measures progress inside SQP.
Example: both methods on the running example. The penalty uses BFGS for each subproblem. The barrier uses Nelder-Mead, which tolerates the $+\infty$ values outside the domain:
import numpy as np
from scipy.optimize import minimize
f = lambda x: (x[0] - 2)**2 + (x[1] - 1)**2
g = lambda x: np.array([x[0]**2 - x[1], x[0] + x[1] - 2]) # g(x) <= 0
print("Quadratic penalty: min f + (rho/2) * sum max(0, g)^2 (approaches from OUTSIDE)")
x = np.array([2.0, 1.0]) # unconstrained minimizer: infeasible
for rho in [1, 10, 100, 1000, 10000]:
P = lambda x: f(x) + rho / 2 * np.sum(np.maximum(0, g(x))**2)
x = minimize(P, x, method="BFGS", options={"gtol": 1e-10}).x # warm start
viol = np.maximum(0, g(x))
print(f" rho = {rho:5d}: x = ({x[0]:.5f}, {x[1]:.5f}) max violation = {viol.max():.1e}"
f" multiplier estimate rho*g+ = {(rho * viol).round(4)}")
print("Log barrier: min f - tau * sum log(-g) (approaches from INSIDE)")
def B(x, tau):
gx = g(x)
return np.inf if np.any(gx >= 0) else f(x) - tau * np.sum(np.log(-gx)) # +inf outside
x = np.array([0.0, 0.5]) # strictly feasible start
for tau in [1, 0.1, 0.01, 0.001, 0.0001]:
x = minimize(B, x, args=(tau,), method="Nelder-Mead", # derivative-free, copes with +inf
options={"xatol": 1e-10, "fatol": 1e-14, "maxiter": 10000}).x
print(f" tau = {tau:6.4f}: x = ({x[0]:.5f}, {x[1]:.5f}) slack = {(-g(x)).round(5)}"
f" multiplier estimate tau/(-g) = {(tau / -g(x)).round(4)}")
# Quadratic penalty: min f + (rho/2) * sum max(0, g)^2 (approaches from OUTSIDE)
# rho = 1: x = (1.24594, 1.07661) max violation = 4.8e-01 multiplier estimate rho*g+ = [0.4758 0.3226]
# rho = 10: x = (1.04053, 1.01917) max violation = 6.4e-02 multiplier estimate rho*g+ = [0.6353 0.5969]
# rho = 100: x = (1.00440, 1.00219) max violation = 6.6e-03 multiplier estimate rho*g+ = [0.6632 0.6589]
# rho = 1000: x = (1.00044, 1.00022) max violation = 6.7e-04 multiplier estimate rho*g+ = [0.6663 0.6659]
# rho = 10000: x = (1.00004, 1.00002) max violation = 6.7e-05 multiplier estimate rho*g+ = [0.6666 0.6664]
# Log barrier: min f - tau * sum log(-g) (approaches from INSIDE)
# tau = 1.0000: x = (0.47863, 0.91189) slack = [0.6828 0.60947] multiplier estimate tau/(-g) = [1.4646 1.6408]
# tau = 0.1000: x = (0.91186, 0.96549) slack = [0.134 0.12265] multiplier estimate tau/(-g) = [0.7463 0.8153]
# tau = 0.0100: x = (0.99015, 0.99523) slack = [0.01483 0.01462] multiplier estimate tau/(-g) = [0.6745 0.684 ]
# tau = 0.0010: x = (0.99900, 0.99950) slack = [0.0015 0.0015] multiplier estimate tau/(-g) = [0.6674 0.6684]
# tau = 0.0001: x = (0.99990, 0.99995) slack = [0.00015 0.00015] multiplier estimate tau/(-g) = [0.6667 0.6668]
The penalty iterates have $x_1 > 1$ (outside the parabola) and the barrier iterates have $x_1 < 1$ (inside). Both multiplier estimates converge to $\tfrac{2}{3}$. The price is that each tenfold improvement in accuracy needs a tenfold change in $\rho$ or $\tau$, which makes the subproblems increasingly ill-conditioned. Practical interior point methods avoid this with primal-dual Newton steps and only a few dozen iterations in total (Section 2.6).
Key takeaways
- Penalty methods add a cost for violation and approach the optimum from outside as $\rho \to \infty$.
- Barrier methods add an infinite wall at the boundary and approach from inside as $\tau \to 0$, following the central path.
- Both give multiplier estimates for free: $\rho g_i^+$ and $\tau / (-g_i)$.
- Both become ill-conditioned in the limit. The augmented Lagrangian and primal-dual interior point methods are the practical fixes.
Practice
- (Beginner) Solve $\min x^2$ s.t. $x \ge 1$ with the quadratic penalty in closed form: show that $x(\rho) = \rho / (2 + \rho)$. How large must $\rho$ be for a violation below $10^{-3}$?
- (Intermediate) For the same problem, solve the barrier subproblem $\min x^2 - \tau \ln(x - 1)$ in closed form and check that $x(\tau) - 1 \approx \tau / 2$ for small $\tau$. Compute the multiplier estimate $\tau / (x(\tau) - 1)$ and its limit.
- (Interview) "Why don't we just set $\rho = 10^{12}$ in a single penalty solve instead of increasing it gradually?" Hint: look at the condition number of $\nabla^2 P$ and what it does to BFGS or Newton steps, and at the value of warm starts.
Augmented Lagrangian Methods
Why it matters. The augmented Lagrangian (method of multipliers) keeps the simplicity of the penalty method while removing its fatal flaw: it reaches the exact solution with a moderate, fixed penalty parameter. It underlies classical NLP codes (LANCELOT, ALGENCAN), many large-scale and GPU solvers, and the ADMM algorithm used for distributed optimization.
Theory. Start with equality constraints. The augmented Lagrangian is the Lagrangian plus a quadratic penalty:
$$ L_\rho(x, \lambda) = f(x) + \lambda^\top h(x) + \frac{\rho}{2} \lVert h(x) \rVert^2 . $$
The method of multipliers alternates two steps:
- $x_{k+1} = \arg\min_x L_\rho(x, \lambda_k)$ (an unconstrained problem).
- $\lambda_{k+1} = \lambda_k + \rho, h(x_{k+1})$ (a price update).
Why it works: at $\lambda = \lambda^$, the term $\lambda^{\top} h$ already makes $x^$ a stationary point, and the penalty only adds curvature that makes $x^$ a strict local minimum of $L_\rho$ for every $\rho$ above some finite threshold. So the penalty parameter does not need to go to infinity. It only needs to be large enough while the multipliers do the rest. Analogy: the pure penalty method is a stick that has to get bigger and bigger. The augmented Lagrangian is a carrot and a stick. A moderate stick, plus a price that keeps adjusting until the constraint is respected voluntarily.
The update in step 2 is a gradient ascent step on the dual function with step size $\rho$ (Section 5.3), the same price-adjustment logic as Lagrangian relaxation (Section 4.2).
Inequalities use the Powell-Hestenes-Rockafellar (PHR) form:
$$ L_\rho(x, \mu) = f(x) + \frac{1}{2\rho} \sum_i \Big( \max{0,\ \mu_i + \rho, g_i(x)}^2 - \mu_i^2 \Big), \qquad \mu_i \leftarrow \max{0,\ \mu_i + \rho, g_i(x)}. $$
A constraint that is comfortably satisfied makes its term vanish and its multiplier drop to 0, which is complementary slackness emerging automatically.
ADMM. For $\min f(x) + g(z)$ s.t. $Ax + Bz = c$, the alternating direction method of multipliers minimizes $L_\rho$ over $x$, then over $z$, then updates the multiplier. Each block is often simple or decomposes into independent pieces that can run in parallel. That is why ADMM is popular for large distributed problems and for first-order QP solvers such as OSQP.
Example: the running example with a fixed $\rho = 10$. With the pure quadratic penalty, $\rho = 10$ left a violation of $6.4 \times 10^{-2}$. The augmented Lagrangian with the same $\rho$:
import numpy as np
from scipy.optimize import minimize
f = lambda x: (x[0] - 2)**2 + (x[1] - 1)**2
g = lambda x: np.array([x[0]**2 - x[1], x[0] + x[1] - 2]) # g(x) <= 0
def aug_lagrangian(x, mu, rho):
# PHR form for inequalities: f + 1/(2 rho) * sum( max(0, mu + rho g)^2 - mu^2 )
return f(x) + np.sum(np.maximum(0, mu + rho * g(x))**2 - mu**2) / (2 * rho)
x, mu, rho = np.array([2.0, 1.0]), np.zeros(2), 10.0 # rho stays FIXED at a moderate value
for k in range(1, 8):
x = minimize(aug_lagrangian, x, args=(mu, rho), method="BFGS", options={"gtol": 1e-12}).x
mu = np.maximum(0, mu + rho * g(x)) # multiplier (price) update
print(f"outer {k}: x = ({x[0]:.8f}, {x[1]:.8f}) mu = {mu.round(6)}"
f" max violation = {np.maximum(0, g(x)).max():.1e}")
# outer 1: x = (1.04052570, 1.01916728) mu = [0.635265 0.59693 ] max violation = 6.4e-02
# outer 2: x = (1.00307100, 1.00314268) mu = [0.665352 0.659067] max violation = 6.2e-03
# outer 3: x = (1.00026819, 1.00040764) mu = [0.66664 0.665825] max violation = 6.8e-04
# outer 4: x = (1.00002592, 1.00004884) mu = [0.66667 0.666573] max violation = 7.5e-05
# outer 5: x = (1.00000268, 1.00000566) mu = [0.666667 0.666656] max violation = 8.3e-06
# outer 6: x = (1.00000029, 1.00000064) mu = [0.666667 0.666665] max violation = 9.3e-07
# outer 7: x = (1.00000003, 1.00000007) mu = [0.666667 0.666666] max violation = 1.0e-07
The first outer iteration ($\mu = 0$) is exactly the penalty solve with $\rho = 10$. After that, every multiplier update cuts the violation by about a factor of 10, and $\mu$ converges to $(\tfrac{2}{3}, \tfrac{2}{3})$ while $\rho$ never changes. Each subproblem stays as well-conditioned as the first.
Key takeaways
- The augmented Lagrangian adds a quadratic penalty to the Lagrangian and updates multipliers with $\lambda \leftarrow \lambda + \rho h(x)$.
- Exact solutions are reached with a finite, moderate $\rho$, which avoids the ill-conditioning of pure penalties.
- The multiplier update is dual ascent, the same price-adjustment idea as Lagrangian relaxation.
- ADMM splits the minimization into blocks and is widely used for large, distributed and first-order methods.
Practice
- (Beginner) For $\min x^2$ s.t. $x = 1$, write $L_\rho(x, \lambda)$, minimize it in closed form, and do three multiplier updates by hand with $\rho = 1$ from $\lambda_0 = 0$. What are $x^$ and $\lambda^$?
- (Intermediate) Rerun the example with $\rho = 1$ and with $\rho = 100$. Count the outer iterations needed for a violation below $10^{-7}$ and the BFGS iterations of each inner solve. Explain the trade-off.
- (Interview) "Explain ADMM to a colleague who knows Lagrangian relaxation." Hint: both price the coupling constraints. ADMM adds a quadratic term so the price updates converge reliably, and it alternates between blocks instead of solving them jointly.
Sequential Quadratic Programming
Why it matters. SQP is Newton's method for constrained optimization. It converges in few iterations, it handles equalities and inequalities naturally, and it is the method of choice for small and medium smooth NLPs with expensive functions (SNOPT, filterSQP, KNITRO's active-set mode, MATLAB's fmincon with the sqp option, SciPy's SLSQP).
Theory. Newton's method in Section 5.1 minimizes a quadratic model of $f$ at each step. SQP does the same for constrained problems: at the current point $(x_k, \mu_k, \lambda_k)$, solve the quadratic programming subproblem
$$ \begin{aligned} \min_d \quad & \nabla f(x_k)^\top d + \tfrac{1}{2}, d^\top B_k, d \ \text{s.t.} \quad & g_i(x_k) + \nabla g_i(x_k)^\top d \le 0, \qquad h_j(x_k) + \nabla h_j(x_k)^\top d = 0, \end{aligned} $$
with $B_k = \nabla^2_{xx} L(x_k, \mu_k, \lambda_k)$ or an approximation. Then set $x_{k+1} = x_k + d$ and take the new multipliers from the QP's own multipliers.
Two details matter:
- Linearized constraints, curved objective. The constraints are replaced by their tangent planes, so their curvature must be captured somewhere. It goes into $B_k$, which is the Hessian of the Lagrangian, not of $f$. In the running example, $\nabla^2 L = 2I + \mu_1 \nabla^2 g_1$, and the term $\mu_1 \nabla^2 g_1$ carries the curvature of the parabola. Using $\nabla^2 f$ alone would converge slowly or not at all.
- Newton on KKT. With equality constraints only, the QP's optimality conditions are the linear system
$$ \begin{pmatrix} B_k & \nabla h(x_k)^\top \ \nabla h(x_k) & 0 \end{pmatrix} \begin{pmatrix} d \ \lambda_{k+1} \end{pmatrix} = -\begin{pmatrix} \nabla f(x_k) \ h(x_k) \end{pmatrix}, $$
which is exactly a Newton step on the KKT equations. That is why SQP inherits Newton's quadratic local convergence. With inequalities, the QP also picks the active set, and near the solution it identifies the correct one.
Making it robust (globalization):
| Ingredient | Purpose |
|---|---|
| Merit function $\phi(x) = f(x) + \rho \sum_i \max{0, g_i} + \rho \sum_j \lvert h_j \rvert$ with a line search | decide whether a step is progress when $f$ and feasibility trade off |
| Filter methods | accept a step if it improves either $f$ or the constraint violation, without a penalty parameter |
| Damped BFGS for $B_k$ | avoid second derivatives and keep the QP convex |
| Second-order correction | avoid the Maratos effect, where a good step is rejected by the merit function near the solution |
SQP is an infeasible-path method: intermediate iterates may violate the constraints, unlike barrier methods. If the model cannot be evaluated outside the feasible region (a square root of a negative number, say), add bounds or use an interior method.
Example: SQP from $(0, 0)$, with cvxpy solving each QP subproblem. Full steps are taken, without a merit function, to show the raw local behavior. Then SciPy's production code solves the same problem.
import numpy as np
import cvxpy as cp
from scipy.optimize import minimize
f = lambda x: (x[0] - 2)**2 + (x[1] - 1)**2
f_grad = lambda x: np.array([2 * (x[0] - 2), 2 * (x[1] - 1)])
g = lambda x: np.array([x[0]**2 - x[1], x[0] + x[1] - 2])
g_jac = lambda x: np.array([[2 * x[0], -1.0], [1.0, 1.0]])
def hess_lagrangian(x, mu): # grad^2 f + sum mu_i grad^2 g_i
return 2 * np.eye(2) + mu[0] * np.array([[2.0, 0.0], [0.0, 0.0]])
x, mu = np.array([0.0, 0.0]), np.zeros(2)
for k in range(1, 7):
d = cp.Variable(2)
lin = g(x) + g_jac(x) @ d <= 0 # linearized constraints
qp = cp.Problem(cp.Minimize(f_grad(x) @ d + 0.5 * cp.quad_form(d, hess_lagrangian(x, mu))),
[lin])
qp.solve()
x, mu = x + d.value, lin.dual_value # full step; QP multipliers -> new mu
print(f"SQP iter {k}: x = ({x[0]:.10f}, {x[1]:.10f}) mu = {mu.round(8)}"
f" |x - x*| = {np.linalg.norm(x - 1):.1e}")
# Production code: scipy's SLSQP (expects constraints as c(x) >= 0)
res = minimize(f, [0.0, 0.0], jac=f_grad, method="SLSQP",
constraints={"type": "ineq", "fun": lambda x: -g(x), "jac": lambda x: -g_jac(x)})
print(f"SLSQP: x = {res.x.round(8)}, f = {res.fun:.6f}, iterations = {res.nit}")
# SQP iter 1: x = (1.5000000000, 0.5000000000) mu = [0. 1.] |x - x*| = 7.1e-01 <- infeasible intermediate iterate
# SQP iter 2: x = (1.0625000000, 0.9375000000) mu = [0.4375 0.5625] |x - x*| = 8.8e-02
# SQP iter 3: x = (1.0012500000, 0.9987500000) mu = [0.65555 0.65805] |x - x*| = 1.8e-03
# SQP iter 4: x = (1.0000005204, 0.9999994796) mu = [0.66665649 0.66665753] |x - x*| = 7.4e-07
# SQP iter 5: x = (1.0000000000, 1.0000000000) mu = [0.66666667 0.66666667] |x - x*| = 1.3e-13
# SQP iter 6: x = (1.0000000000, 1.0000000000) mu = [0.66666667 0.66666667] |x - x*| = 1.1e-16
# SLSQP: x = [1. 1.], f = 1.000000, iterations = 6
The first QP linearizes the parabola at $x_1 = 0$, where it looks flat, and lands at the infeasible point $(1.5, 0.5)$. From then on the error is roughly squared at each iteration ($10^{-1} \to 10^{-3} \to 10^{-7} \to 10^{-13}$), and the multipliers converge with $x$. Compare the six SQP iterations with the five penalty solves or seven augmented Lagrangian solves, each of which ran a full BFGS minimization.
Choosing a constrained NLP method:
| Method | Strengths | Weaknesses | Representative codes |
|---|---|---|---|
| Penalty / barrier (classical) | simple, reuses unconstrained solvers | ill-conditioned, slow | teaching, quick prototypes |
| Interior point (primal-dual barrier) | very large sparse problems, polynomial for convex | hard to warm-start | IPOPT, KNITRO, MOSEK |
| Augmented Lagrangian | moderate penalties, parallelizes, first-order friendly | linear convergence | LANCELOT, ALGENCAN, OSQP (via ADMM) |
| SQP | few iterations, excellent warm starts, active set identified | each iteration solves a QP, which is expensive for very large problems | SNOPT, filterSQP, SLSQP, fmincon (sqp) |
Key takeaways
- SQP solves a QP at each iteration: a quadratic model of the Lagrangian with linearized constraints.
- The QP Hessian must be the Hessian of the Lagrangian, so that constraint curvature is represented.
- With equalities, an SQP step is a Newton step on the KKT system, so local convergence is quadratic (superlinear with BFGS).
- Merit functions or filters, with a line search or trust region, make SQP converge from distant starting points.
Practice
- (Beginner) Write the first SQP subproblem for the running example at $x_0 = (0, 0)$ with $\mu_0 = 0$, and solve it by hand. Confirm that the step lands at $(1.5, 0.5)$.
- (Intermediate) Replace $B_k$ in the example by $\nabla^2 f = 2I$ (ignoring constraint curvature) and rerun. How do the iteration count and the convergence pattern change? Explain using the Hessian of the Lagrangian.
- (Interview) "You have a model with 300 variables whose objective requires a 2-second simulation per evaluation. Would you choose an interior point method or SQP?" Hint: count function evaluations per iteration and in total, consider warm starts when the model is re-solved, and how derivatives are obtained.
<a id="53-convex-optimization"></a>
5.3 Convex Optimization
This subchapter studies the problems that keep the good properties of LP: every local optimum is global, duality gives bounds and certificates, and interior point methods solve them in polynomial time with a few dozen Newton steps. It covers convex duality and then the standard problem classes, each of which strictly contains the previous one:
convex programs any convex objective and constraints
+-- SDP matrix variable X, constraint X positive semidefinite
+-- SOCP constraints ||A x + b||_2 <= c'x + d
+-- convex QCQP convex quadratic objective and constraints
+-- convex QP convex quadratic objective, linear constraints
+-- LP linear objective, linear constraints (Chapter 2)
Modeling tools such as cvxpy, CVX, JuMP and YALMIP implement disciplined convex programming: if a model is built from known convex atoms combined with convexity-preserving rules (Section 1.2), the tool verifies convexity, identifies the class and passes the model to a suitable conic solver (MOSEK, Clarabel, SCS, ECOS, Gurobi). The examples below use cvxpy.
Convex Programs and Duality
Why it matters. For convex problems, the dual problem provides a lower bound that matches the optimum exactly (under a mild condition), so every solution comes with a certificate of optimality. Dual variables are prices, and dualizing a coupling constraint splits a large problem into small independent ones (dual decomposition), which is how electricity markets, network bandwidth allocation and many distributed algorithms work.
Theory. A convex program has the form
$$ \min\ f(x) \quad \text{s.t.} \quad g_i(x) \le 0 \ \ (i = 1, \dots, m), \qquad Ax = b, $$
with $f$ and all $g_i$ convex. Equalities must be affine: $h(x) = 0$ with nonlinear $h$ generally defines a non-convex set (a circle, for example). Its feasible set is convex, and every local minimum is global (Section 1.3).
The Lagrangian is $L(x, \mu, \lambda) = f(x) + \mu^\top g(x) + \lambda^\top (Ax - b)$, and the dual function is
$$ \theta(\mu, \lambda) = \inf_x\ L(x, \mu, \lambda). $$
Two properties hold for any problem, convex or not:
- $\theta$ is concave, because it is a pointwise infimum of functions that are affine in $(\mu, \lambda)$. So the dual problem $\max_{\mu \ge 0,, \lambda} \theta(\mu, \lambda)$ is always a convex problem, even when the primal is not.
- Weak duality: for every feasible $x$ and every $\mu \ge 0$, $\theta(\mu, \lambda) \le L(x, \mu, \lambda) = f(x) + \mu^\top g(x) \le f(x)$. Every dual-feasible point gives a lower bound on the optimum.
The difference $f^* - \theta^*$ is the duality gap. For non-convex problems it is often positive. The Lagrangian bounds of Section 4.2 stay below the integer optimum for exactly this reason.
Strong duality (Slater's theorem). If the problem is convex and some point is strictly feasible ($g_i(\bar x) < 0$ for all $i$, $A\bar x = b$), then $f^* = \theta^*$ and the dual optimum is attained. Then:
- $x^$ and $(\mu^, \lambda^*)$ are optimal if and only if they satisfy the KKT conditions.
- $\mu^$ and $\lambda^$ are the sensitivities of $f^*$ to the constraints, as in LP.
- For affine constraints, strict feasibility is not needed: plain feasibility suffices. LP strong duality is the special case where everything is affine.
Example in one line. For $\min x^2$ s.t. $1 - x \le 0$: $L = x^2 + \mu(1 - x)$ is minimized at $x = \mu/2$, so $\theta(\mu) = \mu - \mu^2/4$. It is maximized at $\mu^* = 2$ with $\theta^* = 1 = f^*$. Every other $\mu \ge 0$ gives a lower bound.
theta(mu) = mu - mu^2/4 f* = 1 (primal optimum, x* = 1)
1 +- - - - - - - * - - - - - - - <- the dual optimum touches f*: zero gap
| .' '.
| .' '. every theta(mu) is a lower bound
| .' '. (weak duality)
0 *----+------+------+-----*-----> mu
0 1 2 3 4
Dual decomposition. If $f$ and the constraints are separable except for a few coupling constraints, dualizing only the coupling constraints splits the inner minimization into independent subproblems, one per block. A coordinator adjusts the multipliers (prices) until the coupling constraints hold. This is the convex counterpart of Lagrangian relaxation.
Example: economic dispatch. Three generators with costs $a_i p_i + b_i p_i^2$ and capacities must meet a demand of 400 MW. Dualizing the demand constraint with price $\lambda$ gives each generator a private problem: produce where its marginal cost $a_i + 2b_i p_i$ equals the price, within its capacity. Bisection on $\lambda$ finds the price at which supply equals demand:
import numpy as np
import cvxpy as cp
# Economic dispatch: cost_i(p) = a_i p + b_i p^2, 0 <= p_i <= cap_i, sum p = D
a = np.array([10.0, 12.0, 20.0]) # $/MWh linear cost
b = np.array([0.05, 0.02, 0.01]) # $/MWh^2 quadratic cost
cap = np.array([150.0, 200.0, 250.0]) # MW
D = 400.0 # MW demand
def best_response(lam):
"""Each generator, offered price lam, maximizes its own profit lam p - cost(p)."""
return np.clip((lam - a) / (2 * b), 0, cap)
def dual_function(lam): # g(lam) = min_p sum cost(p) - lam (sum p - D)
p = best_response(lam)
return np.sum(a * p + b * p**2) - lam * (p.sum() - D)
lo, hi = 0.0, 100.0 # bisection on the price: sum p(lam) is increasing
for _ in range(60):
lam = (lo + hi) / 2
lo, hi = (lam, hi) if best_response(lam).sum() < D else (lo, lam)
p_dual = best_response(lam)
print(f"dual decomposition: price = {lam:.4f} $/MWh, p = {p_dual.round(2)}, g(lam) = {dual_function(lam):.4f}")
p = cp.Variable(3)
balance = cp.sum(p) == D
prob = cp.Problem(cp.Minimize(a @ p + b @ cp.square(p)), [balance, p >= 0, p <= cap])
prob.solve()
print(f"primal (cvxpy) : cost = {prob.value:.4f} $/h, p = {p.value.round(2)}, "
f"dual of balance = {abs(balance.dual_value):.4f}")
for lam_try in [20, 25, 30]:
print(f" g({lam_try}) = {dual_function(lam_try):9.2f} <= primal optimum (weak duality)")
# dual decomposition: price = 21.6667 $/MWh, p = [116.67 200. 83.33], g(lam) = 6783.3333
# primal (cvxpy) : cost = 6783.3333 $/h, p = [116.67 200. 83.33], dual of balance = 21.6667
# g(20) = 6700.00 <= primal optimum (weak duality)
# g(25) = 6450.00 <= primal optimum (weak duality)
# g(30) = 5450.00 <= primal optimum (weak duality)
The dual decomposition and the centralized solver agree on the cost and the dispatch, and the dual value of the demand constraint is the clearing price $\lambda^* = 21.67$ $/MWh: the cost of serving one more MW. Generator 2 runs at its capacity, where its marginal cost is $12 + 0.04 \cdot 200 = 20 < 21.67$. Its capacity constraint has a multiplier of $1.67$, the value of one more MW of capacity at that generator. Every trial price gives a lower bound, as weak duality says.
Conic form. Modern solvers accept convex problems in the conic form $\min c^\top x$ s.t. $Ax = b$, $x \in \mathcal{K}$, where $\mathcal{K}$ is a product of convex cones:
| Cone $\mathcal{K}$ | Members | Problem class |
|---|---|---|
| Nonnegative orthant | $x \ge 0$ | LP |
| Second-order (Lorentz) cone | $\lVert u \rVert_2 \le t$ | SOCP (including convex QP and QCQP) |
| Positive semidefinite cone | symmetric $X \succeq 0$ | SDP |
| Exponential and power cones | $y e^{x/y} \le z$, $x^\alpha y^{1-\alpha} \ge \lvert z \rvert$ | entropy, log-sum-exp, geometric programming |
Conic duality mirrors LP duality, and primal-dual interior point methods (Section 2.6) extend to all these cones.
Key takeaways
- Convex program: convex objective, convex inequality constraints, affine equalities. Local optima are global.
- The dual function is always concave and gives lower bounds (weak duality), for convex and non-convex problems alike.
- Under Slater's condition, strong duality holds, KKT conditions are necessary and sufficient, and multipliers are sensitivities.
- Dualizing coupling constraints splits a problem into independent subproblems coordinated by prices.
Practice
- (Beginner) Derive the dual function of $\min x_1^2 + x_2^2$ s.t. $x_1 + x_2 \ge 2$, maximize it, and check that the optimal values match. Which point is strictly feasible?
- (Intermediate) In the dispatch example, raise the demand to 500 MW. Predict which generators will be at capacity, then rerun both the bisection and cvxpy. How much does the clearing price rise?
- (Interview) "How can a solver prove that a solution is optimal without searching the whole feasible region?" Hint: a dual-feasible point gives a lower bound. Discuss duality gaps, and when they are zero.
Quadratic Programming
Why it matters. Quadratic programs are the first step beyond LP and appear wherever cost or risk grows with the square of a decision: portfolio variance, least-squares fitting with constraints, model predictive control, support vector machines, and the subproblems of SQP. Convex QPs can be solved almost as fast as LPs of the same size.
Theory. A QP has a quadratic objective and linear constraints:
$$ \min\ \tfrac{1}{2} x^\top Q x + c^\top x \quad \text{s.t.} \quad Ax \le b, \quad Ex = d . $$
| $Q$ | Problem | Difficulty |
|---|---|---|
| positive semidefinite | convex QP | polynomial time: local = global |
| positive definite | strictly convex QP | unique optimum |
| indefinite | non-convex QP | NP-hard. Local optima may be many |
Where the optimum lies. Unlike LP, the optimum of a QP does not need to be at a vertex. It can be in the interior of an edge or face, or strictly inside the region if the unconstrained minimum is feasible. The simplex method's vertex search therefore does not apply directly.
Solution methods.
- Equality constraints only: the KKT conditions are a linear system, $\begin{pmatrix} Q & E^\top \ E & 0 \end{pmatrix} \begin{pmatrix} x \ \lambda \end{pmatrix} = \begin{pmatrix} -c \ d \end{pmatrix}$, solved in one step.
- Active-set methods: guess which inequalities are active (the working set), treat them as equalities, and solve the KKT system. If a multiplier is negative, drop that constraint. If the step hits a new constraint, add it. This is the QP analogue of the simplex method. It is excellent for small and medium problems and for warm starts (as in SQP).
- Interior point methods: the same primal-dual barrier approach as for LP (Section 2.6). The best choice for large sparse QPs (Gurobi, CPLEX, MOSEK).
- First-order methods: ADMM-based solvers such as OSQP, which are fast for moderate accuracy and embedded applications.
Example: Markowitz portfolio selection. Choose weights $w$ for four assets with expected returns $\bar r$ and covariance matrix $\Sigma$ to minimize the variance $w^\top \Sigma w$ while reaching a target return $r$:
$$ \min\ w^\top \Sigma w \quad \text{s.t.} \quad \mathbf{1}^\top w = 1, \quad \bar r^\top w \ge r, \quad w \ge 0 . $$
Without the sign constraints, the KKT conditions are a linear system. With them, a QP solver must find which assets are left out (the active set).
import numpy as np
import cvxpy as cp
rbar = np.array([0.04, 0.07, 0.10, 0.12]) # expected annual returns
vol = np.array([0.05, 0.12, 0.18, 0.25]) # volatilities
corr = np.array([[1.0, 0.2, 0.1, 0.0],
[0.2, 1.0, 0.5, 0.3],
[0.1, 0.5, 1.0, 0.6],
[0.0, 0.3, 0.6, 1.0]])
Sigma = np.outer(vol, vol) * corr # covariance matrix (PSD)
# (1) Equality constraints only: the KKT conditions are a LINEAR system
# min 1/2 w'(2 Sigma)w s.t. 1'w = 1, rbar'w = r
r = 0.10
A = np.vstack([np.ones(4), rbar]); rhs = np.array([1.0, r])
KKT = np.block([[2 * Sigma, A.T], [A, np.zeros((2, 2))]])
sol = np.linalg.solve(KKT, np.concatenate([np.zeros(4), rhs]))
w_eq = sol[:4]
print(f"KKT solve (shorting allowed): w = {w_eq.round(4)}, risk = {np.sqrt(w_eq @ Sigma @ w_eq):.4f}")
# (2) Long-only: add w >= 0 and let a QP solver find the active set
w = cp.Variable(4)
ret = rbar @ w >= r
prob = cp.Problem(cp.Minimize(cp.quad_form(w, Sigma)), [cp.sum(w) == 1, ret, w >= 0])
prob.solve()
print(f"QP solve (long only) : w = {w.value.round(4)}, risk = {np.sqrt(prob.value):.4f}")
# (3) Efficient frontier: re-solve for several return targets
print("target risk weights")
for r_t in [0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.11]:
prob = cp.Problem(cp.Minimize(cp.quad_form(w, Sigma)), [cp.sum(w) == 1, rbar @ w >= r_t, w >= 0])
prob.solve()
print(f" {r_t:.2f} {np.sqrt(prob.value):.4f} {np.round(w.value, 3) + 0.0}")
# KKT solve (shorting allowed): w = [-0.1091 0.4225 0.3803 0.3064], risk = 0.1583
# QP solve (long only) : w = [0. 0.232 0.42 0.348], risk = 0.1597
# target risk weights
# 0.05 0.0508 [0.799 0.105 0.04 0.055]
# 0.06 0.0644 [0.617 0.169 0.108 0.105]
# 0.07 0.0847 [0.436 0.232 0.176 0.156]
# 0.08 0.1080 [0.254 0.296 0.244 0.206]
# 0.09 0.1328 [0.073 0.359 0.312 0.256]
# 0.10 0.1597 [0. 0.232 0.42 0.348]
# 0.11 0.1929 [0. 0. 0.5 0.5]
At a 10% target, the unconstrained KKT solution shorts asset 1 (negative weight). Forbidding short sales makes $w_1 \ge 0$ active, and the risk rises from 15.83% to 15.97%. Re-solving for several targets traces the efficient frontier: the minimum risk for each return level. Asset 1 leaves the portfolio at 10%, and asset 2 at 11%.
return
12% | xA4 * efficient portfolios (QP optima)
11% | * x single assets
10% | * xA3
9% | * diversification: every frontier
8% | * point has lower risk than the
7% | * xA2 single asset with the same return
6% | *
5% | *
4% | xA1
+--------------------------------------------> risk (st. dev.)
5% 10% 15% 20% 25%
Key takeaways
- A QP has a quadratic objective and linear constraints. It is convex, and efficiently solvable, when $Q$ is positive semidefinite.
- The optimum can lie inside a face, not only at a vertex.
- Equality-constrained QPs reduce to one linear KKT system. Inequalities are handled by active-set, interior point or first-order (ADMM) methods.
- Markowitz portfolio selection is the classical QP: minimizing variance for each return target traces the efficient frontier.
Practice
- (Beginner) Solve $\min x_1^2 + x_2^2$ s.t. $x_1 + x_2 = 1$ with the KKT linear system. Then add $x_1 \le 0.2$ and determine the new active set by reasoning.
- (Intermediate) Add a constraint to the portfolio model limiting any single asset to at most 40% of the portfolio. Recompute the frontier, and identify for which targets the new constraint is active.
- (Interview) "Why can't we solve a QP with the simplex method directly, and why do QP solvers still feel like LP solvers?" Hint: optima off the vertices, KKT conditions as a linear complementarity problem, and active-set methods as simplex-like pivoting on the KKT system.
Second-Order Cone Programming
Why it matters. Euclidean distances, norms, robust constraints and some probabilistic constraints are not linear or quadratic in the LP/QP sense, yet they are convex. Second-order cone programming handles all of them, and modern interior point solvers solve SOCPs almost as fast as LPs. Many problems that look non-linear and hard turn out to be SOCPs.
Theory. The second-order cone (Lorentz cone, or "ice-cream cone") in $\mathbb{R}^{k+1}$ is
$$ \mathcal{Q}^{k+1} = { (t, u) \in \mathbb{R} \times \mathbb{R}^k : \lVert u \rVert_2 \le t }. $$
t
^
\ | / the cone ||u|| <= t: every horizontal slice at height t
\ | / is a disc of radius t
\ | /
\ | / an SOC constraint ||A x + b|| <= c'x + d says that an
\|/ affine image of x lies inside this cone
------+------> u
An SOCP minimizes a linear objective subject to constraints of the form
$$ \lVert A_i x + b_i \rVert_2 \le c_i^\top x + d_i . $$
With $A_i = 0$ this is a linear constraint, so LP is a special case. Convex quadratic constraints are SOC constraints too: $x^\top Q x \le t$ with $Q = L^\top L$ is the same as $\lVert (2Lx,\ 1 - t) \rVert \le 1 + t$, which is how QPs and QCQPs fit inside SOCP.
Modeling patterns:
| Model element | SOC form |
|---|---|
| Euclidean distance $\lVert x - a \rVert$ in the objective | $\min t$ s.t. $\lVert x - a \rVert \le t$ |
| Sum of distances (location problems) | one SOC constraint per term |
| Hyperbolic constraint $z^2 \le xy$, $x, y \ge 0$ | $\lVert (2z,\ x - y) \rVert \le x + y$ |
| Robust LP: $a^\top x \le b$ for all $a$ in the ellipsoid ${\bar a + Pu : \lVert u \rVert \le 1}$ | $\bar a^\top x + \lVert P^\top x \rVert \le b$ |
| Chance constraint $\Pr(a^\top x \le b) \ge 1 - \varepsilon$, $a \sim N(\bar a, \Sigma)$, $\varepsilon < 0.5$ | $\bar a^\top x + \Phi^{-1}(1 - \varepsilon) \lVert \Sigma^{1/2} x \rVert \le b$ |
The robust row deserves a word: the worst case of $a^\top x$ over the ellipsoid is $\bar a^\top x + \max_{\lVert u \rVert \le 1} u^\top P^\top x = \bar a^\top x + \lVert P^\top x \rVert$. An infinite family of linear constraints collapses into a single SOC constraint (Section 9.4 develops robust optimization further).
Example: warehouse location (Weber problem) and a robust plan. A warehouse serves five stores with weekly truckloads $w_i$, and transport cost is proportional to straight-line distance. Minimizing $\sum_i w_i \lVert x - a_i \rVert$ is an SOCP. Weiszfeld's classical fixed-point iteration cross-checks the answer. The second part compares a nominal LP plan with a robust plan when the resource coefficients are only known to lie in an ellipsoid.
import numpy as np
import cvxpy as cp
# Weber problem: place one warehouse to minimize weighted straight-line transport
stores = np.array([[0.0, 0.0], [10.0, 0.0], [10.0, 8.0], [2.0, 9.0], [5.0, 3.0]])
demand = np.array([30.0, 10.0, 25.0, 15.0, 5.0]) # truckloads per week
x = cp.Variable(2)
t = cp.Variable(5) # t_i >= ||x - a_i|| (one SOC each)
cons = [cp.SOC(t[i], x - stores[i]) for i in range(5)]
prob = cp.Problem(cp.Minimize(demand @ t), cons)
prob.solve()
print(f"SOCP : warehouse at {x.value.round(4)}, cost = {prob.value:.4f}")
# Cross-check with Weiszfeld's fixed-point iteration (a classic special-purpose method)
z = demand @ stores / demand.sum() # start at the center of gravity
print(f"centroid : {z.round(4)}, cost = {demand @ np.linalg.norm(stores - z, axis=1):.4f}")
for _ in range(500):
w = demand / np.linalg.norm(stores - z, axis=1) # inverse-distance weights
z = w @ stores / w.sum()
print(f"Weiszfeld: warehouse at {z.round(4)}, cost = {demand @ np.linalg.norm(stores - z, axis=1):.4f}")
# Robust capacity constraint: a'x <= b for every a in an ellipsoid {a0 + P u : ||u|| <= 1}
# is equivalent to the single SOC constraint a0'x + ||P'x|| <= b
a0, P, budget = np.array([2.0, 3.0]), np.diag([0.4, 0.9]), 30.0
y = cp.Variable(2)
nominal = cp.Problem(cp.Maximize(cp.sum(y)), [a0 @ y <= budget, y >= 0, y <= 10])
robust = cp.Problem(cp.Maximize(cp.sum(y)), [a0 @ y + cp.norm(P.T @ y, 2) <= budget, y >= 0, y <= 10])
nominal.solve(); yn = y.value.copy()
robust.solve(); yr = y.value.copy()
print(f"nominal plan {yn.round(3)}, robust plan {yr.round(3)}")
print(f"worst-case use: nominal {a0 @ yn + np.linalg.norm(P.T @ yn):.3f}, robust {a0 @ yr + np.linalg.norm(P.T @ yr):.3f} (budget {budget})")
# SOCP : warehouse at [4.6284 3.7295], cost = 507.7081
# centroid : [4.7647 4.1176], cost = 508.3444
# Weiszfeld: warehouse at [4.6283 3.7295], cost = 507.7081
# nominal plan [10. 3.333], robust plan [10. 1.885]
# worst-case use: nominal 35.000, robust 30.000 (budget 30.0)
The optimal site is not the center of gravity (the weighted average, which minimizes squared distances instead). Weiszfeld's iteration agrees to four decimals, but it is special-purpose and fails if an iterate lands exactly on a store, whereas the SOCP formulation extends directly to multiple facilities, capacity limits or forbidden zones. The nominal plan uses the full budget of 30 for the expected coefficients, but in the worst case it needs 35. The robust plan gives up some output to remain feasible for every coefficient vector in the ellipsoid.
Key takeaways
- The second-order cone is ${(t, u) : \lVert u \rVert_2 \le t}$. SOCP minimizes a linear function subject to affine images lying in such cones.
- SOCP contains LP, convex QP and convex QCQP.
- Euclidean distances, ellipsoidal robust constraints and Gaussian chance constraints are SOC constraints.
- Interior point solvers handle SOCPs nearly as efficiently as LPs, so recognizing SOC structure is valuable.
Practice
- (Beginner) Write the Weber problem for three stores at $(0, 0)$, $(4, 0)$ and $(0, 3)$ with equal weights as an SOCP, with the $t_i$ variables and cone constraints spelled out.
- (Intermediate) Replace the Euclidean distance in the example by the Manhattan distance $\lvert x_1 - a_{i1} \rvert + \lvert x_2 - a_{i2} \rvert$. Show that the problem becomes an LP, solve it, and compare the locations. Why does the Manhattan solution have coordinates that are weighted medians?
- (Interview) "A logistics manager asks for a delivery plan that respects a fuel budget '95% of the time', with normally distributed consumption. How would you model it?" Hint: the chance constraint becomes an SOC constraint with $\Phi^{-1}(0.95) \approx 1.645$. Discuss the normality assumption and the alternatives (scenarios, robust sets).
Semidefinite Programming
Why it matters. Semidefinite programming optimizes over matrices that must be positive semidefinite. It produces the strongest tractable relaxations of many combinatorial problems (the Goemans-Williamson max-cut algorithm is the classic case), and it is central to robust control, sensor network localization and polynomial optimization. It is the most general of the standard convex classes.
Theory. A symmetric matrix $X$ is positive semidefinite, written $X \succeq 0$, if $v^\top X v \ge 0$ for all $v$, or equivalently if all its eigenvalues are nonnegative. The set of PSD matrices is a convex cone. An SDP in standard form is
$$ \min\ \langle C, X \rangle \quad \text{s.t.} \quad \langle A_i, X \rangle = b_i \ \ (i = 1, \dots, m), \qquad X \succeq 0, $$
where $\langle C, X \rangle = \sum_{ij} C_{ij} X_{ij}$. It is an LP over matrices, with "nonnegative" replaced by "positive semidefinite". If all matrices are diagonal, it reduces to an LP. Its dual is
$$ \max\ b^\top y \quad \text{s.t.} \quad C - \sum_i y_i A_i \succeq 0, $$
with weak duality always, and strong duality under a Slater condition (a strictly feasible, positive definite point).
SDP relaxations of combinatorial problems. The standard recipe is lifting:
- Write binary decisions as $x_i \in {-1, +1}$.
- Replace the products $x_i x_j$ by entries of a matrix $X = x x^\top$. That matrix is PSD, has rank 1 and has $X_{ii} = 1$.
- Drop the non-convex rank-1 condition, keeping $X \succeq 0$ and $\operatorname{diag}(X) = 1$. The result is an SDP whose value bounds the combinatorial optimum.
For max-cut (split the nodes of a graph into two sides to cut as many edges as possible), the number of cut edges is $\sum_{(i,j) \in E} \frac{1 - x_i x_j}{2}$, and the relaxation maximizes $\sum_{(i,j)} \frac{1 - X_{ij}}{2}$. Since $X \succeq 0$, $X = V^\top V$: each node gets a unit vector $v_i$ and $X_{ij} = v_i^\top v_j$. Goemans-Williamson rounding cuts the vectors with a random hyperplane through the origin, putting the nodes on either side of it on opposite sides of the cut. The expected cut is at least $0.878$ times the SDP bound, which is the best ratio achievable in polynomial time under the Unique Games Conjecture (Section 4.4).
Example: max-cut on the 5-cycle and the Petersen graph.
import numpy as np
import cvxpy as cp
from itertools import product
def max_cut_sdp(n, edges, rounds=1000, seed=0):
X = cp.Variable((n, n), symmetric=True) # X_ij = v_i . v_j
cut = sum((1 - X[i, j]) / 2 for i, j in edges)
prob = cp.Problem(cp.Maximize(cut), [X >> 0, cp.diag(X) == 1]) # X PSD, unit vectors
prob.solve()
# Recover vectors v_i from X = V'V (eigendecomposition), then random-hyperplane rounding
eigval, eigvec = np.linalg.eigh(X.value)
V = (eigvec * np.sqrt(np.maximum(eigval, 0))).T # column i is v_i
rng = np.random.default_rng(seed)
best = 0
for _ in range(rounds):
side = np.sign(rng.standard_normal(n) @ V) # which side of a random hyperplane
best = max(best, sum(side[i] != side[j] for i, j in edges))
return prob.value, best
def max_cut_brute(n, edges):
return max(sum(s[i] != s[j] for i, j in edges) for s in product([0, 1], repeat=n))
graphs = {"5-cycle": (5, [(i, (i + 1) % 5) for i in range(5)]),
"Petersen": (10, [(i, (i + 1) % 5) for i in range(5)] + [(i, i + 5) for i in range(5)]
+ [(5 + i, 5 + (i + 2) % 5) for i in range(5)])}
for name, (n, E) in graphs.items():
bound, rounded = max_cut_sdp(n, E)
print(f"{name:9s}: SDP bound = {bound:.4f}, rounded cut = {rounded}, true max cut = {max_cut_brute(n, E)}")
# 5-cycle : SDP bound = 4.5225, rounded cut = 4, true max cut = 4 <- bound (25 + 5 sqrt 5)/8
# Petersen : SDP bound = 12.5000, rounded cut = 12, true max cut = 12
For the 5-cycle, the SDP bound is $\tfrac{25 + 5\sqrt{5}}{8} \approx 4.52$, while the best cut is 4 (an odd cycle cannot be fully cut). The rounding recovers an optimal cut in both graphs. The bound itself is valuable: it proves that no cut of the Petersen graph exceeds 12, the kind of certificate a branch-and-bound search uses for pruning (Section 4.2).
Practical limits. Interior point methods for SDP form and factor dense matrices of size $m \times m$ over $n \times n$ variables, so they handle $n$ up to a few thousand. Larger problems use first-order solvers (SCS), low-rank factorizations $X = VV^\top$ (Burer-Monteiro) and exploitation of sparsity.
Key takeaways
- An SDP optimizes a linear function of a symmetric matrix subject to linear constraints and $X \succeq 0$. It contains LP and SOCP as special cases.
- SDP duality mirrors LP duality, with strong duality under strict feasibility.
- Lifting $x x^\top$ to $X$ and dropping the rank constraint gives strong relaxations of binary quadratic problems.
- Goemans-Williamson rounding of the max-cut SDP guarantees 0.878 of the optimum. SDP scale is limited, so large problems need first-order or low-rank methods.
Practice
- (Beginner) Check which of these matrices are PSD using eigenvalues: $\begin{pmatrix} 2 & 1 \ 1 & 2 \end{pmatrix}$, $\begin{pmatrix} 1 & 2 \ 2 & 1 \end{pmatrix}$, $\begin{pmatrix} 1 & 1 \ 1 & 1 \end{pmatrix}$. For which value of $a$ is $\begin{pmatrix} 1 & a \ a & 1 \end{pmatrix}$ PSD?
- (Intermediate) Run the max-cut code on the complete graph $K_5$ and on a random graph with 12 nodes. Report the SDP bound, the rounded cut and the true maximum, and compute the ratio rounded/bound.
- (Interview) "Your MILP for a binary quadratic problem has a weak LP relaxation. How could an SDP help, and what would worry you?" Hint: the SDP bound can be much tighter, but solve time and problem size are the concern. Think of using it at the root node only or for a few nodes.
Subgradient Methods
Why it matters. Many convex functions in OR have kinks: absolute values, maxima of linear functions, piecewise-linear costs, and every Lagrangian dual function of an integer program (Section 4.2). Gradients do not exist at the kinks, and the optimum is usually at a kink. Subgradient methods are the simplest algorithms that work anyway, and they scale to huge problems because each iteration is very cheap.
Theory. For a convex function $f$, a vector $s$ is a subgradient at $x$ if the linear function through $(x, f(x))$ with slope $s$ never lies above $f$:
$$ f(y) \ge f(x) + s^\top (y - x) \quad \text{for all } y . $$
The set of all subgradients is the subdifferential $\partial f(x)$. Where $f$ is differentiable it contains only $\nabla f(x)$. At a kink it contains many vectors: for $f(x) = \lvert x \rvert$, $\partial f(0) = [-1, 1]$.
f(x) = |x| and two of its subgradients at the kink x = 0
\ / . <- slope 0.5: passes through (0, 0)
\ / . and stays below the graph
\ / .
\ / .
\ / .
\ / .
----------------*---------------------------- <- slope 0: also a supporting line
.
.
. any slope in [-1, 1] works: subdifferential of |x| at 0 = [-1, 1]
The optimality condition becomes $0 \in \partial f(x^*)$: a horizontal supporting line exists at the minimum. The subgradient method iterates
$$ x_{k+1} = x_k - \alpha_k s_k, \qquad s_k \in \partial f(x_k), $$
with a projection onto the feasible set if there are constraints. Its key difference from gradient descent: it is not a descent method. A step along $-s_k$ can increase $f$, so the method tracks the best point found so far, and step sizes are fixed in advance instead of chosen by a line search.
| Step rule | Behavior |
|---|---|
| Constant $\alpha$ | converges to within about $\alpha G^2 / 2$ of $f^*$ ($G$ bounds $\lVert s \rVert$), then oscillates |
| Diminishing: $\alpha_k \to 0$, $\sum_k \alpha_k = \infty$ (e.g. $1/k$, $1/\sqrt{k}$) | converges to $f^*$ |
| Polyak: $\alpha_k = \frac{f(x_k) - f^*}{\lVert s_k \rVert^2}$ | best theoretical rule. Needs $f^*$, estimated by a feasible solution value in Lagrangian relaxation (Section 4.2) |
Convergence is slow: reaching accuracy $\varepsilon$ needs $O(1/\varepsilon^2)$ iterations in the worst case. Faster alternatives for nonsmooth convex problems include bundle methods (which keep a model built from many past subgradients), cutting-plane (Kelley) methods, smoothing techniques and proximal methods for objectives with known structure.
Example: the depot location of Section 1.3. Minimizing the total distance $f(x) = \lvert x - 1 \rvert + \lvert x - 4 \rvert + \lvert x - 9 \rvert$ has its optimum at the median, $x^* = 4$, with $f^* = 8$. At $x = 4$ the subdifferential is $1 + [-1, 1] - 1 = [-1, 1]$, which contains 0. The code compares three step rules from $x_0 = 0$:
import numpy as np
# Depot location from Section 1.3: minimize total distance to customers at 1, 4 and 9
a = np.array([1.0, 4.0, 9.0])
f = lambda x: np.abs(x - a).sum()
subgrad = lambda x: np.sign(x - a).sum() # sign(0) = 0 is a valid choice at a kink
def subgradient_method(x, step_rule, iters=200):
best_x, best_f, history = x, f(x), []
for k in range(1, iters + 1):
x = x - step_rule(k) * subgrad(x)
if f(x) < best_f: # not a descent method: track the best
best_x, best_f = x, f(x)
history.append(x)
return best_x, best_f, history
for name, rule in [("constant 0.7 ", lambda k: 0.7),
("diminishing 2/k ", lambda k: 2.0 / k),
("diminishing 2/sqrt(k)", lambda k: 2.0 / np.sqrt(k))]:
bx, bf, hist = subgradient_method(0.0, rule)
print(f"{name:22s}: best x = {bx:.4f}, best f = {bf:.4f}, last iterates {np.round(hist[-4:], 3)}")
# constant 0.7 : best x = 4.2000, best f = 8.2000, last iterates [3.5 4.2 3.5 4.2]
# diminishing 2/k : best x = 4.0000, best f = 8.0000, last iterates [4.01 4. 3.99 4. ]
# diminishing 2/sqrt(k) : best x = 4.0324, best f = 8.0324, last iterates [3.891 4.033 3.891 4.032]
The constant step ends up bouncing between $3.5$ and $4.2$ and never settles. The diminishing steps converge, slowly and without monotone decrease, which is why the best iterate is recorded.
Key takeaways
- A subgradient defines a supporting line or plane below a convex function. The subdifferential generalizes the gradient to kinks.
- Optimality for nonsmooth convex functions is $0 \in \partial f(x^*)$.
- The subgradient method is not a descent method: track the best iterate and use diminishing or Polyak steps.
- It is slow ($O(1/\varepsilon^2)$) but extremely cheap per iteration. It is the standard tool for Lagrangian duals, and bundle methods improve on it.
Practice
- (Beginner) Compute $\partial f(x)$ for $f(x) = \max(2x + 1,\ -x + 4)$ at $x = 0$, $x = 1$ and $x = 2$, and find the minimizer using $0 \in \partial f(x)$.
- (Intermediate) Extend the example to two dimensions: place a depot in the plane to minimize the sum of Manhattan distances to five customers, using the subgradient method with step $1/k$. Compare the result with the coordinate-wise weighted median.
- (Interview) "Your subgradient method for a Lagrangian dual shows the bound oscillating and improving very slowly after 1,000 iterations. What would you try?" Hint: Polyak steps with a better upper bound, reducing $\theta$ when progress stalls, bundle methods, or switching to an LP-based approach such as column generation.
<a id="6-dynamic-programming"></a>