feat(#162): use exponential action for relaxation trajectories

CHANGELOG.md

@@ -35,6 +35,24 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   eigensolver paths use, so one quantity is not measured two different ways.
 
 ### Changed
+- **Relaxation propagation now computes the action `exp(tL) rho0` with
+  `scipy.sparse.linalg.expm_multiply` instead of materialising the full dense
+  matrix exponential (issue #162).** Exact NumPy-`linspace` grids use SciPy's
+  interval action API; arbitrary grids use independent actions from the same
+  initial state, avoiding stepwise error accumulation. The low-level propagator
+  accepts dense and SciPy sparse matrices and preserves the exact `t=0` state.
+  The #156 fail-closed contract remains in force: non-finite `L*t`, an
+  unrepresentable 1-norm, action failures, or non-finite output states raise
+  `UnrepresentableTrajectoryError` instead of flowing into fitting. The
+  relaxation result now records the backend plus trajectory-wide trace error,
+  Hermiticity defect and minimum Hermitian eigenvalue; positivity drift is
+  measured, never clipped. Dense-reference, sparse-vs-dense, non-uniform-grid,
+  and `L -> cL, t -> t/c` metamorphic tests bound the change. This is a
+  NUMERICAL METHODOLOGY change: floating-point trajectories may differ at
+  round-off level from the old full-`expm` route, although the mathematical
+  evolution is unchanged. The top-level orchestrator remains dense-only; sparse
+  support here is scoped to the propagation primitive. No manifest/schema
+  field changes.
 - **The Gaussian likelihood behind AICc is evaluated in log-RSS space, and an
   exact-zero RSS is now an explicit abstention (issue #135).** This is a
   METHODOLOGY change with user-visible consequences: it can reorder AICc,

CITATION.cff

@@ -153,6 +153,18 @@ abstract: >-
 # than reduction_rtol * (eps * ||L||_F + d**3 * 2**-1074), reduction_rtol
 # defaulting to the existing ZERO_MODE_EPS_FACTOR. Like the restriction itself
 # this moves no reported number: no diagnostic consumes it yet.
+# Also pending for the next cut: relaxation trajectories use the
+# Al-Mohy--Higham exponential-action route exposed by SciPy
+# `sparse.linalg.expm_multiply` instead of materialising `expm(tL)` and then
+# multiplying by the initial state (issue #162). Exact linspace grids use the
+# interval API; non-uniform grids use independent actions from the same initial
+# state. The mathematical semigroup is unchanged, but floating-point values may
+# differ at round-off level, so this is a NUMERICAL METHODOLOGY change rather
+# than a pure refactor. Backend identity plus trace/Hermiticity/positivity-drift
+# measurements are recorded on RelaxationResult. The low-level action path is
+# dense/sparse; this must NOT be described as making diagnose() sparse-native.
+# The #156 representability guard remains a prerequisite and fails closed before
+# non-finite or unrepresentable action inputs can enter downstream fits.
 keywords:
   - open quantum systems
   - Lindblad

README.md

@@ -237,6 +237,17 @@ LiouScope is built for paper-grade reproducibility:
 - **Anchor tests.** `tests/test_anchors.py` locks the numerical anchors that paper figures depend on;
   changes to physics code that move these values are caught in CI.
 - **Paper-figure pipeline.** `figures/generate_all.py` regenerates Fig 1-3 deterministically.
+- **Relaxation trajectories use exponential action, not a materialised propagator.**
+  `diagnostics.relaxation` evolves `exp(tL) rho0` with
+  `scipy.sparse.linalg.expm_multiply`: exact `linspace` grids use its interval
+  API, while arbitrary/non-uniform grids evaluate each requested time directly
+  from the same initial state. `RelaxationResult.trajectory_backend` records
+  which path was used, together with maximum trace error, maximum Hermiticity
+  defect, and the minimum eigenvalue of the Hermitian part seen along the
+  trajectory. These are audit measurements, not silent positivity repairs.
+  The low-level propagation primitive accepts dense and SciPy sparse matrices;
+  the top-level `diagnose()` orchestrator remains dense-only until its separate
+  sparse solver path is wired. Unrepresentable propagation fails closed.
 
 ---
 

benchmarks/issue162_expm_action.py

@@ -0,0 +1,81 @@
+"""Issue #162 benchmark: full dense expm vs exponential action.
+
+Run manually:
+    python benchmarks/issue162_expm_action.py
+
+The script is evidence generation, not a performance gate. It prints elapsed
+time, peak RSS (where ``resource`` is available), and dense/action agreement.
+Each backend runs in a fresh child process so peak RSS is not contaminated by
+the backend that ran first.
+"""
+
+from __future__ import annotations
+
+import multiprocessing as mp
+import time
+
+import numpy as np
+import scipy.linalg as sla
+import scipy.sparse.linalg as spla
+
+
+def _matrix(n: int, seed: int = 162) -> tuple[np.ndarray, np.ndarray]:
+    rng = np.random.default_rng(seed)
+    A = rng.standard_normal((n, n)) + 1j * rng.standard_normal((n, n))
+    # Shift left so exp(A) is not dominated by explosive growth.
+    A -= (float(n) + 1.0) * np.eye(n)
+    b = rng.standard_normal(n) + 1j * rng.standard_normal(n)
+    return np.asarray(A, dtype=complex), np.asarray(b, dtype=complex)
+
+
+def _peak_rss_mb() -> float:
+    try:
+        import resource
+
+        rss = float(resource.getrusage(resource.RUSAGE_SELF).ru_maxrss)
+        # Linux reports KiB; macOS reports bytes. This repo CI/benchmark target
+        # is Linux, but keep the output approximately useful on macOS.
+        return rss / (1024.0 if rss < 1.0e8 else 1024.0**2)
+    except ImportError:
+        return float("nan")
+
+
+def _worker(kind: str, n: int, q: mp.Queue) -> None:
+    A, b = _matrix(n)
+    started = time.perf_counter()
+    if kind == "dense":
+        out = sla.expm(A) @ b
+    elif kind == "action":
+        out = spla.expm_multiply(A, b)
+    else:
+        raise ValueError(kind)
+    q.put((kind, time.perf_counter() - started, _peak_rss_mb(), out))
+
+
+def _run(kind: str, n: int) -> tuple[float, float, np.ndarray]:
+    q: mp.Queue = mp.Queue()
+    proc = mp.Process(target=_worker, args=(kind, n, q))
+    proc.start()
+    proc.join()
+    if proc.exitcode != 0:
+        raise RuntimeError(f"{kind} child failed with exit code {proc.exitcode}")
+    _kind, elapsed, rss, out = q.get()
+    return float(elapsed), float(rss), np.asarray(out)
+
+
+def main() -> None:
+    for n in (64, 128, 256):
+        dense_t, dense_rss, dense = _run("dense", n)
+        action_t, action_rss, action = _run("action", n)
+        rel = float(
+            np.linalg.norm(action - dense)
+            / max(np.linalg.norm(dense), np.finfo(float).tiny)
+        )
+        print(
+            f"n={n:4d}  dense={dense_t:9.4f}s/{dense_rss:9.1f}MB  "
+            f"action={action_t:9.4f}s/{action_rss:9.1f}MB  relerr={rel:.3e}"
+        )
+
+
+if __name__ == "__main__":
+    main()

src/liouscope/_types.py

@@ -223,6 +223,13 @@ class RelaxationResult:
     # uses. Additive + defaulted so older callers / serialised reports stay valid.
     beta_D_linear: float = float("nan")
     linear_fit_model: str = "none"
+    # Issue #162: trajectory backend and numerical-state audit. Additive +
+    # defaulted so older synthetic callers and serialised report objects remain
+    # valid. These fields are measurements, not physicality gates.
+    trajectory_backend: str = "legacy_dense_expm"
+    trajectory_max_trace_error: float = float("nan")
+    trajectory_max_hermiticity_defect: float = float("nan")
+    trajectory_min_eigenvalue: float = float("nan")
 
 
 @dataclass(frozen=True, slots=True, kw_only=True)