5 Constraints and Duality

Most real problems come with conditions that must hold exactly: a budget that cannot be exceeded, weights that must sum to one, a physical quantity that cannot be negative. Chapter 4 handled the simplest such condition, a box, by clipping. This chapter handles general constraints. The central fact is that at a constrained optimum the gradient of the objective usually does not vanish: the boundary is pushing back, and the force it pushes with, the Lagrange multiplier, turns out to be one of the most useful numbers in optimization. It gives the optimality conditions (KKT), a second, “dual” problem whose answer bounds and certifies the first, a price for every constraint, and two very different algorithms: the barrier method, which keeps the iterates strictly inside the feasible region, and ADMM, which splits a problem into pieces that are each easy. The chapter ends with problems where the unknown is a whole function.

5.1 The problem

Definition — Constrained optimization problem

\[ \min_{x \in \mathbb{R}^n} f(x) \quad \text{subject to} \quad g_i(x) \le 0, \;\; i = 1, \dots, m, \qquad Cx = d . \] The set of points satisfying all constraints is the feasible set. An inequality is active at $x$ if $g_i(x) = 0$ and inactive if $g_i(x) < 0$.

The problem is convex when $f$ and every $g_i$ are convex and the equality constraints are linear; the feasible set is then convex, and Chapter 1’s guarantee applies: every local minimum is global.

5.2 Equality constraints and Lagrange multipliers

Start with a single equality constraint, a curve $\{x : c(x) = 0\}$ in the plane. Walk along the curve. At a constrained minimum, $f$ cannot decrease in either direction along it, so the directional derivative of $f$ along the curve’s tangent is zero: $\nabla f$ is perpendicular to the curve. But $\nabla c$ is also perpendicular to the curve (it is perpendicular to its own level set, by Chapter 1). Two vectors perpendicular to the same line are parallel, so \[ \nabla f(x^\star) + \nu\, \nabla c(x^\star) = 0 \quad \text{for some number } \nu . \] The Lagrangian $\mathcal{L}(x, \nu) = f(x) + \nu\, c(x)$ packages this: its gradient in $x$ is the condition above, its gradient in $\nu$ is the constraint itself. A constrained problem in $x$ becomes the search for a stationary point of $\mathcal{L}$ in $(x, \nu)$.

The multiplier has a meaning beyond bookkeeping. If the constraint is relaxed to $c(x) = \epsilon$, the optimal value changes at rate $-\nu$: \[ \frac{d f^\star}{d\epsilon} = -\nu . \] The multiplier is the price of the constraint, the marginal value of loosening it.

Equality-constrained quadratic programs

For a quadratic objective $\tfrac12 x^\top Q x + q^\top x$ with linear constraints $Cx = d$, the stationarity conditions are linear, and the whole problem reduces to one linear system, the KKT system: \[ \begin{bmatrix} Q & C^\top \\ C & 0 \end{bmatrix} \begin{bmatrix} x \\ \nu \end{bmatrix} = \begin{bmatrix} -q \\ d \end{bmatrix} . \] Least squares with linear constraints is a special case. Forcing the slope and intercept of a fitted line to sum to one:

Show the code

import numpy as np

from optimlab.linalg import equality_constrained_least_squares, least_squares
from optimlab.viz import markdown_table

rng = np.random.default_rng(0)
t = np.linspace(0, 10, 40)
y = 2.0 * t + 1.0 + rng.normal(scale=2.0, size=t.size)
A = np.column_stack([t, np.ones_like(t)])
C, d = np.array([[1.0, 1.0]]), np.array([1.0])  # slope + intercept = 1

free = least_squares(A, y).x
constrained = equality_constrained_least_squares(A, y, C, d)
markdown_table(
    ["fit", "slope", "intercept", "slope + intercept"],
    [("unconstrained", f"{free[0]:.3f}", f"{free[1]:.3f}", f"{free.sum():.3f}"),
     ("constrained", f"{constrained[0]:.3f}", f"{constrained[1]:.3f}",
      f"{constrained.sum():.3f}")],
)

fit	slope	intercept	slope + intercept
unconstrained	2.062	0.570	2.632
constrained	2.325	-1.325	1.000

The constraint holds exactly, and the fit gives up as little squared error as it can to satisfy it. Chapter 9 uses the same KKT system to compute an efficient investment portfolio.

Implementation: linalg/qp.py · Checked by: test_linalg.py

5.3 Inequality constraints: the KKT conditions

An inequality $g_i(x) \le 0$ behaves like an equality when it is active and like nothing at all when it is not. It can also push in only one direction: it can stop $x$ from leaving the feasible set, never pull it outward. These two facts become the last two of the four conditions below.

Definition — Karush–Kuhn–Tucker (KKT) conditions

With Lagrangian $\mathcal{L}(x, \lambda, \nu) = f(x) + \sum_i \lambda_i g_i(x) + \nu^\top (Cx - d)$, a point $x^\star$ with multipliers $\lambda^\star, \nu^\star$ satisfies the KKT conditions if

Stationarity: $\nabla f(x^\star) + \sum_i \lambda_i^\star \nabla g_i(x^\star) + C^\top \nu^\star = 0$;
Primal feasibility: $g_i(x^\star) \le 0$ and $Cx^\star = d$;
Dual feasibility: $\lambda_i^\star \ge 0$;
Complementary slackness: $\lambda_i^\star\, g_i(x^\star) = 0$ for every $i$.

Theorem — KKT conditions are necessary, and for convex problems sufficient

Under a mild regularity assumption on the constraints (a constraint qualification; for convex problems, Slater’s condition, the existence of a strictly feasible point, suffices), every local minimum satisfies the KKT conditions for some multipliers. For a convex problem, any point satisfying them is a global minimum (Boyd and Vandenberghe 2004, sec. 5.5).

Read together, the conditions say: at the optimum, the direction that would improve $f$ fastest, $-\nabla f(x^\star)$, is exactly balanced by the outward normals of the active constraints, each scaled by a nonnegative multiplier. If it were not balanced, some feasible direction would still improve $f$. Complementary slackness says an inactive constraint has a zero multiplier and plays no part in the balance; a constraint that is not binding has no price. These conditions generalize $\nabla f(x^\star) = 0$, which is what they reduce to when no constraint is active.

5.4 Duality

The multipliers do more than describe the optimum. Each choice of them defines a lower bound on the optimal value.

Definition — Lagrange dual problem

The dual function is $q(\lambda, \nu) = \inf_x \mathcal{L}(x, \lambda, \nu)$, and the dual problem is \[ \max_{\lambda \ge 0,\, \nu} \; q(\lambda, \nu) . \]

