P-projgd writes when the script runsDictionary of Applied Machine Learning · projection
Numerical companion to the entry projection: it recomputes what the entry states and prints one line per check
One block per paragraph of the entry (marked [P...]): each block verifies numerically what the corresponding statement asserts. Self-contained (numpy/matplotlib only), fixed seed.
Run it without installing anything:uv run https://dictionaryofml.org/terms/projection.py
uv downloads this script and the pinned NumPy and Matplotlib it needs, then runs it; the script fetches any input file it uses. To keep the output files, download projection.py into a folder and run uv run projection.py there. With NumPy and Matplotlib already installed, python3 projection.py, from any directory — it writes its output files into the current directory. Fixed seeds, so the printed numbers reproduce exactly. Download projection.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""
projection.py — numerical companion to the glossary entry 'projection'.
One block per paragraph of the entry (marked [P...]): each block verifies
numerically what the corresponding statement asserts. Self-contained
(numpy/matplotlib only), fixed seed.
Blocks
------
[P-def] The projection onto a closed set is a closest point: the exact
l1-ball projection (via sorting-based soft threshold) beats a
dense sample of other points of the set in Euclidean
distance; for the convex l1-ball it is unique, and the
segment w - proj(w) meets the face of the ball at a right
angle (the marker in the entry's figure). On the two-point
set {(-1,0), (1,0)} the closest point is not unique for
(0, 0.7) — nor for any other point of the vertical axis,
which consists of exactly the points equidistant from the
two elements. For a subspace, the projection map is linear
(orthogonal projection matrix P = B (B^T B)^{-1} B^T with
P^2 = P), while the l1-ball is not a subspace ((tau, 0) is
in the ball, 2 (tau, 0) is not) and its projection violates
additivity — it is not linear (the entry's closing remark).
[P-fund] Idempotence and self-adjointness: every projection is
idempotent (also the nonlinear l1-ball projection); the
subspace projection matrix satisfies P^2 = P and P^T = P;
an oblique projection (idempotent, not self-adjoint) sends
a vector to a point of its range that is farther away than
the orthogonal projection foot.
[P-projgd] Projected GD for Lasso: gradient steps on the training error
followed by l1-ball projections converge to a feasible
iterate whose training error is (near-)optimal among
feasible points, enforcing ||w||_1 <= tau in every
iteration. The projection step is what enforces the
constraint: the same gradient steps without it produce
iterates with ||w||_1 > tau.
Outputs
-------
projection.png : preview figure (checking only).
Data generated by pythondemos/projection.py.
"""
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from pathlib import Path
OUT_DIR = Path(__file__).parent
rng = np.random.default_rng(42)
report = []
def check(name, ok):
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
def proj_l1(w, tau):
"""Exact Euclidean projection onto the l1-ball of radius tau."""
if np.abs(w).sum() <= tau:
return w.copy()
u = np.sort(np.abs(w))[::-1]
css = np.cumsum(u)
rho = np.nonzero(u * np.arange(1, len(w) + 1) > css - tau)[0][-1]
theta = (css[rho] - tau) / (rho + 1.0)
return np.sign(w) * np.maximum(np.abs(w) - theta, 0.0)
The projection onto a closed set is a closest point: the exact l1-ball projection (via sorting-based soft threshold) beats a dense sample of other points of the set in Euclidean distance; for the convex l1-ball it is unique, and the segment w - proj(w) meets the face of the ball at a right angle (the marker in the entry's figure). On the two-point set {(-1,0), (1,0)} the closest point is not unique for (0, 0.7) — nor for any other point of the vertical axis, which consists of exactly the points equidistant from the two elements. For a subspace, the projection map is linear (orthogonal projection matrix P = B (B^T B)^{-1} B^T with P^2 = P), while the l1-ball is not a subspace ((tau, 0) is in the ball, 2 (tau, 0) is not) and its projection violates additivity — it is not linear (the entry's closing remark).
print("[P-def] projection = closest point; linear only on subspaces")
tau = 1.0
w0 = np.array([1.6, 1.1])
p = proj_l1(w0, tau)
check("projection is feasible (||p||_1 <= tau)",
np.abs(p).sum() <= tau + 1e-12)
# closest point: sample the l1-ball densely, none is closer
angles = rng.uniform(0, 2 * np.pi, 20000)
radii = rng.uniform(0, 1, 20000)
raw = np.stack([np.cos(angles), np.sin(angles)], axis=1)
ball = radii[:, None] * raw / np.abs(raw).sum(axis=1, keepdims=True)
dists = np.linalg.norm(ball - w0, axis=1)
check("no sampled point of the ball is closer than the projection",
np.all(dists >= np.linalg.norm(p - w0) - 1e-9))
# right angle (figure marker): w - proj(w) is orthogonal to the face
# x_1 + x_2 = tau of the ball that contains the projection
face = np.array([1.0, -1.0])
check("w - proj(w) is orthogonal to the face containing proj(w)",
abs((w0 - p) @ face) < 1e-9)
# subspace: orthogonal projection matrix, linear and idempotent
B = rng.normal(size=(4, 2)) # columns span a 2-d subspace
P = B @ np.linalg.solve(B.T @ B, B.T)
u1, u2 = rng.normal(size=4), rng.normal(size=4)
check("subspace projection is linear: P(u + u') = P u + P u'",
np.allclose(P @ (u1 + u2), P @ u1 + P @ u2))
check("idempotent: P^2 = P", np.allclose(P @ P, P))
check("residual orthogonal to the subspace: B^T (u - P u) = 0",
np.max(np.abs(B.T @ (u1 - P @ u1))) < 1e-10)
# uniqueness on the convex ball: every near-minimizer is near p
near = ball[dists <= np.linalg.norm(p - w0) + 1e-3]
check("convex set: every near-closest point lies near the projection",
np.all(np.linalg.norm(near - p, axis=1) < 0.15))
# existence on a nonconvex closed set (two points): minimum attained,
# but NOT unique for the midpoint
S = np.array([[1.0, 0.0], [-1.0, 0.0]])
mid = np.array([0.0, 0.7])
d_mid = np.linalg.norm(S - mid, axis=1)
check("nonconvex closed set: a closest point exists but is not unique",
np.isclose(d_mid[0], d_mid[1]))
# ... and likewise for every point of the vertical axis, which consists
# of exactly the points equidistant from the two elements
axis_pts = np.stack([np.zeros(5), np.linspace(-2.0, 2.0, 5)], axis=1)
check("every point of the vertical axis is equidistant from both elements",
np.allclose(np.linalg.norm(axis_pts - S[0], axis=1),
np.linalg.norm(axis_pts - S[1], axis=1)))
off_axis = np.array([0.3, 0.7])
check("a point off the vertical axis has a unique closest element",
not np.isclose(np.linalg.norm(off_axis - S[0]),
np.linalg.norm(off_axis - S[1])))
# the l1-ball is not a subspace ((tau,0) is in the ball, 2(tau,0) is
# not), and its projection is NOT linear
edge = np.array([tau, 0.0])
check("l1-ball is not a subspace: (tau,0) inside, 2(tau,0) outside",
np.abs(edge).sum() <= tau and np.abs(2 * edge).sum() > tau)
a1, a2 = np.array([1.5, 0.0]), np.array([0.0, 1.5])
check("l1-ball projection violates additivity (not a linear map)",
not np.allclose(proj_l1(a1 + a2, tau),
proj_l1(a1, tau) + proj_l1(a2, tau)))
[P-def] projection = closest point; linear only on subspaces [ok] projection is feasible (||p||_1 <= tau) [ok] no sampled point of the ball is closer than the projection [ok] w - proj(w) is orthogonal to the face containing proj(w) [ok] subspace projection is linear: P(u + u') = P u + P u' [ok] idempotent: P^2 = P [ok] residual orthogonal to the subspace: B^T (u - P u) = 0 [ok] convex set: every near-closest point lies near the projection [ok] nonconvex closed set: a closest point exists but is not unique [ok] every point of the vertical axis is equidistant from both elements [ok] a point off the vertical axis has a unique closest element [ok] l1-ball is not a subspace: (tau,0) inside, 2(tau,0) outside [ok] l1-ball projection violates additivity (not a linear map)
Idempotence and self-adjointness: every projection is idempotent (also the nonlinear l1-ball projection); the subspace projection matrix satisfies P^2 = P and P^T = P; an oblique projection (idempotent, not self-adjoint) sends a vector to a point of its range that is farther away than the orthogonal projection foot.
print("[P-fund] idempotent + self-adjoint characterizes orthogonal projection")
check("l1-ball projection is idempotent: proj(proj(w)) = proj(w)",
np.allclose(proj_l1(proj_l1(w0, tau), tau), proj_l1(w0, tau)))
check("subspace projection matrix: P^2 = P and P^T = P",
np.allclose(P @ P, P) and np.allclose(P.T, P))
# near miss: an oblique projection is idempotent but not self-adjoint,
# and its output is not a closest point of its range (the x-axis)
A = np.array([[1.0, 1.0], [0.0, 0.0]])
u = np.array([0.3, 0.8])
foot = np.array([u[0], 0.0])
check("oblique projection: idempotent but not self-adjoint",
np.allclose(A @ A, A) and not np.allclose(A.T, A))
check("oblique output is farther from u than the orthogonal foot",
np.linalg.norm(A @ u - u) > np.linalg.norm(foot - u) + 1e-9)
[P-fund] idempotent + self-adjoint characterizes orthogonal projection [ok] l1-ball projection is idempotent: proj(proj(w)) = proj(w) [ok] subspace projection matrix: P^2 = P and P^T = P [ok] oblique projection: idempotent but not self-adjoint [ok] oblique output is farther from u than the orthogonal foot
Projected GD for Lasso: gradient steps on the training error followed by l1-ball projections converge to a feasible iterate whose training error is (near-)optimal among feasible points, enforcing ||w||_1 <= tau in every iteration. The projection step is what enforces the constraint: the same gradient steps without it produce iterates with ||w||_1 > tau.
print("[P-projgd] projected GD solves the Lasso constraint form")
m, d = 40, 2
X = rng.normal(size=(m, d))
y = X @ np.array([1.2, -0.3]) + 0.05 * rng.normal(size=m)
L = 2 * np.linalg.eigvalsh(X.T @ X / m).max()
w = np.zeros(d)
feasible_all = True
for _ in range(400):
w = proj_l1(w - (1 / L) * (2 / m) * X.T @ (X @ w - y), tau)
feasible_all &= np.abs(w).sum() <= tau + 1e-10
trainerr = lambda v: np.mean((y - X @ v) ** 2)
# compare against a dense sample of the feasible set
cand = tau * ball / np.maximum(np.abs(ball).sum(axis=1, keepdims=True), 1e-12)
trainerrs = np.mean((y[None, :] - cand @ X.T) ** 2, axis=1)
check("every iterate stayed feasible", feasible_all)
# the cost of omitting the projection: plain GD violates the constraint
w_plain = np.zeros(d)
for _ in range(400):
w_plain = w_plain - (1 / L) * (2 / m) * X.T @ (X @ w_plain - y)
check("without the projection step, iterates violate ||w||_1 <= tau",
np.abs(w_plain).sum() > tau + 1e-10)
check("projected-GD training error <= best sampled feasible one + 1e-3",
trainerr(w) <= trainerrs.min() + 1e-3)
# ------------------------------------------------------------ preview
fig, ax = plt.subplots(figsize=(4.2, 4.0))
square = np.array([[1, 0], [0, 1], [-1, 0], [0, -1], [1, 0]]) * tau
ax.plot(square[:, 0], square[:, 1], "k-")
ax.plot(*w0, "ko"); ax.annotate("w", w0)
ax.plot(*p, "rs"); ax.annotate("proj(w)", p)
ax.plot([w0[0], p[0]], [w0[1], p[1]], "k--")
ax.plot(*w, "b^"); ax.annotate("projected GD", w)
ax.set_xlabel("$w_1$"); ax.set_ylabel("$w_2$")
ax.set_aspect("equal"); ax.set_title("[P-def] projection onto the l1-ball")
fig.tight_layout()
fig.savefig(OUT_DIR / "projection.png", dpi=110)
print(f"\n{sum(ok for _, ok in report)}/{len(report)} checks passed")
assert all(ok for _, ok in report)
[P-projgd] projected GD solves the Lasso constraint form [ok] every iterate stayed feasible [ok] without the projection step, iterates violate ||w||_1 <= tau [ok] projected-GD training error <= best sampled feasible one + 1e-3 19/19 checks passed
P-projgd writes when the script runs