Convex Optimization

Why Convex Optimization Matters

Convex optimization chooses the best feasible parameters for an objective whose graph has no hidden non-global local minima. Examples include constrained least squares, sparse linear models, and resource allocation over a convex feasible set. The useful structure is mathematical: convexity connects local optimality conditions to global solutions. It does not mean every step rule is fast, every feasible point is good, or every numerical output is certified optimal.

This lesson covers four complementary methods. Subgradient descent handles convex kinks using supporting slopes. Projected gradient handles a smooth objective with a constraint set that is easy to project onto. Coordinate descent changes one coordinate or block at a time. Proximal gradient separates a smooth loss from a nonsmooth term whose proximal subproblem is easy.

The programs use a small diagonal quadratic, a box, and an optional L1 penalty. This lets us trace all four methods and derive the exact solution independently. Projected gradient is run on the smooth zero-penalty member of the family; the other three can handle the nonsmooth penalty. Their output labels both the objective and the update unit, so different losses and different amounts of derivative work are not mistaken for a performance ranking.

Core Concepts and Terminology

A set $C$ is convex when every segment between its points lies in the set. A function $\Phi$ is convex when, for $\mathbf x,\mathbf z$ in its domain and $0\leq\theta\leq1$,

$$\Phi(\theta\mathbf x+(1-\theta)\mathbf z)\leq\theta\Phi(\mathbf x)+(1-\theta)\Phi(\mathbf z).$$

We study $\min_{\mathbf x\in C}\Phi(\mathbf x)$, with $\Phi=f+r$. The smooth part $f$ has an $L$-Lipschitz gradient when $\lVert\nabla f(\mathbf x)-\nabla f(\mathbf z)\rVert_2\leq L\lVert\mathbf x-\mathbf z\rVert_2$. The nonsmooth part $r$ may have kinks. Our finite, nonempty box is $C=\prod_{j=1}^d[\ell_j,u_j]$.

The reference family is

$$f(\mathbf x)=\frac12\sum_{j=1}^d a_j(x_j-c_j)^2,\qquad r(\mathbf x)=\lambda\lVert\mathbf x\rVert_1,\qquad a_j\gt0,\quad\lambda\geq0.$$

Here $(\nabla f(\mathbf x))_j=a_j(x_j-c_j)$, $L=\max_j a_j$, and $\mu=\min_j a_j\gt0$ is a strong-convexity parameter. The compact feasible box guarantees a minimizer exists, and strict convexity makes it unique. Fixed coordinates with $\ell_j=u_j$ are allowed. Although this family is easy enough to solve analytically, it exposes the operations used in larger problems.

Supporting slopes and constraints. A vector $\mathbf g$ is a subgradient of $\Phi$ at $\mathbf x$ if $\Phi(\mathbf z)\geq\Phi(\mathbf x)+\mathbf g^{\mathsf T}(\mathbf z-\mathbf x)$ for all $\mathbf z$. The subdifferential of the absolute value is $\partial|x|=\{1\}$ for $x\gt0$, $\{-1\}$ for $x\lt0$, and $[-1,1]$ at $0$. The code chooses the valid value $0$ at the kink. This choice does not mean the entire composite objective has zero subgradient there.

Euclidean projection is $\Pi_C(\mathbf y)=\arg\min_{\mathbf z\in C}\tfrac12\lVert\mathbf z-\mathbf y\rVert_2^2$. On a box it is coordinatewise clipping:

$$[\Pi_C(\mathbf y)]_j=\min(u_j,\max(\ell_j,y_j)).$$

The indicator $I_C$ is zero on $C$ and $+\infty$ outside it. A proximal map solves a regularized subproblem,

$$\operatorname{prox}_{\eta h}(\mathbf y)=\arg\min_{\mathbf z}\left(h(\mathbf z)+\frac1{2\eta}\lVert\mathbf z-\mathbf y\rVert_2^2\right),\qquad\eta\gt0.$$

Projection is the special case $h=I_C$. For $h(\mathbf x)=\lambda\lVert\mathbf x\rVert_1+I_C(\mathbf x)$, separability reduces the prox to independent scalar minimizations. Define soft thresholding $S_\tau(v)=\operatorname{sign}(v)\max(|v|-\tau,0)$, with $S_\tau(0)=0$. Then

$$[\operatorname{prox}_{\eta h}(\mathbf y)]_j=\operatorname{clip}_{[\ell_j,u_j]}(S_{\eta\lambda}(y_j)).$$

This composition is exact for L1 plus a product of intervals; it is not a generic identity for arbitrary penalties and constraint sets. Completing the square in each original quadratic coordinate gives an independent solution check:

$$x_j^\star=\operatorname{clip}_{[\ell_j,u_j]}(S_{\lambda/a_j}(c_j)).$$

For $\lambda=0$, the formula reduces to clipping the center. A constrained minimum can have a nonzero ordinary gradient: the active boundary can prevent every feasible downhill move.

How Convex Optimization Works, Step by Step

The common inputs are the initial point, center, positive curvature vector, lower and upper bounds, base step $\eta$, penalty $\lambda$, update cap $K$, and a mode. All modes validate the same rate guard $0\lt\eta\leq1/L$; the coordinate mode validates but does not use this base rate.

1. Validate and project the initial point. Check finite vectors and controls, positive curvatures, matching dimensions, and $\ell_j\leq u_j$. Clip the input point into $C$ before recording step zero. A starting point outside the box is supported; the recorded starting loss belongs to the projected point.
2. Record the starting objective. Evaluate $\Phi$ once, with empty gradient and candidate fields. No derivative work has yet been performed.
3. Choose one update. The four branches below all read the old point. The one-based index is $t=1,2,\ldots,K$.
4. Commit a feasible point and record. Replace the vector only after the branch has constructed its candidate. Evaluate the full objective and store independent copies of all diagnostic vectors. Repeat until exactly $K$ updates have completed.

