1import json
  2import math
  3from pathlib import Path
  4
  5import numpy as np
  6
  7
  8SEED = 2977
  9
 10
 11def gaussian_kernel(X, Y=None):
 12    if Y is None:
 13        Y = X
 14    d2 = ((X[:, None, :] - Y[None, :, :]) ** 2).sum(axis=2)
 15    return np.exp(-0.5 * d2)
 16
 17
 18def hermite_rule_2d(q):
 19    # Gauss-Hermite integrates exp(-t^2) f(t); t=sqrt(2) w gives
 20    # standard-normal integration, with weights normalized to sum to one.
 21    roots, weights = np.polynomial.hermite.hermgauss(q)
 22    nodes_1d = np.sqrt(2.0) * roots
 23    weights_1d = weights / np.sqrt(np.pi)
 24    nodes = np.array([(a, b) for a in nodes_1d for b in nodes_1d])
 25    weights_2d = np.array([wa * wb for wa in weights_1d for wb in weights_1d])
 26    weights_2d /= weights_2d.sum()
 27    return nodes, weights_2d
 28
 29
 30def quadrature_features(X, q):
 31    W, a = hermite_rule_2d(q)
 32    # G(w,x)=exp(i w^T x), Z_mi=sqrt(a_i)G(w_i,x_m).
 33    Z = np.exp(1j * (X @ W.T)) * np.sqrt(a)[None, :]
 34    return Z
 35
 36
 37def random_features(X, n, rng):
 38    W = rng.normal(size=(n, X.shape[1]))
 39    Z = np.exp(1j * (X @ W.T)) / np.sqrt(n)
 40    return Z
 41
 42
 43def eig_metrics(K, Kh, top=20):
 44    lam = np.linalg.eigvalsh(K)[::-1]
 45    lhat = np.linalg.eigvalsh((Kh + Kh.conj().T) / 2.0)[::-1]
 46    j = min(top, len(lam))
 47    top_abs = float(np.mean(np.abs(lam[:j] - lhat[:j])))
 48    top_rel = float(np.mean(np.abs(lam[:j] - lhat[:j]) / np.maximum(lam[:j], 1e-12)))
 49    op = float(np.linalg.norm(K - Kh, ord=2))
 50    # Weyl's inequality is checked for every ordered eigenvalue.
 51    weyl_max = float(np.max(np.abs(lam - lhat)))
 52    return {"top20_abs_eigen_error": top_abs, "top20_relative_eigen_error": top_rel,
 53            "operator_error": op, "weyl_bound_holds": bool(weyl_max <= op + 2e-8),
 54            "largest_exact_eigenvalue": float(lam[0]), "largest_approx_eigenvalue": float(lhat[0])}
 55
 56
 57def krr_mse(Xtr, ytr, Xte, yte, Ztr, Zte, ridge=1e-3):
 58    # Complex feature-space ridge regression; same representation for both methods.
 59    A = Ztr.conj().T @ Ztr + ridge * np.eye(Ztr.shape[1])
 60    theta = np.linalg.solve(A, Ztr.conj().T @ ytr)
 61    pred = np.real(Zte @ theta)
 62    return float(np.mean((pred - yte) ** 2))
 63
 64
 65def main():
 66    rng = np.random.default_rng(SEED)
 67    # Fixed points make all methods directly comparable. A smooth target is
 68    # deliberately compatible with the RBF kernel used by the Fourier integral.
 69    X = rng.normal(size=(120, 2))
 70    K = gaussian_kernel(X)
 71    rows = []
 72    for q in (3, 4, 6, 8):
 73        Zq = quadrature_features(X, q)
 74        mq = eig_metrics(K, Zq @ Zq.conj().T)
 75        random_runs = []
 76        for s in range(20):
 77            Zr = random_features(X, q * q, np.random.default_rng(SEED + 1000 * q + s))
 78            random_runs.append(eig_metrics(K, Zr @ Zr.conj().T))
 79        rr = {key: float(np.mean([r[key] for r in random_runs]))
 80              for key in ("top20_abs_eigen_error", "top20_relative_eigen_error", "operator_error",
 81                          "largest_approximate_eigenvalue") if key in random_runs[0]}
 82        rows.append({"features": q * q, "quadrature": mq, "random_mean": rr,
 83                     "random_std_top20_relative": float(np.std([r["top20_relative_eigen_error"] for r in random_runs]))})
 84
 85    # Equal-width learning check at 64 features, with a held-out split.
 86    Xtr, Xte = X[:80], X[80:]
 87    ytr = np.sin(Xtr[:, 0]) + 0.25 * np.cos(1.5 * Xtr[:, 1])
 88    yte = np.sin(Xte[:, 0]) + 0.25 * np.cos(1.5 * Xte[:, 1])
 89    Zqtr, Zqte = quadrature_features(Xtr, 8), quadrature_features(Xte, 8)
 90    q_mse = krr_mse(Xtr, ytr, Xte, yte, Zqtr, Zqte)
 91    random_mses = []
 92    for s in range(20):
 93        # Generate the same frequencies for train and test in each replicate.
 94        W = np.random.default_rng(SEED + 7000 + s).normal(size=(64, 2))
 95        Ztr = np.exp(1j * (Xtr @ W.T)) / 8.0
 96        Zte = np.exp(1j * (Xte @ W.T)) / 8.0
 97        random_mses.append(krr_mse(Xtr, ytr, Xte, yte, Ztr, Zte))
 98
 99    # Independent cheap algebra sanity check of Z Z* and the probability rule.
100    W, a = hermite_rule_2d(4)
101    Z = np.exp(1j * (X[:12] @ W.T)) * np.sqrt(a)[None, :]
102    # Independent matrix-form evaluation of the same quadrature sum.
103    direct = np.exp(1j * (X[:12] @ W.T)) @ np.diag(a) @ np.exp(-1j * (W @ X[:12].T))
104    algebra_error = float(np.max(np.abs(Z @ Z.conj().T - direct)))
105
106    result = {"seed": SEED, "kernel": "exp(-||x-y||^2/2)", "rows": rows,
107              "learning_64_features": {"quadrature_mse": q_mse,
108                                        "random_mean_mse": float(np.mean(random_mses)),
109                                        "random_std_mse": float(np.std(random_mses))},
110              "sanity": {"weights_sum": float(a.sum()), "weights_nonnegative": bool(np.all(a >= 0)),
111                         "ZZh_identity_max_error": algebra_error}}
112    Path("results.json").write_text(json.dumps(result, indent=2))
113    print(json.dumps(result, indent=2))
114
115
116if __name__ == "__main__":
117    main()