Add spectral reconstruction from sequences to ε-machines

README.rst

@@ -172,6 +172,7 @@ read off computational-mechanics quantities:
 
    eps = EpsilonMachine.from_hmm(gm)          # minimize an HMM presentation
    # eps = EpsilonMachine.from_sequence(data, method="cssr", Lmax=4)  # infer
+   # eps = EpsilonMachine.from_sequence(data, method="spectral", prefix_length=3, rank=2)
 
    eps.statistical_complexity()   # 0.9183 bits  (C_mu)
    eps.entropy_rate()             # 0.6667 bits/symbol  (h_mu)

docs/generators/epsilon_inference.rst

@@ -9,10 +9,11 @@ Sample-based reconstruction of ε-machines from observed symbol sequences.
 This complements the **oracle** path :func:`~sofic.generators.epsilon_construction.build_epsilon_machine`,
 which merges probabilistically equivalent states in a *given* generator.
 
-Two algorithms are implemented:
+Three algorithms are implemented:
 
 * **CSSR** (Causal-State Splitting Reconstruction) :cite:`Shalizi2002,Shalizi2004`
 * **Subtree merging** (Crutchfield--Young batch reconstruction) :cite:`CrutchfieldYoung1989,Crutchfield1994`
+* **Spectral** (Hankel SVD, then mixed-state extraction) :cite:`Balle2014,Hsu2012,Ellison2009`
 
 Quick start
 ===========
@@ -55,16 +56,50 @@ numerical tolerance for finite-sample estimates.
 
 .. autofunction:: subtree_merge
 
+Spectral reconstruction
+=======================
+
+Spectral reconstruction learns a weighted finite automaton from the Hankel matrix
+of block probabilities :cite:`Balle2014,Hsu2012`, then extracts causal states as
+mixed states of the learned operators :cite:`Ellison2009`. When the operators
+are non-negative in the learned basis this is a Mealy projection followed by
+:meth:`~sofic.generators.epsilon_machine.EpsilonMachine.from_hmm`; signed
+operators use mixed-state enumeration of the observable operators (not a
+clustering heuristic). The same extraction is
+:func:`~sofic.inference.spectral.project_to_epsilon_machine`.
+
+Pass ``rank`` when the model order is known (two for the golden mean and even
+process). Otherwise the Hankel singular-value gap selects the rank.
+
+.. code-block:: python
+
+   inferred = EpsilonMachine.from_sequence(
+       observations, method="spectral", prefix_length=3, rank=2
+   )
+
+.. autofunction:: spectral
+
 Unified entry point
 ===================
 
 Use :meth:`~sofic.generators.epsilon_machine.EpsilonMachine.from_sequence` to
-dispatch to CSSR or subtree merging (see :doc:`epsilon_machine`).
+dispatch to CSSR, subtree merging, or spectral reconstruction
+(see :doc:`epsilon_machine`).
+
+Related inference methods
+=========================
+
+**transCSSR** — input/output ε-transducers; see :doc:`epsilon_transducer_inference`.
 
-Related inference methods (not yet implemented)
-===============================================
+**Bayesian structural inference** — conjugate Dirichlet–multinomial evidence over
+candidate unifilar topologies :cite:`Strelioff2014`; see :doc:`../inference/epsilon`.
+This is not a Gibbs clustering heuristic over history labels.
 
-**transCSSR** — input/output ε-transducers (Darmon, 2014).
+Subtree merging is the literature reconstruction by morph clustering; a separate
+agglomerative "causal-state merging" procedure is not provided. k-means on
+history-morph embeddings (sometimes labelled "neural state discovery" despite
+involving no neural network) is also not provided — mixed-state extraction is
+the causal-state construction used after spectral learning.
 
 **VLMC / context algorithm** — sparse Markov trees; not causally minimal in general.
 
@@ -76,4 +111,5 @@ Related inference methods (not yet implemented)
 
 **RKHS ε-machines** — continuous-time extension (arXiv:2011.14821).
 
-See also :doc:`epsilon_machine`, :doc:`constructions`, and :doc:`hmm_inference`.
+See also :doc:`epsilon_machine`, :doc:`constructions`, :doc:`hmm_inference`,
+and :doc:`../inference/spectral`.

docs/inference/inference.rst

