#!/usr/bin/env python3 """ PROVENANCE: SURROGATE Computation 129 -- a spectral gap reaches the relaxation exponent, NOT the stationary-mean coefficient: the A6 bias residual is immune to mixing rates ============================================================================ Backs the open_research.md ยง1 A6 bias-half reduction. The (R) residual binds at the LEADING COEFFICIENT eta_{K,R} < d_0^2/4 of the discrepancy | A_n(u) - (d_0^2/4)||grad u||^2 | <= eta_{K,R} ||u||^2_{H^1}, NOT at an exponent. The natural-looking route flagged in Remark 2 -- a covariance-decay / spectral-gap mixing estimate -- controls the STATISTICAL FLUCTUATION of A_n(u) about its mean (the O(|D_n|^{-1/2}) term, already settled by McDiarmid, Computation 125), but is STRUCTURALLY BLIND to the stationary-mean VALUE, i.e. the systematic equidistribution defect that sets the coefficient. A spectral gap is a relaxation RATE; the stationary distribution is the top (kernel) eigenvector, decoupled from the gap. So no mixing exponent, however sharp, reaches eta_{K,R}: the residual lives in ker L, not its complement. The split mirrors the paper's own systematic (O(d_0^4), deterministic) vs statistical (O(d_0^2 |D_n|^{-1/2}), mixing-controlled) decomposition. This script demonstrates the decoupling two ways. CHECKS ------ (1) GAP BLIND TO THE MEAN. A one-parameter family of reversible 2-state Markov chains with FIXED spectral gap g but stationary mean swept across (0,1): identical gap, different mean. The gap does not pin the stationary mean. (2) FLUCTUATION vs SYSTEMATIC BIAS in A_n(u). Points drawn with a controlled equidistribution defect eps (density 1 + eps*cos 2pi x) plus sampling noise: the std of A_n(u) decays as |D_n|^{-1/2} (the mixing/gap-controlled fluctuation), while the systematic mean-deviation stays ~ eps*G INDEPENDENT of |D_n|. Sharper mixing (larger N) shrinks the fluctuation but never the systematic coefficient. (3) Consequence: eta_{K,R} = systematic (deterministic equidistribution) part + statistical (mixing) part; only the latter is gap-accessible, so the binding coefficient half is foreclosed to spectral-gap / mixing methods. A "pass": (1) gap == g to machine precision across a >=0.8-wide mean sweep; (2) mean-deviation ~ eps*G flat in N while std halves per 4x N (slope ~ -1/2 in log-log); both confirm the coefficient is decoupled from any relaxation rate. """ from __future__ import annotations import numpy as np def two_state(g, p): """Reversible 2-state chain, stationary (1-p, p), built to have gap exactly g. P = [[1-gp, gp],[g(1-p), 1-g(1-p)]] has eigenvalues 1 and 1-g, stationary (1-p,p).""" a, b = g * p, g * (1.0 - p) P = np.array([[1 - a, a], [b, 1 - b]]) lam = np.sort(np.linalg.eigvals(P).real)[::-1] gap = 1.0 - lam[1] mean = float(np.array([1 - p, p]) @ np.array([0.0, 1.0])) # state values {0,1} return gap, mean, P def sample_defect(N, eps, rng): """N points on [0,1] from density 1 + eps*cos(2 pi x) by rejection (envelope 1+eps).""" out = np.empty(0) while out.size < N: x = rng.uniform(0.0, 1.0, size=2 * N) keep = rng.uniform(0.0, 1.0, size=2 * N) < (1 + eps * np.cos(2 * np.pi * x)) / (1 + eps) out = np.concatenate([out, x[keep]]) return out[:N] def main(): print("=" * 100) print(" Computation 129 -- spectral gap reaches the exponent, not the stationary-mean coefficient (A6)") print("=" * 100) print() # ---- (1) gap blind to the stationary mean ------------------------------- print(" (1) FIXED spectral gap g = 0.5, sweep the stationary mean p over (0,1):") g = 0.5 ps = (0.1, 0.3, 0.5, 0.7, 0.9) gaps, means = [], [] for p in ps: gp, m, _ = two_state(g, p) gaps.append(gp); means.append(m) print(f" p = {p:.1f}: spectral gap = {gp:.6f} stationary mean = {m:.4f}") p1 = (max(abs(x - g) for x in gaps) < 1e-9) and (max(means) - min(means) >= 0.8) print(f" -> identical gap {g} across a {max(means)-min(means):.1f}-wide mean range: {p1}") print(f" a relaxation rate does NOT determine the stationary mean (the coefficient).") print() # ---- (2) fluctuation (gap-controlled) vs systematic bias (coefficient) --- eps = 0.3 G = 0.25 # G = \int_0^1 cos^2(pi x) cos(2 pi x) dx = 1/4 target = 0.5 # \int_0^1 cos^2(pi x) dx bias_pred = eps * G print(f" (2) A_n(u) with u' s.t. |grad u|^2 = cos^2(pi x); density defect eps = {eps}:") print(f" predicted systematic bias eps*G = {bias_pred:.4f} (n-INDEPENDENT); fluctuation ~ N^(-1/2)") print(f" {'N':>8} {'mean(A_n)-target':>18} {'std(A_n)':>12}") rng = np.random.default_rng(20260628) gfun = lambda x: np.cos(np.pi * x) ** 2 Ns = (250, 1000, 4000, 16000) biases, stds = [], [] for N in Ns: vals = np.array([gfun(sample_defect(N, eps, rng)).mean() for _ in range(300)]) biases.append(vals.mean() - target); stds.append(vals.std()) print(f" {N:>8} {vals.mean()-target:>+18.4f} {vals.std():>12.5f}") # systematic bias flat in N; std halves per 4x N bias_flat = max(abs(b - bias_pred) for b in biases) < 0.01 ratios = [stds[i] / stds[i + 1] for i in range(len(stds) - 1)] std_decay = abs(np.mean(ratios) - 2.0) < 0.4 # ~2x per 4x N <=> slope -1/2 print(f" systematic bias flat at ~{bias_pred:.3f} (max dev {max(abs(b-bias_pred) for b in biases):.4f}): {bias_flat}") print(f" std ratio per 4x N = {[round(r,2) for r in ratios]} (~2.0 => slope -1/2): {std_decay}") print() # ---- assessment --------------------------------------------------------- print("=" * 100) print(" ASSESSMENT") print("=" * 100) print(f" (1) gap blind to the stationary mean (same gap, 0.8-wide mean sweep) : {p1}") print(f" (2) systematic coefficient flat in N while the fluctuation decays : {bias_flat and std_decay}") ok = p1 and bias_flat and std_decay print() msg = ("CONFIRMS -- a spectral gap reaches only the fluctuation exponent; the systematic-mean " "coefficient eta_{K,R} is decoupled, so the A6 bias half is foreclosed to mixing/spectral-gap " "methods") if ok else "MISMATCH -- revise" print(" RESULT: " + msg) if __name__ == "__main__": main()