SUBGRAD: projected composite subgradient. Form $g_{t,j}=a_j(x_{t-1,j}-c_j)+\lambda\,\operatorname{sign}_0(x_{t-1,j})$, where the chosen sign at zero is $0$. Use a diminishing step $\alpha_t=\eta/t$ and

$$\mathbf x_t=\Pi_C(\mathbf x_{t-1}-\alpha_t\mathbf g_t).$$

The squared-summable, nonsummable harmonic schedule is intentional. Subgradient updates need not decrease $\Phi$ on every step, even with a convex objective.

PG: projected gradient for the smooth problem. The implementation requires $\lambda=0$ and rejects a nonzero penalty:

$$\mathbf x_t=\Pi_C(\mathbf x_{t-1}-\eta\nabla f(\mathbf x_{t-1})).$$

Replacing the gradient by a subgradient would instead produce the preceding method. Ignoring the L1 term does not solve the penalized problem.

CD: cyclic exact coordinate descent. Choose $j=(t-1)\bmod d$, leave all other entries unchanged, and solve that scalar subproblem:

$$x_{t,j}=\operatorname{clip}_{[\ell_j,u_j]}(S_{\lambda/a_j}(c_j)).$$

The code evaluates one old-point partial derivative for the trace, then uses this direct scalar minimizer. Its recorded forward field is $c_j$, its thresholded field is $S_{\lambda/a_j}(c_j)$, and its step-size field is $1/a_j$. The threshold is computed as $\lambda/a_j$ directly. One sweep means $d$ coordinate updates, not one update. Exact completion in one sweep is a special property of this separable diagonal objective.

PROX: proximal gradient for the composite problem. Make the smooth forward step $\mathbf y_t=\mathbf x_{t-1}-\eta\nabla f(\mathbf x_{t-1})$, then solve the combined penalty-and-box proximal subproblem:

$$x_{t,j}=\operatorname{clip}_{[\ell_j,u_j]}(S_{\eta\lambda}(y_{t,j})).$$

All smooth gradients are evaluated at the same old vector, even if a larger replacement objective couples coordinates. A history row numbered $t$ records the new point and its objective, but its gradient belongs to the update from $\mathbf x_{t-1}$. For CD, gradient and intermediate fields have one entry identified by coordinate; full-vector branches use $d$ entries. Coordinate $-1$ denotes a full-vector update or the initial row.

The solver returns its last point. It does not stop on zero gradients, no-ops, or a previously optimal state, and it does not replace the live state by a best-so-far point.

flowchart TD
    A["Validate vectors, box, positive curvature, rate, mode, and penalty"] --> B["Project initial point; evaluate starting objective; record"]
    B --> C{"Completed K updates?"}
    C -->|Yes| Z["Return final feasible point, counts, and history"]
    C -->|No| D{"Coordinate descent?"}
    D -->|Yes| E["Pick cyclic coordinate; record one partial derivative"]
    E --> F["Soft-threshold center by penalty over curvature; clip one coordinate"]
    D -->|No| G{"Subgradient mode?"}
    G -->|Yes| H["Composite subgradient; step = base rate divided by update index"]
    G -->|No| I["Smooth full gradient; fixed step; PG requires penalty zero"]
    H --> J["Compute whole forward candidate from old point"]
    I --> J
    J --> L{"Proximal mode?"}
    L -->|Yes| M["Soft-threshold by step times penalty"]
    L -->|No| N["Keep the forward candidate"]
    M --> P["Clip candidate to the box"]
    N --> P
    F --> R["Commit; evaluate full objective; copy state; count work"]
    P --> R
    R --> C

Worked Example with a Full Trace

Use initial point $\mathbf x_0=(0,1)$, center $\mathbf c=(3,-1)$, curvatures $(2,4)$, box $[-1,1]^2$, base step $\eta=1/4$, and $K=3$. The smooth part is $f(\mathbf x)=(x_1-3)^2+2(x_2+1)^2$, with $L=4$. The starting point is already feasible. There is no randomness.

For SUBGRAD, CD, and PROX, set $\lambda=1$. The starting smooth value is $9+8=17$, the L1 term is $1$, and $\Phi(\mathbf x_0)=18$. The exact solution formula gives $x_1^\star=\operatorname{clip}(S_{1/2}(3))=1$ and $x_2^\star=\operatorname{clip}(S_{1/4}(-1))=-3/4$, with $\Phi^\star=4+1/8+7/4=47/8=5.875$.

Projected subgradient, $\lambda=1$.

1. At $(0,1)$, smooth gradient $(-6,8)$ plus chosen L1 subgradient $(0,1)$ gives $\mathbf g_1=(-6,9)$. With $\alpha_1=1/4$, the forward point is $(3/2,-5/4)$; clipping gives $(1,-1)$. Its smooth value is $4$, its L1 term is $2$, and $\Phi=6$.
2. At $(1,-1)$, the smooth gradient is $(-4,0)$ and the L1 subgradient is $(1,-1)$, so $\mathbf g_2=(-3,-1)$. With $\alpha_2=1/8$, the forward point is $(11/8,-7/8)$; clipping gives $(1,-7/8)$. Smooth loss is $4+1/32$, penalty is $15/8$, and $\Phi=189/32=5.90625$.
3. At $(1,-7/8)$, smooth gradient $(-4,1/2)$ plus $(1,-1)$ gives $\mathbf g_3=(-3,-1/2)$. With $\alpha_3=1/12$, the forward point is $(5/4,-5/6)$; clipping gives $(1,-5/6)$. The smooth part is $4+1/18$, penalty is $11/6$, and $\Phi=53/9\approx5.8888888889$. Fractions describe the real-arithmetic recurrence; the program uses rounded binary64 division for the third step.

