Dictionary of Applied Machine Learning · Sobolev space
Numerical companion to the entry Sobolev space: it recomputes what the entry states and prints one line per check
Backs every claim of the entry: the defining integration-by-parts identity of the weak derivative, the equivalence "bounded weak gradient <=> Lipschitz", the spectral-norm bound on the weak gradient of a ReLU network, the certified radius that a bounded weak gradient buys, the failure of a piecewise-constant hypothesis to have any such radius, and the graph counterpart used by GTVMin. Self-contained (numpy only), fixed seeds. Run from the repo root or from pythondemos/ — no files are written.
Run it with python3 sobolevspace.py, from any directory — it writes no files and only prints its checks. Requires NumPy only, and uses fixed seeds, so the printed numbers reproduce exactly. Download sobolevspace.py
One cell per block of the script: the code, and what that code printed when it last ran here
"""
sobolevspace.py — numerical companion to the glossary entry 'Sobolev space'.
Purpose
-------
Backs every claim of the entry: the defining integration-by-parts
identity of the weak derivative, the equivalence "bounded weak
gradient <=> Lipschitz", the spectral-norm bound on the weak gradient
of a ReLU network, the certified radius that a bounded weak gradient
buys, the failure of a piecewise-constant hypothesis to have any such
radius, and the graph counterpart used by GTVMin. Self-contained
(numpy only), fixed seeds. Run from the repo root or from
pythondemos/ — no files are written.
Blocks
------
[P-weakderiv] For the kinked function f(x) = max(0, x) on (-2, 2),
the step g = 1_{x>0} satisfies the defining identity
int f(x) phi'(x) dx = - int g(x) phi(x) dx for smooth
test functions phi vanishing at the boundary, even
though f is not differentiable at x = 0. g is not a
weak derivative of the jump function 1_{x>0}: there the
same identity fails for every locally integrable g,
which the block shows by the jump term it produces.
[P-lipschitz] For a 2-layer ReLU network on R^2 the supremum of the
weak gradient norm and the largest difference quotient
|f(x) - f(x')| / ||x - x'|| agree to three decimals:
a bounded weak gradient IS the Lipschitz constant.
[P-spectral] The product of the spectral norms of the weight
matrices bounds that constant (valid but loose:
4.94 against 2.12 here).
[P-certified] If |f(x)| = gamma and f is L-Lipschitz, no perturbation
of norm < gamma / L changes sign(f(x)): 0 sign flips in
60000 perturbations drawn inside the certified radius.
For a linear hypothesis the weak gradient is the
constant w, so gamma / L is the distance of the feature
vector from the decision boundary (the SVM statement).
[P-jump] A decision stump is piecewise constant with a jump, so
it has no bounded weak gradient: a perturbation of norm
2e-9 flips its prediction, and no certified radius
exists at any margin.
[P-tvh1] Sharpening a transition of width eps: the squared H^1
seminorm grows like 4 / (3 eps) while the total variation
stays at the height of the limiting jump — why the L^1
gradient penalty tolerates jumps and the H^1 penalty does
not.
[P-graph] Graph counterpart: for a graph signal sampled from a
smooth function the edge-difference energy (the squared
graph gradient summed over edges, i.e. the Laplacian
quadratic form used by GTVMin) is small, while a signal
with one jump makes it large — the discrete analogue of
a bounded weak derivative.
"""
import numpy as np # the only numerical dependency
report = [] # collects (check name, pass/fail) pairs
def check(name, ok): # records and prints one verification
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
For the kinked function f(x) = max(0, x) on (-2, 2), the step g = 1_{x>0} satisfies the defining identity int f(x) phi'(x) dx = - int g(x) phi(x) dx for smooth test functions phi vanishing at the boundary, even though f is not differentiable at x = 0. g is not a weak derivative of the jump function 1_{x>0}: there the same identity fails for every locally integrable g, which the block shows by the jump term it produces.
# f(x) = max(0, x) has no derivative at x = 0, but the step g = 1_{x>0}
# satisfies the definition of the weak derivative: for every smooth phi
# vanishing near the boundary, int f phi' = - int g phi.
grid = np.linspace(-2.0, 2.0, 400001)
dx = grid[1] - grid[0]
f_kink = np.maximum(grid, 0.0) # continuous, kinked
g_step = (grid > 0).astype(float) # candidate weak derivative
f_jump = g_step.copy() # discontinuous at 0
# test functions phi(x) = (4 - x^2)^3 * x^p vanish (with derivative) at +-2
for p in (0, 1, 2):
phi = (4.0 - grid ** 2) ** 3 * grid ** p
dphi = np.gradient(phi, dx)
lhs = np.trapezoid(f_kink * dphi, grid) # int f phi'
rhs = -np.trapezoid(g_step * phi, grid) # - int g phi
check(f"[P-weakderiv] integration-by-parts identity holds for phi_{p} "
f"(lhs {lhs:+.4f} = rhs {rhs:+.4f})", abs(lhs - rhs) < 1e-3)
# the jump function admits no locally integrable weak derivative: the
# identity would force int f_jump phi' = -int g phi for all phi, but the
# left-hand side carries the value -phi(0) of the jump, which no integrable
# g can reproduce (its contribution shrinks with the width of phi's support).
defects = []
for eps in (0.4, 0.2, 0.1): # bumps of shrinking width
bump = np.where(np.abs(grid) < eps,
np.exp(-1.0 / np.maximum(1e-12, 1 - (grid / eps) ** 2)), 0.0)
bump /= bump.max()
dbump = np.gradient(bump, dx)
defects.append(np.trapezoid(f_jump * dbump, grid)) # -> -phi(0) = -1
check("[P-weakderiv] jump function: int f phi' stays at the jump height "
f"(-1) as the bump narrows ({', '.join(f'{d:+.3f}' for d in defects)}) "
"— no integrable weak derivative",
all(abs(d + 1.0) < 0.02 for d in defects))
[ok] [P-weakderiv] integration-by-parts identity holds for phi_0 (lhs -58.5143 = rhs -58.5140) [ok] [P-weakderiv] integration-by-parts identity holds for phi_1 (lhs -32.0000 = rhs -32.0000) [ok] [P-weakderiv] integration-by-parts identity holds for phi_2 (lhs -26.0063 = rhs -26.0063) [ok] [P-weakderiv] jump function: int f phi' stays at the jump height (-1) as the bump narrows (-1.000, -1.000, -1.000) — no integrable weak derivative
For a 2-layer ReLU network on R^2 the supremum of the weak gradient norm and the largest difference quotient |f(x) - f(x')| / ||x - x'|| agree to three decimals: a bounded weak gradient IS the Lipschitz constant.
rng = np.random.default_rng(42)
W1 = rng.normal(size=(6, 2))
b1 = rng.normal(size=6)
w2 = rng.normal(size=6)
b2 = 0.3
def f_net(X): # 2-layer ReLU network on R^2
return np.maximum(X @ W1.T + b1, 0.0) @ w2 + b2
X = rng.normal(size=(4000, 2)) * 2.0
i = rng.integers(0, len(X), 20000)
j = rng.integers(0, len(X), 20000)
sep = np.linalg.norm(X[i] - X[j], axis=1)
ok = sep > 1e-9
lip_emp = np.max(np.abs(f_net(X[i][ok]) - f_net(X[j][ok])) / sep[ok])
h = 1e-6 # weak gradient, evaluated off the kinks
grad = np.stack([(f_net(X + h * np.eye(2)[k]) - f_net(X - h * np.eye(2)[k]))
/ (2 * h) for k in range(2)], axis=1)
sup_grad = np.linalg.norm(grad, axis=1).max()
check(f"[P-lipschitz] sup ||weak gradient|| ({sup_grad:.4f}) equals the "
f"largest difference quotient ({lip_emp:.4f})",
abs(sup_grad - lip_emp) < 1e-3)
[ok] [P-lipschitz] sup ||weak gradient|| (2.1249) equals the largest difference quotient (2.1247)
The product of the spectral norms of the weight matrices bounds that constant (valid but loose: 4.94 against 2.12 here).
L_spec = np.linalg.norm(W1, 2) * np.linalg.norm(w2)
check(f"[P-spectral] product of spectral norms bounds it: {lip_emp:.4f} "
f"<= {L_spec:.4f}", lip_emp <= L_spec + 1e-9)
[ok] [P-spectral] product of spectral norms bounds it: 2.1247 <= 4.9447
If |f(x)| = gamma and f is L-Lipschitz, no perturbation of norm < gamma / L changes sign(f(x)): 0 sign flips in 60000 perturbations drawn inside the certified radius. For a linear hypothesis the weak gradient is the constant w, so gamma / L is the distance of the feature vector from the decision boundary (the SVM statement).
flips, tested = 0, 0
for x in X[:300]:
gamma = abs(f_net(x[None])[0])
radius = gamma / lip_emp # certified radius
D = rng.normal(size=(200, 2))
D /= np.linalg.norm(D, axis=1, keepdims=True)
D *= rng.uniform(0.0, 0.999 * radius, size=(200, 1))
s0 = np.sign(f_net(x[None])[0])
flips += int(np.sum(np.sign(f_net(x + D)) != s0))
tested += len(D)
check(f"[P-certified] no sign flip inside the radius gamma/L "
f"({flips} flips in {tested} perturbations)", flips == 0)
# linear hypothesis: the weak gradient is the constant w, so gamma / L is
# the distance of the feature vector from the decision boundary
w_lin = np.array([1.5, -2.0])
b_lin = 0.4
x_test = np.array([0.7, 1.1])
gamma_lin = abs(w_lin @ x_test + b_lin)
dist = gamma_lin / np.linalg.norm(w_lin)
check(f"[P-certified] linear case: gamma/L ({gamma_lin / np.linalg.norm(w_lin):.4f}) "
f"is the distance to the decision boundary ({dist:.4f})",
abs(gamma_lin / np.linalg.norm(w_lin) - dist) < 1e-12)
[ok] [P-certified] no sign flip inside the radius gamma/L (0 flips in 60000 perturbations) [ok] [P-certified] linear case: gamma/L (0.3000) is the distance to the decision boundary (0.3000)
A decision stump is piecewise constant with a jump, so it has no bounded weak gradient: a perturbation of norm 2e-9 flips its prediction, and no certified radius exists at any margin.
stump = lambda X: np.where(X[:, 0] > 0.0, 1.0, -1.0)
eps = 1e-9
left = stump(np.array([[-eps, 0.5]]))[0]
right = stump(np.array([[+eps, 0.5]]))[0]
check(f"[P-jump] stump flips from {left:+.0f} to {right:+.0f} across a "
f"perturbation of norm {2 * eps:.0e} — no certified radius",
left != right)
[ok] [P-jump] stump flips from -1 to +1 across a perturbation of norm 2e-09 — no certified radius
Sharpening a transition of width eps: the squared H^1 seminorm grows like 4 / (3 eps) while the total variation stays at the height of the limiting jump — why the L^1 gradient penalty tolerates jumps and the H^1 penalty does not.
# Why the L^1 gradient penalty (total variation) tolerates jumps while the
# squared H^1 seminorm does not: sharpen a transition of width eps and watch
# the two energies. For f(x) = tanh(x / eps) the total variation stays at
# the height 2 of the limiting jump, while int |f'|^2 grows like 4 / (3 eps).
xs = np.linspace(-3.0, 3.0, 600001)
h_x = xs[1] - xs[0]
h1, tv = [], []
for eps in (0.2, 0.05, 0.0125):
df = (1.0 / eps) / np.cosh(xs / eps) ** 2 # derivative of tanh(x/eps)
h1.append(np.trapezoid(df ** 2, xs)) # squared H^1 seminorm
tv.append(np.trapezoid(np.abs(df), xs)) # total variation
check("[P-tvh1] H^1 energy grows without bound as the transition sharpens "
f"({', '.join(f'{v:.1f}' for v in h1)})",
all(b > 3.5 * a for a, b in zip(h1, h1[1:])))
check("[P-tvh1] total variation stays at the height of the jump "
f"({', '.join(f'{v:.4f}' for v in tv)})",
all(abs(v - 2.0) < 1e-3 for v in tv))
[ok] [P-tvh1] H^1 energy grows without bound as the transition sharpens (6.7, 26.7, 106.7) [ok] [P-tvh1] total variation stays at the height of the jump (2.0000, 2.0000, 2.0000)
Graph counterpart: for a graph signal sampled from a smooth function the edge-difference energy (the squared graph gradient summed over edges, i.e. the Laplacian quadratic form used by GTVMin) is small, while a signal with one jump makes it large — the discrete analogue of a bounded weak derivative. """ import numpy as np # the only numerical dependency report = [] # collects (check name, pass/fail) pairs def check(name, ok): # records and prints one verification report.append((name, bool(ok))) print(f" [{'ok' if ok else 'FAIL'}] {name}") # ------------------------------------------------------------ [P-weakderiv] # f(x) = max(0, x) has no derivative at x = 0, but the step g = 1_{x>0} # satisfies the definition of the weak derivative: for every smooth phi # vanishing near the boundary, int f phi' = - int g phi. grid = np.linspace(-2.0, 2.0, 400001) dx = grid[1] - grid[0] f_kink = np.maximum(grid, 0.0) # continuous, kinked g_step = (grid > 0).astype(float) # candidate weak derivative f_jump = g_step.copy() # discontinuous at 0 # test functions phi(x) = (4 - x^2)^3 * x^p vanish (with derivative) at +-2 for p in (0, 1, 2): phi = (4.0 - grid ** 2) ** 3 * grid ** p dphi = np.gradient(phi, dx) lhs = np.trapezoid(f_kink * dphi, grid) # int f phi' rhs = -np.trapezoid(g_step * phi, grid) # - int g phi check(f"[P-weakderiv] integration-by-parts identity holds for phi_{p} " f"(lhs {lhs:+.4f} = rhs {rhs:+.4f})", abs(lhs - rhs) < 1e-3) # the jump function admits no locally integrable weak derivative: the # identity would force int f_jump phi' = -int g phi for all phi, but the # left-hand side carries the value -phi(0) of the jump, which no integrable # g can reproduce (its contribution shrinks with the width of phi's support). defects = [] for eps in (0.4, 0.2, 0.1): # bumps of shrinking width bump = np.where(np.abs(grid) < eps, np.exp(-1.0 / np.maximum(1e-12, 1 - (grid / eps) ** 2)), 0.0) bump /= bump.max() dbump = np.gradient(bump, dx) defects.append(np.trapezoid(f_jump * dbump, grid)) # -> -phi(0) = -1 check("[P-weakderiv] jump function: int f phi' stays at the jump height " f"(-1) as the bump narrows ({', '.join(f'{d:+.3f}' for d in defects)}) " "— no integrable weak derivative", all(abs(d + 1.0) < 0.02 for d in defects)) # ------------------------------------------------------------- [P-lipschitz] rng = np.random.default_rng(42) W1 = rng.normal(size=(6, 2)) b1 = rng.normal(size=6) w2 = rng.normal(size=6) b2 = 0.3 def f_net(X): # 2-layer ReLU network on R^2 return np.maximum(X @ W1.T + b1, 0.0) @ w2 + b2
n = 60
nodes = np.linspace(0.0, 1.0, n)
edges = [(k, k + 1) for k in range(n - 1)] # path graph, unit weights
smooth = np.sin(2 * np.pi * nodes) # samples of a smooth function
jumpy = np.where(nodes < 0.5, -1.0, 1.0) # one jump of height 2
energy = lambda s: sum((s[a] - s[b]) ** 2 for a, b in edges)
check(f"[P-graph] edge-difference energy: smooth signal {energy(smooth):.4f} "
f"<< signal with one jump {energy(jumpy):.4f}",
energy(smooth) < 0.1 * energy(jumpy))
n_ok = sum(ok for _, ok in report)
print(f"\n{n_ok}/{len(report)} checks pass")
if n_ok != len(report):
raise SystemExit(1)
[ok] [P-graph] edge-difference energy: smooth signal 0.3342 << signal with one jump 4.0000 12/12 checks pass