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 with python3 projection.py, from any directory — it writes its output files into the current directory. Requires NumPy and Matplotlib only, and uses fixed seeds, so the printed numbers reproduce exactly. Download projection.py
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. 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 projection violates additivity — it is not
linear (the entry's closing remark).
[P-projgd] Projected GD for Lasso: gradient steps followed by l1-ball
projections converge to a feasible iterate whose empirical
risk is (near-)optimal among feasible points, enforcing the
constraint ||w||_1 <= tau in every iteration.
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
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. 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 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))
# 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: all near-closest points cluster at 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]))
# l1-ball projection is NOT linear
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] 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: all near-closest points cluster at the projection [ok] nonconvex closed set: a closest point exists but is not unique [ok] l1-ball projection violates additivity (not a linear map)
Projected GD for Lasso: gradient steps followed by l1-ball projections converge to a feasible iterate whose empirical risk is (near-)optimal among feasible points, enforcing the constraint ||w||_1 <= tau in every iteration.
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
risk = 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)
risks = np.mean((y[None, :] - cand @ X.T) ** 2, axis=1)
check("every iterate stayed feasible", feasible_all)
check("projected-GD risk <= best sampled feasible risk + 1e-3",
risk(w) <= risks.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_aspect("equal"); ax.set_title("[P-def] projection onto the l1-ball")
fig.tight_layout()
fig.savefig("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] projected-GD risk <= best sampled feasible risk + 1e-3 10/10 checks passed
P-projgd writes when the script runs