diff --git a/CHANGELOG.md b/CHANGELOG.md index 24889a0..c08d72b 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/docs/generators/bidirectional_epsilon_machine.rst b/docs/generators/bidirectional_epsilon_machine.rst index f87a872..23f3b14 100644 --- a/docs/generators/bidirectional_epsilon_machine.rst +++ b/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 diff --git a/docs/generators/generative_models.rst b/docs/generators/generative_models.rst index 5959c27..59a01cf 100644 --- a/docs/generators/generative_models.rst +++ b/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 ------ diff --git a/sofic/generators/bidirectional_epsilon_machine.py b/sofic/generators/bidirectional_epsilon_machine.py index d265fa0..c126b47 100644 --- a/sofic/generators/bidirectional_epsilon_machine.py +++ b/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. diff --git a/sofic/generators/lifted_common_information.py b/sofic/generators/lifted_common_information.py index e31b130..7d9ab19 100644 --- a/sofic/generators/lifted_common_information.py +++ b/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 diff --git a/sofic/generators/minimal_generative_model.py b/sofic/generators/minimal_generative_model.py index 272830a..3358bf7 100644 --- a/sofic/generators/minimal_generative_model.py +++ b/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", diff --git a/tests/test_lifted_common_information.py b/tests/test_lifted_common_information.py index 54b4d3e..f902b2a 100644 --- a/tests/test_lifted_common_information.py +++ b/tests/test_lifted_common_information.py @@ -107,6 +107,32 @@ def test_rll_two_infinity_constraints_bind() -> None: assert _reproduces(model, bidir) +def test_generative_gap_rll_two_infinity() -> None: + bidir = _rll_two_infinity() + exact = bidir.generative_gap(rng=np.random.default_rng(0)) + assert exact.measure == "exact" + assert exact.lower == pytest.approx(1.4056, abs=2e-3) + assert exact.upper == pytest.approx(1.5, abs=1e-4) + assert exact.gap == pytest.approx(exact.upper - exact.lower) + assert exact.gap > 0.09 + assert exact.upper == pytest.approx(exact.model.generative_complexity(), abs=1e-6) + wyner = bidir.generative_gap(measure="wyner", rng=np.random.default_rng(0)) + assert wyner.lower == pytest.approx(1.129, abs=5e-3) + assert wyner.gap > 0.1 + + +def test_generative_gap_vanishes_for_markov_chains() -> None: + bidir = golden_mean_bidirectional(0.5) + for measure in ("exact", "wyner"): + result = bidir.generative_gap(measure=measure, rng=np.random.default_rng(0)) + assert 0.0 <= result.gap < 1e-3 + + +def test_generative_gap_rejects_unknown_measure() -> None: + with pytest.raises(ValueError, match="measure"): + golden_mean_bidirectional(0.5).generative_gap(measure="functional") + + def test_functional_model_always_generates() -> None: """dit's functional partition fails the lifting conditions here; the lifted partition search does not.""" bidir = _rll_two_infinity()