import json, math, time from pathlib import Path import numpy as np SEED = 2784 rng = np.random.default_rng(SEED) def envelope(A_factors): A = np.eye(A_factors[0].shape[0], dtype=np.complex128) P = np.eye(A.shape[0], dtype=np.float64) for F in A_factors: A = F @ A P = np.abs(F) @ P return A, P def power_radius(P, iters=100): v = np.ones(P.shape[0]) / math.sqrt(P.shape[0]) for _ in range(iters): z = P @ v n = np.linalg.norm(z) if n == 0: return 0.0 v = z / n return float(np.linalg.norm(P @ v) / (np.linalg.norm(v) + 1e-15)) def verify(): # Prediction 1: |F_L...F_1| <= product |F_k| entrywise for every draw. violations, ratios = [], [] for _ in range(5000): L = int(rng.integers(1, 5)) fs = [(rng.normal(size=(4,4)) + 1j*rng.normal(size=(4,4))) / math.sqrt(8) for _ in range(L)] A, P = envelope(fs) violations.append(float(np.max(np.abs(A)-P))) ratios.append(float(np.max(np.abs(A)/(P+1e-15)))) max_violation = max(violations) # Prediction 2: two equal paths have cancellation ratio |cos(theta/2)|. thetas = np.linspace(0, math.pi, 9) cancel_obs = [] cancel_pred = [] for th in thetas: z = 1.0 + np.exp(1j*th) cancel_obs.append(abs(z)/2.0) cancel_pred.append(abs(math.cos(th/2))) # Prediction 3: scalar repeated dynamics. Envelope threshold is gamma=1/2, # while actual threshold is gamma=1/(2|cos(theta/2)|). gamma_grid = np.linspace(.1, 4.0, 391) threshold_rows = [] for th in [0.0, math.pi/2, math.pi*0.9, math.pi]: c = abs(math.cos(th/2)) env_factor = 2.0 * gamma_grid actual_factor = 2.0 * c * gamma_grid env_cross = gamma_grid[np.argmax(env_factor >= 1.0)] act_cross = gamma_grid[np.argmax(actual_factor >= 1.0)] if c > 1e-12 and actual_factor.max() >= 1.0 else None threshold_rows.append({ 'theta': float(th), 'predicted_envelope_gamma': .5, 'observed_envelope_gamma': float(env_cross), 'predicted_actual_gamma': None if c == 0 else float(1/(2*c)), 'observed_actual_gamma': None if act_cross is None else float(act_cross), 'cancellation_factor': float(c), }) return { 'entrywise_bound': {'trials': 5000, 'max_violation': max_violation, 'max_absA_over_P': max(ratios)}, 'phase_cancellation': [{'theta': float(t), 'predicted_ratio': float(p), 'observed_ratio': float(o)} for t,p,o in zip(thetas,cancel_pred,cancel_obs)], 'growth_thresholds': threshold_rows, } def torch_benchmark(): # Optional torch benchmark: learn a stable scalar-output complex recurrence. # Both models have identical factorized transition parameterization and steps. try: import torch torch.manual_seed(SEED) device = 'cuda' if torch.cuda.is_available() else 'cpu' try: torch.zeros(1, device=device).sum().item() except Exception: device = 'cpu' dtype = torch.complex64 n, T, B, steps = 4, 35, 64, 350 g = torch.Generator(device=device).manual_seed(SEED) # Stable teacher with long memory; inputs are real but state is complex. F1t = (0.34*torch.eye(n, device=device, dtype=dtype) + 0.025*(torch.randn(n,n, generator=g, device=device)+1j*torch.randn(n,n,generator=g,device=device))) F2t = (0.90*torch.eye(n, device=device, dtype=dtype) + 0.025*(torch.randn(n,n, generator=g, device=device)+1j*torch.randn(n,n,generator=g,device=device))) At = F2t @ F1t x = torch.randn(B,T,1, generator=g, device=device, dtype=torch.float32).to(dtype) with torch.no_grad(): h = torch.zeros(B,n,device=device,dtype=dtype); ys=[] for t in range(T): h = h @ At.T + x[:,t] @ torch.ones(1,n,device=device,dtype=dtype)*.12 ys.append(h[:,0].real) y = torch.stack(ys,1) def run(kind): torch.manual_seed(SEED+ (0 if kind=='baseline' else 1)) F1 = torch.nn.Parameter((.72*torch.eye(n,dtype=dtype,device=device) + .11*(torch.randn(n,n,device=device)+1j*torch.randn(n,n,device=device))) ) F2 = torch.nn.Parameter((.72*torch.eye(n,dtype=dtype,device=device) + .11*(torch.randn(n,n,device=device)+1j*torch.randn(n,n,device=device))) ) read = torch.nn.Parameter(torch.randn(n,device=device,dtype=dtype)*.1) opt = torch.optim.Adam([F1,F2,read], lr=.012) t0=time.perf_counter(); unstable=0 for _ in range(steps): opt.zero_grad(); A=F2@F1; P=torch.abs(F2)@torch.abs(F1) h=torch.zeros(B,n,device=device,dtype=dtype); pred=[] for t in range(T): h=h@A.T + x[:,t]@torch.ones(1,n,device=device,dtype=dtype)*.12 pred.append((h*read).sum(1).real) loss=torch.stack(pred,1).sub(y).pow(2).mean() if kind=='envelope': # Smooth differentiable conservative penalty: row-sum envelope above target. env=P.sum(1).amax() loss=loss + .08*torch.relu(env-.96).pow(2) + .01*(P.sum()-torch.abs(A).sum())/(P.sum()+1e-6) loss.backward() gn=torch.nn.utils.clip_grad_norm_([F1,F2,read], 10.0) if not torch.isfinite(loss) or not torch.isfinite(gn): unstable += 1 opt.step() with torch.no_grad(): A=F2@F1; P=torch.abs(F2)@torch.abs(F1) # long rollout against teacher from fresh sequence xx=torch.randn(B,120,1,generator=g,device=device,dtype=torch.float32).to(dtype) hp=torch.zeros(B,n,device=device,dtype=dtype); ht=torch.zeros_like(hp) err=[] for t in range(120): hp=hp@A.T+xx[:,t]@torch.ones(1,n,device=device,dtype=dtype)*.12 ht=ht@At.T+xx[:,t]@torch.ones(1,n,device=device,dtype=dtype)*.12 err.append(((hp*read).sum(1).real-(ht[:,0].real))**2) longerr=torch.stack(err,1).mean().sqrt().item() rhoP=power_radius(P.detach().cpu().numpy()) rhoA=max(abs(np.linalg.eigvals(A.detach().cpu().numpy()))) return {'final_train_loss':float(loss.item()),'long_rollout_rmse':longerr, 'unstable_steps':unstable,'rho_A':float(rhoA),'rho_P':float(rhoP), 'seconds':time.perf_counter()-t0,'device':device} return {'baseline':run('baseline'),'envelope':run('envelope')} except Exception as e: return {'error':repr(e), 'note':'torch benchmark failed; algebraic verification remains valid'} if __name__ == '__main__': out={'seed':SEED,'verification':verify(),'benchmark':torch_benchmark()} Path('results.json').write_text(json.dumps(out, indent=2)) print(json.dumps(out, indent=2))