1import json, math, time
 2from pathlib import Path
 3import numpy as np
 4
 5SEED = 1303
 6
 7def block_matrix(n_left, h, delta=1.0):
 8    """Uniformly sampled discretization of 1/(s+r), one dyadic r block."""
 9    s = np.arange(1, n_left + 1, dtype=float)[:, None] + delta - 1.0
10    r = h + np.arange(max(1, int(h)), dtype=float) + delta
11    return 1.0 / (s + r[None, :])
12
13
14def svd_observations():
15    rows = []
16    for h in [4, 8, 16, 32, 64]:
17        H = block_matrix(4 * h, h)
18        sv = np.linalg.svd(H / np.linalg.norm(H, 2), compute_uv=False)
19        # The paper predicts s_(N+1) <= 4^-N. This normalized check is
20        # conservative: normalization removes irrelevant scale factors.
21        bounds = [{"N": n, "observed": float(sv[n]), "predicted_bound": 4.0 ** (-n),
22                   "satisfies": bool(sv[n] <= 4.0 ** (-n) * (1 + 1e-10))}
23                  for n in range(min(5, len(sv)))]
24        rows.append({"h": h, "singular_values": [float(x) for x in sv[:6]],
25                     "geometric_bound_check": bounds,
26                     "ratios": [float(sv[i+1] / sv[i]) for i in range(min(4, len(sv)-1))]})
27    return rows
28
29
30def rank_scaling():
31    rows = []
32    for h in [8, 16, 32, 64]:
33        sv = np.linalg.svd(block_matrix(4*h, h), compute_uv=False)
34        sv /= sv[0]
35        for eps in [1e-2, 1e-3, 1e-4]:
36            observed = int(np.sum(sv > eps))
37            predicted = int(math.ceil(math.log(1 / eps, 4)))
38            rows.append({"h": h, "epsilon": eps, "observed_rank": observed,
39                         "predicted_rank_bound": predicted,
40                         "bound_holds": observed <= predicted})
41    return rows
42
43
44def application_scaling():
45    """Compare cached SVD application with dense multiplication over sizes."""
46    rng = np.random.default_rng(SEED)
47    rows = []
48    for h in [64, 128, 256, 512]:
49        nl, nr, d, rank = 4*h, h, 16, 4
50        H = block_matrix(nl, h)
51        U, sv, VT = np.linalg.svd(H, full_matrices=False)
52        A, B = U[:, :rank] * sv[:rank][None, :], VT[:rank, :].T
53        V = rng.standard_normal((nr, d))
54        dense = H @ V
55        approx = A @ (B.T @ V)
56        err = np.linalg.norm(dense - approx) / np.linalg.norm(dense)
57        # Small fixed repetitions; warmup avoids first-call allocation bias.
58        H @ V; A @ (B.T @ V)
59        reps = 8 if h <= 256 else 4
60        t0 = time.perf_counter()
61        for _ in range(reps): H @ V
62        td = (time.perf_counter() - t0) / reps
63        t0 = time.perf_counter()
64        for _ in range(reps): A @ (B.T @ V)
65        tf = (time.perf_counter() - t0) / reps
66        dense_work = nl * nr * d
67        fact_work = (nl + nr) * rank * d
68        rows.append({"h": h, "n_left": nl, "n_right": nr, "rank": rank,
69                     "relative_output_error": float(err), "dense_ms": td*1000,
70                     "factorized_ms": tf*1000, "measured_speedup": td/max(tf,1e-12),
71                     "ideal_work_ratio": dense_work / fact_work})
72    return rows
73
74
75def main():
76    out = {"seed": SEED, "singular_decay": svd_observations(),
77           "rank_scaling": rank_scaling(), "application_scaling": application_scaling()}
78    Path("results.json").write_text(json.dumps(out, indent=2))
79    print(json.dumps(out, indent=2))
80
81if __name__ == "__main__": main()