generative_gap, and describe the lifted search as restricted to lifted generators

CHANGELOG.md

@@ -9,15 +9,19 @@ brute-force reference or a closed-form value and has a regression test.
 
 - `minimal_generative_model`, `wyner_generative_model`,
   `functional_generative_model` and `generative_complexity()` optimize only
-  over generator states: besides `I[S+ : S- | G] = 0` they impose
+  over lifted generator states: besides `I[S+ : S- | G] = 0` they impose
   `I[G_0 X_0 : S-_1 | G_1] = 0` and `I[S+_0 : X_0 G_1 | G_0] = 0`, the
-  conditions under which the HMM built from the channel generates the process.
+  conditions under which the HMM built from the channel (applied at both ends
+  of every bidirectional edge) generates the process. Lifted generators are a
+  proper subclass, so the exact and Wyner values are upper bounds over all
+  generators; the functional search is exact over deterministic ones.
   The reported value is the returned model's own (no functional fallback), the
   exact value upper-bounds the generative complexity, and models carry
   `markov_residual`. The functional model, previously returned even when it
   did not reproduce the process, now always does. The unconstrained values are
   `BidirectionalEpsilonMachine.wyner_common_information()` and
-  `exact_common_information()`.
+  `exact_common_information()`, and `generative_gap()` reports how much the
+  lifted-generator requirement raises each (`G[S+ : S-] <= C_g <= H[G]`).
 - `correction="bonferroni"` is the default for every CSSR learner (process,
   subtree, transducer, stack), and process CSSR counts its length-`(L+1)`
   resolution tests; pass `correction=None` for the old behavior.

docs/generators/bidirectional_epsilon_machine.rst

@@ -137,4 +137,4 @@ API
 ===
 
 .. autoclass:: BidirectionalEpsilonMachine
-   :members: from_forward, from_pair, joint_distribution, step_distribution, forward_epsilon_machine, reverse_epsilon_machine, entropy_rate, statistical_complexity, excess_entropy, bidirectional_crypticity, minimal_generative_model, wyner_generative_model, functional_generative_model, gacs_korner_generative_model, generative_complexity, wyner_common_information, exact_common_information, predicted_information, bound_information, ephemeral_information, structural_ephemeral_information, parallel_edge_information, bound_structural_information, bound_parallel_edge_information, information_anatomy, caekl_causal_information
+   :members: from_forward, from_pair, joint_distribution, step_distribution, forward_epsilon_machine, reverse_epsilon_machine, entropy_rate, statistical_complexity, excess_entropy, bidirectional_crypticity, minimal_generative_model, wyner_generative_model, functional_generative_model, gacs_korner_generative_model, generative_complexity, wyner_common_information, exact_common_information, generative_gap, predicted_information, bound_information, ephemeral_information, structural_ephemeral_information, parallel_edge_information, bound_structural_information, bound_parallel_edge_information, information_anatomy, caekl_causal_information

docs/generators/generative_models.rst

@@ -12,7 +12,8 @@ these minimal generators from a
 :class:`~sofic.generators.bidirectional_epsilon_machine.BidirectionalEpsilonMachine`,
 each realizing a different **common information** between the forward causal
 state :math:`S^+` and the reverse causal state :math:`S^-`, restricted to
-auxiliaries that are the states of a generator (see `Lifted generators`_).
+auxiliaries that are the states of a *lifted* generator (see
+`Lifted generators`_).
 
 The generative complexity :math:`C_g = \H{G}` of a presentation with generative
 states :math:`G` obeys :math:`C_g \le C_\mu`. The three generator
@@ -27,12 +28,12 @@ where the Wyner value is :math:`\op{I}{(S^+, S^-) : G}` and the others are
 :math:`\H{G}` (the first two up to optimizer accuracy).
 
 * :func:`minimal_generative_model` — the **exact common information**
-  :cite:`Kumar2014`: the minimum-entropy generator state :math:`G` rendering
-  :math:`S^+` and :math:`S^-` conditionally independent, found by a
+  :cite:`Kumar2014`: the minimum-entropy lifted generator state :math:`G`
+  rendering :math:`S^+` and :math:`S^-` conditionally independent, found by a
   stochastic optimizer. Its entropy bounds :math:`C_g` from above.
 * :func:`wyner_generative_model` — the **Wyner common information**
   :cite:`Wyner1975`, minimizing :math:`\op{I}{(S^+, S^-) : G}` over the same