Weak duality takes one line: for any feasible $x$ and any $\lambda \ge 0$, the added terms $\lambda_i g_i(x)$ are $\le 0$ and $\nu^\top(Cx - d) = 0$, so $q(\lambda, \nu) \le \mathcal{L}(x, \lambda, \nu) \le f(x)$. Every dual-feasible point bounds the primal optimum from below, whether or not the problem is convex. The difference between a primal objective value and a dual one is the duality gap, and a small gap is a certificate that a candidate is nearly optimal. Strong duality, a zero gap at the optimum, holds for convex problems that satisfy Slater’s condition. The dual function is always concave, so the dual problem is convex even when the primal is not.

The dual of a linear program, and shadow prices

For a linear program the dual is again a linear program. Take the workshop problem of Chapter 3, $\max\, 3x_1 + 5x_2$ subject to three machine-hour limits $Ax \le b$ with $b = (4, 12, 18)$. Its dual is \[ \min_{y \ge 0} \; 4y_1 + 12y_2 + 18y_3 \quad \text{s.t.} \quad y_1 + 3y_3 \ge 3, \quad 2y_2 + 2y_3 \ge 5 . \] Read $y_i$ as the price of one hour on machine $i$. The constraints say that the hours needed to make one unit of each product are worth at least the profit the product earns, and the objective is the total value of all available hours. Strong duality says that at the optimum, the total value of the hours equals the maximum profit.

Show the code

from optimlab.optimizers.linear_programming import LinearProgram, dual, simplex

def workshop(hours):
    return LinearProgram(c=[-3.0, -5.0], A_ub=[[1, 0], [0, 2], [3, 2]], b_ub=hours)

hours = [4.0, 12.0, 18.0]
primal = simplex(workshop(hours))
prices = simplex(dual(workshop(hours))).x
rows = []
for i in range(3):
    more = list(hours)
    more[i] += 1
    gain = -simplex(workshop(more)).objective - (-primal.objective)
    rows.append((f"machine {i + 1}", f"{prices[i]:.2f}", f"{gain:.2f}"))
markdown_table(["constraint", "dual variable (shadow price)", "profit gained from +1 hour"], rows)

Table 5.1: The workshop problem’s dual solution next to a direct experiment: add one hour to each machine and re-solve. The dual variables are exactly the marginal profit of each machine’s hours. Machine 1 is not fully used at the optimum, so its price is zero (complementary slackness).

constraint	dual variable (shadow price)	profit gained from +1 hour
machine 1	0.00	0.00
machine 2	1.50	1.50
machine 3	1.00	1.00

The maximum profit is 36, and so is the value of the hours at their shadow prices, 36: zero duality gap.

Implementation: optimizers/linear_programming.py · Checked by: test_linear_programming.py

5.5 The barrier method

The KKT conditions say what an optimum looks like. The barrier method is one way to find it, for problems with inequality constraints: replace each hard wall with a smooth penalty that becomes infinite at the wall, solve the resulting unconstrained problem with Newton’s method, and gradually make the penalty weaker relative to the objective.

Definition — Logarithmic barrier and central path

For $t > 0$, let \[ x^\star(t) = \arg\min_x \; t\, f(x) - \sum_{i=1}^m \log\big(-g_i(x)\big). \] The curve $t \mapsto x^\star(t)$ is the central path. As $t \to \infty$ it converges to a solution of the constrained problem.

The logarithm is defined only where every $g_i(x) < 0$, so each barrier problem lives strictly inside the feasible set and its minimizer stays away from the walls. Increasing $t$ shifts the balance toward $f$, and the minimizer slides toward the boundary.

Why the central path converges, and how fast

Setting the gradient of the barrier objective to zero gives \[ t\, \nabla f(x) + \sum_i \frac{1}{-g_i(x)}\, \nabla g_i(x) = 0 \quad\Longleftrightarrow\quad \nabla f(x) + \sum_i \lambda_i \nabla g_i(x) = 0 \;\;\text{with}\;\; \lambda_i = \frac{1}{-t\, g_i(x)} > 0 . \] That is KKT stationarity, with positive multipliers. Complementary slackness does not hold exactly: $\lambda_i\,(-g_i(x)) = 1/t$ instead of 0. The duality gap of this primal–dual pair is therefore exactly $\sum_i \lambda_i\,(-g_i) = m/t$. Each point on the central path certifies its own accuracy: it is within $m/t$ of optimal.

Algorithm — Barrier method

Given a strictly feasible $x$, an initial $t > 0$, a growth factor $\gamma > 1$ (typically 10 to 20) and a tolerance $\epsilon$:

Centering: starting from the current $x$, minimize $t f(x) - \sum_i \log(-g_i(x))$ with Newton’s method. The line search rejects any trial point that leaves the feasible set before testing for sufficient decrease.
Stop if $m/t < \epsilon$.
Increase $t \leftarrow \gamma t$ and go to step 1, warm-starting from the current point.

The warm start is what makes this efficient: the previous center is close to the next one, so each centering step needs only a few Newton iterations. The number of outer iterations is about $\log\big(m/(t_0\epsilon)\big)/\log \gamma$ (Boyd and Vandenberghe 2004, sec. 11.3). Interior-point methods for linear and convex programming refine this idea, and they solve such problems in polynomial time.

A small example: minimize $(x_1 - 3)^2 + (x_2 - 3)^2$, a bowl centered at $(3, 3)$, subject to $x_1 \le 1$ and $x_2 \le 1$. The constrained optimum is the corner $(1, 1)$, where both constraints are active.

Show the code

from optimlab.optimizers.barrier_method import ConstrainedProblem, barrier_method

problem = ConstrainedProblem(
    f=lambda x: (x[0] - 3.0) ** 2 + (x[1] - 3.0) ** 2,
    x0=np.array([0.5, 0.5]),
    inequality_constraints=[lambda x: x[0] - 1.0, lambda x: x[1] - 1.0],
    name="corner",
)
barrier = barrier_method(problem)

Show the code

from optimlab.viz import central_path_figure

fig = central_path_figure(problem, barrier)
fig.update_layout(title="The central path")
fig.show()

Figure 5.1: The central path for the corner problem. Each marker is the minimizer of one barrier problem; as t grows they move from the interior toward the corner (1, 1), reaching it only in the limit.

Show the code

from optimlab.viz import duality_gap_figure

duality_gap_figure(barrier.grad_norm_trajectory, tol=1e-8).show()

Figure 5.2: The duality gap m/t after each outer iteration, on a log scale. Multiplying t by a fixed factor each round divides the gap by the same factor: a straight line, a geometric rate.

The method reaches (1.000000, 1.000000) in 10 outer iterations. At that point the multipliers can be read off the formula above, $\lambda_i = 1/(-t\, g_i(x))$, and stationarity can be checked directly by drawing the vectors:

