Continuous-time Pauli propagation: real-space adaptive Lindbladian evolution

Cargo.toml

@@ -14,6 +14,7 @@ members = [
     "crates/ppvm-pauli-word",
     "crates/ppvm-pauli-sum",
     "crates/ppvm-sym",
+    "crates/ppvm-lindblad",
     "crates/ppvm-python-native",
     "crates/ppvm-tableau",
     "crates/ppvm-stim", "crates/stim-parser",

crates/ppvm-lindblad/Cargo.toml

@@ -0,0 +1,24 @@
+[package]
+name = "ppvm-lindblad"
+version = "0.1.0"
+edition = "2024"
+description = "Direct Heisenberg-picture Lindbladian evolution on an adaptive Pauli-string basis."
+
+[dependencies]
+fxhash = "0.2.1"
+ndarray = "0.17"
+num = "0.4.3"
+ppvm-traits     = { version = "0.1.0", path = "../ppvm-traits" }
+ppvm-pauli-word = { version = "0.1.0", path = "../ppvm-pauli-word" }
+ppvm-pauli-sum  = { version = "0.1.0", path = "../ppvm-pauli-sum" }
+rayon = "1.11"
+# Matrix-exponential action (Al-Mohy & Higham). QuSpin-rust is MIT-licensed;
+# the pinned rev is the commit that added the LICENSE file.
+quspin-expm = { git = "https://github.com/QuSpin/QuSpin-rust", rev = "a0ad6c9fe2e8063208f9ba1c6677150c993bb554" }
+# `QuSpinError` (the error type returned by the `LinearOperator` trait methods
+# we implement in `mf_expm.rs`) is not re-exported from `quspin-expm`'s root,
+# so we depend on `quspin-types` directly. Same git rev as `quspin-expm`.
+quspin-types = { git = "https://github.com/QuSpin/QuSpin-rust", rev = "a0ad6c9fe2e8063208f9ba1c6677150c993bb554" }
+
+[dev-dependencies]
+approx = "0.5.1"

crates/ppvm-lindblad/src/config.rs

@@ -0,0 +1,61 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+//! Configuration objects for the predictor-corrector stepper.
+
+/// Truncation and execution policy for a single predictor-corrector step
+/// ([`crate::LindbladSpec::pc_step`], [`crate::LindbladSpec::pc_step_timed`],
+/// [`crate::orbit_rep::pc_step_orbit_rep`]).
+///
+/// These are the per-run *tuning knobs*, kept separate from the per-call data
+/// (`basis`, `coeffs`, `dt`, `protected`, and — on the orbit path — the
+/// translation group and momentum).
+///
+/// `max_basis` is the primary accuracy/cost dial; `admit_basis` selects the
+/// displacement scheme; `drop_tol` is the churn valve of the admission-bound
+/// scheme; `tau_add` is a wall optimization at most.
+#[derive(Debug, Clone, Copy)]
+pub struct PcStepConfig {
+    /// Hard rank cap on the retained basis: after the corrector, only the
+    /// top-`max_basis` strings by `|coeff|` are kept (protected words always
+    /// survive). The primary convergence dial — verify by re-running at 2×.
+    pub max_basis: usize,
+    /// Working-set (admission) bound. When `Some(a)` with `a > max_basis`,
+    /// enrichment may grow the live basis to `a` and the final cap performs a
+    /// genuine top-`max_basis`-of-union rank displacement (the analog of
+    /// two-site TDVP truncation at `χ_max`); `drop_tol` is then not needed
+    /// for membership turnover. `None` bounds admission by `max_basis`
+    /// itself — the valve scheme, which requires `drop_tol > 0` to keep the
+    /// basis adapting once it fills.
+    pub admit_basis: Option<usize>,
+    /// Magnitude prune applied after the corrector: basis entries whose
+    /// `|coeff|` is below `drop_tol` are discarded (protected words are
+    /// always kept). `<= 0.0` disables pruning — valid only with
+    /// `admit_basis` set, otherwise the basis freezes once it fills the cap.
+    pub drop_tol: f64,
+    /// Optional absolute rate threshold on leakage admission: a candidate is
+    /// admitted only if its inflow rate exceeds `tau_add`. This is the
+    /// natural (dt- and drop_tol-independent) parameterization — the
+    /// admission accuracy cliff sits at a fixed `tau_add`. `None` = no
+    /// filter, the recommended default with cap-based truncation.
+    pub tau_add: Option<f64>,
+    /// When `Some(n)`, run the entire step inside a freshly built rayon
+    /// thread pool of `n` threads (useful for benchmarking parallel
+    /// scaling). When `None`, the global rayon pool is used.
+    pub num_threads: Option<usize>,
+}
+
+impl Default for PcStepConfig {
+    /// Uncapped, unfiltered, no pruning: the near-exact reference
+    /// configuration. Production runs should set `max_basis` (and usually
+    /// `admit_basis ≈ 2-3×` it).
+    fn default() -> Self {
+        Self {
+            max_basis: usize::MAX,
+            admit_basis: None,
+            drop_tol: 0.0,
+            tau_add: None,
+            num_threads: None,
+        }
+    }
+}

crates/ppvm-lindblad/src/error.rs

@@ -0,0 +1,76 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+//! Error type for [`crate::LindbladSpec`] construction and stepping.
+
+use crate::MAX_QUBITS;
+use std::fmt;
+
+/// Errors raised when constructing a [`crate::LindbladSpec`].
+#[derive(Debug, Clone)]
+pub enum Error {
+    TooManyQubits {
+        got: usize,
+    },
+    LengthMismatch {
+        what: &'static str,
+        a: usize,
+        b: usize,
+    },
+    InvalidPauliCode {
+        code: u8,
+    },
+    InvalidPauliChar {
+        c: char,
+    },
+    WrongLength {
+        expected: usize,
+        got: usize,
+    },
+    NegativeRate {
+        index: usize,
+        rate: f64,
+    },
+    EmptyLincomb {
+        index: usize,
+    },
+    Internal(String),
+}
+
+impl fmt::Display for Error {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            Error::TooManyQubits { got } => {
+                write!(
+                    f,
+                    "LindbladSpec supports n_qubits ≤ {MAX_QUBITS}; got {got}"
+                )
+            }
+            Error::LengthMismatch { what, a, b } => {
+                write!(f, "{what}: expected matching lengths, got {a} and {b}")
+            }
+            Error::InvalidPauliCode { code } => write!(
+                f,
+                "Pauli code must be 0 (I), 1 (X), 2 (Z), or 3 (Y); got {code}"
+            ),
+            Error::InvalidPauliChar { c } => {
+                write!(f, "invalid Pauli character '{c}'; expected I, X, Y, or Z")
+            }
+            Error::WrongLength { expected, got } => {
+                write!(f, "Pauli string has length {got} but n_qubits = {expected}")
+            }
+            Error::NegativeRate { index, rate } => {
+                write!(f, "jump rate must be non-negative; got γ_{index} = {rate}")
+            }
+            Error::EmptyLincomb { index } => {
+                write!(
+                    f,
+                    "jump {index}: lincomb must contain at least one Pauli term"
+                )
+            }
+            Error::Internal(msg) => write!(f, "internal error: {msg}"),
+        }
+    }
+}
+
+impl std::error::Error for Error {}

