import json, math, time, random from pathlib import Path import numpy as np import torch SEED = 1620 np.random.seed(SEED); random.seed(SEED); torch.manual_seed(SEED) if torch.cuda.is_available(): try: torch.cuda.manual_seed_all(SEED) except Exception: pass try: device = torch.device("cuda" if torch.cuda.is_available() else "cpu") # A tiny allocation makes CUDA failures visible before the experiment. if device.type == "cuda": torch.zeros(1, device=device) except Exception: device = torch.device("cpu") def rand_u2(rng): z = rng.normal(size=(2,2)) + 1j*rng.normal(size=(2,2)) q, r = np.linalg.qr(z) q = q @ np.diag(np.conj(np.diag(r))/np.abs(np.diag(r))) return q.astype(np.complex128) def apply_butterfly(x, gates, scale=1.0): """Apply log2(N) stages of local 2x2 gates to vectors or a matrix.""" n = x.shape[-1] y = np.array(x, dtype=np.complex128, copy=True) levels = int(math.log2(n)) for s in range(levels): half = 2**s step = 2*half out = y.copy() for base in range(0, n, step): for j in range(half): inds = [base+j, base+j+half] out[..., inds] = np.einsum('ab,...b->...a', scale*gates[s][base//step*half+j], y[..., inds]) y = out return y def full_transform(n, rng, scale=1.0): gates = [[rand_u2(rng) for _ in range(n//2)] for _ in range(int(math.log2(n)))] # apply_butterfly expects one gate per pair at every stage cols = np.eye(n, dtype=np.complex128) return apply_butterfly(cols, gates, scale=scale), gates def unitary_checks(): rng = np.random.default_rng(SEED) rows=[] for n in [4, 8, 16, 32, 64]: T, gates = full_transform(n, rng) gram_err = float(np.linalg.norm(T.conj().T @ T - np.eye(n), ord='fro')) x = rng.normal(size=n)+1j*rng.normal(size=n) ratio = float(np.linalg.norm(T@x)/np.linalg.norm(x)) rows.append({'N':n, 'levels':int(math.log2(n)), 'unitarity_fro_error':gram_err, 'norm_ratio':ratio, 'predicted_norm_ratio':1.0}) # Controlled violation: multiplying each local gate by (1+eps) gives exactly # (1+eps)^levels for every vector, because each path crosses one gate per level. drift=[] n=32; x=rng.normal(size=n)+1j*rng.normal(size=n) # Sweep both perturbation strength and repeated depth. The exact prediction is # (1+eps)^(log2(N)*repetitions), since each path crosses one gate per level. for eps in [0.0, 1e-3, 1e-2, 5e-2]: rng_eps=np.random.default_rng(SEED+1) Tp, _ = full_transform(n, rng_eps, scale=1+eps) for depth in [1,4,8]: observed=float(np.linalg.norm(np.linalg.matrix_power(Tp, depth)@x)/np.linalg.norm(x)) predicted=float((1+eps)**(int(math.log2(n))*depth)) drift.append({'epsilon':eps, 'repetitions':depth, 'observed':observed, 'predicted':predicted, 'relative_error':abs(observed-predicted)/(predicted if predicted else 1)}) # Since T is unitary, its conjugate transpose is its exact inverse. rng_inv=np.random.default_rng(SEED+2); T,_=full_transform(32,rng_inv) inverse_err=float(np.linalg.norm(T.conj().T@T@np.ones(32)-np.ones(32))) return rows, drift, inverse_err class StructuredMixer(torch.nn.Module): # Four real coordinates parameterize every U(2), using a standard SU(2) chart # plus a global phase. Gates are shared over samples but not over butterfly pairs. def __init__(self, n, init=None): super().__init__(); self.n=n; self.levels=int(math.log2(n)) self.p=torch.nn.Parameter(torch.zeros(self.levels, n//2, 4)) if init is not None: with torch.no_grad(): self.p.copy_(init) def gates(self): a,b,c,d=self.p.unbind(-1) ca=torch.cos(c); sa=torch.sin(c) phase=torch.exp(1j*a) row0=torch.stack([torch.exp(1j*b)*ca, torch.exp(1j*d)*sa], -1) row1=torch.stack([-torch.exp(-1j*d)*sa, torch.exp(-1j*b)*ca], -1) return phase[...,None,None]*torch.stack([row0,row1],-2) def forward(self, x): y=x gs=self.gates() for s in range(self.levels): half=2**s; step=2*half; out=y.clone() for base in range(0,self.n,step): for j in range(half): ids=[base+j,base+j+half] out[:,ids]=torch.einsum('ab,nb->na',gs[s,base//step*half+j],y[:,ids]) y=out return y def fit_toy(): torch.manual_seed(SEED); n=16; samples=256 # Generate a target exactly in the proposed topology, making identifiability fair. true=StructuredMixer(n).to(device) with torch.no_grad(): true.p.normal_(0,0.35) x=(torch.randn(samples,n)+1j*torch.randn(samples,n)).to(device) with torch.no_grad(): y=true(x).detach() structured=StructuredMixer(n).to(device) dense=torch.nn.Parameter((0.05*(torch.randn(n,n)+1j*torch.randn(n,n))).to(device)) os=torch.optim.Adam(structured.parameters(),lr=0.06) od=torch.optim.Adam([dense],lr=0.06) losses_s=[]; losses_d=[]; t0=time.time() for step in range(301): os.zero_grad(); pred=structured(x); ls=(pred-y).abs().pow(2).mean(); ls.backward(); os.step() od.zero_grad(); pred2=x@dense.T; ld=(pred2-y).abs().pow(2).mean(); ld.backward(); od.step() if step in [0,50,100,200,300]: losses_s.append(float(ls.detach().cpu())); losses_d.append(float(ld.detach().cpu())) elapsed=time.time()-t0 return {'N':n,'samples':samples,'steps':300,'structured_params':int(sum(p.numel() for p in structured.parameters())), 'dense_complex_params':int(dense.numel()),'structured_final_mse':losses_s[-1], 'dense_final_mse':losses_d[-1],'structured_curve':losses_s,'dense_curve':losses_d, 'seconds':elapsed,'device':str(device)} def main(): checks, drift, inverse_err=unitary_checks() fit=fit_toy() result={'seed':SEED,'predictions':[ 'For every N=2^k, ||T x||/||x|| = 1 exactly up to floating point.', 'T^*T-I is at machine precision and does not grow materially with N.', 'If every local gate is scaled by 1+eps, repeated-pass norm drift is (1+eps)^(k*r), k=log2(N).' ],'unitarity_sweep':checks,'controlled_drift':drift,'inverse_reconstruction_error':inverse_err,'toy_fit':fit} Path('results.json').write_text(json.dumps(result,indent=2)) print(json.dumps(result,indent=2)) if __name__=='__main__': main()