Projected gradient, $\lambda=0$. This run solves the smooth constrained problem, whose starting value is $17$ and optimum value is $4$.

1. Smooth gradient at $(0,1)$ is $(-6,8)$. The forward point is $(3/2,-1)$; clipping gives $(1,-1)$ with objective $4$.
2. The new smooth gradient is $(-4,0)$, which is nonzero. Its forward point is $(2,-1)$, and projection returns $(1,-1)$ again with objective $4$.
3. The same gradient, forward point, and projected point repeat. The point is globally optimal for the constrained smooth problem because the first coordinate cannot increase beyond its upper bound, not because the ordinary gradient vanishes.

Cyclic exact coordinate descent, $\lambda=1$.

1. Choose zero-based coordinate $0$. Its recorded old partial derivative is $-6$. Its curvature inverse is $1/2$, threshold is $1/2$, and soft-thresholded center is $3-1/2=5/2$. Clipping gives $x_1=1$, while the other coordinate stays $1$. At $(1,1)$, smooth loss is $4+8=12$, penalty is $2$, and $\Phi=14$.
2. Choose coordinate $1$. Its old partial derivative is $4(1+1)=8$. Curvature inverse and threshold are both $1/4$. Soft-thresholded center is $-1+1/4=-3/4$, already inside the box. The new point is $(1,-3/4)$ with $\Phi=47/8$. One complete two-coordinate sweep has solved this separable problem.
3. Choose coordinate $0$ again. Its old partial derivative is $-4$; the same center, threshold, and clipping keep it at $1$. The other coordinate stays $-3/4$, and the objective remains $47/8$.

Proximal gradient, $\lambda=1$.

1. Smooth gradient at $(0,1)$ is $(-6,8)$. The forward point is $(3/2,-1)$. Threshold $\eta\lambda=1/4$ gives $(5/4,-3/4)$, and box projection gives $(1,-3/4)$. Smooth loss is $4+1/8=4.125$, the penalty is $7/4=1.75$, and $\Phi=47/8$.
2. At $(1,-3/4)$ the smooth gradient is $(-4,1)$. The forward point is $(2,-1)$. Thresholding gives $(7/4,-3/4)$, and projection returns $(1,-3/4)$ with the same loss.
3. The same gradient, forward point, thresholded point, and accepted point repeat. This is a proximal fixed point with a nonzero smooth gradient.

The three full-vector methods each compute three full gradients, or six derivative components in dimension two. CD computes three partial derivatives. All four runs also evaluate four complete objectives, including the initial value. PG's loss is for a different objective; neither smaller printed loss nor equal update counts alone make a fair comparison.

Implementation in Python, C++, and Java

Each panel is a complete copy-paste-runnable program. The shared callable is optimize(initial, center, weights, lower, upper, rate, penalty, steps, mode), with modes SUBGRAD, PG, CD, and PROX. The printed sample is the entire three-update trace above. Parameters are projected before the initial objective is recorded.

The teaching API supports $1\leq d\leq32$, integer $0\leq K\leq2000$, finite initial coordinates, centers, and bounds of magnitude at most $10^6$, positive curvatures $10^{-6}\leq a_j\leq10^6$, $0\lt\eta\leq10^6$ with $\eta\leq1/L$, and $0\leq\lambda\leq10^6$. Empty boxes and mismatched shapes are rejected. PG additionally requires $\lambda=0$. CD checks the base rate for a uniform input contract but uses its own exact scalar subproblem rather than that rate. Python rejects booleans as numeric controls and requires an actual integer update count; C++ and Java encode the count type in their signatures.

The code uses scalar branches for clipping and soft thresholding, with a zero result exactly at the threshold. It keeps the multiplication, addition, and subtraction order the same in all languages. For full-vector methods, every derivative is computed from the same old vector before any candidate is committed. CD records exactly one partial derivative for diagnostics, even though the direct minimizer formula does not need that diagnostic.

The result contains the last feasible point, its objective, counts of full gradients, counts of coordinate partials, derivative-component count, objective calls, and complete state history. Initial candidate fields are empty. SUBGRAD and PG leave thresholded empty; CD's diagnostic vectors contain one scalar, and its coordinate field identifies which parameter was updated. A full-vector update uses coordinate $-1$.

Run Python 3 with python main.py. Compile C++17 with g++ -std=c++17 -ffp-contract=off main.cpp -o convex and run ./convex; avoid fast-math. Save Java as Main.java, use Java 17 or later, compile with javac Main.java, and run java Main. Java's decimal helper only formats stored binary64 values to eight decimal places with half-even rounding.

These are reference implementations for the stated diagonal family, not generic constrained solvers. There is no objective callback, line search, random coordinate selection, active-set shortcut, or convergence-tolerance exit. A tiny subgradient step can underflow or round to a no-op; it is still counted. Replacing the objective or feasible set requires new gradients, scalar subproblems or proximal maps, and appropriate smoothness and numerical assumptions.

from dataclasses import dataclass
import math


@dataclass
class State:
    iteration: int
    point: list
    coordinate: int
    step_size: float
    gradient: list
    forward: list
    thresholded: list
    value: float


@dataclass
class Result:
    point: list
    value: float
    steps: int
    full_gradients: int
    coordinate_gradients: int
    gradient_components: int
    value_evaluations: int
    history: list


def clip(x, lo, hi):
    return lo if x < lo else hi if x > hi else x


def soft(x, threshold):
    if x > threshold:
        return x - threshold
    if x < -threshold:
        return x + threshold
    return 0.0


def objective(point, center, weights, penalty):
    total, absolute_sum = 0.0, 0.0
    for x, c, a in zip(point, center, weights):
        delta = x - c
        square = delta * delta
        term = a * square
        total = total + term
        absolute_sum = absolute_sum + abs(x)
    smooth = 0.5 * total
    regularizer = penalty * absolute_sum
    return smooth + regularizer


