Generator verification: equations are necessary only; exact word-process check

CHANGELOG.md

@@ -5,7 +5,20 @@
 A second correctness review: every fix below was reproduced against a
 brute-force reference or a closed-form value and has a regression test.
 
-### Bidirectional builds
+### Generator verification
+
+- The one-step generator equations (edge law, forward stationarity, summed future-side consistency) are
+  necessary but not sufficient: a kernel with kernel-LP violation 0 can emit the wrong word process,
+  because nothing forces a state's belief `R_w` to be the generated future's actual conditional. On survey
+  process `((1,·),(2,0),(1,0))` an 11-state example matches words up to length 3 and is off by 0.024 at
+  length 8. `generator_kernel`, the result flags (`n_g_optima_generating`, `separating_variables_generate`,
+  `cg_exceeds_g`) and the docs now say so; positive violations still rule a variable out.
+- Candidate verification (`_reproduction_error`, used by every search that returns a generator) adds a
+  difference-automaton test (`_equivalence_error`): breadth-first search over words keeping those whose
+  difference vectors are independent (at most `n1 + n2`, Schützenberger / Tzeng), with the largest
+  word-probability deviation over the visited words. The fixed cutoff (at most length 8) alone could
+  accept a candidate that matched short words and drifted later; two HMMs are only guaranteed equal once
+  they agree up to length `n1 + n2 - 1`.
 
 - `stationary_distribution_hmm` solves chains with more than 2,000 states sparsely (sparse transition
   matrix from the edges, one balance equation replaced by the normalization, `spsolve`). The dense

docs/generators/generative_models.rst

@@ -168,12 +168,20 @@ Two narrower tools share the same equations:
   interleaved process's minimal generator is the memory-1 lift of a single
   flag, and :func:`~sofic.generators.coupled_generators.memory_lift_search`
   finds it.
-* :func:`~sofic.generators.coupled_generators.generator_kernel` decides whether
-  a given channel's variable is the state of *any* generator: for a fixed channel
-  the equations are linear in :math:`K`, and a linear program returns the
-  smallest violation and the kernel. The exact-common-information optimum of the
+* :func:`~sofic.generators.coupled_generators.generator_kernel` tests a given
+  channel's variable against the generator equations: for a fixed channel the
+  equations are linear in :math:`K`, and a linear program returns the smallest
+  violation and the kernel. The equations are necessary, so a positive violation
+  rules the variable out: the exact-common-information optimum of the
   :math:`(2, \infty)` process has violation 0.0625, so no generator has it as
-  its state. :func:`~sofic.generators.coupled_generators.bilevel_search`
+  its state. They are not sufficient: future-side consistency only constrains
+  the belief :math:`R_w` through its one-step average, and a kernel with
+  violation 0 can still emit the wrong word distribution (on a survey process an
+  11-state example matches words up to length 3 and is off by about 0.02 at
+  length 8). Every candidate the pipeline reports is therefore checked against
+  the source's word process: word probabilities up to a fixed length, and a
+  difference-automaton test (Schützenberger, Tzeng) over words that span all word
+  vectors, so a candidate that matches short words but drifts later is rejected. :func:`~sofic.generators.coupled_generators.bilevel_search`
   (experimental) searches channels against this violation.
 
 All three searches are local and return upper bounds on :math:`C_g` for the

sofic/generators/coupled_generators.py

@@ -464,9 +464,10 @@ def value(z: np.ndarray, mu: float) -> float:
 def generator_kernel(bidir: BidirectionalEpsilonMachine | _GeneratorProblem, c: np.ndarray) -> tuple[float, np.ndarray]:
     """Best kernel for a fixed channel ``c[pair, w]``: the generator equations are linear in ``K``.
 
-    Solves the linear program minimizing the largest violation ``t`` over stochastic kernels; ``t = 0`` means
-    the channel's variable is the state of a generator, whatever the coupling. ``c`` is indexed by the
-    causal-state pairs in sorted order. Returns ``(t, K)``.
+    Solves the linear program minimizing the largest violation ``t`` over stochastic kernels. The equations
+    are necessary, so ``t > 0`` rules the channel's variable out as a generator state, whatever the coupling;
+    they are not sufficient, so ``t = 0`` still needs a word-level check (the HMM built from ``K`` can emit
+    the wrong word distribution). ``c`` is indexed by the causal-state pairs in sorted order. Returns ``(t, K)``.
     """
     problem = _as_problem(bidir)
     pi, Q, R, _ = problem.channel_seed(c)

