Convex Optimization, in depth: recognising convexity, the LP to SDP ladder, proximal and Newton methods with certificates

Why convexity is the line that matters

A set C is convex if the segment between any two of its points stays inside it. A function f is convex if its domain is convex and f(θx + (1−θ)y) ≤ θf(x) + (1−θ)f(y) for θ in [0, 1]. In words, the chord lies on or above the graph. A convex optimisation problem minimises a convex function subject to convex inequality constraints g_i(x) ≤ 0 and affine equality constraints Ax = b.

The central fact has a short proof. Suppose x is a local minimum but some feasible y has f(y) < f(x). The segment from x to y is feasible, and by the chord inequality f on it is at most θf(y) + (1−θ)f(x), strictly below f(x) for every θ > 0. So better points exist arbitrarily close to x, a contradiction. Convexity removes the question every nonconvex practitioner lives with: is this the answer, or just an answer?

Optimality can also be checked: for differentiable f, x is optimal exactly when ∇f(x)^T(y − x) ≥ 0 for every feasible y, which is ∇f(x) = 0 without constraints. And solvers can return a lower bound with the answer. The difference, the duality gap, is a certificate, and the worked example below stops on one.

Recognising convex functions

Proving convexity from the definition is slow. In practice you combine a small set of known convex atoms using rules that preserve convexity. For twice-differentiable functions there is also a direct test: f is convex exactly when its Hessian is positive semidefinite everywhere on the domain. For example, the Hessian of ½‖Ax − b‖² is A^TA, which is always positive semidefinite.

Rules that preserve convexity:

• Nonnegative weighted sums of convex functions are convex.
• Composition with an affine map, f(Ax + b), keeps convexity. This is why every regression loss of a linear model is convex in the weights.
• The pointwise maximum of convex functions is convex. Hinge loss max(0, 1 − y w^Tx) is one example.
• Composing h(g(x)) is convex when h is convex and nondecreasing and g is convex. For example, exp of a convex function is convex.
• Intersections of convex sets are convex, so you can stack constraints freely.

The problem ladder: LP, QP, SOCP, SDP

Solvers do not take arbitrary convex functions. They take problems in one of a few standard forms, each a special case of the next. A linear program (LP) has a linear objective and linear constraints. A convex quadratic program (QP) adds a positive semidefinite quadratic objective. A second-order cone program (SOCP) allows constraints of the form ‖Ax + b‖₂ ≤ c^Tx + d, which covers robust and norm-bounded problems. A semidefinite program (SDP) constrains a symmetric matrix variable to be positive semidefinite.

Model your problem in the lowest class that fits. Inner-class solvers are faster, more accurate and scale further: an LP with millions of variables is routine, while a 5,000 by 5,000 SDP variable already has about 12.5 million free entries.

Modelling with disciplined convex programming

Modelling languages such as CVXPY, CVX and Convex.jl implement disciplined convex programming (DCP): they track curvature through every operation, accept the problem only if convexity follows from the rules, and rewrite it as a cone program. A long-only minimum-variance portfolio with a return target:

import cvxpy as cp
import numpy as np

mu = np.array([0.08, 0.12, 0.10, 0.07])          # expected annual returns
Sigma = np.array([[0.10, 0.02, 0.04, 0.00],
                  [0.02, 0.20, 0.05, 0.01],
                  [0.04, 0.05, 0.15, 0.02],
                  [0.00, 0.01, 0.02, 0.05]])       # covariance, must be PSD

w = cp.Variable(4)
risk = cp.quad_form(w, Sigma)
budget = cp.sum(w) == 1
target = mu @ w >= 0.09
prob = cp.Problem(cp.Minimize(risk), [budget, target, w >= 0])
prob.solve()

print(prob.status)          # check this before trusting w.value
print(w.value, risk.value)
print(target.dual_value)    # marginal variance per unit of required return

Three habits matter here. Always read prob.status before using the answer. The values infeasible, unbounded and optimal_inaccurate all appear in practice. Read the dual value of a constraint you care about: here it prices one more unit of required return in variance. And when DCP rejects a model, either the problem is not convex or its convexity cannot be proven by the rules. cp.sqrt(x)**2 is rejected even though it equals x, while cp.norm(x, 2) is accepted where cp.sqrt(cp.sum_squares(x)) is not. Rewrite using the atom that states your intent.

First-order and proximal methods

First-order methods use only gradients, so each iteration is cheap: one or two matrix-vector products. If f is L-smooth (its gradient is L-Lipschitz), gradient descent with step 1/L reaches f(x_k) − f* ≤ L‖x₀ − x*‖²/(2k). If f is also μ-strongly convex, the error shrinks by about a factor (1 − μ/L) each step, so the condition number κ = L/μ sets the speed. Nesterov acceleration improves the general rate to O(1/k²) and the strongly convex dependence to roughly √κ.

Plain gradients fail on nonsmooth terms such as the ℓ₁ penalty that produces sparse solutions. Proximal gradient takes a gradient step on the smooth part, then applies the closed-form proximal operator of the nonsmooth part. For λ‖x‖₁ that is soft-thresholding: shrink each coordinate toward zero by λ/L, clipping at zero. That is ISTA; adding Nesterov momentum gives FISTA.

import numpy as np

def soft_threshold(v, t):
    return np.sign(v) * np.maximum(np.abs(v) - t, 0.0)

def lasso_objective(A, b, lam, x):
    return 0.5 * np.sum((A @ x - b) ** 2) + lam * np.sum(np.abs(x))

def lasso_gap(A, b, lam, x):
    """Duality gap: a certificate, not a heuristic."""
    r = A @ x - b
    scale = min(1.0, lam / np.max(np.abs(A.T @ r)))
    nu = scale * r                        # dual-feasible point
    dual = -0.5 * nu @ nu - nu @ b
    return lasso_objective(A, b, lam, x) - dual

