Heaviest-term lower bound on C_g over generator states

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
 

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

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,40 +126,87 @@ 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.
 
     Branches best first on the widest coordinate of the heaviest term's box and stops once every open box
     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()
     na, nb = P.shape
     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"]

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,

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