import json, math from pathlib import Path import numpy as np from sklearn.linear_model import LogisticRegression from sklearn.model_selection import train_test_split from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler from sklearn.metrics import accuracy_score SEED = 1461 rng = np.random.default_rng(SEED) def laplacian_path(n): A = np.zeros((n, n)) for i in range(n - 1): A[i, i + 1] = A[i + 1, i] = 1.0 return np.diag(A.sum(1)) - A def spectral_data(q): n = len(q) H = laplacian_path(n) + np.diag(q) lam, phi = np.linalg.eigh(H) S = phi ** 2 sums = np.array([lam[i] + lam[j] for i in range(n) for j in range(i, n)]) diffs = np.array([lam[i] - lam[j] for i in range(n) for j in range(i)]) gap_sum = np.min(np.abs(sums[:, None] - sums[None, :] + np.eye(len(sums)) * 1e9)) if len(sums) > 1 else np.inf gap_diff = np.min(np.abs(diffs[:, None] - diffs[None, :] + np.eye(len(diffs)) * 1e9)) if len(diffs) > 1 else np.inf overlap = min(np.max(np.abs(phi[:, i] * phi[:, j])) for i in range(n) for j in range(i)) return H, lam, phi, S, float(gap_sum), float(gap_diff), float(np.linalg.svd(S, compute_uv=False)[-1]), float(overlap) def observations(phi, lam, u, times): c = phi.T @ u return np.array([np.abs(phi @ (np.exp(-1j * t * lam) * c)) for t in times]).ravel() def phase_aligned_jacobian(phi, lam, u, times): # finite differences in 2n real coordinates; remove the global-phase tangent n = len(u) x = np.r_[u.real, u.imag] f0 = observations(phi, lam, u, times) J = np.empty((len(f0), 2*n)) eps = 2e-6 for k in range(2*n): xp = x.copy(); xp[k] += eps xm = x.copy(); xm[k] -= eps J[:, k] = (observations(phi, lam, xp[:n] + 1j*xp[n:], times) - observations(phi, lam, xm[:n] + 1j*xm[n:], times)) / (2*eps) # tangent to global phase is (-Im u, Re u), which is exactly a null direction tangent = np.r_[-u.imag, u.real] tangent /= np.linalg.norm(tangent) Q, _ = np.linalg.qr(np.column_stack([tangent, rng.normal(size=(2*n, 2*n-1))])) B = Q[:, 1:] # QR above may not preserve tangent as first column due nonorthogonal random columns; # use a robust null-space basis from SVD instead. _, _, vh = np.linalg.svd(tangent[None, :]) B = vh[1:].T return np.linalg.svd(J @ B, compute_uv=False) class SpectralPhaselessEncoder: """Fixed graph Schrodinger evolution exposed only through coordinate magnitudes.""" def __init__(self, graph_laplacian, potential, dt=1.0, steps=6): self.H = np.asarray(graph_laplacian, float) + np.diag(np.asarray(potential, float)) self.lam, self.phi = np.linalg.eigh(self.H) self.dt, self.steps = float(dt), int(steps) self.U = self.phi @ np.diag(np.exp(-1j * self.dt * self.lam)) @ self.phi.T def transform(self, h): x = np.asarray(h, complex).copy() out = [] for _ in range(self.steps): out.append(np.abs(x)) x = self.U @ x return np.concatenate(out) def diagnostics(self): n = len(self.lam) sums = np.array([self.lam[i]+self.lam[j] for i in range(n) for j in range(i,n)]) gap = np.min(np.abs(sums[:,None]-sums[None,:] + np.eye(len(sums))*1e9)) S = self.phi**2 overlap = min(np.max(np.abs(self.phi[:,i]*self.phi[:,j])) for i in range(n) for j in range(i)) return {"sigma_min_S": float(np.linalg.svd(S, compute_uv=False)[-1]), "pair_sum_gap": float(gap), "minimum_pair_overlap": float(overlap)} def math_checks(n=6): rows = [] for scale in [0.02, 0.1, 0.3, 1.0, 3.0]: vals = [] for _ in range(40): q = scale * rng.normal(size=n) *_, gap_sum, gap_diff, smin, overlap = spectral_data(q) vals.append((gap_sum, gap_diff, smin, overlap)) a = np.array(vals) rows.append({"potential_scale": scale, "median_pair_sum_gap": float(np.median(a[:,0])), "median_difference_gap": float(np.median(a[:,1])), "median_sigma_min_S": float(np.median(a[:,2])), "median_overlap": float(np.median(a[:,3]))}) q = rng.normal(size=n) _, lam, phi, S, gs, gd, sm, ov = spectral_data(q) u = rng.normal(size=n) + 1j*rng.normal(size=n) ranks = [] mins = [] for K in range(1, 9): times = np.linspace(0, 7.0, K) sv = phase_aligned_jacobian(phi, lam, u, times) ranks.append(int(np.sum(sv > 1e-5 * sv[0]))) mins.append(float(sv[-1])) # Noise amplification: perturb observations and solve a local least-squares correction. noise_rows = [] for K in [1, 2, 3, 4, 6, 8]: times = np.linspace(0, 7.0, K) sv = phase_aligned_jacobian(phi, lam, u, times) noise_rows.append({"K": K, "jacobian_rank": int(np.sum(sv > 1e-5*sv[0])), "smallest_nonphase_singular": float(sv[-1]), "predicted_noise_gain_1_over_sigma": float(1/sv[-1])}) return {"selected_operator": {"eigenvalues": lam.tolist(), "pair_sum_gap": gs, "difference_gap": gd, "sigma_min_S": sm, "overlap_min": ov}, "conditioning_sweep": rows, "time_sample_sweep": {"K": list(range(1,9)), "local_rank": ranks, "smallest_singular": mins}, "noise_sweep": noise_rows, "global_phase_max_abs_error": float(max(np.max(np.abs(observations(phi, lam, u, np.linspace(0,7,6)) - observations(phi, lam, u*np.exp(1j*1.234), np.linspace(0,7,6)))), 0.0))} def classification(n=6, N=2400): q = rng.normal(size=n) _, lam, phi, *_ = spectral_data(q) # Label is a phase-invariant spectral energy contrast; random global phases cannot affect it. U = rng.normal(size=(N,n)) + 1j*rng.normal(size=(N,n)) c = U @ phi label = ((np.abs(c[:,0])**2 + np.abs(c[:,1])**2) > (np.abs(c[:,2])**2 + np.abs(c[:,3])**2)).astype(int) phase = rng.uniform(0, 2*np.pi, N) U *= np.exp(1j*phase)[:,None] def enc(K): ts = np.linspace(0, 7.0, K) return np.concatenate([np.abs((phi @ (np.exp(-1j*t*lam)[:,None] * (U @ phi).T)).T) for t in ts], axis=1) Xidea = enc(6) Xbase = np.abs(U) tr, te = train_test_split(np.arange(N), test_size=.3, random_state=SEED, stratify=label) def fit(X): clf = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000, random_state=SEED)) clf.fit(X[tr], label[tr]); return accuracy_score(label[te], clf.predict(X[te])) # Add sensor noise at test time, matching the requested robustness check. noise = .05 * rng.normal(size=Xidea[te].shape) return {"clean_accuracy": {"phaseless_K6": float(fit(Xidea)), "single_time_coordinate_magnitude": float(fit(Xbase))}, "noisy_test_accuracy_phaseless_K6": float(accuracy_score(label[te], make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000, random_state=SEED)).fit(Xidea[tr], label[tr]).predict(Xidea[te] + noise))), "N": N, "train_test": [len(tr), len(te)]} if __name__ == "__main__": out = {"seed": SEED, "math": math_checks(), "classification": classification()} Path("results.json").write_text(json.dumps(out, indent=2)) print(json.dumps(out, indent=2))