From 5a5cfc723decc515c30eba625db42bc1e9b9e642 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Sat, 14 Nov 2026 17:49:44 -0700 Subject: [PATCH] Heaviest-term lower bound on C_g with the heavy term restricted to generator states heaviest_term_bound takes a weight_cap per box; generator_heaviest_term_bound supplies the one-explicit-state generator-cap LP. Opt-in pipeline stage (generator_heaviest_term=True). Co-authored-by: Cursor --- CHANGELOG.md | 5 ++ docs/generators/generative_models.rst | 11 +++++ sofic/generators/exact_ci_bounds.py | 68 +++++++++++++++++++++++---- sofic/generators/generative_search.py | 27 ++++++++++- tests/test_exact_ci_bounds.py | 50 ++++++++++++++++++++ 5 files changed, 152 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ce47d10..857ae6e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -54,6 +54,11 @@ brute-force reference or a closed-form value and has a regression test. joints `((0,1),(2,·),(0,1))`, `((0,1),(0,2),(0,1))` and `((0,1),(2,1),(0,·))` it certifies `G` to within 1e-5 of the estimate (0.9335, 0.9243, 0.4717; was 0.9118, 0.9056, 0.4436), which closes the first bracket. +- `generator_heaviest_term_bound` bounds `C_g` the same way with the heavy term restricted to generator + states: `heaviest_term_bound(weight_cap=...)` caps each box's weight by the one-explicit-state + generator-cap LP over the box of beliefs, and prunes boxes the LP rules out. On `((0,1),(0,2),(0,1))` + it certifies `C_g >= 0.928 > G = 0.924`, but only after ~220k boxes, so the pipeline runs it on + request (`generator_heaviest_term=True`). ### Generator caps diff --git a/docs/generators/generative_models.rst b/docs/generators/generative_models.rst index 44df4d3..71845c2 100644 --- a/docs/generators/generative_models.rst +++ b/docs/generators/generative_models.rst @@ -301,6 +301,16 @@ the certified lower bound is still below the exact-common-information estimate a :math:`|S^+| + |S^-| \le 7`. On three survey processes with :math:`3 \times 3` joints it raises the lower bound to the estimate, closing ``((0,1),(2,·),(0,1))`` at :math:`0.9335` bits. +A generator's states are an exact decomposition whose heaviest term is its heaviest state, and that +state satisfies the generator equations, so the same branching bounds :math:`C_g` when each box's +weight is also capped by the one-explicit-state generator-cap LP over the box of beliefs +(:func:`~sofic.generators.exact_ci_bounds.generator_heaviest_term_bound`). Where the LP excludes the +:math:`\op{G}{}`-optimal heavy term the bound exceeds :math:`\op{G}{S^+ : S^-}`: on +``((0,1),(0,2),(0,1))`` it certifies :math:`C_g \ge 0.928 > \op{G}{} = 0.924`, but only after about +220,000 boxes (half an hour), because the LP excludes that term only in boxes of radius about +:math:`10^{-4}`. The pipeline therefore runs it only on request (``generator_heaviest_term=True``). +This combination is not from the literature. + :func:`~sofic.generators.generator_caps.certified_generator_weight_cap` removes the LP tolerance: every branch-and-bound leaf gets an exact bound from rationalized HiGHS duals or an exact Farkas certificate of infeasibility (:mod:`sofic.generators.certificates`, after @@ -357,6 +367,7 @@ Coupled generators and bounds .. autofunction:: sofic.generators.generator_caps.generator_weight_cap .. autofunction:: sofic.generators.generator_caps.generator_coverage_bound .. autofunction:: sofic.generators.exact_ci_bounds.heaviest_term_bound +.. autofunction:: sofic.generators.exact_ci_bounds.generator_heaviest_term_bound .. autofunction:: sofic.generators.exact_ci_bounds.capped_coverage_entropy .. autoclass:: sofic.generators.exact_ci_bounds.HeaviestTermBoundResult .. autofunction:: sofic.generators.generator_caps.certified_generator_weight_cap diff --git a/sofic/generators/exact_ci_bounds.py b/sofic/generators/exact_ci_bounds.py index d4a2a6d..a4844b2 100644 --- a/sofic/generators/exact_ci_bounds.py +++ b/sofic/generators/exact_ci_bounds.py @@ -21,13 +21,18 @@ import heapq import itertools import time +from collections.abc import Callable from dataclasses import dataclass +from typing import TYPE_CHECKING import numpy as np from sofic.generators.coverage_bound import mu1_vertex, mu_k_spectral_cap from sofic.generators.generator_caps import _tighten +if TYPE_CHECKING: + from sofic.generators.coupled_generators import _GeneratorProblem + Box = tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, float, float] @@ -121,6 +126,7 @@ def heaviest_term_bound( target: float | None = None, max_nodes: int = 200_000, deadline: float | None = None, + weight_cap: Callable[[Box], float | None] | None = None, ) -> HeaviestTermBoundResult: """Certified lower bound on :math:`\\op{G}{X : Y}` for the joint ``P``; see the module docstring. @@ -128,6 +134,10 @@ def heaviest_term_bound( is at least ``target``, after ``max_nodes`` boxes, or at ``deadline`` (a :func:`time.monotonic` value). The caps are :func:`~sofic.generators.coverage_bound.mu1_vertex` and the spectral caps, as in the certified coverage bound. The work grows quickly with ``|X| + |Y|``; it is meant for small joints. + + ``weight_cap(box)``, if given, is an extra upper bound on the heavy term's weight over a box (``None`` + when no admissible term lies in it, which prunes the box); + :func:`generator_heaviest_term_bound` uses it to restrict the heavy term to generator states. """ P = np.asarray(P, dtype=float) P = P / P.sum() @@ -135,26 +145,68 @@ def heaviest_term_bound( mu1 = mu1_vertex(P) t_hi = min(1.0, mu1.value if mu1 is not None else 1.0) t_lo = 1.0 / (2 ** min(na, nb) - 1) - root: Box = (np.zeros(na), np.ones(na), np.zeros(nb), np.ones(nb), t_lo, t_hi) - first = _box_bound(P, root) + + def bound(box: Box) -> tuple[float | None, Box]: + if weight_cap is not None: + cap = weight_cap(box) + if cap is None or cap < box[4]: + return None, box + box = (*box[:5], min(box[5], cap)) + return _box_bound(P, box), box + + first, root = bound((np.zeros(na), np.ones(na), np.zeros(nb), np.ones(nb), t_lo, t_hi)) heap: list[tuple[float, int, Box]] = [] if first is None else [(first, 0, root)] order = itertools.count(1) nodes = 0 goal = np.inf if target is None else target while heap and nodes < max_nodes and (deadline is None or time.monotonic() < deadline): - bound, _, box = heap[0] - if bound >= goal: + lowest, _, box = heap[0] + if lowest >= goal: break heapq.heappop(heap) nodes += 1 qlo, qhi, rlo, rhi, tlo, thi = box widths = np.concatenate([qhi - qlo, rhi - rlo, [thi - tlo]]) - for child in _split(box, int(np.argmax(widths))): - value = _box_bound(P, child) + for split in _split(box, int(np.argmax(widths))): + value, child = bound(split) if value is not None: - heapq.heappush(heap, (max(bound, value), next(order), child)) + heapq.heappush(heap, (max(lowest, value), next(order), child)) value = heap[0][0] if heap else float("inf") return HeaviestTermBoundResult(value, nodes, len(heap), target is not None and value >= target) -__all__ = ["HeaviestTermBoundResult", "capped_coverage_entropy", "heaviest_term_bound"] +def generator_heaviest_term_bound( + problem: _GeneratorProblem, + *, + target: float | None = None, + max_nodes: int = 200_000, + deadline: float | None = None, +) -> HeaviestTermBoundResult: + """Certified lower bound on the generative complexity :math:`C_g` by branching on the heaviest *state*. + + A generator's states give an exact decomposition of the causal-state joint + :math:`p(s^+, s^-) = \\sum_w \\pi_w Q_w R_w^\\top`, whose heaviest term is its heaviest state; that state also + satisfies the generator equations. So :func:`heaviest_term_bound` applies with an extra cap on the heavy + weight: the one-explicit-state generator-cap LP over the box of the state's beliefs + (:func:`~sofic.generators.generator_caps.generator_weight_cap`). Where that LP excludes the + :math:`\\op{G}{S^+ : S^-}`-optimal heavy term the bound exceeds :math:`\\op{G}{S^+ : S^-}`. + + The root's weight floor :math:`1/(2^m - 1)`, :math:`m = \\min(|S^+|, |S^-|)`, is the cardinality bound on + optimal decompositions, and a generator may have more states. It stays valid: a generator whose + heaviest state weighs less has :math:`\\H{\\pi} \\ge \\log_2(2^m - 1) \\ge \\log_2 m \\ge \\min(C_\\mu^+, + C_\\mu^-) \\ge C_g`, so it cannot lower :math:`C_g`. + """ + from sofic.generators.generator_caps import _HeaviestStatesRelaxation + + relaxation = _HeaviestStatesRelaxation(problem, 1) + P = problem.lifted.plus.T @ (problem.p_pair[:, None] * problem.lifted.minus) + return heaviest_term_bound( + P, + target=target, + max_nodes=max_nodes, + deadline=deadline, + weight_cap=lambda box: relaxation.bound((box[:4],)), + ) + + +__all__ = ["HeaviestTermBoundResult", "capped_coverage_entropy", "generator_heaviest_term_bound", "heaviest_term_bound"] diff --git a/sofic/generators/generative_search.py b/sofic/generators/generative_search.py index 4534e84..cc4c2a2 100644 --- a/sofic/generators/generative_search.py +++ b/sofic/generators/generative_search.py @@ -60,7 +60,7 @@ targeted_line_split_seeds, ) from sofic.generators.epsilon_machine import EpsilonMachine -from sofic.generators.exact_ci_bounds import heaviest_term_bound +from sofic.generators.exact_ci_bounds import generator_heaviest_term_bound, heaviest_term_bound from sofic.generators.generator_bounds import ( flow_bound_result, maximal_rectangles, @@ -340,6 +340,7 @@ def _bracket_search( coverage_mu_starts: int = 25, coverage_k_max: int = 8, coverage_certified_only: bool = False, + generator_heaviest_term: bool = False, generator_caps: int = 2, targeted: bool = True, entropy_sweep_starts: int = 80, @@ -623,6 +624,26 @@ def maybe_return_early() -> GenerativeModelResult | None: ) if htb.value > lower + 1e-9: lower, lower_source, lower_certified = htb.value, "heaviest-term bound (certified)", True + best_upper = min((c[0] for c in candidates), default=np.inf) + if ( + generator_heaviest_term + and not wyner + and lower < best_upper - _CERT_TOL + and problem.n_plus + problem.n_minus <= _HEAVIEST_TERM_MAX_DIMS + and not out_of_time("generator heaviest-term bound") + ): + gtb = generator_heaviest_term_bound(problem, target=best_upper - 0.5 * _CERT_TOL, deadline=share(0.3)) + log.append( + { + "stage": "generator heaviest-term bound (certified)", + "value": gtb.value, + "nodes": gtb.nodes, + "reached": gtb.reached, + } + ) + if gtb.value > lower + 1e-9: + lower, lower_source = gtb.value, "generator heaviest-term bound (certified)" + lower_certified = True if bilevel and n_generating == 0 and optima and not out_of_time("bilevel search"): for ent, viol, c in bilevel_search(bidir, optima[:8], deadline=share(0.15))[:4]: log.append({"stage": "bilevel search", "entropy": ent, "kernel violation": viol}) @@ -841,6 +862,7 @@ def generative_model_from_bidirectional( coverage_mu_starts: int = 25, coverage_k_max: int = 8, coverage_certified_only: bool = False, + generator_heaviest_term: bool = False, generator_caps: int = 2, targeted: bool = True, entropy_sweep_starts: int = 80, @@ -870,6 +892,7 @@ def generative_model_from_bidirectional( coverage_mu_starts=coverage_mu_starts, coverage_k_max=coverage_k_max, coverage_certified_only=coverage_certified_only, + generator_heaviest_term=generator_heaviest_term, generator_caps=generator_caps, targeted=targeted, entropy_sweep_starts=entropy_sweep_starts, @@ -965,6 +988,7 @@ def generative_model( coverage_mu_starts: int = 25, coverage_k_max: int = 8, coverage_certified_only: bool = False, + generator_heaviest_term: bool = False, generator_caps: int = 2, targeted: bool = True, entropy_sweep_starts: int = 80, @@ -1013,6 +1037,7 @@ def generative_model( coverage_mu_starts=coverage_mu_starts, coverage_k_max=coverage_k_max, coverage_certified_only=coverage_certified_only, + generator_heaviest_term=generator_heaviest_term, generator_caps=generator_caps, targeted=targeted, entropy_sweep_starts=entropy_sweep_starts, diff --git a/tests/test_exact_ci_bounds.py b/tests/test_exact_ci_bounds.py index 6fbc22c..67237d7 100644 --- a/tests/test_exact_ci_bounds.py +++ b/tests/test_exact_ci_bounds.py @@ -61,3 +61,53 @@ def test_survey_case_closes_against_its_exact_common_information() -> None: result = heaviest_term_bound(P, target=best - 1e-2, max_nodes=100_000) assert result.reached assert best - 1e-2 <= result.value <= best + 1e-9 + + +def test_weight_cap_that_admits_no_term_prunes_every_box() -> None: + D = np.array([[0.3, 0.2], [0.2, 0.3]]) + result = heaviest_term_bound(D, max_nodes=10, weight_cap=lambda box: None) + assert result.value == np.inf + assert result.open_boxes == 0 + + +def _survey_case_a(): + """Survey process ``((0,1),(0,2),(0,1))`` (v11 seed 54), where C_g >= 0.928 > G = 0.92433.""" + from sofic.examples.epsilon_machines import from_symbol_matrices + from sofic.generators.bidirectional_epsilon_machine import BidirectionalEpsilonMachine + from sofic.generators.coupled_generators import _GeneratorProblem + from sofic.generators.epsilon_machine import EpsilonMachine + + a, b, c = 0.06761675137492751, 0.44247275706537287, 0.45789456421998476 + mats = { + 0: np.array([[a, 0, 0], [b, 0, 0], [c, 0, 0]]), + 1: np.array([[0, 1 - a, 0], [0, 0, 1 - b], [0, 1 - c, 0]]), + } + machine = from_symbol_matrices((0, 1, 2), (0, 1), mats) + forward = machine if isinstance(machine, EpsilonMachine) else EpsilonMachine.from_hmm(machine) + return _GeneratorProblem(BidirectionalEpsilonMachine.from_forward(forward)) + + +def test_generator_equations_exclude_the_g_optimal_heavy_term() -> None: + from sofic.generators.generator_caps import _HeaviestStatesRelaxation + + problem = _survey_case_a() + # heaviest term of the G-optimal decomposition (entropy 0.92433) + t = 0.6943850024064206 + q = np.array([0.06761675137400137, 0.5912539107207103, 0.3411293379052883]) + r = np.array([0.18172612379896527, 0.46738361190932987, 0.3508902642917049]) + relaxation = _HeaviestStatesRelaxation(problem, 1) + wide = relaxation.bound( + ((np.clip(q - 1e-2, 0, 1), np.clip(q + 1e-2, 0, 1), np.clip(r - 1e-2, 0, 1), np.clip(r + 1e-2, 0, 1)),) + ) + tight = relaxation.bound( + ((np.clip(q - 1e-4, 0, 1), np.clip(q + 1e-4, 0, 1), np.clip(r - 1e-4, 0, 1), np.clip(r + 1e-4, 0, 1)),) + ) + assert wide is not None and wide > t + assert tight is not None and tight < t - 1e-3 + + +def test_generator_heaviest_term_bound_is_below_a_known_generator() -> None: + from sofic.generators.exact_ci_bounds import generator_heaviest_term_bound + + result = generator_heaviest_term_bound(_survey_case_a(), max_nodes=300) + assert 0.85 < result.value <= 0.9697134 + 1e-9 # entropy of a word-verified generator