Projected and Proximal Gradient Methods for Composite Objectives

Two kinds of structure defeat plain gradient descent: a hard constraint that the solution must satisfy exactly, and a term in the objective that has no gradient at the point where the optimum sits. Both are handled by the same manoeuvre, taking a gradient step with respect to the differentiable part and then applying an operator that enforces the remaining structure, but the operators are different objects with different computational profiles. Both converge under conditions that are close to the descent lemma’s, and both fail in ways that are specific to the composition rather than to the gradient step.

Clipping a step is not projecting it

A tempting implementation of “stay inside the box [0, 1]ⁿ” is to compute the unconstrained step and then clamp each coordinate into the box. The result is a feasible point, but not the projected gradient step. Clamping the coordinates of x − η∇f(x) is the projection of the gradient step only when the box is axis-aligned and the projection is applied as the composition of a step and a projection, which it is: x⁺ = Π_C(x − η∇f(x)). What differs is the direction. The projected point is not generally x minus a step along any direction derived from the gradient alone; the descent direction is the vector x − Π_C(x − ηg) rescaled, which is the gradient projected onto the tangent cone of C at x. Treating this as a step of clipped gradient and then applying a line search to the underlying gradient gives a direction that fails the descent test, and the accepted steps stop decreasing the objective.

The projection itself is cheap for the sets that arise most often. A box or a halfspace costs one pass over the coordinates, so O(n). The Euclidean ball of radius r costs one norm evaluation and one scaling. The probability simplex, defined by non-negative entries that sum to one, costs a sort, so O(n log n) naively, or O(n) expected with a linear-time selection algorithm for the threshold value; the projection has the closed form xᵢ = max(vᵢ − θ, 0) where θ is chosen so the entries sum to one. The ℓ₁ ball of radius r reduces to the simplex case by a sign flip and a translation. Sets with polyhedral constraints expressed as Ax ≤ b for a general matrix A require a quadratic program per projection, which changes the character of the algorithm entirely.

Convergence of the projected iteration

With f convex and L-smooth and C convex and closed, the projected iteration with η ≤ 1/L satisfies

f(x_{k+1}) ≤ f(x_k) − (1/(2η))‖x_{k+1} − x_k‖²,

an exact analogue of the unconstrained decrease bound with the projected step length in place of the gradient norm. Convexity then yields O(1/k) in function value, and μ-strong convexity yields the same (1 − μ/L) contraction per iteration as the unconstrained case. Three requirements are doing the work: C must be convex, so that the projection is unique and the tangent-cone characterisation holds; the projection must be exact, so that a round-off-level error is the only error; and the step size must be bounded by the smoothness of f alone, with the constraint contributing nothing to L. That last point is a frequent source of confusion, because a heavily constrained problem can have a global curvature far larger than that of f on the whole space, and the relevant constant is still the one for f.

Inexact projection changes the result qualitatively. If each projection has error bounded by δ, the iteration converges to a neighbourhood of the constrained optimum whose radius is proportional to δ/η. Reducing the step size shrinks that neighbourhood at the cost of slower progress, which is the same trade as the noise floor under stochastic gradients. Iterative projections, such as alternating projections onto two sets, fall into this category and usually need the tolerance tied to the step size rather than fixed.

The proximal operator and composite objectives

The second structure is a term in the objective that has no gradient where the optimum is. Write the objective as f(x) + g(x) with f smooth and g convex and possibly non-smooth: an ℓ₁ norm, an indicator function of a convex set, a total-variation term. The proximal operator of g with step η is

prox_{ηg}(v) = argmin_u { g(u) + (1/(2η))‖u − v‖² },

