Week 6: Interior-Point Methods

Why Not Projected Gradient Alone#

Projected gradient descent (intermediate Optim Week 9) works for constrained convex problems but has a fundamental weakness: the projection operation itself can be expensive. For a polyhedron defined by $m$ linear inequalities, projection is a QP that may be as hard as the original problem. For the PSD cone, projection requires an eigendecomposition (as in Week 5's lab) — feasible but $O(n^3)$ per step.

More critically, projected gradient produces iterates on the boundary — it bounces along constraint faces. Boundary iterates mean some constraints are active, which means the gradient may be discontinuous when a constraint joins or leaves the active set. This active-set chattering slows convergence and makes the method unreliable for high-accuracy solutions.

Interior-point methods (IPMs) take the opposite approach: stay strictly inside the feasible region, and gently approach the boundary only in the limit. They solve a sequence of unconstrained (or equality-constrained) smooth approximations, each of which is solved by Newton's method. The result is polynomial-time convergence to high accuracy for LP, QP, SOCP, and SDP.

Logarithmic Barrier#

The main idea of an IPM is to replace each inequality constraint $h_j(x) \leq 0$ with a barrier term that penalises approaching the boundary:

$\phi(x) = -\sum_{j=1}^p \log(-h_j(x))$

As $h_j(x) \to 0^-$ (approaching the boundary from inside), $-\log(-h_j(x)) \to +\infty$ — the barrier explodes, preventing the iterates from ever touching or crossing the boundary. The feasible set's interior is the domain of $\phi$.

The barrier problem for a fixed parameter $t > 0$ is:

$\min_x \; t f(x) + \phi(x) \quad \text{s.t.} \quad g_i(x) = 0$

where the equality constraints remain explicit. For a pure inequality-constrained problem:

$\min_x \; t f(x) - \sum_{j=1}^p \log(-h_j(x))$

The parameter $t$ controls the trade-off: when $t$ is small, the barrier dominates and the minimiser is forced toward the "analytic centre" of the feasible set (the point that maximises the product of distances to all constraint boundaries). When $t$ is large, the $t f(x)$ term dominates and the minimiser approaches the true constrained optimum.

Worked 1D example. $\min \; x$ subject to $x \geq 0$ (i.e., $h(x) = -x \leq 0$).

Barrier: $\phi(x) = -\log(x)$ for $x > 0$. The barrier problem is $\min_{x > 0} \; t x - \log(x)$.

Set derivative to zero: $t - 1/x = 0 \Rightarrow x^*(t) = 1/t$. As $t \to \infty$, $x^*(t) \to 0$, the true optimum. The path $x^*(t) = 1/t$ for $t > 0$ is the central path for this trivial problem.

Central Path#

The central path is the set $\{x^*(t) \mid t > 0\}$ where $x^*(t)$ is the minimiser of the barrier problem with parameter $t$. Under mild conditions (convexity, Slater, self-concordance), the central path is a smooth curve that converges to the primal optimum $x^*$ as $t \to \infty$. Simultaneously, the dual variables $\lambda_j^*(t) = 1/(t \cdot (-h_j(x^*(t))))$ converge to the optimal dual variables $\lambda_j^*$.

Central path in KKT terms. The barrier subproblem's KKT conditions give a perturbed complementary slackness:

$\lambda_j(t) \cdot h_j(x(t)) = -1/t$

As $t \to \infty$, the right-hand side $\to 0$, recovering exact complementary slackness. The IPM follows the central path by solving a sequence of Newton systems that track this perturbation toward zero.

Primal-dual IPMs. Modern IPMs solve the primal and dual simultaneously, taking Newton steps on the perturbed KKT system. At each iteration, the barrier parameter $t$ is increased (e.g., $t \leftarrow \mu t$ with $\mu > 1$), and one or a few Newton steps solve the updated subproblem. The method terminates when the duality gap $m/t$ (where $m$ is the number of inequalities) is below tolerance — typically $10^{-8}$ or tighter.

Newton on the Barrier#

Why does Newton work well on the barrier subproblem? The log barrier has a remarkable property called self-concordance: its third derivative is bounded by a constant times its second derivative to the power $3/2$:

$|\phi'''(x)| \leq 2 (\phi''(x))^{3/2}$

This property, discovered by Nesterov and Nemirovskii, ensures that Newton's method on self-concordant functions converges in a predictable number of iterations independent of the problem conditioning — the Newton decrement provides a reliable step-size and stopping criterion. This is the mathematical engine behind IPM's polynomial-time guarantees.

Practical Newton step for IPM. The barrier gradient and Hessian are:

$\nabla \phi(x) = \sum_j \frac{1}{-h_j(x)} \nabla h_j(x)$ $\nabla^2 \phi(x) = \sum_j \frac{1}{(-h_j(x))^2} \nabla h_j(x) \nabla h_j(x)^T + \frac{1}{-h_j(x)} \nabla^2 h_j(x)$

The Newton direction solves $(\nabla^2 \phi - \text{equality projections}) \Delta x = -\nabla \phi$. For LP ($\nabla^2 h_j = 0$), the Hessian simplifies and the system can be reduced to the normal equations form $A D A^T$, where $D$ is a diagonal scaling derived from the barrier.

Solver Practice#

When you call a production IPM solver (e.g., commercial LP/QP/SDP solvers or open-source equivalents), you get:

1. Approximate primal solution $x$: feasible to within tolerance.
2. Approximate dual solution $(\lambda, \nu)$: provides sensitivity/shadow price information.
3. Duality gap: $p^* - d^*$ (or its estimate), converging to zero.
4. Status flags: Optimal, Infeasible, Unbounded, or numerical difficulties.

IPM vs first-order for ML. IPMs excel for moderate-sized problems (thousands to low millions of variables) where high accuracy matters. For large-scale ML training with millions of parameters and billions of data points, first-order methods (SGD/AdamAdaptive Moment Estimation (optimizer)) dominate because:

• IPM requires solving a linear system at each iteration — $O(d^3)$ for dense, intractable for $d = 10^7$ deep net parameters.
• ML objectives are typically nonconvex, where IPM's global guarantees do not apply.
• High accuracy (8+ digits) is irrelevant when generalisation is the real goal.

IPMs are the tool for structured, moderate-scale convex problems — the "structure path" in Week 14's decision flow.

Browser lab#

Consider $\min_{x,y} \; x + y$ subject to $x \geq 0, y \geq 0, x + y \geq 1$. Trace the barrier subproblem minimisers for increasing $t$ via gradient descent on the barrier objective, and visualise the central path.

import numpy as np
import matplotlib.pyplot as plt

def barrier_obj(xy, t):
    """t*(x + y) - log(x) - log(y) - log(x + y - 1)"""
    x, y = xy
    if x <= 0 or y <= 0 or x + y <= 1:
        return np.inf
    return t * (x + y) - np.log(x) - np.log(y) - np.log(x + y - 1)

def barrier_grad(xy, t):
    """Gradient of the barrier objective."""
    x, y = xy
    if x <= 0 or y <= 0 or x + y <= 1:
        return np.array([np.inf, np.inf])
    d_barrier_x = -1/x - 1/(x + y - 1)
    d_barrier_y = -1/y - 1/(x + y - 1)
    return np.array([t + d_barrier_x, t + d_barrier_y])

# Trace central path for increasing t
t_values = [0.5, 1, 2, 4, 8, 16, 32, 64]
path = []

for t in t_values:
    # Minimise barrier via GD from a feasible start
    xy = np.array([2.0, 2.0])  # feasible: x>0, y>0, x+y>1
    for _ in range(2000):
        g = barrier_grad(xy, t)
        step = 0.01 * xy  # heuristic: step size proportional to position
        # Clamp to avoid crossing boundary
        xy_new = xy - np.minimum(step, xy * 0.5) * np.sign(g) * np.minimum(np.abs(g), 1.0)
        # Simple: just do GD with small fixed step after checking domain
        xy_new_feasible = xy - 0.005 * g / (np.linalg.norm(g) + 1e-8)
        if np.all(xy_new_feasible > 0) and (xy_new_feasible[0] + xy_new_feasible[1] > 1):
            xy = xy_new_feasible
        else:
            xy = xy - 0.001 * g / (np.linalg.norm(g) + 1e-8)
    path.append(xy.copy())
    print(f"t={t:4.0f}: minimiser = ({xy[0]:.4f}, {xy[1]:.4f}), obj = {xy[0]+xy[1]:.4f}")

path = np.array(path)

# Plot the central path on the feasible set
fig, ax = plt.subplots(figsize=(7, 6))
# Feasible region: x >= 0, y >= 0, x + y >= 1
x_plot = np.linspace(-0.1, 3, 300)
y_plot = np.linspace(-0.1, 3, 300)
X, Y = np.meshgrid(x_plot, y_plot)
feasible = (X >= 0) & (Y >= 0) & (X + Y >= 1)
ax.contourf(X, Y, feasible.astype(float), levels=1, colors=["lightgreen"], alpha=0.3)

# Plot objective contours
Z = X + Y
ax.contour(X, Y, Z, levels=8, colors="gray", alpha=0.4, linestyles="--")

ax.plot(path[:, 0], path[:, 1], "o-", color="tab:red", linewidth=2, markersize=6, label="Central path")
for i, t in enumerate(t_values):
    ax.annotate(f"t={t}", (path[i, 0] + 0.05, path[i, 1] + 0.05), fontsize=8)

# True optimum: x* = 0.5, y* = 0.5 (symmetric problem)
ax.plot(0.5, 0.5, "*", color="black", markersize=14, label="Optimum (0.5, 0.5)")

ax.set_xlim(-0.1, 3)
ax.set_ylim(-0.1, 3)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("Central path for $\min x+y$ s.t. $x,y \geq 0,\; x+y \geq 1$")
ax.legend()
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)
plt.show()