Capacity bounds for epsilon-transducers; trapdoor and Ising example channels

CHANGELOG.md

@@ -131,6 +131,16 @@ brute-force reference or a closed-form value and has a regression test.
 
 ### New
 
+- `sofic.generators.capacity`: capacity bounds for finite ε-transducers
+  (unifilar finite-state channels). `EpsilonTransducer.feedback_capacity()`
+  brackets the feedback capacity with Q-graph bounds (convex upper bound,
+  certified graph-based-encoder lower bound), the Q-graph duality bound, or
+  belief-state value iteration; `EpsilonTransducer.capacity()` brackets the
+  capacity without feedback with an optimized Markov input and the Q-graph
+  duality dynamic program. Results are `CapacityBounds(lower, upper, ...)`.
+  The optional `capacity` extra installs cvxpy for the convex step.
+- `trapdoor()` / `trapdoor_et()` and `ising_channel()` / `ising_channel_et()`
+  example channels, each with a uniform initial state.
 - `sofic.inference.Trials`: `learn_epsilon_machine_cssr`,
   `learn_epsilon_transducer_cssr`, and `suggest_max_history` learn from
   independent trials, pooling counts without any history spanning two trials

docs/examples.rst

@@ -230,7 +230,14 @@ parametrized tests and sweeps:
 
 A parallel set of transducer factories (``bit_flip``, ``parity``, ``delay``,
 ``binary_channel``, …) lives alongside the processes and produces
-:class:`~sofic.automata.transducers.MealyMachine` instances.
+:class:`~sofic.automata.transducers.MealyMachine` instances. Two of them are
+classic channels with memory whose feedback capacities are known:
+``trapdoor`` (Blackwell's trapdoor channel, :math:`C_{\mathrm{fb}} = \log_2\varphi`
+:cite:`Permuter2008`) and ``ising_channel`` (:math:`C_{\mathrm{fb}} \approx 0.575522`
+:cite:`Elishco2014`), both with a uniform initial state. ``trapdoor_et`` and
+``ising_channel_et`` return their two-state
+:class:`~sofic.generators.epsilon_transducer.EpsilonTransducer` presentations, which
+feed the :doc:`channel-capacity bounds <generators/channel_capacity>`.
 
 Symbolic-shift examples (:mod:`sofic.examples.shifts`) provide the
 sofic-Dyck factories listed above; access them via

docs/generators/channel_capacity.rst

