diff --git a/CHANGELOG.md b/CHANGELOG.md index 19e0ce9..c844fef 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,17 @@ A second correctness review: every fix below was reproduced against a brute-force reference or a closed-form value and has a regression test. +### Finite-state rounded HMM filters + +- `quantized_filter` builds the reachable log-odds-grid observer of Feng & Meng + (arXiv:2611.13760) for positive-transition Moore HMMs, with bit-valued event, + uniform-word, and total-forward-KL budgets. It preserves the exact initial + label, permits zero emissions, exposes label-capacity/error metadata, and + separates recursively rounded generation from posterior-initialized private + simulation. Bounds are theorem-derived in real arithmetic, explicitly not + interval-certified (`bounds.rigorous=False`). Tests reproduce the paper's + log-grid state counts and check finite-horizon losses against exact enumeration. + ### ε-machine construction - Processes whose hidden state no finite word reveals (e.g. a period-2 process whose phase stays diff --git a/docs/generators/generators.rst b/docs/generators/generators.rst index 54129b3..ccffc3a 100644 --- a/docs/generators/generators.rst +++ b/docs/generators/generators.rst @@ -34,6 +34,7 @@ Computational mechanics block_convergence bidirectional_epsilon_machine mixed_state_presentation + quantized_filter nonunifilar_entropy_rate edge_machine information_anatomy diff --git a/docs/generators/mixed_state_presentation.rst b/docs/generators/mixed_state_presentation.rst index 79ba4c1..0c6b2b0 100644 --- a/docs/generators/mixed_state_presentation.rst +++ b/docs/generators/mixed_state_presentation.rst @@ -10,6 +10,12 @@ probability simplices over hidden states of an HMM. The construction is the mixed-state presentation used for nonunifilar hidden processes and time-symmetric computational mechanics :cite:`Ellison2009`. +This constructor enumerates beliefs and raises ``MixedStateExplosionError`` +if they do not close within the state cap. Numeric equality tolerances are +not approximation-error budgets. For an opt-in finite approximation of a +known positive-transition Moore HMM, see :doc:`quantized_filter` +:cite:`feng2026recursive`. + Use :meth:`MixedStatePresentation.to_recurrent` to drop transient belief states and work only with the recurrent component. Pure recurrent components are returned as :class:`~sofic.generators.epsilon_machine.EpsilonMachine`; diff --git a/docs/generators/predictive_rate_distortion.rst b/docs/generators/predictive_rate_distortion.rst index eea3596..72d8bdc 100644 --- a/docs/generators/predictive_rate_distortion.rst +++ b/docs/generators/predictive_rate_distortion.rst @@ -28,6 +28,12 @@ distortion reported is the predictive information lost, Algorithm ========= +This is a static encoder from causal states to features, not a construction +of a recursively updated finite-state observer. Its information rate is not +an allocated label-capacity bound. For an explicit rounded filter of a known +positive-transition Moore HMM, with separate query and sequential-generation +losses, see :doc:`quantized_filter` :cite:`feng2026recursive`. + * **Future morphs** :math:`P(X_{0:L} = w \mid S = s) = \langle \delta_s | T^{(w)} | \mathbf{1} \rangle` are exact. Futures are built by prepending symbols to the column vectors :math:`T^{(w)}\mathbf{1}`. Words whose columns are diff --git a/docs/generators/quantized_filter.rst b/docs/generators/quantized_filter.rst new file mode 100644 index 0000000..e382de7 --- /dev/null +++ b/docs/generators/quantized_filter.rst @@ -0,0 +1,200 @@ +.. py:module:: sofic.generators.quantized_filter + +**************************************** +Finite-state quantized filtering +**************************************** + +A known HMM can have infinitely many reachable Bayesian beliefs. Instead of +trying to enumerate them exactly, :func:`quantized_filter` rounds the filter +onto a finite log-odds grid and builds only the labels reachable from its +specified initial belief. The construction follows Feng and Meng +:cite:`feng2026recursive`, Theorem 3.1, Theorem 4.3, and Appendix A.4. + +This is an opt-in approximation, not a replacement for +:doc:`mixed_state_presentation`. It does not learn the HMM, minimize the number +of states, or solve the :doc:`predictive_rate_distortion` bottleneck. + +.. warning:: + + The bounds below are theorem-derived bounds in real arithmetic. This + implementation uses float64 and checks rounding distances numerically; + ``result.bounds.rigorous`` is **False**. It is not an outward-rounded + interval certificate. Precision too fine to resolve, nonpositive rounded + beliefs, and reachable-state cap overruns raise ``ValueError``. + + Only numeric :class:`~sofic.generators.moore.MooreHMM` models with strictly + positive hidden transitions are accepted. Emissions may be zero. A positive + *summed* transition matrix of an arbitrary Mealy HMM does not imply the + required contraction, so converting an arbitrary Mealy model is not a way + to obtain these guarantees. + +Three different operations +========================== + +The observer's retained label identifies a representative predictive belief +:math:`z`. It can be used in three ways :cite:`feng2026recursive`: + +* **Event query:** supply the likelihood vector :math:`v_B(i)=P(B\mid S_t=i)` + for a future event, and decode :math:`zv_B` with ``event_probability``. + Its log error is bounded independently of the event's horizon. Constructing + the likelihood vector is the caller's responsibility. +* **Recursive generation:** ``to_generator()`` returns a fixed-start, + unifilar :class:`~sofic.generators.mealy.MealyHMM`. It emits from the + representative belief and rounds again after each output. Conditional + errors accumulate, so a fixed *total* loss requires finer precision at + longer horizons. +* **Private-state simulation:** ``private_generator(state=...)`` initializes + a copy of the original HMM from a retained belief. It samples a hidden + state once and then follows the original dynamics. Only initialization + uses the approximation; its event-law error is horizon-independent + (Proposition 4.4). The paper's :math:`K+m` label count assumes the observation + and simulation phases reuse storage. Keeping both live concurrently is a + different memory model. + +``state_after(history)`` begins at the exact initial label 0. ``step(label, +next_symbol)`` maintains the observer online while retaining only a label. +Unknown or impossible observations are rejected. Bounds concern histories +starting from the specified initial law, not arbitrary resets of the source +alone. The arrays exposed by the result are read-only, and source models are +copied so later model edits cannot change the observer's meaning. + +Accuracy and units +================== + +All public tolerances, decoded log errors, and KL losses use **bits**. The +paper instead uses nats for losses. A requested bit budget :math:`\varepsilon` +is converted internally to :math:`\epsilon=(\log 2)\varepsilon` nats. + +For the departure-emission Moore convention, let :math:`P` be the hidden +transition matrix and :math:`D_x` the diagonal emission matrix. The predictive +belief :math:`u_t=P(S_t\mid X_{0:t})` evolves as + +.. math:: + + u_{t+1}=\Psi_{X_t}(u_t),\qquad + \Psi_x(u)=\frac{uD_xP}{uD_x\mathbf{1}}. + +Positive transitions make this a contraction in the natural-log Hilbert +metric :math:`d_H`, with Birkhoff coefficient :math:`\kappa<1`: + +.. math:: + + d_H(\Psi_x(u),\Psi_x(z))\leq\kappa d_H(u,z). + +If each rounding introduces at most :math:`\eta`, then +:math:`d_H(u_t,z_t)\leq\eta/(1-\kappa)`. The constructor chooses a grid spacing +at most :math:`(1-\kappa)\xi`, where + +.. math:: + + \xi=\begin{cases} + \epsilon, & \text{event queries},\\ + \epsilon/T, & \text{uniform word log error},\\ + \log(1+\sqrt{\epsilon/T}), & \text{total forward KL}. + \end{cases} + +Use ``loss="event"`` without a horizon, or ``loss="word"`` / ``loss="kl"`` +with a positive ``horizon=T``. The latter KL bound is uniform over supplied +observation histories, which is stronger than averaging over histories. +``bounds.error_bound`` reports the bound for the requested criterion; +``word_log_error_bound(T)`` and ``total_kl_bound(T)`` can also inspect a +constructed observer at another horizon. + +With :math:`L=\log(p_{\max}/p_{\min})`, the full grid has at most + +.. math:: + + 1+\left(1+\left\lceil\frac{2L}{\eta}\right\rceil\right)^{m-1} + +labels, including the exact initial label. Only reachable labels are built; +identical transition rows have an exact one- or two-label shortcut. The +constructor honors an explicit model initial law or ``start=...``; stationarity +is used only when no initial law is supplied. ``grid_states`` reports the +unpruned construction, ``n_states`` the reachable count, and ``memory_bits`` +is :math:`\lceil\log_2 K\rceil`, **not** stationary state entropy, :math:`C_\mu`, +:math:`C_g`, or an information-bottleneck rate. + +A per-symbol KL allowance :math:`a` corresponds to total budget :math:`aT`. +Then the required grid precision need not change with :math:`T`. The horizon +cost is not a universal barrier to fixed-memory models with a fixed perplexity +tolerance (Section 6.1 of :cite:`feng2026recursive`). + +Example: queries versus generation +======================================== + +The two-regime model of the paper has hidden flip probability 0.2 and emission +noise 0.1. The following uses existing HMM word-probability and sampling APIs: + +.. code-block:: python + + import math + from itertools import product + + import numpy as np + from sofic import MooreHMM, quantized_filter + from sofic.generators.matrices import emission_tensors + + source = MooreHMM( + observation_alphabet=frozenset("01"), + initial_distribution={0: 0.5, 1: 0.5}, + ) + for state in (0, 1): + source.add_transition(state, state, 0.8) + source.add_transition(state, 1 - state, 0.2) + source.set_emission_distribution(0, {"0": 0.9, "1": 0.1}) + source.set_emission_distribution(1, {"0": 0.1, "1": 0.9}) + + # Match the paper's 0.1-nat event-error budget explicitly. + observer = quantized_filter(source, 0.1 / math.log(2)) + assert (observer.grid_states, observer.n_states, observer.memory_bits) == (72, 29, 5) + assert not observer.bounds.rigorous + history = tuple("0110100111") + label = observer.state_after(history) + + # Query the event "no alarm in the next 12 observations". + _, matrices = emission_tensors(source) + likelihood = np.linalg.matrix_power(matrices["0"], 12) @ np.ones(2) + event_probability = observer.event_probability(likelihood, state=label) + + # The same approximate future law can be sampled with private hidden states. + private = observer.private_generator(state=label) + observations, hidden_states = private.sample(100, rng=np.random.default_rng(0)) + assert np.isclose(private.word_probability("0" * 12), event_probability) + + # In contrast, this generator rounds again after every visible output. + rounded = observer.to_generator() + exact_belief = source.forward(history, normalize=True)[-1] + words = list(product("01", repeat=8)) + p = np.array([source.word_probability(w, start=exact_belief) for w in words]) + q = np.array([rounded.word_probability(w, start=label) for w in words]) + total_kl_bits = float(p @ np.log2(p / q)) + assert total_kl_bits <= observer.bounds.total_kl_bound(8) + 1e-12 + + # Choose precision for a particular total-KL forecast budget instead. + forecast = quantized_filter(source, 0.1 / math.log(2), loss="kl", horizon=100) + assert (forecast.grid_states, forecast.n_states) == (225, 51) + assert forecast.bounds.error_bound <= 0.1 / math.log(2) + 1e-12 + +These reachable-state counts reproduce the log-grid rows of Table 3 using +float64 comparisons, not the paper's interval checks. ``to_generator()`` +returns an ordinary serializable Mealy HMM; it does not claim to preserve the +source exactly or to be a minimal epsilon-machine. + +The explicit finite-word KL calculation grows exponentially in the forecast +horizon. For larger problems, the KL chain rule :cite:`Cover2006` gives the +same total by summing expected conditional KLs along **source** prefixes, +running the candidate's own updates along each prefix. Monte Carlo estimates +of those expectations are not uniform error certificates. Existing +:doc:`relative_entropy_rate` functions instead concern stationary, asymptotic +bits per symbol. + +API +=== + +.. autofunction:: quantized_filter + +.. autoclass:: QuantizedFilter + :members: + +.. autoclass:: QuantizedFilterBounds + :members: diff --git a/docs/references.bib b/docs/references.bib index beb4a3c..7befadb 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -1239,6 +1239,17 @@ @article{riechers2019fraudulent doi = {10.1103/PhysRevResearch.3.013170}, } +@article{feng2026recursive, + author = {Feng, Pengbin and Meng, Chunlei}, + title = {Recursive Memory for Prediction and Generation}, + journal = {arXiv preprint arXiv:2611.13760}, + year = {2026}, + eprint = {2611.13760}, + archivePrefix = {arXiv}, + primaryClass = {stat.AP}, + url = {https://arxiv.org/abs/2611.13760}, +} + @article{Marzen2016, author = {Marzen, Sarah E. and Crutchfield, James P.}, title = {Predictive Rate-Distortion for Infinite-Order {Markov} Processes}, diff --git a/sofic/__init__.py b/sofic/__init__.py index 6b86cf4..68cff62 100644 --- a/sofic/__init__.py +++ b/sofic/__init__.py @@ -73,6 +73,8 @@ MooreHMM, NMachine, ProbabilisticFiniteAutomaton, + QuantizedFilter, + QuantizedFilterBounds, QuasiRealization, QuasiStochasticModel, StochasticModel, @@ -82,6 +84,7 @@ is_lumpable, lump, minimal_generative_model, + quantized_filter, wheeler_presentation, wheeler_statistical_complexity, wyner_generative_model, @@ -146,6 +149,9 @@ "NMachine", "ObservationTable", "ProbabilisticFiniteAutomaton", + "QuantizedFilter", + "QuantizedFilterBounds", + "quantized_filter", "QuasiRealization", "QuasiStochasticModel", "RegularLanguage", diff --git a/sofic/generators/__init__.py b/sofic/generators/__init__.py index ec18a96..914d1dc 100644 --- a/sofic/generators/__init__.py +++ b/sofic/generators/__init__.py @@ -49,6 +49,7 @@ from sofic.generators.nmachine import NMachine from sofic.generators.pfa import ProbabilisticFiniteAutomaton from sofic.generators.predictive_rd import PredictiveRateDistortionCurve, predictive_rate_distortion +from sofic.generators.quantized_filter import QuantizedFilter, QuantizedFilterBounds, quantized_filter from sofic.generators.quasi_realization import QuasiRealization from sofic.generators.relative_entropy_rate import ( RelativeEntropyRateBounds, @@ -119,6 +120,9 @@ "MooreHMM", "NMachine", "ProbabilisticFiniteAutomaton", + "QuantizedFilter", + "QuantizedFilterBounds", + "quantized_filter", "QuasiRealization", "QuasiStochasticModel", "StochasticModel", diff --git a/sofic/generators/quantized_filter.py b/sofic/generators/quantized_filter.py new file mode 100644 index 0000000..435810c --- /dev/null +++ b/sofic/generators/quantized_filter.py @@ -0,0 +1,388 @@ +r"""Finite-state rounded Bayesian filters for positive-transition Moore HMMs. + +The log-odds grid and loss-dependent precision follow Feng & Meng, +:cite:`feng2026recursive`, Theorems 3.1 and 4.3 and Appendix A.4. Bounds are +analytic real-arithmetic bounds; this float64 implementation does not provide +an outward-rounded numerical certificate. +""" + +from __future__ import annotations + +import math +from collections.abc import Hashable, Iterable, Sequence +from dataclasses import dataclass, field +from numbers import Integral +from typing import Any, Literal + +import numpy as np + +from sofic.generators.matrices import StartSpec, start_vector +from sofic.generators.mealy import MealyHMM +from sofic.generators.moore import MooreHMM +from sofic.graph import ATTR_EMISSION_DIST +from sofic.properties import transition_matrix + +__all__ = ["QuantizedFilter", "QuantizedFilterBounds", "quantized_filter"] + +_LOG2 = math.log(2.0) +Loss = Literal["event", "word", "kl"] + + +def _positive_integer(value: object, name: str) -> int: + if isinstance(value, bool) or not isinstance(value, Integral) or value < 1: + raise ValueError(f"{name} must be a positive integer") + return int(value) + + +@dataclass(frozen=True) +class QuantizedFilterBounds: + r"""Theorem-derived error bounds for a rounded filter :cite:`feng2026recursive`. + + ``tolerance`` and all decoded losses are in **bits**. ``rounding_error`` + and ``hilbert_error`` are natural-log Hilbert distances; ``contraction`` + is Birkhoff's coefficient. ``max_rounding_error`` is the largest float64 + distance checked over the reachable transitions, not an interval bound. + + ``rigorous`` is always false: the theorems hold in real arithmetic, but + grid rounding, normalization, and the transition checks use float64. + This record is not a machine-verified certificate or an optimal-memory + claim. Bounds hold along histories starting from the constructor's exact + initial label, not after an arbitrary reset of the source alone. + """ + + loss: Loss + tolerance: float + horizon: int | None + contraction: float + rounding_error: float + hilbert_error: float + max_rounding_error: float + rigorous: bool = field(default=False, init=False) + + @property + def event_log_error_bound(self) -> float: + """Absolute log2 error for any supplied positive-probability future event.""" + return self.hilbert_error / _LOG2 + + def word_log_error_bound(self, horizon: int) -> float: + """Uniform absolute log2 word error through ``horizon`` generated symbols.""" + return _positive_integer(horizon, "horizon") * self.event_log_error_bound + + def total_kl_bound(self, horizon: int) -> float: + """Total forward KL in bits, uniformly over observed histories. + + Each conditional KL is at most both the log-ratio bound and the + chi-square bound ``expm1(hilbert_error)**2`` (Theorem 4.3). + """ + horizon = _positive_integer(horizon, "horizon") + error = self.hilbert_error + one_step = min(error, math.expm1(error) ** 2) if error < 1.0 else error + return horizon * (one_step / _LOG2) + + @property + def error_bound(self) -> float: + """Bound for the loss and horizon requested from :func:`quantized_filter`.""" + if self.loss == "event": + return self.event_log_error_bound + assert self.horizon is not None + if self.loss == "word": + return self.word_log_error_bound(self.horizon) + return self.total_kl_bound(self.horizon) + + +@dataclass(frozen=True) +class QuantizedFilter: + r"""Reachable finite-state observer and its approximate beliefs. + + Construct with :func:`quantized_filter`, following :cite:`feng2026recursive`. + Label 0 stores the exact initial belief. ``beliefs[z]`` is a distribution + over ``basis_states``; ``successors[z, a]`` is the next label on + ``symbols[a]`` or -1 for an impossible observation. Arrays are read-only. + ``grid_states`` counts the full construction before reachable-state + pruning, including the separate initial label. + + Event queries and private simulation use the original HMM from a retained + belief. :meth:`to_generator` instead rounds again after each output, so its + errors can accumulate. These are deliberately separate operations. + """ + + basis_states: tuple[Hashable, ...] + symbols: tuple[Any, ...] + beliefs: np.ndarray + successors: np.ndarray + output_probabilities: np.ndarray + grid_states: int + bounds: QuantizedFilterBounds + _source: MooreHMM = field(repr=False, compare=False) + + @property + def n_states(self) -> int: + """Number of reachable labels, including the exact initial label.""" + return len(self.beliefs) + + @property + def memory_bits(self) -> int: + """Allocated label capacity ``ceil(log2(n_states))``, not state entropy.""" + return (self.n_states - 1).bit_length() + + def _state(self, state: int) -> int: + if isinstance(state, bool) or not isinstance(state, Integral) or not 0 <= state < self.n_states: + raise ValueError(f"state must be a label in [0, {self.n_states})") + return int(state) + + def step(self, state: int, symbol: Any) -> int: + """Update one retained label; reject unknown or impossible observations.""" + state = self._state(state) + try: + column = self.symbols.index(symbol) + except ValueError as error: + raise ValueError(f"unknown observation {symbol!r}") from error + target = int(self.successors[state, column]) + if target < 0: + raise ValueError(f"observation {symbol!r} is impossible from state {state}") + return target + + def state_after(self, observations: Iterable[Any], *, state: int = 0) -> int: + """Read observations from label 0, or continue a previous observation phase.""" + state = self._state(state) + for symbol in observations: + state = self.step(state, symbol) + return state + + def event_probability(self, likelihood: Sequence[float], *, state: int = 0) -> float: + r"""Decode ``beliefs[state] @ likelihood`` for a supplied future event. + + ``likelihood[i]`` must be the event's probability given hidden state + ``basis_states[i]`` at the next emission. Theorem 3.1 of + :cite:`feng2026recursive` bounds its relative/log error at every horizon; + constructing the event-likelihood vector is the caller's responsibility. + """ + state = self._state(state) + vector = np.asarray(likelihood, dtype=float) + if vector.shape != (len(self.basis_states),) or not np.all(np.isfinite(vector)): + raise ValueError("likelihood must have one finite probability per hidden state") + if np.any(vector < 0.0) or np.any(vector > 1.0): + raise ValueError("likelihood values must lie in [0, 1]") + return float(self.beliefs[state] @ vector) + + def to_generator(self) -> MealyHMM: + r"""Return the rounded, fixed-start unifilar generator (Appendix A.4). + + This is an approximate generator, not an epsilon-machine or an exact + mixed-state presentation :cite:`feng2026recursive`. Returned models + can be sampled, queried, or serialized with the ordinary HMM APIs. + """ + generator = MealyHMM(initial_distribution={0: 1.0}, observation_alphabet=frozenset(self.symbols)) + for state, row in enumerate(self.successors): + generator.graph.add_state(state) + for column, target in enumerate(row): + if target >= 0: + generator.add_transition( + state, int(target), self.symbols[column], float(self.output_probabilities[state, column]) + ) + return generator + + def private_generator(self, *, state: int = 0) -> MooreHMM: + r"""Initialize a copy of the original hidden generator from a retained belief. + + Only initialization uses the approximation (Proposition 4.4 of + :cite:`feng2026recursive`). Its event-law error is horizon-independent. + The paper's ``K + m`` label accounting assumes observation and simulation + phases reuse storage; keeping this observer live concurrently is not + that memory model. Subsequent edits to the copy cannot alter the filter. + """ + state = self._state(state) + generator = self._source.copy() + generator.initial_distribution = dict(zip(self.basis_states, self.beliefs[state].tolist(), strict=True)) + return generator + + +def _model_data(model: MooreHMM, start: StartSpec) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + model.validate() + states = model.reindex().states + if not states: + raise ValueError("quantized_filter requires a nonempty MooreHMM") + transition, _states = transition_matrix(model, states=states) + transition = np.asarray(transition, dtype=float) + if not np.all(np.isfinite(transition)) or np.any(transition <= 0.0): + raise ValueError("quantized_filter requires strictly positive hidden transition probabilities") + symbols = tuple(sorted(model.observation_alphabet, key=repr)) + emissions = np.array( + [[model.graph.state_attrs(s).get(ATTR_EMISSION_DIST, {}).get(y, 0.0) for y in symbols] for s in states], + dtype=float, + ) + initial = np.asarray(start_vector(model, start), dtype=float) + for name, array in (("transition", transition), ("emission", emissions), ("initial", initial[None, :])): + if ( + not np.all(np.isfinite(array)) + or np.any(array < 0.0) + or not np.allclose(array.sum(axis=1), 1.0, rtol=0.0, atol=1e-12) + ): + raise ValueError(f"{name} probabilities must be finite, nonnegative, and normalized") + return transition, emissions, initial + + +def _precision(tolerance: float, loss: Loss, horizon: int | None) -> tuple[float, int | None]: + if isinstance(tolerance, bool) or not math.isfinite(tolerance) or tolerance <= 0.0: + raise ValueError("tolerance must be finite and strictly positive (in bits)") + if loss not in ("event", "word", "kl"): + raise ValueError("loss must be 'event', 'word', or 'kl'") + if loss == "event": + if horizon is not None: + raise ValueError("event loss is horizon-independent; omit horizon") + return tolerance * _LOG2, None + horizon = _positive_integer(horizon, "horizon") + log_budget = math.log(tolerance) + math.log(_LOG2) - math.log(horizon) + if loss == "word": + return math.exp(log_budget), horizon + # log(1 + sqrt(epsilon_nats / T)), stable for very large horizons. + return float(np.logaddexp(0.0, log_budget / 2.0)), horizon + + +def quantized_filter( + model: MooreHMM, + tolerance: float, + *, + loss: Loss = "event", + horizon: int | None = None, + start: StartSpec = None, + max_states: int = 10_000, +) -> QuantizedFilter: + r"""Build a finite rounded filter for a known positive-transition Moore HMM. + + Implements the log-odds construction of :cite:`feng2026recursive`, Theorem + 3.1 (event queries), Appendix A.4 (uniform words), and Theorem 4.3 (total + forward KL). No stationary initialization, positive emissions, or finite + exact mixed-state presentation is required. Arbitrary Mealy HMMs are not + accepted: positivity of their summed transition matrix is insufficient. + + Parameters + ---------- + model + Numeric state-emitting HMM with strictly positive hidden transitions. + The departure-emission convention is that of ``MooreHMM.to_mealy``. + tolerance + Positive error budget in **bits** (log2 probability error or KL). + loss + ``'event'``: uniform future-event log error, independent of horizon. + ``'word'``: uniform word log error for the rounded generator through + ``horizon``. ``'kl'``: total forward KL for that generator through + ``horizon``, uniformly over observed histories, not just on average. + horizon + Positive integer required for ``'word'`` and ``'kl'``; omitted for + ``'event'``. A per-symbol budget ``a`` corresponds to total ``a*horizon``. + start + Initial law in the usual HMM start-vector convention. Defaults to the + model's initial law (stationary only when the model has no explicit law). + Label 0 stores it exactly, including boundary beliefs. + max_states + Safety cap on reachable labels. The full Cartesian grid is never + allocated, but its size is exponential in hidden-belief dimension. + + Notes + ----- + With ``L = log(p_max/p_min)`` and Birkhoff coefficient ``kappa``, round + log-odds coordinates at spacing at most ``(1-kappa)*xi``. Here ``xi`` is + ``epsilon``, ``epsilon/T``, or ``log(1+sqrt(epsilon/T))`` for the three + losses, after converting the bit budget to nats. Hilbert contraction + bounds repeated rounding, and mixture likelihood ratios control events. + + The returned bounds describe this construction in real arithmetic; + ``bounds.rigorous`` is false because numerical operations are float64, + not interval-certified. Unresolvable grids or nonpositive rounded beliefs + are rejected rather than given a spurious bound. This is not an optimizer + of allocated memory, statistical complexity, or the information bottleneck. + + Raises + ------ + TypeError + If ``model`` is not a Moore HMM. + ValueError + For invalid probabilities/options, a reachable-state cap, or precision + too fine for float64. + """ + if not isinstance(model, MooreHMM): + raise TypeError("quantized_filter requires a MooreHMM, not a general edge-emitting HMM") + precision, horizon = _precision(tolerance, loss, horizon) + max_states = _positive_integer(max_states, "max_states") + source = model.copy() + transition, emissions, initial = _model_data(source, start) + basis = source.reindex().states + symbols = tuple(sorted(source.observation_alphabet, key=repr)) + constant = bool(np.all(transition == transition[0])) + log_p = np.log(transition) + diameter = max(float(np.ptp(log_p - row, axis=1).max()) for row in log_p) + contraction = math.tanh(diameter / 4.0) + gap = 1.0 - contraction + width = float(log_p.max() - log_p.min()) + if constant: + intervals = 0 + spacing = 0.0 + grid_states = 1 if np.array_equal(initial, transition[0]) else 2 + else: + radius = gap * precision + ratio = 2.0 * width / radius if radius > 0.0 else math.inf + if not math.isfinite(ratio) or ratio > 2**48: + raise ValueError("requested filter precision is too fine for float64") + intervals = max(1, math.ceil(ratio)) + spacing = 2.0 * width / intervals + grid_states = 1 + (intervals + 1) ** (len(basis) - 1) + + beliefs = [initial] + labels: dict[tuple[int, ...], int] = {} + if constant and np.array_equal(initial, transition[0]): + labels[()] = 0 + successors: list[list[int]] = [] + probabilities: list[np.ndarray] = [] + max_rounding_error = 0.0 + cursor = 0 + while cursor < len(beliefs): + belief = beliefs[cursor] + output = belief @ emissions + row = [-1] * len(symbols) + for column, probability in enumerate(output): + likelihood = emissions[:, column] + possible = bool(np.any((belief > 0.0) & (likelihood > 0.0))) + if not possible: + continue + if probability <= 0.0: + raise ValueError("observation probability underflows in float64") + # Scaling the emission law avoids underflow for globally rare symbols. + posterior = belief * (likelihood / likelihood.max()) + posterior /= posterior.sum() + updated = posterior @ transition + updated /= updated.sum() + if constant: + key = () + rounded = transition[0] + else: + coordinates = np.log(updated[:-1]) - math.log(updated[-1]) + position = (coordinates + width) / (2.0 * width) * intervals + key = tuple(int(i) for i in np.clip(np.floor(position + 0.5), 0, intervals)) + coordinates = np.r_[-width + (2.0 * width / intervals) * np.asarray(key), 0.0] + rounded = np.exp(coordinates - coordinates.max()) + rounded /= rounded.sum() + if np.any(rounded <= 0.0): + raise ValueError("rounded belief underflows in float64") + error = float(np.ptp(np.log(rounded) - np.log(updated))) + max_rounding_error = max(max_rounding_error, error) + if error > spacing: + raise ValueError("grid rounding cannot resolve the requested precision in float64") + target = labels.get(key) + if target is None: + if len(beliefs) >= max_states: + raise ValueError(f"quantized filter exceeds max_states={max_states}; increase the cap or tolerance") + target = labels[key] = len(beliefs) + beliefs.append(rounded) + row[column] = target + successors.append(row) + probabilities.append(output) + cursor += 1 + + arrays = [np.array(beliefs), np.array(successors, dtype=int), np.array(probabilities)] + for array in arrays: + array.setflags(write=False) + bounds = QuantizedFilterBounds( + loss, float(tolerance), horizon, contraction, spacing, spacing / gap, max_rounding_error + ) + return QuantizedFilter(basis, symbols, arrays[0], arrays[1], arrays[2], grid_states, bounds, source) diff --git a/tests/test_quantized_filter.py b/tests/test_quantized_filter.py new file mode 100644 index 0000000..192a4da --- /dev/null +++ b/tests/test_quantized_filter.py @@ -0,0 +1,345 @@ +"""Rounded filters and loss guarantees of Feng & Meng, arXiv:2611.13760.""" + +from __future__ import annotations + +import math +from itertools import product + +import numpy as np +import pytest + +from sofic import MooreHMM, QuantizedFilter, QuantizedFilterBounds, quantized_filter +from sofic.exceptions import StochasticValidationError +from sofic.generators.matrices import emission_tensors + + +def _moore(transition=None, emissions=None, initial=(0.5, 0.5)): + transition = np.array([[0.8, 0.2], [0.2, 0.8]] if transition is None else transition, dtype=float) + emissions = np.array([[0.9, 0.1], [0.1, 0.9]] if emissions is None else emissions, dtype=float) + symbols = tuple(str(i) for i in range(emissions.shape[1])) + model = MooreHMM(observation_alphabet=frozenset(symbols)) + if initial is not None: + model.initial_distribution = dict(enumerate(initial)) + for i, row in enumerate(transition): + for j, probability in enumerate(row): + if probability: + model.add_transition(i, j, float(probability)) + model.set_emission_distribution(i, dict(zip(symbols, emissions[i], strict=True))) + return model + + +def _words(max_length): + for length in range(max_length + 1): + yield from product("01", repeat=length) + + +def _likelihood(matrices, word): + vector = np.ones(len(next(iter(matrices.values())))) + for symbol in reversed(word): + vector = matrices[symbol] @ vector + return vector + + +def _forecast_laws(model, observer, history, horizon): + """Independent enumeration of the two laws and the source-prefix KL chain rule.""" + _, matrices = emission_tensors(model) + initial = model.forward(history, normalize=True)[-1] + nodes = [(initial, 1.0, observer.state_after(history))] + chain_kl = 0.0 + for _ in range(horizon): + following = [] + for mass, q_mass, state in nodes: + p_mass = float(mass.sum()) + conditional = np.array([(mass @ matrices[y]).sum() / p_mass for y in observer.symbols]) + q_row = observer.output_probabilities[state] + chain_kl += p_mass * float(np.sum(conditional * np.log2(conditional / q_row))) + for column, symbol in enumerate(observer.symbols): + following.append((mass @ matrices[symbol], q_mass * q_row[column], observer.successors[state, column])) + nodes = following + p = np.array([mass.sum() for mass, _, _ in nodes]) + q = np.array([mass for _, mass, _ in nodes]) + return p, q, chain_kl + + +def test_public_exports(): + from sofic.generators import QuantizedFilter as Filter + from sofic.generators import QuantizedFilterBounds as Bounds + from sofic.generators import quantized_filter as construct + + assert (Filter, Bounds, construct) == (QuantizedFilter, QuantizedFilterBounds, quantized_filter) + + +@pytest.mark.parametrize( + ("loss", "horizon", "grid_states", "reachable"), + [("event", None, 72, 29), ("kl", 100, 225, 51), ("kl", 10000, 2198, 183)], +) +def test_paper_table_3_log_grid_counts(loss, horizon, grid_states, reachable): + # The paper's epsilon=0.1 is in nats; the public API uses bits. + observer = quantized_filter(_moore(), 0.1 / math.log(2), loss=loss, horizon=horizon) + assert observer.grid_states == grid_states + assert observer.n_states == reachable + assert observer.memory_bits == (reachable - 1).bit_length() + assert observer.bounds.error_bound <= observer.bounds.tolerance * (1 + 1e-12) + assert observer.bounds.max_rounding_error <= observer.bounds.rounding_error + assert observer.bounds.contraction == pytest.approx(0.6) + assert not observer.bounds.rigorous + + +@pytest.mark.parametrize("loss", ["event", "word", "kl"]) +def test_finite_horizon_guarantees_and_kl_chain_rule(loss): + model = _moore() + horizon = 8 + observer = quantized_filter(model, 0.1, loss=loss, horizon=None if loss == "event" else horizon) + for history in [(), ("1",), tuple("011010")]: + p, q, chain_kl = _forecast_laws(model, observer, history, horizon) + log_ratios = np.log2(p / q) + direct_kl = float(p @ log_ratios) + assert direct_kl == pytest.approx(chain_kl, abs=2e-14) + assert direct_kl <= observer.bounds.total_kl_bound(horizon) + 2e-14 + assert np.max(np.abs(log_ratios)) <= observer.bounds.word_log_error_bound(horizon) + 2e-14 + assert observer.bounds.error_bound <= 0.1 * (1 + 1e-12) + + +@pytest.mark.parametrize("n_states", [2, 3, 4]) +def test_hilbert_and_event_bounds_after_all_short_histories(n_states): + rng = np.random.default_rng(n_states) + transition = rng.uniform(0.2, 1.0, (n_states, n_states)) + transition /= transition.sum(axis=1, keepdims=True) + emissions = rng.uniform(0.2, 1.0, (n_states, 2)) + emissions /= emissions.sum(axis=1, keepdims=True) + model = _moore(transition, emissions, initial=np.arange(1, n_states + 1) / sum(range(1, n_states + 1))) + observer = quantized_filter(model, 0.5) + _, matrices = emission_tensors(model) + likelihoods = [_likelihood(matrices, word) for word in product("01", repeat=5)] + for history in _words(5): + exact = model.forward(history, normalize=True)[-1] + state = observer.state_after(history) + retained = observer.beliefs[state] + distance = np.ptp(np.log(exact) - np.log(retained)) + assert distance <= observer.bounds.hilbert_error + 2e-14 + for likelihood in likelihoods: + actual = observer.event_probability(likelihood, state=state) + expected = float(exact @ likelihood) + assert abs(math.log2(actual / expected)) <= observer.bounds.event_log_error_bound + 2e-14 + + +def test_query_and_private_simulation_do_not_accumulate_rounding(): + model = _moore() + observer = quantized_filter(model, 0.5) + history = tuple("0110100111") + state = observer.state_after(history) + exact = model.forward(history, normalize=True)[-1] + private = observer.private_generator(state=state) + rounded = observer.to_generator() + word = tuple("1" * 1024) + true_log = model.log_word_probability(word, start=exact) + private_error = abs(private.log_word_probability(word) - true_log) + rounded_error = abs(rounded.log_word_probability(word, start=state) - true_log) + assert private_error <= observer.bounds.event_log_error_bound + 1e-11 + assert rounded_error > 10 * observer.bounds.event_log_error_bound + _, matrices = emission_tensors(model) + likelihood = _likelihood(matrices, tuple("10110")) + assert private.word_probability(tuple("10110")) == pytest.approx( + observer.event_probability(likelihood, state=state) + ) + observations, hidden = private.sample(20, np.random.default_rng(7)) + assert len(observations) == len(hidden) == 20 + + +def test_generator_is_fixed_start_unifilar_and_serializable(): + model = _moore() + observer = quantized_filter(model, 0.2, loss="word", horizon=20) + generator = observer.to_generator() + generator.validate() + assert generator.is_unifilar() + assert generator.initial_distribution == {0: 1.0} + restored = type(generator).from_yaml(generator.to_yaml()) + for word in _words(4): + state, probability = 0, 1.0 + for symbol in word: + probability *= observer.output_probabilities[state, observer.symbols.index(symbol)] + state = observer.step(state, symbol) + assert generator.word_probability(word) == pytest.approx(probability) + assert restored.word_probability(word) == pytest.approx(probability) + + +def test_initial_law_is_exact_and_explicit_start_overrides_model(): + model = _moore(initial=(0.1, 0.9)) + observer = quantized_filter(model, 0.1) + np.testing.assert_array_equal(observer.beliefs[0], [0.1, 0.9]) + np.testing.assert_array_equal(observer.private_generator().forward([], normalize=True)[0], [0.1, 0.9]) + boundary = quantized_filter(model, 0.1, start=0) + np.testing.assert_array_equal(boundary.beliefs[0], [1.0, 0.0]) + vector = quantized_filter(model, 0.1, start={0: 0.7, 1: 0.3}) + np.testing.assert_array_equal(vector.beliefs[0], [0.7, 0.3]) + assert observer.state_after("0110") == observer.state_after("10", state=observer.state_after("01")) + + +def test_unspecified_initial_uses_stationary_law(): + model = _moore([[0.7, 0.3], [0.1, 0.9]], initial=None) + observer = quantized_filter(model, 0.5) + np.testing.assert_allclose(observer.beliefs[0], [0.25, 0.75]) + + +@pytest.mark.parametrize("initial", [(0.25, 0.75), (1.0, 0.0)]) +def test_identical_transition_rows_are_exact(initial): + model = _moore([[0.25, 0.75], [0.25, 0.75]], initial=initial) + observer = quantized_filter(model, 0.1, loss="word", horizon=100) + expected_states = 1 if initial == (0.25, 0.75) else 2 + assert observer.n_states == observer.grid_states == expected_states + assert observer.bounds.error_bound == observer.bounds.hilbert_error == 0 + for word in _words(4): + assert observer.to_generator().word_probability(word) == pytest.approx(model.word_probability(word)) + + +def test_single_hidden_state(): + observer = quantized_filter(_moore([[1.0]], [[0.3, 0.7]], initial=(1.0,)), 0.1, max_states=1) + assert observer.n_states == observer.grid_states == 1 + assert observer.memory_bits == 0 + assert observer.bounds.error_bound == 0 + + +def test_zero_emissions_and_impossible_observations(): + model = _moore(emissions=[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], initial=(1.0, 0.0)) + observer = quantized_filter(model, 0.2) + assert observer.event_probability([0.0, 0.0]) == 0 + with pytest.raises(ValueError, match="impossible"): + observer.state_after("1") + for label in range(observer.n_states): + with pytest.raises(ValueError, match="impossible"): + observer.step(label, "2") + assert observer.state_after("010") >= 0 + generator = observer.to_generator() + generator.validate() + assert generator.word_probability("1") == 0 + assert generator.word_probability("02") == 0 + + +def test_globally_rare_symbol_does_not_underflow_the_update(): + observer = quantized_filter(_moore(emissions=[[1 - 1e-300, 1e-300], [1 - 2e-300, 2e-300]]), 0.2) + state = observer.state_after("11111") + assert np.all(np.isfinite(observer.beliefs[state])) + assert np.all(observer.beliefs[state] > 0) + + +def test_per_symbol_kl_allowance_has_horizon_independent_grid(): + short = quantized_filter(_moore(), 0.001 * 10, loss="kl", horizon=10) + long = quantized_filter(_moore(), 0.001 * 1000, loss="kl", horizon=1000) + assert short.grid_states == long.grid_states + np.testing.assert_array_equal(short.successors, long.successors) + np.testing.assert_allclose(short.beliefs, long.beliefs) + + +def test_large_tolerance_uses_finite_linear_kl_bound(): + observer = quantized_filter(_moore(), 10.0) + assert observer.bounds.hilbert_error > 1.0 + assert observer.bounds.total_kl_bound(10) == observer.bounds.word_log_error_bound(10) + assert observer.bounds.error_bound <= 10.0 + + +def test_underflowing_positive_observation_is_rejected(): + model = _moore(emissions=[[1e-300, 1.0], [0.0, 1.0]], initial=(1e-300, 1.0)) + with pytest.raises(ValueError, match="underflows"): + quantized_filter(model, 0.2) + + +def test_source_and_returned_models_do_not_mutate_observer(): + model = _moore() + observer = quantized_filter(model, 0.2) + expected = observer.private_generator().word_probability("011") + model.set_emission_distribution(0, {"0": 1.0, "1": 0.0}) + private = observer.private_generator() + private.set_emission_distribution(1, {"0": 1.0, "1": 0.0}) + observer.to_generator().initial_distribution = {1: 1.0} + assert observer.private_generator().word_probability("011") == expected + for array in (observer.beliefs, observer.successors, observer.output_probabilities): + assert not array.flags.writeable + + +@pytest.mark.parametrize("tolerance", [0.0, -1.0, math.nan, math.inf, True]) +def test_invalid_tolerance(tolerance): + with pytest.raises(ValueError, match="tolerance"): + quantized_filter(_moore(), tolerance) + + +@pytest.mark.parametrize( + "options", + [ + {"loss": "unknown"}, + {"horizon": 1}, + {"loss": "kl"}, + {"loss": "word", "horizon": 0}, + {"loss": "kl", "horizon": True}, + {"loss": "kl", "horizon": 1.5}, + {"max_states": 0}, + ], +) +def test_invalid_options(options): + with pytest.raises(ValueError): + quantized_filter(_moore(), 0.1, **options) + + +def test_reachable_cap_does_not_allocate_full_grid(): + transition = 0.4 * np.eye(5) + 0.6 / 5 + emissions = np.column_stack((np.linspace(0.1, 0.9, 5), np.linspace(0.9, 0.1, 5))) + with pytest.raises(ValueError, match="max_states=3"): + quantized_filter(_moore(transition, emissions, initial=np.ones(5) / 5), 0.01, max_states=3) + + +def test_general_mealy_is_rejected_even_with_positive_total_transition(): + with pytest.raises(TypeError, match="MooreHMM"): + quantized_filter(_moore().to_mealy(), 0.1) + + +@pytest.mark.parametrize("transition", [[[1.0, 0.0], [0.2, 0.8]], [[math.nan, 0.2], [0.2, 0.8]]]) +def test_invalid_hidden_transitions(transition): + with pytest.raises((ValueError, StochasticValidationError)): + quantized_filter(_moore(transition), 0.1) + + +@pytest.mark.parametrize("emissions", [[[0.7, 0.2], [0.2, 0.8]], [[math.nan, 0.1], [0.2, 0.8]]]) +def test_invalid_emissions(emissions): + with pytest.raises((ValueError, StochasticValidationError)): + quantized_filter(_moore(emissions=emissions), 0.1) + + +def test_empty_model_and_invalid_start(): + with pytest.raises(ValueError, match="nonempty"): + quantized_filter(MooreHMM(), 0.1) + with pytest.raises(ValueError, match="initial"): + quantized_filter(_moore(), 0.1, start=[0.2, 0.3]) + + +@pytest.mark.parametrize("state", [-1, 100000, 1.5, True]) +def test_invalid_state(state): + observer = quantized_filter(_moore(), 0.5) + with pytest.raises(ValueError, match="state"): + observer.state_after([], state=state) + with pytest.raises(ValueError, match="state"): + observer.private_generator(state=state) + + +@pytest.mark.parametrize("likelihood", [[0.3], [0.2, math.inf], [-0.1, 0.5], [1.1, 0.5]]) +def test_invalid_event_likelihood(likelihood): + with pytest.raises(ValueError, match="likelihood"): + quantized_filter(_moore(), 0.5).event_probability(likelihood) + + +def test_unknown_symbol_and_invalid_bound_horizon(): + observer = quantized_filter(_moore(), 0.5) + with pytest.raises(ValueError, match="unknown"): + observer.step(0, "unknown") + for method in (observer.bounds.word_log_error_bound, observer.bounds.total_kl_bound): + with pytest.raises(ValueError, match="horizon"): + method(0) + + +@pytest.mark.parametrize("options", [{}, {"loss": "word", "horizon": 10**100}, {"loss": "kl", "horizon": 10**100}]) +def test_unresolvable_precision_is_rejected(options): + with pytest.raises(ValueError, match="float64"): + quantized_filter(_moore(), 1e-20, **options) + + +def test_numerically_unresolvable_contraction_is_rejected(): + with pytest.raises(ValueError, match="float64"): + quantized_filter(_moore([[1.0, 1e-100], [1e-100, 1.0]]), 0.1)