and the composite iteration is x⁺ = prox_{ηg}(x − η∇f(x)). The indicator function of a convex set reproduces projected gradient descent exactly, since its proximal operator is the projection. For the ℓ₁ norm the proximal operator is soft thresholding, prox_{ηλ‖·‖₁}(v)ᵢ = sign(vᵢ)·max(|vᵢ| − ηλ, 0), a single pass over the coordinates. That formula is the reason ℓ₁-regularised problems are tractable: the penalty removes a fixed amount from every coefficient and sends small ones exactly to zero, so the sparsity is produced by an inexpensive element-wise operation rather than by a search over active sets.

def soft_threshold(lam, v):
    return [0.0 if abs(vi) <= lam else (vi - lam if vi > 0 else vi + lam) for vi in v]

def ista(A, b, lam, eta, iters):
    """Minimise 0.5*||Ax-b||^2 + lam*||x||_1 with A as a list of rows."""
    n = len(A[0])
    x = [0.0] * n
    for _ in range(iters):
        # gradient of the smooth part: A^T (A x - b)
        g = [sum(A[i][j] * (sum(A[i][k] * x[k] for k in range(n)) - b[i])
                 for i in range(len(A))) for j in range(n)]
        x = soft_threshold(eta * lam, [x[j] - eta * g[j] for j in range(n)])
    return x

The resulting algorithm is known as ISTA or the proximal gradient method, and its function value decreases at O(1/k), the same rate as plain gradient descent on a smooth problem. The accelerated version composes the same two operators with the momentum structure of Nesterov’s accelerated scheme and reaches O(1/k²). The composition is not free: the accelerated composite method requires the extrapolation to be applied to the sequence of proximal points, and the resulting objective sequence is not monotone, so implementations that check for an increase and restart recover monotonicity at the cost of losing part of the theoretical rate. The step size again depends only on the smooth part: η ≤ 1/L where L is the smoothness constant of f; for a least-squares data term with design matrix A, that constant is the largest eigenvalue of AᵀA, which can be estimated by power iteration rather than computed.

Choosing between a constraint and a penalty

Framed as f plus an indicator, a hard constraint is exact and non-negotiable. Framed as a penalty λ·max(0, c(x))², the constraint is soft: the returned point violates it by an amount that shrinks as λ grows, while the smoothness constant of the penalised objective grows without bound and the conditioning degrades with it. The two formulations also produce different optima. For a constrained problem with an optimal value p*, the penalised problem returns a point with a smaller objective value and a violated constraint, and no single λ reproduces the constrained solution exactly; matching the KKT conditions of the constrained problem requires the multiplier as well as the penalty, which is the idea behind augmented Lagrangian methods. Proximal methods sit between the two: the penalty is applied through an exact operator rather than added to the smooth part, so ℓ₁-penalised objectives with any λ > 0 have exact zeros in their solutions, whereas adding a smooth surrogate for the same penalty produces solutions that are small but never exactly zero. The λ parameter then plays the role of the multiplier from the KKT conditions, and the trade-off it controls is a genuine modelling choice about sparsity rather than a numerical knob.

Where the composition degrades

Non-convex feasible sets break the projection machinery in a specific way. The projection is still well defined for a closed set, and the iteration still has a stationary-point guarantee if the projection error is small, but the guarantee is stated with respect to the projected-gradient fixed-point condition Π_C(x − η∇f(x)) = x, which on a non-convex set is satisfiable at points that are neither local minima nor even local stationary points of the constrained problem. A closed-form proximal operator is not always available: the proximal operator of the ℓ₁ norm over overlapping groups, for instance, requires solving a small optimization problem per group, and the total cost per iteration can exceed the cost of the gradient step it is meant to accompany. Group and structured penalties also change the sparsity pattern in ways that are easy to get wrong: overlapping groups produce solutions that are not sparse at the group level in any simple sense, and the support of the solution is determined by the subdifferential evaluated at zero, so the model’s structure is a property of the operator rather than of any threshold applied afterwards. Finally, warm-starting across a sequence of decreasing values of λ is almost always necessary for the solver to be competitive, since a cold start at a small λ begins from a dense point far from the sparse solution; the solver’s iteration count at each λ is what determines whether the whole path is affordable.