@@ -0,0 +1,144 @@
+.. channel_capacity.rst
+.. py:module:: sofic.generators.capacity
+
+****************
+Channel Capacity
+****************
+
+A finite :doc:`ε-transducer <epsilon_transducer>` is a *unifilar finite-state
+channel*: the pair ``(x, y)`` observed in causal state ``s`` determines the next
+causal state. :mod:`sofic.generators.capacity` brackets the capacity of such a
+channel, with feedback (the encoder sees past outputs) and without.
+
+Every method returns a :class:`CapacityBounds` ``lower <= C <= upper`` in bits
+per channel use. Exact answers are not available in general: the feedback
+capacity of finite-state channels is not Turing-computable
+:cite:`Grigorescu2024`, and the capacity without feedback is open even for
+the trapdoor channel. ``CapacityBounds.value`` returns the midpoint once the
+bracket is narrower than ``tol``.
+
+The methods assume a strongly connected channel (checked), so the capacity
+does not depend on the initial state. Inputs undefined at a state act as input
+constraints for the feedback bounds; the no-feedback bounds require every
+input to be admissible everywhere.
+
+Feedback capacity
+=================
+
+With feedback, the encoder knows the causal state, and the feedback capacity
+is the limit of normalized directed information :cite:`Permuter2008,Yang2005`.
+
+Q-graphs
+--------
+
+A :class:`QGraph` maps output histories onto a finite set of contexts by
+walking a deterministic graph with one edge per output symbol
+:cite:`Sabag2017`. For any Q-graph,
+
+.. math::
+
+   C_{\mathrm{fb}} \le \max_{P_{X \mid S, Q}} I(X, S; Y \mid Q),
+
+with ``(S, Q)`` stationary :cite:`Sabag2017`. This is a convex program in the
+joint ``P(s, q, x)`` :cite:`Sabag2020`, solved with cvxpy when the optional
+``capacity`` extra is installed (``pip install sofic[capacity]``) and with
+scipy otherwise. On the lower side, any graph-based encoder
+:math:`P(x \mid s, q)` is a valid feedback input; its directed-information rate
+:math:`\bar h(Y) - \mathbb{E}_\pi H(Y \mid X, S)` is bracketed with block
+entropies :cite:`Cover2006` and optimized. When the encoder is also
+BCJR-invariant the single-letter rate :math:`I(X, S; Y \mid Q)` is achievable
+:cite:`Sabag2017,Sabag2020`. :func:`qgraph_bounds` tries Markov Q-graphs of the
+last ``k`` outputs plus any graphs you pass.
+
+Markov Q-graphs are not always enough. The trapdoor channel
+:cite:`Permuter2008` is closed by a three-node graph:
+
+.. ipython::
+
+   In [1]: import numpy as np
+
+   In [2]: from sofic.examples.processes import trapdoor_et
+
+   In [3]: from sofic.generators.capacity import QGraph, qgraph_bounds
+
+   In [4]: trapdoor = trapdoor_et()
+
+   In [5]: graph = QGraph(("0", "1"), np.array([[0, 1], [0, 2], [1, 2]]), ("a", "b", "c"))
+
+   In [6]: bounds = qgraph_bounds(trapdoor, qgraphs=[graph], max_order=1)
+
+   @doctest
+   In [7]: round(bounds.upper, 6), round(np.log2((1 + np.sqrt(5)) / 2), 6)
+   Out[7]: (0.694242, 0.694242)
+
+   @doctest
+   In [8]: bounds.gap < 1e-4
+   Out[8]: True
+
+For the Ising channel :cite:`Elishco2014` (``ising_channel_et``), the four-node
+graph with table ``[[1, 2], [0, 2], [0, 3], [0, 2]]`` closes the bracket at
+:math:`C_{\mathrm{fb}} \approx 0.575522`.
+
+Duality bound
+-------------
+
+:func:`duality_upper_bound` fixes a test distribution :math:`R(y \mid q)` on a
+Q-graph. For a unifilar channel the resulting bound is the optimal average
+reward of a *finite* Markov decision process with state ``(s, q)``, action
+``x`` and reward :math:`D(P_{Y \mid x, s} \| R_{Y \mid q})`
+:cite:`Huleihel2025`. Every test distribution gives a valid bound; the
+starting one is the output law of the Q-graph upper-bound optimizer, refined
+by Nelder--Mead.
+
+Value iteration
+---------------
+
+:func:`value_iteration_estimate` solves the belief-state dynamic program of
+:cite:`Permuter2008,Yang2005` (state :math:`P(S \mid y^t)`, action
+:math:`P(x \mid s)`, reward :math:`I(X, S; Y)`) on a gridded simplex. Its
+bracket bounds the gain of the discretized problem only, so it is returned
+with ``certified=False``; ``details["policy"]`` holds the greedy input law at
+every grid belief.
+
+:func:`feedback_capacity` dispatches on ``method`` (``"qgraph"``,
+``"duality"``, ``"value_iteration"``, or ``"all"``).
+
+Capacity without feedback
+=========================
+
+* :func:`markov_input_lower_bound` optimizes an order-``k`` Markov input and
+  reports a block-entropy lower bound on its information rate. That rate is
+  the objective of the generalized Blahut--Arimoto algorithm
+  :cite:`Vontobel2008,Arnold2006`, maximized here directly.
+* :func:`duality_upper_bound_no_feedback` evaluates the dual bound with a
+  Q-graph test distribution, a deterministic dynamic program over
+  :math:`P(S, Q \mid x^t)` :cite:`Huleihel2021`. Its value function is convex,
+  so interpolating on a simplex grid over-estimates it and the gridded bound
+  stays valid at any resolution.
+
+For the trapdoor channel, an order-3 Markov input certifies ``C >= 0.5727``
+and the order-1 Markov Q-graph at ``belief_resolution=24`` gives
+``C <= 0.5936``. The best published bounds are ``0.572`` and ``0.584``
+:cite:`Huleihel2021`. :func:`capacity` dispatches on ``method``
+(``"markov_input"``, ``"duality"``, or ``"all"``).
+
+API
+===
+
+.. autoclass:: CapacityBounds
+   :members: gap, value
+
+.. autoclass:: QGraph
+   :members: markov, from_dfa, aligned
+
+.. autofunction:: feedback_capacity
+.. autofunction:: qgraph_bounds
+.. autofunction:: qgraph_upper_bound
+.. autofunction:: qgraph_lower_bound
+.. autofunction:: encoder_rate_bounds
+.. autofunction:: duality_upper_bound
+.. autofunction:: value_iteration_estimate
+.. autofunction:: capacity
+.. autofunction:: markov_input_lower_bound
+.. autofunction:: duality_upper_bound_no_feedback
+.. autofunction:: channel_tensors

