"""Heavy-tail path-adaptive optimizer pool: small NumPy reference implementation.""" import numpy as np class RestartedAdaGrad: def __init__(self, dim, horizon, alpha=0.25, epsilon=1e-8, x0=None, radius=10.0): self.dim, self.horizon = dim, horizon self.alpha, self.epsilon, self.radius = alpha, epsilon, radius self.x = np.zeros(dim) if x0 is None else np.array(x0, dtype=float).copy() self.age = 0 self.acc = np.full(dim, epsilon, dtype=float) def step(self, grad): grad = np.asarray(grad, dtype=float) # B_k(t) is the current block; reset the diagonal accumulator at restart. if self.age >= self.horizon: self.age, self.acc = 0, np.full(self.dim, self.epsilon) self.acc += grad * grad self.x -= self.alpha * grad / np.sqrt(self.acc) self.x = np.clip(self.x, -self.radius, self.radius) self.age += 1 return self.x class HeavyTailPool: def __init__(self, dim=1, horizons=(1, 2, 4, 8, 16), alpha=.25, beta=2.0, epsilon=1e-8, x0=None, radius=10.): self.horizons = tuple(horizons) self.experts = [RestartedAdaGrad(dim, h, alpha, epsilon, x0, radius) for h in self.horizons] self.q = np.ones(len(horizons)) / len(horizons) self.v = np.zeros(len(horizons)) self.beta, self.epsilon = beta, epsilon def step(self, grad_fn, validation_loss_fn): # Each expert receives its own gradient, as in the proposed construction. for e in self.experts: e.step(grad_fn(e.x)) losses = np.array([validation_loss_fn(e.x) for e in self.experts]) mean_loss = float(self.q @ losses) z = losses - mean_loss self.v += z*z logits = np.log(self.q) - self.beta*z/np.sqrt(self.epsilon + self.v) logits -= logits.max() self.q = np.exp(logits); self.q /= self.q.sum() x = sum(q*e.x for q, e in zip(self.q, self.experts)) return x, losses, self.q.copy() def fixed_adagrad(x0, gradients, alpha=.25, epsilon=1e-8, radius=10.): x = np.array(x0, dtype=float).copy(); acc = np.full_like(x, epsilon) out=[] for g in gradients: acc += g*g; x -= alpha*g/np.sqrt(acc); x=np.clip(x,-radius,radius); out.append(x.copy()) return np.array(out) def run_shift(T=240, shift=80, seed=0, horizons=(1,2,4,8,16,32), alpha=.32, beta=3.): rng=np.random.default_rng(seed); target=np.zeros(T); target[shift:]=2.0 pool=HeavyTailPool(1,horizons,alpha,beta,x0=[0.],radius=4.) xs=[]; qs=[]; losses=[] for t in range(T): # convex quadratic observation with moderate heavy-tailed gradient noise def gf(x, t=t): return (x-target[t:t+1]) + .06*rng.standard_t(1.7, size=1) def lf(x, t=t): return .5*float((x[0]-target[t])**2) x,l,q=pool.step(gf,lf); xs.append(float(x[0])); qs.append(q); losses.append(float(.5*(x[0]-target[t])**2)) return np.array(xs),np.array(losses),np.array(qs),target def main(): # 1) Exact constant-gradient scaling: displacement is alpha sum 1/sqrt(eps+t g^2). # Prediction: for fixed g and large H, displacement ~ 2 alpha sqrt(H)/|g|. alpha=.7; g=2.; eps=1e-8 hs=np.array([16,32,64,128,256,512,1024]) observed=[]; predicted=[] for h in hs: a=eps+np.arange(1,h+1)*g*g observed.append(alpha*np.sum(g/np.sqrt(a))) predicted.append(2*alpha*np.sqrt(h)) # / sign(g), because |g| cancels in 1-D slope=np.polyfit(np.log(hs),np.log(observed),1)[0] print('PREDICTION displacement ~ sqrt(H):') print('H observed predicted ratio') for h,o,p in zip(hs,observed,predicted): print(h, f'{o:.6f}', f'{p:.6f}', f'{o/p:.4f}') print('observed_loglog_slope',slope,'predicted',.5) # 2) Meta mechanism: on fixed losses, log q_i/q_j follows sqrt(T) for a constant gap. # Construct two experts with losses [0, delta], initially equal. z=[-delta/2,+delta/2]. beta=2.; delta=.2; Ts=np.array([25,100,400,1600]); logodds=[]; pred=[] for T in Ts: v=np.zeros(2); lo=0. for _ in range(T): z=np.array([-delta/2,delta/2]); v += z*z lo += float((-beta*z[0]/np.sqrt(eps+v[0])) - (-beta*z[1]/np.sqrt(eps+v[1]))) logodds.append(lo) # Each expert has v_k=t*delta^2/4, so the two log-weight # increments sum asymptotically to 2*beta/sqrt(t). pred.append(4*beta*np.sqrt(T)) meta_slope=np.polyfit(np.log(Ts),np.log(logodds),1)[0] print('\nPREDICTION meta log-odds ~ gap*sqrt(T):') print('T observed predicted ratio') for t,o,p in zip(Ts,logodds,pred): print(t, f'{o:.6f}', f'{p:.6f}', f'{o/p:.4f}') print('observed_loglog_slope',meta_slope,'predicted',.5) # 3) Dynamic prediction: after a jump, best expert's recovery lag scales with H. # Recovery lag is the first time its squared error is below 0.25 after shift. print('\nPREDICTION dynamic recovery lag is O(H):') for h in [1,2,4,8,16,32]: x,loss,q,target=run_shift(T=180,shift=60,seed=11,horizons=(h,),alpha=.32,beta=3.) post=np.where(np.abs(x[60:]-2.)<.5)[0] lag=int(post[0]) if len(post) else 180-60 print(h,lag, 'lag_over_H',f'{lag/h:.3f}') # Mini comparison: same shifted quadratic, pool versus fixed AdaGrad, equal one gradient/expert seeds=range(10); pool_scores=[]; base_scores=[]; pool_recovery=[]; base_recovery=[] for seed in seeds: x,loss,q,target=run_shift(seed=seed) rng=np.random.default_rng(seed); gs=[] for t in range(240): gs.append(np.array([0.0])) # use deterministic baseline trajectory below # baseline gets the same target sequence and no noisy gradient, isolating adaptation behavior tg=np.zeros(240); tg[80:]=2.; xx=fixed_adagrad([0.], [np.array([0.])-tg[t:t+1] for t in range(240)], alpha=.32) bl=.5*(xx[:,0]-tg)**2 pool_scores.append(loss[80:].mean()); base_scores.append(bl[80:].mean()) pool_recovery.append(np.where(np.abs(x[80:]-2.)<.5)[0][0] if np.any(np.abs(x[80:]-2.)<.5) else 160) base_recovery.append(np.where(np.abs(xx[80:,0]-2.)<.5)[0][0] if np.any(np.abs(xx[80:,0]-2.)<.5) else 160) print('\nCOMPARISON (pool has noisy heavy-tail gradients; baseline is fixed AdaGrad):') print('pool post-shift loss',np.mean(pool_scores),'+/-',np.std(pool_scores)/np.sqrt(10)) print('baseline post-shift loss',np.mean(base_scores),'+/-',np.std(base_scores)/np.sqrt(10)) print('pool recovery',np.mean(pool_recovery),'baseline recovery',np.mean(base_recovery)) if __name__=='__main__': main()