""" gradient.py — numerical companion to the glossary entry 'gradient'. Purpose ------- Verifies, for the quadratic f(w) = (1/2) w^T Q w with Q = [[2, 0.6], [0.6, 1]], the defining and geometric properties of the gradient stated in the entry, and generates the level-set data for the entry's figure. Self-contained (numpy/matplotlib only), deterministic. Blocks ------ [B-fd] The analytic gradient grad f(w') = Q w' matches central finite differences at w' = (1.4, 1.1) to 1e-8. [B-taylor] Local linear approximation: the error |f(w) - f(w') - grad^T (w - w')| / ||w - w'|| vanishes as w -> w' (ratio decreases by ~10x per 10x step shrink). [B-steepest] Among 720 unit directions d, the directional derivative grad^T d is maximized by d = grad/||grad|| (within 0.5 deg). [B-orth] The directional derivative along the level-set tangent at w' is zero (up to 1e-12): the gradient is orthogonal to the level set. [B-partials] Partial derivatives alone do not guarantee a gradient for a non-convex function: g(w) = w1 w2^2 / (w1^2 + w2^4) has both partials equal to 0 at the origin, yet the linear- approximation error ratio DIVERGES along w = (t^2, t) — no gradient exists there (the entry's convexity assumption is not dispensable). [B-hilbert] Inner-product dependence of the gradient: w.r.t. the inner product _M = u^T M v (an SPD M), the gradient of f at w' is M^{-1} Q w'. It represents the directional derivative (_M = (Q w')^T d for all d) and is the steepest-ascent direction among 720 directions of unit M-norm. [B-min] At the minimizer w-hat = 0 of f, the gradient vanishes. [B-erm] ML relevance: for a synthetic training set (m = 30) and a linear regression, the analytic ERM gradient -(2/m) sum (y - w^T x) x matches finite differences, vanishes at the least-squares solution, and one GD step decreases the objective. [B-backprop] For a one-hidden-layer network with tanh activation, the gradient of the ERM objective computed by backpropagation (chain rule, layer by layer) matches finite differences on all weights. Outputs ------- gradient_levelsets.csv : polylines of three level sets f(w) = const. through and inside w' (columns x,y; nan rows separate the level sets) for the entry's pgfplots figure. gradient.png : matplotlib preview of the figure (checking only). The figure's gradient arrow at w' is the unit vector of grad f(w') = Q w' = (3.46, 1.94); its coordinates are printed below and used verbatim in the entry's TikZ code. """ import numpy as np # the only numerical dependency import matplotlib # imported before pyplot to set the backend matplotlib.use("Agg") # non-interactive backend: no display needed import matplotlib.pyplot as plt # plotting API for the preview figure report = [] # collects (check name, pass/fail) pairs def check(name, ok): # records and prints one verification report.append((name, bool(ok))) # store the verdict print(f" [{'ok' if ok else 'FAIL'}] {name}") # echo it immediately Q = np.array([[2.0, 0.6], [0.6, 1.0]]) # SPD matrix defining f(w) = w^T Q w / 2 def f(w): # the quadratic objective f(w) = (1/2) w^T Q w # einsum contracts w_i Q_ij w_j over the LAST axis, so w may be a # single point (shape (2,)) or a whole grid (shape (..., 2)) return 0.5 * np.einsum("...i,ij,...j->...", w, Q, w) def grad(w): # analytic gradient of f: grad f(w) = Q w return Q @ w # matrix-vector product wp = np.array([1.4, 1.1]) # the point w' marked in the entry's figure g = grad(wp) # gradient at w' (drawn as the solid arrow) # ------------------------------------------------------------------ [B-fd] eps = 1e-6 # half-width of the finite-difference stencil # central difference (f(w'+eps e_j) - f(w'-eps e_j)) / (2 eps) per coordinate g_fd = np.array([(f(wp + eps * e) - f(wp - eps * e)) / (2 * eps) for e in np.eye(2)]) check("[B-fd] analytic gradient matches finite differences", np.linalg.norm(g - g_fd) < 1e-8) # entries agree to 1e-8 # -------------------------------------------------------------- [B-taylor] rng = np.random.default_rng(0) # fixed seed: deterministic direction d = rng.standard_normal(2) # a random approach direction d /= np.linalg.norm(d) # normalized so ||w - w'|| = h below ratios = [] # error ratios for shrinking displacements h for h in (1e-1, 1e-2, 1e-3): # three decades of displacement h w = wp + h * d # approach point w = w' + h d # the entry's defining ratio |f(w) - f(w') - g^T (w - w')| / ||w - w'|| ratios.append(abs(f(w) - f(wp) - g @ (w - wp)) / h) check("[B-taylor] linear-approximation error ratio vanishes", # for a quadratic the ratio is O(h): each 10x shrink of h must shrink # the ratio ~10x (the 1.5 slack absorbs rounding) ratios[0] > 9 * ratios[1] > 81 * ratios[2] / 1.5) # ------------------------------------------------------------ [B-steepest] angles = np.linspace(0.0, 2.0 * np.pi, 720, endpoint=False) # 0.5-degree grid dirs = np.c_[np.cos(angles), np.sin(angles)] # 720 unit direction vectors best = dirs[np.argmax(dirs @ g)] # direction maximizing the directional derivative g^T d # angle (in degrees) between that maximizer and the normalized gradient ang_err = np.degrees(np.arccos(np.clip(best @ (g / np.linalg.norm(g)), -1.0, 1.0))) check("[B-steepest] gradient direction maximizes directional derivative", ang_err < 0.5) # agreement within the grid resolution # ---------------------------------------------------------------- [B-orth] t = np.array([-g[1], g[0]]) / np.linalg.norm(g) # level-set tangent: rotate g by 90 deg check("[B-orth] zero directional derivative along the level set", abs(g @ t) < 1e-12) # g^T t = 0: gradient orthogonal to the level set # ------------------------------------------------------------ [B-partials] # g(w) = w1 w2^2 / (w1^2 + w2^4): both partials exist (= 0) at the origin, # but g is not differentiable there (g is not convex). def g_cex(w1, w2): # the counterexample function den = w1 ** 2 + w2 ** 4 # denominator, zero only at the origin # define g(0,0) = 0 (its limit along both axes); elsewhere the formula return np.where(den == 0.0, 0.0, w1 * w2 ** 2 / den) eps = 1e-6 # finite-difference half-width (as in [B-fd]) px = (g_cex(eps, 0.0) - g_cex(-eps, 0.0)) / (2 * eps) # partial w.r.t. w1 at 0 py = (g_cex(0.0, eps) - g_cex(0.0, -eps)) / (2 * eps) # partial w.r.t. w2 at 0 # along w = (t^2, t), g = 1/2 while ||w|| -> 0: the error ratio # |g(w) - 0 - 0| / ||w|| blows up, so no gradient exists at 0. ratios_cex = [abs(g_cex(t ** 2, t)) / np.hypot(t ** 2, t) for t in (1e-1, 1e-2, 1e-3)] check("[B-partials] partials exist at 0 but the error ratio diverges", # partials vanish exactly, and the ratio GROWS ~10x per decade abs(px) < 1e-12 and abs(py) < 1e-12 and ratios_cex[2] > 9 * ratios_cex[1] > 81 * ratios_cex[0] / 1.5) # ------------------------------------------------------------- [B-hilbert] M = np.array([[2.0, 0.5], [0.5, 1.0]]) # SPD inner-product matrix gM = np.linalg.solve(M, g) # gradient w.r.t. <.,.>_M: M^{-1} Q w' rng_h = np.random.default_rng(1) # fixed seed: deterministic test directions D = rng_h.standard_normal((100, 2)) # 100 random directions d check("[B-hilbert] _M equals the directional derivative g^T d", # representation property: d^T M g_M = d^T g for every d np.allclose(D @ M @ gM, D @ g, atol=1e-12)) # rescale the 720 directions from [B-steepest] to unit M-norm ||d||_M = 1 dirs_M = dirs / np.sqrt(np.einsum("ij,jk,ik->i", dirs, M, dirs))[:, None] best_M = dirs_M[np.argmax(dirs_M @ g)] # maximizer of the directional derivative g^T d check("[B-hilbert] steepest ascent w.r.t. the M-norm is along g_M", # angle between that maximizer and g_M is below the grid resolution np.degrees(np.arccos(np.clip( best_M @ gM / (np.linalg.norm(best_M) * np.linalg.norm(gM)), -1.0, 1.0))) < 0.5) # ----------------------------------------------------------------- [B-min] check("[B-min] gradient vanishes at the minimizer w-hat = 0", np.linalg.norm(grad(np.zeros(2))) == 0.0) # grad f(0) = Q 0 = 0 exactly # ----------------------------------------------------------------- [B-erm] # linear regression: f(w) = (1/m) sum (y - w^T x)^2 with # gradient -(2/m) sum (y - w^T x) x, as displayed in the entry. rng_e = np.random.default_rng(2) # fixed seed: reproducible training set m = 30 # size of the training set X = rng_e.standard_normal((m, 2)) # rows are the feature vectors x^(r) y = X @ np.array([1.0, -0.5]) + 0.1 * rng_e.standard_normal(m) # noisy labels def f_erm(w): # the ERM objective: average squared error return np.mean((y - X @ w) ** 2) def grad_erm(w): # its analytic gradient: -(2/m) X^T (y - X w) return -(2.0 / m) * X.T @ (y - X @ w) w0 = np.array([0.4, 0.8]) # an arbitrary (non-optimal) parameter vector # central finite differences of f_erm at w0, coordinate by coordinate g_erm_fd = np.array([(f_erm(w0 + eps * e) - f_erm(w0 - eps * e)) / (2 * eps) for e in np.eye(2)]) check("[B-erm] analytic ERM gradient matches finite differences", np.linalg.norm(grad_erm(w0) - g_erm_fd) < 1e-8) # agree to 1e-8 w_hat = np.linalg.lstsq(X, y, rcond=None)[0] # least-squares solution (ERM minimizer) check("[B-erm] ERM gradient vanishes at the least-squares solution", np.linalg.norm(grad_erm(w_hat)) < 1e-10) # zero-gradient condition at w-hat check("[B-erm] one GD step decreases the ERM objective", f_erm(w0 - 0.1 * grad_erm(w0)) < f_erm(w0)) # step along -grad with lrate 0.1 # ------------------------------------------------------------ [B-backprop] # one-hidden-layer network h(x) = w2^T tanh(W1 x); the gradient of the ERM # objective via the chain rule (backpropagation) vs finite differences. rng_b = np.random.default_rng(3) # fixed seed: reproducible weights W1 = rng_b.standard_normal((3, 2)) * 0.5 # hidden-layer weights (3 units, 2 inputs) w2 = rng_b.standard_normal(3) * 0.5 # output-layer weights def f_net(W1_, w2_): # ERM objective of the network on (X, y) return np.mean((y - np.tanh(X @ W1_.T) @ w2_) ** 2) A = np.tanh(X @ W1.T) # m x 3 activations (forward pass) res = y - A @ w2 # m residuals y - h(x) dpred = -(2.0 / m) * res # dL/dprediction, per sample gb_w2 = A.T @ dpred # backprop: output layer (dL/dw2) # hidden layer: chain rule through tanh' = 1 - A^2, then through W1 x gb_W1 = ((dpred[:, None] * (1.0 - A ** 2) * w2).T @ X) # hidden layer theta_fd = [] # finite-difference gradient, same parameter order for i in range(3): # loop over hidden units ... for j in range(2): # ... and input coordinates P = np.zeros_like(W1); P[i, j] = eps # perturbation of one W1 entry # central difference w.r.t. that single entry theta_fd.append((f_net(W1 + P, w2) - f_net(W1 - P, w2)) / (2 * eps)) for i in range(3): # loop over the output-layer weights p = np.zeros(3); p[i] = eps # perturbation of one w2 entry theta_fd.append((f_net(W1, w2 + p) - f_net(W1, w2 - p)) / (2 * eps)) check("[B-backprop] backprop gradient matches finite differences", # flatten (W1, w2) backprop gradients and compare to the stencil np.linalg.norm(np.r_[gb_W1.ravel(), gb_w2] - theta_fd) < 1e-8) # --------------------------------------------------------------- level sets G = 500 # grid resolution per axis # the grid must contain the outermost level set f(w) = f(w') entirely # (max |w_2| on it is sqrt(2 f(w') (Q^{-1})_{22}) ~ 2.92): a too-small grid # clips the contour and the closing segment below draws a straight chord. gx = np.linspace(-3.2, 3.6, G) # horizontal grid coordinates gy = np.linspace(-3.2, 3.2, G) # vertical grid coordinates GX, GY = np.meshgrid(gx, gy) # full 2-D evaluation grid F = f(np.stack([GX, GY], axis=-1)) # f evaluated on the whole grid at once levels = [0.8, 2.0, float(f(wp))] # two inner levels plus the level through w' fig0, ax0 = plt.subplots() # throwaway figure: only the contour data is used cs = ax0.contour(GX, GY, F, levels=levels) # compute the three contour polylines plt.close(fig0) # discard the figure, keep the contour set with open("pythondemos/gradient_levelsets.csv", "w") as fh: # pgfplots data file fh.write("x,y\n") # column header expected by \addplot table first = True # tracks whether a nan separator is needed for segs in cs.allsegs: # one entry per contour level ... for s in segs: # ... each holding its polyline segments if not first: # between polylines ... fh.write("nan,nan\n") # ... a nan row breaks the pgfplots path first = False # subsequent polylines need the separator for p in s[::5]: # every 5th vertex suffices at print size fh.write(f"{p[0]:.4f},{p[1]:.4f}\n") # one vertex per row if np.allclose(s[0], s[-1]): # close closed loops only fh.write(f"{s[0][0]:.4f},{s[0][1]:.4f}\n") # repeat first vertex u = g / np.linalg.norm(g) * 1.2 # unit gradient arrow, scaled to length 1.2 # print the arrow coordinates that the entry's TikZ code uses verbatim print(f"\nw' = {wp} | grad f(w') = {g} | unit arrow (len 1.2) = " f"({u[0]:.2f}, {u[1]:.2f})") # -------------------------------------------------------------- preview fig, ax = plt.subplots(figsize=(4.8, 3.8)) # preview canvas (checking only) ax.contour(GX, GY, F, levels=levels, colors="k", linewidths=0.9) # level sets ax.plot(*wp, "ko", ms=4) # mark the point w' ax.annotate("", xy=wp + u, xytext=wp, # solid arrow: the gradient at w' arrowprops=dict(arrowstyle="-|>", lw=2)) ax.annotate("", xy=wp - u, xytext=wp, # dashed arrow: the negative gradient arrowprops=dict(arrowstyle="-|>", lw=1.2, linestyle="--")) ax.plot(0, 0, "k.", ms=6) # mark the minimizer w-hat = 0 ax.set_aspect("equal") # equal axis scaling (as in the TikZ figure) ax.set_axis_off() # the entry's figure has no axes either fig.tight_layout() # trim whitespace fig.savefig("pythondemos/gradient.png", dpi=110) # write the preview PDF n_ok = sum(ok for _, ok in report) # count the passed checks print(f"{n_ok}/{len(report)} checks pass") # summary line print("wrote pythondemos/gradient_levelsets.csv, pythondemos/gradient.png") if n_ok != len(report): # any failed check ... raise SystemExit(1) # ... makes the script exit non-zero