Constraints and projected/penalized intuition

Compare the proposal, projection and penalty

Parameter space: θ₁ and θ₂

Loss slice: θ₂ = 0

Current parameter vector: (0, 1)
Current objective value: 4
Actual gradient components: (−2, 4)
Learning rate: 1
Unconstrained update vector: (2, −4)
Next parameter vector: (1, −1)
New objective value: 2.500000
Unconstrained proposed parameters: (2, −3)
Projection correction: (−1, 2)
Actual projected displacement: (1, −2)
Feasible region: box [-1, 1]²
Requested start before feasibility projection: (0, 1)
Evaluated feasible start: (0, 1)
Configured L2 penalty strength: 0

≈ binary64 calculation · = analytic classification. Replay selects the same evaluated step in every view.

Iteration 0: at θ ≈ (0, 1), L ≈ 4 and gradient ≈ (−2, 4). Rate η ≈ 1 gives update ≈ (2, −4). Before projection the proposal is ≈ (2, −3). Projection gives θ next ≈ (1, −1) and L next ≈ 2.5. The box projection corrects the proposal by ≈ (−1, 2); actual displacement is ≈ (1, −2). Feasibility and a soft objective penalty are different mechanisms. Stop: iteration budget. Parameters outside the displayed window are retained in the table; the view never changes the path. Feasible region: box [-1, 1]². λ = 0 adds λ‖θ‖² to the objective; a finite penalty does not enforce the region.

L=\tfrac a2(\theta_1-c)^2+2\theta_2^2

z_k=\theta_k-\eta\nabla J(\theta_k),\quad\theta_{k+1}=P_C(z_k)

z_{0}\approx \begin{bmatrix}0\\1\end{bmatrix}+\begin{bmatrix}2\\-4\end{bmatrix},\quad\theta_{1}=P_C(z_{0})\approx \begin{bmatrix}1\\-1\end{bmatrix}

J(\theta)=L(\theta)+\lambda\|\theta\|^2,\quad\nabla J=\nabla L+2\lambda\theta

