1import json
  2import numpy as np
  3
  4
  5def opnorm(a):
  6    return np.linalg.norm(a, 2)
  7
  8
  9def gram(stack, p=None):
 10    if p is None:
 11        p = stack[0].shape[1]
 12    return sum((w.T @ w for w in stack), np.zeros((p, p)))
 13
 14
 15def theorem_bound(true_stack, errors):
 16    return sum(2 * opnorm(w) * e + e * e for w, e in zip(true_stack, errors))
 17
 18
 19def computable_bound(estimated_stack, errors):
 20    # Since ||Omega|| <= ||Omega_hat|| + ||Delta||, this is a valid
 21    # estimated-state-only inflation of the theorem bound.
 22    return sum(2 * (opnorm(h) + e) * e + e * e
 23               for h, e in zip(estimated_stack, errors))
 24
 25
 26def greedy_stack(candidates, max_n=8):
 27    selected, selected_ids, current = [], [], -np.inf
 28    p = candidates[0].shape[1]
 29    for i, w in enumerate(candidates):
 30        trial = selected + [w]
 31        score = np.linalg.eigvalsh(gram(trial, p))[0]
 32        if len(selected) < p or score > current + 1e-10:
 33            selected, selected_ids, current = trial, selected_ids + [i], score
 34        if len(selected) >= max_n:
 35            break
 36    return selected, selected_ids
 37
 38
 39def mechanism_checks(rng):
 40    p, m = 2, 18
 41    true = [rng.normal(size=(3, p)) for _ in range(m)]
 42    true = [w / opnorm(w) for w in true]
 43
 44    rows = []
 45    violations = 0
 46    conservative_violations = 0
 47    for e in np.geomspace(1e-5, 0.8, 12):
 48        estimated, errors = [], []
 49        for w in true:
 50            d = rng.normal(size=w.shape); d *= e / opnorm(d)
 51            estimated.append(w + d); errors.append(e)
 52        diff = opnorm(gram(estimated) - gram(true))
 53        b = theorem_bound(true, errors)
 54        bc = computable_bound(estimated, errors)
 55        violations += int(diff > b * (1 + 1e-10))
 56        conservative_violations += int(diff > bc * (1 + 1e-10))
 57        rows.append((e, diff, b, bc, diff / b))
 58
 59    data = np.array(rows)
 60    # Prediction 1: exact perturbation is bounded for every mismatch size.
 61    # Prediction 2: theorem bound has O(e) scaling at small e.
 62    small_slope = np.median(data[:5, 1] / data[:5, 0])
 63    # Prediction 3: the quadratic term becomes non-negligible at large e;
 64    # report observed quadratic/theorem-bound contribution fraction.
 65    large_fraction = np.median((m * data[-4:, 0] ** 2) / data[-4:, 2])
 66
 67    # Certificate threshold with one fixed perturbation direction.  This
 68    # isolates the predicted monotone loss of the conservative certificate.
 69    gamma = 0.18
 70    directions = []
 71    for w in true:
 72        d = rng.normal(size=w.shape)
 73        directions.append(d / opnorm(d))
 74    threshold_rows = []
 75    for e in np.linspace(0, 0.55, 56):
 76        estimated = [w + e*d for w, d in zip(true, directions)]
 77        gh, gt = gram(estimated), gram(true)
 78        rho = theorem_bound(true, [e] * m)
 79        q = np.linalg.eigvalsh(gh)[0] - rho
 80        actual = np.linalg.eigvalsh(gt)[0]
 81        threshold_rows.append((e, q, actual))
 82    cert = np.array(threshold_rows)
 83    positive = cert[cert[:, 1] > gamma]
 84    observed_e = positive[-1, 0] if len(positive) else np.nan
 85    pred_e = np.nan
 86    for a, b in zip(cert[:-1], cert[1:]):
 87        if (a[1] - gamma) * (b[1] - gamma) <= 0:
 88            pred_e = a[0] + (gamma-a[1]) * (b[0]-a[0]) / (b[1]-a[1])
 89            break
 90    # Directly verify Weyl at every sweep point, not only q-positive points.
 91    cert_margin = float(np.min(cert[:, 2] - cert[:, 1]))
 92    return {
 93        'bound_violations': int(violations),
 94        'computable_bound_violations': int(conservative_violations),
 95        'bound_max_ratio': float(max(r[4] for r in rows)),
 96        'small_epsilon_slope_median': float(small_slope),
 97        'large_epsilon_quadratic_fraction': float(large_fraction),
 98        'gamma': gamma,
 99        'certificate_crossing_epsilon_observed': float(observed_e),
100        'certificate_crossing_epsilon_interpolated': float(pred_e),
101        'minimum_actual_minus_q': cert_margin,
102        'bound_table': [[float(x) for x in r] for r in rows],
103    }
104
105
106def identification(rng):
107    p, n_candidates = 2, 45
108    theta = np.array([1.25, -0.8])
109    candidates = []
110    for i in range(n_candidates):
111        angle = 0.12 * i
112        base = np.array([np.cos(angle), np.sin(angle)])
113        w = np.tile(base, (4, 1)) + 0.08 * rng.normal(size=(4, p))
114        candidates.append(w / opnorm(w))
115    noise = 0.035
116    y = [w @ theta + noise * rng.normal(size=4) for w in candidates]
117    greedy, idx = greedy_stack(candidates, max_n=8)
118    g = gram(greedy)
119
120    def estimate(ids):
121        H = gram([candidates[i] for i in ids], p)
122        b = sum((candidates[i].T @ y[i] for i in ids), np.zeros(p))
123        return np.linalg.solve(H + 1e-5*np.eye(p), b)
124
125    errors_gated, errors_all, activated = [], [], False
126    for t in range(1, n_candidates + 1):
127        errors_all.append(float(np.linalg.norm(estimate(range(t)) - theta)))
128        gs = [i for i in idx if i < t]
129        if len(gs) >= p and np.linalg.eigvalsh(gram([candidates[i] for i in gs], p))[0] > 0.05:
130            activated = True
131            est = estimate(gs)
132        else:
133            est = estimate(gs) if activated else np.zeros(p)
134        errors_gated.append(float(np.linalg.norm(est-theta)))
135    return {
136        'greedy_stack_size': len(greedy),
137        'greedy_lambda_min': float(np.linalg.eigvalsh(g)[0]),
138        'baseline_final_parameter_error': errors_all[-1],
139        'idea_final_parameter_error': errors_gated[-1],
140        'baseline_best_error': min(errors_all),
141        'idea_best_error': min(errors_gated),
142        'activation_step': next((t for t in range(1, n_candidates + 1) if len([i for i in idx if i < t]) >= p and np.linalg.eigvalsh(gram([candidates[i] for i in idx if i < t], p))[0] > 0.05), None),
143        'errors_baseline': errors_all,
144        'errors_idea': errors_gated,
145    }
146
147
148def main():
149    rng = np.random.default_rng(2914)
150    result = {'seed': 2914, 'mechanism': mechanism_checks(rng), 'identification': identification(rng)}
151    with open('results.json', 'w') as f: json.dump(result, f, indent=2)
152    print(json.dumps(result, indent=2))
153
154
155if __name__ == '__main__':
156    main()