import json, math, os import numpy as np SEED = 1195 rng = np.random.default_rng(SEED) def delayed_quadratic(a, delay, steps=5000, x0=1.0): # x[t+1] = x[t] - a*x[t-delay], with zero-padded history x = np.ones(steps + delay + 1) * x0 for t in range(steps): x[t+delay+1] = x[t+delay] - a*x[t] tail = np.max(np.abs(x[-500:])) return tail, np.max(np.abs(x)) def theoretical_delay_limit(d): # Exact Schur stability boundary for z[t+1]=z[t]-a z[t-d]. return 2.0 * math.sin(math.pi / (4*d + 2)) def numerical_boundary(d): lo, hi = 0.0, 3.0 # binary search using a long trajectory; margin avoids classifying boundary oscillations as stable for _ in range(38): mid = (lo + hi) / 2 tail, peak = delayed_quadratic(mid, d) if np.isfinite(tail) and tail < 0.8 and peak < 20: lo = mid else: hi = mid return lo def lipschitz_check(): # F(x)=.5*lambda*x^2 has L=lambda, so the stated inequality should be exact. lam = 3.7 xs = np.linspace(-2.0, 2.0, 101) ys = np.linspace(-2.0, 2.0, 101) ratios = [] for x in xs: for y in ys: if abs(x-y) > 1e-10: ratios.append((abs(lam*x-lam*y)**2)/(lam**2*(x-y)**2)) return float(np.max(ratios)), float(np.min(ratios)) def delay_error_scaling(): # A slowly moving parameter trajectory x_t=v*t makes displacement exactly proportional to d. # For a quadratic, ||grad(x_t)-grad(x_{t-d})||^2 = L^2 ||x_t-x_{t-d}||^2, # hence absolute error contribution is predicted to scale as d^2. lam, v, t = 2.0, 0.013, 1000 ds = np.arange(1, 9) errs = np.array([(lam*v*d)**2 for d in ds]) slope = np.polyfit(np.log(ds), np.log(errs), 1)[0] return ds.tolist(), errs.tolist(), float(slope), float(errs[-1]/errs[0]) def async_stream(seed=SEED, n_parts=8, workers=8, c=4, steps=180, eta=0.16): """Small event-driven coded stream on separable logistic regression shards. Each worker has a replica assignment and random compute delays. The coordinator uses the newest available gradient per shard once every shard is at most c versions old. """ rg = np.random.default_rng(seed) dim, shard_size = 10, 48 true_w = rg.normal(size=dim) shards = [] for i in range(n_parts): X = rg.normal(size=(shard_size, dim)) y = (X @ true_w + .35*rg.normal(size=shard_size) > 0).astype(float) shards.append((X, y)) # two replicas, assigned round-robin with distinct worker speeds assignments = [[i for i in range(n_parts) if i % workers == w or i % workers == (w-1) % workers] for w in range(workers)] speed = np.exp(rg.normal(np.log(1.0), .38, size=workers)) # queue entries: (completion time, worker, partition, version, gradient) queue = [] table = [None] * n_parts beta_hist = [np.zeros(dim)] clock, t, updates = 0.0, 0, 0 ages, wall = [], [] def grad(beta, i): X, y = shards[i] z = np.clip(X @ beta, -30, 30) p = 1/(1+np.exp(-z)) return X.T @ (p-y) / len(y) def schedule(w, part, version, beta): # occasional stragglers, but no unrealistic seconds are needed in normalized time delay = speed[w] * float(np.exp(rg.normal(0, .35))) if rg.random() < .07: delay *= 4.5 queue.append((clock+delay, w, part, version, grad(beta, part))) # initial replicas on version 0 for w in range(workers): for i in assignments[w]: schedule(w, i, 0, beta_hist[0]) while t < steps and queue: queue.sort(key=lambda z: z[0]) clock = queue[0][0] ready = [q for q in queue if q[0] <= clock + 1e-12] queue = [q for q in queue if q[0] > clock + 1e-12] for _, w, i, v, g in ready: if table[i] is None or v >= table[i][0]: table[i] = (v, g) if all(x is not None and 0 <= t-x[0] < c for x in table): used = [x[0] for x in table] gmean = np.mean([x[1] for x in table], axis=0) beta = beta_hist[-1] - eta*gmean beta_hist.append(beta); ages.extend([t-v for v in used]); wall.append(clock) t += 1; updates += 1 # each worker begins one new partition evaluation on the new snapshot for w in range(workers): schedule(w, (t+w) % n_parts, t, beta) # objective on fresh data Xv = rg.normal(size=(800, dim)); yv=(Xv@true_w+.35*rg.normal(size=800)>0).astype(float) z=np.clip(Xv@beta_hist[-1],-30,30) loss=float(np.mean(np.logaddexp(0,z)-yv*z)) return dict(loss=loss, updates=updates, wall=float(wall[-1]) if wall else float('inf'), mean_age=float(np.mean(ages)) if ages else float('inf'), max_age=int(max(ages)) if ages else -1) def sync_stream(seed=SEED, steps=180, eta=.16): # Same stochastic problem and delay distribution, but one round waits for all workers. # Approximate synchronous wall time by max of eight independent worker times per round. rg=np.random.default_rng(seed); dim=10; n=8; m=48 tw=rg.normal(size=dim); shards=[] for i in range(n): X=rg.normal(size=(m,dim)); y=(X@tw+.35*rg.normal(size=m)>0).astype(float); shards.append((X,y)) beta=np.zeros(dim); total=0. def g(i): X,y=shards[i]; z=np.clip(X@beta,-30,30); p=1/(1+np.exp(-z)); return X.T@(p-y)/m for t in range(steps): gs=[g(i) for i in range(n)]; beta-=eta*np.mean(gs,axis=0) delays=rg.lognormal(0,.38,size=8); delays[rg.random(8)<.07]*=4.5; total+=float(np.max(delays)) X=rg.normal(size=(800,dim)); y=(X@tw+.35*rg.normal(size=800)>0).astype(float); z=np.clip(X@beta,-30,30) return dict(loss=float(np.mean(np.logaddexp(0,z)-y*z)), updates=steps, wall=total, mean_age=0., max_age=0) def main(): stability=[] for d in [0,1,2,3,5,8]: pred=theoretical_delay_limit(d); obs=numerical_boundary(d) stability.append(dict(delay=d, predicted=pred, observed=obs, relative_error=abs(obs-pred)/pred)) lip=lipschitz_check(); ds, errs, slope, growth=delay_error_scaling() base=sync_stream(); idea={str(c):async_stream(c=c) for c in [2,4,8]} out={'seed':SEED, 'stability_boundary':stability, 'lipschitz_bound_max_ratio':lip[0], 'lipschitz_bound_min_ratio':lip[1], 'delay_error_delays':ds, 'delay_error_values':errs, 'delay_error_loglog_slope':slope, 'delay_error_growth_d1_to_d8':growth, 'baseline_sync':base, 'pipelined':idea} with open('results.json','w') as f: json.dump(out,f,indent=2) print(json.dumps(out,indent=2)) if __name__=='__main__': main()