Evaluated binary64 history. Replay changes the selected row, never the stored mathematics.
k	θ(k)	L(k)	gradient	η	proposed update	θ(k+1)	L(k+1)	unconstrained proposal	projection correction	base loss
0	(0, 1)	4	(−2, 4)	1	(2, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	4
1	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5
2	(1, 1)	2.5	(−1, 4)	1	(1, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	2.5
3	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5
4	(1, 1)	2.5	(−1, 4)	1	(1, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	2.5
5	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5
6	(1, 1)	2.5	(−1, 4)	1	(1, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	2.5
7	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5
8	(1, 1)	2.5	(−1, 4)	1	(1, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	2.5
9	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5
10	(1, 1)	2.5	(−1, 4)	1	(1, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	2.5
11	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5

Evaluated iteration 0

z_{0}\approx \begin{bmatrix}0\\1\end{bmatrix}+\begin{bmatrix}2\\-4\end{bmatrix},\quad\theta_{1}=P_C(z_{0})\approx \begin{bmatrix}1\\-1\end{bmatrix}

SpeedStep 1 of 12

Parameter θ₁

from −3 to 3, up to 2 decimal places

Parameter θ₂

from −3 to 3, up to 2 decimal places

Quadratic curvature a

from 0.25 to 4, up to 2 decimal places

Quadratic center c

from −2 to 2, up to 2 decimal places

Learning rate η

from 0.01 to 2, up to 2 decimal places

Iteration budget

from 1 to 40, whole numbers only

L2 penalty strength λ

from 0 to 4, up to 2 decimal places

Inspect the requested and feasible starts, then distinguish each unconstrained proposal, projected point, projection correction, and actual displacement. Switch between the fixed box and unit disk.

At the canonical start $(0,1)$ , the gradient is $(-2,4)$ and the rate is 1. The unconstrained proposal is $(2,-3)$ . Predict where the box rule sends it before advancing: the optimizer reaches $(1,-1)$ , not the proposal. Now switch to the disk and predict which direction radial projection preserves.

Intuition: a boundary and a preference are different rules

A hard constraint declares some parameter points inadmissible. Projection repairs an update by choosing the nearest admissible point to its proposal. The proposal remains part of the explanation even when it is never used as the next optimizer state.

A penalty instead changes how expensive parameter points are. An L2 penalty makes large distance from the origin more costly, balancing that preference against the base objective. It does not reject a point merely because it lies outside a box or disk. A finite penalized minimum can remain outside either region.

Feasibility and progress are also different. Repeated projection can keep every point allowed while the trajectory oscillates. The canonical box example deliberately shows this: a boundary is a geometric rule, not an automatic convergence mechanism.

Definition: the feasible set and the changed objective

Write $\theta=(x,y)$ and use the bounded quadratic family

L(\theta)=\frac a2(x-c)^2+2y^2,\qquad \nabla L=(a(x-c),4y),\qquad H_L=\begin{bmatrix}a&0\\0&4\end{bmatrix}.

The explorer has $1/4\leq a\leq4$ and $-2\leq c\leq2$ . Its hard-region choices are no constraint, the fixed box $C=[-1,1]^2$ , and the unit disk $C=\{\theta:\norm{\theta}_2\leq1\}$ . For a closed convex region, Euclidean projection is

P_C(z)=\argmin_{u\in C}\frac12\norm{u-z}_2^2.

The separate penalty mechanism uses $0\leq\lambda\leq4$ :

J(\theta)=L(\theta)+\lambda(x^2+y^2),\qquad \nabla J=\nabla L+2\lambda\theta,\qquad H_J=H_L+2\lambda I.

At $\lambda=0$ , J equals L. When a penalty is active, its gradient supplies the proposal and both the base loss and penalized objective remain meaningful quantities to inspect. Choosing a hard region independently applies its projection to that proposal.

Mechanics: record initialization and both parts of a step

For a requested initial point $\theta_{\mathrm{requested}}$ , the feasible initial state is explicitly

\theta_0=P_C(\theta_{\mathrm{requested}}),\qquad r_0=\theta_0-\theta_{\mathrm{requested}}.

Record the request, feasible start, and initialization correction before evaluating the first gradient. With no constraint, projection is the identity. With a box request $(2,-3)$ , the feasible start is $(1,-1)$ and its correction is $(-1,2)$ . This adjustment is initialization, not iteration 0’s gradient step.

For each iteration, let g be the gradient of the selected objective and retain

z_k=\theta_k-\eta g_k,\qquad \theta_{k+1}=P_C(z_k),\qquad q_k=\theta_{k+1}-z_k,\qquad d_k=\theta_{k+1}-\theta_k.

Here $q_k$ is the projection correction; $d_k$ is the actual parameter displacement. They satisfy $d_k=-\eta g_k+q_k$ . Drawing the proposal vector as the actual displacement would hide the constraint’s effect.

Reset restores $a=1,c=2$ , requested start $(0,1)$ , $\eta=1$ , $K=12$ , $\lambda=0$ , and the box. Replay selects stored updates. Compare the unconstrained preset at $\eta=0.1$ , the disk at the canonical rate, and the no-constraint penalty presets at $\lambda=1$ and $\lambda=3$ , both with rate $0.1$ .

Derivation: nearest points and the penalized minimum

The box distance objective separates into two coordinate squares. Each coordinate’s nearest allowed value is its clamp:

P_{\mathrm{box}}(z)= \bigl(\min(1,\max(-1,z_x)),\min(1,\max(-1,z_y))\bigr).

For the disk, an interior point already has zero projection distance. If $\norm z_2>1$ , the triangle inequality gives $\norm{z-u}_2\geq\norm z_2-1$ for every feasible u. Equality holds at $u=z/\norm z_2$ , so

P_{\mathrm{disk}}(z)=\frac{z}{\max(1,\norm z_2)}.

This also has a direct KKT check. Use Lagrangian $\frac12\norm{u-z}_2^2+\frac\mu2(\norm u_2^2-1)$ . For an exterior z, choose $u=z/\norm z_2$ and $\mu=\norm z_2-1$ . Then

u-z+\mu u=0,\qquad \norm u_2^2=1,\qquad \mu\geq0,\qquad \mu(\norm u_2^2-1)=0.

Feasibility, stationarity, and complementarity identify the unique nearest point for this convex projection problem. They concern distance to z, not minimizing L in one operation.

On the whole plane, the penalized quadratic has gradient $((a+2\lambda)x-ac,(4+2\lambda)y)$ . Its positive-definite Hessian gives the unique minimum

\theta_J^*=\left(\frac{ac}{a+2\lambda},0\right),\qquad J_{\min}=\frac{a\lambda c^2}{a+2\lambda}.

For no hard constraint, errors around this changed target multiply by $1-\eta(a+2\lambda)$ and $1-\eta(4+2\lambda)$ . Increasing the penalty increases curvature as well as changing the target; it need not rescue an excessive rate.

Worked examples: box repair, radial repair, and a soft preference

At $(0,1)$ with $a=1,c=2$ , the base loss is 4. The first box step is

g_0=(-2,4),\qquad z_0=(2,-3),\qquad \theta_1=(1,-1),\qquad L_1=5/2.

Its correction and actual displacement are different:

q_0=(1,-1)-(2,-3)=(-1,2),\qquad d_0=(1,-1)-(0,1)=(1,-2).

Worked example: Compare the proposal, projection and penalty

Parameter space: θ₁ and θ₂

Loss slice: θ₂ = 0

Current parameter vector: (0, 1)
Current objective value: 4
Actual gradient components: (−2, 4)
Learning rate: 1
Unconstrained update vector: (2, −4)
Next parameter vector: (1, −1)
New objective value: 2.500000
Unconstrained proposed parameters: (2, −3)
Projection correction: (−1, 2)
Actual projected displacement: (1, −2)
Feasible region: box [-1, 1]²
Requested start before feasibility projection: (0, 1)
Evaluated feasible start: (0, 1)
Configured L2 penalty strength: 0

≈ binary64 calculation · = analytic classification. Replay selects the same evaluated step in every view.

L=\tfrac a2(\theta_1-c)^2+2\theta_2^2

z_k=\theta_k-\eta\nabla J(\theta_k),\quad\theta_{k+1}=P_C(z_k)

z_{0}\approx \begin{bmatrix}0\\1\end{bmatrix}+\begin{bmatrix}2\\-4\end{bmatrix},\quad\theta_{1}=P_C(z_{0})\approx \begin{bmatrix}1\\-1\end{bmatrix}

J(\theta)=L(\theta)+\lambda\|\theta\|^2,\quad\nabla J=\nabla L+2\lambda\theta

Evaluated binary64 history. Replay changes the selected row, never the stored mathematics.
k	θ(k)	L(k)	gradient	η	proposed update	θ(k+1)	L(k+1)	unconstrained proposal	projection correction	base loss
0	(0, 1)	4	(−2, 4)	1	(2, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	4
1	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5
2	(1, 1)	2.5	(−1, 4)	1	(1, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	2.5
3	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5
4	(1, 1)	2.5	(−1, 4)	1	(1, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	2.5
5	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5
6	(1, 1)	2.5	(−1, 4)	1	(1, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	2.5
7	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5
8	(1, 1)	2.5	(−1, 4)	1	(1, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	2.5
9	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5
10	(1, 1)	2.5	(−1, 4)	1	(1, −4)	(1, −1)	2.5	(2, −3)	(−1, 2)	2.5
11	(1, −1)	2.5	(−1, −4)	1	(1, 4)	(1, 1)	2.5	(2, 3)	(−1, −2)	2.5

Evaluated iteration 0

z_{0}\approx \begin{bmatrix}0\\1\end{bmatrix}+\begin{bmatrix}2\\-4\end{bmatrix},\quad\theta_{1}=P_C(z_{0})\approx \begin{bmatrix}1\\-1\end{bmatrix}

SpeedStep 1 of 12

Replay the canonical box step: loss 4, gradient (−2,4), proposal (2,−3), projected point (1,−1), correction (−1,2), displacement (1,−2), and next base loss 5/2.

The next gradient is $(-1,-4)$ , proposal $(2,3)$ , and projected point $(1,1)$ . Subsequent y values alternate between $-1$ and 1 at this rate, with loss $5/2$ . The constrained minimizer is $(1,0)$ with loss $1/2$ . Every iterate is feasible, yet this run does not approach that minimum. At the minimizer the raw gradient is $(-1,0)$ , while projection makes the actual displacement zero.

With the same first proposal and the disk instead of the box,

\theta_1=\frac{(2,-3)}{\sqrt{13}},\qquad q_0=\left(\frac2{\sqrt{13}}-2,-\frac3{\sqrt{13}}+3\right),\qquad d_0=\left(\frac2{\sqrt{13}},-\frac3{\sqrt{13}}-1\right).

Its norm is 1, with projection multiplier $\sqrt{13}-1$ . Its next base loss is $46/13-4/\sqrt{13}\approx2.42906$ . Coordinate clamping to $(1,-1)$ would leave the disk.

For the penalty preset, use no hard constraint, $\lambda=1$ , and $\eta=1/10$ . At the same start, $J_0=5$ and $\nabla J=(-2,6)$ , giving

\theta_1=(1/5,2/5),\qquad L_1=97/50,\qquad J_1=107/50=2.14.

Here $\theta_J^*=(2/3,0)$ , $J_{\min}=4/3$ , and $H_J=\operatorname{diag}(3,6)$ . The stronger penalty $\lambda=3$ starts at $J_0=7$ , uses gradient $(-2,10)$ , and reaches $(1/5,0)$ with $J_1=87/50$ . Its preferred x coordinate is $2/7$ . These are changes to the optimization problem, not projections onto a hard region.

Executable reference code

The standalone standard-library reference uses 80-digit Decimal arithmetic for the bounded two-coordinate recurrence and disk square roots. It records the requested start, feasible start, initialization correction, and each step’s base and penalized evaluations, proposal, projection, correction, and displacement. Literal assertions check the complete box step, disk feasibility and KKT equations, and the penalty’s first loss and minimum.

For unconstrained penalty runs it independently compares the computed path with powers of the two error factors around the changed minimizer. The NumPy reference advances actual two-component arrays, then compares every finite recorded value against a separate scalar high-precision recurrence, including both projection mechanisms. Its absolute tolerance is $10^{-12}$ on the bounded checked runs.

code/python/optimization_constraints.py

"""Bounded two-coordinate projection and L2 penalty, with 80-digit Decimal checks."""
from decimal import Decimal as D, localcontext


def project(point, constraint):
    x, y = point
    if constraint == "none":
        return point
    if constraint == "box":
        return (min(D(1), max(D(-1), x)), min(D(1), max(D(-1), y)))
    if constraint == "disk":
        divisor = max(D(1), (x*x+y*y).sqrt())
        return (x/divisor, y/divisor)
    raise ValueError("Constraint must be none, box, or disk")


def evaluate(point, a, c, lam):
    x, y = point
    loss = a*(x-c)*(x-c)/2+2*y*y
    gradient = (a*(x-c), 4*y)
    return {"loss": loss, "penalized_loss": loss+lam*(x*x+y*y),
            "gradient": gradient, "penalized_gradient": (gradient[0]+2*lam*x, gradient[1]+2*lam*y),
            "hessian": ((a, D(0)), (D(0), D(4))),
            "penalized_hessian": ((a+2*lam, D(0)), (D(0), 4+2*lam))}


def trajectory(start=(0, 1), eta="1", iterations=12, a="1", c="2", lam="0", constraint="box"):
    if len(start) != 2 or type(iterations) is not int or not 1 <= iterations <= 40:
        raise ValueError("Require two coordinates and an integer budget in [1,40]")
    with localcontext() as context:
        context.prec = 80
        requested = tuple(D(str(value)) for value in start)
        eta, a, c, lam = (D(str(value)) for value in (eta, a, c, lam))
        if not all(value.is_finite() for value in (*requested, eta, a, c, lam)):
            raise ValueError("Require finite inputs")
        if eta <= 0 or not D("0.25") <= a <= 4 or not -2 <= c <= 2 or not 0 <= lam <= 4:
            raise ValueError("Require eta>0, a in [0.25,4], c in [-2,2], and lambda in [0,4]")
        point = feasible_start = project(requested, constraint)
        frames = []
        for k in range(iterations):
            current = evaluate(point, a, c, lam)
            gx, gy = current["penalized_gradient"]
            proposed = (point[0]-eta*gx, point[1]-eta*gy)
            projected = project(proposed, constraint)
            correction = (projected[0]-proposed[0], projected[1]-proposed[1])
            displacement = (projected[0]-point[0], projected[1]-point[1])
            frames.append({"k": k, "point": point, "evaluation": current,
                           "proposed": proposed, "projected": projected, "correction": correction,
                           "displacement": displacement, "next_evaluation": evaluate(projected, a, c, lam)})
            point = projected
        return {"requested": requested, "start": feasible_start,
                "initial_correction": (feasible_start[0]-requested[0], feasible_start[1]-requested[1]),
                "frames": frames}


def main():
    with localcontext() as context:
        context.prec = 80
        tolerance = D("1e-70")
        box = trajectory()
        first = box["frames"][0]
        assert first["evaluation"]["loss"] == 4
        assert first["evaluation"]["gradient"] == (-2, 4)
        assert first["proposed"] == (2, -3) and first["projected"] == (1, -1)
        assert first["correction"] == (-1, 2) and first["displacement"] == (1, -2)
        assert first["next_evaluation"]["loss"] == D("2.5")
        assert box["frames"][1]["proposed"] == (2, 3)
        for k, frame in enumerate(box["frames"]):
            assert frame["projected"] == (1, D(-1)**(k+1))
            assert frame["next_evaluation"]["loss"] == D("2.5")
            assert all(-1 <= value <= 1 for value in frame["projected"])
        initialized = trajectory((2, -3), iterations=1)
        assert initialized["requested"] == (2, -3)
        assert initialized["start"] == (1, -1) and initialized["initial_correction"] == (-1, 2)
        assert trajectory((2, -3), iterations=1, constraint="none")["start"] == (2, -3)
        boundary = trajectory((1, 0), iterations=1)["frames"][0]
        assert boundary["evaluation"]["gradient"] == (-1, 0)
        assert boundary["displacement"] == (0, 0) and boundary["evaluation"]["loss"] == D("0.5")

        disk = trajectory(constraint="disk")
        radial = disk["frames"][0]
        root = D(13).sqrt()
        projected = radial["projected"]
        assert projected == (2/root, -3/root)
        assert radial["correction"] == (2/root-2, -3/root+3)
        assert radial["displacement"] == (2/root, -3/root-1)
        squared_norm = projected[0]**2+projected[1]**2
        multiplier = root-1
        assert multiplier > 0 and abs(squared_norm-1) < tolerance
        assert all(abs(p-z+multiplier*p) < tolerance for p, z in zip(projected, (2, -3)))
        assert abs(multiplier*(squared_norm-1)) < tolerance
        distance = (projected[0]-2)**2+(projected[1]+3)**2
        assert abs(distance-(root-1)**2) < tolerance
        assert abs(radial["next_evaluation"]["loss"]-(D(46)/13-4/root)) < tolerance
        for frame in disk["frames"]:
            x, y = frame["projected"]
            assert x*x+y*y <= 1+tolerance
        assert project((D(0), D(0)), "disk") == (0, 0)
        assert project((D("0.2"), D("0.4")), "disk") == (D("0.2"), D("0.4"))
        disk_start = trajectory((2, -3), iterations=1, constraint="disk")
        assert disk_start["start"] == projected
        assert disk_start["initial_correction"] == radial["correction"]

        for lam in (0, 1, 3):
            run = trajectory(eta="0.1", lam=lam, constraint="none")
            target = D(2)/(1+2*lam)
            rx, ry = 1-D("0.1")*(1+2*lam), 1-D("0.1")*(4+2*lam)
            for k, frame in enumerate(run["frames"], 1):
                expected = (target-target*rx**k, ry**k)
                assert all(abs(actual-value) < tolerance for actual, value in zip(frame["projected"], expected))
        plain = trajectory(eta="0.1", constraint="none")["frames"][0]
        assert plain["projected"] == (D("0.2"), D("0.6"))
        assert plain["next_evaluation"]["loss"] == D("2.34")
        penalty = trajectory(eta="0.1", lam=1, constraint="none")["frames"][0]
        assert penalty["evaluation"]["penalized_loss"] == 5
        assert penalty["evaluation"]["penalized_gradient"] == (-2, 6)
        assert penalty["evaluation"]["penalized_hessian"] == ((3, 0), (0, 6))
        assert penalty["projected"] == (D("0.2"), D("0.4"))
        assert penalty["next_evaluation"]["loss"] == D("1.94")
        assert penalty["next_evaluation"]["penalized_loss"] == D("2.14")
        minimum = evaluate((D(2)/3, D(0)), D(1), D(2), D(1))
        assert abs(minimum["penalized_loss"]-D(4)/3) < tolerance
        assert all(abs(value) < tolerance for value in minimum["penalized_gradient"])
        strong = trajectory(eta="0.1", lam=3, constraint="none")["frames"][0]
        assert strong["evaluation"]["penalized_loss"] == 7
        assert strong["evaluation"]["penalized_gradient"] == (-2, 10)
        assert strong["projected"] == (D("0.2"), 0)
        assert strong["next_evaluation"]["penalized_loss"] == D("1.74")
        assert D(4)*2/(4+2*1) > 1  # Finite lambda=1 need not enforce either hard region.

        for kwargs in ({"a": 0}, {"c": 3}, {"lam": -1}, {"eta": 0},
                       {"iterations": 41}, {"start": ("NaN", 1)}, {"constraint": "triangle"}):
            try:
                trajectory(**kwargs)
            except ValueError:
                pass
            else:
                raise AssertionError("Invalid bounded input accepted")
        print("Decimal box: proposal (2,-3), projection (1,-1), correction (-1,2), displacement (1,-2), L=2.5.")
        print("Decimal disk: (2,-3)/sqrt(13); feasibility, nearest-point distance, and KKT residuals checked.")
        print("Penalty: first J=2.14; minimum at (2/3,0), J=4/3; closed-form paths checked.")
        print("Feasibility does not imply convergence; finite L2 penalty does not enforce a hard region.")


if __name__ == "__main__":
    main()

code/python/optimization_constraints_numpy.py

"""Actual two-coordinate NumPy recurrence, checked against independent 80-digit scalars."""
from decimal import Decimal as D, localcontext
import numpy as np


def project(point, constraint):
    if constraint == "none":
        return point.copy()
    if constraint == "box":
        return np.clip(point, -1.0, 1.0)
    if constraint == "disk":
        return point/max(1.0, float(np.hypot(point[0], point[1])))
    raise ValueError("Constraint must be none, box, or disk")


def evaluate(point, a, c, lam):
    x, y = point
    loss = a*(x-c)*(x-c)/2+2*y*y
    gradient = np.array([a*(x-c), 4*y])
    penalized_loss = loss+lam*np.dot(point, point)
    penalized_gradient = gradient+2*lam*point
    if not np.all(np.isfinite([loss, penalized_loss, *gradient, *penalized_gradient])):
        raise OverflowError("Non-finite binary64 objective or gradient")
    return {"loss": float(loss), "penalized_loss": float(penalized_loss), "gradient": gradient,
            "penalized_gradient": penalized_gradient, "hessian": np.diag([a, 4.0]),
            "penalized_hessian": np.diag([a+2*lam, 4+2*lam])}


def trajectory(start=(0, 1), eta=1.0, iterations=12, a=1.0, c=2.0, lam=0.0, constraint="box"):
    requested = np.asarray(start, dtype=np.float64).copy()
    if requested.shape != (2,) or not np.all(np.isfinite(requested)):
        raise ValueError("Require a finite two-coordinate start")
    if type(iterations) is not int or not 1 <= iterations <= 40:
        raise ValueError("Require an integer budget in [1,40]")
    if not np.all(np.isfinite([eta, a, c, lam])) or eta <= 0 or not 0.25 <= a <= 4 or not -2 <= c <= 2 or not 0 <= lam <= 4:
        raise ValueError("Require finite parameters, eta>0, a in [0.25,4], c in [-2,2], lambda in [0,4]")
    point = feasible_start = project(requested, constraint)
    frames = []
    with np.errstate(over="ignore", invalid="ignore"):
        for k in range(iterations):
            current = evaluate(point, a, c, lam)
            proposed = point-eta*current["penalized_gradient"]
            if not np.all(np.isfinite(proposed)):
                raise OverflowError("Non-finite binary64 proposal")
            projected = project(proposed, constraint)
            frames.append({"k": k, "point": point, "evaluation": current, "proposed": proposed,
                           "projected": projected, "correction": projected-proposed,
                           "displacement": projected-point, "next_evaluation": evaluate(projected, a, c, lam)})
            point = projected
    return {"requested": requested, "start": feasible_start,
            "initial_correction": feasible_start-requested, "frames": frames}


def main():
    def check(actual, expected):
        np.testing.assert_allclose(actual, np.asarray(expected, dtype=np.float64), rtol=0, atol=1e-12)

    first = trajectory()["frames"][0]
    np.testing.assert_array_equal(first["evaluation"]["gradient"], [-2, 4])
    for key, literal in (("proposed", [2, -3]), ("projected", [1, -1]),
                         ("correction", [-1, 2]), ("displacement", [1, -2])):
        np.testing.assert_array_equal(first[key], literal)
    assert first["evaluation"]["loss"] == 4 and first["next_evaluation"]["loss"] == 2.5
    initialized = trajectory((2, -3), iterations=1)
    np.testing.assert_array_equal(initialized["start"], [1, -1])
    np.testing.assert_array_equal(initialized["initial_correction"], [-1, 2])
    disk = trajectory(constraint="disk")["frames"][0]
    root = np.sqrt(13.0)
    check(disk["projected"], np.array([2, -3])/root)
    assert abs(np.dot(disk["projected"], disk["projected"])-1) < 1e-14
    check(disk["projected"]-disk["proposed"]+(root-1)*disk["projected"], [0, 0])
    check(np.dot(disk["correction"], disk["correction"]), (root-1)**2)
    check(project(np.array([0.0, 0.0]), "disk"), [0, 0])

    cases = [
        {"constraint": "box"}, {"constraint": "disk"}, {"constraint": "none", "eta": 0.1},
        {"constraint": "none", "lam": 1, "eta": 0.1},
        {"constraint": "none", "lam": 3, "eta": 0.1},
        {"constraint": "disk", "start": (2, -3)},
        {"constraint": "box", "a": 0.25, "c": -2, "lam": 4, "eta": 0.05},
        {"constraint": "none", "a": 4, "c": 2, "lam": 1, "eta": 0.1},
    ]
    with localcontext() as context:
        context.prec = 80
        for settings in cases:
            result = trajectory(**settings)
            a, c, lam, eta = (D(str(settings.get(key, default))) for key, default in
                              (("a", 1), ("c", 2), ("lam", 0), ("eta", 1)))
            constraint = settings["constraint"]
            requested = tuple(D(str(value)) for value in settings.get("start", (0, 1)))

            def scalar_project(x, y):
                if constraint == "box":
                    return min(D(1), max(D(-1), x)), min(D(1), max(D(-1), y))
                if constraint == "disk" and x*x+y*y > 1:
                    radius = (x*x+y*y).sqrt()
                    return x/radius, y/radius
                return x, y

            x, y = scalar_project(*requested)
            check(result["requested"], requested)
            check(result["start"], (x, y))
            check(result["initial_correction"], (x-requested[0], y-requested[1]))
            for frame in result["frames"]:
                check(frame["point"], (x, y))
                gx, gy = (a+2*lam)*x-a*c, (4+2*lam)*y
                zx, zy = x-eta*gx, y-eta*gy
                nx, ny = scalar_project(zx, zy)
                check(frame["proposed"], (zx, zy))
                check(frame["projected"], (nx, ny))
                check(frame["correction"], (nx-zx, ny-zy))
                check(frame["displacement"], (nx-x, ny-y))
                for evaluation, px, py in ((frame["evaluation"], x, y), (frame["next_evaluation"], nx, ny)):
                    expanded_loss = a*px*px/2-a*c*px+a*c*c/2+2*py*py
                    check(evaluation["loss"], expanded_loss)
                    check(evaluation["penalized_loss"], expanded_loss+lam*(px*px+py*py))
                    check(evaluation["gradient"], (a*px-a*c, 4*py))
                    check(evaluation["penalized_gradient"], ((a+2*lam)*px-a*c, (4+2*lam)*py))
                    check(evaluation["hessian"], [[a, 0], [0, 4]])
                    check(evaluation["penalized_hessian"], [[a+2*lam, 0], [0, 4+2*lam]])
                if constraint == "box":
                    assert np.all(np.abs(frame["projected"]) <= 1)
                if constraint == "disk":
                    assert np.hypot(*frame["projected"]) <= 1+1e-14
                x, y = nx, ny

    penalty = trajectory(eta=0.1, lam=1, constraint="none")["frames"][0]
    check(penalty["projected"], [0.2, 0.4])
    check(penalty["next_evaluation"]["penalized_loss"], 2.14)
    optimum = evaluate(np.array([2/3, 0.0]), 1, 2, 1)
    check(optimum["penalized_gradient"], [0, 0])
    check(optimum["penalized_loss"], 4/3)
    check(optimum["penalized_hessian"], [[3, 0], [0, 6]])
    outside = np.array([4/3, 0.0])
    check(evaluate(outside, 4, 2, 1)["penalized_gradient"], [0, 0])
    assert np.hypot(*outside) > 1 and outside[0] > 1
    for settings in ({"start": [0, np.inf]}, {"a": 0}, {"lam": -1}, {"constraint": "triangle"},
                     {"iterations": 41}, {"eta": 0}, {"c": 3}):
        try:
            trajectory(**settings)
        except ValueError:
            pass
        else:
            raise AssertionError("Invalid bounded input accepted")
    print("NumPy box, disk, unconstrained, and penalty recurrences passed all recorded-value checks.")
    print("Independent 80-digit scalar reference: absolute tolerance 1e-12 on eight bounded runs.")
    print("Requested starts, nearest-point/KKT residuals, and the finite-penalty counterexample checked.")


if __name__ == "__main__":
    main()

Explore: predict the correction before comparing paths

Request a start outside the box and locate the initialization correction before reading iteration 0. Repeat with the disk: an exterior request preserves its direction under radial scaling. An interior request needs no correction. Switching to no constraint preserves the original request.

For the canonical first proposal, compare box and disk endpoints without changing the current gradient. Identify which correction component comes from a coordinate boundary and which comes from radial scaling. Next reduce the box rate and compare feasible progress with the canonical alternating run.

Remove the hard constraint and compare penalties 0, 1, and 3 at rate $0.1$ . Track J and L separately, and predict each preferred x coordinate from $ac/(a+2\lambda)$ . As a counterexample to guaranteed hard feasibility, set $a=4,c=2,\lambda=1$ : the penalized minimizer is $(4/3,0)$ , outside both fixed regions.

Common mistakes and counterexamples

Projection acts on the proposal, not the gradient coordinates or the original requested start forever. The box clamp is not the disk rule. Projection minimizes Euclidean distance to its input; it does not solve the constrained loss problem in a single update.

The penalty contributes $2\lambda\theta$ to the gradient and $2\lambda I$ to the Hessian. Adding it only to a displayed loss would leave the wrong optimizer recurrence. Likewise, calling J the base loss would conceal the preference being traded against data fit.

A feasible path can oscillate, and a constrained minimum can have nonzero raw gradient. At $(1,0)$ in this example, the negative gradient points beyond the allowed region and projection returns the same point. An unconstrained raw-gradient stopping rule would miss that constrained minimum.

Numerical considerations

The stage and NumPy path use binary64 approximations. Disk normalization generally involves an irrational norm, so displayed coordinates use approximation signs. A computed boundary norm can differ slightly from 1 through rounding; the references check feasibility and KKT residuals with stated numerical tolerances.

Retain the proposed and projected values even when rounded labels appear identical. A small projection correction or small actual displacement can disappear at display precision. A large raw gradient at a boundary need not imply a large feasible displacement.

The high-precision reference supplies an independent numerical check, not symbolic exactness for every square root. Its positive curvature, nonnegative penalty, finite inputs, and bounded iteration count are validated. Comparing J between different penalty settings also changes the quantity being compared.

Engineering and ML connections

Hard parameter limits can encode an admissible physical range or a norm budget. L2 regularization instead adds a cost for parameter size, changing the balance of the modeled objective. The coefficient and parameter units matter; the same numeric penalty after rescaling parameters describes a different preference.

Learning-rate behavior explains how the added curvature affects a plain update. Least squares separates fit from residual geometry, while convexity supplies the structure behind this toy problem’s unique minima. Numerical diagnostics monitor a run’s arithmetic and stopping reasons; regularizing a model and protecting a numerical update are separate decisions.

Exercises

Keep the requested start, feasible state, unconstrained proposal, projection correction, and actual displacement distinct. Use exact fractions for rational hand examples and the radial expression for the disk.

Summary and checkpoints

Projection enforces an allowed region by finding the nearest feasible point to each proposal. A penalty changes the objective, its gradient, its curvature, and its minimizer. Check that you can reconstruct the canonical box step, justify radial disk projection, record an infeasible requested start honestly, and recover J’s minimum $4/3$ at $(2/3,0)$ . Neither feasible iterates nor a finite penalty alone establish convergence.

Do we restrict the allowed points or change their cost?

Compare the proposal, projection and penalty

Parameter space: θ₁ and θ₂

Loss slice: θ₂ = 0

Evaluated iteration 0

Intuition: a boundary and a preference are different rules

Definition: the feasible set and the changed objective

Mechanics: record initialization and both parts of a step

Derivation: nearest points and the penalized minimum

Worked examples: box repair, radial repair, and a soft preference

Worked example: Compare the proposal, projection and penalty

Parameter space: θ₁ and θ₂

Loss slice: θ₂ = 0

Evaluated iteration 0

Executable reference code

Explore: predict the correction before comparing paths

Common mistakes and counterexamples

Numerical considerations

Engineering and ML connections

Exercises

Separate the proposal from the feasible step

Record an infeasible requested start

Differentiate the changed objective

Locate the penalized minimum

Distinguish nearest-point projection from a penalty

Summary and checkpoints

The map of the book

Mathematical Language

Linear Algebra

Calculus

Probability

Statistics

Optimization

Information Theory

Numerical Computing