From 2014c5c5c6223a6819a681d015c34d090945cce9 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Sat, 14 Nov 2026 17:28:21 -0700 Subject: [PATCH] =?UTF-8?q?=CE=B5-machines=20of=20asymptotically=20synchro?= =?UTF-8?q?nizing=20processes:=20seed=20the=20mixed=20states=20from=20a=20?= =?UTF-8?q?limit=20belief?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Mixed-state presentations of processes whose hidden state no finite word reveals either closed on hundreds of near-boundary beliefs (rounding reflects the posterior's random walk) or exceeded max_states. Seed from a limit belief instead and accept the result only if it generates the source process exactly. Co-authored-by: Cursor --- CHANGELOG.md | 13 ++++ sofic/generators/epsilon_construction.py | 91 ++++++++++++++++++++++-- tests/test_asymptotic_synchronization.py | 34 +++++++++ 3 files changed, 134 insertions(+), 4 deletions(-) create mode 100644 tests/test_asymptotic_synchronization.py diff --git a/CHANGELOG.md b/CHANGELOG.md index ce47d10..c81f776 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,19 @@ A second correctness review: every fix below was reproduced against a brute-force reference or a closed-form value and has a regression test. +### ε-machine construction + +- Processes whose hidden state no finite word reveals (e.g. a period-2 process whose phase stays + uncertain) synchronize only asymptotically. Their mixed-state presentation either closed on hundreds of + beliefs within round-off of the simplex boundary (rounding acts as a reflecting barrier for the + posterior's random walk) or exceeded `max_states`, so `EpsilonMachine.from_hmm` returned a huge + machine or raised `MixedStateExplosionError` although the ε-machine is finite. It now seeds the + presentation from a limit belief (the heaviest near-boundary belief, or the belief after filtering a + long sample) with sub-`1e-8` entries set to zero, and accepts the result only if it generates the + source process exactly (difference-automaton test). On the 138 minimal binary 3-state topologies, + 9 reverse-machine "explosions" are now 3-state reverse ε-machines, and finite/infinite now agrees + with `reverse_epsilon_machine_is_finite` on every topology. + ### Generator verification - The one-step generator equations (edge law, forward stationarity, summed future-side consistency) are diff --git a/sofic/generators/epsilon_construction.py b/sofic/generators/epsilon_construction.py index 5cdc27f..e42cad4 100644 --- a/sofic/generators/epsilon_construction.py +++ b/sofic/generators/epsilon_construction.py @@ -14,11 +14,11 @@ import numpy as np -from sofic.exceptions import StochasticValidationError +from sofic.exceptions import MixedStateExplosionError, StochasticValidationError from sofic.generators.base import HiddenMarkovModel from sofic.generators.epsilon_machine import EpsilonMachine from sofic.generators.mealy import MealyHMM -from sofic.generators.mixed_state import MixedStatePresentation +from sofic.generators.mixed_state import MixedState, MixedStatePresentation from sofic.generators.prob import ( array_sum, as_prob, @@ -34,19 +34,102 @@ TransitionSignature = tuple[tuple[Any, int, Any], ...] +_SUPPORT_ATOL = 1e-8 +_LIMIT_VERIFY_ATOL = 1e-9 + def build_epsilon_machine(hmm: HiddenMarkovModel, *, max_states: int = 10_000) -> EpsilonMachine: hmm = hmm.to_mealy() - presentation = _unifilar_presentation(hmm, max_states=max_states) + try: + presentation = _unifilar_presentation(hmm, max_states=max_states) + except MixedStateExplosionError: + limit = _limit_machine(hmm, _filtered_belief(hmm), max_states=max(64, max_states // 10)) + if limit is None: + raise + return limit + machine = _minimize(presentation, constraints=getattr(hmm, "symbol_constraints", None)) + seed = _boundary_belief(presentation) + if seed is None: + return machine + limit = _limit_machine(hmm, seed, max_states=len(list(presentation.states())) - 1) + if limit is None or len(list(limit.states())) >= len(list(machine.states())): + return machine + return limit + + +def _minimize(presentation: MealyHMM, *, constraints: Any = None) -> EpsilonMachine: stationary = presentation.stationary_distribution() if is_zero(array_sum(stationary)): raise StochasticValidationError("generator must have a stationary distribution") - constraints = getattr(presentation, "symbol_constraints", None) or getattr(hmm, "symbol_constraints", None) + constraints = getattr(presentation, "symbol_constraints", None) or constraints partitions = _refine_probabilistic_partitions(presentation, constraints=constraints) return _quotient_machine(presentation, partitions, stationary, constraints=constraints) +def _boundary_belief(presentation: MealyHMM) -> np.ndarray | None: + """The heaviest recurrent belief, if the recurrent beliefs sit within round-off of the simplex boundary. + + When no finite word reveals the hidden state (e.g. the phase of a period-2 process), beliefs + converge to the causal states without reaching them. Rounding then closes the mixed-state + presentation on a reflecting boundary: a large recurrent class of beliefs with positive entries + below ``_SUPPORT_ATOL``, carrying geometrically decaying mass. Returns ``None`` for any other + presentation. + """ + states = [state for state in presentation.states() if isinstance(state, MixedState)] + if not states or any(state.is_symbolic() for state in states): + return None + beliefs = np.array([state.as_array() for state in states]) + if not ((beliefs > 0.0) & (beliefs < _SUPPORT_ATOL)).any(): + return None + stationary = np.asarray(presentation.stationary_distribution(), dtype=float) + mass = dict(zip(presentation.reindex().states, stationary, strict=True)) + return max(states, key=lambda state: mass[state]).as_array() + + +def _filtered_belief(hmm: MealyHMM, *, length: int = 4096) -> np.ndarray | None: + """Belief over ``hmm``'s states after Bayes-filtering a sampled word of ``length`` symbols.""" + from sofic.generators.matrices import symbol_matrices + + if has_symbolic(hmm.stationary_distribution().ravel()): + return None + matrices = {symbol: np.asarray(matrix, dtype=float) for symbol, matrix in symbol_matrices(hmm).items()} + symbols, _ = hmm.sample(length, rng=np.random.default_rng(0)) + belief = np.asarray(hmm.stationary_distribution(), dtype=float) + for symbol in symbols: + belief = belief @ matrices[symbol] + belief /= belief.sum() + return belief + + +def _limit_machine(hmm: MealyHMM, seed: np.ndarray | None, *, max_states: int) -> EpsilonMachine | None: + """The ε-machine generated from a limit belief, if it is finite and generates ``hmm``'s process. + + ``seed`` approximates a causal state of a process that synchronizes only asymptotically; entries + below ``_SUPPORT_ATOL`` are set to zero, and the mixed-state presentation it generates under exact + Bayes updates (at most ``max_states`` beliefs) is minimized. The result is accepted only if it + generates the same process as ``hmm`` + (:func:`~sofic.generators.minimal_generative_model._equivalence_error`): a finite unifilar + presentation of the process minimizes to its ε-machine. + """ + if seed is None: + return None + seed = np.where(seed < _SUPPORT_ATOL, 0.0, seed) + from sofic.generators.minimal_generative_model import _equivalence_error + from sofic.generators.mixed_state_construction import build_mixed_state_presentation + + try: + limit = build_mixed_state_presentation( + hmm, initial_mixed_state=MixedState.from_vector(seed), max_states=max(1, max_states) + ) + except MixedStateExplosionError: + return None + candidate = _minimize(_unifilar_presentation(limit)) + if _equivalence_error(candidate, hmm) > _LIMIT_VERIFY_ATOL: + return None + return candidate + + def _unifilar_presentation(hmm: MealyHMM, *, max_states: int = 10_000) -> MealyHMM: """Return a row-unifilar presentation generating the same process.""" if not hmm.is_unifilar(): diff --git a/tests/test_asymptotic_synchronization.py b/tests/test_asymptotic_synchronization.py new file mode 100644 index 0000000..2412bee --- /dev/null +++ b/tests/test_asymptotic_synchronization.py @@ -0,0 +1,34 @@ +"""ε-machines of processes whose hidden state no finite word reveals.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from sofic.examples.epsilon_machines import from_symbol_matrices +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.minimal_generative_model import _equivalence_error +from sofic.generators.reversal import time_reverse_stochastic + + +def _paired_bits(e: float): + """Period 2: a fair bit, then a copy of it flipped with probability ``e``; the phase is never certain.""" + mats = { + 0: np.array([[0, 0.5, 0], [1 - e, 0, 0], [e, 0, 0]]), + 1: np.array([[0, 0, 0.5], [e, 0, 0], [1 - e, 0, 0]]), + } + return from_symbol_matrices(("A", "B0", "B1"), (0, 1), mats) + + +# 0.4: rounding closes the mixed-state presentation on hundreds of near-limit beliefs; +# 0.1: the mixed-state presentation exceeds ``max_states`` first. +@pytest.mark.parametrize("e", [0.4, 0.1]) +def test_reverse_epsilon_machine_of_asymptotically_synchronizing_process(e): + forward = _paired_bits(e) + reversed_hmm = time_reverse_stochastic(forward) + + reverse = EpsilonMachine.from_hmm(reversed_hmm) + + assert len(list(reverse.states())) == 3 + assert float(reverse.statistical_complexity()) == pytest.approx(1.5) + assert _equivalence_error(reverse, reversed_hmm) < 1e-12