diff --git a/CHANGELOG.md b/CHANGELOG.md index 19e0ce9..e42a0df 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/generators/generative_models.rst b/docs/generators/generative_models.rst index 1f84a43..7ee60b3 100644 --- a/docs/generators/generative_models.rst +++ b/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 diff --git a/sofic/generators/exact_ci_bounds.py b/sofic/generators/exact_ci_bounds.py index a4844b2..7af4fbb 100644 --- a/sofic/generators/exact_ci_bounds.py +++ b/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", +] diff --git a/tests/test_exact_ci_bounds.py b/tests/test_exact_ci_bounds.py index 67237d7..464de94 100644 --- a/tests/test_exact_ci_bounds.py +++ b/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