Show the code

from optimlab.viz import kkt_geometry_figure

kkt_geometry_figure(problem, barrier).show()

Figure 5.3: KKT stationarity at the corner optimum: -∇f(x*) and the two active constraints’ gradients, each scaled by its multiplier. The scaled constraint gradients add up (dashed) to exactly -∇f. The arrows are drawn with a small sideways offset only so that they do not overlap.

Both multipliers are positive (both constraints are holding the optimum in place) and, by the problem’s symmetry, equal.

Implementation: optimizers/barrier_method.py · Checked by: test_barrier_method.py

5.6 ADMM: splitting a problem into easy pieces

The barrier method needs derivatives of everything. Many useful problems are sums of pieces that are each easy to handle on their own but awkward together: a smooth data-fit term plus a nonsmooth penalty, or terms that live on different machines. The alternating direction method of multipliers (Boyd et al. 2011) handles such sums by giving each piece its own copy of the variable and using a multiplier to make the copies agree.

From the augmented Lagrangian to ADMM

Write $\min_x f(x) + h(x)$ as $\min_{x, z} f(x) + h(z)$ subject to $x = z$. The augmented Lagrangian adds a quadratic penalty to the ordinary one: \[ \mathcal{L}_\rho(x, z, \nu) = f(x) + h(z) + \nu^\top (x - z) + \frac{\rho}{2}\lVert x - z \rVert^2 . \] The classical method of multipliers minimizes $\mathcal{L}_\rho$ over $(x, z)$ jointly, then moves the multiplier uphill on the dual, $\nu \leftarrow \nu + \rho(x - z)$. ADMM makes one change: instead of minimizing over $x$ and $z$ jointly, which is as hard as the original problem, it minimizes over $x$ and then over $z$, once each. With the scaled multiplier $u = \nu/\rho$, both minimizations become proximal steps (Chapter 4):

Definition — ADMM (scaled form)

\[ \begin{aligned} x_{k+1} &= \operatorname{prox}_{f/\rho}\big(z_k - u_k\big), \\ z_{k+1} &= \operatorname{prox}_{h/\rho}\big(x_{k+1} + u_k\big), \\ u_{k+1} &= u_k + x_{k+1} - z_{k+1} . \end{aligned} \]

Each step touches only one piece of the objective, and neither needs to be differentiable. For convex $f$ and $h$, ADMM converges for every $\rho > 0$. The penalty $\rho$ is a speed setting, not a stability condition.

Stopping: two residuals, not one

Two things must go to zero. The primal residual $\lVert x_k - z_k \rVert$ measures whether the copies agree. The dual residual $\rho \lVert z_k - z_{k-1} \rVert$ measures whether $z$, and with it the multiplier, has stopped moving. It is tempting to stop on the primal residual alone. That is unsafe: $x$ and $z$ can agree by coincidence in the middle of a run, long before either is correct. On a lasso problem with $\rho = 50$ it happens within a few iterations:

Show the code

import plotly.graph_objects as go

from optimlab.optimizers.admm import ADMMProblem, admm
from optimlab.optimizers.proximal_gradient import CompositeProblem, proximal_gradient, soft_threshold
from optimlab.viz.theme import contrasting_categorical, layout_template

rng = np.random.default_rng(0)
m, n = 50, 20
A = rng.standard_normal((m, n))
x_true = np.zeros(n)
x_true[[2, 5, 9]] = [3.0, -2.0, 1.5]
b = A @ x_true + 0.01 * rng.standard_normal(m)
alpha = 1.0
AtA, Atb = A.T @ A, A.T @ b

def lasso_admm(rho):
    return ADMMProblem(
        prox_f=lambda v, t: np.linalg.solve(AtA + np.eye(n) / t, Atb + v / t),
        prox_g=lambda v, t: soft_threshold(v, alpha * t),
        x0=np.zeros(n),
        f_obj=lambda x: 0.5 * np.sum((A @ x - b) ** 2),
        g_obj=lambda x: alpha * np.sum(np.abs(x)),
    )

pg = proximal_gradient(
    CompositeProblem(grad_smooth=lambda x: A.T @ (A @ x - b),
                     prox_nonsmooth=lambda v, t: soft_threshold(v, alpha * t), x0=np.zeros(n)),
    lr=1.0 / np.linalg.eigvalsh(AtA).max(), max_iter=5000, tol=1e-10,
)
f_star = 0.5 * np.sum((A @ pg.x - b) ** 2) + alpha * np.sum(np.abs(pg.x))
run = admm(lasso_admm(50.0), rho=50.0, max_iter=5000, tol=1e-6, record_history=True)
primal = np.array(run.history["primal_residual"])
dual_res = np.array(run.history["dual_residual"])
first_dip = int(np.argmax(primal < 1e-6)) + 1
gap_at_dip = (run.f_trajectory[first_dip] - f_star) / f_star

Show the code

colors = contrasting_categorical()
iters = np.arange(1, run.n_iter + 1)
fig = go.Figure([
    go.Scatter(x=iters, y=np.maximum(primal, 1e-17), name="primal residual ‖x − z‖",
               line={"color": colors[0], "width": 2}),
    go.Scatter(x=iters, y=np.maximum(dual_res, 1e-17), name="dual residual ρ‖zₖ − zₖ₋₁‖",
               line={"color": colors[1], "width": 2}),
    go.Scatter(x=iters, y=np.maximum(np.array(run.f_trajectory[1:]) - f_star, 1e-17),
               name="objective − optimum", line={"color": colors[2], "width": 2, "dash": "dot"}),
])
fig.update_layout(**layout_template(xaxis_title="iteration", yaxis_title="value",
                                    yaxis={"type": "log", "exponentformat": "power"},
                                    legend={"y": 1, "x": 1.02, "xanchor": "left"},
                                    margin={"r": 230}))
fig.show()

Figure 5.4: ADMM on a lasso problem with ρ = 50 (log scale). The primal residual collapses to numerical zero early, a coincidence of x and z agreeing, while the dual residual and the objective’s distance from the optimum are still large. The primal residual then rises again. Only when both residuals are small is the run actually finished.

At iteration 5 the primal residual is below $10^{-6}$ while the objective is still 50% above its optimal value. A primal-only rule would have stopped there. The standard rule requires both residuals to be small (Boyd et al. 2011, sec. 3.3.1), and the dual residual is far from small at that point.

ADMM and proximal gradient agree

The lasso was solved in Chapter 4 by proximal gradient, which takes an explicit gradient step. ADMM never takes one, yet the two land on the same answer:

Show the code

