1import json
 2import time
 3import numpy as np
 4from coreset import residual_coreset, fit_min_norm
 5
 6rng = np.random.default_rng(2878)
 7
 8def make_data(n, d, r, m, noise=0.25):
 9    q, _ = np.linalg.qr(rng.normal(size=(d, r)))
10    x = rng.normal(size=(n, r)) @ q.T
11    theta = rng.normal(size=(m, d))
12    y = x @ theta.T + noise * rng.normal(size=(n, m))
13    return x, y
14
15def one(n, d, r, m):
16    x, y = make_data(n, d, r, m)
17    Wfull = fit_min_norm(x, y)
18    t0 = time.perf_counter()
19    idx, ww, _, rr, it = residual_coreset(x, y)
20    build_sec = time.perf_counter() - t0
21    wc = np.zeros(n); wc[idx] = ww
22    residual = y - x @ Wfull.T
23    full_normal = residual.T @ x
24    core_normal = residual.T @ (wc[:, None] * x)
25    inv_err = np.linalg.norm(core_normal-full_normal) / max(np.linalg.norm(y.T@x), 1e-12)
26    t0 = time.perf_counter(); Wc = fit_min_norm(x[idx], y[idx], ww); core_fit_sec = time.perf_counter()-t0
27    sol_err = np.linalg.norm(Wc-Wfull) / max(np.linalg.norm(Wfull), 1e-12)
28    full_loss = np.sum((x @ Wfull.T-y)**2)
29    core_at_full = np.sum((x @ Wc.T-y)**2)
30    ridx = rng.choice(n, size=len(idx), replace=False)
31    rw = np.full(len(ridx), n/len(ridx))
32    Wr = fit_min_norm(x[ridx], y[ridx], rw)
33    random_err = np.linalg.norm(Wr-Wfull) / max(np.linalg.norm(Wfull), 1e-12)
34    random_full_loss = np.sum((x @ Wr.T-y)**2)
35    return dict(n=n, d=d, r=r, m=m, support=len(idx), bound=(m+1)*r,
36      compression=n/len(idx), iterations=it, rank=rr, invariant_rel=float(inv_err),
37      solution_rel=float(sol_err), full_loss=float(full_loss), coreset_solution_full_loss=float(core_at_full),
38      full_loss_rel_error=float(abs(core_at_full-full_loss)/max(full_loss,1e-12)),
39      random_solution_rel=float(random_err), random_full_loss=float(random_full_loss),
40      build_sec=float(build_sec), core_fit_sec=float(core_fit_sec))
41
42def main():
43    sweep = [one(160,12,r,m) for r in [1,2,3,4] for m in [1,2,4]]
44    scaling = [one(n,12,3,2) for n in [40,80,160,320]]
45    allrows=sweep+scaling
46    max_inv=max(z['invariant_rel'] for z in allrows)
47    max_sol=max(z['solution_rel'] for z in allrows)
48    bound_ok=bool(all(z['support'] <= z['bound'] for z in allrows))
49    scaling_summary=[{'n':z['n'], 'support':z['support'], 'predicted_support':9,
50                      'compression':z['compression']} for z in scaling]
51    summary = {'prediction_checks': {
52      'P1_normal_equation_invariance': {'prediction':'relative error ~ machine precision', 'observed_max':max_inv, 'pass':bool(max_inv<1e-10)},
53      'P2_support_bound': {'prediction':'support <= (m+1)r', 'observed_all_pass':bound_ok,
54                           'observed_supports':[z['support'] for z in sweep]},
55      'P3_rank_fixed_scaling': {'prediction':'support remains 9 and compression scales as n/9', 'observed':scaling_summary,
56                                'compression_over_n':[z['compression']/z['n'] for z in scaling]},
57      'P4_solution_exactness': {'prediction':'full-data loss at coreset fit equals full optimum', 'observed_max_relative_parameter_error':max_sol,
58                                'observed_max_full_loss_relative_error':max(z['full_loss_rel_error'] for z in allrows), 'pass':bool(max_sol<1e-10)}
59    }, 'random_control_mean_parameter_error':float(np.mean([z['random_solution_rel'] for z in sweep])), 'sweep':sweep, 'scaling':scaling}
60    with open('results.json','w') as f: json.dump(summary,f,indent=2)
61    print(json.dumps(summary['prediction_checks'], indent=2))
62    print('random_control_mean_parameter_error', summary['random_control_mean_parameter_error'])
63
64if __name__=='__main__': main()