@@ -31,7 +31,8 @@ The historical names ``InferMC`` and ``InferEM`` are retained as aliases for
 .. note::
 
    For *non-Bayesian* reconstruction — Causal-State Splitting Reconstruction
-   (CSSR) and subtree merging — see the point-estimate routines in
+   (CSSR), subtree merging, and spectral mixed-state extraction — see the
+   point-estimate routines in
    :doc:`../generators/epsilon_inference`,
    :doc:`../generators/hmm_inference`, and
    :doc:`../generators/stack_inference`.

docs/inference/spectral.rst

@@ -61,6 +61,25 @@ in the learned basis and raises :class:`SpectralInferenceError` otherwise.
 
    machine = project_to_nmachine(model)   # observable-operator generator with a graph
 
+Projection to an ε-machine
+==========================
+
+:func:`project_to_epsilon_machine` extracts the causal presentation: a
+non-negative Mealy projection when one exists in the learned basis, otherwise
+mixed-state enumeration of the observable operators
+:cite:`Ellison2009`. The same path is
+:func:`~sofic.generators.epsilon_inference.spectral` /
+``EpsilonMachine.from_sequence(..., method="spectral")``.
+
+.. code-block:: python
+
+   from sofic.generators.epsilon_inference import spectral
+   from sofic.examples import golden_mean
+
+   process = golden_mean(0.5)
+   eps = spectral(word_probability=process.word_probability, alphabet=(0, 1), prefix_length=3, rank=2)
+   len(list(eps.states()))  # 2
+
 API
 ===
 
@@ -74,4 +93,6 @@ API
 
 .. autofunction:: project_to_mealy
 
+.. autofunction:: project_to_epsilon_machine
+
 .. autoexception:: SpectralInferenceError

sofic/generators/__init__.py

@@ -16,7 +16,7 @@
     synergistic_information_flow,
     transfer_entropy,
 )
-from sofic.generators.epsilon_inference import cssr, subtree_merge
+from sofic.generators.epsilon_inference import cssr, spectral, subtree_merge
 from sofic.generators.epsilon_machine import EpsilonMachine
 from sofic.generators.epsilon_transducer import EpsilonTransducer
 from sofic.generators.lumping import LumpabilityError, is_lumpable, lump, normalize_partition
@@ -71,6 +71,7 @@
     "is_lumpable",
     "lump",
     "normalize_partition",
+    "spectral",
     "subtree_merge",
     "fit_stack_hmm_mle",
     "learn_stack_hmm_papni",

sofic/generators/epsilon_inference.py

@@ -1,13 +1,15 @@
-"""Sample-based ε-machine reconstruction (CSSR and subtree merging).
+"""Sample-based ε-machine reconstruction (CSSR, subtree merging, and spectral).
 
 CSSR follows Shalizi, Shalizi & Crutchfield (arXiv:cs/0210025). Subtree merging
