1import json
 2import time
 3from pathlib import Path
 4import numpy as np
 5
 6SEED = 2820
 7A, B, S1, S2 = 1.0, 0.8, 0.5, 2.0
 8# Under K(u)=integral exp(i omega u) Khat(omega)d omega,
 9# a exp(-s^2 u^2/2) has nonnegative spectral mass a and omega~N(0,s^2).
10C = A + B
11NEG_MASS = B / C
12
13
14def exact_kernel(u, bb=B):
15    u = np.asarray(u)
16    return A * np.exp(-0.5 * (S1 * u) ** 2) - bb * np.exp(-0.5 * (S2 * u) ** 2)
17
18
19def draw_freqs(m, rng, bb=B):
20    c = A + bb
21    neg = rng.random(m) < (bb / c) if c else np.zeros(m, dtype=bool)
22    omega = np.where(neg, rng.normal(0.0, S2, m), rng.normal(0.0, S1, m))
23    signs = np.where(neg, -1.0, 1.0)
24    return omega, signs, c
25
26
27def features(x, omega):
28    c, s = np.cos(np.outer(x, omega)), np.sin(np.outer(x, omega))
29    return np.stack((c, s), axis=2).reshape(len(x), -1)
30
31
32def rff_kernel(x, y, m, signed, seed):
33    omega, signs, c = draw_freqs(m, np.random.default_rng(seed))
34    zx, zy = features(x, omega), features(y, omega)
35    d = np.repeat(signs if signed else np.ones(m), 2)
36    return (c / m) * ((zx * d) @ zy.T)
37
38
39def factor_aggregate(x, v, m, signed, seed):
40    omega, signs, c = draw_freqs(m, np.random.default_rng(seed))
41    z = features(x, omega)
42    d = np.repeat(signs if signed else np.ones(m), 2)
43    return (c / m) * ((z * d) @ (z.T @ v))
44
45
46def slope(ms, errors):
47    return float(np.polyfit(np.log(np.asarray(ms)), np.log(np.asarray(errors)), 1)[0])
48
49
50def main():
51    rng = np.random.default_rng(SEED)
52    n_pairs, trials = 180, 100
53    u = rng.uniform(-3.0, 3.0, n_pairs)
54    truth = exact_kernel(u)
55    ms = [16, 32, 64, 128, 256, 512]
56    signed_rmse, positive_rmse, signed_bias = [], [], []
57    for m in ms:
58        se, pe, estimates = [], [], []
59        for t in range(trials):
60            omega, signs, c = draw_freqs(m, np.random.default_rng(SEED + 10000*m + t))
61            vals = np.cos(np.outer(u, omega))
62            est = (c / m) * (vals * signs).sum(axis=1)
63            pos = (c / m) * vals.sum(axis=1)
64            se.append(np.mean((est - truth) ** 2)); pe.append(np.mean((pos - truth) ** 2))
65            estimates.append(est)
66        signed_rmse.append(float(np.sqrt(np.mean(se))))
67        positive_rmse.append(float(np.sqrt(np.mean(pe))))
68        signed_bias.append(float(np.mean(np.concatenate(estimates) - np.tile(truth, trials))))
69    observed_slope = slope(ms, signed_rmse)
70    positive_limit = float(np.sqrt(np.mean((A*np.exp(-0.5*(S1*u)**2) + B*np.exp(-0.5*(S2*u)**2) - truth)**2)))
71
72    mass_rows = []
73    for bb in [0.0, 0.2, 0.5, 0.8, 1.2]:
74        signed_truth = exact_kernel(u, bb)
75        positive_kernel = A*np.exp(-0.5*(S1*u)**2) + bb*np.exp(-0.5*(S2*u)**2)
76        discrepancy = float(np.sqrt(np.mean((positive_kernel - signed_truth)**2)))
77        predicted = float(2.0 * bb * np.sqrt(np.mean(np.exp(-(S2*u)**2))))
78        mass_rows.append({'negative_amplitude': bb, 'negative_mass_fraction': bb/(A+bb),
79                          'observed_positive_limit_rmse': discrepancy, 'predicted_linear_rmse': predicted})
80
81    n, r, m = 220, 8, 128
82    x, v = rng.uniform(-3, 3, n), rng.normal(size=(n, r))
83    exact = exact_kernel(x[:, None] - x[None, :]) @ v
84    errs, times = [], {'exact': [], 'signed': []}
85    for t in range(20):
86        st = time.perf_counter(); out = factor_aggregate(x, v, m, True, SEED+50000+t); times['signed'].append(time.perf_counter()-st)
87        errs.append(float(np.linalg.norm(out-exact)/np.linalg.norm(exact)))
88        st = time.perf_counter(); _ = exact_kernel(x[:, None]-x[None, :]) @ v; times['exact'].append(time.perf_counter()-st)
89    result = {'kernel': 'exp(-0.5*(0.5u)^2) - 0.8 exp(-0.5*(2u)^2)', 'C': C, 'negative_mass_fraction': NEG_MASS,
90      'prediction_1_M_minus_half': {'M': ms, 'signed_rmse': signed_rmse, 'slope': observed_slope, 'predicted_slope': -0.5, 'signed_bias_at_M512': signed_bias[-1]},
91      'prediction_2_positive_features': {'positive_rmse': positive_rmse, 'predicted_asymptotic_rmse': positive_limit},
92      'prediction_3_negative_mass_sweep': mass_rows,
93      'aggregation': {'n': n, 'r': r, 'M': m, 'relative_errors': errs, 'median_exact_seconds': float(np.median(times['exact'])), 'median_signed_factor_seconds': float(np.median(times['signed']))}}
94    Path('results.json').write_text(json.dumps(result, indent=2)); print(json.dumps(result, indent=2))
95
96if __name__ == '__main__': main()