crates/ppvm-lindblad/src/expm.rs

@@ -0,0 +1,145 @@
+// SPDX-FileCopyrightText: 2026 The PPVM Authors
+// SPDX-License-Identifier: Apache-2.0
+
+//! Taylor-partition selection for the `quspin-expm`-backed `exp(t·A)·b`
+//! engine (driven from [`crate::mf_expm`]): the `(m, s)` selection tables
+//! [`THETA`] / [`THETA_LOOSE`] from Al-Mohy & Higham (2011), used to pick
+//! the partition handed to `quspin-expm`'s `from_parts`.
+
+/// `θ_m` table from Al-Mohy & Higham (2011), Table A.3, for double
+/// precision (unit roundoff `u = 2^{-53}`).
+///
+/// `θ_m` bounds `‖A‖₁` such that the degree-`m` Taylor polynomial
+/// approximates `exp(A)` to within `u`. We pick `(m, s)` with
+/// `s ≥ ⌈‖tA‖₁ / θ_m⌉` and minimise `m·s` (total SpMV count).
+pub(crate) const THETA: &[(u32, f64)] = &[
+    (1, 2.29e-16),
+    (2, 2.58e-8),
+    (3, 1.39e-5),
+    (4, 3.40e-4),
+    (5, 2.40e-3),
+    (6, 9.07e-3),
+    (7, 2.38e-2),
+    (8, 5.00e-2),
+    (9, 8.96e-2),
+    (10, 1.44e-1),
+    (11, 2.14e-1),
+    (12, 3.00e-1),
+    (13, 4.00e-1),
+    (14, 5.14e-1),
+    (15, 6.41e-1),
+    (16, 7.81e-1),
+    (17, 9.31e-1),
+    (18, 1.09),
+    (19, 1.26),
+    (20, 1.44),
+    (21, 1.62),
+    (22, 1.82),
+    (23, 2.01),
+    (24, 2.22),
+    (25, 2.43),
+    (26, 2.64),
+    (27, 2.86),
+    (28, 3.08),
+    (29, 3.31),
+    (30, 3.54),
+];
+
+/// `θ_m` table for a relaxed backward-error tolerance `tol = 1e-6`, computed
+/// with the same Al-Mohy & Higham (2011) construction as [`THETA`] (the
+/// backward-error series `h_{m+1}(x) = log(e^{-x} T_m(x))`; validated by
+/// reproducing the `u = 2^{-53}` table above to ~2 significant figures).
+///
+/// The predictor-corrector truncates the Pauli basis at `drop_tol` (typically
+/// 1e-3), so computing `exp` to double-precision backward error (~1e-16) is
+/// ~10 orders more accurate than the state it acts on. Using `tol = 1e-6`
+/// (still ~1000x tighter than the truncation) admits a lower-degree Taylor
+/// polynomial for the same `‖tA‖`, cutting the SpMV count (e.g. 23 -> 13 at
+/// `‖tA‖ ≈ 2`) with no measurable effect on the truncated result.
+pub(crate) const THETA_LOOSE: &[(u32, f64)] = &[
+    (1, 2.000e-06),
+    (2, 2.447e-03),
+    (3, 2.863e-02),
+    (4, 1.025e-01),
+    (5, 2.262e-01),
+    (6, 3.911e-01),
+    (7, 5.866e-01),
+    (8, 8.045e-01),
+    (9, 1.039),
+    (10, 1.285),
+    (11, 1.539),
+    (12, 1.801),
+    (13, 2.067),
+    (14, 2.337),
+    (15, 2.610),
+    (16, 2.885),
+    (17, 3.162),
+    (18, 3.441),
+    (19, 3.721),
+    (20, 4.001),
+    (21, 4.282),
+    (22, 4.564),
+    (23, 4.847),
+    (24, 5.129),
+    (25, 5.412),
+    (26, 5.696),
+    (27, 5.979),
+    (28, 6.263),
+    (29, 6.546),
+    (30, 6.830),
+];
+
+/// Pick `(m, s)` minimising `s·m` subject to `s ≥ ⌈t_norm / θ_m⌉, s ≥ 1`,
+/// using the `θ_m` table `theta`. Restricted to the table's `m` range; for
+/// larger norms `s` simply grows linearly.
+fn select_ms_with(t_norm: f64, theta: &[(u32, f64)]) -> (u32, u32) {
+    if t_norm <= 0.0 {
+        return (1, 1);
+    }
+    let mut best_m = 1u32;
+    let mut best_s = 1u32;
+    let mut best_cost = u64::MAX;
+    for &(m, th) in theta {
+        let s_f = (t_norm / th).ceil();
+        let s = if s_f >= 1.0 { s_f as u32 } else { 1 };
+        let cost = (m as u64) * (s as u64);
+        if cost < best_cost {
+            best_cost = cost;
+            best_m = m;
+            best_s = s;
+        }
+    }
+    (best_m, best_s)
+}
+
+/// `(m, s)` selection at double-precision backward error ([`THETA`]).
+pub(crate) fn select_ms(t_norm: f64) -> (u32, u32) {
+    select_ms_with(t_norm, THETA)
+}
+
+/// `(m, s)` selection at the relaxed `tol = 1e-6` backward error
+/// ([`THETA_LOOSE`]) — fewer SpMVs, used on the truncated PC expm path.
+pub(crate) fn select_ms_loose(t_norm: f64) -> (u32, u32) {
+    select_ms_with(t_norm, THETA_LOOSE)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn ms_selection_sane() {
+        // tiny norm → small m, s = 1
+        let (m, s) = select_ms(1e-9);
+        assert!(m <= 5, "expected small m for tiny norm, got m={m}");
+        assert_eq!(s, 1);
+
+        // moderate norm → m·s should be ~10-50
+        let (m, s) = select_ms(1.0);
+        assert!((m * s) <= 50, "moderate norm cost too high: m={m} s={s}");
+
+        // large norm → s grows
+        let (_m, s) = select_ms(100.0);
+        assert!(s >= 20, "large norm should require many steps, got s={s}");
+    }
+}