Lower bound on lifted generators

CHANGELOG.md

@@ -78,6 +78,11 @@ brute-force reference or a closed-form value and has a regression test.
   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`).
+- `lifted_heaviest_term_bound` bounds the state entropy of every *lifted* generator from below. A lifted
+  generator's heaviest state has channel `t q_a r_b / p(a, b)` at both ends of each edge; coarse-graining
+  the other end to "is that state" turns each lifting condition into a condition on `(q, r)` alone, checked
+  on boxes by interval arithmetic inside `heaviest_term_bound`. On `((0,1),(0,2),(0,2))` it shows every
+  lifted generator has entropy >= 1.4815 while `C_g = G = 1.4225`.
 
 ### Generator caps
 

docs/generators/generative_models.rst

@@ -316,6 +316,19 @@ weight is also capped by the one-explicit-state generator-cap LP over the box of
 :math:`10^{-4}`. The pipeline therefore runs it only on request (``generator_heaviest_term=True``).
 This combination is not from the literature.
 
+Lifted generators get their own bound.
+:func:`~sofic.generators.exact_ci_bounds.lifted_heaviest_term_bound` uses that a lifted generator draws
+consecutive states from one channel, so its heaviest state :math:`e` has channel
+:math:`c_e(a, b) = t\, q_a r_b / p(a, b)` at both ends of every edge. Coarse-graining the other end to the
+event that it is :math:`e` turns each lifting condition into a condition on :math:`(q, r)` alone: a
+polynomial in the beliefs must be the same for every supported forward (or reverse) state. Its
+coefficients are nonnegative, so it is checked on boxes by interval arithmetic and prunes the
+heaviest-term branching. A coupled generator need not pass the test. On dataset process
+``((0,1),(0,2),(0,2))`` (a :math:`3 \times 3` joint) every lifted generator has entropy at least
+:math:`1.4815`, with a lifted generator at :math:`1.4817`, while :math:`C_g = \op{G}{S^+ : S^-} = 1.4225`:
+the minimal generator is not lifted even though an exact-common-information optimum generates. This
+bound 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
@@ -373,6 +386,7 @@ Coupled generators and bounds
 .. 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.lifted_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

@@ -23,7 +23,7 @@
 import time
 from collections.abc import Callable
 from dataclasses import dataclass
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Any
 
 import numpy as np
 
@@ -209,4 +209,95 @@ def generator_heaviest_term_bound(
     )
 
 
-__all__ = ["HeaviestTermBoundResult", "capped_coverage_entropy", "generator_heaviest_term_bound", "heaviest_term_bound"]
+class _LiftedTermTest:
+    """Necessary conditions for a decomposition term to be a state of a *lifted* generator (see below)."""
+
+    def __init__(self, lifted: Any) -> None:
+        self.pa = lifted.plus.argmax(axis=1)
+        self.pb = lifted.minus.argmax(axis=1)
+        self.p = lifted.p
+        self.edges = lifted.edges
+        self.na, self.nb = lifted.plus.shape[1], lifted.minus.shape[1]
+
+    def _forward(self, q: np.ndarray, r: np.ndarray) -> np.ndarray:
+        u = q[self.pa] * r[self.pb] / self.p
+        inner = np.einsum("pxq,q->px", self.edges, u)
+        out = np.zeros((self.na, inner.shape[1]))
+        np.add.at(out, self.pa, (r[self.pb] / self.p)[:, None] * inner)
+        return out
+
+    def _reverse(self, q: np.ndarray, r: np.ndarray) -> np.ndarray:
+        u = q[self.pa] * r[self.pb] / self.p
+        inner = np.einsum("pxq,p->qx", self.edges, u)
+        out = np.zeros((self.nb, inner.shape[1]))
+        np.add.at(out, self.pb, (q[self.pa] / self.p)[:, None] * inner)
+        return out
+
+    def violation(self, q: np.ndarray, r: np.ndarray, *, support_tol: float = 1e-12) -> float:
+        """Largest spread of ``L_B`` over ``supp q`` or of ``L_A`` over ``supp r`` (0 for a lifted state)."""
+        worst = 0.0
+        for values, support in ((self._forward(q, r), q > support_tol), (self._reverse(q, r), r > support_tol)):
+            rows = values[support]
+            if len(rows) > 1:
+                worst = max(worst, float((rows.max(axis=0) - rows.min(axis=0)).max()))
+        return worst
+
+    def excludes(self, box: Box) -> bool:
+        """Whether no term with beliefs in ``box`` passes: two surely supported states have disjoint intervals."""
+        qlo, qhi, rlo, rhi = box[:4]
+        for f, support in ((self._forward, qlo > 0.0), (self._reverse, rlo > 0.0)):
+            lo, hi = f(qlo, rlo)[support], f(qhi, rhi)[support]
+            if len(lo) > 1 and (lo.max(axis=0) > hi.min(axis=0) + 1e-12).any():
+                return True
+        return False
+
+
+def lifted_heaviest_term_bound(
+    problem: _GeneratorProblem,
+    *,
+    target: float | None = None,
+    max_nodes: int = 200_000,
+    deadline: float | None = None,
+) -> HeaviestTermBoundResult:
+    """Certified lower bound on the state entropy of every *lifted* generator, by branching on its heaviest state.
+
+    A lifted generator draws :math:`W_0` and :math:`W_1` from one channel :math:`c(w \\mid s^+, s^-)`, independently
+    given the bidirectional edge (see :mod:`sofic.generators.lifted_common_information`). Its states form an exact
+    decomposition :math:`p = \\sum_w t_w q_w r_w^\\top`, so state :math:`e` has channel
+    :math:`c_e(a, b) = t\\, q_a r_b / p(a, b)`. Coarse-graining :math:`W_1` (or :math:`W_0`) to the event that it
+    equals :math:`e` cannot raise the lifting conditions' conditional mutual informations, so with :math:`E` the
+    joint edge law and :math:`u(a, b) = q_a r_b / p(a, b)`, for every symbol :math:`x`
+
+    .. math::
+
+       L_B(a, x) = \\sum_b \\frac{r_b}{p(a, b)} \\sum_{a', b'} E[(a, b), x, (a', b')]\\, u(a', b')
+       \\quad\\text{is constant over } a \\in \\operatorname{supp} q,
+
+       L_A(b', x) = \\sum_{a'} \\frac{q_{a'}}{p(a', b')} \\sum_{a, b} E[(a, b), x, (a', b')]\\, u(a, b)
+       \\quad\\text{is constant over } b' \\in \\operatorname{supp} r
+
+    (the factor :math:`t` cancels). A general generator need not satisfy them: it may couple consecutive states.
+    The coefficients are nonnegative, so a box's lower and upper corners bound :math:`L_A, L_B`; a box is pruned
+    when two surely supported states get disjoint intervals, inside :func:`heaviest_term_bound`. The root's weight
+    floor stays valid because the forward ε-machine is a lifted generator of entropy
+    :math:`C_\\mu^+ \\le \\log_2 |S^+|`, so a lifted generator whose heaviest state weighs less than
+    :math:`1/(2^m - 1)` cannot be minimal.
+    """
+    test = _LiftedTermTest(problem.lifted)
+    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: None if test.excludes(box) else box[5],
+    )
+
+
+__all__ = [
+    "HeaviestTermBoundResult",
+    "capped_coverage_entropy",
+    "generator_heaviest_term_bound",
+    "heaviest_term_bound",
+    "lifted_heaviest_term_bound",
+]

tests/test_exact_ci_bounds.py

@@ -111,3 +111,42 @@ def test_generator_heaviest_term_bound_is_below_a_known_generator() -> None:
 
     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
+
+
+def _dataset_case_t055():
+    """Dataset process ``n3-t055-d2`` (code ((0,1),(0,2),(0,2))): C_g = G = 1.42254, minimal lifted model 1.4815+."""
+    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.29501936538151036, 0.7629418805405772, 0.1416369489322688
+    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, 0, 1 - c]]),
+    }
+    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_epsilon_machine_states_pass_the_lifted_term_test() -> None:
+    from sofic.generators.exact_ci_bounds import _LiftedTermTest
+
+    problem = _dataset_case_t055()
+    test = _LiftedTermTest(problem.lifted)
+    _, Q, R, _ = problem.epsilon_seed()
+    assert all(test.violation(q, r) < 1e-12 for q, r in zip(Q, R, strict=True))
+
+
+def test_lifted_bound_separates_lifted_generators_from_g() -> None:
+    from sofic.generators.exact_ci_bounds import lifted_heaviest_term_bound
+
+    problem = _dataset_case_t055()
+    P = problem.lifted.plus.T @ (problem.p_pair[:, None] * problem.lifted.minus)
+    lifted = lifted_heaviest_term_bound(problem, target=1.43, max_nodes=50_000)
+    static = heaviest_term_bound(P, target=1.43, max_nodes=50_000)
+    assert lifted.reached  # every lifted generator has entropy >= 1.43
+    assert static.value <= 1.42255 + 1e-6  # a coupled generator of entropy 1.4225448 exists, so G <= 1.42255
+    c_mu = float(-(problem.p_plus * np.log2(problem.p_plus)).sum())
+    assert lifted.value <= c_mu  # the epsilon-machine is a lifted generator