import json import math from pathlib import Path import numpy as np def clip_norm(x, tau): n = np.linalg.norm(x) return x if n <= tau else x * (tau / n) def k_formula(alpha, beta, p): if alpha <= p * beta: return beta return (alpha / p) * (alpha * (p - 1) / (p * (alpha - beta))) ** (p - 1) def r_star_formula(alpha, beta, p): if alpha <= p * beta: return 1.0 return p * (alpha - beta) / (alpha * (p - 1)) def normalized_cost(r, alpha, beta): return beta * r * r if r <= 1 else alpha * (r - 1) + beta def verify_math(): rows = [] # Dense radial sweep is the exact one-dimensional reduction in the theorem. rs = np.unique(np.r_[np.linspace(0, 1, 20001), np.geomspace(1.000001, 1e4, 50000)]) for p in (1.25, 1.5, 2.0): for alpha, beta in ((0.8, 1.0), (2.0, 1.0), (1.0, 1.0), (3.0, .5), (0.25, .5)): vals = np.array([normalized_cost(r, alpha, beta) / max(r ** p, 1e-300) for r in rs]) j = int(np.argmax(vals)) pred_r = r_star_formula(alpha, beta, p) # At p=2 in the energy phase the theorem predicts a whole tie interval r<=1. tie_valid = (p == 2.0 and alpha <= p*beta and rs[j] <= 1.0 + 1e-8) radius_error = 0.0 if tie_valid else abs(float(rs[j])-pred_r) rows.append({"p":p, "alpha":alpha, "beta":beta, "predicted_phase":"energy" if alpha <= p*beta else "bias", "observed_phase":"energy" if rs[j] <= 1.0001 else "bias", "boundary_margin":alpha-p*beta, "predicted_r":pred_r, "observed_r":float(rs[j]), "radius_error":radius_error, "p2_tie_valid":bool(tie_valid), "predicted_K":k_formula(alpha,beta,p), "observed_K":float(vals[j])}) # Explicit transition sweep: the maximizer should remain at r=1 through alpha=p beta, # then move continuously outside with the exact r_star formula. transition = [] p0, beta0 = 1.5, 1.0 for alpha0 in (1.40, 1.49, 1.50, 1.51, 1.60, 2.00): vv = np.array([normalized_cost(r, alpha0, beta0) / max(r**p0, 1e-300) for r in rs]) jj = int(np.argmax(vv)) transition.append({"alpha":alpha0, "predicted_r":r_star_formula(alpha0,beta0,p0), "observed_r":float(rs[jj]), "predicted_phase":"energy" if alpha0 <= p0*beta0 else "bias", "observed_K":float(vv[jj]), "predicted_K":k_formula(alpha0,beta0,p0)}) # Pure residual: max (r-1)_+ / r^p = c_p at p/(p-1), and scaling is tau-independent. residual_rows = [] for p in (1.25, 1.5, 2.0): rs2 = np.geomspace(1.0000001, 1e5, 100000) q = (rs2 - 1) / rs2**p j = int(np.argmax(q)) cp = (p-1)**(p-1) / p**p residual_rows.append({"p":p, "predicted_r":p/(p-1), "observed_r":float(rs2[j]), "predicted_c":cp, "observed_c":float(q[j])}) # Direct vector scaling check over several tau and dimensions. scaling = [] rng = np.random.default_rng(7) for p in (1.5, 2.0): cp = (p-1)**(p-1) / p**p for tau in (.2, 1., 5.): r = p/(p-1) x = np.array([r*tau, 0.0, 0.0]) resid = np.linalg.norm(x-clip_norm(x,tau)) ratio = resid / (tau**(1-p) * np.linalg.norm(x)**p) scaling.append({"p":p,"tau":tau,"ratio":float(ratio),"predicted":cp}) return rows, residual_rows, scaling, transition class PhaseAwareController: def __init__(self, p=2.0, alpha=1.0, beta=1.0, tau=1.0, eta=.035): self.p, self.alpha, self.beta = p, alpha, beta self.tau, self.eta = tau, eta self.phase = "energy" if alpha <= p*beta else "bias" self.history = [] def step(self, gradients): norms = np.linalg.norm(gradients, axis=1) clipped = gradients.copy() factors = np.minimum(1., self.tau / np.maximum(norms, 1e-12)) clipped *= factors[:, None] residual = np.mean(np.linalg.norm(gradients-clipped, axis=1)) energy = np.mean(np.sum(clipped*clipped, axis=1)) / self.tau # Targets are deliberately fixed relative to the initial trust region. target = .22 if self.phase == "energy" else .10 signal = energy if self.phase == "energy" else residual delta = np.clip(self.eta * (signal-target), -.08, .08) self.tau *= math.exp(delta) self.tau = float(np.clip(self.tau, .08, 8.0)) self.history.append((self.tau, residual, energy, np.mean(factors < 1))) return clipped def percentile_clip(g, tau_state, percentile=90): norms = np.linalg.norm(g, axis=1) threshold = np.percentile(norms, percentile) tau = min(tau_state, threshold) factors = np.minimum(1., tau / np.maximum(norms, 1e-12)) return g * factors[:,None], tau def mini_experiment(seed=11, steps=700, batch=64, dim=8): rng = np.random.default_rng(seed) methods = ["fixed", "percentile", "phase_aware"] results = {} for method in methods: w = np.ones(dim) * 3.0 tau, taus = 1.0, [] ctl = PhaseAwareController(alpha=1., beta=1., tau=1.) if method == "phase_aware" else None losses, updates, clips, residuals, energies = [], [], [], [], [] for t in range(steps): # Quadratic objective with occasional Pareto-like, direction-random outlier batches. g = w[None,:] + rng.normal(0, .18, (batch,dim)) if t % 70 == 35: g += rng.normal(size=(batch,dim)) * 12.0 if method == "fixed": factors = np.minimum(1., tau / np.maximum(np.linalg.norm(g,axis=1),1e-12)) cg = g * factors[:,None] elif method == "percentile": cg, tau = percentile_clip(g, tau) factors = np.minimum(1., tau / np.maximum(np.linalg.norm(g,axis=1),1e-12)) else: old_tau = ctl.tau cg = ctl.step(g) tau = ctl.tau factors = np.minimum(1., old_tau / np.maximum(np.linalg.norm(g,axis=1),1e-12)) update = .045 * np.mean(cg, axis=0) w -= update losses.append(.5*np.dot(w,w)); updates.append(np.linalg.norm(update)) clips.append(np.mean(factors < 1)); residuals.append(np.mean(np.linalg.norm(g-cg,axis=1))) energies.append(np.mean(np.sum(cg*cg,axis=1))/max(tau,1e-9)); taus.append(tau) # recovery is steps needed after final outlier to get below .1 loss, or sentinel tail = np.array(losses[596:]) hit = np.where(tail < .1)[0] results[method] = {"final_loss":float(losses[-1]), "median_update":float(np.median(updates)), "update_std":float(np.std(updates)), "clip_fraction":float(np.mean(clips)), "mean_residual":float(np.mean(residuals)), "mean_energy":float(np.mean(energies)), "recovery_steps_after_last_outlier":int(hit[0]) if len(hit) else None, "final_tau":float(taus[-1])} return results def main(): math_rows, residual_rows, scaling, transition = verify_math() mini = mini_experiment() out = {"math_sweep":math_rows,"transition_sweep":transition,"residual_sweep":residual_rows,"scaling_check":scaling,"mini_experiment":mini} Path("results.json").write_text(json.dumps(out, indent=2)) print(json.dumps(out, indent=2)) if __name__ == "__main__": main()