def optimize(initial, center, weights, lower, upper, rate, penalty, steps, mode):
    if type(steps) is not int or not 0 <= steps <= 2000:
        raise ValueError("steps must be an integer in [0,2000]")
    if mode not in ("SUBGRAD", "PG", "CD", "PROX"):
        raise ValueError("invalid method")

    def number(x):
        return type(x) is int or (type(x) is float and math.isfinite(x))

    d = len(initial)
    if not 1 <= d <= 32 or any(len(v) != d for v in [center, weights, lower, upper]):
        raise ValueError("invalid vector sizes")
    if any(not number(x) or abs(x) > 1000000
           for row in [initial, center, lower, upper] for x in row):
        raise ValueError("invalid coordinates or bounds")
    if any(not number(a) or not 1e-6 <= a <= 1000000 for a in weights):
        raise ValueError("invalid positive curvature")
    if any(lo > hi for lo, hi in zip(lower, upper)):
        raise ValueError("empty box")
    if not number(rate) or not 0 < rate <= 1000000:
        raise ValueError("invalid learning rate")
    if not number(penalty) or not 0 <= penalty <= 1000000:
        raise ValueError("invalid L1 penalty")
    if rate > 1.0 / float(max(weights)):
        raise ValueError("learning rate exceeds 1/L")
    if mode == "PG" and penalty != 0:
        raise ValueError("PG requires a smooth objective: penalty must be zero")
    point, c, a, lo, hi = [list(map(float, v)) for v in [initial, center, weights, lower, upper]]
    rate, penalty = float(rate), float(penalty)
    point = [clip(x, l, u) for x, l, u in zip(point, lo, hi)]
    value = objective(point, c, a, penalty)
    history = [State(0, point[:], -1, 0.0, [], [], [], value)]
    full, partial = 0, 0
    for k in range(1, steps + 1):
        candidate = point[:]
        coordinate = -1
        gradient, forward, thresholded = [], [], []
        if mode == "CD":
            coordinate = j = (k - 1) % d
            delta = point[j] - c[j]
            gradient = [a[j] * delta]
            step_size = 1.0 / a[j]
            threshold = penalty / a[j]
            forward = [c[j]]
            thresholded = [soft(c[j], threshold)]
            candidate[j] = clip(thresholded[0], lo[j], hi[j])
            partial += 1
        else:
            step_size = rate / float(k) if mode == "SUBGRAD" else rate
            threshold = step_size * penalty
            for x, cj, aj in zip(point, c, a):
                delta = x - cj
                g = aj * delta
                if mode == "SUBGRAD":
                    sign = 1.0 if x > 0 else -1.0 if x < 0 else 0.0
                    regularizer = penalty * sign
                    g = g + regularizer
                gradient.append(g)
                displacement = step_size * g
                forward.append(x - displacement)
            for j in range(d):
                y = forward[j]
                if mode == "PROX":
                    y = soft(y, threshold)
                    thresholded.append(y)
                candidate[j] = clip(y, lo[j], hi[j])
            full += 1
        point = candidate
        value = objective(point, c, a, penalty)
        history.append(State(k, point[:], coordinate, step_size, gradient[:],
                             forward[:], thresholded[:], value))
    return Result(point[:], value, steps, full, partial, d * full + partial,
                  steps + 1, history)


def vector(xs):
    return "(" + ",".join(f"{x:.8f}" for x in xs) + ")"


def main():
    for mode in ("SUBGRAD", "PG", "CD", "PROX"):
        penalty = 0.0 if mode == "PG" else 1.0
        r = optimize([0, 1], [3, -1], [2, 4], [-1, -1], [1, 1],
                     0.25, penalty, 3, mode)
        print(f"{mode} lambda={penalty:.8f}")
        for s in r.history:
            print(f"k={s.iteration} coord={s.coordinate} x={vector(s.point)} F={s.value:.8f}")
        print(f"steps={r.steps} full={r.full_gradients} partial={r.coordinate_gradients} "
              f"components={r.gradient_components} values={r.value_evaluations}")


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
using namespace std;

struct State {
    int iteration;
    vector<double> point;
    int coordinate;
    double step_size;
    vector<double> gradient, forward, thresholded;
    double value;
};
struct Result {
    vector<double> point;
    double value;
    int steps, full_gradients, coordinate_gradients, gradient_components, value_evaluations;
    vector<State> history;
};

double clip(double x, double lo, double hi) {
    return x < lo ? lo : x > hi ? hi : x;
}
double soft(double x, double threshold) {
    if (x > threshold) return x - threshold;
    if (x < -threshold) return x + threshold;
    return 0.0;
}
double objective(const vector<double>& point, const vector<double>& center,
                 const vector<double>& weights, double penalty) {
    double total = 0.0, absolute_sum = 0.0;
    for (size_t j = 0; j < point.size(); ++j) {
        double delta = point[j] - center[j];
        double square = delta * delta;
        double term = weights[j] * square;
        total = total + term;
        absolute_sum = absolute_sum + abs(point[j]);
    }
    double smooth = 0.5 * total;
    double regularizer = penalty * absolute_sum;
    return smooth + regularizer;
}