docs/generators/epsilon_transducer.rst

@@ -57,11 +57,16 @@ quantity, reusing the directional-flow estimators of :doc:`directional_flow`.
 * :func:`directed_information` / :func:`transfer_entropy` — input-to-output
   information flow of the driven joint process.
 
+These measures fix the input process. Optimizing over inputs gives the channel
+capacity; :meth:`~EpsilonTransducer.feedback_capacity` and
+:meth:`~EpsilonTransducer.capacity` bracket it with and without feedback
+(:doc:`channel_capacity`).
+
 API
 ===
 
 .. autoclass:: EpsilonTransducer
-   :members: from_channel, from_iohmm, from_joint_generator, from_paired_sequences, is_unifilar, causal_states, statistical_complexity, driven_entropy_rate, directed_information, transfer_entropy, from_wfst
+   :members: from_channel, from_iohmm, from_joint_generator, from_paired_sequences, is_unifilar, causal_states, statistical_complexity, driven_entropy_rate, directed_information, transfer_entropy, feedback_capacity, capacity, from_wfst
 
 .. currentmodule:: sofic.generators.channel_measures
 

docs/generators/generators.rst

@@ -29,6 +29,7 @@ Computational mechanics
 
    epsilon_machine
    epsilon_transducer
+   channel_capacity
    block_convergence
    bidirectional_epsilon_machine
    mixed_state_presentation

docs/references.bib

@@ -1311,3 +1311,113 @@ @incollection{blackwell1957entropy
   address   = {Prague},
   year      = {1957},
 }
