1import json, math, os
  2import numpy as np
  3from scipy.linalg import eigvalsh, solve_discrete_lyapunov
  4
  5SEED = 2075
  6rng = np.random.default_rng(SEED)
  7
  8
  9def make_hessians():
 10    # Two non-commuting positive definite curvature vertices.
 11    r = np.array([[math.cos(.42), -math.sin(.42)], [math.sin(.42), math.cos(.42)]])
 12    h1 = np.diag([1.0, 8.0])
 13    h2 = r @ np.diag([1.5, 5.5]) @ r.T
 14    return [h1, h2]
 15
 16
 17def lmi_margin(Hs, K, alpha, sigma, P=None, lam=1.0):
 18    """Return largest eigenvalue of supplied S-procedure matrices.
 19    Negative means the fixed-P certificate is strict. P is normalized if absent.
 20    """
 21    n = Hs[0].shape[0]
 22    if P is None:
 23        # A valid Lyapunov candidate for the no-staleness closed loop.
 24        Abar = -K @ sum(Hs) / len(Hs)
 25        P = solve_discrete_lyapunov(np.eye(n) + .02*Abar, np.eye(n))
 26        P = P / np.trace(P) * n
 27    worst = -np.inf
 28    for H in Hs:
 29        A = -K @ H
 30        B = A.copy()
 31        tl = A.T@P + P@A + 2*alpha*P + lam*sigma*sigma*np.eye(n)
 32        M = np.block([[tl, P@B], [B.T@P, -lam*np.eye(n)]])
 33        worst = max(worst, eigvalsh((M+M.T)/2).max())
 34    return float(worst), P
 35
 36
 37def run(h, K, sigma, steps=2500, seed=0, max_interval=10):
 38    rng = np.random.default_rng(seed)
 39    Hs = make_hessians()
 40    x = np.array([2.0, -1.5])
 41    last = x.copy()
 42    total = 0.0
 43    events = 0
 44    max_norm = 0.0
 45    losses = []
 46    intervals = []
 47    since = 0
 48    # Alternating vertex is a deliberately adversarial but bounded curvature sequence.
 49    for t in range(steps):
 50        H = Hs[t % len(Hs)]
 51        g_stale = H @ last
 52        x = x - h * (K @ g_stale)
 53        since += 1
 54        # e proxy is x itself because x*=0, matching the proposed trigger proxy.
 55        d = last - x
 56        if (d @ d > sigma*sigma*(x @ x + 1e-12)) or since >= max_interval:
 57            last = x.copy()
 58            events += 1
 59            intervals.append(since)
 60            since = 0
 61        q = .5*x @ H @ x
 62        total += q
 63        losses.append(q)
 64        max_norm = max(max_norm, float(np.linalg.norm(x)))
 65        if not np.isfinite(q) or max_norm > 1e12:
 66            return dict(stable=False, final_loss=float('inf'), events=events,
 67                        rate=events/steps, mean_interval=np.mean(intervals) if intervals else steps,
 68                        max_norm=max_norm, losses=losses)
 69    return dict(stable=True, final_loss=float(losses[-1]), events=events,
 70                rate=events/steps, mean_interval=np.mean(intervals) if intervals else steps,
 71                max_norm=max_norm, losses=losses)
 72
 73
 74def main():
 75    Hs = make_hessians()
 76    # Dense K is aligned to average Hessian; diagonal baseline has same trace/average scale.
 77    Hbar = sum(Hs)/len(Hs)
 78    K_dense = np.linalg.inv(Hbar) * 0.72
 79    K_diag = np.diag(0.72/np.diag(Hbar))
 80    # Prediction 1: no-stale discrete quadratic stability boundary h*rho(KH)<2.
 81    print('HESSIANS', [h.tolist() for h in Hs])
 82    for name,K in [('diag',K_diag),('dense',K_dense)]:
 83        rho=max(abs(np.linalg.eigvals(K@H)) .max() for H in Hs)
 84        predicted=2/rho
 85        rows=[]
 86        for mult in [0.90, 1.02]:
 87            h=mult*predicted
 88            rr=run(h,K,0.0,steps=160,seed=1,max_interval=1)
 89            rows.append({'mult':mult,'h':h,'stable':rr['stable'],'max_norm':rr['max_norm']})
 90        print('STABILITY',name,json.dumps({'rho':float(rho),'predicted_h_boundary':predicted,'observed':rows}))
 91    # Prediction 2: trigger threshold zero gives stale-free updates; larger sigma permits longer intervals
 92    # and should increase error relative to a synchronized reference.
 93    h=.16
 94    print('TRIGGER_SWEEP')
 95    for name,K in [('diag',K_diag),('dense',K_dense)]:
 96        ref=run(h,K,0.0,steps=800,seed=3,max_interval=1)
 97        for sigma in [0.0,.01,.03,.08,.16,.30]:
 98            rr=run(h,K,sigma,steps=800,seed=3,max_interval=50)
 99            # average loss over final 20%, normalized to stale-free final loss floor
100            tail=float(np.mean(rr['losses'][-160:]))
101            ref_tail=float(np.mean(ref['losses'][-160:]))
102            print(json.dumps({'kind':'trigger','method':name,'sigma':sigma,
103                'events':rr['events'],'event_fraction':rr['rate'],'mean_interval':rr['mean_interval'],
104                'final_tail_loss':tail,'loss_ratio_vs_sync':tail/(ref_tail+1e-30),
105                'stable':rr['stable']}))
106    # Prediction 3: local first-order communication interval scales approximately sigma/(h ||KH||)
107    # for small sigma, so event fraction should scale approximately linearly with h/sigma.
108    print('SCALING_SWEEP')
109    K=K_dense
110    for h in [.04,.08,.16]:
111        for sigma in [.02,.04,.08]:
112            rr=run(h,K,sigma,steps=1000,seed=4,max_interval=1000)
113            print(json.dumps({'kind':'scaling','h':h,'sigma':sigma,'events':rr['events'],
114                'event_fraction':rr['rate'],'mean_interval':rr['mean_interval'],'stable':rr['stable']}))
115    # Fixed-K LMI sanity check: report monotone margin as sigma increases.
116    print('LMI_SWEEP')
117    for name,K in [('diag',K_diag),('dense',K_dense)]:
118        for sigma in [0.0,.02,.05,.1,.2]:
119            margin,_=lmi_margin(Hs,K,alpha=.02,sigma=sigma,lam=1.0)
120            print(json.dumps({'kind':'lmi','method':name,'sigma':sigma,'worst_eigenvalue':margin,'certified':margin<0}))
121
122if __name__ == '__main__':
123    main()