Result optimize(const vector<double>& initial, const vector<double>& center,
                const vector<double>& weights, const vector<double>& lower,
                const vector<double>& upper, double rate, double penalty,
                int steps, const string& mode) {
    if (steps < 0 || steps > 2000 ||
        (mode != "SUBGRAD" && mode != "PG" && mode != "CD" && mode != "PROX"))
        throw invalid_argument("invalid steps or method");
    size_t d = initial.size();
    if (d < 1 || d > 32 || center.size() != d || weights.size() != d ||
        lower.size() != d || upper.size() != d)
        throw invalid_argument("invalid vector sizes");
    for (size_t j = 0; j < d; ++j) {
        if (!isfinite(initial[j]) || abs(initial[j]) > 1000000 ||
            !isfinite(center[j]) || abs(center[j]) > 1000000 ||
            !isfinite(lower[j]) || abs(lower[j]) > 1000000 ||
            !isfinite(upper[j]) || abs(upper[j]) > 1000000 || lower[j] > upper[j] ||
            !isfinite(weights[j]) || weights[j] < 1e-6 || weights[j] > 1000000)
            throw invalid_argument("invalid coordinates, box, or curvature");
    }
    if (!isfinite(rate) || rate <= 0 || rate > 1000000 ||
        !isfinite(penalty) || penalty < 0 || penalty > 1000000)
        throw invalid_argument("invalid rate or penalty");
    if (rate > 1.0 / *max_element(weights.begin(), weights.end()))
        throw invalid_argument("learning rate exceeds 1/L");
    if (mode == "PG" && penalty != 0)
        throw invalid_argument("PG requires penalty zero");
    vector<double> point = initial;
    for (size_t j = 0; j < d; ++j) point[j] = clip(point[j], lower[j], upper[j]);
    double value = objective(point, center, weights, penalty);
    vector<State> history;
    history.push_back(State{0, point, -1, 0.0, {}, {}, {}, value});
    int full = 0, partial = 0;
    for (int k = 1; k <= steps; ++k) {
        vector<double> candidate = point, gradient, forward, thresholded;
        int coordinate = -1;
        double step_size;
        if (mode == "CD") {
            coordinate = (k - 1) % int(d);
            size_t j = size_t(coordinate);
            double delta = point[j] - center[j];
            gradient.push_back(weights[j] * delta);
            step_size = 1.0 / weights[j];
            double threshold = penalty / weights[j];
            forward.push_back(center[j]);
            thresholded.push_back(soft(center[j], threshold));
            candidate[j] = clip(thresholded[0], lower[j], upper[j]);
            ++partial;
        } else {
            step_size = mode == "SUBGRAD" ? rate / double(k) : rate;
            double threshold = step_size * penalty;
            for (size_t j = 0; j < d; ++j) {
                double delta = point[j] - center[j];
                double g = weights[j] * delta;
                if (mode == "SUBGRAD") {
                    double sign = point[j] > 0 ? 1.0 : point[j] < 0 ? -1.0 : 0.0;
                    double regularizer = penalty * sign;
                    g = g + regularizer;
                }
                gradient.push_back(g);
                double displacement = step_size * g;
                forward.push_back(point[j] - displacement);
            }
            for (size_t j = 0; j < d; ++j) {
                double y = forward[j];
                if (mode == "PROX") {
                    y = soft(y, threshold);
                    thresholded.push_back(y);
                }
                candidate[j] = clip(y, lower[j], upper[j]);
            }
            ++full;
        }
        point = candidate;
        value = objective(point, center, weights, penalty);
        history.push_back(State{k, point, coordinate, step_size, gradient, forward, thresholded, value});
    }
    return Result{point, value, steps, full, partial, int(d) * full + partial, steps + 1, history};
}

string vector_text(const vector<double>& xs) {
    ostringstream out;
    out << fixed << setprecision(8) << '(';
    for (size_t j = 0; j < xs.size(); ++j) {
        if (j) out << ',';
        out << xs[j];
    }
    out << ')';
    return out.str();
}
int main() {
    for (const string& mode : {string("SUBGRAD"), string("PG"), string("CD"), string("PROX")}) {
        double penalty = mode == "PG" ? 0.0 : 1.0;
        Result r = optimize({0, 1}, {3, -1}, {2, 4}, {-1, -1}, {1, 1}, 0.25, penalty, 3, mode);
        cout << fixed << setprecision(8) << mode << " lambda=" << penalty << '\n';
        for (const auto& s : r.history)
            cout << "k=" << s.iteration << " coord=" << s.coordinate << " x=" << vector_text(s.point)
                 << " F=" << s.value << '\n';
        cout << "steps=" << r.steps << " full=" << r.full_gradients
             << " partial=" << r.coordinate_gradients << " components=" << r.gradient_components
             << " values=" << r.value_evaluations << '\n';
    }
    return 0;
}

import java.math.BigDecimal;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.List;