+
+@article{Permuter2008,
+  author  = {Permuter, Haim and Cuff, Paul and Van Roy, Benjamin and Weissman, Tsachy},
+  title   = {Capacity of the Trapdoor Channel With Feedback},
+  journal = {IEEE Transactions on Information Theory},
+  volume  = {54},
+  number  = {7},
+  pages   = {3150--3165},
+  year    = {2008},
+  doi     = {10.1109/TIT.2008.924681},
+}
+
+@article{Elishco2014,
+  author  = {Elishco, Ohad and Permuter, Haim},
+  title   = {Capacity and Coding for the {Ising} Channel With Feedback},
+  journal = {IEEE Transactions on Information Theory},
+  volume  = {60},
+  number  = {9},
+  pages   = {5138--5149},
+  year    = {2014},
+  doi     = {10.1109/TIT.2014.2331951},
+}
+
+@article{Yang2005,
+  author  = {Yang, Shaohua and Kav{\v{c}}i{\'c}, Aleksandar and Tatikonda, Sekhar},
+  title   = {Feedback Capacity of Finite-State Machine Channels},
+  journal = {IEEE Transactions on Information Theory},
+  volume  = {51},
+  number  = {3},
+  pages   = {799--810},
+  year    = {2005},
+  doi     = {10.1109/TIT.2004.842626},
+}
+
+@article{Sabag2017,
+  author  = {Sabag, Oron and Permuter, Haim H. and Pfister, Henry D.},
+  title   = {A Single-Letter Upper Bound on the Feedback Capacity of Unifilar Finite-State Channels},
+  journal = {IEEE Transactions on Information Theory},
+  volume  = {63},
+  number  = {3},
+  pages   = {1392--1409},
+  year    = {2017},
+  doi     = {10.1109/TIT.2016.2636851},
+}
+
+@article{Sabag2020,
+  author  = {Sabag, Oron and Huleihel, Bashar and Permuter, Haim H.},
+  title   = {Graph-Based Encoders and Their Performance for Finite-State Channels With Feedback},
+  journal = {IEEE Transactions on Communications},
+  volume  = {68},
+  number  = {4},
+  pages   = {2106--2117},
+  year    = {2020},
+  doi     = {10.1109/TCOMM.2020.2965454},
+}
+
+@article{Huleihel2021,
+  author  = {Huleihel, Bashar and Sabag, Oron and Permuter, Haim H. and Kashyap, Navin and Shamai (Shitz), Shlomo},
+  title   = {Computable Upper Bounds on the Capacity of Finite-State Channels},
+  journal = {IEEE Transactions on Information Theory},
+  volume  = {67},
+  number  = {9},
+  pages   = {5674--5692},
+  year    = {2021},
+  doi     = {10.1109/TIT.2021.3091691},
+}
+
+@article{Huleihel2025,
+  author  = {Huleihel, Bashar and Sabag, Oron and Aharoni, Ziv and Permuter, Haim H.},
+  title   = {The Duality Upper Bound for Finite-State Channels With Feedback},
+  journal = {IEEE Transactions on Information Theory},
+  volume  = {71},
+  number  = {5},
+  pages   = {3255--3270},
+  year    = {2025},
+  doi     = {10.1109/TIT.2025.3545688},
+}
+
+@article{Vontobel2008,
+  author  = {Vontobel, Pascal O. and Kav{\v{c}}i{\'c}, Aleksandar and Arnold, Dieter M. and Loeliger, Hans-Andrea},
+  title   = {A Generalization of the {Blahut}--{Arimoto} Algorithm to Finite-State Channels},
+  journal = {IEEE Transactions on Information Theory},
+  volume  = {54},
+  number  = {5},
+  pages   = {1887--1918},
+  year    = {2008},
+  doi     = {10.1109/TIT.2008.920243},
+}
+
+@article{Arnold2006,
+  author  = {Arnold, Dieter M. and Loeliger, Hans-Andrea and Vontobel, Pascal O. and Kav{\v{c}}i{\'c}, Aleksandar and Zeng, Wei},
+  title   = {Simulation-Based Computation of Information Rates for Channels With Memory},
+  journal = {IEEE Transactions on Information Theory},
+  volume  = {52},
+  number  = {8},
+  pages   = {3498--3508},
+  year    = {2006},
+  doi     = {10.1109/TIT.2006.878110},
+}
+
+@article{Grigorescu2024,
+  author  = {Grigorescu, Andrea and Boche, Holger and Schaefer, Rafael F. and Poor, H. Vincent},
+  title   = {Capacity of Finite State Channels With Feedback: Algorithmic and Optimization Theoretic Properties},
+  journal = {IEEE Transactions on Information Theory},
+  volume  = {70},
+  number  = {8},
+  pages   = {5413--5426},
+  year    = {2024},
+  doi     = {10.1109/TIT.2024.3411919},
+}

pyproject.toml

@@ -59,7 +59,7 @@ docs = [
     "sofic[symbolic]",
 ]
 dev = [
-    "sofic[test,viz]",
+    "sofic[test,viz,capacity]",
     "sofic[docs]",
     "ruff",
     "ty",
@@ -75,6 +75,9 @@ bayes = [
 sat = [
     "python-sat>=1.8",
 ]
+capacity = [
+    "cvxpy>=1.4",
+]
 symbolic = [
     "sympy>=1.12",
     "dit[symbolic]>=2.3",