1#!/usr/bin/env python3
  2"""Schur-Riesz Greedy Adapter Expansion numerical MVP."""
  3import json
  4import numpy as np
  5
  6SEED = 2145
  7
  8
  9def weighted_projector(B, y, rtol=1e-11):
 10    y = np.asarray(y, float)
 11    n = len(y)
 12    if B.size == 0:
 13        return np.zeros((n, n)), np.diag(np.sqrt(y))
 14    Y = np.diag(y)
 15    gram = B.T @ Y @ B
 16    P = B @ np.linalg.pinv(gram, rcond=rtol) @ B.T @ Y
 17    return P, np.diag(np.sqrt(y))
 18
 19
 20def residual_response(B, C, y):
 21    P, Ys = weighted_projector(B, y)
 22    Q = C - P @ C
 23    WQ = Ys @ Q
 24    sv = np.linalg.svd(WQ, compute_uv=False)
 25    scale = max(1.0, float(sv[0]) if len(sv) else 0.0)
 26    positive = sv[sv > 1e-9 * scale]
 27    return Q, float(np.sum(WQ * WQ)), positive
 28
 29
 30def fit_loss(X, target, y):
 31    sw = np.sqrt(y)
 32    Xw, tw = sw[:, None] * X, sw * target
 33    residual = tw - Xw @ np.linalg.pinv(Xw) @ tw
 34    return 0.5 * float(residual @ residual)
 35
 36
 37def exact_increment(X, C, target, y):
 38    before = fit_loss(X, target, y)
 39    after = fit_loss(np.column_stack([X, C]), target, y)
 40    Q, _, _ = residual_response(X, C, y)
 41    sw = np.sqrt(y)
 42    Qw, rw = sw[:, None] * Q, sw * target
 43    predicted = 0.5 * float(rw @ Qw @ np.linalg.pinv(Qw.T @ Qw) @ Qw.T @ rw)
 44    return before - after, predicted
 45
 46
 47def make_problem(seed=SEED, n=180, d=3, candidates=18):
 48    r = np.random.default_rng(seed)
 49    y = np.exp(r.normal(0, .35, n))
 50    base = r.normal(size=(n, d))
 51    # Several candidates are redundant, several are useful, and several are noise.
 52    blocks = []
 53    labels = []
 54    for j in range(candidates):
 55        if j < 5:
 56            C = base[:, [j % d]] + .015 * r.normal(size=(n, 1))
 57            labels.append("redundant")
 58        elif j in (5, 7, 11):
 59            C = r.normal(size=(n, 1))
 60            labels.append("useful")
 61        else:
 62            C = .07 * r.normal(size=(n, 1))
 63            labels.append("weak")
 64        blocks.append(C)
 65    target = 1.8 * blocks[5][:, 0] - 1.2 * blocks[7][:, 0] + .25 * base[:, 0] + .35 * r.normal(size=n)
 66    return y, base, blocks, labels, target
 67
 68
 69def math_checks():
 70    r = np.random.default_rng(SEED)
 71    n = 80
 72    y = np.exp(r.normal(size=n))
 73    B = r.normal(size=(n, 3))
 74    C = r.normal(size=(n, 2))
 75    P, Ys = weighted_projector(B, y)
 76    Q, gain, sv = residual_response(B, C, y)
 77    orth = np.linalg.norm(B.T @ (y[:, None] * Q))
 78    idem = np.linalg.norm(P @ P - P)
 79
 80    # Prediction 1: residual gain vanishes as candidate approaches incumbent.
 81    scales = [0.0, .25, .5, 1.0, 2.0]
 82    # candidate = incumbent direction + scale * novel direction
 83    novel = r.normal(size=(n, 1))
 84    gain_scale = []
 85    for a in scales:
 86        _, g, _ = residual_response(B, B[:, [0]] + a * novel, y)
 87        gain_scale.append(g)
 88    # Prediction 2: gain scales quadratically with candidate amplitude.
 89    amps = [.25, .5, 1., 2.]
 90    amp_gain = [residual_response(B, a * novel, y)[1] for a in amps]
 91    ratios = [amp_gain[i] / amp_gain[2] for i in range(4)]
 92    # Prediction 3: exact least-squares improvement equals projected residual energy.
 93    target = r.normal(size=n)
 94    observed, predicted = exact_increment(B, C, target, y)
 95    return {
 96        "projection_orthogonality_abs": orth,
 97        "projector_idempotence_abs": idem,
 98        "scale_sweep": {str(a): g for a, g in zip(scales, gain_scale)},
 99        "amplitude_sweep_gain": {str(a): g for a, g in zip(amps, amp_gain)},
100        "amplitude_ratios_vs_a1": {str(a): ratios[i] for i, a in enumerate(amps)},
101        "expected_amplitude_ratios": {str(a): a*a for a in amps},
102        "exact_loss_drop": observed,
103        "predicted_loss_drop": predicted,
104        "relative_prediction_error": abs(observed-predicted) / max(1e-12, abs(observed)),
105        "positive_singular_values": sv.tolist(),
106    }
107
108
109def greedy_experiment():
110    y, base, blocks, labels, target = make_problem()
111    # Baseline: fixed order; idea: largest projected energy, with a sensible
112    # lower bound that rejects nearly-null blocks. Both use one scalar block.
113    def run(order):
114        X = base.copy()
115        losses = [fit_loss(X, target, y)]
116        selected = []
117        for j in order:
118            Q, g, sv = residual_response(X, blocks[j], y)
119            if len(sv) and sv[0] >= .15 and sv[0] <= 30:
120                selected.append(j)
121                X = np.column_stack([X, blocks[j]])
122                losses.append(fit_loss(X, target, y))
123            if len(selected) == 3: break
124        return selected, losses
125    fixed, fixed_losses = run(list(range(len(blocks))))
126    # Greedy recomputes residual gains each round.
127    remaining = set(range(len(blocks)))
128    X = base.copy(); greedy = [fit_loss(X, target, y)]; accepted = []; gains = []
129    for _ in range(3):
130        scored = []
131        for j in remaining:
132            _, g, sv = residual_response(X, blocks[j], y)
133            if len(sv) and .15 <= sv[0] <= 30:
134                scored.append((g, j, float(sv[0])))
135        if not scored: break
136        g, j, s = max(scored)
137        accepted.append(j); gains.append(g)
138        X = np.column_stack([X, blocks[j]]); remaining.remove(j)
139        greedy.append(fit_loss(X, target, y))
140    return {
141        "fixed_order_selected": fixed,
142        "fixed_order_labels": [labels[j] for j in fixed],
143        "fixed_losses": fixed_losses,
144        "greedy_selected": accepted,
145        "greedy_labels": [labels[j] for j in accepted],
146        "greedy_projected_gains": gains,
147        "greedy_losses": greedy,
148        "final_loss_ratio_greedy_over_fixed": greedy[-1] / fixed_losses[-1],
149    }
150
151
152def main():
153    out = {"seed": SEED, "math_checks": math_checks(), "toy_experiment": greedy_experiment()}
154    print(json.dumps(out, indent=2))
155
156
157if __name__ == "__main__":
158    main()