public class Main {
    static class State {
        int iteration;
        double[] point;
        int coordinate;
        double step_size;
        double[] gradient, forward, thresholded;
        double value;
        State(int k, double[] x, int j, double rate, double[] g, double[] y,
              double[] shrunk, double f) {
            iteration = k; point = x.clone(); coordinate = j; step_size = rate;
            gradient = g.clone(); forward = y.clone(); thresholded = shrunk.clone(); value = f;
        }
    }
    static class Result {
        double[] point;
        double value;
        int steps, full_gradients, coordinate_gradients, gradient_components, value_evaluations;
        List<State> history;
        Result(double[] x, double f, int k, int full, int partial, int d, List<State> rows) {
            point = x.clone(); value = f; steps = k; full_gradients = full;
            coordinate_gradients = partial; gradient_components = d * full + partial;
            value_evaluations = k + 1; history = rows;
        }
    }
    static double clip(double x, double lo, double hi) {
        return x < lo ? lo : x > hi ? hi : x;
    }
    static double soft(double x, double threshold) {
        if (x > threshold) return x - threshold;
        if (x < -threshold) return x + threshold;
        return 0.0;
    }
    static double objective(double[] point, double[] center, double[] weights, double penalty) {
        double total = 0.0, absoluteSum = 0.0;
        for (int j = 0; j < point.length; ++j) {
            double delta = point[j] - center[j];
            double square = delta * delta;
            double term = weights[j] * square;
            total = total + term;
            absoluteSum = absoluteSum + Math.abs(point[j]);
        }
        double smooth = 0.5 * total;
        double regularizer = penalty * absoluteSum;
        return smooth + regularizer;
    }
    static Result optimize(double[] initial, double[] center, double[] weights,
                           double[] lower, double[] upper, double rate, double penalty,
                           int steps, String mode) {
        if (steps < 0 || steps > 2000 ||
            (!"SUBGRAD".equals(mode) && !"PG".equals(mode) &&
             !"CD".equals(mode) && !"PROX".equals(mode)))
            throw new IllegalArgumentException("invalid steps or method");
        int d = initial.length;
        if (d < 1 || d > 32 || center.length != d || weights.length != d ||
            lower.length != d || upper.length != d)
            throw new IllegalArgumentException("invalid vector sizes");
        double largest = 0.0;
        for (int j = 0; j < d; ++j) {
            if (!Double.isFinite(initial[j]) || Math.abs(initial[j]) > 1000000 ||
                !Double.isFinite(center[j]) || Math.abs(center[j]) > 1000000 ||
                !Double.isFinite(lower[j]) || Math.abs(lower[j]) > 1000000 ||
                !Double.isFinite(upper[j]) || Math.abs(upper[j]) > 1000000 || lower[j] > upper[j] ||
                !Double.isFinite(weights[j]) || weights[j] < 1e-6 || weights[j] > 1000000)
                throw new IllegalArgumentException("invalid coordinates, box, or curvature");
            largest = Math.max(largest, weights[j]);
        }
        if (!Double.isFinite(rate) || rate <= 0 || rate > 1000000 ||
            !Double.isFinite(penalty) || penalty < 0 || penalty > 1000000)
            throw new IllegalArgumentException("invalid rate or penalty");
        if (rate > 1.0 / largest) throw new IllegalArgumentException("learning rate exceeds 1/L");
        if (mode.equals("PG") && penalty != 0)
            throw new IllegalArgumentException("PG requires penalty zero");
        double[] point = initial.clone();
        for (int j = 0; j < d; ++j) point[j] = clip(point[j], lower[j], upper[j]);
        double value = objective(point, center, weights, penalty);
        List<State> history = new ArrayList<>();
        history.add(new State(0, point, -1, 0.0, new double[0], new double[0], new double[0], value));
        int full = 0, partial = 0;
        for (int k = 1; k <= steps; ++k) {
            double[] candidate = point.clone(), gradient, forward, thresholded;
            int coordinate = -1;
            double stepSize;
            if (mode.equals("CD")) {
                int j = (k - 1) % d;
                coordinate = j;
                double delta = point[j] - center[j];
                gradient = new double[]{weights[j] * delta};
                stepSize = 1.0 / weights[j];
                double threshold = penalty / weights[j];
                forward = new double[]{center[j]};
                thresholded = new double[]{soft(center[j], threshold)};
                candidate[j] = clip(thresholded[0], lower[j], upper[j]);
                ++partial;
            } else {
                stepSize = mode.equals("SUBGRAD") ? rate / (double) k : rate;
                double threshold = stepSize * penalty;
                gradient = new double[d]; forward = new double[d];
                thresholded = mode.equals("PROX") ? new double[d] : new double[0];
                for (int j = 0; j < d; ++j) {
                    double delta = point[j] - center[j];
                    double g = weights[j] * delta;
                    if (mode.equals("SUBGRAD")) {
                        double sign = point[j] > 0 ? 1.0 : point[j] < 0 ? -1.0 : 0.0;
                        double regularizer = penalty * sign;
                        g = g + regularizer;
                    }
                    gradient[j] = g;
                    double displacement = stepSize * g;
                    forward[j] = point[j] - displacement;
                }
                for (int j = 0; j < d; ++j) {
                    double y = forward[j];
                    if (mode.equals("PROX")) {
                        y = soft(y, threshold);
                        thresholded[j] = y;
                    }
                    candidate[j] = clip(y, lower[j], upper[j]);
                }
                ++full;
            }
            point = candidate;
            value = objective(point, center, weights, penalty);
            history.add(new State(k, point, coordinate, stepSize, gradient, forward, thresholded, value));
        }
        return new Result(point, value, steps, full, partial, d, history);
    }
    static String fixed8(double x) {
        String s = new BigDecimal(x).setScale(8, RoundingMode.HALF_EVEN).toPlainString();
        return Double.doubleToRawLongBits(x) < 0 && s.charAt(0) != '-' ? "-" + s : s;
    }
    static String vector(double[] xs) {
        StringBuilder s = new StringBuilder("(");
        for (int j = 0; j < xs.length; ++j) {
            if (j > 0) s.append(',');
            s.append(fixed8(xs[j]));
        }
        return s.append(')').toString();
    }
    public static void main(String[] args) {
        for (String mode : new String[]{"SUBGRAD", "PG", "CD", "PROX"}) {
            double penalty = mode.equals("PG") ? 0.0 : 1.0;
            Result r = optimize(new double[]{0, 1}, new double[]{3, -1}, new double[]{2, 4},
                new double[]{-1, -1}, new double[]{1, 1}, 0.25, penalty, 3, mode);
            System.out.println(mode + " lambda=" + fixed8(penalty));
            for (State s : r.history)
                System.out.println("k=" + s.iteration + " coord=" + s.coordinate + " x=" +
                    vector(s.point) + " F=" + fixed8(s.value));
            System.out.println("steps=" + r.steps + " full=" + r.full_gradients +
                " partial=" + r.coordinate_gradients + " components=" + r.gradient_components +
                " values=" + r.value_evaluations);
        }
    }
}

Expected standard output, identical in all three languages:

