1import json, math, time
  2import numpy as np
  3from scipy.linalg import eigh
  4
  5SEED = 1226
  6
  7
  8def unitary(A, t):
  9    w, V = eigh(A)
 10    return (V * np.exp(-1j * t * w)) @ V.T
 11
 12
 13def dyadic_exact(A, lam, K):
 14    n = A.shape[0]
 15    I = np.eye(n, dtype=complex)
 16    R = 1j * np.linalg.inv(I - math.exp(-lam) * unitary(A, 1.0))
 17    for k in range(1, K + 1):
 18        t = 2.0 ** (-k)
 19        R -= 1j * t * np.linalg.inv(I + math.exp(-lam * t) * unitary(A, t))
 20    return R
 21
 22
 23def dyadic_neumann(A, X, lam, K, J):
 24    Y = np.zeros_like(X, dtype=complex)
 25    U = unitary(A, 1.0)
 26    Z = X.astype(complex)
 27    for j in range(J + 1):
 28        Y += 1j * math.exp(-j * lam) * Z
 29        Z = U @ Z
 30    for k in range(1, K + 1):
 31        t = 2.0 ** (-k)
 32        c = math.exp(-lam * t)
 33        U = unitary(A, t)
 34        Z = X.astype(complex)
 35        for j in range(J + 1):
 36            Y -= 1j * t * ((-c) ** j) * Z
 37            Z = U @ Z
 38    return Y
 39
 40
 41def scalar_verification():
 42    vals = np.linspace(-2.0, 2.0, 17)
 43    A = np.diag(vals)
 44    out = {}
 45    for lam in [0.2, 0.5, 1.0, 2.0]:
 46        target = np.diag(1.0 / (vals - 1j * lam))
 47        errs = []
 48        for K in [2, 4, 6, 8, 10, 12]:
 49            errs.append(float(np.max(np.abs(dyadic_exact(A, lam, K) - target))))
 50        # Exact resolvent tail is bounded by sum_{k>K} 2^-k = 2^-K,
 51        # since ||(I+cU)^-1|| <= 1/(1-c) is not uniformly useful; empirically
 52        # the measured tail should decrease with K and remain below a constant*2^-K.
 53        ratios = [errs[i] / errs[i + 1] for i in range(len(errs)-1)]
 54        # Neumann prediction: at scale t, remainder is at most c^(J+1)/(1-c).
 55        nerrs = []
 56        for J in [1, 2, 4, 8, 12, 16]:
 57            X = np.eye(len(vals))
 58            nerrs.append(float(np.max(np.abs(dyadic_neumann(A, X, lam, 10, J) - dyadic_exact(A, lam, 10)))))
 59        # Unitary norm prediction over random vectors and times.
 60        x = np.random.default_rng(7).normal(size=len(vals))
 61        norm_errors = [abs(np.linalg.norm(unitary(A, t) @ x) / np.linalg.norm(x) - 1) for t in [1, .5, .125, .01]]
 62        out[str(lam)] = {"K": [2,4,6,8,10,12], "exact_errors": errs,
 63                         "K_error_ratios": ratios, "J": [1,2,4,8,12,16],
 64                         "neumann_errors": nerrs, "unitary_relative_norm_errors": norm_errors}
 65    return out
 66
 67
 68def graph_data(n=80, p_in=.16, p_out=.035, d=8):
 69    rng = np.random.default_rng(SEED)
 70    y = np.repeat([0, 1], n // 2)
 71    G = np.zeros((n, n), float)
 72    for i in range(n):
 73        for j in range(i+1, n):
 74            p = p_in if y[i] == y[j] else p_out
 75            if rng.random() < p: G[i,j] = G[j,i] = 1
 76    deg = G.sum(1)
 77    L = np.eye(n) - (G / np.sqrt(np.maximum(deg, 1)[:,None] * np.maximum(deg, 1)[None,:]))
 78    # Features include weak signal; held-out nodes use same fixed feature matrix.
 79    X = rng.normal(0, 1, (n, d))
 80    X[:, 0] += (2*y-1) * .65
 81    perm = rng.permutation(n)
 82    tr, va = perm[:n//2], perm[n//2:]
 83    return L, X, y, tr, va
 84
 85
 86def train_models():
 87    # Tiny ridge classifier compares standard low-degree polynomial filtering to dyadic filtering.
 88    L, X, y, tr, va = graph_data()
 89    # Rescale normalized Laplacian to [-1,1] (spectrum is approximately [0,2]).
 90    A = L - np.eye(len(L))
 91    rng = np.random.default_rng(SEED)
 92    results = {}
 93    configs = [("baseline_poly", 6), ("dyadic_K2_J2", 6), ("dyadic_K4_J2", 6)]
 94    for name, budget in configs:
 95        t0 = time.perf_counter()
 96        if name == "baseline_poly":
 97            Z = X.astype(complex); F = [Z]
 98            U = A @ Z
 99            for _ in range(budget):
100                F.append(U); U = A @ U
101            H = np.concatenate([z.real for z in F], axis=1)
102        else:
103            K = int(name.split('K')[1].split('_')[0]); J = int(name.split('J')[1])
104            Hc = dyadic_neumann(A, X, .8, K, J)
105            H = np.concatenate([X, Hc.real, Hc.imag], axis=1)
106        # closed-form ridge, no training randomness
107        lam = 1.0
108        W = np.linalg.solve(H[tr].T @ H[tr] + lam*np.eye(H.shape[1]), H[tr].T @ (2*y[tr]-1))
109        pred = (H @ W > 0).astype(int)
110        acc = float(np.mean(pred[va] == y[va]))
111        norm = float(np.linalg.norm(H) / np.linalg.norm(X))
112        results[name] = {"val_accuracy": acc, "feature_norm_ratio": norm,
113                         "seconds": time.perf_counter()-t0, "sparse_op_budget": budget}
114    return results
115
116
117def main():
118    np.random.seed(SEED)
119    report = {"seed": SEED, "verification": scalar_verification(), "mini_experiment": train_models()}
120    with open("results.json", "w") as f: json.dump(report, f, indent=2)
121    print(json.dumps(report, indent=2))
122
123if __name__ == "__main__": main()