Markov Spectral Equivariant Layer / markov_filter_experiment.py
Mechanism confirmed, baseline not beaten
1#!/usr/bin/env python3
2"""MVP verification of a positive finite-rank Markov spectral layer on Z_M.
3
4The cyclic grid is a finite compact-group analogue of SO(2). DFT convolution
5is exactly equivariant to cyclic shifts. The Fejer kernel is constructed as a
6squared Dirichlet kernel, so positivity and unit mass are not assumptions.
7"""
8import json
9import math
10from pathlib import Path
11import numpy as np
12
13SEED = 1203
14rng = np.random.default_rng(SEED)
15M = 256
16
17
18def fejer_kernel(M, L):
19 """Nonnegative unit-mass kernel with 2L-1 retained Fourier modes."""
20 j = np.arange(M)
21 # squared finite geometric sum; works stably at j=0
22 z = np.exp(2j * np.pi * j / M)
23 d = np.zeros(M, dtype=np.complex128)
24 for r in range(L):
25 d += z ** r
26 k = (np.abs(d) ** 2) / (M * L)
27 return k.real
28
29
30def centered_mask(M, L):
31 freqs = np.fft.fftfreq(M) * M
32 return (np.abs(freqs) < L).astype(float)
33
34
35def signed_kernel(M, L):
36 return np.fft.ifft(centered_mask(M, L)).real
37
38
39def conv(x, k):
40 return np.fft.ifft(np.fft.fft(x) * np.fft.fft(k)).real
41
42
43def shift(x, s):
44 return np.roll(x, s)
45
46
47def induced_inf_norm(k):
48 # For a circulant convolution, this is exactly its l_infinity operator norm.
49 return float(np.sum(np.abs(k)))
50
51
52def spectral_coefficients(k):
53 return np.fft.fft(k).real
54
55
56def run():
57 Ls = [2, 4, 8, 16, 32, 64]
58 kernel_rows = []
59 attenuation_rows = []
60 random_rows = []
61 equivariance_errors = []
62 order_errors = []
63
64 # Prediction 1: positivity + mass imply Markov properties for every L.
65 # Prediction 2: the induced sup norm is exactly one for Fejer, while signed
66 # truncation has a growing norm (and can amplify a worst-case input).
67 for L in Ls:
68 kf = fejer_kernel(M, L)
69 ks = signed_kernel(M, L)
70 af = spectral_coefficients(kf)
71 raw_norm = induced_inf_norm(ks)
72 fejer_norm = induced_inf_norm(kf)
73 worst = np.sign(np.roll(ks, np.argmax(np.abs(ks))))
74 # A direct worst-case witness for the first convolution row.
75 witness_ratio = np.max(np.abs(conv(worst, ks))) / np.max(np.abs(worst))
76 f = rng.uniform(-1, 1, M)
77 q = rng.uniform(-1, 1, M)
78 markov_ratio = np.max(np.abs(conv(f, kf) - conv(q, kf))) / np.max(np.abs(f-q))
79 # Deterministic order test q=f+r, r>=0: positivity predicts F(q)-F(f)>=0.
80 f0 = rng.normal(size=M)
81 r0 = rng.random(M)
82 order_errors.append(float(max(0.0, -np.min(conv(f0 + r0, kf) - conv(f0, kf)))))
83 kernel_rows.append({
84 "L": L, "min_fejer_kernel": float(kf.min()),
85 "mass_error": float(abs(kf.sum() - 1)),
86 "fejer_inf_norm": fejer_norm, "signed_inf_norm": raw_norm,
87 "signed_worst_case_ratio": float(witness_ratio),
88 "fejer_random_ratio": float(markov_ratio),
89 })
90
91 # Prediction 3: Fejer eigenvalue at frequency k is max(1-|k|/L,0).
92 checks = []
93 for k in [0, max(1, L // 2), L - 1, L, min(M // 2, L + 3)]:
94 expected = max(1.0 - abs(k) / L, 0.0)
95 actual = af[k % M]
96 checks.append({"k": k, "observed": float(actual), "predicted": expected,
97 "abs_error": float(abs(actual - expected))})
98 attenuation_rows.append({"L": L, "checks": checks})
99
100 # Exact equivariance test: F(T_s x) = T_s F(x).
101 x = rng.normal(size=M)
102 s = int(rng.integers(0, M))
103 equivariance_errors.append(float(np.max(np.abs(conv(shift(x, s), kf) - shift(conv(x, kf), s)))))
104
105 # Average bounded random-input behavior, plus 20 repeated layers.
106 vals = []
107 vals_s = []
108 for _ in range(100):
109 x = rng.uniform(-1, 1, M)
110 vals.append(np.max(np.abs(conv(x, kf))) / np.max(np.abs(x)))
111 vals_s.append(np.max(np.abs(conv(x, ks))) / np.max(np.abs(x)))
112 xf = rng.uniform(-1, 1, M)
113 xs = xf.copy()
114 for _ in range(20):
115 xf, xs = conv(xf, kf), conv(xs, ks)
116 random_rows.append({"L": L, "mean_fejer_ratio": float(np.mean(vals)),
117 "mean_signed_ratio": float(np.mean(vals_s)),
118 "20_layer_signed_norm_ratio": float(np.max(np.abs(xs)) / np.max(np.abs(xf)) if np.max(np.abs(xf)) else 0.0)})
119
120 # Correct, reproducible stacked-layer comparison with a common initial signal.
121 stack_rows = []
122 for L in Ls:
123 kf, ks = fejer_kernel(M, L), signed_kernel(M, L)
124 x0 = rng.uniform(-1, 1, M)
125 xf, xs = x0.copy(), x0.copy()
126 for _ in range(20):
127 xf, xs = conv(xf, kf), conv(xs, ks)
128 stack_rows.append({"L": L, "fejer_20_layer_sup_ratio": float(np.max(abs(xf))/np.max(abs(x0))),
129 "signed_20_layer_sup_ratio": float(np.max(abs(xs))/np.max(abs(x0)))})
130
131 result = {
132 "seed": SEED, "M": M,
133 "predictions": [
134 "For every L, min(K)>=0 and sum(K)=1, hence induced sup norm and random error ratios are <=1.",
135 "Matched-rank signed truncation has induced sup norm sum(abs(K)); its worst-case ratio equals that value and generally grows with rank.",
136 "Fejer Fourier eigenvalue at integer frequency k is max(1-|k|/L,0), giving linear attenuation to the cutoff and zero beyond it."
137 ],
138 "kernel_sweep": kernel_rows,
139 "attenuation_sweep": attenuation_rows,
140 "equivariance_max_errors": equivariance_errors,
141 "order_violation_max": max(order_errors),
142 "random_and_stack": random_rows,
143 "stacked_20_layers": stack_rows,
144 "summary": {
145 "max_fejer_mass_error": max(r["mass_error"] for r in kernel_rows),
146 "min_fejer_kernel": min(r["min_fejer_kernel"] for r in kernel_rows),
147 "max_equivariance_error": max(equivariance_errors),
148 "max_attenuation_error": max(c["abs_error"] for row in attenuation_rows for c in row["checks"]),
149 "signed_norms": [r["signed_inf_norm"] for r in kernel_rows],
150 "fejer_norms": [r["fejer_inf_norm"] for r in kernel_rows]
151 }
152 }
153 Path("results.json").write_text(json.dumps(result, indent=2))
154 print(json.dumps(result["summary"], indent=2))
155 print("kernel_sweep")
156 for r in kernel_rows:
157 print(r)
158 print("stacked_20_layers")
159 for r in stack_rows:
160 print(r)
161
162if __name__ == "__main__":
163 run()