• SUBGRAD lambda=1.00000000
• k=0 coord=-1 x=(0.00000000,1.00000000) F=18.00000000
• k=1 coord=-1 x=(1.00000000,-1.00000000) F=6.00000000
• k=2 coord=-1 x=(1.00000000,-0.87500000) F=5.90625000
• k=3 coord=-1 x=(1.00000000,-0.83333333) F=5.88888889
• steps=3 full=3 partial=0 components=6 values=4
• PG lambda=0.00000000
• k=0 coord=-1 x=(0.00000000,1.00000000) F=17.00000000
• k=1 coord=-1 x=(1.00000000,-1.00000000) F=4.00000000
• k=2 coord=-1 x=(1.00000000,-1.00000000) F=4.00000000
• k=3 coord=-1 x=(1.00000000,-1.00000000) F=4.00000000
• steps=3 full=3 partial=0 components=6 values=4
• CD lambda=1.00000000
• k=0 coord=-1 x=(0.00000000,1.00000000) F=18.00000000
• k=1 coord=0 x=(1.00000000,1.00000000) F=14.00000000
• k=2 coord=1 x=(1.00000000,-0.75000000) F=5.87500000
• k=3 coord=0 x=(1.00000000,-0.75000000) F=5.87500000
• steps=3 full=0 partial=3 components=3 values=4
• PROX lambda=1.00000000
• k=0 coord=-1 x=(0.00000000,1.00000000) F=18.00000000
• k=1 coord=-1 x=(1.00000000,-0.75000000) F=5.87500000
• k=2 coord=-1 x=(1.00000000,-0.75000000) F=5.87500000
• k=3 coord=-1 x=(1.00000000,-0.75000000) F=5.87500000
• steps=3 full=3 partial=0 components=6 values=4

Complexity, Correctness, and Invariants

Let $d$ be the dimension and $K$ the requested number of updates. Validation, initial projection, and the starting objective cost $O(d)$. A full gradient, thresholding, projection, and evaluation of the reference objective each cost $O(d)$. The bare exact coordinate update is $O(1)$ for this separable family, but these programs also copy the whole parameter vector, evaluate the whole objective, and store a whole snapshot on every coordinate update.

Consequently best, average, and worst running time for every mode in these trace-producing programs are $\Theta((K+1)d)$ for fixed valid $K,d$. There is no early termination. Working vectors require $O(d)$ space, while retained history uses $\Theta((K+1)d)$ total space. A streamlined CD implementation without full monitoring or snapshots can reduce its core update to $O(1)$; it must also change the reported observation policy. The full-gradient modes use $Kd$ derivative components, whereas CD records $K$. Every run uses $K+1$ full objective calls.

For a general objective, separate the costs of a full gradient $C_g$, a chosen partial derivative or scalar subproblem $C_j$, full objective monitoring $C_f$, and a projection or proximal map $C_{\rm map}$. General projection is not necessarily a cheap clip. A unit labeled “iteration” can mean a coordinate update, a sweep, or a full-vector step; state which one before comparing complexities.

Feasibility and recurrence invariant. At every recorded state, $\mathbf x_t\in C$. Initialization establishes this by projection. Each full-vector branch computes its prescribed old-point gradient or subgradient, constructs a whole candidate, and clips every coordinate. CD changes only its selected coordinate to a clipped scalar solution. Thus feasibility and the chosen recurrence are preserved by induction. The value stored in each row is evaluated at that same committed point. The diagnostic derivative in that row belongs to the incoming update. Copying snapshots prevents later updates from rewriting earlier states. The invariant describes the rounded code recurrence in binary64, not exact real-arithmetic values for arbitrary inputs.

Optimality is not a zero ordinary gradient test. For this composite box problem, the convex optimality condition is

$$\mathbf0\in\nabla f(\mathbf x^\star)+\lambda\partial\lVert\mathbf x^\star\rVert_1+N_C(\mathbf x^\star),$$

where $N_C(\mathbf x)=\{\mathbf w:\mathbf w^{\mathsf T}(\mathbf z-\mathbf x)\leq0\text{ for all }\mathbf z\in C\}$ is the normal cone. For example, at $(1,-3/4)$, the smooth gradient is $(-4,1)$, the L1 subgradient is $(1,-1)$, and the remaining vector $(-3,0)$ is balanced by an outward normal $(3,0)$ at the active upper bound.

An equivalent residual is the proximal gradient mapping

$$G_\eta(\mathbf x)=\frac{\mathbf x-\operatorname{prox}_{\eta h}(\mathbf x-\eta\nabla f(\mathbf x))}{\eta},\qquad h=\lambda\lVert\cdot\rVert_1+I_C.$$

In exact arithmetic $G_\eta(\mathbf x)=\mathbf0$ characterizes a minimizer of this convex problem. For PG, set $\lambda=0$. A floating-point zero or a small residual requires an appropriate numerical tolerance and scaling; the programs deliberately make no convergence claim from a rounded no-op.

Projected/proximal descent. Proximal gradient minimizes the quadratic surrogate

$$f(\mathbf x)+\nabla f(\mathbf x)^{\mathsf T}(\mathbf z-\mathbf x)+\frac1{2\eta}\lVert\mathbf z-\mathbf x\rVert_2^2+\lambda\lVert\mathbf z\rVert_1+I_C(\mathbf z).$$

For an $L$-Lipschitz smooth gradient and $0\lt\eta\leq1/L$, this surrogate majorizes the true composite objective and equals it at the old feasible point. Its exact minimizer therefore cannot increase the objective. PG is the same argument with the L1 coefficient zero. Standard convex proximal-gradient bounds give an $O(1/K)$ objective gap with the appropriate fixed step and exact operations. The implementation's diagonal formula evaluates the required prox exactly as a mathematical map, subject to floating-point rounding. See Parikh and Boyd's proximal algorithms for the general framework.

Subgradient progress is different. With a valid subgradient $\mathbf g_t$, nonexpansive projection and the supporting inequality give

