"""Checks Proposition 6.4: the gap formula of part (a) and the reduction of part (b). Part (b): a 3-CNF on n0 variables becomes a Set Cover instance with Y = 2*n0 sets and t = n0, and then samples in R^(Y+2) with E(z) = sum ([]_+ - y)^2. Checked with exact rationals: the formula is satisfiable iff at most t sets cover; a cover gives E(z) <= beta^2 t with ||z||^2 <= Y+2; the integer normalization keeps samples and weight in the unit ball and labels in [0,1], and divides the error by Lambda^2; and 1/eps = O(n0^9). Step <2>3 is tested numerically by gradient descent from cover-like and random starts: a point with E(z) < beta^2 (t+1) must give a family of at most t sets that covers the universe, so on unsatisfiable formulas no such point may appear. This is evidence for that step, which the paper proves, and not a proof. Part (a): the three-term form of the gap of [GKMR21] and the bounds on its terms for k >= lambda N. Run: python -I check_relu_setcover.py """ import itertools import math import random from fractions import Fraction as Fr random.seed(20261030) def ceil_sqrt(a): r = math.isqrt(a) return r if r * r == a else r + 1 def literals(n0): return [(v, s) for v in range(n0) for s in (True, False)] def random_cnf(n0, m0): """A clause is a frozenset of literals of width at most 3, without complementary literals.""" clauses = set() for _ in range(4 * m0): if len(clauses) == m0: break vs = random.sample(range(n0), random.randint(1, min(3, n0))) clauses.add(frozenset((v, random.random() < 0.5) for v in vs)) return sorted(clauses, key=sorted) def satisfiable(n0, clauses): for a in itertools.product((True, False), repeat=n0): if all(any(a[v] == s for v, s in c) for c in clauses): return True return False def set_cover(n0, clauses): """Universe: the variables 0..n0-1, then the clauses. One set for each literal.""" sets = [{v} | {n0 + i for i, c in enumerate(clauses) if (v, s) in c} for v, s in literals(n0)] return sets, n0 + len(clauses) def find_cover(sets, size_w, t): for r in range(t + 1): for fam in itertools.combinations(range(len(sets)), r): if len(set().union(*[sets[i] for i in fam])) == size_w: return fam return None def samples(sets, size_w, beta): """Coordinates: 0..Y-1 for the sets, Y for z_1, Y+1 for z_*. Vectors as lists of integers.""" Y = len(sets) S = [] for w in range(size_w): x = [1 if w in T else 0 for T in sets] + [1, 0] S.append((x, Fr(0))) for i in range(Y): x = [0] * (Y + 2) x[i] = x[Y + 1] = 1 S.append((x, beta)) S.append(([0] * Y + [1, 0], Fr(1))) S.append(([0] * Y + [0, 1], beta)) return S def energy(z, S): return sum((max(sum(a * b for a, b in zip(z, x)), 0) - y) ** 2 for x, y in S) def check_exact(trials=300): sat = unsat = 0 for _ in range(trials): n0 = random.randint(1, 4) clauses = random_cnf(n0, random.randint(0, 6 * n0)) assert len(clauses) < 8 * n0 ** 3 sets, size_w = set_cover(n0, clauses) Y, t = 2 * n0, n0 fam = find_cover(sets, size_w, t) assert (fam is not None) == satisfiable(n0, clauses) beta = Fr(1, 100 * Y * Y) S = samples(sets, size_w, beta) assert len(S) == 3 * n0 + len(clauses) + 2 assert all(sum(a * a for a in x) <= Y + 1 for x, _ in S) L1, L2 = ceil_sqrt(Y + 1), ceil_sqrt(Y + 2) L = L1 * L2 if fam is None: unsat += 1 continue sat += 1 z = [Fr(0)] * (Y + 2) for i in fam: z[i] = Fr(-1) z[Y], z[Y + 1] = Fr(1), beta E = energy(z, S) assert E <= beta * beta * t assert sum(a * a for a in z) <= Y + 2 <= L2 * L2 Sn = [([Fr(a, L1) for a in x], y / L) for x, y in S] zn = [a / L2 for a in z] assert all(sum(a * a for a in x) <= 1 and 0 <= y <= 1 for x, y in Sn) assert sum(a * a for a in zn) <= 1 assert energy(zn, Sn) == E / (L * L) eps = beta ** 2 / (2 * len(S) * L * L) assert E / (len(S) * L * L) + eps < beta ** 2 * (t + 1) / (len(S) * L * L) assert sat and unsat print(f"Exact checks of the reduction and the normalization on {sat} satisfiable and {unsat} unsatisfiable formulas") def check_count(): worst = 0.0 for n0 in (1, 2, 3, 5, 10, 100, 1000, 10000): Y = 2 * n0 size_s = 3 * n0 + (8 * n0 ** 3 - 1) + 2 L = ceil_sqrt(Y + 1) * ceil_sqrt(Y + 2) assert L * L <= 16 * (Y + 1) * (Y + 2) worst = max(worst, 2 * size_s * L * L * (100 * Y * Y) ** 2 / n0 ** 9) assert worst < 1e8 print(f"1/eps <= {worst:.3g} * n0^9 for the tested n0, with M0 = 8 n0^3 - 1 clauses") def descend(z, sparse, Y, t, beta, sets, size_w, lr, steps): """Gradient descent on E. Returns the smallest E seen and the number of points below the threshold.""" thr = beta * beta * (t + 1) * (1 - 1e-9) best, hits = math.inf, 0 for _ in range(steps): grad = [0.0] * len(z) E = 0.0 for idx, y in sparse: p = sum(z[i] for i in idx) r = max(p, 0.0) - y E += r * r if p > 0: for i in idx: grad[i] += 2 * r best = min(best, E) if E < thr: hits += 1 fam = [i for i in range(Y) if z[i] + z[Y + 1] <= 0] assert len(fam) <= t assert len(set().union(*[sets[i] for i in fam])) == size_w z = [a - lr * g for a, g in zip(z, grad)] return best, hits def check_soundness(formulas=24, starts=3, steps=600): low_unsat, sat_hits, n_unsat, n_sat = math.inf, 0, 0, 0 for _ in range(formulas): n0 = random.randint(1, 3) clauses = random_cnf(n0, random.randint(1, 6 * n0)) sets, size_w = set_cover(n0, clauses) Y, t = 2 * n0, n0 beta = 0.01 / (Y * Y) S = samples(sets, size_w, Fr(1, 100 * Y * Y)) sparse = [([i for i, a in enumerate(x) if a], float(y)) for x, y in S] lr = 1 / (2 * sum(len(idx) for idx, _ in sparse)) fam = find_cover(sets, size_w, t) for s in range(starts): if s == 0: chosen = fam if fam is not None else random.sample(range(Y), min(t + 1, Y)) z = [-1.0 if i in chosen else 0.0 for i in range(Y)] + [1.0, beta] z = [a + random.gauss(0, 1e-3 * beta) for a in z] else: z = [random.uniform(-1, 0.2) for _ in range(Y)] + [random.uniform(0, 1), random.uniform(0, 0.1)] best, hits = descend(z, sparse, Y, t, beta, sets, size_w, lr, steps) if fam is None: assert hits == 0 low_unsat = min(low_unsat, best / (beta * beta * (t + 1))) else: sat_hits += hits n_sat += fam is not None n_unsat += fam is None assert n_unsat and sat_hits print(f"Descent: {n_unsat} unsatisfiable formulas, smallest E / (beta^2 (t+1)) = {low_unsat:.4f}; " f"{sat_hits} points below the threshold on {n_sat} satisfiable formulas all gave covers (numeric evidence)") def check_gap_formula(): rt2 = math.sqrt(2) count = 0 for d in (1, 3, 10): gamma, zeta = 1 / (1000 * d), 1 / (1e10 * d * d) for lam in (1.0, 0.1, 0.001): n_min = math.ceil(5e6 * d ** 4 / lam) for N in (n_min, 10 * n_min, 1000 * n_min): k = max(2, math.ceil(lam * N)) lo_first = gamma * zeta * math.sqrt(lam) / (12 * rt2 * d * math.sqrt(N)) hi_first = gamma * zeta / (8 * math.sqrt(2 * lam * N)) for M in {math.ceil(k / 3), d * N // 2}: for xi in {math.ceil(k / 3), M}: if not k / 3 <= xi <= M <= d * N / 2: continue delta = 1 / (2 * math.sqrt(k)) eps = (gamma * zeta * xi * delta / (4 * rt2 * M) - (1 - gamma) * zeta * k * delta ** 2 / (8 * N) - gamma * zeta * xi * delta ** 2 / (32 * M)) first = gamma * zeta * xi / (8 * math.sqrt(2 * k) * M) second = (1 - gamma) * zeta / (32 * N) third = gamma * zeta * xi / (128 * k * M) assert math.isclose(eps, first - second - third, rel_tol=1e-9) assert lo_first * (1 - 1e-9) <= first <= hi_first * (1 + 1e-9) assert third <= first / 2 assert lo_first / 4 * (1 - 1e-9) <= eps <= hi_first * (1 + 1e-9) count += 1 print(f"Part (a): gap formula and Theta(1/sqrt N) bounds checked on {count} parameter settings, for N >= 5e6 d^4 / lambda") if __name__ == "__main__": check_gap_formula() check_exact() check_count() check_soundness()