"""Brute-force checks of Theorem 3.1(b) (posterior predictive from two error profiles), Theorem 5.4 (local Fourier features), Proposition 5.1 (halfspace forms) and Theorem 3.3 (Bayes risk, monotonicity). Run: python -I check_inference_linear_risk.py """ import itertools import random from fractions import Fraction from math import comb, log2 random.seed(7) # ---------------- random local classes ---------------- def random_local_class(N, kappa, nx): """Domain of nx inputs; input x has a set R(x) of <= kappa bits and a random predicate f_x.""" cls = [] for _ in range(nx): r = random.randint(0, min(kappa, N)) R = tuple(sorted(random.sample(range(N), r))) table = {beta: random.randint(0, 1) for beta in itertools.product((0, 1), repeat=r)} cls.append((R, table)) return cls def label(cls, w, x): R, table = cls[x] return table[tuple(w[i] for i in R)] def profile(cls, N, S): """Error profile coefficients c_j, j = 0..len(S).""" c = [0] * (len(S) + 1) for w in itertools.product((0, 1), repeat=N): e = sum(label(cls, w, x) != y for x, y in S) c[e] += 1 return c def check_profile_recursion(trials=200): for _ in range(trials): N = random.randint(2, 7) kappa = random.randint(1, 3) cls = random_local_class(N, kappa, 6) m = random.randint(1, 6) S = [(random.randrange(6), random.randint(0, 1)) for _ in range(m)] x = random.randrange(6) eta = Fraction(random.randint(1, 9), 20) # in (0, 1/2) rho = eta / (1 - eta) c = profile(cls, N, S) cp = profile(cls, N, S + [(x, 1)]) # recursion of step <2>2 a, b = [0] * (m + 1), [0] * (m + 1) prev_b = 0 for j in range(m + 1): a[j] = cp[j] - prev_b b[j] = c[j] - a[j] prev_b = b[j] assert cp[m + 1] == b[m] # direct counts for j in range(m + 1): aj = bj = 0 for w in itertools.product((0, 1), repeat=N): if sum(label(cls, w, xx) != y for xx, y in S) == j: if label(cls, w, x) == 1: aj += 1 else: bj += 1 assert (aj, bj) == (a[j], b[j]) p_formula = sum(a[j] * rho ** j for j in range(m + 1)) / sum(c[j] * rho ** j for j in range(m + 1)) num = den = Fraction(0) for w in itertools.product((0, 1), repeat=N): wt = rho ** sum(label(cls, w, xx) != y for xx, y in S) den += wt num += wt * label(cls, w, x) assert p_formula == num / den print("Posterior predictive from two error profiles (Theorem 3.1(b)): OK on", trials, "instances") # ---------------- local Fourier features ---------------- def chi(T, bits): return (-1) ** sum(bits[i] for i in T) def fourier_features(R, table): """psi_T for T subset of R (as tuples of global indices), by the definition.""" r = len(R) feats = {} for size in range(r + 1): for T in itertools.combinations(R, size): pos = [R.index(i) for i in T] s = sum(table[beta] * (-1) ** sum(beta[p] for p in pos) for beta in itertools.product((0, 1), repeat=r)) feats[T] = Fraction(s, 2 ** r) return feats def check_linearization(trials=200): for _ in range(trials): N = random.randint(1, 8) kappa = random.randint(1, min(3, N)) cls = random_local_class(N, kappa, 5) for R, table in cls: feats = fourier_features(R, table) for f in feats.values(): assert (f * 2 ** kappa).denominator == 1 assert len(feats) <= 2 ** len(R) for w in itertools.product((0, 1), repeat=N): val = sum(f * chi(T, w) for T, f in feats.items()) assert val == table[tuple(w[i] for i in R)] # bound D_kappa <= (eN/kappa)^kappa D = sum(comb(N, j) for j in range(kappa + 1)) assert D <= (2.718281828459045 * N / kappa) ** kappa + 1e-9 print("Local Fourier features (Theorem 5.4): OK on", trials, "random local classes") # ---------------- halfspace forms ---------------- def check_halfspaces(): for n in range(1, 6): for x in itertools.product((-1, 0, 1), repeat=n): r = sum(1 for v in x if v != 0) for sigma in itertools.product((-1, 1), repeat=n): ip = sum(s * v for s, v in zip(sigma, x)) ld = int(any(v != 0 and v == s for s, v in zip(sigma, x))) assert ld == int(ip + r - 2 >= 0) if ld == 0: assert ip + r - 2 == -2 bp = int(ip >= 0) if bp == 0: assert ip <= -1 print("Halfspace forms (Proposition 5.1): OK for n <= 5") # ---------------- Bayes risk ---------------- def bayes_risk(cls, N, D, eta, m): """Exact Bayes risk R*_m for the uniform prior over {0,1}^N, input distribution D (list of probs).""" nx = len(D) W = list(itertools.product((0, 1), repeat=N)) labels = [[label(cls, w, x) for x in range(nx)] for w in W] total = Fraction(0) for xs in itertools.product(range(nx), repeat=m + 1): px = Fraction(1) for x in xs: px *= D[x] if px == 0: continue q = xs[-1] for ys in itertools.product((0, 1), repeat=m): # joint probability of labels and posterior predictive of target 1 jw = [] for wi in range(len(W)): pr = Fraction(1, len(W)) for x, y in zip(xs[:-1], ys): pr *= (1 - eta) if labels[wi][x] == y else eta jw.append(pr) z = sum(jw) if z == 0: continue p1 = sum(jw[wi] * ((1 - eta) if labels[wi][q] == 1 else eta) for wi in range(len(W))) total += px * min(p1, z - p1) return total def binary_entropy(e): e = float(e) return 0.0 if e == 0 else -e * log2(e) - (1 - e) * log2(1 - e) def check_risk(trials=12): for t in range(trials): N = random.randint(2, 3) cls = random_local_class(N, 2, 4) weights = [random.randint(1, 4) for _ in range(4)] D = [Fraction(wt, sum(weights)) for wt in weights] eta = Fraction(0) if t % 2 == 0 else Fraction(random.choice([1, 2, 3]), 10) risks = [bayes_risk(cls, N, D, eta, m) for m in range(0, 5)] for m in range(len(risks) - 1): assert risks[m + 1] <= risks[m], (risks, eta) for m, r in enumerate(risks): bound = float(eta) + (1 - 2 * float(eta)) * N / ((1 - binary_entropy(eta)) * (m + 1)) assert float(r) <= bound + 1e-12, (m, float(r), bound) print("Bayes risk (Theorem 3.3): monotone and below the bound on", trials, "random task families, m <= 4") if __name__ == "__main__": check_profile_recursion() check_linearization() check_halfspaces() check_risk()