def lasso_fista(A, b, lam, iters=20000, tol=1e-6, accelerate=True):
    """minimise 0.5*||Ax - b||^2 + lam*||x||_1"""
    L = np.linalg.norm(A, 2) ** 2          # Lipschitz constant of the smooth part
    x = np.zeros(A.shape[1]); y = x.copy(); t = 1.0
    for k in range(1, iters + 1):
        grad = A.T @ (A @ y - b)
        x_new = soft_threshold(y - grad / L, lam / L)
        if accelerate:
            t_new = (1 + np.sqrt(1 + 4 * t * t)) / 2
            y = x_new + ((t - 1) / t_new) * (x_new - x)
            t = t_new
        else:
            y = x_new
        if k % 10 == 0 and lasso_gap(A, b, lam, x_new) <= tol:
            return x_new, k
        x = x_new
    return x, iters

The gap function is where duality earns its place. The lasso dual is to maximise −½‖ν‖² − ν^Tb subject to ‖A^Tν‖_∞ ≤ λ, and at the optimum ν equals the residual Ax − b. Scaling the current residual into the constraint gives a feasible dual point whose value bounds the optimum from below. A gap under 10⁻⁶ proves the objective is within 10⁻⁶ of the best possible.

Worked example: lasso and conditioning

To see what conditioning does, run the code on a 200 by 500 lasso with a 10-sparse true signal and small noise. Make the columns of A share a common component with correlation ρ, and stop when the duality gap falls below 10⁻⁶, checked every ten iterations.

Conditioning dominates: correlated columns cost two orders of magnitude more iterations with the same code. Acceleration is not a free win either: on the easy problem FISTA took more iterations than ISTA. FISTA is not monotone, and a common cause is momentum overshooting once the set of nonzero coordinates has settled. Practical implementations add adaptive restarts that reset the momentum whenever the objective goes up. And a smaller λ means a less sparse, harder problem. If you need many values of λ, solve from large to small and warm-start each solve from the previous answer. That technique, the regularisation path, is how production lasso solvers stay fast.

Newton and interior-point methods

Newton's method solves ∇²f(x) Δx = −∇f(x) and backtracks along Δx until the objective drops enough. Near the optimum convergence is quadratic, roughly doubling the correct digits each iteration. The natural stopping test is the Newton decrement λ² = −∇f^TΔx, since λ²/2 estimates the remaining suboptimality.

def logistic_newton(X, y, mu=1e-3, tol=1e-10, max_iter=50):
    """minimise mean log(1 + exp(-y * Xw)) + mu/2 ||w||^2, labels y in {-1, +1}"""
    m, n = X.shape
    w = np.zeros(n)
    f = lambda w: np.mean(np.logaddexp(0, -y * (X @ w))) + 0.5 * mu * w @ w
    for k in range(max_iter):
        s = 1 / (1 + np.exp(y * (X @ w)))            # sigmoid(-margin)
        g = -(X.T @ (y * s)) / m + mu * w
        H = (X.T * (s * (1 - s))) @ X / m + mu * np.eye(n)
        step = np.linalg.solve(H, -g)
        if -g @ step / 2 <= tol:                       # half the Newton decrement squared
            return w, k
        t, fw = 1.0, f(w)
        while f(w + t * step) > fw + 0.25 * t * (g @ step):   # Armijo backtracking
            t *= 0.5
        w = w + t * step
    return w, max_iter

On a 1,000 by 20 synthetic logistic regression this stopped after 7 iterations. Gradient descent with step 1/L needed 871 iterations just to bring the gradient norm below 10⁻⁶. The price is per iteration: forming and factoring the Hessian costs O(mn² + n³), which is trivial at n = 20 and impossible at n = 10⁶.

Interior-point methods extend Newton to constraints by replacing each inequality with a logarithmic barrier and following the central path as the barrier weight shrinks. They typically converge in a few tens of Newton steps, which is why most conic solvers behind modelling tools use them. Their limit is the sparse factorisation each step needs.

Choosing a method

For repeated solves, declare cp.Parameter objects so CVXPY canonicalises the problem once and re-solves it with new data.

Failure modes

• A nonconvex term in disguise. A product of two variables x·y, a ratio, or a quadratic with an indefinite matrix makes the problem nonconvex. Every guarantee above is then gone. DCP catches this; hand-written solvers do not.
• Step size above 1/L. Gradient and proximal methods can diverge or oscillate. Estimate L with a power iteration, or use backtracking.
• Bad scaling. Features in very different units inflate κ. Standardise columns first.
• Trusting optimal_inaccurate. This status means the solver stopped short of its tolerance. Rescale the data, try a different solver, or tighten settings before you ship the answer.
• Infeasible by accident. Contradictory constraints, such as a return target above the best single asset, produce infeasible. Relax constraints one at a time, or add slack variables with a penalty to see which one binds.
• Stopping on step size. Small steps can mean slow progress, not optimality. Stop on a duality gap or Newton decrement when one is available.

What to do next

1. Take one loss you already minimise and prove it is convex with the atom and rule tables above, or find the term that breaks convexity.
2. Run the lasso code on your own data with both settings of accelerate, stopping on the duality gap. Then add a restart rule and measure again.
3. Rewrite one hand-tuned heuristic in your system, such as an allocation or schedule, as a CVXPY model, and compare objective values.
4. Read the dual values of your binding constraints and write down what each one prices.
5. Learn why the gap closes in Lagrangian duality and how LPs specialise it in LP duality.
6. Compare against the simplex method in the simplex article, and see why nonconvex training behaves differently in the optimisation landscape article and gradient descent analysis.