-  generator states.
+  lifted generator states.
 * :func:`functional_generative_model` — the **functional common information**:
   the smallest *deterministic* function of :math:`(S^+, S^-)` that separates
   the two causal states and is a generator state, found by a partition search.
@@ -91,13 +92,27 @@ returns, bounds the generative complexity of Löhr & Ay :cite:`LohrAy2009` from
 above. ``markov_residual`` reports the largest constraint information (bits) at
 the returned channel.
 
+Lifted generators are a proper subclass of generators: the lift draws
+:math:`G_0` and :math:`G_1` independently given the bidirectional edge, while a
+general generator may couple them. For
+:math:`Z_t = Z_{t-2} \oplus \mathrm{Bernoulli}(0.1)` (two interleaved
+first-order chains), the product of the two chains' minimal generators has
+state :math:`G_t = (U_t, U_{t+1})`, which shares :math:`U_{t+1}` with
+:math:`G_{t+1}`. It is a generator state of entropy 1.982, essentially
+:math:`\op{G}{S^+ : S^-}`, but violates both lifting conditions (by 0.019 bits),
+so the lifted search cannot reach it and returns a larger value. The exact and
+Wyner values here are therefore upper bounds over lifted generators, not
+minima over all generators. For a deterministic channel the coupling is
+forced, so the functional search is exact over all deterministic generator
+states.
+
 For first-order Markov processes the constrained and unconstrained values
 coincide. They differ already for the second-order Markov process in which
 every 1 is followed by at least two 0s: its unconstrained exact (and functional)
 optimum, of entropy 1.406, is a deterministic variable whose value
 :math:`(C, R_{01})` ("the past ends in 10 and the future starts with 01") can
 only be preceded by a value that does not know the past ended in 10, so no HMM
-has it as its state; the constrained value is 1.5. The stochastic optimizer is
+has it as its state; the lifted search finds 1.5. The stochastic optimizer is
 local and seeded (``starts``, ``rng``); on harder processes more starts can
 find lower values.
 
@@ -109,6 +124,23 @@ common informations of the past and future — and are available as measures:
    bidir.wyner_common_information()   # C[S+ : S-] = C[past : future]
    bidir.exact_common_information()   # G[S+ : S-] = G[past : future]
 
+Because every generator state separates the causal states, they bracket the
+generative complexity: :math:`\op{G}{S^+ : S^-} \le C_g \le \H{G}` for the
+state :math:`G` of :func:`minimal_generative_model`. :func:`generative_gap`
+reports both ends and their difference, an upper bound on
+:math:`C_g - \op{G}{S^+ : S^-}`, and the analogous Wyner comparison:
+
+.. code-block:: python
+
+   gap = bidir.generative_gap()            # measure="exact"
+   gap.lower, gap.upper, gap.gap           # 1.406, 1.5, 0.094 for the process above
+   bidir.generative_gap(measure="wyner")   # 1.129, 1.311, 0.182
+
+Both ends are optimizer estimates, and the upper end is over lifted generators
+only, so ``gap`` can be positive even when some non-lifted generator attains
+:math:`\op{G}{S^+ : S^-}` (the interleaved process above). The gap vanishes for
+first-order Markov processes.
+
 API
 ===
 
@@ -125,6 +157,9 @@ Measures
 
 .. autofunction:: wyner_common_information
 .. autofunction:: exact_common_information
+.. autofunction:: generative_gap
+.. autoclass:: GenerativeGap
+   :members: gap
 
 Models
 ------

sofic/generators/bidirectional_epsilon_machine.py

@@ -620,6 +620,12 @@ def exact_common_information(self, **kwargs: Any) -> float:
 
         return exact_common_information(self, **kwargs)
 