rows = []
for rho in [0.5, 1.0, 5.0, 20.0, 50.0]:
    r = admm(lasso_admm(rho), rho=rho, max_iter=5000, tol=1e-6)
    rows.append((rho, r.n_iter, f"{r.f:.4f}", f"{np.max(np.abs(r.x - pg.x)):.1e}"))
markdown_table(["ρ", "iterations", "objective", "max difference from proximal gradient"], rows)

Table 5.2: ADMM on the same lasso problem for several values of ρ. All converge to the proximal-gradient solution; ρ only changes how many iterations it takes.

ρ	iterations	objective	max difference from proximal gradient
0.5	1877	6.4663	4.5e-07
1.0	939	6.4663	4.5e-07
5.0	190	6.4663	4.4e-07
20.0	62	6.4663	3.3e-07
50.0	43	6.4663	3.7e-08

Implementation: optimizers/admm.py · Checked by: test_admm.py

5.7 When the unknown is a function: calculus of variations

So far the unknown has been a vector. In the calculus of variations it is a whole function $y(x)$, and the objective is an integral, a functional: \[ J[y] = \int_a^b F\big(x, y(x), y'(x)\big)\, dx, \qquad y(a), y(b) \text{ fixed}. \] The first-order condition is derived the same way as $\nabla f = 0$: perturb the candidate and require the first-order change to vanish. Let $y + \epsilon\eta$ be a perturbed curve, with $\eta(a) = \eta(b) = 0$ so the endpoints stay fixed. Then \[ \frac{d}{d\epsilon} J[y + \epsilon \eta]\Big|_{\epsilon = 0} = \int_a^b \big( F_y\, \eta + F_{y'}\, \eta' \big)\, dx = \int_a^b \Big( F_y - \frac{d}{dx} F_{y'} \Big)\, \eta\, dx , \] after integrating the second term by parts (the boundary term vanishes because $\eta$ does at both ends). For this to be zero for every perturbation $\eta$, the bracket must vanish everywhere:

Definition — Euler–Lagrange equation

