1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2486
  6
  7def gale(X):
  8    """Return rows spanning the left nullspace of X (X is a x b)."""
  9    u, s, vh = np.linalg.svd(X.T, full_matrices=True)
 10    a, b = X.shape
 11    return u[:, a:].T, s
 12
 13def normalized_penalty(Y, X, D=None, eps=1e-12):
 14    if D is None:
 15        D = np.eye(X.shape[1])
 16    r = Y @ D @ X.T
 17    return np.linalg.norm(r, 'fro')**2 / (np.linalg.norm(Y, 'fro')**2 * np.linalg.norm(X, 'fro')**2 + eps)
 18
 19def prediction_sweeps(rng):
 20    # P1: exact SVD Gale residual is floating-point roundoff for every valid size.
 21    exact = []
 22    for a, b in [(2, 5), (3, 7), (5, 9), (8, 13)]:
 23        X = rng.normal(size=(a, b))
 24        Y, _ = gale(X)
 25        exact.append({'a': a, 'b': b, 'residual_norm': float(np.linalg.norm(Y @ X.T)),
 26                      'relative_penalty': float(normalized_penalty(Y, X)),
 27                      'rank': int(np.linalg.matrix_rank(X))})
 28
 29    # P2: perturbing Y by gamma Z predicts residual norm O(gamma), penalty O(gamma^2).
 30    a, b = 4, 9
 31    X = rng.normal(size=(a, b)); Y, _ = gale(X)
 32    Z = rng.normal(size=Y.shape); Z /= np.linalg.norm(Z)
 33    gammas = np.logspace(-6, -1, 8)
 34    rows = []
 35    for g in gammas:
 36        Yn = Y + g * Z
 37        rn = np.linalg.norm(Yn @ X.T)
 38        p = normalized_penalty(Yn, X)
 39        rows.append({'gamma': float(g), 'residual_norm': float(rn), 'penalty': float(p),
 40                     'norm_over_gamma': float(rn / g), 'penalty_over_gamma2': float(p / g**2)})
 41    slope_norm = float(np.polyfit(np.log(gammas), np.log([r['residual_norm'] for r in rows]), 1)[0])
 42    slope_pen = float(np.polyfit(np.log(gammas), np.log([r['penalty'] for r in rows]), 1)[0])
 43
 44    # P3: D=cI preserves the null relation; a nonconstant diagonal gauge generally breaks it.
 45    Drows = []
 46    for spread in [0., .01, .03, .1, .3, 1.0]:
 47        logd = rng.normal(0, spread, size=b)
 48        D = np.diag(np.exp(logd))
 49        Drows.append({'spread': spread, 'residual_norm': float(np.linalg.norm(Y @ D @ X.T)),
 50                      'penalty': float(normalized_penalty(Y, X, D))})
 51    positive = [r for r in Drows if r['spread'] > 0]
 52    gauge_slope = float(np.polyfit(np.log([r['spread'] for r in positive]), np.log([r['residual_norm'] for r in positive]), 1)[0])
 53    return exact, rows, slope_norm, slope_pen, Drows, gauge_slope
 54
 55def fit_ridge(A, y, lam=1e-6):
 56    # Deterministic least-squares ridge, with intercept.
 57    A1 = np.column_stack([A, np.ones(len(A))])
 58    reg = lam * np.eye(A1.shape[1]); reg[-1, -1] = 0
 59    return np.linalg.solve(A1.T @ A1 + reg, A1.T @ y)
 60
 61def toy_comparison(rng):
 62    # A tiny representation probe: primary summary has a dimensions; adding
 63    # YX^T supplies an orthogonal (b-a)-dimensional feature summary.
 64    a, b, d, n = 3, 8, 12, 400
 65    W = rng.normal(size=(a, d)); W /= np.linalg.norm(W, axis=1, keepdims=True)
 66    primary, combined, targets = [], [], []
 67    for _ in range(n):
 68        H = rng.normal(size=(d, b))
 69        P = W @ H
 70        Y, _ = gale(P)
 71        dual = Y @ H.T                         # (b-a) x d
 72        p = P.mean(axis=1)
 73        q = dual.mean(axis=1)
 74        primary.append(p)
 75        combined.append(np.concatenate([p, q]))
 76        # Target depends on both channels; this tests whether the complement
 77        # carries useful information, not whether it changes the null identity.
 78        targets.append(1.2 * p[0] - .7 * p[1] + .5 * q[0] + .8 * q[1] + .05 * rng.normal())
 79    primary, combined, targets = map(np.asarray, (primary, combined, targets))
 80    cut = 280
 81    wb = fit_ridge(primary[:cut], targets[:cut]); wg = fit_ridge(combined[:cut], targets[:cut])
 82    pb = np.column_stack([primary[cut:], np.ones(n-cut)]) @ wb
 83    pg = np.column_stack([combined[cut:], np.ones(n-cut)]) @ wg
 84    return {'primary_test_mse': float(np.mean((pb-targets[cut:])**2)),
 85            'gale_combined_test_mse': float(np.mean((pg-targets[cut:])**2)),
 86            'feature_width_primary': a, 'feature_width_combined': a + (b-a)}
 87
 88def main():
 89    rng = np.random.default_rng(SEED)
 90    exact, rows, sn, sp, drows, gs = prediction_sweeps(rng)
 91    result = {'seed': SEED, 'predictions': {
 92        'P1_exact_nullspace': '||YX^T|| remains at floating-point roundoff',
 93        'P2_noise_scaling': '||YX^T|| ~ gamma^1 and normalized penalty ~ gamma^2',
 94        'P3_diagonal_gauge': 'D=cI preserves zero; nonconstant D produces residual'},
 95        'exact_sweep': exact, 'noise_sweep': rows,
 96        'noise_loglog_slopes': {'residual_norm': sn, 'penalty': sp},
 97        'gauge_loglog_slope': gs,
 98        'diagonal_sweep': drows, 'toy_comparison': toy_comparison(rng)}
 99    Path('results.json').write_text(json.dumps(result, indent=2))
100    print(json.dumps(result, indent=2))
101
102if __name__ == '__main__':
103    main()