-follows Crutchfield & Young (PRL 1989; PRE 1994).
+follows Crutchfield & Young (PRL 1989; PRE 1994). Spectral reconstruction learns
+a weighted finite automaton by Hankel SVD :cite:`Balle2014,Hsu2012` and extracts
+causal states as mixed states of the learned operators :cite:`Ellison2009`.
 """
 
 from __future__ import annotations
 
 from collections import Counter, defaultdict
-from collections.abc import Callable, Sequence
+from collections.abc import Callable, Iterable, Sequence
 from dataclasses import dataclass, field
 from typing import Any, ClassVar, Literal
 
@@ -641,3 +643,61 @@ def subtree_merge(
         history: state_id for state_id, histories_in_state in states.items() for history in histories_in_state
     }
     return _counts_to_mealy(states, counts, history_to_state, seq, length=L)
+
+
+def spectral(
+    sequences: Iterable[Any] | None = None,
+    *,
+    word_probability: Callable[[Sequence[Any]], float] | None = None,
+    alphabet: Sequence[Any] | None = None,
+    rank: int | None = None,
+    prefix_length: int = 3,
+    suffix_length: int | None = None,
+    singular_value_threshold: float = 1e-3,
+    min_singular_value: float = 1e-12,
+    max_states: int = 10_000,
+) -> EpsilonMachine:
+    """Reconstruct an ε-machine by spectral learning then mixed-state extraction.
+
+    Learns a weighted finite automaton / observable-operator model from block
+    statistics :cite:`Balle2014,Hsu2012`, then extracts causal states as the
+    mixed states of those operators :cite:`Ellison2009`. When the learned
+    operators are non-negative this is a Mealy projection followed by
+    :meth:`~sofic.generators.epsilon_machine.EpsilonMachine.from_hmm`; signed
+    operators use mixed-state enumeration rather than a clustering heuristic.
+
+    Parameters
+    ----------
+    sequences
+        A single observed realization or an iterable of realizations. Ignored
+        when ``word_probability`` is given.
+    word_probability
+        Optional exact block-probability function ``f(word) -> float``.
+        ``alphabet`` is then required.
+    alphabet
+        Observation alphabet. Inferred from ``sequences`` when omitted.
+    rank
+        Number of latent states. When ``None`` the rank is chosen from the
+        Hankel singular-value spectrum.
+    prefix_length, suffix_length
+        Maximum lengths of the prefix and suffix bases. ``suffix_length``
+        defaults to ``prefix_length``.
+    singular_value_threshold, min_singular_value
+        Cutoffs for automatic rank selection; see
+        :func:`~sofic.inference.spectral.learn_spectral_wfa`.
+    max_states
+        Safety cap on enumerated mixed states.
+    """
+    from sofic.inference.spectral import learn_spectral_wfa, project_to_epsilon_machine
+
+    model = learn_spectral_wfa(
+        sequences,
+        word_probability=word_probability,
+        alphabet=alphabet,
+        rank=rank,
+        prefix_length=prefix_length,
+        suffix_length=suffix_length,
+        singular_value_threshold=singular_value_threshold,
+        min_singular_value=min_singular_value,
+    )
+    return project_to_epsilon_machine(model, max_states=max_states)

sofic/generators/epsilon_machine.py

@@ -51,7 +51,7 @@ def from_sequence(
         cls,
         sequence: Sequence[Any],
         *,
-        method: Literal["cssr", "subtree"] = "cssr",
+        method: Literal["cssr", "subtree", "spectral"] = "cssr",
         **kwargs: Any,
     ) -> EpsilonMachine:
         """Reconstruct an ε-machine from an observed symbol sequence.
@@ -61,11 +61,14 @@ def from_sequence(
         sequence
             Observed process realization.
         method
-            ``"cssr"`` for Causal-State Splitting Reconstruction, or
-            ``"subtree"`` for depth-``L`` subtree merging (pass ``L=...``).
+            ``"cssr"`` for Causal-State Splitting Reconstruction,
+            ``"subtree"`` for depth-``L`` subtree merging (pass ``L=...``), or
+            ``"spectral"`` for Hankel-SVD learning followed by mixed-state
+            extraction.
         **kwargs
-            Forwarded to :func:`~sofic.generators.epsilon_inference.cssr` or
-            :func:`~sofic.generators.epsilon_inference.subtree_merge`.
+            Forwarded to :func:`~sofic.generators.epsilon_inference.cssr`,
+            :func:`~sofic.generators.epsilon_inference.subtree_merge`, or
+            :func:`~sofic.generators.epsilon_inference.spectral`.
         """
         if method == "cssr":
             from sofic.generators.epsilon_inference import cssr
@@ -75,6 +78,10 @@ def from_sequence(
             from sofic.generators.epsilon_inference import subtree_merge
 
             return subtree_merge(sequence, **kwargs)
+        if method == "spectral":
+            from sofic.generators.epsilon_inference import spectral
+
+            return spectral(sequence, **kwargs)
         raise ValueError(f"unknown inference method {method!r}")
 
     def copy(self) -> Self:

sofic/inference/__init__.py

@@ -17,6 +17,7 @@
     SpectralInferenceError,
     hankel_matrices,
     learn_spectral_wfa,
+    project_to_epsilon_machine,
     project_to_mealy,
     project_to_nmachine,
     spectral_singular_values,
@@ -37,6 +38,7 @@
     "SpectralInferenceError",
     "hankel_matrices",
     "learn_spectral_wfa",
+    "project_to_epsilon_machine",
     "project_to_mealy",
     "project_to_nmachine",
     "spectral_singular_values",