\[ \frac{\partial F}{\partial y} - \frac{d}{dx}\frac{\partial F}{\partial y'} = 0 . \]

It plays exactly the role of $\nabla f(x^\star) = 0$, now as a differential equation (Gelfand and Fomin 2000). A classic instance is the soap film stretched between two coaxial rings of radius $y_0$ a distance $\ell$ apart. The film minimizes its area, $J[y] = 2\pi \int_0^\ell y \sqrt{1 + y'^2}\, dx$. Because $F$ does not depend on $x$ explicitly, the Euler–Lagrange equation has a first integral, $F - y' F_{y'} = \text{const}$, which here reads $y / \sqrt{1 + y'^2} = c$, and its solution is a catenary, \[ y(x) = c \cosh\!\Big(\frac{x - \ell/2}{c}\Big), \qquad y_0 = c\cosh\!\Big(\frac{\ell}{2c}\Big). \]

The same problem can be solved without any of this calculus. Represent $y$ by its values on a grid of $N + 1$ points, replace the integral with a sum, and the functional becomes an ordinary function of $N - 1$ unknowns (the interior heights), which BFGS (Chapter 2) can minimize:

Show the code

import jax.numpy as jnp
from scipy.optimize import brentq

from optimlab.core import Problem
from optimlab.optimizers import bfgs

y0, ell, N = 1.0, 1.0, 40
xs = np.linspace(0.0, ell, N + 1)
dx = xs[1] - xs[0]

def area(y_inner):  # the 2π factor does not change the minimizer
    y = jnp.concatenate([jnp.array([y0]), y_inner, jnp.array([y0])])
    return jnp.sum((y[:-1] + y[1:]) / 2 * jnp.sqrt(1.0 + (jnp.diff(y) / dx) ** 2)) * dx

film = bfgs(Problem(f=area, x0=np.full(N - 1, y0)), max_iter=500)
y_opt = np.concatenate([[y0], film.x, [y0]])
c = brentq(lambda c: c * np.cosh(ell / (2 * c)) - y0, 0.5, 3.0)  # the stable (larger) root
y_catenary = c * np.cosh((xs - ell / 2) / c)
max_gap = np.max(np.abs(y_opt - y_catenary))

colors = contrasting_categorical()
fig = go.Figure([
    go.Scatter(x=xs, y=np.full(N + 1, y0), mode="lines", name="flat starting guess",
               line={"color": "#898781", "width": 1.5, "dash": "dot"}),
    go.Scatter(x=xs, y=y_catenary, mode="lines", name="catenary (Euler–Lagrange)",
               line={"color": "#52514e", "width": 4}),
    go.Scatter(x=xs, y=y_opt, mode="markers", name="BFGS on the discretized area",
               marker={"size": 6, "color": colors[0]}),
])
fig.update_layout(**layout_template(title="A soap film between two rings",
                                    xaxis_title="x", yaxis_title="radius y(x)"))
fig.show()

Figure 5.5: Discretize, then optimize: starting from a flat film, BFGS on 39 interior heights bows the curve inward until it lies on the analytic catenary from the Euler–Lagrange equation.

The largest gap between the optimized heights and the catenary is 1.4e-05. Solving the differential equation and minimizing a fine discretization of the functional are, in the limit, the same computation. Direct discretization is how most optimal-control problems are solved in practice, as Chapter 9 shows.

5.8 What these examples cannot show

Scale. Real interior-point solvers exploit sparsity in the Newton systems and handle equality constraints and infeasible starting points; the two-variable example shows the idea, not the engineering.
Non-convex constraints. The KKT conditions remain necessary, but duality gaps can be positive, and a KKT point need not be a global optimum.
ADMM’s speed. ADMM reaches moderate accuracy quickly and high accuracy slowly. Its advantages, decomposition and parallelism, appear only on problems far larger than these.

5.9 Connections

Chapter 1: the KKT conditions generalize $\nabla f(x^\star) = 0$, and convexity is what makes them sufficient.
Chapter 3: duality gives every linear program a twin, and its variables are the constraints’ shadow prices.
Chapter 4: proximal operators are ADMM’s building blocks, and box constraints are the simplest constraints of all.
Chapter 9 applies the KKT system to portfolio selection, the barrier method to fair resource allocation, and discretize-then-optimize to optimal control.

Boyd, Stephen, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein. 2011. “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.” Foundations and Trends in Machine Learning 3 (1). https://web.stanford.edu/~boyd/papers/admm_distr_stats.html.

Boyd, Stephen, and Lieven Vandenberghe. 2004. Convex Optimization. Cambridge University Press. https://web.stanford.edu/~boyd/cvxbook/.

Gelfand, Israel M., and Sergei V. Fomin. 2000. Calculus of Variations. Dover Publications.

# Constraints and Duality {#sec-constraints} Most real problems come with conditions that must hold exactly: a budget that cannot be exceeded, weights that must sum to one, a physical quantity that cannot be negative. @sec-least-squares handled the simplest such condition, a box, by clipping. This chapter handles general constraints. The central fact is that at a constrained optimum the gradient of the objective usually does *not* vanish: the boundary is pushing back, and the force it pushes with, the **Lagrange multiplier**, turns out to be one of the most useful numbers in optimization. It gives the optimality conditions (KKT), a second, "dual" problem whose answer bounds and certifies the first, a price for every constraint, and two very different algorithms: the barrier method, which keeps the iterates strictly inside the feasible region, and ADMM, which splits a problem into pieces that are each easy. The chapter ends with problems where the unknown is a whole function. ## The problem ::: {.callout-note title="Definition — Constrained optimization problem"} $$ \min_{x \in \mathbb{R}^n} f(x) \quad \text{subject to} \quad g_i(x) \le 0, \;\; i = 1, \dots, m, \qquad Cx = d . $$ The set of points satisfying all constraints is the **feasible set**. An inequality is **active** at $x$ if $g_i(x) = 0$ and **inactive** if $g_i(x) < 0$. ::: The problem is **convex** when $f$ and every $g_i$ are convex and the equality constraints are linear; the feasible set is then convex, and @sec-geometry's guarantee applies: every local minimum is global. ## Equality constraints and Lagrange multipliers Start with a single equality constraint, a curve $\{x : c(x) = 0\}$ in the plane. Walk along the curve. At a constrained minimum, $f$ cannot decrease in either direction along it, so the directional derivative of $f$ along the curve's tangent is zero: $\nabla f$ is perpendicular to the curve. But $\nabla c$ is also perpendicular to the curve (it is perpendicular to its own level set, by @sec-geometry). Two vectors perpendicular to the same line are parallel, so $$ \nabla f(x^\star) + \nu\, \nabla c(x^\star) = 0 \quad \text{for some number } \nu . $$ The **Lagrangian** $\mathcal{L}(x, \nu) = f(x) + \nu\, c(x)$ packages this: its gradient in $x$ is the condition above, its gradient in $\nu$ is the constraint itself. A constrained problem in $x$ becomes the search for a stationary point of $\mathcal{L}$ in $(x, \nu)$. The multiplier has a meaning beyond bookkeeping. If the constraint is relaxed to $c(x) = \epsilon$, the optimal value changes at rate $-\nu$: $$ \frac{d f^\star}{d\epsilon} = -\nu . $$ The multiplier is the **price** of the constraint, the marginal value of loosening it. ### Equality-constrained quadratic programs For a quadratic objective $\tfrac12 x^\top Q x + q^\top x$ with linear constraints $Cx = d$, the stationarity conditions are linear, and the whole problem reduces to one linear system, the **KKT system**: $$ \begin{bmatrix} Q & C^\top \\ C & 0 \end{bmatrix} \begin{bmatrix} x \\ \nu \end{bmatrix} = \begin{bmatrix} -q \\ d \end{bmatrix} . $$ Least squares with linear constraints is a special case. Forcing the slope and intercept of a fitted line to sum to one: ```{python} #| label: kkt-qp import numpy as np from optimlab.linalg import equality_constrained_least_squares, least_squares from optimlab.viz import markdown_table rng = np.random.default_rng(0) t = np.linspace(0, 10, 40) y = 2.0 * t + 1.0 + rng.normal(scale=2.0, size=t.size) A = np.column_stack([t, np.ones_like(t)]) C, d = np.array([[1.0, 1.0]]), np.array([1.0]) # slope + intercept = 1 free = least_squares(A, y).x constrained = equality_constrained_least_squares(A, y, C, d) markdown_table( ["fit", "slope", "intercept", "slope + intercept"], [("unconstrained", f"{free[0]:.3f}", f"{free[1]:.3f}", f"{free.sum():.3f}"), ("constrained", f"{constrained[0]:.3f}", f"{constrained[1]:.3f}", f"{constrained.sum():.3f}")], ) ``` The constraint holds exactly, and the fit gives up as little squared error as it can to satisfy it. @sec-applications uses the same KKT system to compute an efficient investment portfolio. [Implementation: [`linalg/qp.py`]({{< var src >}}/linalg/qp.py) · Checked by: [`test_linalg.py`]({{< var tests >}}/test_linalg.py)]{.impl} ## Inequality constraints: the KKT conditions An inequality $g_i(x) \le 0$ behaves like an equality when it is active and like nothing at all when it is not. It can also push in only one direction: it can stop $x$ from leaving the feasible set, never pull it outward. These two facts become the last two of the four conditions below. ::: {.callout-note title="Definition — Karush–Kuhn–Tucker (KKT) conditions"} With Lagrangian $\mathcal{L}(x, \lambda, \nu) = f(x) + \sum_i \lambda_i g_i(x) + \nu^\top (Cx - d)$, a point $x^\star$ with multipliers $\lambda^\star, \nu^\star$ satisfies the KKT conditions if 1. **Stationarity**: $\nabla f(x^\star) + \sum_i \lambda_i^\star \nabla g_i(x^\star) + C^\top \nu^\star = 0$; 2. **Primal feasibility**: $g_i(x^\star) \le 0$ and $Cx^\star = d$; 3. **Dual feasibility**: $\lambda_i^\star \ge 0$; 4. **Complementary slackness**: $\lambda_i^\star\, g_i(x^\star) = 0$ for every $i$. ::: ::: {.callout-tip title="Theorem — KKT conditions are necessary, and for convex problems sufficient"} Under a mild regularity assumption on the constraints (a *constraint qualification*; for convex problems, Slater's condition, the existence of a strictly feasible point, suffices), every local minimum satisfies the KKT conditions for some multipliers. For a convex problem, any point satisfying them is a global minimum [@boyd2004convex, §5.5]. ::: Read together, the conditions say: at the optimum, the direction that would improve $f$ fastest, $-\nabla f(x^\star)$, is exactly balanced by the outward normals of the active constraints, each scaled by a nonnegative multiplier. If it were not balanced, some feasible direction would still improve $f$. Complementary slackness says an inactive constraint has a zero multiplier and plays no part in the balance; a constraint that is not binding has no price. These conditions generalize $\nabla f(x^\star) = 0$, which is what they reduce to when no constraint is active. ## Duality The multipliers do more than describe the optimum. Each choice of them defines a lower bound on the optimal value. ::: {.callout-note title="Definition — Lagrange dual problem"} The **dual function** is $q(\lambda, \nu) = \inf_x \mathcal{L}(x, \lambda, \nu)$, and the **dual problem** is $$ \max_{\lambda \ge 0,\, \nu} \; q(\lambda, \nu) . $$ ::: **Weak duality** takes one line: for any feasible $x$ and any $\lambda \ge 0$, the added terms $\lambda_i g_i(x)$ are $\le 0$ and $\nu^\top(Cx - d) = 0$, so $q(\lambda, \nu) \le \mathcal{L}(x, \lambda, \nu) \le f(x)$. Every dual-feasible point bounds the primal optimum from below, whether or not the problem is convex. The difference between a primal objective value and a dual one is the **duality gap**, and a small gap is a certificate that a candidate is nearly optimal. **Strong duality**, a zero gap at the optimum, holds for convex problems that satisfy Slater's condition. The dual function is always concave, so the dual problem is convex even when the primal is not. ### The dual of a linear program, and shadow prices For a linear program the dual is again a linear program. Take the workshop problem of @sec-linear-programming, $\max\, 3x_1 + 5x_2$ subject to three machine-hour limits $Ax \le b$ with $b = (4, 12, 18)$. Its dual is $$ \min_{y \ge 0} \; 4y_1 + 12y_2 + 18y_3 \quad \text{s.t.} \quad y_1 + 3y_3 \ge 3, \quad 2y_2 + 2y_3 \ge 5 . $$ Read $y_i$ as the price of one hour on machine $i$. The constraints say that the hours needed to make one unit of each product are worth at least the profit the product earns, and the objective is the total value of all available hours. Strong duality says that at the optimum, the total value of the hours equals the maximum profit. ```{python} #| label: tbl-shadow-prices #| tbl-cap: "The workshop problem's dual solution next to a direct experiment: add one hour to each machine and re-solve. The dual variables are exactly the marginal profit of each machine's hours. Machine 1 is not fully used at the optimum, so its price is zero (complementary slackness)." from optimlab.optimizers.linear_programming import LinearProgram, dual, simplex def workshop(hours): return LinearProgram(c=[-3.0, -5.0], A_ub=[[1, 0], [0, 2], [3, 2]], b_ub=hours) hours = [4.0, 12.0, 18.0] primal = simplex(workshop(hours)) prices = simplex(dual(workshop(hours))).x rows = [] for i in range(3): more = list(hours) more[i] += 1 gain = -simplex(workshop(more)).objective - (-primal.objective) rows.append((f"machine {i + 1}", f"{prices[i]:.2f}", f"{gain:.2f}")) markdown_table(["constraint", "dual variable (shadow price)", "profit gained from +1 hour"], rows) ``` The maximum profit is `{python} f"{-primal.objective:.0f}"`, and so is the value of the hours at their shadow prices, `{python} f"{np.dot(prices, hours):.0f}"`: zero duality gap. [Implementation: [`optimizers/linear_programming.py`]({{< var src >}}/optimizers/linear_programming.py) · Checked by: [`test_linear_programming.py`]({{< var tests >}}/test_linear_programming.py)]{.impl} ## The barrier method The KKT conditions say what an optimum looks like. The barrier method is one way to find it, for problems with inequality constraints: replace each hard wall with a smooth penalty that becomes infinite at the wall, solve the resulting unconstrained problem with Newton's method, and gradually make the penalty weaker relative to the objective. ::: {.callout-note title="Definition — Logarithmic barrier and central path"} For $t > 0$, let $$ x^\star(t) = \arg\min_x \; t\, f(x) - \sum_{i=1}^m \log\big(-g_i(x)\big). $$ The curve $t \mapsto x^\star(t)$ is the **central path**. As $t \to \infty$ it converges to a solution of the constrained problem. ::: The logarithm is defined only where every $g_i(x) < 0$, so each barrier problem lives strictly inside the feasible set and its minimizer stays away from the walls. Increasing $t$ shifts the balance toward $f$, and the minimizer slides toward the boundary. ### Why the central path converges, and how fast Setting the gradient of the barrier objective to zero gives $$ t\, \nabla f(x) + \sum_i \frac{1}{-g_i(x)}\, \nabla g_i(x) = 0 \quad\Longleftrightarrow\quad \nabla f(x) + \sum_i \lambda_i \nabla g_i(x) = 0 \;\;\text{with}\;\; \lambda_i = \frac{1}{-t\, g_i(x)} > 0 . $$ That is KKT stationarity, with positive multipliers. Complementary slackness does not hold exactly: $\lambda_i\,(-g_i(x)) = 1/t$ instead of 0. The duality gap of this primal–dual pair is therefore exactly $\sum_i \lambda_i\,(-g_i) = m/t$. Each point on the central path certifies its own accuracy: it is within $m/t$ of optimal. ::: {.callout-note title="Algorithm — Barrier method"} Given a strictly feasible $x$, an initial $t > 0$, a growth factor $\gamma > 1$ (typically 10 to 20) and a tolerance $\epsilon$: 1. **Centering**: starting from the current $x$, minimize $t f(x) - \sum_i \log(-g_i(x))$ with Newton's method. The line search rejects any trial point that leaves the feasible set before testing for sufficient decrease. 2. **Stop** if $m/t < \epsilon$. 3. **Increase** $t \leftarrow \gamma t$ and go to step 1, warm-starting from the current point. ::: The warm start is what makes this efficient: the previous center is close to the next one, so each centering step needs only a few Newton iterations. The number of outer iterations is about $\log\big(m/(t_0\epsilon)\big)/\log \gamma$ [@boyd2004convex, §11.3]. Interior-point methods for linear and convex programming refine this idea, and they solve such problems in polynomial time. A small example: minimize $(x_1 - 3)^2 + (x_2 - 3)^2$, a bowl centered at $(3, 3)$, subject to $x_1 \le 1$ and $x_2 \le 1$. The constrained optimum is the corner $(1, 1)$, where both constraints are active. ```{python} #| label: barrier-setup from optimlab.optimizers.barrier_method import ConstrainedProblem, barrier_method problem = ConstrainedProblem( f=lambda x: (x[0] - 3.0) ** 2 + (x[1] - 3.0) ** 2, x0=np.array([0.5, 0.5]), inequality_constraints=[lambda x: x[0] - 1.0, lambda x: x[1] - 1.0], name="corner", ) barrier = barrier_method(problem) ``` ```{python} #| label: fig-central-path #| fig-cap: "The central path for the corner problem. Each marker is the minimizer of one barrier problem; as t grows they move from the interior toward the corner (1, 1), reaching it only in the limit." from optimlab.viz import central_path_figure fig = central_path_figure(problem, barrier) fig.update_layout(title="The central path") fig.show() ``` ```{python} #| label: fig-duality-gap #| fig-cap: "The duality gap m/t after each outer iteration, on a log scale. Multiplying t by a fixed factor each round divides the gap by the same factor: a straight line, a geometric rate." from optimlab.viz import duality_gap_figure duality_gap_figure(barrier.grad_norm_trajectory, tol=1e-8).show() ``` The method reaches `{python} f"({barrier.x[0]:.6f}, {barrier.x[1]:.6f})"` in `{python} barrier.n_iter` outer iterations. At that point the multipliers can be read off the formula above, $\lambda_i = 1/(-t\, g_i(x))$, and stationarity can be checked directly by drawing the vectors: ```{python} #| label: fig-kkt-geometry #| fig-cap: "KKT stationarity at the corner optimum: -∇f(x*) and the two active constraints' gradients, each scaled by its multiplier. The scaled constraint gradients add up (dashed) to exactly -∇f. The arrows are drawn with a small sideways offset only so that they do not overlap." from optimlab.viz import kkt_geometry_figure kkt_geometry_figure(problem, barrier).show() ``` Both multipliers are positive (both constraints are holding the optimum in place) and, by the problem's symmetry, equal. [Implementation: [`optimizers/barrier_method.py`]({{< var src >}}/optimizers/barrier_method.py) · Checked by: [`test_barrier_method.py`]({{< var tests >}}/test_barrier_method.py)]{.impl} ## ADMM: splitting a problem into easy pieces The barrier method needs derivatives of everything. Many useful problems are sums of pieces that are each easy to handle on their own but awkward together: a smooth data-fit term plus a nonsmooth penalty, or terms that live on different machines. The **alternating direction method of multipliers** [@boyd2011distributed] handles such sums by giving each piece its own copy of the variable and using a multiplier to make the copies agree. ### From the augmented Lagrangian to ADMM Write $\min_x f(x) + h(x)$ as $\min_{x, z} f(x) + h(z)$ subject to $x = z$. The **augmented Lagrangian** adds a quadratic penalty to the ordinary one: $$ \mathcal{L}_\rho(x, z, \nu) = f(x) + h(z) + \nu^\top (x - z) + \frac{\rho}{2}\lVert x - z \rVert^2 . $$ The classical *method of multipliers* minimizes $\mathcal{L}_\rho$ over $(x, z)$ jointly, then moves the multiplier uphill on the dual, $\nu \leftarrow \nu + \rho(x - z)$. ADMM makes one change: instead of minimizing over $x$ and $z$ jointly, which is as hard as the original problem, it minimizes over $x$ and then over $z$, once each. With the scaled multiplier $u = \nu/\rho$, both minimizations become proximal steps (@sec-least-squares): ::: {.callout-note title="Definition — ADMM (scaled form)"} $$ \begin{aligned} x_{k+1} &= \operatorname{prox}_{f/\rho}\big(z_k - u_k\big), \\ z_{k+1} &= \operatorname{prox}_{h/\rho}\big(x_{k+1} + u_k\big), \\ u_{k+1} &= u_k + x_{k+1} - z_{k+1} . \end{aligned} $$ ::: Each step touches only one piece of the objective, and neither needs to be differentiable. For convex $f$ and $h$, ADMM converges for every $\rho > 0$. The penalty $\rho$ is a speed setting, not a stability condition. ### Stopping: two residuals, not one Two things must go to zero. The **primal residual** $\lVert x_k - z_k \rVert$ measures whether the copies agree. The **dual residual** $\rho \lVert z_k - z_{k-1} \rVert$ measures whether $z$, and with it the multiplier, has stopped moving. It is tempting to stop on the primal residual alone. That is unsafe: $x$ and $z$ can agree by coincidence in the middle of a run, long before either is correct. On a lasso problem with $\rho = 50$ it happens within a few iterations: ```{python} #| label: admm-setup import plotly.graph_objects as go from optimlab.optimizers.admm import ADMMProblem, admm from optimlab.optimizers.proximal_gradient import CompositeProblem, proximal_gradient, soft_threshold from optimlab.viz.theme import contrasting_categorical, layout_template rng = np.random.default_rng(0) m, n = 50, 20 A = rng.standard_normal((m, n)) x_true = np.zeros(n) x_true[[2, 5, 9]] = [3.0, -2.0, 1.5] b = A @ x_true + 0.01 * rng.standard_normal(m) alpha = 1.0 AtA, Atb = A.T @ A, A.T @ b def lasso_admm(rho): return ADMMProblem( prox_f=lambda v, t: np.linalg.solve(AtA + np.eye(n) / t, Atb + v / t), prox_g=lambda v, t: soft_threshold(v, alpha * t), x0=np.zeros(n), f_obj=lambda x: 0.5 * np.sum((A @ x - b) ** 2), g_obj=lambda x: alpha * np.sum(np.abs(x)), ) pg = proximal_gradient( CompositeProblem(grad_smooth=lambda x: A.T @ (A @ x - b), prox_nonsmooth=lambda v, t: soft_threshold(v, alpha * t), x0=np.zeros(n)), lr=1.0 / np.linalg.eigvalsh(AtA).max(), max_iter=5000, tol=1e-10, ) f_star = 0.5 * np.sum((A @ pg.x - b) ** 2) + alpha * np.sum(np.abs(pg.x)) run = admm(lasso_admm(50.0), rho=50.0, max_iter=5000, tol=1e-6, record_history=True) primal = np.array(run.history["primal_residual"]) dual_res = np.array(run.history["dual_residual"]) first_dip = int(np.argmax(primal < 1e-6)) + 1 gap_at_dip = (run.f_trajectory[first_dip] - f_star) / f_star ``` ```{python} #| label: fig-admm-residuals #| fig-cap: "ADMM on a lasso problem with ρ = 50 (log scale). The primal residual collapses to numerical zero early, a coincidence of x and z agreeing, while the dual residual and the objective's distance from the optimum are still large. The primal residual then rises again. Only when both residuals are small is the run actually finished." colors = contrasting_categorical() iters = np.arange(1, run.n_iter + 1) fig = go.Figure([ go.Scatter(x=iters, y=np.maximum(primal, 1e-17), name="primal residual ‖x − z‖", line={"color": colors[0], "width": 2}), go.Scatter(x=iters, y=np.maximum(dual_res, 1e-17), name="dual residual ρ‖zₖ − zₖ₋₁‖", line={"color": colors[1], "width": 2}), go.Scatter(x=iters, y=np.maximum(np.array(run.f_trajectory[1:]) - f_star, 1e-17), name="objective − optimum", line={"color": colors[2], "width": 2, "dash": "dot"}), ]) fig.update_layout(**layout_template(xaxis_title="iteration", yaxis_title="value", yaxis={"type": "log", "exponentformat": "power"}, legend={"y": 1, "x": 1.02, "xanchor": "left"}, margin={"r": 230})) fig.show() ``` At iteration `{python} first_dip` the primal residual is below $10^{-6}$ while the objective is still `{python} f"{100 * gap_at_dip:.0f}"`% above its optimal value. A primal-only rule would have stopped there. The standard rule requires both residuals to be small [@boyd2011distributed, §3.3.1], and the dual residual is far from small at that point. ### ADMM and proximal gradient agree The lasso was solved in @sec-least-squares by proximal gradient, which takes an explicit gradient step. ADMM never takes one, yet the two land on the same answer: ```{python} #| label: tbl-admm-rho #| tbl-cap: "ADMM on the same lasso problem for several values of ρ. All converge to the proximal-gradient solution; ρ only changes how many iterations it takes." rows = [] for rho in [0.5, 1.0, 5.0, 20.0, 50.0]: r = admm(lasso_admm(rho), rho=rho, max_iter=5000, tol=1e-6) rows.append((rho, r.n_iter, f"{r.f:.4f}", f"{np.max(np.abs(r.x - pg.x)):.1e}")) markdown_table(["ρ", "iterations", "objective", "max difference from proximal gradient"], rows) ``` [Implementation: [`optimizers/admm.py`]({{< var src >}}/optimizers/admm.py) · Checked by: [`test_admm.py`]({{< var tests >}}/test_admm.py)]{.impl} ## When the unknown is a function: calculus of variations So far the unknown has been a vector. In the calculus of variations it is a whole function $y(x)$, and the objective is an integral, a *functional*: $$ J[y] = \int_a^b F\big(x, y(x), y'(x)\big)\, dx, \qquad y(a), y(b) \text{ fixed}. $$ The first-order condition is derived the same way as $\nabla f = 0$: perturb the candidate and require the first-order change to vanish. Let $y + \epsilon\eta$ be a perturbed curve, with $\eta(a) = \eta(b) = 0$ so the endpoints stay fixed. Then $$ \frac{d}{d\epsilon} J[y + \epsilon \eta]\Big|_{\epsilon = 0} = \int_a^b \big( F_y\, \eta + F_{y'}\, \eta' \big)\, dx = \int_a^b \Big( F_y - \frac{d}{dx} F_{y'} \Big)\, \eta\, dx , $$ after integrating the second term by parts (the boundary term vanishes because $\eta$ does at both ends). For this to be zero for *every* perturbation $\eta$, the bracket must vanish everywhere: ::: {.callout-note title="Definition — Euler–Lagrange equation"} $$ \frac{\partial F}{\partial y} - \frac{d}{dx}\frac{\partial F}{\partial y'} = 0 . $$ ::: It plays exactly the role of $\nabla f(x^\star) = 0$, now as a differential equation [@gelfand2000calculus]. A classic instance is the soap film stretched between two coaxial rings of radius $y_0$ a distance $\ell$ apart. The film minimizes its area, $J[y] = 2\pi \int_0^\ell y \sqrt{1 + y'^2}\, dx$. Because $F$ does not depend on $x$ explicitly, the Euler–Lagrange equation has a first integral, $F - y' F_{y'} = \text{const}$, which here reads $y / \sqrt{1 + y'^2} = c$, and its solution is a **catenary**, $$ y(x) = c \cosh\!\Big(\frac{x - \ell/2}{c}\Big), \qquad y_0 = c\cosh\!\Big(\frac{\ell}{2c}\Big). $$ The same problem can be solved without any of this calculus. Represent $y$ by its values on a grid of $N + 1$ points, replace the integral with a sum, and the functional becomes an ordinary function of $N - 1$ unknowns (the interior heights), which BFGS (@sec-gradient-methods) can minimize: ```{python} #| label: fig-catenary #| fig-cap: "Discretize, then optimize: starting from a flat film, BFGS on 39 interior heights bows the curve inward until it lies on the analytic catenary from the Euler–Lagrange equation." import jax.numpy as jnp from scipy.optimize import brentq from optimlab.core import Problem from optimlab.optimizers import bfgs y0, ell, N = 1.0, 1.0, 40 xs = np.linspace(0.0, ell, N + 1) dx = xs[1] - xs[0] def area(y_inner): # the 2π factor does not change the minimizer y = jnp.concatenate([jnp.array([y0]), y_inner, jnp.array([y0])]) return jnp.sum((y[:-1] + y[1:]) / 2 * jnp.sqrt(1.0 + (jnp.diff(y) / dx) ** 2)) * dx film = bfgs(Problem(f=area, x0=np.full(N - 1, y0)), max_iter=500) y_opt = np.concatenate([[y0], film.x, [y0]]) c = brentq(lambda c: c * np.cosh(ell / (2 * c)) - y0, 0.5, 3.0) # the stable (larger) root y_catenary = c * np.cosh((xs - ell / 2) / c) max_gap = np.max(np.abs(y_opt - y_catenary)) colors = contrasting_categorical() fig = go.Figure([ go.Scatter(x=xs, y=np.full(N + 1, y0), mode="lines", name="flat starting guess", line={"color": "#898781", "width": 1.5, "dash": "dot"}), go.Scatter(x=xs, y=y_catenary, mode="lines", name="catenary (Euler–Lagrange)", line={"color": "#52514e", "width": 4}), go.Scatter(x=xs, y=y_opt, mode="markers", name="BFGS on the discretized area", marker={"size": 6, "color": colors[0]}), ]) fig.update_layout(**layout_template(title="A soap film between two rings", xaxis_title="x", yaxis_title="radius y(x)")) fig.show() ``` The largest gap between the optimized heights and the catenary is `{python} f"{max_gap:.1e}"`. Solving the differential equation and minimizing a fine discretization of the functional are, in the limit, the same computation. Direct discretization is how most optimal-control problems are solved in practice, as @sec-applications shows. ## What these examples cannot show - **Scale.** Real interior-point solvers exploit sparsity in the Newton systems and handle equality constraints and infeasible starting points; the two-variable example shows the idea, not the engineering. - **Non-convex constraints.** The KKT conditions remain necessary, but duality gaps can be positive, and a KKT point need not be a global optimum. - **ADMM's speed.** ADMM reaches moderate accuracy quickly and high accuracy slowly. Its advantages, decomposition and parallelism, appear only on problems far larger than these. ## Connections - **@sec-geometry**: the KKT conditions generalize $\nabla f(x^\star) = 0$, and convexity is what makes them sufficient. - **@sec-linear-programming**: duality gives every linear program a twin, and its variables are the constraints' shadow prices. - **@sec-least-squares**: proximal operators are ADMM's building blocks, and box constraints are the simplest constraints of all. - **@sec-applications** applies the KKT system to portfolio selection, the barrier method to fair resource allocation, and discretize-then-optimize to optimal control.