$$\lVert\mathbf x_t-\mathbf x^\star\rVert_2^2\leq\lVert\mathbf x_{t-1}-\mathbf x^\star\rVert_2^2-2\alpha_t(\Phi(\mathbf x_{t-1})-\Phi^\star)+\alpha_t^2\lVert\mathbf g_t\rVert_2^2.$$

If all chosen subgradients have norm at most $G$ and the initial distance is at most $R$, summing yields

$$\min_{0\leq t\lt K}(\Phi(\mathbf x_t)-\Phi^\star)\leq\frac{R^2+G^2\sum_{t=1}^K\alpha_t^2}{2\sum_{t=1}^K\alpha_t}.$$

Our box bounds the subgradients, and $\alpha_t=\eta/t$ has a divergent sum and a finite sum of squares. This proves convergence of the best objective gap in exact arithmetic as the horizon grows; it does not prove that every update is a descent step or that the finite last state is optimal. This general bound for the harmonic schedule is only of order $1/\log K$. The familiar $O(1/\sqrt K)$ general convex subgradient bound uses a suitable different schedule, such as a horizon-dependent constant step. The code returns its last state; a best-so-far or averaged reporting policy would need to be added explicitly. Boyd and Park's subgradient notes explain these distinctions.

Coordinate correctness and limits. For the diagonal family, completing the square in a single coordinate gives the clipped soft-threshold formula. Updating one coordinate cannot increase the exact objective. Since no coordinate depends on another, visiting all $d$ coordinates once gives the unique global solution. For coupled objectives, this one-sweep conclusion disappears, and convergence depends on coordinate selection, smoothness, convexity, and separable penalty/constraint assumptions. Wright's coordinate-descent survey develops that broader setting.

Variants, Pitfalls, and Optimizations

Match the method to the objective. Ordinary projected gradient is for a smooth objective. Adding an L1 penalty creates a kink and requires a suitable subgradient or proximal treatment. Our PG branch refuses a nonzero penalty rather than silently ignoring it. At a kink, choosing one valid subgradient does not make it a descent direction. Convexity guarantees a global interpretation of valid optimality conditions, not monotone behavior of every first-order rule.

Do not interchange projection and shrinkage. For L1 plus a box, solve the scalar proximal subproblem to justify clipping after soft thresholding. Reversing these operations gives a different map. On the first coordinate of the worked step, forward value $1.5$ and threshold $0.25$ give the correct result $\operatorname{clip}(1.25)=1$; clipping first then shrinking gives $0.75$. Both are feasible, but the latter does not minimize the required proximal subproblem. With a box that excludes zero, shrinking after projection can even move outside the box. Proximal maps of arbitrary sums generally cannot be obtained by composing their separate maps.

A boundary optimum need not have zero gradient. In the smooth PG sample, $(1,-1)$ has gradient $(-4,0)$; its projected gradient mapping is zero. In the penalized sample, the solution balances smooth gradient, L1 subgradient, and the box normal. Use a residual adapted to the actual constrained composite problem. One small displacement can also be a rounding artifact or an overly small rate, and a fixed number of coordinate updates smaller than a sweep can leave untouched coordinates far from their optima.

Coordinate descent requires the right scalar problem. The direct center-threshold formula in the code relies on a diagonal quadratic. For ordinary lasso with $f(\mathbf x)=\tfrac12\lVert A\mathbf x-\mathbf b\rVert_2^2$, a nonzero column $A_j$ has coordinate curvature $\lVert A_j\rVert_2^2$ and a scalar update based on

$$S_{\lambda/\lVert A_j\rVert_2^2}\left(x_j-\frac{A_j^{\mathsf T}(A\mathbf x-\mathbf b)}{\lVert A_j\rVert_2^2}\right),$$

followed by interval clipping if the feasible set is a box. Maintain the residual to avoid recomputing a full matrix product for every coordinate. Zero columns need their own scalar treatment. A full-vector method instead needs a global gradient smoothness bound such as the largest eigenvalue of $A^{\mathsf T}A$, not merely its largest diagonal entry.

Arbitrary nonsmooth nonseparable objectives can have coordinatewise traps: on $[-1,1]^2$, $\max(x_1,x_2)$ has no strictly improving one-coordinate move from $(0,0)$ while moving both coordinates to $(-1,-1)$ decreases the objective. Thus coordinatewise lack of progress alone does not certify a global minimizer outside suitable structural assumptions.

General constraints change the map cost. A box clips in $O(d)$. A simplex projection usually needs more work, and a generic convex projection may itself require optimization. For intersections, one pass of sequential clipping or projection is not generally the Euclidean projection onto the intersection. Use an appropriate projection, proximal splitting, or a solver designed for those constraints.

Schedules, scaling, and numerical state. The harmonic schedule here is a convergence-oriented teaching choice, not a universally fastest subgradient schedule. Backtracking, acceleration, randomized or block coordinates, warm starts along a regularization path, and active-set strategies can improve appropriate problems while changing the state and work accounting. Keep the objective and selection policy fixed when comparing them.

Changing the scaling of data changes curvature, step-size limits, and the relative meaning of $\lambda$. Preserve a coordinate index or iteration count when resuming cyclic CD or a diminishing-step run. Record feasibility error, a suitable residual, objective change, and actual work in larger applications. The reference programs use bounded inputs, a finite box, and fixed caps; they deliberately reject zero curvature instead of hiding degenerate-coordinate rules.

Practice Problems and Real-World Use

Begin with the common two-dimensional trace, then change one structural feature at a time: remove the penalty, activate a bound, create a kink optimum, or couple coordinates. Compare the same objective when evaluating solution quality, and compare derivative components, objective calls, and elapsed work when evaluating cost.

For broader foundations, Boyd and Vandenberghe's Convex Optimization develops convexity, duality, and optimality conditions. The method-specific primary references above explain why different update rules require different convergence arguments.