+    def generative_gap(self, measure: str = "exact", **kwargs: Any) -> Any:
+        """Exact (or Wyner) common information vs. its value over generator states; see ``generative_gap``."""
+        from sofic.generators.minimal_generative_model import generative_gap
+
+        return generative_gap(self, measure=measure, **kwargs)
+
     def information_anatomy(self) -> dict[str, Any]:
         """Return ρ_μ, b_μ, r_μ, h_μ, E, χ, χ±, and the structural/gauge refinement.
 

sofic/generators/lifted_common_information.py

@@ -40,14 +40,26 @@
 Every feasible point is a generator state, so the lifted exact value bounds
 the generative complexity :math:`C_g` of Löhr & Ay :cite:`LohrAy2009` from
 above. Both lifted values bound the corresponding unconstrained common
-informations from above; they coincide for first-order Markov processes. They
-can differ because only generators whose kernel is the lifted one are
-searched, and because the optimizer is local.
+informations from above; they coincide for first-order Markov processes.
+
+Lifted generators are a proper subclass of generators. The lift draws
+:math:`G_0` and :math:`G_1` independently given the bidirectional edge; a
+general generator may couple them. For
+:math:`Z_t = Z_{t-2} \oplus \mathrm{Bernoulli}(0.1)`, the interleaved product of
+two first-order generators has state :math:`G_t = (U_t, U_{t+1})`, sharing
+:math:`U_{t+1}` with :math:`G_{t+1}`. It is a generator state of entropy 1.982,
+essentially :math:`\op{G}{S^+ : S^-}`, yet violates both lifting conditions by
+0.019 bits, so the lifted search cannot reach it. For a deterministic channel
+the coupling is forced, so lifting is no restriction: the functional search
+optimizes over all deterministic generator states.
+
+The lifted values can also exceed the optimum over lifted generators because the
+optimizer is local.
 """
 
 from __future__ import annotations
 
-from collections.abc import Hashable
+from collections.abc import Hashable, Sequence
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -77,6 +89,9 @@
 _SUPPORT_CUTOFF = 1e-8
 _REFINE_TOL = 1e-12
 _EXACT_ATOL = 1e-10
+_PERTURB_SCALES = (0.05, 0.2)
+_PERTURB_COPIES = 2
+_MIX_WEIGHTS = (0.5, 0.8)
 _EXHAUSTIVE_PAIRS = 9
 
 
@@ -245,7 +260,11 @@ def _optimize(
     tol: float,
     rng: np.random.Generator,
     extra: list[np.ndarray],
+    more_starts: Sequence[np.ndarray] = (),
+    generator_seeds: bool = True,
 ):