sofic/generators/generative_search.py

@@ -99,6 +99,14 @@ class GenerativeModelResult:
     (state entropy over generators whose state is a function of :math:`(S^+, S^-)`).
     ``budget_exhausted`` is set when a ``time_budget`` ran out and later stages were skipped; the bracket
     is still valid, just possibly wider. Each ``stages`` entry records ``elapsed`` seconds.
+
+    ``n_g_optima_generating``, ``separating_variables_generate`` and ``cg_exceeds_g`` report the kernel LP
+    (:func:`~sofic.generators.coupled_generators.generator_kernel`), whose one-step generator equations
+    are necessary but not sufficient. A failed kernel LP rules an exact-common-information optimum out as a
+    generator state, so ``cg_exceeds_g=True`` (no optimum found passes) is evidence that :math:`C_g` exceeds
+    :math:`\\op{G}{S^+ : S^-}`, though only over the optima found. A passing kernel LP is only a candidate:
+    ``cg_exceeds_g=False`` does not prove :math:`C_g = \\op{G}{S^+ : S^-}`. ``model`` and ``upper`` always come
+    from a candidate whose word process was checked against the source (``_reproduction_error``).
     """
 
     model: CoupledGenerativeModel
@@ -128,7 +136,7 @@ def gap(self) -> float:
 
     @property
     def g_optima_generate(self) -> bool:
-        """Whether some :math:`\\op{G}{S^+ : S^-}` optimum passed the generator kernel LP."""
+        """Whether some :math:`\\op{G}{S^+ : S^-}` optimum passed the generator kernel LP (a necessary condition)."""
         return self.n_g_optima_generating > 0
 
 

sofic/generators/minimal_generative_model.py

@@ -888,13 +888,67 @@ def _reproduction_length(reference: Any) -> int:
 
 
 def _reproduction_error(model: Any, reference: Any, *, max_length: int) -> float:
-    """Max abs word-probability deviation of ``model`` from ``reference`` up to ``max_length``."""
+    """Deviation of ``model``'s word process from ``reference``'s, over all word lengths.
+
+    The larger of the max abs word-probability deviation up to ``max_length`` and
+    :func:`_equivalence_error`, which bounds the deviation of every word. Agreement up to a fixed
+    length is not enough: a model can match all short words and still differ later (models that satisfy
+    the one-step generator equations do this), while two HMMs with ``n1`` and ``n2`` states are only
+    guaranteed equal once they agree up to length ``n1 + n2 - 1``.
+    """
     error = 0.0
     for length in range(max_length + 1):
         model_words = model.word_probabilities(length, sparse=False)
         reference_words = reference.word_probabilities(length, sparse=False)
         for word in set(model_words) | set(reference_words):
             error = max(error, abs(model_words.get(word, 0.0) - reference_words.get(word, 0.0)))
+    return max(error, _equivalence_error(model, reference))
+
+
+def _equivalence_error(model: Any, reference: Any, *, rank_tol: float = 1e-7) -> float:
+    """Word-probability deviation of two HMMs on a spanning set of words (difference-automaton test).
+
+    The difference automaton has start vector ``[s1, -s2]`` and block-diagonal symbol matrices, so a word
+    ``w``'s probability difference is ``v_w . 1`` with ``v_w = [s1, -s2] M_w``. Breadth-first search
+    expands a word only when its ``v_w`` is (relatively, ``rank_tol``) outside the span of the words kept
+    so far; at most ``n1 + n2`` are kept, and every word's ``v_w`` lies in their span (Schutzenberger,
+    Tzeng). In exact arithmetic the processes are therefore equal iff every visited word has deviation 0.
+    Returns the largest ``|v_w . 1|`` over the visited words, which include words far longer than any fixed
+    cutoff along the directions where the two processes could differ.
+    """
+    from sofic.generators.matrices import start_vector, symbol_matrices
+
+    m1, m2 = model.to_mealy(), reference.to_mealy()
+    A, B = symbol_matrices(m1), symbol_matrices(m2)
+    s1 = np.asarray(start_vector(m1, None), dtype=float)
+    s2 = np.asarray(start_vector(m2, None), dtype=float)
+    n1, n2 = len(s1), len(s2)
+    blocks = {}
+    for x in set(A) | set(B):
+        M = np.zeros((n1 + n2, n1 + n2))
+        if x in A:
+            M[:n1, :n1] = np.asarray(A[x], dtype=float)
+        if x in B:
+            M[n1:, n1:] = np.asarray(B[x], dtype=float)
+        blocks[x] = M
+    ones = np.ones(n1 + n2)
+    basis: list[np.ndarray] = []
+    queue = [np.concatenate([s1, -s2])]
+    error = 0.0
+    while queue and len(basis) < n1 + n2:
+        v = queue.pop(0)
+        error = max(error, abs(float(v @ ones)))
+        residual = v.copy()
+        for _ in range(2):  # re-orthogonalize: one Gram-Schmidt pass loses orthogonality on near-dependent vectors
+            for b in basis:
+                residual -= (b @ residual) * b
+        norm = float(np.linalg.norm(residual))
+        if norm <= rank_tol * max(1.0, float(np.linalg.norm(v))):
+            continue
+        basis.append(residual / norm)
+        queue.extend(v @ M for M in blocks.values())
+    for v in queue:  # one-step extensions of the last kept words
+        error = max(error, abs(float(v @ ones)))
     return error
 
 

tests/test_reproduction_check.py

@@ -0,0 +1,45 @@
+"""Word-process verification of candidate generators beyond a fixed word length."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from sofic.examples import golden_mean
+from sofic.examples.epsilon_machines import from_symbol_matrices
+from sofic.generators.mealy import MealyHMM
+from sofic.generators.minimal_generative_model import _equivalence_error, _reproduction_error
+
+
+def _fair_coin():
+    return from_symbol_matrices((0,), (0, 1), {0: np.array([[0.5]]), 1: np.array([[0.5]])})
+
+
+def _alternating():
+    return from_symbol_matrices(
+        (0, 1), (0, 1), {0: np.array([[0.0, 1.0], [0.0, 0.0]]), 1: np.array([[0.0, 0.0], [1.0, 0.0]])}
+    )
+
+
+def test_processes_equal_on_single_symbols_are_told_apart() -> None:
+    # both emit each symbol with probability 1/2; they differ first on words of length 2
+    coin, alternating = _fair_coin(), _alternating()
+    assert _equivalence_error(coin, alternating) >= 0.25 - 1e-12
+    assert _reproduction_error(coin, alternating, max_length=1) >= 0.25 - 1e-12
+
+
+def test_presentations_of_one_process_agree() -> None:
+    g = golden_mean(0.5)
+    assert _equivalence_error(g, g) <= 1e-12
+    # a non-unifilar presentation of the golden mean: state A split into two copies
+    split = MealyHMM(
+        observation_alphabet=frozenset({"0", "1"}), initial_distribution={"A1": 1 / 3, "A2": 1 / 3, "B": 1 / 3}
+    )
+    for state in ("A1", "A2", "B"):
+        split.graph.add_state(state)
+    for source in ("A1", "A2"):
+        split.add_transition(source, "A1", "0", 0.25)
+        split.add_transition(source, "A2", "0", 0.25)
+        split.add_transition(source, "B", "1", 0.5)
+    split.add_transition("B", "A1", "0", 0.5)
+    split.add_transition("B", "A2", "0", 0.5)
+    assert _equivalence_error(split, g) <= 1e-12