+    """Continuation from the generator seeds (if ``generator_seeds``), the ``extra`` channels, ``more_starts``
+    and ``starts`` random channels; the ``extra`` channels are also kept and refined as they are."""
     n = len(problem.pairs)
     shape = (n, k)
     row_sums = {
@@ -283,9 +302,9 @@ def solve(c0: np.ndarray, closure_tol: float, screen_tol: float, support: np.nda
     # Continuation to a coarse tolerance: the plain common-information problem first, then ever tighter closure.
     coarse = max(tol, _COARSE_TOL)
     schedule = [1e3, *np.geomspace(1e-1, coarse, max(2, int(np.ceil(np.log10(1e-1 / coarse) / 2)) + 1))]
-    fixed = problem.seeds(k) + extra
+    fixed = (problem.seeds(k) if generator_seeds else []) + extra
     candidates = [(problem.objective(c0, kind)[0], problem.markov_residual(c0), c0) for c0 in fixed]
-    for c0 in fixed + [rng.dirichlet(np.ones(k), size=n) for _ in range(starts)]:
+    for c0 in fixed + list(more_starts) + [rng.dirichlet(np.ones(k), size=n) for _ in range(starts)]:
         c = c0
         for closure_tol in schedule:
             c = solve(c, closure_tol, coarse)
@@ -366,11 +385,28 @@ def lifted_model(
                 seed /= seed.sum(axis=1, keepdims=True)
             k_seed = max(seed.shape[1], k if bound is not None else 0)
             seed = np.pad(seed, ((0, 0), (0, k_seed - seed.shape[1])))
+            perturbed = [
+                (1 - eps) * seed + eps * rng.dirichlet(np.ones(k_seed), size=len(problem.pairs))
+                for eps in _PERTURB_SCALES
+                for _ in range(_PERTURB_COPIES)
+            ]
             seed_refined, seed_candidates = _optimize(
-                problem, objective, k_seed, 0, maxiter, markov_atol / 2, rng, [seed]
+                problem, objective, k_seed, 0, maxiter, markov_atol / 2, rng, [seed], perturbed, generator_seeds=False
             )
             refined += seed_refined
             candidates += seed_candidates
+            # Time-sharing mixtures: an independent coin picks the unconstrained channel or an eps-machine
+            # channel; disjoint state labels keep the mixture separating.
+            for machine in (problem.plus, problem.minus, np.eye(len(problem.pairs))):
+                width = k_seed + machine.shape[1]
+                if bound is not None and width > k:
+                    continue
+                mixtures = [np.hstack([lam * seed, (1 - lam) * machine]) for lam in _MIX_WEIGHTS]
+                mix_refined, mix_candidates = _optimize(
+                    problem, objective, width, 0, maxiter, markov_atol / 2, rng, [], mixtures, generator_seeds=False
+                )
+                refined += mix_refined
+                candidates += mix_candidates
         # Prefer refined channels: looser candidates can undercut the true optimum by about sqrt(residual).
         feasible = [cand for cand in refined if cand[1] <= markov_atol] or [
             cand for cand in candidates if cand[1] <= markov_atol

sofic/generators/minimal_generative_model.py

@@ -4,6 +4,7 @@
 
 from collections import defaultdict
 from collections.abc import Hashable
+from dataclasses import dataclass
 from typing import TYPE_CHECKING, Any, cast
 
 import numpy as np
@@ -163,12 +164,16 @@ def minimal_generative_model(
     is the model's own state entropy, an upper bound on the generative complexity
     :cite:`LohrAy2009`; the unconstrained value is
     :meth:`~sofic.generators.bidirectional_epsilon_machine.BidirectionalEpsilonMachine.exact_common_information`.
+    Lifted generators draw consecutive states independently given the
+    bidirectional edge; a generator that couples them can have smaller state
+    entropy, so this is not necessarily the minimum over all generators.
 
     The starts are ``S+``, ``S-``, ``(S+, S-)``, the unconstrained optimum (found
     with dit; ``niter``, ``maxiter``, ``polish`` and ``backend`` as in
-    :class:`dit.multivariate.common_informations.ExactCommonInformation`), and
-    ``starts`` random channels over ``bound`` states (default: the number of
-    joint states). Constraint informations up to ``markov_atol`` bits are
+    :class:`dit.multivariate.common_informations.ExactCommonInformation`),
+    perturbations of it (blended with random channels), time-sharing mixtures of
+    it with each of the three ε-machine channels, and ``starts`` random channels
+    over ``bound`` states (default: the number of joint states). Constraint informations up to ``markov_atol`` bits are
     accepted and the best candidates refined toward exact closure; a residual
     ``r`` allows the value and word probabilities to be off by roughly
     ``sqrt(r)``. The model is checked against the source up to
@@ -214,8 +219,9 @@ def wyner_generative_model(
     The auxiliary state ``G`` minimizes ``I[(S+, S-) : G]`` :cite:`Wyner1975`
     subject to rendering ``S+`` and ``S-`` conditionally independent and to the
     two lifting conditions of :func:`minimal_generative_model`, whose parameters
-    it shares. The returned model's state entropy ``H[G]`` is not generally equal
-    to that optimized mutual information. The unconstrained value is
+    it shares, and like it searches lifted generators only. The returned model's
+    state entropy ``H[G]`` is not generally equal to that optimized mutual
+    information. The unconstrained value is
     :meth:`~sofic.generators.bidirectional_epsilon_machine.BidirectionalEpsilonMachine.wyner_common_information`.
     """
     from sofic.generators.lifted_common_information import lifted_model
@@ -335,6 +341,95 @@ def _unconstrained_common_information(
     return value
 
 
+@dataclass(frozen=True)
+class GenerativeGap:
+    """Gap between a common information of ``(S+, S-)`` and its value over lifted generator states.
+
+    Attributes
+    ----------
+    measure
+        ``"exact"`` or ``"wyner"``.
+    lower
+        The unconstrained common information: ``G[S+ : S-]`` or ``C[S+ : S-]``.
+    upper
+        The same quantity over lifted generator states: ``H[G]`` of
+        :func:`minimal_generative_model`, or ``I[(S+, S-) : G]`` of
+        :func:`wyner_generative_model`.
+    model
+        The generator achieving ``upper``.
+    """
+
+    measure: str
+    lower: float
+    upper: float
+    model: MinimalGenerativeModel | WynerGenerativeModel
+
+    @property
+    def gap(self) -> float:
+        """``upper - lower``, in bits."""
+        return self.upper - self.lower
+
+
+def generative_gap(
+    bidir: BidirectionalEpsilonMachine,
+    *,
+    measure: str = "exact",
+    bound: int | None = None,
+    niter: int | None = 10,
+    maxiter: int = 1000,
+    polish: float | bool = 1e-8,
+    backend: str = "numpy",
+    starts: int = 8,
+    markov_atol: float = 1e-8,
+    cutoff: float = 1e-10,
+    reproduction_atol: float = 1e-3,
+    rng: np.random.Generator | None = None,
+) -> GenerativeGap:
+    """How much requiring a lifted generator raises the exact (or Wyner) common information.
+
+    Every generator state renders ``S+`` and ``S-`` conditionally independent, so
+    for ``measure="exact"``
+
+    .. math:: \\op{G}{S^+ : S^-} \\le C_g \\le \\H{G},
+
+    with ``G`` the state of :func:`minimal_generative_model` and ``C_g`` the
+    generative complexity :cite:`LohrAy2009`: ``gap`` bounds ``C_g - G[S+ : S-]``
+    from above. Because only lifted generators are searched, ``gap`` can be
+    positive even when a non-lifted generator attains the unconstrained value (see
+    :mod:`sofic.generators.lifted_common_information`). ``measure="wyner"``
+    compares ``C[S+ : S-]`` with :func:`wyner_generative_model`. Both sides come
+    from local optimizers, so the gap is an estimate; it vanishes for first-order
+    Markov processes. Since both sides upper-bound the true unconstrained value,
+    ``lower`` is the smaller of the two. Parameters are those of
+    :func:`minimal_generative_model`.
+    """
+    if measure not in ("exact", "wyner"):
+        raise ValueError(f"measure must be 'exact' or 'wyner', not {measure!r}")
+    rng = rng or np.random.default_rng()
+    optimizer_kwargs = {
+        "bound": bound,
+        "niter": niter,
+        "maxiter": maxiter,
+        "polish": polish,
+        "backend": backend,
+        "cutoff": cutoff,
+        "rng": rng,
+    }
+    if measure == "exact":
+        lower = exact_common_information(bidir, **optimizer_kwargs)
+        model = minimal_generative_model(
+            bidir, starts=starts, markov_atol=markov_atol, reproduction_atol=reproduction_atol, **optimizer_kwargs
+        )
+        upper = model.exact_common_information
+    else:
+        lower = wyner_common_information(bidir, **optimizer_kwargs)
+        model = wyner_generative_model(
+            bidir, starts=starts, markov_atol=markov_atol, reproduction_atol=reproduction_atol, **optimizer_kwargs
+        )
+        upper = model.wyner_common_information
+    return GenerativeGap(measure=measure, lower=min(lower, upper), upper=upper, model=model)
+
+
 def functional_generative_model(
     bidir: BidirectionalEpsilonMachine,
     *,
@@ -886,11 +981,13 @@ def _model_from_channel(
 __all__ = [
     "FunctionalGenerativeModel",
     "GacsKornerGenerativeModel",
+    "GenerativeGap",
     "MinimalGenerativeModel",
     "WynerGenerativeModel",
     "exact_common_information",
     "functional_generative_model",
     "gacs_korner_generative_model",
+    "generative_gap",
     "minimal_generative_model",
     "wyner_common_information",
     "wyner_generative_model",