From 7b52cf6b803d0b8ac33acfcca6748fd6136bbc3d Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 16 Sep 2026 09:58:17 +0200 Subject: [PATCH 1/7] feat(ppvm-lindblad): Kossakowski-form dissipator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add the general GKSL form `D*(O) = Σ_nm K_nm (A_n† O A_m − ½{A_n† A_m, O})`, alongside the existing jump-operator form. Each `(n, m)` entry is compiled once into a `kossakowski::Pair`; Hermitian-conjugate pairs are folded into a single upper-triangle entry, since both sandwiches produce the same output words with conjugate phase and the action keeps only the real part. `PairShape` distinguishes the diagonal `(n,n)` entry from a folded off-diagonal one and carries the term list only the folded case needs, so the off-diagonal-only path cannot be reached on a diagonal pair. `precompute_ldagger_l` generalizes to `precompute_adag_b(a, b)` and moves to `algebra` next to `PauliTerm`, which the new module shares. Co-Authored-By: Claude Opus 5 --- Cargo.lock | 188 ++++++++++++++ crates/ppvm-lindblad/Cargo.toml | 13 + crates/ppvm-lindblad/src/algebra.rs | 46 +++- crates/ppvm-lindblad/src/error.rs | 19 ++ crates/ppvm-lindblad/src/kossakowski.rs | 330 ++++++++++++++++++++++++ crates/ppvm-lindblad/src/lib.rs | 1 + crates/ppvm-lindblad/src/spec.rs | 65 +++-- 7 files changed, 628 insertions(+), 34 deletions(-) create mode 100644 crates/ppvm-lindblad/src/kossakowski.rs diff --git a/Cargo.lock b/Cargo.lock index 110a5f3c1..59ad47df8 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1343,6 +1343,114 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "glam" +version = "0.14.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "333928d5eb103c5d4050533cec0384302db6be8ef7d3cebd30ec6a35350353da" + +[[package]] +name = "glam" +version = "0.15.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3abb554f8ee44336b72d522e0a7fe86a29e09f839a36022fa869a7dfe941a54b" + +[[package]] +name = "glam" +version = "0.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4126c0479ccf7e8664c36a2d719f5f2c140fbb4f9090008098d2c291fa5b3f16" + +[[package]] +name = "glam" +version = "0.17.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e01732b97afd8508eee3333a541b9f7610f454bb818669e66e90f5f57c93a776" + +[[package]] +name = "glam" +version = "0.18.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "525a3e490ba77b8e326fb67d4b44b4bd2f920f44d4cc73ccec50adc68e3bee34" + +[[package]] +name = "glam" +version = "0.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b8509e6791516e81c1a630d0bd7fbac36d2fa8712a9da8662e716b52d5051ca" + +[[package]] +name = "glam" +version = "0.20.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f43e957e744be03f5801a55472f593d43fabdebf25a4585db250f04d86b1675f" + +[[package]] +name = "glam" +version = "0.21.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "518faa5064866338b013ff9b2350dc318e14cc4fcd6cb8206d7e7c9886c98815" + +[[package]] +name = "glam" +version = "0.22.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "12f597d56c1bd55a811a1be189459e8fad2bbc272616375602443bdfb37fa774" + +[[package]] +name = "glam" +version = "0.23.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e4afd9ad95555081e109fe1d21f2a30c691b5f0919c67dfa690a2e1eb6bd51c" + +[[package]] +name = "glam" +version = "0.24.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b5418c17512bdf42730f9032c74e1ae39afc408745ebb2acf72fbc4691c17945" + +[[package]] +name = "glam" +version = "0.25.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "151665d9be52f9bb40fc7966565d39666f2d1e69233571b71b87791c7e0528b3" + +[[package]] +name = "glam" +version = "0.27.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9e05e7e6723e3455f4818c7b26e855439f7546cf617ef669d1adedb8669e5cb9" + +[[package]] +name = "glam" +version = "0.28.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "779ae4bf7e8421cf91c0b3b64e7e8b40b862fba4d393f59150042de7c4965a94" + +[[package]] +name = "glam" +version = "0.29.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8babf46d4c1c9d92deac9f7be466f76dfc4482b6452fc5024b5e8daf6ffeb3ee" + +[[package]] +name = "glam" +version = "0.30.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19fc433e8437a212d1b6f1e68c7824af3aed907da60afa994e7f542d18d12aa9" + +[[package]] +name = "glam" +version = "0.31.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "556f6b2ea90b8d15a74e0e7bb41671c9bdf38cd9f78c284d750b9ce58a2b5be7" + +[[package]] +name = "glam" +version = "0.32.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f70749695b063ecbf6b62949ccccde2e733ec3ecbbd71d467dca4e5c6c97cca0" + [[package]] name = "gxhash" version = "3.5.0" @@ -1738,6 +1846,51 @@ dependencies = [ "windows-sys 0.61.2", ] +[[package]] +name = "nalgebra" +version = "0.34.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "df76ea0ff5c7e6b88689085804d6132ded0ddb9de5ca5b8aeb9eeadc0508a70a" +dependencies = [ + "approx", + "glam 0.14.0", + "glam 0.15.2", + "glam 0.16.0", + "glam 0.17.3", + "glam 0.18.0", + "glam 0.19.0", + "glam 0.20.5", + "glam 0.21.3", + "glam 0.22.0", + "glam 0.23.0", + "glam 0.24.2", + "glam 0.25.0", + "glam 0.27.0", + "glam 0.28.0", + "glam 0.29.3", + "glam 0.30.10", + "glam 0.31.1", + "glam 0.32.1", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "nalgebra-macros" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "973e7178a678cfd059ccec50887658d482ce16b0aa9da3888ddeab5cd5eb4889" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.118", +] + [[package]] name = "ndarray" version = "0.17.2" @@ -2036,7 +2189,10 @@ dependencies = [ name = "ppvm-lindblad" version = "0.1.0" dependencies = [ + "approx", + "criterion 0.7.0", "fxhash", + "nalgebra", "ndarray", "num", "ppvm-pauli-sum", @@ -2730,6 +2886,15 @@ version = "1.0.20" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "28d3b2b1366ec20994f1fd18c3c594f05c5dd4bc44d8bb0c1c632c8d6829481f" +[[package]] +name = "safe_arch" +version = "0.7.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "96b02de82ddbe1b636e6170c21be622223aea188ef2e139be0a5b219ec215323" +dependencies = [ + "bytemuck", +] + [[package]] name = "same-file" version = "1.0.6" @@ -2893,6 +3058,19 @@ dependencies = [ "libc", ] +[[package]] +name = "simba" +version = "0.9.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c99284beb21666094ba2b75bbceda012e610f5479dfcc2d6e2426f53197ffd95" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "paste", + "wide", +] + [[package]] name = "similar" version = "2.7.0" @@ -3563,6 +3741,16 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "wide" +version = "0.7.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ce5da8ecb62bcd8ec8b7ea19f69a51275e91299be594ea5cc6ef7819e16cd03" +dependencies = [ + "bytemuck", + "safe_arch", +] + [[package]] name = "winapi" version = "0.3.9" diff --git a/crates/ppvm-lindblad/Cargo.toml b/crates/ppvm-lindblad/Cargo.toml index 900388d36..25bd6dd32 100644 --- a/crates/ppvm-lindblad/Cargo.toml +++ b/crates/ppvm-lindblad/Cargo.toml @@ -19,3 +19,16 @@ quspin-expm = { git = "https://github.com/QuSpin/QuSpin-rust", rev = "a0ad6c9fe2 # 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" +criterion = "0.7.0" +nalgebra = "0.34" + +[[bench]] +name = "kossakowski" +harness = false + +[[bench]] +name = "drug_dipolar" +harness = false diff --git a/crates/ppvm-lindblad/src/algebra.rs b/crates/ppvm-lindblad/src/algebra.rs index da44c3042..7a2b8c177 100644 --- a/crates/ppvm-lindblad/src/algebra.rs +++ b/crates/ppvm-lindblad/src/algebra.rs @@ -10,10 +10,54 @@ //! keeps a copy that returns the unpacked `(word, phase)` pair without //! constructing a phased wrapper. -use crate::word::{W_CHUNKS, Word}; +use crate::word::{Chunk, W_CHUNKS, Word}; +use fxhash::FxHashMap; use num::Complex; use ppvm_traits::PauliWordTrait; +/// Magnitude below which an expanded Pauli coefficient is treated as +/// cancellation noise and dropped. +pub(crate) const COEFF_DROP_TOL: f64 = 1e-14; + +/// One Pauli term in a complex linear combination (a single summand of +/// `L = Σ_a λ_a P_a`, or of a precomputed product such as `L†L`). +#[derive(Clone)] +pub(crate) struct PauliTerm { + pub(crate) word: Word, + pub(crate) coeff: Complex, +} + +/// Expand `A†B = (Σ_a λ_a P_a)† (Σ_b μ_b P_b) = Σ_{a,b} λ_a* μ_b P_a P_b` +/// as a Pauli linear combination, dropping FP-noise zeros. For `A = B` +/// (the jump-operator `L†L`) the coefficients are real; in general they +/// are complex. +pub(crate) fn precompute_adag_b(a_terms: &[PauliTerm], b_terms: &[PauliTerm]) -> Vec { + let zero = Complex::new(0.0, 0.0); + let mut acc: FxHashMap> = FxHashMap::default(); + for a in a_terms { + for b in b_terms { + let (word, phase) = pauli_mul(&a.word, &b.word); + let coeff = a.coeff.conj() * b.coeff * phase_factor(phase); + *acc.entry(word).or_insert(zero) += coeff; + } + } + acc.into_iter() + .filter(|(_, c)| c.norm() > COEFF_DROP_TOL) + .map(|(word, coeff)| PauliTerm { word, coeff }) + .collect() +} + +/// Union of the supports (`xbits | zbits`) of every term, as raw chunks. +pub(crate) fn support_mask(terms: &[PauliTerm]) -> [Chunk; W_CHUNKS] { + let mut mask = [0 as Chunk; W_CHUNKS]; + for t in terms { + for (i, slot) in mask.iter_mut().enumerate() { + *slot |= t.word.xbits.data[i] | t.word.zbits.data[i]; + } + } + mask +} + #[inline(always)] pub(crate) fn phase_factor(phase: u8) -> Complex { match phase & 3 { diff --git a/crates/ppvm-lindblad/src/error.rs b/crates/ppvm-lindblad/src/error.rs index e0c1eb5f0..6c2235e86 100644 --- a/crates/ppvm-lindblad/src/error.rs +++ b/crates/ppvm-lindblad/src/error.rs @@ -34,6 +34,17 @@ pub enum Error { EmptyLincomb { index: usize, }, + /// Row `row` of the Kossakowski matrix is not `n_ops` wide. + KMatrixRowLength { + row: usize, + expected: usize, + got: usize, + }, + /// `K_nm ≠ conj(K_mn)`: not a valid GKSL pair matrix. + KMatrixNotHermitian { + n: usize, + m: usize, + }, Internal(String), } @@ -68,6 +79,14 @@ impl fmt::Display for Error { "jump {index}: lincomb must contain at least one Pauli term" ) } + Error::KMatrixRowLength { row, expected, got } => write!( + f, + "kossakowski K row {row} has length {got}; expected {expected} (one per operator)" + ), + Error::KMatrixNotHermitian { n, m } => write!( + f, + "kossakowski K must be Hermitian; K[{n}][{m}] ≠ conj(K[{m}][{n}])" + ), Error::Internal(msg) => write!(f, "internal error: {msg}"), } } diff --git a/crates/ppvm-lindblad/src/kossakowski.rs b/crates/ppvm-lindblad/src/kossakowski.rs new file mode 100644 index 000000000..65aa76723 --- /dev/null +++ b/crates/ppvm-lindblad/src/kossakowski.rs @@ -0,0 +1,330 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Kossakowski-form dissipator. +//! +//! For a family of operators `A_n` and a Hermitian, positive-semidefinite +//! pair matrix `K`, the adjoint dissipator is +//! +//! ```text +//! D*(O) = Σ_{n,m} K_nm ( A_n† O A_m − ½ {A_n† A_m, O} ). +//! ``` +//! +//! This is the general GKSL form; the diagonal `K = diag(γ_k)` case is the +//! jump-operator form handled by [`crate::spec::JumpKind::General`]. +//! +//! Each `(n, m)` pair is compiled once into a [`Pair`]. Hermitian-conjugate +//! pairs `(n,m)` and `(m,n)` are *folded* into a single upper-triangle entry: +//! both sandwiches produce the same output words with conjugate phase and the +//! final action keeps only the real part, so one entry that doubles-and-takes- +//! `Re` suffices, halving the pair count. [`PairShape`] records which of the +//! two a compiled pair is, and carries the extra term list that only the +//! folded case needs. + +use crate::Error; +use crate::algebra::{ + COEFF_DROP_TOL, PauliTerm, comm_product, pauli_mul, phase_factor, precompute_adag_b, + support_mask, +}; +use crate::word::{Chunk, W_CHUNKS, Word, parse_pauli_string}; +use fxhash::FxHashMap; +use num::Complex; +use std::collections::BTreeSet; + +/// Relative tolerance for the Hermiticity check on `K`. +const HERMITICITY_TOL: f64 = 1e-10; + +/// Sandwich table of a pair, grouped by the left word: one +/// `(P_a, [(P_b, coeff), …])` group per distinct `P_a`, so `P_a · p` is +/// computed once per group and reused across its `P_b` partners. +type SandwichGroups = Vec<(Word, Vec<(Word, Complex)>)>; + +/// Which of the two compiled pair forms a [`Pair`] is. +/// +/// The distinction changes the meaning of [`Pair::dd`] and selects the term +/// list used by the one-sided commutator path, so it is modelled as a sum +/// type rather than a flag: the off-diagonal-only term list cannot be +/// reached on a diagonal pair. +pub(crate) enum PairShape { + /// `n == m`. [`Pair::dd`] is `K_nn · A_n†A_n`. + Diagonal, + /// `n < m`, folding in the conjugate `(m, n)` pair. [`Pair::dd`] is the + /// Hermitian sum `2·Re(K_nm·A_n†A_m)` used by the both-sided + /// anticommutator. + OffDiagonal { + /// The anti-Hermitian difference `−2i·Im(K_nm·A_n†A_m)` + /// (pure-imaginary coefficients), used by the one-sided commutator + /// of the folded conjugate pair. + dd_anti: Vec, + }, +} + +/// One compiled `(n, m)` pair of a Kossakowski dissipator. +pub(crate) struct Pair { + sand: SandwichGroups, + /// `A_n†A_m` scaled by `K_nm`; see [`PairShape`] for the exact form. + dd: Vec, + shape: PairShape, + /// Support masks of `A_n` and `A_m`, for the one-sided fast path. + left_mask: [Chunk; W_CHUNKS], + right_mask: [Chunk; W_CHUNKS], +} + +/// Compile a Kossakowski dissipator into one [`Pair`] per non-negligible +/// upper-triangle entry of `K`. +/// +/// Returns the pairs alongside, for each pair, the union support of its two +/// operators, so the caller can index them by qubit. +pub(crate) fn compile( + ops: &[Vec<(String, Complex)>], + k: &[Vec>], + n_qubits: usize, +) -> Result)>, Error> { + let max_abs = validate_k(k, ops.len())?; + let (parsed, op_support) = parse_ops(ops, n_qubits)?; + + let pair_tol = COEFF_DROP_TOL * max_abs; + let mut out = Vec::new(); + for n in 0..ops.len() { + for m in n..ops.len() { + if k[n][m].norm() <= pair_tol { + continue; + } + let pair = compile_pair(&parsed[n], &parsed[m], k[n][m], n != m); + let mut union: BTreeSet = op_support[n].iter().copied().collect(); + union.extend(op_support[m].iter().copied()); + out.push((pair, union.into_iter().collect())); + } + } + Ok(out) +} + +/// Check that `k` is square with side `n_ops` and Hermitian. Returns the +/// largest `|K_nm|`, which sets the scale for the negligible-pair cutoff. +fn validate_k(k: &[Vec>], n_ops: usize) -> Result { + if k.len() != n_ops { + return Err(Error::LengthMismatch { + what: "kossakowski ops and K rows", + a: n_ops, + b: k.len(), + }); + } + for (row, entries) in k.iter().enumerate() { + if entries.len() != n_ops { + return Err(Error::KMatrixRowLength { + row, + expected: n_ops, + got: entries.len(), + }); + } + } + + let max_abs = k + .iter() + .flat_map(|row| row.iter().map(|c| c.norm())) + .fold(0.0_f64, f64::max); + + // A non-Hermitian K is not a valid GKSL pair matrix and would produce an + // action that does not preserve Hermiticity. + let tol = HERMITICITY_TOL * max_abs.max(1.0); + for (n, row_n) in k.iter().enumerate() { + for (m, k_nm) in row_n.iter().enumerate().skip(n) { + if (k_nm - k[m][n].conj()).norm() > tol { + return Err(Error::KMatrixNotHermitian { n, m }); + } + } + } + Ok(max_abs) +} + +/// Parsed operator table: the Pauli terms of each `A_n`, and each `A_n`'s +/// union support. +type ParsedOps = (Vec>, Vec>); + +/// Parse each operator's Pauli lincomb, returning the parsed terms and each +/// operator's union support. +fn parse_ops(ops: &[Vec<(String, Complex)>], n_qubits: usize) -> Result { + let mut parsed = Vec::with_capacity(ops.len()); + let mut supports = Vec::with_capacity(ops.len()); + for (i, op) in ops.iter().enumerate() { + if op.is_empty() { + return Err(Error::EmptyLincomb { index: i }); + } + let mut terms = Vec::with_capacity(op.len()); + let mut union: BTreeSet = BTreeSet::new(); + for (s, c) in op { + let (word, support) = parse_pauli_string(s, n_qubits)?; + union.extend(support.iter().copied()); + terms.push(PauliTerm { word, coeff: *c }); + } + parsed.push(terms); + supports.push(union.into_iter().collect()); + } + Ok((parsed, supports)) +} + +/// Compile the `(n, m)` entry with `A_n = a_terms`, `A_m = b_terms`. +fn compile_pair( + a_terms: &[PauliTerm], + b_terms: &[PauliTerm], + k_nm: Complex, + off_diag: bool, +) -> Pair { + // A_n†A_m as `Σ γ_w W`, then scaled by K_nm. + let adag_b = precompute_adag_b(a_terms, b_terms); + let (dd, shape) = if off_diag { + // Splitting K_nm·γ_w into its Hermitian and anti-Hermitian halves is + // what lets the conjugate (m,n) pair be dropped: the (m,n) sandwich + // contributes the complex conjugate, so the sum is 2·Re on the + // both-sided path and 2i·Im on the one-sided one. + let mut dd = Vec::with_capacity(adag_b.len()); + let mut dd_anti = Vec::with_capacity(adag_b.len()); + for t in &adag_b { + let c = k_nm * t.coeff; + if c.re.abs() > COEFF_DROP_TOL { + dd.push(PauliTerm { + word: t.word, + coeff: Complex::new(2.0 * c.re, 0.0), + }); + } + if c.im.abs() > COEFF_DROP_TOL { + dd_anti.push(PauliTerm { + word: t.word, + coeff: Complex::new(0.0, -2.0 * c.im), + }); + } + } + (dd, PairShape::OffDiagonal { dd_anti }) + } else { + let dd = adag_b + .iter() + .map(|t| PauliTerm { + word: t.word, + coeff: k_nm * t.coeff, + }) + .collect(); + (dd, PairShape::Diagonal) + }; + + let sand = a_terms + .iter() + .map(|a| { + let rights = b_terms + .iter() + .map(|b| (b.word, a.coeff.conj() * b.coeff * k_nm)) + .collect(); + (a.word, rights) + }) + .collect(); + + Pair { + sand, + dd, + shape, + left_mask: support_mask(a_terms), + right_mask: support_mask(b_terms), + } +} + +impl Pair { + /// Accumulate this pair's contribution to `L*(p)` into `local`. + pub(crate) fn accumulate(&self, p: &Word, local: &mut FxHashMap>) { + let mut p_bits = [0 as Chunk; W_CHUNKS]; + for (i, slot) in p_bits.iter_mut().enumerate() { + *slot = p.xbits.data[i] | p.zbits.data[i]; + } + let hits = |mask: &[Chunk; W_CHUNKS]| (0..W_CHUNKS).any(|i| mask[i] & p_bits[i] != 0); + let (hit_l, hit_r) = (hits(&self.left_mask), hits(&self.right_mask)); + + // The pair is only visited when `p` overlaps at least one side, so + // "not both" means exactly one. + if hit_l && hit_r { + self.accumulate_both_sided(p, local); + } else { + self.accumulate_one_sided(p, hit_r, local); + } + } + + /// One-sided fast path: when `p` is disjoint from one of the two + /// operators the sandwich and anticommutator collapse to a commutator. + /// For a diagonal pair with `D = K·A_n†A_m`: + /// + /// ```text + /// p disjoint from A_n (left): C = −½ [D, p] + /// p disjoint from A_m (right): C = +½ [D, p] + /// ``` + /// + /// For a folded off-diagonal pair the two conjugate one-sided + /// contributions combine into `±½ [F, p]` with the anti-Hermitian + /// `F = dd_anti` and the *opposite* sign. With `[P_c, p] = −i·eps·out` + /// from [`comm_product`], the term coefficient is `∓ t_c · (i/2) · eps`. + fn accumulate_one_sided( + &self, + p: &Word, + hit_r: bool, + local: &mut FxHashMap>, + ) { + let zero = Complex::new(0.0, 0.0); + let (terms, half_i) = match &self.shape { + PairShape::Diagonal => ( + &self.dd, + if hit_r { + Complex::new(0.0, 0.5) + } else { + Complex::new(0.0, -0.5) + }, + ), + PairShape::OffDiagonal { dd_anti } => ( + dd_anti, + if hit_r { + Complex::new(0.0, -0.5) + } else { + Complex::new(0.0, 0.5) + }, + ), + }; + for t in terms { + let (out, eps) = comm_product(&t.word, p); + if eps != 0.0 { + *local.entry(out).or_insert(zero) += t.coeff * half_i * eps; + } + } + } + + /// Both sides hit: full sandwich plus anticommutator. The sandwich is + /// grouped by the left word so `P_a · p` is computed once per distinct + /// `P_a` and reused across all its `P_b` partners. For a folded + /// off-diagonal pair the sandwich is doubled and its real part taken + /// (the conjugate `(m,n)` pair supplies the other half). + fn accumulate_both_sided(&self, p: &Word, local: &mut FxHashMap>) { + let zero = Complex::new(0.0, 0.0); + let fold = matches!(self.shape, PairShape::OffDiagonal { .. }); + for (wa, rights) in &self.sand { + let (r_ap, phi1) = pauli_mul(wa, p); + // Hoisted out of the inner loop: the fold is a property of the + // pair, not of the term. + if fold { + for (wb, c0) in rights { + let (s, phi2) = pauli_mul(&r_ap, wb); + let v = c0 * phase_factor(phi1 + phi2); + *local.entry(s).or_insert(zero) += Complex::new(2.0 * v.re, 0.0); + } + } else { + for (wb, c0) in rights { + let (s, phi2) = pauli_mul(&r_ap, wb); + *local.entry(s).or_insert(zero) += c0 * phase_factor(phi1 + phi2); + } + } + } + + // −½{D, p}. For Pauli words, {P_c, p} = 2·sign·R when they commute + // (P_c·p = sign·R) and 0 when they anti-commute; the ½ cancels the 2. + for t in &self.dd { + let (r, phase) = pauli_mul(&t.word, p); + if phase & 1 == 0 { + let sign = if phase == 0 { 1.0 } else { -1.0 }; + *local.entry(r).or_insert(zero) -= t.coeff * Complex::new(sign, 0.0); + } + } + } +} diff --git a/crates/ppvm-lindblad/src/lib.rs b/crates/ppvm-lindblad/src/lib.rs index 57dff8aa2..55c2feb92 100644 --- a/crates/ppvm-lindblad/src/lib.rs +++ b/crates/ppvm-lindblad/src/lib.rs @@ -38,6 +38,7 @@ mod basis; pub mod config; pub mod error; pub(crate) mod expm; +mod kossakowski; mod scalar; pub mod sector; mod spec; diff --git a/crates/ppvm-lindblad/src/spec.rs b/crates/ppvm-lindblad/src/spec.rs index 4a8046eea..78b48be05 100644 --- a/crates/ppvm-lindblad/src/spec.rs +++ b/crates/ppvm-lindblad/src/spec.rs @@ -4,7 +4,10 @@ //! Precompiled Lindbladian: construction and the single-Pauli `L*` kernel. use crate::Error; -use crate::algebra::{anti_commutes, comm_product, pauli_mul, phase_factor}; +use crate::algebra::{ + PauliTerm, anti_commutes, comm_product, pauli_mul, phase_factor, precompute_adag_b, +}; +use crate::kossakowski; use crate::word::{MAX_QUBITS, Word, parse_pauli_string, word_support}; use fxhash::FxHashMap; use num::Complex; @@ -16,17 +19,9 @@ struct HTerm { coeff: f64, } -/// One Pauli term in a complex linear combination (a single summand of -/// `L = Σ_a λ_a P_a` or of the precomputed `L†L`). -#[derive(Clone)] -struct PauliTerm { - word: Word, - coeff: Complex, -} - -/// One jump operator `L_k` with rate `γ_k`. The `HermitianPauli` variant -/// is a fast path; `General` handles arbitrary complex Pauli sums. -#[derive(Clone)] +/// One entry of the dissipator. `HermitianPauli` and `General` are the +/// jump-operator form (`K` diagonal); `Kossakowski` is one compiled pair of +/// the general form. See [`crate::kossakowski`]. enum JumpKind { HermitianPauli { word: Word, @@ -37,26 +32,7 @@ enum JumpKind { dagger_dagger: Vec, // L†L = Σ_c μ_c P_c (μ_c ∈ ℝ) rate: f64, }, -} - -/// Expand `L†L = (Σ_a λ_a P_a)† (Σ_b λ_b P_b) = Σ_{a,b} λ_a* λ_b P_a P_b` -/// as a Pauli linear combination, dropping FP-noise zeros. Coefficients are -/// real because `L†L` is Hermitian; we keep them complex for arithmetic -/// uniformity. -fn precompute_ldagger_l(terms: &[PauliTerm]) -> Vec { - let zero = Complex::new(0.0, 0.0); - let mut acc: FxHashMap> = FxHashMap::default(); - for a in terms { - for b in terms { - let (word, phase) = pauli_mul(&a.word, &b.word); - let coeff = a.coeff.conj() * b.coeff * phase_factor(phase); - *acc.entry(word).or_insert(zero) += coeff; - } - } - acc.into_iter() - .filter(|(_, c)| c.norm() > 1e-14) - .map(|(word, coeff)| PauliTerm { word, coeff }) - .collect() + Kossakowski(kossakowski::Pair), } /// Union of `index[q]` for each `q ∈ p_support`, deduped. @@ -163,7 +139,7 @@ impl LindbladSpec { for q in union_support { j_support_idx[q as usize].push(k as u32); } - let dagger_dagger = precompute_ldagger_l(&terms); + let dagger_dagger = precompute_adag_b(&terms, &terms); j_kinds.push(JumpKind::General { terms, dagger_dagger, @@ -180,6 +156,28 @@ impl LindbladSpec { }) } + /// Add a Kossakowski-form dissipator + /// `D*(O) = Σ_{n,m} K_nm ( A_n† O A_m − ½ {A_n† A_m, O} )`. + /// + /// `ops` lists the operators `A_n` as complex Pauli linear combinations; + /// `k` is the `n_ops × n_ops` Hermitian pair matrix. Contributions are + /// added to any jumps already present. See [`crate::kossakowski`] for + /// how each `(n, m)` entry is compiled. + pub fn add_kossakowski( + &mut self, + ops: &[Vec<(String, Complex)>], + k: &[Vec>], + ) -> Result<(), Error> { + for (pair, support) in kossakowski::compile(ops, k, self.n_qubits)? { + let idx = self.j_kinds.len() as u32; + self.j_kinds.push(JumpKind::Kossakowski(pair)); + for q in support { + self.j_support[q as usize].push(idx); + } + } + Ok(()) + } + pub fn n_qubits(&self) -> usize { self.n_qubits } @@ -263,6 +261,7 @@ impl LindbladSpec { } } } + JumpKind::Kossakowski(pair) => pair.accumulate(p, local), } } From e6266b52e085f481b71ab88222a90970c683085e Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 16 Sep 2026 09:58:50 +0200 Subject: [PATCH 2/7] bench(lindblad): Kossakowski superradiance and dipolar-relaxation workloads `kossakowski` compares the pair-matrix path against the equivalent eigenmode-jump decomposition on a superradiant chain; `drug_dipolar` covers the 2-local rank-2 tensors of dipolar relaxation, where each `A_n` carries several Pauli terms. `drug_profile` is the same model wired to a plain timing loop for use under a profiler. Co-Authored-By: Claude Opus 5 --- crates/ppvm-lindblad/benches/drug_dipolar.rs | 261 ++++++++++++++++++ crates/ppvm-lindblad/benches/kossakowski.rs | 173 ++++++++++++ crates/ppvm-lindblad/examples/drug_profile.rs | 217 +++++++++++++++ 3 files changed, 651 insertions(+) create mode 100644 crates/ppvm-lindblad/benches/drug_dipolar.rs create mode 100644 crates/ppvm-lindblad/benches/kossakowski.rs create mode 100644 crates/ppvm-lindblad/examples/drug_profile.rs diff --git a/crates/ppvm-lindblad/benches/drug_dipolar.rs b/crates/ppvm-lindblad/benches/drug_dipolar.rs new file mode 100644 index 000000000..3e1281e62 --- /dev/null +++ b/crates/ppvm-lindblad/benches/drug_dipolar.rs @@ -0,0 +1,261 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Per-step cost of the Kossakowski-form dissipator on a *molecular dipolar +//! relaxation* workload (the ZULF-NMR drug-FID application), which stresses +//! the path differently from the superradiance chain in `kossakowski.rs`: +//! +//! - operators are **2-local rank-2 tensors** with ~4 Pauli terms each (one +//! channel per site pair, per spatial harmonic `m`), not single-site σ⁻; +//! - `K` is **block-diagonal** over the 5 spatial components `m ∈ −2..2`, +//! each block a dense `P×P` Gram matrix over the `P = C(N,2)` pairs; +//! - the basis strings are dense/high-weight, so most candidate pairs hit +//! the **both-sided sandwich** (12-product) path — the arm the 2026-07-16 +//! ledger flagged as remaining headroom. +//! +//! Both specs (eigenmode jumps vs Kossakowski) generate the identical action; +//! the benchmark measures representation cost on one full `pc_step`. + +use criterion::{Criterion, criterion_group, criterion_main}; +use num::Complex; +use ppvm_lindblad::{JumpInput, LindbladSpec, PcStepConfig, Word, parse_pauli_string}; +use std::hint::black_box; + +const B: usize = 4096; +const GROW_STEPS: usize = 3; +const DT: f64 = 1e-3; +const N_M: usize = 5; + +/// Dipolar-pair geometry: the `(a, b)` site pairs, each pair's coupling +/// magnitude, and its unit separation vector. +type Geometry = (Vec<(usize, usize)>, Vec, Vec<[f64; 3]>); + +/// A Kossakowski dissipator as handed to `LindbladSpec::add_kossakowski`: +/// the operators `A_n` as Pauli lincombs, and the pair matrix `K`. +type KossakowskiModel = (Vec)>>, Vec>>); +// spatial harmonics m = -2..2 + +fn pstr(n: usize, sites: &[(usize, char)]) -> String { + let mut s = vec!['I'; n]; + for &(q, c) in sites { + s[q] = c; + } + s.into_iter().collect() +} + +/// Deterministic pseudo-random 3D unit direction + distance for pair (a,b), +/// so the model is reproducible without an RNG dependency in the bench. +fn hashf(mut x: u64) -> f64 { + x ^= x >> 33; + x = x.wrapping_mul(0xff51afd7ed558ccd); + x ^= x >> 33; + (x >> 11) as f64 / (1u64 << 53) as f64 +} + +/// Dipolar coupling `b` and unit vector for every pair, from placing spins on +/// a jittered chain (real molecules: `b ∝ 1/r³`, generic directions). +fn geometry(n: usize) -> Geometry { + let pos: Vec<[f64; 3]> = (0..n) + .map(|i| { + [ + i as f64 + 0.3 * hashf(i as u64 * 3 + 1), + 0.4 * hashf(i as u64 * 3 + 2), + 0.4 * hashf(i as u64 * 3 + 3), + ] + }) + .collect(); + let mut pairs = Vec::new(); + let mut bmag = Vec::new(); + let mut dir = Vec::new(); + for a in 0..n { + for b in (a + 1)..n { + let d = [ + pos[a][0] - pos[b][0], + pos[a][1] - pos[b][1], + pos[a][2] - pos[b][2], + ]; + let r = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt(); + pairs.push((a, b)); + bmag.push(1.0 / (r * r * r)); + dir.push([d[0] / r, d[1] / r, d[2] / r]); + } + } + (pairs, bmag, dir) +} + +/// Real rank-2 spherical harmonics (up to normalization) of a unit vector, +/// ordered m = -2,-1,0,1,2 — the spatial factors that make `Γ` rank 5. +fn y2(u: &[f64; 3]) -> [f64; N_M] { + let (x, y, z) = (u[0], u[1], u[2]); + [ + x * y, + y * z, + (3.0 * z * z - 1.0) / 2.0, + x * z, + (x * x - y * y) / 2.0, + ] +} + +/// The 2-local rank-2 tensor operator on pair (a,b) for tensor component +/// `mt` — a representative 4-term Pauli lincomb matching the high-field +/// dressed-tensor forms (T^(2,±2): XX∓YY ± i(XY±YX), etc.). The exact +/// coefficients are immaterial to the cost profile; the term *count* and +/// 2-locality are what matter. +fn tensor_op(n: usize, a: usize, b: usize, mt: usize) -> Vec<(String, Complex)> { + let (i, j) = (Complex::new(0.0, 1.0), Complex::new(1.0, 0.0)); + match mt { + 0 | 4 => { + let s = if mt == 0 { i } else { -i }; // ±2 components + vec![ + (pstr(n, &[(a, 'X'), (b, 'X')]), j), + (pstr(n, &[(a, 'Y'), (b, 'Y')]), -j), + (pstr(n, &[(a, 'X'), (b, 'Y')]), s), + (pstr(n, &[(a, 'Y'), (b, 'X')]), s), + ] + } + 1 | 3 => { + let s = if mt == 1 { i } else { -i }; // ±1 components + vec![ + (pstr(n, &[(a, 'X'), (b, 'Z')]), j), + (pstr(n, &[(a, 'Y'), (b, 'Z')]), s), + (pstr(n, &[(a, 'Z'), (b, 'X')]), j), + (pstr(n, &[(a, 'Z'), (b, 'Y')]), s), + ] + } + _ => vec![ + // m = 0 + (pstr(n, &[(a, 'X'), (b, 'X')]), j), + (pstr(n, &[(a, 'Y'), (b, 'Y')]), j), + (pstr(n, &[(a, 'Z'), (b, 'Z')]), Complex::new(2.0, 0.0)), + ], + } +} + +fn hamiltonian_terms(n: usize, pairs: &[(usize, usize)], bmag: &[f64]) -> Vec<(String, f64)> { + let mut h = Vec::new(); + for (k, &(a, b)) in pairs.iter().enumerate() { + let jc = 0.1 * bmag[k]; // scalar J-coupling, XX+YY + h.push((pstr(n, &[(a, 'X'), (b, 'X')]), jc)); + h.push((pstr(n, &[(a, 'Y'), (b, 'Y')]), jc)); + } + h +} + +/// Kossakowski ops (`N_M` blocks of `P` pair tensors) and the block-diagonal +/// `K = blockdiag(Γ_m)`, `Γ_m[μν] = Σ_{m'} c_μ^{m'} c_ν^{m'}` with +/// `c_μ^{m'} = b_μ Y_2^{m'}(r̂_μ)` — a rank-5 Gram block, exactly as the drug +/// pickles decompose. +fn kossakowski_model(n: usize) -> KossakowskiModel { + let (pairs, bmag, dir) = geometry(n); + let p = pairs.len(); + let c: Vec<[f64; N_M]> = (0..p) + .map(|k| { + let y = y2(&dir[k]); + std::array::from_fn(|mp| bmag[k] * y[mp]) + }) + .collect(); + let mut ops = Vec::with_capacity(N_M * p); + for mt in 0..N_M { + for &(a, b) in &pairs { + ops.push(tensor_op(n, a, b, mt)); + } + } + let m_ops = N_M * p; + let mut k = vec![vec![Complex::new(0.0, 0.0); m_ops]; m_ops]; + for mt in 0..N_M { + let off = mt * p; + for mu in 0..p { + for nu in 0..p { + let g: f64 = (0..N_M).map(|mp| c[mu][mp] * c[nu][mp]).sum(); + k[off + mu][off + nu] = Complex::new(g, 0.0); + } + } + } + (ops, k) +} + +/// Eigenmode jumps of the block-diagonal `K` (the dense representation the +/// Kossakowski path replaces): per block, `L_ν = √γ_ν Σ_μ V_μν T_μ`. +fn eigenmode_jumps(ops: &[Vec<(String, Complex)>], k: &[Vec>]) -> Vec { + let m_ops = ops.len(); + let p = m_ops / N_M; + let mut jumps = Vec::new(); + for mt in 0..N_M { + let off = mt * p; + let block = nalgebra::DMatrix::from_fn(p, p, |a, b| k[off + a][off + b].re); + let eig = nalgebra::SymmetricEigen::new(block); + for nu in 0..p { + let g = eig.eigenvalues[nu]; + if g < 1e-12 { + continue; + } + let mut lin = Vec::new(); + for mu in 0..p { + let v = eig.eigenvectors[(mu, nu)]; + if v.abs() > 1e-14 { + for (s, cc) in &ops[off + mu] { + lin.push((s.clone(), cc * Complex::new(v, 0.0))); + } + } + } + jumps.push(JumpInput { + lincomb: lin, + rate: g, + }); + } + } + jumps +} + +/// Initial observable: γ-weighted transverse magnetization Σ_i X_i (the coil +/// quadrature), a sparse single-site sum like the drug FID initial operator. +fn observable(n: usize) -> (Vec, Vec) { + let mut basis = Vec::new(); + let mut coeffs = Vec::new(); + for a in 0..n { + basis.push(parse_pauli_string(&pstr(n, &[(a, 'X')]), n).unwrap().0); + coeffs.push(1.0); + } + (basis, coeffs) +} + +fn bench_drug(c: &mut Criterion) { + let mut group = c.benchmark_group("pc_step_drug_dipolar"); + group.sample_size(10); + for n in [10usize, 20, 32] { + let (pairs, bmag, _) = geometry(n); + let h = hamiltonian_terms(n, &pairs, &bmag); + let (ops, k) = kossakowski_model(n); + + let spec_eig = LindbladSpec::new(n, &h, &eigenmode_jumps(&ops, &k)).unwrap(); + let mut spec_koss = LindbladSpec::new(n, &h, &[]).unwrap(); + spec_koss.add_kossakowski(&ops, &k).unwrap(); + + let cfg = PcStepConfig { + max_basis: B, + admit_basis: Some(3 * B), + ..Default::default() + }; + let (mut basis, mut coeffs) = observable(n); + for _ in 0..GROW_STEPS { + spec_koss + .pc_step(&mut basis, &mut coeffs, DT, &[], &cfg) + .unwrap(); + } + + for (label, spec) in [("eigenmode", &spec_eig), ("kossakowski", &spec_koss)] { + group.bench_function(format!("{label}_n{n}"), |bch| { + bch.iter(|| { + let mut b = basis.clone(); + let mut cf = coeffs.clone(); + spec.pc_step(&mut b, &mut cf, DT, &[], &cfg).unwrap(); + black_box(cf.len()) + }) + }); + } + } + group.finish(); +} + +criterion_group!(benches, bench_drug); +criterion_main!(benches); diff --git a/crates/ppvm-lindblad/benches/kossakowski.rs b/crates/ppvm-lindblad/benches/kossakowski.rs new file mode 100644 index 000000000..22724a6d0 --- /dev/null +++ b/crates/ppvm-lindblad/benches/kossakowski.rs @@ -0,0 +1,173 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Per-step cost of the Kossakowski-form dissipator vs the equivalent +//! eigenmode-jump representation, on the subwavelength superradiance chain +//! (free-space photon-mediated collective σ⁻ decay, d = 0.1 λ₀). +//! +//! Both specs generate the identical adjoint action; the difference is +//! representation cost: eigenmode jumps pay `N · (2N)²` Pauli products per +//! dissipator evaluation, Kossakowski pairs pay `4·nnz(Γ) = 4N²`. +//! +//! The benchmark grows a realistic working basis with a few capped +//! `pc_step` calls, then measures one full `pc_step` (two leakage passes + +//! predictor/corrector expm) from a cloned copy of that basis. + +use criterion::{Criterion, criterion_group, criterion_main}; +use num::Complex; +use ppvm_lindblad::{JumpInput, LindbladSpec, PcStepConfig, Word, parse_pauli_string}; +use std::f64::consts::PI; +use std::hint::black_box; + +const G0: f64 = 1.0; +const D_OVER_LAM: f64 = 0.1; +const B: usize = 4096; +const GROW_STEPS: usize = 3; +const DT: f64 = 0.01; + +/// Free-space couplings `(J, Γ)` of a chain along x with spacing `d·λ₀`, +/// circular polarization `(1, i, 0)/√2`. +fn chain_couplings(n: usize) -> (Vec>, Vec>) { + let k0 = 2.0 * PI; + let mut j = vec![vec![0.0; n]; n]; + let mut gam = vec![vec![0.0; n]; n]; + for (a, (j_row, gam_row)) in j.iter_mut().zip(gam.iter_mut()).enumerate() { + for b in 0..n { + if a == b { + gam_row[b] = G0; + continue; + } + let r = (a as f64 - b as f64).abs() * D_OVER_LAM; + let kr = k0 * r; + let e = Complex::from_polar(1.0, kr); + let pref = e / (4.0 * PI * k0 * k0 * r * r * r); + // p†·G·p with p = (1, i, 0)/√2 and r̂ = x̂: + // p†·(kr²+ikr−1)·1·p = (kr²+ikr−1); p†·r̂r̂·p = 1/2. + let g = pref + * (Complex::new(kr * kr - 1.0, kr) - Complex::new(kr * kr - 3.0, 3.0 * kr) * 0.5); + j_row[b] = -3.0 * PI * G0 / k0 * g.re; + gam_row[b] = 6.0 * PI * G0 / k0 * g.im; + } + } + (j, gam) +} + +fn pstr(n: usize, sites: &[(usize, char)]) -> String { + let mut s = vec!['I'; n]; + for &(q, c) in sites { + s[q] = c; + } + s.into_iter().collect() +} + +fn hamiltonian_terms(n: usize, j: &[Vec]) -> Vec<(String, f64)> { + let mut h = Vec::new(); + for (a, j_row) in j.iter().enumerate() { + for (b, &j_ab) in j_row.iter().enumerate().skip(a + 1) { + if j_ab.abs() > 1e-14 { + h.push((pstr(n, &[(a, 'X'), (b, 'X')]), j_ab / 2.0)); + h.push((pstr(n, &[(a, 'Y'), (b, 'Y')]), j_ab / 2.0)); + } + } + } + h +} + +fn sigma_minus(site: usize, n: usize) -> Vec<(String, Complex)> { + vec![ + (pstr(n, &[(site, 'X')]), Complex::new(0.5, 0.0)), + (pstr(n, &[(site, 'Y')]), Complex::new(0.0, -0.5)), + ] +} + +/// Eigenmode jumps `L_ν = √γ_ν Σ_j V_jν σ⁻_j` from `Γ = V diag(γ) Vᵀ`. +fn eigenmode_jumps(n: usize, gam: &[Vec]) -> Vec { + let mat = nalgebra::DMatrix::from_fn(n, n, |a, b| gam[a][b]); + let eig = nalgebra::SymmetricEigen::new(mat); + let mut jumps = Vec::new(); + for nu in 0..n { + let g = eig.eigenvalues[nu]; + if g < 1e-12 { + continue; + } + let mut lin = Vec::new(); + for j in 0..n { + let v = eig.eigenvectors[(j, nu)]; + if v.abs() > 1e-14 { + for (p, c) in sigma_minus(j, n) { + lin.push((p, c * v)); + } + } + } + jumps.push(JumpInput { + lincomb: lin, + rate: g, + }); + } + jumps +} + +/// `O = Σ_nm Γ_nm σ⁺_n σ⁻_m` as a real Pauli sum. +fn observable(n: usize, gam: &[Vec]) -> (Vec, Vec) { + let mut basis = Vec::new(); + let mut coeffs = Vec::new(); + let mut push = |s: String, c: f64| { + basis.push(parse_pauli_string(&s, n).unwrap().0); + coeffs.push(c); + }; + push(pstr(n, &[]), n as f64 * G0 / 2.0); + for (a, gam_row) in gam.iter().enumerate() { + push(pstr(n, &[(a, 'Z')]), G0 / 2.0); + for (b, &g_ab) in gam_row.iter().enumerate().skip(a + 1) { + push(pstr(n, &[(a, 'X'), (b, 'X')]), g_ab / 2.0); + push(pstr(n, &[(a, 'Y'), (b, 'Y')]), g_ab / 2.0); + } + } + (basis, coeffs) +} + +fn bench_kossakowski(c: &mut Criterion) { + let mut group = c.benchmark_group("pc_step_superradiance"); + group.sample_size(10); + for n in [10usize, 20, 30] { + let (j, gam) = chain_couplings(n); + let h = hamiltonian_terms(n, &j); + + let spec_eig = LindbladSpec::new(n, &h, &eigenmode_jumps(n, &gam)).unwrap(); + let mut spec_koss = LindbladSpec::new(n, &h, &[]).unwrap(); + let ops: Vec<_> = (0..n).map(|q| sigma_minus(q, n)).collect(); + let k: Vec>> = gam + .iter() + .map(|row| row.iter().map(|&v| Complex::new(v, 0.0)).collect()) + .collect(); + spec_koss.add_kossakowski(&ops, &k).unwrap(); + + // Grow a realistic capped working basis once (shared by both). + let cfg = PcStepConfig { + max_basis: B, + admit_basis: Some(3 * B), + ..Default::default() + }; + let (mut basis, mut coeffs) = observable(n, &gam); + for _ in 0..GROW_STEPS { + spec_koss + .pc_step(&mut basis, &mut coeffs, DT, &[], &cfg) + .unwrap(); + } + + for (label, spec) in [("eigenmode", &spec_eig), ("kossakowski", &spec_koss)] { + group.bench_function(format!("{label}_n{n}"), |bch| { + bch.iter(|| { + let mut b = basis.clone(); + let mut cf = coeffs.clone(); + spec.pc_step(&mut b, &mut cf, DT, &[], &cfg).unwrap(); + black_box(cf.len()) + }) + }); + } + } + group.finish(); +} + +criterion_group!(benches, bench_kossakowski); +criterion_main!(benches); diff --git a/crates/ppvm-lindblad/examples/drug_profile.rs b/crates/ppvm-lindblad/examples/drug_profile.rs new file mode 100644 index 000000000..5c8f4e1bf --- /dev/null +++ b/crates/ppvm-lindblad/examples/drug_profile.rs @@ -0,0 +1,217 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Fast single-config profiler for the Kossakowski-form dissipator on the +//! molecular dipolar-relaxation (ZULF drug-FID) workload — the A/B harness +//! for the 2026-07-17-drug-kossakowski autotune campaign. +//! +//! Prints per-phase `pc_step_timed` breakdown (median of N steps) for the +//! Kossakowski path only. The eigenmode representation is *not* measured here +//! (its O(N³) dense-jump blowup is the thing this path removes — see the +//! `drug_dipolar` criterion bench for the documented representation ratio). +//! +//! Usage: `cargo run --release --example drug_profile -- [N] [B] [STEPS]` + +use num::Complex; +use ppvm_lindblad::{LindbladSpec, PcStepConfig, Word, parse_pauli_string}; +use std::time::Instant; + +const N_M: usize = 5; + +/// Dipolar-pair geometry: the `(a, b)` site pairs, each pair's coupling +/// magnitude, and its unit separation vector. +type Geometry = (Vec<(usize, usize)>, Vec, Vec<[f64; 3]>); + +/// A Kossakowski dissipator as handed to `LindbladSpec::add_kossakowski`: +/// the operators `A_n` as Pauli lincombs, and the pair matrix `K`. +type KossakowskiModel = (Vec)>>, Vec>>); + +fn pstr(n: usize, sites: &[(usize, char)]) -> String { + let mut s = vec!['I'; n]; + for &(q, c) in sites { + s[q] = c; + } + s.into_iter().collect() +} + +fn hashf(mut x: u64) -> f64 { + x ^= x >> 33; + x = x.wrapping_mul(0xff51afd7ed558ccd); + x ^= x >> 33; + (x >> 11) as f64 / (1u64 << 53) as f64 +} + +fn geometry(n: usize) -> Geometry { + let pos: Vec<[f64; 3]> = (0..n) + .map(|i| { + [ + i as f64 + 0.3 * hashf(i as u64 * 3 + 1), + 0.4 * hashf(i as u64 * 3 + 2), + 0.4 * hashf(i as u64 * 3 + 3), + ] + }) + .collect(); + let (mut pairs, mut bmag, mut dir) = (Vec::new(), Vec::new(), Vec::new()); + for a in 0..n { + for b in (a + 1)..n { + let d = [ + pos[a][0] - pos[b][0], + pos[a][1] - pos[b][1], + pos[a][2] - pos[b][2], + ]; + let r = (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt(); + pairs.push((a, b)); + bmag.push(1.0 / (r * r * r)); + dir.push([d[0] / r, d[1] / r, d[2] / r]); + } + } + (pairs, bmag, dir) +} + +fn y2(u: &[f64; 3]) -> [f64; N_M] { + let (x, y, z) = (u[0], u[1], u[2]); + [ + x * y, + y * z, + (3.0 * z * z - 1.0) / 2.0, + x * z, + (x * x - y * y) / 2.0, + ] +} + +fn tensor_op(n: usize, a: usize, b: usize, mt: usize) -> Vec<(String, Complex)> { + let (i, j) = (Complex::new(0.0, 1.0), Complex::new(1.0, 0.0)); + match mt { + 0 | 4 => { + let s = if mt == 0 { i } else { -i }; + vec![ + (pstr(n, &[(a, 'X'), (b, 'X')]), j), + (pstr(n, &[(a, 'Y'), (b, 'Y')]), -j), + (pstr(n, &[(a, 'X'), (b, 'Y')]), s), + (pstr(n, &[(a, 'Y'), (b, 'X')]), s), + ] + } + 1 | 3 => { + let s = if mt == 1 { i } else { -i }; + vec![ + (pstr(n, &[(a, 'X'), (b, 'Z')]), j), + (pstr(n, &[(a, 'Y'), (b, 'Z')]), s), + (pstr(n, &[(a, 'Z'), (b, 'X')]), j), + (pstr(n, &[(a, 'Z'), (b, 'Y')]), s), + ] + } + _ => vec![ + (pstr(n, &[(a, 'X'), (b, 'X')]), j), + (pstr(n, &[(a, 'Y'), (b, 'Y')]), j), + (pstr(n, &[(a, 'Z'), (b, 'Z')]), Complex::new(2.0, 0.0)), + ], + } +} + +fn model(n: usize) -> (Vec<(String, f64)>, KossakowskiModel) { + let (pairs, bmag, dir) = geometry(n); + let p = pairs.len(); + let mut h = Vec::new(); + for (k, &(a, b)) in pairs.iter().enumerate() { + let jc = 0.1 * bmag[k]; + h.push((pstr(n, &[(a, 'X'), (b, 'X')]), jc)); + h.push((pstr(n, &[(a, 'Y'), (b, 'Y')]), jc)); + } + let c: Vec<[f64; N_M]> = (0..p) + .map(|k| { + let y = y2(&dir[k]); + std::array::from_fn(|mp| bmag[k] * y[mp]) + }) + .collect(); + let mut ops = Vec::with_capacity(N_M * p); + for mt in 0..N_M { + for &(a, b) in &pairs { + ops.push(tensor_op(n, a, b, mt)); + } + } + let m_ops = N_M * p; + let mut k = vec![vec![Complex::new(0.0, 0.0); m_ops]; m_ops]; + for mt in 0..N_M { + let off = mt * p; + for mu in 0..p { + for nu in 0..p { + let g: f64 = (0..N_M).map(|mp| c[mu][mp] * c[nu][mp]).sum(); + k[off + mu][off + nu] = Complex::new(g, 0.0); + } + } + } + (h, (ops, k)) +} + +fn observable(n: usize) -> (Vec, Vec) { + let mut basis = Vec::new(); + let mut coeffs = Vec::new(); + for a in 0..n { + basis.push(parse_pauli_string(&pstr(n, &[(a, 'X')]), n).unwrap().0); + coeffs.push(1.0); + } + (basis, coeffs) +} + +fn main() { + let args: Vec = std::env::args().collect(); + let n: usize = args.get(1).and_then(|s| s.parse().ok()).unwrap_or(20); + let b: usize = args.get(2).and_then(|s| s.parse().ok()).unwrap_or(4096); + let steps: usize = args.get(3).and_then(|s| s.parse().ok()).unwrap_or(8); + let dt = 1e-3; + + let (h, (ops, k)) = model(n); + let mut spec = LindbladSpec::new(n, &h, &[]).unwrap(); + let t0 = Instant::now(); + spec.add_kossakowski(&ops, &k).unwrap(); + let build_ms = t0.elapsed().as_secs_f64() * 1e3; + + let cfg = PcStepConfig { + max_basis: b, + admit_basis: Some(3 * b), + ..Default::default() + }; + let (mut basis, mut coeffs) = observable(n); + // Grow into a realistic capped basis (not timed). + for _ in 0..3 { + spec.pc_step(&mut basis, &mut coeffs, dt, &[], &cfg) + .unwrap(); + } + + let mut totals = Vec::new(); + let (mut l1, mut e1, mut x1, mut l2, mut e2, mut x2) = (0u64, 0u64, 0u64, 0u64, 0u64, 0u64); + for _ in 0..steps { + let mut bb = basis.clone(); + let mut cf = coeffs.clone(); + let t = spec.pc_step_timed(&mut bb, &mut cf, dt, &[], &cfg).unwrap(); + totals.push(t.total_us()); + l1 += t.leakage1_us; + e1 += t.expand1_us; + x1 += t.expm1_us; + l2 += t.leakage2_us; + e2 += t.expand2_us; + x2 += t.expm2_us; + } + totals.sort_unstable(); + let med = totals[totals.len() / 2] as f64 / 1e3; + let s = steps as f64; + println!( + "N={n} B={b} pairs={} ops={} nnz(K)={}", + n * (n - 1) / 2, + ops.len(), + N_M * (n * (n - 1) / 2) * (n * (n - 1) / 2) + ); + println!(" add_kossakowski build: {build_ms:.0} ms"); + println!(" median total/step: {med:.1} ms (over {steps} steps)"); + println!( + " phase avg (ms): leak1 {:.1} expm1 {:.1} leak2 {:.1} expm2 {:.1} expand {:.1}", + l1 as f64 / s / 1e3, + x1 as f64 / s / 1e3, + l2 as f64 / s / 1e3, + x2 as f64 / s / 1e3, + (e1 + e2) as f64 / s / 1e3, + ); + let diss = (l1 + l2) as f64; + let tot = (l1 + e1 + x1 + l2 + e2 + x2) as f64; + println!(" leakage(action) share: {:.0}%", 100.0 * diss / tot); +} From dc8c6234e049ec5efec974b3d46c1bde90e8febe Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 16 Sep 2026 09:59:07 +0200 Subject: [PATCH 3/7] feat(python): Kossakowski dissipator on Lindbladian MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `Lindbladian(..., kossakowski=(ops, K))` forwards an operator family and its pair matrix to the native spec. Mathematically identical to passing the eigenmode jumps of `K = V diag(γ) V†`, but the per-string action cost scales with the number of nonzero `K_nm` entries rather than carrying an extra factor of `M`. Python validates `K` for Hermiticity and positive-semidefiniteness before the call. The core rejects a non-Hermitian `K` too, but the Python check is not redundant: `eigvalsh` reads only one triangle, so without it a non-Hermitian `K` would pass the PSD test on its symmetric part. Co-Authored-By: Claude Opus 5 --- crates/ppvm-python-native/src/lindblad.rs | 39 +++++++++++--- ppvm-python/src/ppvm/_core.pyi | 2 + ppvm-python/src/ppvm/lindblad.py | 62 ++++++++++++++++++++++- 3 files changed, 96 insertions(+), 7 deletions(-) diff --git a/crates/ppvm-python-native/src/lindblad.rs b/crates/ppvm-python-native/src/lindblad.rs index 477e1e834..792b2e8b2 100644 --- a/crates/ppvm-python-native/src/lindblad.rs +++ b/crates/ppvm-python-native/src/lindblad.rs @@ -48,6 +48,15 @@ use crate::pauli_arr::{ check_coeffs_len, check_group_qubits, check_momentum_len, decode_basis, encode_basis, }; +/// Convert the `(pauli_string, re, im)` triple encoding used across the +/// Python boundary into a complex Pauli linear combination. +fn to_lincomb(terms: Vec<(String, f64, f64)>) -> Vec<(String, Complex)> { + terms + .into_iter() + .map(|(s, re, im)| (s, Complex::new(re, im))) + .collect() +} + /// Pack `Vec<(Word, f64)>` into the standard PyO3 return shape. fn pack_pauli_map<'py>( py: Python<'py>, @@ -72,14 +81,22 @@ impl LindbladSpec { /// `jump_lincombs[k]` is a list of `(pauli_string, real, imag)` triples /// encoding `L_k = Σ_a (re + i·im) P_a`. A length-1 jump with `im == 0` /// is routed to the Hermitian-Pauli fast path (with rate scaled by `re²`). + /// + /// `kossakowski_ops` / `kossakowski_k` optionally add a Kossakowski-form + /// dissipator `D*(O) = Σ_nm K_nm (A_n† O A_m − ½{A_n†A_m, O})`: + /// `kossakowski_ops[i]` is the Pauli lincomb of `A_i` in the same triple + /// encoding, `kossakowski_k[n][m] = (re, im)` the Hermitian pair matrix. #[new] - #[pyo3(signature = (n_qubits, h_terms, h_coeffs, jump_lincombs, jump_rates))] + #[pyo3(signature = (n_qubits, h_terms, h_coeffs, jump_lincombs, jump_rates, + kossakowski_ops = vec![], kossakowski_k = vec![]))] fn new( n_qubits: usize, h_terms: Vec, h_coeffs: Vec, jump_lincombs: Vec>, jump_rates: Vec, + kossakowski_ops: Vec>, + kossakowski_k: Vec>, ) -> PyResult { if h_terms.len() != h_coeffs.len() { return Err(PyValueError::new_err( @@ -96,14 +113,24 @@ impl LindbladSpec { .into_iter() .zip(jump_rates) .map(|(lincomb, rate)| JumpInput { - lincomb: lincomb - .into_iter() - .map(|(s, re, im)| (s, Complex::new(re, im))) - .collect(), + lincomb: to_lincomb(lincomb), rate, }) .collect(); - let inner = CoreSpec::new(n_qubits, &h, &jumps).map_err(map_err)?; + let mut inner = CoreSpec::new(n_qubits, &h, &jumps).map_err(map_err)?; + if !kossakowski_ops.is_empty() || !kossakowski_k.is_empty() { + let ops: Vec)>> = + kossakowski_ops.into_iter().map(to_lincomb).collect(); + let k: Vec>> = kossakowski_k + .into_iter() + .map(|row| { + row.into_iter() + .map(|(re, im)| Complex::new(re, im)) + .collect() + }) + .collect(); + inner.add_kossakowski(&ops, &k).map_err(map_err)?; + } Ok(Self { inner }) } diff --git a/ppvm-python/src/ppvm/_core.pyi b/ppvm-python/src/ppvm/_core.pyi index e9877b7f8..86819da03 100644 --- a/ppvm-python/src/ppvm/_core.pyi +++ b/ppvm-python/src/ppvm/_core.pyi @@ -365,6 +365,8 @@ class LindbladSpec: h_coeffs: list[float], jump_lincombs: list[list[tuple[str, float, float]]], jump_rates: list[float], + kossakowski_ops: list[list[tuple[str, float, float]]] = ..., + kossakowski_k: list[list[tuple[float, float]]] = ..., ) -> None: ... @property def n_qubits(self) -> int: ... diff --git a/ppvm-python/src/ppvm/lindblad.py b/ppvm-python/src/ppvm/lindblad.py index a71dce745..8dc3f9455 100644 --- a/ppvm-python/src/ppvm/lindblad.py +++ b/ppvm-python/src/ppvm/lindblad.py @@ -156,6 +156,21 @@ class Lindbladian: Pauli linear combination such as `sigma_plus` or `sigma_minus`. ``rate`` is the non-negative GKSL rate ``γ_k``. + kossakowski: + Optional ``(ops, K)`` pair adding a Kossakowski-form dissipator + ``D*(O) = Σ_nm K_nm (A_n† O A_m − ½{A_n†A_m, O})``. ``ops`` is a + length-``M`` sequence of Pauli lincombs (same format as a + ``jump_op``, e.g. ``[sigma_minus(j, n) for j in range(n)]``); + ``K`` is an ``(M, M)`` Hermitian positive-semidefinite array + (e.g. the collective-decay pair matrix ``Γ_nm``). Mathematically + identical to passing the eigenmode jumps + ``L_ν = √γ_ν Σ_j V*_jν A_j`` of ``K = V diag(γ) V†``, but the + per-string action cost scales with the number of nonzero + ``K_nm`` entries instead of carrying an extra factor of ``M``. + May coexist with ``jump_terms`` (both contribute). Raises + ``ValueError`` if ``K`` is not Hermitian or has an eigenvalue + below ``−tol·‖K‖`` (a non-PSD ``K`` is not a valid GKSL + generator; no silent clipping). Examples -------- @@ -167,13 +182,23 @@ class Lindbladian: >>> jumps = [(sigma_minus(0, 2), 0.5)] >>> Lindbladian(2, [("XX", 1.0)], jumps) + + Collective decay from a pair matrix: + + >>> ops = [sigma_minus(j, 2) for j in range(2)] + >>> gamma = [[1.0, 0.6], [0.6, 1.0]] + >>> Lindbladian(2, [("XX", 1.0)], kossakowski=(ops, gamma)) """ + #: Relative tolerance for the Hermiticity and PSD checks on ``K``. + _K_TOL = 1e-10 + def __init__( self, n_qubits: int, h_terms: Iterable[tuple[str, float]], jump_terms: Iterable[tuple[str | PauliLincomb, float]] = (), + kossakowski: tuple[Sequence[str | PauliLincomb], npt.ArrayLike] | None = None, ): self.n_qubits = int(n_qubits) h_strs: list[str] = [] @@ -186,7 +211,42 @@ def __init__( for jump_op, rate in jump_terms: j_lincombs.append(_normalize_jump(jump_op)) j_rates.append(float(rate)) - self._spec = _LindbladSpec(self.n_qubits, h_strs, h_coeffs, j_lincombs, j_rates) + k_ops, k_mat = self._normalize_kossakowski(kossakowski) + self._spec = _LindbladSpec( + self.n_qubits, h_strs, h_coeffs, j_lincombs, j_rates, k_ops, k_mat + ) + + @classmethod + def _normalize_kossakowski( + cls, kossakowski: tuple[Sequence[str | PauliLincomb], npt.ArrayLike] | None + ) -> tuple[list[list[tuple[str, float, float]]], list[list[tuple[float, float]]]]: + """Validate and flatten the ``(ops, K)`` pair for the native constructor. + + The core rejects a non-Hermitian ``K`` too, but the check is repeated + here because ``eigvalsh`` reads only one triangle: without it a + non-Hermitian ``K`` would pass the PSD test on its symmetric part. + """ + if kossakowski is None: + return [], [] + ops, k = kossakowski + k_ops = [_normalize_jump(op) for op in ops] + k_arr = np.asarray(k, dtype=np.complex128) + if k_arr.shape != (len(k_ops), len(k_ops)): + raise ValueError( + f"kossakowski K has shape {k_arr.shape} but ops has length {len(k_ops)}" + ) + scale = max(float(np.abs(k_arr).max(initial=0.0)), 1.0) + if float(np.abs(k_arr - k_arr.conj().T).max(initial=0.0)) > cls._K_TOL * scale: + raise ValueError("kossakowski K is not Hermitian") + evals = np.linalg.eigvalsh(k_arr) + if float(evals.min(initial=0.0)) < -cls._K_TOL * scale: + raise ValueError( + f"kossakowski K is not positive semidefinite " + f"(min eigenvalue {evals.min():.3e}); a non-PSD pair matrix " + f"is not a valid GKSL generator" + ) + k_mat = [[(float(v.real), float(v.imag)) for v in row] for row in k_arr] + return k_ops, k_mat @property def num_h_terms(self) -> int: From b52d22cc3120fdb067ffb60a5a4d6f9c0cce359e Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Wed, 16 Sep 2026 09:59:24 +0200 Subject: [PATCH 4/7] test(lindblad): Kossakowski dissipator against dense references MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `test_kossakowski` checks the pair-matrix action against a dense superoperator and against the equivalent eigenmode-jump decomposition; `test_orbit_dissipative` covers the orbit-representative path on a ring of emitters, including the non-unital Z → I flow across orbits of different sizes. The orbit tests reconstruct a translation-invariant sum as `Σ |orbit(rep)| · c_rep`. Rep coefficients are member coefficients on every orbit since the stabilized-orbit fix in d8313e0e, so the previous uniform `|G| ·` factor would over-count orbits with a non-trivial stabilizer. Co-Authored-By: Claude Opus 5 --- .../test/lindblad/_dissipative_refs.py | 195 +++++++++++++++ ppvm-python/test/lindblad/test_kossakowski.py | 170 +++++++++++++ .../test/lindblad/test_orbit_dissipative.py | 232 ++++++++++++++++++ 3 files changed, 597 insertions(+) create mode 100644 ppvm-python/test/lindblad/_dissipative_refs.py create mode 100644 ppvm-python/test/lindblad/test_kossakowski.py create mode 100644 ppvm-python/test/lindblad/test_orbit_dissipative.py diff --git a/ppvm-python/test/lindblad/_dissipative_refs.py b/ppvm-python/test/lindblad/_dissipative_refs.py new file mode 100644 index 000000000..fcf46c528 --- /dev/null +++ b/ppvm-python/test/lindblad/_dissipative_refs.py @@ -0,0 +1,195 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Shared references for the collective-decay (Kossakowski) tests. + +Vendored from the superradiant-burst study (free-space photon-mediated +Lindbladian, arXiv:2309.11376 Eqs. 6-8): emitter couplings from the dyadic +Green's function, the ppvm model builder, and the exact Lindblad reference +via the excitation-number cascade. `exact_rate_sector` is geometry-agnostic +(it takes the J and Gamma matrices), so the same reference covers chains +and rings. + +Conventions: excited = |0> (Z = +1), sigma^- = (X - i Y)/2; +H = sum_{n(t), the photon emission rate. +""" + +import itertools + +import numpy as np + +LAM = 1.0 +K0 = 2 * np.pi / LAM +G0 = 1.0 +POL = np.array([1.0, 1.0j, 0.0]) / np.sqrt(2.0) + + +def greens(r_vec): + """Free-space dyadic Green's function G(r, omega0).""" + r = np.linalg.norm(r_vec) + rh = np.outer(r_vec, r_vec) / r**2 + kr = K0 * r + pref = np.exp(1j * kr) / (4 * np.pi * K0**2 * r**3) + return pref * ((kr**2 + 1j * kr - 1) * np.eye(3) - (kr**2 + 3j * kr - 3) * rh) + + +def couplings(pos): + """J_nm, Gamma_nm from emitter positions; J_nn = 0, Gamma_nn = Gamma_0.""" + n = len(pos) + J = np.zeros((n, n)) + Gam = np.zeros((n, n)) + for a in range(n): + for b in range(n): + if a == b: + Gam[a, b] = G0 + continue + g = POL.conj() @ greens(pos[a] - pos[b]) @ POL + J[a, b] = -3 * np.pi * G0 / K0 * g.real + Gam[a, b] = 6 * np.pi * G0 / K0 * g.imag + return J, Gam + + +def chain_positions(n, d): + return [np.array([j * d, 0.0, 0.0]) for j in range(n)] + + +def ring_positions(n, d): + """n emitters on a circle with nearest-neighbour arc spacing ~d. + + Chord-based radius so that adjacent emitters sit exactly d apart. + The resulting J/Gamma matrices are circulant (exact C_n symmetry). + """ + radius = d / (2 * np.sin(np.pi / n)) + return [ + np.array([radius * np.cos(2 * np.pi * j / n), radius * np.sin(2 * np.pi * j / n), 0.0]) + for j in range(n) + ] + + +def pstr(n, **sites): + s = ["I"] * n + for k, v in sites.items(): + s[int(k[1:])] = v + return "".join(s) + + +def hamiltonian_terms(n, J): + """H = sum_{a 1e-14: + h_terms.append((pstr(n, **{f"q{a}": "X", f"q{b}": "X"}), J[a, b] / 2)) + h_terms.append((pstr(n, **{f"q{a}": "Y", f"q{b}": "Y"}), J[a, b] / 2)) + return h_terms + + +def eigenmode_jumps(ops, K): + """Jump list equivalent to the Kossakowski pair (ops, K). + + K = V diag(g) V^dagger; L_nu = sqrt(g_nu) sum_j conj(V_j_nu) A_j with + rate g_nu (the sqrt is absorbed into the rate as g_nu). + """ + g_nu, V = np.linalg.eigh(np.asarray(K, dtype=complex)) + jumps = [] + for nu in range(len(ops)): + if g_nu[nu] < 1e-12: + continue + lin = [] + for j, op in enumerate(ops): + v = np.conj(V[j, nu]) + if abs(v) > 1e-14: + for p, c in op: + lin.append((p, complex(c) * v)) + jumps.append((lin, float(g_nu[nu]))) + return jumps + + +def rate_observable(n, Gam): + """O = sum_nm Gamma_nm s+_n s-_m as {pauli_string: real_coeff}.""" + obs = {pstr(n): n * G0 / 2} + for a in range(n): + obs[pstr(n, **{f"q{a}": "Z"})] = G0 / 2 + for b in range(a + 1, n): + obs[pstr(n, **{f"q{a}": "X", f"q{b}": "X"})] = Gam[a, b] / 2 + obs[pstr(n, **{f"q{a}": "Y", f"q{b}": "Y"})] = Gam[a, b] / 2 + return obs + + +def exact_rate_sector(n, J, Gam, T_run, dt_out, dt_inner=2e-3): + """Exact Lindblad R(t) via the excitation-number cascade (dense blocks). + + From the fully inverted state, H conserves the excitation number M and + every jump lowers M symmetrically on both sides of rho, so rho(t) is a + direct sum of C(n, M)-sized blocks, evolved here with RK4. + """ + sectors = [] + for M in range(n, -1, -1): + confs = [frozenset(c) for c in itertools.combinations(range(n), M)] + sectors.append({c: i for i, c in enumerate(confs)}) + + def block(i, mat): + idx = sectors[i] + out = np.zeros((len(idx), len(idx)), dtype=complex) + for conf, b in idx.items(): + for m in conf: + for nn in range(n): + if nn == m: + out[b, b] += mat[m, m] + elif nn not in conf: + out[idx[(conf - {m}) | {nn}], b] += mat[nn, m] + return out + + H = [block(i, J) for i in range(n + 1)] + A = [block(i, Gam) for i in range(n + 1)] + add = [] + for i in range(1, n + 1): + idx, idx_up = sectors[i], sectors[i - 1] + amap = np.full((n, len(idx)), -1, dtype=np.int64) + for conf, b in idx.items(): + for m in range(n): + if m not in conf: + amap[m, b] = idx_up[conf | {m}] + add.append(amap) + + def feed(i, rho_up): + amap = add[i - 1] + out = np.zeros((len(sectors[i]), len(sectors[i])), dtype=complex) + padded = np.pad(rho_up, ((0, 1), (0, 1))) + for m in range(n): + for nn in range(n): + if abs(Gam[nn, m]) < 1e-14: + continue + sub = padded[amap[m][:, None], amap[nn][None, :]] + mask = (amap[m][:, None] >= 0) & (amap[nn][None, :] >= 0) + out += Gam[nn, m] * np.where(mask, sub, 0.0) + return out + + def rhs(blocks): + out = [] + for i, r in enumerate(blocks): + d = -1j * (H[i] @ r - r @ H[i]) - 0.5 * (A[i] @ r + r @ A[i]) + if i > 0: + d += feed(i, blocks[i - 1]) + out.append(d) + return out + + blocks = [np.zeros((len(s), len(s)), dtype=complex) for s in sectors] + blocks[0][0, 0] = 1.0 + n_out = round(T_run / dt_out) + sub = max(1, int(np.ceil(dt_out / dt_inner))) + h = dt_out / sub + R = np.zeros(n_out + 1) + R[0] = sum(np.trace(r @ a).real for r, a in zip(blocks, A)) + for k in range(n_out): + for _ in range(sub): + k1 = rhs(blocks) + k2 = rhs([r + h / 2 * d for r, d in zip(blocks, k1)]) + k3 = rhs([r + h / 2 * d for r, d in zip(blocks, k2)]) + k4 = rhs([r + h * d for r, d in zip(blocks, k3)]) + blocks = [ + r + h / 6 * (a + 2 * b + 2 * c + e) for r, a, b, c, e in zip(blocks, k1, k2, k3, k4) + ] + R[k + 1] = sum(np.trace(r @ a).real for r, a in zip(blocks, A)) + return np.arange(n_out + 1) * dt_out, R diff --git a/ppvm-python/test/lindblad/test_kossakowski.py b/ppvm-python/test/lindblad/test_kossakowski.py new file mode 100644 index 000000000..96570f353 --- /dev/null +++ b/ppvm-python/test/lindblad/test_kossakowski.py @@ -0,0 +1,170 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Kossakowski-form dissipator: exact equivalence with the eigenmode-jump +representation, physics regression against the exact excitation-cascade +reference, and leakage coverage. +""" + +import pathlib + +import numpy as np +import pytest + +from ppvm import Lindbladian +from ppvm.lindblad import _basis_to_codes, _codes_to_basis, sigma_minus + +from ._dissipative_refs import ( + chain_positions, + couplings, + eigenmode_jumps, + exact_rate_sector, + hamiltonian_terms, + rate_observable, +) + +BIG = 10_000_000 # uncapped max_basis + +FIG11_H5 = pathlib.Path( + "/Users/alexschuckert/dev/26_ppvm/CTPP Figures/fig11_superradiant_burst/data.h5" +) + + +def run_steps(lind, obs, n, dt, steps): + """Evolve {string: coeff} a few uncapped pc steps; return final dict.""" + strings = list(obs) + basis = _basis_to_codes(strings, n) + coeff = np.array([obs[s] for s in strings], dtype=np.float64) + for _ in range(steps): + basis, coeff = lind.pc_step_arr(basis, coeff, dt, max_basis=BIG, drop_tol=0.0) + return dict(zip(_codes_to_basis(basis), coeff)) + + +def random_psd(rng, m, complex_k=False): + a = rng.standard_normal((m, m)) + if complex_k: + a = a + 1j * rng.standard_normal((m, m)) + return a @ a.conj().T / m + + +def random_lincomb(rng, n): + """A random 2-term Pauli lincomb with complex coefficients.""" + ops = "IXYZ" + out = [] + for _ in range(2): + s = "".join(rng.choice(list(ops)) for _ in range(n)) + if s == "I" * n: + s = "X" + s[1:] + c = complex(rng.standard_normal(), rng.standard_normal()) + out.append((s, c)) + return out + + +@pytest.mark.parametrize("complex_k", [False, True]) +def test_random_model_equivalence(complex_k): + """Eigenmode jumps built from K and kossakowski=(ops, K) produce the + same evolution to near machine precision (same dt, uncapped basis).""" + rng = np.random.default_rng(7 if complex_k else 3) + n, dt, steps = 4, 0.02, 3 + # ops: all single-site sigma^- plus one random 2-term lincomb + ops = [sigma_minus(j, n) for j in range(n)] + [random_lincomb(rng, n)] + K = random_psd(rng, len(ops), complex_k) + h_terms = [("XX" + "I" * (n - 2), 0.9), ("I" + "ZZ" + "I" * (n - 3), -0.4)] + obs = {"Z" + "I" * (n - 1): 1.0, "IXY" + "I" * (n - 3): 0.3} + + out_k = run_steps(Lindbladian(n, h_terms, kossakowski=(ops, K)), obs, n, dt, steps) + out_j = run_steps(Lindbladian(n, h_terms, eigenmode_jumps(ops, K)), obs, n, dt, steps) + + assert set(out_k) == set(out_j) + max_dev = max(abs(out_k[s] - out_j[s]) for s in out_k) + assert max_dev < 1e-12, f"representations diverged: max |dc| = {max_dev:.2e}" + + +def test_kossakowski_coexists_with_jump_terms(): + """kossakowski= and jump_terms may both contribute.""" + n = 2 + ops = [sigma_minus(j, n) for j in range(n)] + K = [[1.0, 0.5], [0.5, 1.0]] + both = Lindbladian(n, [], [("ZI", 0.3)], kossakowski=(ops, K)) + only_k = Lindbladian(n, [], kossakowski=(ops, K)) + only_j = Lindbladian(n, [], [("ZI", 0.3)]) + a_both = both.action("XI") + a_sum = {} + for d in (only_k.action("XI"), only_j.action("XI")): + for s, c in d.items(): + a_sum[s] = a_sum.get(s, 0.0) + c + for s in set(a_both) | set(a_sum): + assert abs(a_both.get(s, 0.0) - a_sum.get(s, 0.0)) < 1e-13 + + +def superradiance_chain(n, d_over_lam=0.1): + J, Gam = couplings(chain_positions(n, d_over_lam)) + ops = [sigma_minus(j, n) for j in range(n)] + return J, Gam, hamiltonian_terms(n, J), ops, rate_observable(n, Gam) + + +def rate_trace(lind, obs, n, dt, steps): + """R(t) on the fully inverted state = sum of {I,Z}-string coefficients.""" + strings = list(obs) + basis = _basis_to_codes(strings, n) + coeff = np.array([obs[s] for s in strings], dtype=np.float64) + R = np.zeros(steps + 1) + for k in range(steps + 1): + iz = np.all((basis == 0) | (basis == 2), axis=1) # codes: I=0, Z=2 + R[k] = coeff[iz].sum() + if k == steps: + break + basis, coeff = lind.pc_step_arr(basis, coeff, dt, max_basis=BIG, drop_tol=0.0) + return R + + +def test_superradiance_physics_regression(): + """N=6 subwavelength chain, full basis, T=1: the Kossakowski path matches + the eigenmode path to ~1e-12 and the exact cascade reference to < 1e-4.""" + n, dt, T = 6, 0.01, 1.0 + steps = round(T / dt) + J, Gam, h_terms, ops, obs = superradiance_chain(n) + + R_k = rate_trace(Lindbladian(n, h_terms, kossakowski=(ops, Gam)), obs, n, dt, steps) + R_j = rate_trace(Lindbladian(n, h_terms, eigenmode_jumps(ops, Gam)), obs, n, dt, steps) + assert np.abs(R_k - R_j).max() < 1e-11, ( + f"kossakowski vs eigenmode R(t): {np.abs(R_k - R_j).max():.2e}" + ) + + _, R_exact = exact_rate_sector(n, J, Gam, T_run=T, dt_out=dt) + err = np.abs(R_k - R_exact).max() + assert err < 1e-4, f"kossakowski vs exact cascade: max |dR| = {err:.2e}" + + +@pytest.mark.skipif(not FIG11_H5.exists(), reason="fig11 data.h5 not present") +def test_superradiance_vs_fig11_reference(): + """Cross-check R(t) against the stored exact reference of the + superradiant-burst study (same model, N=6, first 1/Gamma_0).""" + h5py = pytest.importorskip("h5py") + n, dt, T = 6, 0.01, 1.0 + steps = round(T / dt) + _, Gam, h_terms, ops, obs = superradiance_chain(n) + R_k = rate_trace(Lindbladian(n, h_terms, kossakowski=(ops, Gam)), obs, n, dt, steps) + with h5py.File(FIG11_H5, "r") as h5: + R_ref = h5["n6/exact"][: steps + 1] + err_ref = np.abs(R_k - R_ref).max() + assert err_ref < 1e-4, f"vs fig11 data.h5 n6/exact: {err_ref:.2e}" + + +def test_leakage_covers_kossakowski_terms(): + """Leakage of a Z string under a pure-Kossakowski dissipator is nonzero + and identical to the eigenmode-jump leakage (admission sees the same + physics).""" + n = 3 + ops = [sigma_minus(j, n) for j in range(n)] + _, Gam = couplings(chain_positions(n, 0.1)) + lk = Lindbladian(n, [], kossakowski=(ops, Gam)) + lj = Lindbladian(n, [], eigenmode_jumps(ops, Gam)) + + basis = ["ZII"] + coeffs = np.array([1.0]) + leak_k = lk.leakage(basis, coeffs) + leak_j = lj.leakage(basis, coeffs) + assert leak_k, "pure-Kossakowski dissipator produced empty leakage" + assert set(leak_k) == set(leak_j) + for s in leak_k: + assert abs(leak_k[s] - leak_j[s]) < 1e-12 diff --git a/ppvm-python/test/lindblad/test_orbit_dissipative.py b/ppvm-python/test/lindblad/test_orbit_dissipative.py new file mode 100644 index 000000000..d58aa6b2d --- /dev/null +++ b/ppvm-python/test/lindblad/test_orbit_dissipative.py @@ -0,0 +1,232 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Dissipative generators on the orbit-rep path. + +`pc_step_orbit_rep` evolves canonical translation-orbit representatives +with complex coefficients under the convention set by the momentum +projector `canonicalize_basis_arr_complex`: + + c_rep = coeff of the representative word itself, + +on every orbit, whether or not it has a non-trivial stabilizer. + +Under this convention the phase-aware action is exact for ANY equivariant +generator, including jump and Kossakowski dissipators: transitions between +orbits of different sizes (e.g. the non-unital Z → I flow of σ⁻ decay, +where I is stabilized by the whole group) carry the correct weight. An +orbit contributes ``|orbit| · c_rep`` to a translation-invariant {I,Z} +sum, so the emission rate in rep space is ``R = Σ |orbit(rep)| · c_rep`` +over {I,Z} reps — see `cyclic_orbit_size` below. + +Model: ring of N emitters (positions on a circle → circulant J and Γ, +exact C_N symmetry), collective σ⁻ decay, momentum sector k = 0. +""" + +import numpy as np +import pytest + +from ppvm import Lindbladian +from ppvm._core import TranslationGroup, canonicalize_basis_arr_complex +from ppvm.lindblad import _basis_to_codes, _codes_to_basis, sigma_minus + +from ._dissipative_refs import ( + couplings, + eigenmode_jumps, + exact_rate_sector, + hamiltonian_terms, + rate_observable, + ring_positions, +) + +BIG = 10_000_000 + + +def ring_model(n, d_over_lam=0.1): + J, Gam = couplings(ring_positions(n, d_over_lam)) + assert np.allclose(Gam, np.roll(np.roll(Gam, 1, 0), 1, 1)), "Gamma not circulant" + assert np.linalg.eigvalsh(Gam).min() > -1e-12, "Gamma not PSD" + ops = [sigma_minus(j, n) for j in range(n)] + return J, Gam, hamiltonian_terms(n, J), ops, rate_observable(n, Gam) + + +def to_rep(basis, coeff, group, mom): + b, c = canonicalize_basis_arr_complex(basis, np.asarray(coeff, dtype=np.complex128), group, mom) + return dict(zip(_codes_to_basis(b), c)) + + +@pytest.mark.parametrize("representation", ["kossakowski", "eigenmode"]) +def test_orbit_matches_full_basis(representation): + """N=6 ring, uncapped: orbit-rep evolution equals the full-basis + real-space evolution projected to rep space, coefficient by + coefficient.""" + n, dt, steps = 6, 0.02, 4 + _, Gam, h_terms, ops, obs = ring_model(n) + if representation == "kossakowski": + lind = Lindbladian(n, h_terms, kossakowski=(ops, Gam)) + else: + lind = Lindbladian(n, h_terms, eigenmode_jumps(ops, Gam)) + group = TranslationGroup.chain_1d(n) + mom = np.array([0], dtype=np.int32) + + strings = list(obs) + basis0 = _basis_to_codes(strings, n) + coeff0 = np.array([obs[s] for s in strings]) + + b, c = basis0.copy(), coeff0.copy() + for _ in range(steps): + b, c = lind.pc_step_arr(b, c, dt, max_basis=BIG, drop_tol=0.0) + full = to_rep(b, c, group, mom) + + br, cr = canonicalize_basis_arr_complex(basis0, coeff0.astype(np.complex128), group, mom) + for _ in range(steps): + br, cr = lind.pc_step_orbit_rep( + br, cr, dt, max_basis=BIG, group=group, momentum=mom, drop_tol=0.0 + ) + orbit = dict(zip(_codes_to_basis(br), cr)) + + assert set(full) == set(orbit) + max_dev = max(abs(full[s] - orbit[s]) for s in full) + assert max_dev < 1e-12, f"orbit vs full-basis: max |dc| = {max_dev:.2e}" + + +def cyclic_orbit_size(word: str) -> int: + """Number of distinct cyclic rotations of `word` — its orbit size under + `TranslationGroup.chain_1d`. Words fixed by a non-trivial shift (e.g. + the identity, or `ZIZI`) have an orbit smaller than `|G|`.""" + return len({word[i:] + word[:i] for i in range(len(word))}) + + +def orbit_rate_trace(lind, obs, n, dt, steps, group, mom, max_basis=BIG, admit=None): + """R(t) from the orbit-rep evolution: the {I,Z} sum over all real-space + words, reassembled as `Σ_reps |orbit(rep)| · c_rep`.""" + strings = list(obs) + basis0 = _basis_to_codes(strings, n) + coeff0 = np.array([obs[s] for s in strings]) + br, cr = canonicalize_basis_arr_complex(basis0, coeff0.astype(np.complex128), group, mom) + R = np.zeros(steps + 1) + peak = 0 + for k in range(steps + 1): + iz = np.all((br == 0) | (br == 2), axis=1) + sizes = np.array([cyclic_orbit_size(w) for w in _codes_to_basis(br[iz])]) + R[k] = (sizes * cr[iz]).sum().real + peak = max(peak, len(cr)) + if k == steps: + break + br, cr = lind.pc_step_orbit_rep( + br, + cr, + dt, + max_basis=max_basis, + group=group, + momentum=mom, + drop_tol=0.0, + admit_basis=admit, + ) + return R, peak + + +def test_orbit_rate_vs_exact_cascade(): + """N=6 ring, T=1, full rep basis: R(t) traced down from the orbit-rep + evolution matches the excitation-cascade ED to < 1e-4.""" + n, dt, T = 6, 0.01, 1.0 + steps = round(T / dt) + J, Gam, h_terms, ops, obs = ring_model(n) + lind = Lindbladian(n, h_terms, kossakowski=(ops, Gam)) + group = TranslationGroup.chain_1d(n) + mom = np.array([0], dtype=np.int32) + + R_orbit, _ = orbit_rate_trace(lind, obs, n, dt, steps, group, mom) + _, R_exact = exact_rate_sector(n, J, Gam, T_run=T, dt_out=dt) + err = np.abs(R_orbit - R_exact).max() + assert err < 1e-4, f"orbit-rep R(t) vs exact cascade: max |dR| = {err:.2e}" + + +def test_orbit_truncated_sanity(): + """N=10 ring, genuinely truncated: no NaNs or blowup, and the error + against the (matched-capacity) real-space run improves monotonically + with the rep budget.""" + n, dt, steps = 10, 0.02, 20 + _, Gam, h_terms, ops, obs = ring_model(n) + lind = Lindbladian(n, h_terms, kossakowski=(ops, Gam)) + group = TranslationGroup.chain_1d(n) + mom = np.array([0], dtype=np.int32) + + # Real-space reference at matched effective capacity B_full = n * B_reps. + b = _basis_to_codes(list(obs), n) + c = np.array([obs[s] for s in obs]) + R_full = np.zeros(steps + 1) + for k in range(steps + 1): + iz = np.all((b == 0) | (b == 2), axis=1) + R_full[k] = c[iz].sum() + if k == steps: + break + b, c = lind.pc_step_arr( + b, c, dt, max_basis=n * 2048, drop_tol=0.0, admit_basis=3 * n * 2048 + ) + + errs = {} + for b_reps in (512, 2048): + R, peak = orbit_rate_trace( + lind, + obs, + n, + dt, + steps, + group, + mom, + max_basis=b_reps, + admit=3 * b_reps, + ) + assert np.all(np.isfinite(R)), f"non-finite R(t) at B_reps={b_reps}" + assert np.abs(R).max() < 5 * n, f"R(t) blowup at B_reps={b_reps}" + assert peak <= 3 * b_reps, "admission bound violated" + errs[b_reps] = np.abs(R - R_full).max() + + assert errs[2048] <= errs[512], ( + f"no improvement with rep budget: err(2048)={errs[2048]:.3e} > err(512)={errs[512]:.3e}" + ) + # At matched capacity the two representations should agree closely. + assert errs[2048] < 0.05 * np.abs(R_full).max(), ( + f"orbit-rep tracks real-space poorly: {errs[2048]:.3e}" + ) + + +def test_identity_bookkeeping_closed_form(): + """Uniform single-site decay (K = Γ0·1): from O = Σ_j Z_j the exact + solution is coeff_{Z_j}(t) = e^{-Γ0 t} per site and + coeff_I(t) = -n(1 - e^{-Γ0 t}) (each site pours into the identity). + Rep-space coefficients are member coefficients on every orbit, + stabilized or not, so c_Z = e^{-Γ0 t} and c_I = coeff_I — the identity + orbit has a single member and carries its full weight. This pins the + non-unital bookkeeping across orbits of different sizes.""" + n, dt, steps = 6, 0.01, 40 + ops = [sigma_minus(j, n) for j in range(n)] + lind = Lindbladian(n, [], kossakowski=(ops, np.eye(n))) + group = TranslationGroup.chain_1d(n) + mom = np.array([0], dtype=np.int32) + + # O = Σ_j Z_j → one Z rep with c = 1 (k = 0 eigenstate); + # canonicalize_first rewrites the row to the lex-min representative. + strings = ["Z" + "I" * (n - 1)] + br = _basis_to_codes(strings, n) + cr = np.array([1.0 + 0.0j]) + for _ in range(steps): + br, cr = lind.pc_step_orbit_rep( + br, + cr, + dt, + max_basis=BIG, + group=group, + momentum=mom, + drop_tol=0.0, + canonicalize_first=True, + ) + out = dict(zip(_codes_to_basis(br), cr)) + + t = steps * dt + (z_rep,) = [s for s in out if s.count("Z") == 1 and set(s) <= {"I", "Z"}] + c_z = out[z_rep] + c_i = out["I" * n] + assert abs(c_z - np.exp(-t)) < 1e-9, f"c_Z = {c_z} vs {np.exp(-t)}" + expected_i = -n * (1 - np.exp(-t)) + assert abs(c_i - expected_i) < 1e-9, f"c_I = {c_i} vs coeff_I = {expected_i}" From 235caf15196e408d30e72f5a2ff7b75d07dea40b Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Mon, 28 Sep 2026 13:42:11 +0200 Subject: [PATCH 5/7] feat(lindblad): variable Pauli-word width, lifting the 128-qubit cap MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `LindbladSpec` was hard-capped at 128 qubits by the compile-time `W_CHUNKS`, so e.g. a 130-qubit Lindbladian failed with ValueError: LindbladSpec supports n_qubits ≤ 128; got 130 which `wide-pauli-words` (4c079995) had already fixed; the fix did not make it into the split. Port it onto the current spec/basis/step layout, including the Kossakowski dissipator added in this PR. - `ppvm-lindblad` is const-generic in the chunk count: `Word`, `LindbladSpec`, and every function touching words. Both default to `C = W_CHUNKS` (128 qubits), so existing Rust callers are unchanged. `max_qubits::()`, `chunks_for(n)`, `WIDTHS` (128/256/512-qubit chunk counts, in `Chunk` units so the wasm32 u32 layout keeps working) and `MAX_SUPPORTED_QUBITS = 512` are exported. `Error::TooManyQubits` now carries the width's capacity. - PyO3: `LindbladSpec` holds an `AnySpec { W128, W256, W512 }` picked at construction and dispatches with `with_spec!`; the symmetry bindings use `with_width!`. Registers of ≤128 qubits monomorphize to exactly the old layout. More than 512 qubits is a `ValueError`. - The orbit path's masked-shift generators cover up to 1024 bits, so all three widths take the fast canonicalizer. Tests: `tests/word_width.rs` — capacities, the 130/512-qubit cases, the width error, bit-identical real-path results for the same problem padded into 128/256/512-qubit words (Kossakowski and orbit-rep to rounding: their accumulation order follows word hashes, which depend on the width), and an orbit-rep step at 130 qubits. Python: `test_word_width.py` builds and steps 130–512-qubit Lindbladians, rejects 513, and runs the momentum-orbit step across the 128-qubit boundary. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01TmnGNkxQqpmFmGszAV7dm5 --- crates/ppvm-lindblad/benches/drug_dipolar.rs | 12 +- crates/ppvm-lindblad/benches/kossakowski.rs | 8 +- crates/ppvm-lindblad/examples/drug_profile.rs | 10 +- crates/ppvm-lindblad/src/algebra.rs | 27 +- crates/ppvm-lindblad/src/basis.rs | 56 +-- crates/ppvm-lindblad/src/error.rs | 10 +- crates/ppvm-lindblad/src/kossakowski.rs | 51 +- crates/ppvm-lindblad/src/lib.rs | 12 +- crates/ppvm-lindblad/src/mf_expm.rs | 32 +- crates/ppvm-lindblad/src/sector.rs | 9 +- crates/ppvm-lindblad/src/spec.rs | 48 +- crates/ppvm-lindblad/src/step.rs | 24 +- crates/ppvm-lindblad/src/tests.rs | 24 +- crates/ppvm-lindblad/src/truncate.rs | 33 +- crates/ppvm-lindblad/src/word.rs | 88 +++- crates/ppvm-lindblad/tests/word_width.rs | 274 +++++++++++ crates/ppvm-python-native/src/lindblad.rs | 435 +++++++++++------- crates/ppvm-python-native/src/pauli_arr.rs | 39 +- crates/ppvm-python-native/src/symmetry.rs | 75 +-- ppvm-python/src/ppvm/lindblad.py | 5 +- ppvm-python/test/lindblad/test_word_width.py | 86 ++++ 21 files changed, 951 insertions(+), 407 deletions(-) create mode 100644 crates/ppvm-lindblad/tests/word_width.rs create mode 100644 ppvm-python/test/lindblad/test_word_width.py diff --git a/crates/ppvm-lindblad/benches/drug_dipolar.rs b/crates/ppvm-lindblad/benches/drug_dipolar.rs index 3e1281e62..ba17caa65 100644 --- a/crates/ppvm-lindblad/benches/drug_dipolar.rs +++ b/crates/ppvm-lindblad/benches/drug_dipolar.rs @@ -18,7 +18,7 @@ use criterion::{Criterion, criterion_group, criterion_main}; use num::Complex; -use ppvm_lindblad::{JumpInput, LindbladSpec, PcStepConfig, Word, parse_pauli_string}; +use ppvm_lindblad::{JumpInput, LindbladSpec, PcStepConfig, W_CHUNKS, Word, parse_pauli_string}; use std::hint::black_box; const B: usize = 4096; @@ -213,7 +213,11 @@ fn observable(n: usize) -> (Vec, Vec) { let mut basis = Vec::new(); let mut coeffs = Vec::new(); for a in 0..n { - basis.push(parse_pauli_string(&pstr(n, &[(a, 'X')]), n).unwrap().0); + basis.push( + parse_pauli_string::(&pstr(n, &[(a, 'X')]), n) + .unwrap() + .0, + ); coeffs.push(1.0); } (basis, coeffs) @@ -227,8 +231,8 @@ fn bench_drug(c: &mut Criterion) { let h = hamiltonian_terms(n, &pairs, &bmag); let (ops, k) = kossakowski_model(n); - let spec_eig = LindbladSpec::new(n, &h, &eigenmode_jumps(&ops, &k)).unwrap(); - let mut spec_koss = LindbladSpec::new(n, &h, &[]).unwrap(); + let spec_eig = ::new(n, &h, &eigenmode_jumps(&ops, &k)).unwrap(); + let mut spec_koss = ::new(n, &h, &[]).unwrap(); spec_koss.add_kossakowski(&ops, &k).unwrap(); let cfg = PcStepConfig { diff --git a/crates/ppvm-lindblad/benches/kossakowski.rs b/crates/ppvm-lindblad/benches/kossakowski.rs index 22724a6d0..f47c08799 100644 --- a/crates/ppvm-lindblad/benches/kossakowski.rs +++ b/crates/ppvm-lindblad/benches/kossakowski.rs @@ -15,7 +15,7 @@ use criterion::{Criterion, criterion_group, criterion_main}; use num::Complex; -use ppvm_lindblad::{JumpInput, LindbladSpec, PcStepConfig, Word, parse_pauli_string}; +use ppvm_lindblad::{JumpInput, LindbladSpec, PcStepConfig, W_CHUNKS, Word, parse_pauli_string}; use std::f64::consts::PI; use std::hint::black_box; @@ -112,7 +112,7 @@ fn observable(n: usize, gam: &[Vec]) -> (Vec, Vec) { let mut basis = Vec::new(); let mut coeffs = Vec::new(); let mut push = |s: String, c: f64| { - basis.push(parse_pauli_string(&s, n).unwrap().0); + basis.push(parse_pauli_string::(&s, n).unwrap().0); coeffs.push(c); }; push(pstr(n, &[]), n as f64 * G0 / 2.0); @@ -133,8 +133,8 @@ fn bench_kossakowski(c: &mut Criterion) { let (j, gam) = chain_couplings(n); let h = hamiltonian_terms(n, &j); - let spec_eig = LindbladSpec::new(n, &h, &eigenmode_jumps(n, &gam)).unwrap(); - let mut spec_koss = LindbladSpec::new(n, &h, &[]).unwrap(); + let spec_eig = ::new(n, &h, &eigenmode_jumps(n, &gam)).unwrap(); + let mut spec_koss = ::new(n, &h, &[]).unwrap(); let ops: Vec<_> = (0..n).map(|q| sigma_minus(q, n)).collect(); let k: Vec>> = gam .iter() diff --git a/crates/ppvm-lindblad/examples/drug_profile.rs b/crates/ppvm-lindblad/examples/drug_profile.rs index 5c8f4e1bf..c82415d61 100644 --- a/crates/ppvm-lindblad/examples/drug_profile.rs +++ b/crates/ppvm-lindblad/examples/drug_profile.rs @@ -13,7 +13,7 @@ //! Usage: `cargo run --release --example drug_profile -- [N] [B] [STEPS]` use num::Complex; -use ppvm_lindblad::{LindbladSpec, PcStepConfig, Word, parse_pauli_string}; +use ppvm_lindblad::{LindbladSpec, PcStepConfig, W_CHUNKS, Word, parse_pauli_string}; use std::time::Instant; const N_M: usize = 5; @@ -147,7 +147,11 @@ fn observable(n: usize) -> (Vec, Vec) { let mut basis = Vec::new(); let mut coeffs = Vec::new(); for a in 0..n { - basis.push(parse_pauli_string(&pstr(n, &[(a, 'X')]), n).unwrap().0); + basis.push( + parse_pauli_string::(&pstr(n, &[(a, 'X')]), n) + .unwrap() + .0, + ); coeffs.push(1.0); } (basis, coeffs) @@ -161,7 +165,7 @@ fn main() { let dt = 1e-3; let (h, (ops, k)) = model(n); - let mut spec = LindbladSpec::new(n, &h, &[]).unwrap(); + let mut spec = ::new(n, &h, &[]).unwrap(); let t0 = Instant::now(); spec.add_kossakowski(&ops, &k).unwrap(); let build_ms = t0.elapsed().as_secs_f64() * 1e3; diff --git a/crates/ppvm-lindblad/src/algebra.rs b/crates/ppvm-lindblad/src/algebra.rs index 7a2b8c177..5fd3f4750 100644 --- a/crates/ppvm-lindblad/src/algebra.rs +++ b/crates/ppvm-lindblad/src/algebra.rs @@ -22,8 +22,8 @@ pub(crate) const COEFF_DROP_TOL: f64 = 1e-14; /// One Pauli term in a complex linear combination (a single summand of /// `L = Σ_a λ_a P_a`, or of a precomputed product such as `L†L`). #[derive(Clone)] -pub(crate) struct PauliTerm { - pub(crate) word: Word, +pub(crate) struct PauliTerm { + pub(crate) word: Word, pub(crate) coeff: Complex, } @@ -31,9 +31,12 @@ pub(crate) struct PauliTerm { /// as a Pauli linear combination, dropping FP-noise zeros. For `A = B` /// (the jump-operator `L†L`) the coefficients are real; in general they /// are complex. -pub(crate) fn precompute_adag_b(a_terms: &[PauliTerm], b_terms: &[PauliTerm]) -> Vec { +pub(crate) fn precompute_adag_b( + a_terms: &[PauliTerm], + b_terms: &[PauliTerm], +) -> Vec> { let zero = Complex::new(0.0, 0.0); - let mut acc: FxHashMap> = FxHashMap::default(); + let mut acc: FxHashMap, Complex> = FxHashMap::default(); for a in a_terms { for b in b_terms { let (word, phase) = pauli_mul(&a.word, &b.word); @@ -48,8 +51,8 @@ pub(crate) fn precompute_adag_b(a_terms: &[PauliTerm], b_terms: &[PauliTerm]) -> } /// Union of the supports (`xbits | zbits`) of every term, as raw chunks. -pub(crate) fn support_mask(terms: &[PauliTerm]) -> [Chunk; W_CHUNKS] { - let mut mask = [0 as Chunk; W_CHUNKS]; +pub(crate) fn support_mask(terms: &[PauliTerm]) -> [Chunk; C] { + let mut mask = [0 as Chunk; C]; for t in terms { for (i, slot) in mask.iter_mut().enumerate() { *slot |= t.word.xbits.data[i] | t.word.zbits.data[i]; @@ -73,9 +76,9 @@ pub(crate) fn phase_factor(phase: u8) -> Complex { /// Two Pauli strings anti-commute iff /// `popcount(a.x & b.z) + popcount(a.z & b.x)` is odd. #[inline(always)] -pub(crate) fn anti_commutes(a: &Word, b: &Word) -> bool { +pub(crate) fn anti_commutes(a: &Word, b: &Word) -> bool { let mut bits: u32 = 0; - for i in 0..W_CHUNKS { + for i in 0..C { bits += (a.xbits.data[i] & b.zbits.data[i]).count_ones(); bits += (a.zbits.data[i] & b.xbits.data[i]).count_ones(); } @@ -88,7 +91,7 @@ pub(crate) fn anti_commutes(a: &Word, b: &Word) -> bool { /// - `eps = -2.0` if `h·p` has phase `+i` (so `i·[h,p] = -2·out`), /// - `eps = +2.0` if `h·p` has phase `-i` (so `i·[h,p] = +2·out`). #[inline(always)] -pub(crate) fn comm_product(h: &Word, p: &Word) -> (Word, f64) { +pub(crate) fn comm_product(h: &Word, p: &Word) -> (Word, f64) { let (out, phase) = pauli_mul(h, p); let eps = match phase { 1 => -2.0, @@ -101,11 +104,11 @@ pub(crate) fn comm_product(h: &Word, p: &Word) -> (Word, f64) { /// Full Pauli product `p · q`: returns `(out, phase)` where the product /// is `ω · out` with `ω = i^phase`. #[inline(always)] -pub(crate) fn pauli_mul(p: &Word, q: &Word) -> (Word, u8) { - let mut out = Word::new(p.n_qubits()); +pub(crate) fn pauli_mul(p: &Word, q: &Word) -> (Word, u8) { + let mut out = Word::::new(p.n_qubits()); let mut sign_count: u32 = 0; let mut imag_count: u32 = 0; - for i in 0..W_CHUNKS { + for i in 0..C { let a = p.xbits.data[i]; let b = p.zbits.data[i]; let c = q.xbits.data[i]; diff --git a/crates/ppvm-lindblad/src/basis.rs b/crates/ppvm-lindblad/src/basis.rs index 64cbeeec8..09e162177 100644 --- a/crates/ppvm-lindblad/src/basis.rs +++ b/crates/ppvm-lindblad/src/basis.rs @@ -18,8 +18,8 @@ const CHUNK_SIZE: usize = 4096; /// Build a `word → row` map for a basis assumed to contain unique Pauli /// words; debug-asserts the uniqueness invariant. -pub fn build_basis_index(basis: &[Word]) -> FxHashMap { - let mut index: FxHashMap = FxHashMap::default(); +pub fn build_basis_index(basis: &[Word]) -> FxHashMap, u32> { + let mut index: FxHashMap, u32> = FxHashMap::default(); for (i, w) in basis.iter().enumerate() { let prev = index.insert(*w, i as u32); debug_assert!( @@ -32,15 +32,15 @@ pub fn build_basis_index(basis: &[Word]) -> FxHashMap { index } -impl LindbladSpec { +impl LindbladSpec { /// Off-basis component of `L*( Σ_j coeffs[j] · basis[j] )`. Output /// strings that lie in `basis` or in `protected` are dropped. pub fn leakage( &self, - basis: &[Word], + basis: &[Word], coeffs: &[f64], - protected: &[Word], - ) -> Result, Error> { + protected: &[Word], + ) -> Result, f64)>, Error> { self.leakage_with_prune(basis, coeffs, protected, usize::MAX, 0.0) } @@ -56,12 +56,12 @@ impl LindbladSpec { /// `room ≥ all candidates`, nothing is dropped — the near-exact case. pub fn leakage_with_prune( &self, - basis: &[Word], + basis: &[Word], coeffs: &[f64], - protected: &[Word], + protected: &[Word], max_basis: usize, tau_add: f64, - ) -> Result, Error> { + ) -> Result, f64)>, Error> { if basis.len() != coeffs.len() { return Err(Error::LengthMismatch { what: "basis and coeffs", @@ -80,16 +80,16 @@ impl LindbladSpec { let order = order_by_desc_mag(coeffs); let room = max_basis.saturating_sub(basis.len()); let n_qubits = self.n_qubits(); - let mut merged: FxHashMap = FxHashMap::default(); + let mut merged: FxHashMap, f64> = FxHashMap::default(); for chunk_indices in order.chunks(CHUNK_SIZE) { - let local: Vec> = chunk_indices + let local: Vec, f64)>> = chunk_indices .par_iter() .map_init( || { ( Vec::::with_capacity(n_qubits), Vec::::with_capacity(128), - FxHashMap::>::with_capacity_and_hasher( + FxHashMap::, Complex>::with_capacity_and_hasher( 128, FxBuildHasher::default(), ), @@ -132,7 +132,7 @@ impl LindbladSpec { /// /// Precondition: `basis` must not contain duplicate Pauli words /// (asserted in debug builds). - pub fn generator(&self, basis: &[Word]) -> Vec<(usize, usize, f64)> { + pub fn generator(&self, basis: &[Word]) -> Vec<(usize, usize, f64)> { let index = build_basis_index(basis); let n_qubits = self.n_qubits(); @@ -146,7 +146,7 @@ impl LindbladSpec { ( Vec::::with_capacity(n_qubits), Vec::::with_capacity(128), - FxHashMap::>::with_capacity_and_hasher( + FxHashMap::, Complex>::with_capacity_and_hasher( 128, FxBuildHasher::default(), ), @@ -178,10 +178,10 @@ impl LindbladSpec { /// component of `L*( Σ_j coeffs[j] · basis[j] )` with complex `coeffs`. pub fn leakage_complex( &self, - basis: &[Word], + basis: &[Word], coeffs: &[Complex], - protected: &[Word], - ) -> Result)>, Error> { + protected: &[Word], + ) -> Result, Complex)>, Error> { if basis.len() != coeffs.len() { return Err(Error::LengthMismatch { what: "basis and coeffs", @@ -194,12 +194,12 @@ impl LindbladSpec { protected.iter().map(|w| (word_hash(w), ())).collect(); let n_qubits = self.n_qubits(); - let mut merged: FxHashMap> = FxHashMap::default(); + let mut merged: FxHashMap, Complex> = FxHashMap::default(); for chunk_start in (0..basis.len()).step_by(CHUNK_SIZE) { let chunk_end = (chunk_start + CHUNK_SIZE).min(basis.len()); let chunk_basis = &basis[chunk_start..chunk_end]; let chunk_coeffs = &coeffs[chunk_start..chunk_end]; - let local: Vec)>> = chunk_basis + let local: Vec, Complex)>> = chunk_basis .par_iter() .zip(chunk_coeffs.par_iter()) .map_init( @@ -207,7 +207,7 @@ impl LindbladSpec { ( Vec::::with_capacity(n_qubits), Vec::::with_capacity(128), - FxHashMap::>::with_capacity_and_hasher( + FxHashMap::, Complex>::with_capacity_and_hasher( 128, FxBuildHasher::default(), ), @@ -260,12 +260,12 @@ impl LindbladSpec { /// (room ≥ all candidates) disables the cap — the near-exact case. pub fn leakage_orbit_rep( &self, - basis: &[Word], + basis: &[Word], coeffs: &[Complex], - protected: &[Word], + protected: &[Word], sector: &Sector<'_>, max_basis: usize, - ) -> Result)>, Error> { + ) -> Result, Complex)>, Error> { if basis.len() != coeffs.len() { return Err(Error::LengthMismatch { what: "basis and coeffs", @@ -275,22 +275,22 @@ impl LindbladSpec { } // Membership is tested on the canonical rep `r_q`, so unlike the // real path these are full-Word sets, not `word_hash` tables. - let in_basis: FxHashSet<&Word> = basis.iter().collect(); - let protected_set: FxHashSet<&Word> = protected.iter().collect(); + let in_basis: FxHashSet<&Word> = basis.iter().collect(); + let protected_set: FxHashSet<&Word> = protected.iter().collect(); let order = order_by_desc_mag(coeffs); let room = max_basis.saturating_sub(basis.len()); let n_qubits = self.n_qubits(); - let mut merged: FxHashMap> = FxHashMap::default(); + let mut merged: FxHashMap, Complex> = FxHashMap::default(); for chunk_indices in order.chunks(CHUNK_SIZE) { - let local: Vec)>> = chunk_indices + let local: Vec, Complex)>> = chunk_indices .par_iter() .map_init( || { ( Vec::::with_capacity(n_qubits), Vec::::with_capacity(128), - FxHashMap::>::with_capacity_and_hasher( + FxHashMap::, Complex>::with_capacity_and_hasher( 128, FxBuildHasher::default(), ), diff --git a/crates/ppvm-lindblad/src/error.rs b/crates/ppvm-lindblad/src/error.rs index 6c2235e86..1e10c53dd 100644 --- a/crates/ppvm-lindblad/src/error.rs +++ b/crates/ppvm-lindblad/src/error.rs @@ -3,14 +3,15 @@ //! 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 { + /// `got` qubits do not fit the word width in use, which holds `max`. TooManyQubits { got: usize, + max: usize, }, LengthMismatch { what: &'static str, @@ -51,11 +52,8 @@ pub enum Error { 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::TooManyQubits { got, max } => { + write!(f, "LindbladSpec supports n_qubits ≤ {max}; got {got}") } Error::LengthMismatch { what, a, b } => { write!(f, "{what}: expected matching lengths, got {a} and {b}") diff --git a/crates/ppvm-lindblad/src/kossakowski.rs b/crates/ppvm-lindblad/src/kossakowski.rs index 65aa76723..07c6220af 100644 --- a/crates/ppvm-lindblad/src/kossakowski.rs +++ b/crates/ppvm-lindblad/src/kossakowski.rs @@ -37,7 +37,7 @@ const HERMITICITY_TOL: f64 = 1e-10; /// Sandwich table of a pair, grouped by the left word: one /// `(P_a, [(P_b, coeff), …])` group per distinct `P_a`, so `P_a · p` is /// computed once per group and reused across its `P_b` partners. -type SandwichGroups = Vec<(Word, Vec<(Word, Complex)>)>; +type SandwichGroups = Vec<(Word, Vec<(Word, Complex)>)>; /// Which of the two compiled pair forms a [`Pair`] is. /// @@ -45,7 +45,7 @@ type SandwichGroups = Vec<(Word, Vec<(Word, Complex)>)>; /// list used by the one-sided commutator path, so it is modelled as a sum /// type rather than a flag: the off-diagonal-only term list cannot be /// reached on a diagonal pair. -pub(crate) enum PairShape { +pub(crate) enum PairShape { /// `n == m`. [`Pair::dd`] is `K_nn · A_n†A_n`. Diagonal, /// `n < m`, folding in the conjugate `(m, n)` pair. [`Pair::dd`] is the @@ -55,19 +55,19 @@ pub(crate) enum PairShape { /// The anti-Hermitian difference `−2i·Im(K_nm·A_n†A_m)` /// (pure-imaginary coefficients), used by the one-sided commutator /// of the folded conjugate pair. - dd_anti: Vec, + dd_anti: Vec>, }, } /// One compiled `(n, m)` pair of a Kossakowski dissipator. -pub(crate) struct Pair { - sand: SandwichGroups, +pub(crate) struct Pair { + sand: SandwichGroups, /// `A_n†A_m` scaled by `K_nm`; see [`PairShape`] for the exact form. - dd: Vec, - shape: PairShape, + dd: Vec>, + shape: PairShape, /// Support masks of `A_n` and `A_m`, for the one-sided fast path. - left_mask: [Chunk; W_CHUNKS], - right_mask: [Chunk; W_CHUNKS], + left_mask: [Chunk; C], + right_mask: [Chunk; C], } /// Compile a Kossakowski dissipator into one [`Pair`] per non-negligible @@ -75,11 +75,11 @@ pub(crate) struct Pair { /// /// Returns the pairs alongside, for each pair, the union support of its two /// operators, so the caller can index them by qubit. -pub(crate) fn compile( +pub(crate) fn compile( ops: &[Vec<(String, Complex)>], k: &[Vec>], n_qubits: usize, -) -> Result)>, Error> { +) -> Result, Vec)>, Error> { let max_abs = validate_k(k, ops.len())?; let (parsed, op_support) = parse_ops(ops, n_qubits)?; @@ -139,11 +139,14 @@ fn validate_k(k: &[Vec>], n_ops: usize) -> Result { /// Parsed operator table: the Pauli terms of each `A_n`, and each `A_n`'s /// union support. -type ParsedOps = (Vec>, Vec>); +type ParsedOps = (Vec>>, Vec>); /// Parse each operator's Pauli lincomb, returning the parsed terms and each /// operator's union support. -fn parse_ops(ops: &[Vec<(String, Complex)>], n_qubits: usize) -> Result { +fn parse_ops( + ops: &[Vec<(String, Complex)>], + n_qubits: usize, +) -> Result, Error> { let mut parsed = Vec::with_capacity(ops.len()); let mut supports = Vec::with_capacity(ops.len()); for (i, op) in ops.iter().enumerate() { @@ -164,12 +167,12 @@ fn parse_ops(ops: &[Vec<(String, Complex)>], n_qubits: usize) -> Result( + a_terms: &[PauliTerm], + b_terms: &[PauliTerm], k_nm: Complex, off_diag: bool, -) -> Pair { +) -> Pair { // A_n†A_m as `Σ γ_w W`, then scaled by K_nm. let adag_b = precompute_adag_b(a_terms, b_terms); let (dd, shape) = if off_diag { @@ -226,14 +229,14 @@ fn compile_pair( } } -impl Pair { +impl Pair { /// Accumulate this pair's contribution to `L*(p)` into `local`. - pub(crate) fn accumulate(&self, p: &Word, local: &mut FxHashMap>) { - let mut p_bits = [0 as Chunk; W_CHUNKS]; + pub(crate) fn accumulate(&self, p: &Word, local: &mut FxHashMap, Complex>) { + let mut p_bits = [0 as Chunk; C]; for (i, slot) in p_bits.iter_mut().enumerate() { *slot = p.xbits.data[i] | p.zbits.data[i]; } - let hits = |mask: &[Chunk; W_CHUNKS]| (0..W_CHUNKS).any(|i| mask[i] & p_bits[i] != 0); + let hits = |mask: &[Chunk; C]| (0..C).any(|i| mask[i] & p_bits[i] != 0); let (hit_l, hit_r) = (hits(&self.left_mask), hits(&self.right_mask)); // The pair is only visited when `p` overlaps at least one side, so @@ -260,9 +263,9 @@ impl Pair { /// from [`comm_product`], the term coefficient is `∓ t_c · (i/2) · eps`. fn accumulate_one_sided( &self, - p: &Word, + p: &Word, hit_r: bool, - local: &mut FxHashMap>, + local: &mut FxHashMap, Complex>, ) { let zero = Complex::new(0.0, 0.0); let (terms, half_i) = match &self.shape { @@ -296,7 +299,7 @@ impl Pair { /// `P_a` and reused across all its `P_b` partners. For a folded /// off-diagonal pair the sandwich is doubled and its real part taken /// (the conjugate `(m,n)` pair supplies the other half). - fn accumulate_both_sided(&self, p: &Word, local: &mut FxHashMap>) { + fn accumulate_both_sided(&self, p: &Word, local: &mut FxHashMap, Complex>) { let zero = Complex::new(0.0, 0.0); let fold = matches!(self.shape, PairShape::OffDiagonal { .. }); for (wa, rights) in &self.sand { diff --git a/crates/ppvm-lindblad/src/lib.rs b/crates/ppvm-lindblad/src/lib.rs index 55c2feb92..683e34751 100644 --- a/crates/ppvm-lindblad/src/lib.rs +++ b/crates/ppvm-lindblad/src/lib.rs @@ -28,8 +28,11 @@ //! FP noise). //! //! Pauli strings are stored as [`ppvm_pauli_word::word::PauliWord`] backed by -//! two 64-bit chunks (≤128 qubits; four 32-bit chunks on 32-bit targets) -//! with cached hashes for fast HashMap lookup. The hot-path commutator/ +//! a fixed array of `C` chunks (64-bit, or 32-bit on 32-bit targets) with +//! cached hashes for fast HashMap lookup. The crate is const-generic in `C` +//! ([`Word`], [`LindbladSpec`]); the default is the 128-qubit width, +//! and [`chunks_for`] picks the narrowest of the 128/256/512-qubit widths +//! for a register. The hot-path commutator/ //! product loops bypass the higher-level word API and operate directly on //! the raw chunks for speed. @@ -55,7 +58,10 @@ pub use error::Error; pub use sector::{Sector, canonicalize_basis_to_rep}; pub use spec::{JumpInput, LindbladSpec}; pub use step::PcStepTimings; -pub use word::{MAX_QUBITS, Word, codes_from_word, parse_pauli_string, word_from_codes}; +pub use word::{ + CHUNK_BITS, MAX_QUBITS, MAX_SUPPORTED_QUBITS, W_CHUNKS, WIDTHS, Word, chunks_for, + codes_from_word, max_qubits, parse_pauli_string, word_from_codes, +}; #[cfg(test)] mod tests; diff --git a/crates/ppvm-lindblad/src/mf_expm.rs b/crates/ppvm-lindblad/src/mf_expm.rs index 7f9f58323..caca82c5a 100644 --- a/crates/ppvm-lindblad/src/mf_expm.rs +++ b/crates/ppvm-lindblad/src/mf_expm.rs @@ -50,10 +50,10 @@ type PerCol = Vec<(f64, T)>; /// outputs (in- and out-of-basis, an upper bound on the column 1-norm) and /// `diag` the coefficient of the output Word equal to the input Word. The /// cache is reused by [`CscOp`] across every Krylov/Taylor matvec. -fn build_mf_cols( - spec: &LindbladSpec, - basis: &[Word], - index: &FxHashMap, +fn build_mf_cols( + spec: &LindbladSpec, + basis: &[Word], + index: &FxHashMap, u32>, ) -> (Cols, PerCol) { basis .par_iter() @@ -62,7 +62,7 @@ fn build_mf_cols( ( Vec::::with_capacity(spec.n_qubits()), Vec::::with_capacity(128), - FxHashMap::>::with_capacity_and_hasher( + FxHashMap::, Complex>::with_capacity_and_hasher( 128, FxBuildHasher::default(), ), @@ -113,10 +113,10 @@ fn build_mf_cols( /// upper bound: several distinct outputs `q` can share one rep, so the /// out-of-basis magnitudes are not attributable to a column of `M`. `diag` /// accumulates for the same reason. -fn build_orbit_rep_cols( - spec: &LindbladSpec, - basis: &[Word], - index: &FxHashMap, +fn build_orbit_rep_cols( + spec: &LindbladSpec, + basis: &[Word], + index: &FxHashMap, u32>, sector: &Sector<'_>, ) -> (Cols>, PerCol>) { basis @@ -127,7 +127,7 @@ fn build_orbit_rep_cols( ( Vec::::with_capacity(spec.n_qubits()), Vec::::with_capacity(128), - FxHashMap::>::with_capacity_and_hasher( + FxHashMap::, Complex>::with_capacity_and_hasher( 128, FxBuildHasher::default(), ), @@ -348,9 +348,9 @@ where /// ONE action pass builds the CSC cache `cols` (reused across every matvec) /// and, in the same pass, the `(raw, diag)` data the `μ`/1-norm selection /// needs; [`expm_apply_cached`] does the rest. -pub(crate) fn expm_apply_mf( - spec: &LindbladSpec, - basis: &[Word], +pub(crate) fn expm_apply_mf( + spec: &LindbladSpec, + basis: &[Word], dt: f64, coeffs: &[f64], drop_tol: f64, @@ -386,9 +386,9 @@ pub(crate) fn expm_apply_mf( /// The expensive phase-aware action is computed ONCE here (via /// [`build_orbit_rep_cols`]) and reused, CSC-style, across every /// Krylov–Taylor matvec, exactly as on the real path. -pub(crate) fn expm_apply_orbit_rep( - spec: &LindbladSpec, - basis: &[Word], +pub(crate) fn expm_apply_orbit_rep( + spec: &LindbladSpec, + basis: &[Word], sector: &Sector<'_>, dt: f64, coeffs: &[Complex], diff --git a/crates/ppvm-lindblad/src/sector.rs b/crates/ppvm-lindblad/src/sector.rs index 3973de327..181d9359b 100644 --- a/crates/ppvm-lindblad/src/sector.rs +++ b/crates/ppvm-lindblad/src/sector.rs @@ -88,7 +88,10 @@ impl<'a> Sector<'a> { /// is incompatible with `k`): the coefficient of such a rep is /// identically zero, so the term is dropped. #[inline] - pub fn canonicalize_phase(&self, q: &Word) -> Option<(Word, Complex, usize)> { + pub fn canonicalize_phase( + &self, + q: &Word, + ) -> Option<(Word, Complex, usize)> { let (rep, idx, orbit_size) = self .group .canonicalize_in_sector_indexed(q, &self.characters)?; @@ -105,7 +108,7 @@ impl<'a> Sector<'a> { /// `ĉ_r = |orbit_r| · c_r` (what `momentum_merge_pauli_sum_pair` /// uses). It is `|G|` only for free orbits. #[inline] - pub fn orbit_size(&self, w: &Word) -> Option { + pub fn orbit_size(&self, w: &Word) -> Option { self.group .canonicalize_in_sector_indexed(w, &self.characters) .map(|(_, _, orbit_size)| orbit_size) @@ -119,7 +122,7 @@ impl<'a> Sector<'a> { /// /// Does NOT deduplicate — if multiple input entries collapse to the /// same rep, both are kept (caller should run a merge afterwards). -pub fn canonicalize_basis_to_rep(basis: &mut [Word], group: &TranslationGroup) { +pub fn canonicalize_basis_to_rep(basis: &mut [Word], group: &TranslationGroup) { for w in basis.iter_mut() { *w = group.canonicalize(w); } diff --git a/crates/ppvm-lindblad/src/spec.rs b/crates/ppvm-lindblad/src/spec.rs index 78b48be05..013884d0a 100644 --- a/crates/ppvm-lindblad/src/spec.rs +++ b/crates/ppvm-lindblad/src/spec.rs @@ -8,31 +8,31 @@ use crate::algebra::{ PauliTerm, anti_commutes, comm_product, pauli_mul, phase_factor, precompute_adag_b, }; use crate::kossakowski; -use crate::word::{MAX_QUBITS, Word, parse_pauli_string, word_support}; +use crate::word::{W_CHUNKS, Word, check_width, parse_pauli_string, word_support}; use fxhash::FxHashMap; use num::Complex; /// Parsed Hamiltonian term. #[derive(Clone)] -struct HTerm { - word: Word, +struct HTerm { + word: Word, coeff: f64, } /// One entry of the dissipator. `HermitianPauli` and `General` are the /// jump-operator form (`K` diagonal); `Kossakowski` is one compiled pair of /// the general form. See [`crate::kossakowski`]. -enum JumpKind { +enum JumpKind { HermitianPauli { - word: Word, + word: Word, rate: f64, }, General { - terms: Vec, // L = Σ_a λ_a P_a - dagger_dagger: Vec, // L†L = Σ_c μ_c P_c (μ_c ∈ ℝ) + terms: Vec>, // L = Σ_a λ_a P_a + dagger_dagger: Vec>, // L†L = Σ_c μ_c P_c (μ_c ∈ ℝ) rate: f64, }, - Kossakowski(kossakowski::Pair), + Kossakowski(kossakowski::Pair), } /// Union of `index[q]` for each `q ∈ p_support`, deduped. @@ -52,10 +52,10 @@ fn candidate_terms(p_support: &[u32], index: &[Vec], scratch: &mut Vec /// call rather than cached: for sparse-local Hamiltonians a per-word cache /// costs more than the recompute (hash lookup ≳ recompute) and its several /// KB per cached word dominate memory at large basis sizes. -pub struct LindbladSpec { +pub struct LindbladSpec { n_qubits: usize, - h_terms: Vec, - j_kinds: Vec, + h_terms: Vec>, + j_kinds: Vec>, /// `h_support[q]` = indices of Hamiltonian terms acting on qubit `q`. h_support: Vec>, /// `j_support[q]` = indices of jumps whose support contains qubit `q`. @@ -72,7 +72,7 @@ pub struct JumpInput { pub rate: f64, } -impl LindbladSpec { +impl LindbladSpec { /// Construct a Lindbladian spec from Hamiltonian terms and jump operators. /// /// `h_terms` are `(pauli_string, coefficient)` pairs forming the Hermitian @@ -84,11 +84,9 @@ impl LindbladSpec { h_terms: &[(String, f64)], jumps: &[JumpInput], ) -> Result { - if n_qubits > MAX_QUBITS { - return Err(Error::TooManyQubits { got: n_qubits }); - } + check_width::(n_qubits)?; - let mut h_parsed: Vec = Vec::with_capacity(h_terms.len()); + let mut h_parsed: Vec> = Vec::with_capacity(h_terms.len()); let mut h_support_idx: Vec> = vec![Vec::new(); n_qubits]; for (i, (s, c)) in h_terms.iter().enumerate() { let (word, support) = parse_pauli_string(s, n_qubits)?; @@ -98,7 +96,7 @@ impl LindbladSpec { h_parsed.push(HTerm { word, coeff: *c }); } - let mut j_kinds: Vec = Vec::with_capacity(jumps.len()); + let mut j_kinds: Vec> = Vec::with_capacity(jumps.len()); let mut j_support_idx: Vec> = vec![Vec::new(); n_qubits]; for (k, jump) in jumps.iter().enumerate() { if jump.rate < 0.0 { @@ -126,7 +124,7 @@ impl LindbladSpec { } // General path: parse all terms, precompute L†L, record union support. - let mut terms: Vec = Vec::with_capacity(jump.lincomb.len()); + let mut terms: Vec> = Vec::with_capacity(jump.lincomb.len()); let mut union_support: std::collections::BTreeSet = std::collections::BTreeSet::new(); for (s, c) in &jump.lincomb { @@ -192,8 +190,8 @@ impl LindbladSpec { /// Apply `L*` to a single Pauli string `p`. Returns the output Pauli /// strings and their real coefficients (zero entries omitted). - pub fn action(&self, p: &Word) -> Vec<(Word, f64)> { - let mut out: FxHashMap = FxHashMap::default(); + pub fn action(&self, p: &Word) -> Vec<(Word, f64)> { + let mut out: FxHashMap, f64> = FxHashMap::default(); let mut s1 = Vec::new(); let mut s2 = Vec::new(); self.accumulate_action(p, 1.0, &mut out, &mut s1, &mut s2); @@ -204,11 +202,11 @@ impl LindbladSpec { /// `L*(p)` contributes (without the input coefficient). pub(crate) fn compute_action_terms( &self, - p: &Word, + p: &Word, scratch_support: &mut Vec, scratch_cands: &mut Vec, - scratch_local: &mut FxHashMap>, - ) -> Vec<(Word, f64)> { + scratch_local: &mut FxHashMap, Complex>, + ) -> Vec<(Word, f64)> { word_support(p, scratch_support); let zero = Complex::new(0.0, 0.0); scratch_local.clear(); @@ -283,9 +281,9 @@ impl LindbladSpec { /// Accumulate `scale · L*(p)` into `out`. fn accumulate_action( &self, - p: &Word, + p: &Word, scale: f64, - out: &mut FxHashMap, + out: &mut FxHashMap, f64>, scratch_support: &mut Vec, scratch_cands: &mut Vec, ) { diff --git a/crates/ppvm-lindblad/src/step.rs b/crates/ppvm-lindblad/src/step.rs index 97e593ead..e5dce345e 100644 --- a/crates/ppvm-lindblad/src/step.rs +++ b/crates/ppvm-lindblad/src/step.rs @@ -50,7 +50,7 @@ impl Phase { } } -impl LindbladSpec { +impl LindbladSpec { /// One predictor-corrector step `O ← exp(dt·L*) O` in the adaptive /// real-coefficient Pauli basis: first-hop leakage admission, predictor /// exponential, second-hop admission from the predicted state, corrector @@ -61,10 +61,10 @@ impl LindbladSpec { /// `protected` words are never dropped. All tuning knobs live in `cfg`. pub fn pc_step( &self, - basis: &mut Vec, + basis: &mut Vec>, coeffs: &mut Vec, dt: f64, - protected: &[Word], + protected: &[Word], cfg: &PcStepConfig, ) -> Result<(), Error> { self.run_in_pool(cfg, |this| { @@ -78,10 +78,10 @@ impl LindbladSpec { /// spots. pub fn pc_step_timed( &self, - basis: &mut Vec, + basis: &mut Vec>, coeffs: &mut Vec, dt: f64, - protected: &[Word], + protected: &[Word], cfg: &PcStepConfig, ) -> Result { self.run_in_pool(cfg, |this| { @@ -107,10 +107,10 @@ impl LindbladSpec { fn pc_step_inner( &self, - basis: &mut Vec, + basis: &mut Vec>, coeffs: &mut Vec, dt: f64, - protected: &[Word], + protected: &[Word], cfg: &PcStepConfig, timed: bool, ) -> Result { @@ -175,7 +175,7 @@ impl LindbladSpec { /// Compute `exp(dt · M) · b` for the in-basis-restricted generator /// `M`, matrix-free, via `quspin-expm` (see [`crate::mf_expm`]). - fn expm_step(&self, basis: &[Word], dt: f64, b: &[f64], drop_tol: f64) -> Vec { + fn expm_step(&self, basis: &[Word], dt: f64, b: &[f64], drop_tol: f64) -> Vec { mf_expm::expm_apply_mf(self, basis, dt, b, drop_tol) } @@ -202,10 +202,10 @@ impl LindbladSpec { /// Honours `cfg.num_threads` the same way [`Self::pc_step`] does. pub fn pc_step_orbit_rep( &self, - basis: &mut Vec, + basis: &mut Vec>, coeffs: &mut Vec>, dt: f64, - protected: &[Word], + protected: &[Word], sector: &Sector<'_>, cfg: &PcStepConfig, ) -> Result<(), Error> { @@ -216,10 +216,10 @@ impl LindbladSpec { fn pc_step_orbit_rep_inner( &self, - basis: &mut Vec, + basis: &mut Vec>, coeffs: &mut Vec>, dt: f64, - protected: &[Word], + protected: &[Word], sector: &Sector<'_>, cfg: &PcStepConfig, ) -> Result<(), Error> { diff --git a/crates/ppvm-lindblad/src/tests.rs b/crates/ppvm-lindblad/src/tests.rs index 7a1e4a81a..8d4651ef0 100644 --- a/crates/ppvm-lindblad/src/tests.rs +++ b/crates/ppvm-lindblad/src/tests.rs @@ -45,13 +45,13 @@ fn jump_hpauli(s: &str, rate: f64) -> JumpInput { #[test] fn z_dephasing_action_on_x() { // L = Z on a single qubit; L*(X) = γ(ZXZ - X) = γ(-X - X) = -2γ X. - let spec = LindbladSpec::new( + let spec = ::new( 1, &[("X".to_string(), 0.0)], // no Hamiltonian &[jump_hpauli("Z", 0.5)], ) .unwrap(); - let (x, _) = parse_pauli_string("X", 1).unwrap(); + let (x, _) = parse_pauli_string::("X", 1).unwrap(); let terms = spec.action(&x); assert_eq!(terms.len(), 1); assert!((terms[0].1 - (-1.0)).abs() < 1e-12); // -2·0.5 = -1 @@ -68,10 +68,10 @@ fn amplitude_damping_action_on_z() { ], rate: 1.0, }; - let spec = LindbladSpec::new(1, &[], &[sigma_minus]).unwrap(); - let (z, _) = parse_pauli_string("Z", 1).unwrap(); + let spec = ::new(1, &[], &[sigma_minus]).unwrap(); + let (z, _) = parse_pauli_string::("Z", 1).unwrap(); let terms = spec.action(&z); - let (i_word, _) = parse_pauli_string("I", 1).unwrap(); + let (i_word, _) = parse_pauli_string::("I", 1).unwrap(); let mut i_coeff = 0.0; let mut z_coeff = 0.0; for (w, c) in &terms { @@ -88,7 +88,7 @@ fn amplitude_damping_action_on_z() { #[test] fn word_codec_roundtrip() { let codes = [0u8, 1, 2, 3, 1, 0, 3, 2]; - let w = word_from_codes(&codes).unwrap(); + let w = word_from_codes::(&codes).unwrap(); let mut out = vec![0u8; codes.len()]; codes_from_word(&w, &mut out); assert_eq!(out.as_slice(), &codes); @@ -127,11 +127,11 @@ fn assert_orbit_rep_matches_projection( ) { use ppvm_pauli_sum::symmetry::canonicalize_pauli_sum_complex; - let spec = LindbladSpec::new(n, h_terms, &[]).unwrap(); + let spec = ::new(n, h_terms, &[]).unwrap(); let group = ppvm_pauli_sum::symmetry::TranslationGroup::chain_1d(n); let basis_full: Vec = seed .iter() - .map(|(s, _)| parse_pauli_string(s, n).unwrap().0) + .map(|(s, _)| parse_pauli_string::(s, n).unwrap().0) .collect(); let coeffs_full: Vec> = seed.iter().map(|(_, c)| *c).collect(); @@ -252,14 +252,14 @@ fn complex_full_matches_real_at_kzero() { h_terms.push((s.into_iter().collect(), 1.0)); } } - let spec = LindbladSpec::new(n, &h_terms, &[]).unwrap(); + let spec = ::new(n, &h_terms, &[]).unwrap(); let mut basis_r: Vec = (0..n) .map(|j| { let mut s = vec!['I'; n]; s[j] = 'Z'; let st: String = s.into_iter().collect(); - let (w, _) = parse_pauli_string(&st, n).unwrap(); + let (w, _) = parse_pauli_string::(&st, n).unwrap(); w }) .collect(); @@ -344,7 +344,7 @@ fn pc_step_matches_symmetry_merged_on_small_chain() { } } // No dissipation. - let spec = LindbladSpec::new(n, &h_terms, &[]).unwrap(); + let spec = ::new(n, &h_terms, &[]).unwrap(); let group = TranslationGroup::chain_1d(n); // Initial: O(0) = Σ_j Z_j (translation-invariant). @@ -353,7 +353,7 @@ fn pc_step_matches_symmetry_merged_on_small_chain() { let mut s = vec!['I'; n]; s[j] = 'Z'; let st: String = s.into_iter().collect(); - let (w, _) = parse_pauli_string(&st, n).unwrap(); + let (w, _) = parse_pauli_string::(&st, n).unwrap(); w }) .collect(); diff --git a/crates/ppvm-lindblad/src/truncate.rs b/crates/ppvm-lindblad/src/truncate.rs index c96a927fb..1fe7a9fb3 100644 --- a/crates/ppvm-lindblad/src/truncate.rs +++ b/crates/ppvm-lindblad/src/truncate.rs @@ -24,7 +24,10 @@ pub(crate) fn order_by_desc_mag(coeffs: &[T]) -> Vec { /// candidate map — `room` being the number of strings we could actually /// admit to the basis, so there is no point tracking more. Applied after /// each accumulation chunk. -pub(crate) fn cap_map_to_room(merged: &mut FxHashMap, room: usize) { +pub(crate) fn cap_map_to_room( + merged: &mut FxHashMap, T>, + room: usize, +) { if merged.len() <= room { return; } @@ -41,17 +44,17 @@ pub(crate) fn cap_map_to_room(merged: &mut FxHashMap, room: u /// Compact `basis` / `coeffs` in place: drop entries whose coefficient /// magnitude is below `drop_tol` unless the word appears in `protected`. /// No-op when `drop_tol ≤ 0`. -pub(crate) fn prune_basis( - basis: &mut Vec, +pub(crate) fn prune_basis( + basis: &mut Vec>, coeffs: &mut Vec, drop_tol: f64, - protected: &[Word], + protected: &[Word], ) { if drop_tol <= 0.0 { return; } debug_assert_eq!(basis.len(), coeffs.len()); - let protected_set: FxHashSet<&Word> = protected.iter().collect(); + let protected_set: FxHashSet<&Word> = protected.iter().collect(); retain_in_place(basis, coeffs, |w, c| { c.mag() >= drop_tol || protected_set.contains(w) }); @@ -61,16 +64,16 @@ pub(crate) fn prune_basis( /// `max_basis` largest-magnitude terms (protected strings always kept), /// dropping the rest. Rank-based total-basis bound; dual of `drop_tol`. /// A `max_basis` large enough to cover the whole basis is a no-op. -pub(crate) fn cap_basis( - basis: &mut Vec, +pub(crate) fn cap_basis( + basis: &mut Vec>, coeffs: &mut Vec, max_basis: usize, - protected: &[Word], + protected: &[Word], ) { if basis.len() <= max_basis { return; } - let protected_set: FxHashSet<&Word> = protected.iter().collect(); + let protected_set: FxHashSet<&Word> = protected.iter().collect(); let n_prot = basis.iter().filter(|w| protected_set.contains(w)).count(); let slots = max_basis.saturating_sub(n_prot); let mut mags: Vec = basis @@ -96,10 +99,10 @@ pub(crate) fn cap_basis( /// expm/leakage peak memory) never exceeds `max_basis`. New strings get /// coefficient 0; the surrounding expm fills them. No magnitude filter: the /// top-`room` by `|leakage|` are added (a large `max_basis` adds them all). -pub(crate) fn add_leakage_capped( - basis: &mut Vec, +pub(crate) fn add_leakage_capped( + basis: &mut Vec>, coeffs: &mut Vec, - mut leak: Vec<(Word, T)>, + mut leak: Vec<(Word, T)>, max_basis: usize, ) { let room = max_basis.saturating_sub(basis.len()); @@ -117,10 +120,10 @@ pub(crate) fn add_leakage_capped( /// Keep the `basis`/`coeffs` entries satisfying `keep`, preserving order, /// by swapping survivors down and truncating. -fn retain_in_place( - basis: &mut Vec, +fn retain_in_place( + basis: &mut Vec>, coeffs: &mut Vec, - mut keep: impl FnMut(&Word, &T) -> bool, + mut keep: impl FnMut(&Word, &T) -> bool, ) { let mut write = 0; for read in 0..basis.len() { diff --git a/crates/ppvm-lindblad/src/word.rs b/crates/ppvm-lindblad/src/word.rs index 15c762d9b..d79e9fe57 100644 --- a/crates/ppvm-lindblad/src/word.rs +++ b/crates/ppvm-lindblad/src/word.rs @@ -16,31 +16,68 @@ pub(crate) type Chunk = u64; #[cfg(not(target_pointer_width = "64"))] pub(crate) type Chunk = u32; -/// Chunks per word; words pack up to 128 qubits on every target. -#[cfg(target_pointer_width = "64")] -pub(crate) const W_CHUNKS: usize = 2; -#[cfg(not(target_pointer_width = "64"))] -pub(crate) const W_CHUNKS: usize = 4; +/// Bits per storage chunk. +pub const CHUNK_BITS: usize = Chunk::BITS as usize; + +/// Chunks for a 128-qubit word: the default width, and the only one +/// instantiated for registers of at most 128 qubits. +pub const W_CHUNKS: usize = 128 / CHUNK_BITS; + +/// Chunk counts of the word widths the Python layer instantiates, narrowest +/// first: 128, 256 and 512 qubits. [`chunks_for`] picks among them. +pub const WIDTHS: [usize; 3] = [128 / CHUNK_BITS, 256 / CHUNK_BITS, 512 / CHUNK_BITS]; + +/// Largest register any instantiated width supports. +pub const MAX_SUPPORTED_QUBITS: usize = 512; -/// Maximum number of qubits supported by [`Word`]. -pub const MAX_QUBITS: usize = 128; +/// Maximum number of qubits of the default-width [`Word`] (128). +pub const MAX_QUBITS: usize = max_qubits::(); -/// The Pauli-word storage type used throughout this crate. +/// Capacity of a `C`-chunk [`Word`], in qubits. +pub const fn max_qubits() -> usize { + C * CHUNK_BITS +} + +/// The narrowest entry of [`WIDTHS`] that holds `n_qubits`, or `None` above +/// [`MAX_SUPPORTED_QUBITS`]. +pub const fn chunks_for(n_qubits: usize) -> Option { + let mut i = 0; + while i < WIDTHS.len() { + if n_qubits <= WIDTHS[i] * CHUNK_BITS { + return Some(WIDTHS[i]); + } + i += 1; + } + None +} + +/// The Pauli-word storage type used throughout this crate, `C` chunks wide. /// -/// `[Chunk; W_CHUNKS]` covers up to 128 qubits; the `FxBuildHasher` -/// matches the hash used by the `FxHashMap` keys we wrap with; -/// `REHASH=true` means `set()` keeps the cached hash in sync. -pub type Word = PauliWord<[Chunk; W_CHUNKS], FxBuildHasher, true>; +/// The crate is const-generic in the chunk count so a register of any size +/// up to [`MAX_SUPPORTED_QUBITS`] gets a fixed-size word; the default +/// `C = W_CHUNKS` covers 128 qubits, byte-for-byte the historical layout. +/// The `FxBuildHasher` matches the hash used by the `FxHashMap` keys we +/// wrap with; `REHASH=true` means `set()` keeps the cached hash in sync. +pub type Word = PauliWord<[Chunk; C], FxBuildHasher, true>; + +/// `Err(TooManyQubits)` unless a `C`-chunk word holds `n_qubits`. +pub(crate) fn check_width(n_qubits: usize) -> Result<(), Error> { + if n_qubits > max_qubits::() { + return Err(Error::TooManyQubits { + got: n_qubits, + max: max_qubits::(), + }); + } + Ok(()) +} /// Build a [`Word`] from a length-`n_qubits` slice of Pauli labels /// (`0=I, 1=X, 2=Z, 3=Y` — the [`ppvm_traits::char::Pauli`] discriminants). /// Sets all bits and rehashes once. -pub fn word_from_codes(codes: &[u8]) -> Result { +pub fn word_from_codes(codes: &[u8]) -> Result, Error> { let n_qubits = codes.len(); - if n_qubits > MAX_QUBITS { - return Err(Error::TooManyQubits { got: n_qubits }); - } - let mut w = Word::new(n_qubits); + check_width::(n_qubits)?; + let mut w = Word::::new(n_qubits); for (q, &b) in codes.iter().enumerate() { if b > 3 { return Err(Error::InvalidPauliCode { code: b }); @@ -57,7 +94,7 @@ pub fn word_from_codes(codes: &[u8]) -> Result { } /// Inverse of [`word_from_codes`]: write `n_qubits` Pauli labels into `out`. -pub fn codes_from_word(w: &Word, out: &mut [u8]) { +pub fn codes_from_word(w: &Word, out: &mut [u8]) { debug_assert_eq!(out.len(), w.n_qubits()); for (q, slot) in out.iter_mut().enumerate() { let xb = w.xbits[q] as u8; @@ -68,10 +105,11 @@ pub fn codes_from_word(w: &Word, out: &mut [u8]) { /// Parse a `"IXYZ..."` string into a [`Word`] together with the list of /// qubits where the Pauli is non-identity (the term's support). -pub fn parse_pauli_string(s: &str, n_qubits: usize) -> Result<(Word, Vec), Error> { - if n_qubits > MAX_QUBITS { - return Err(Error::TooManyQubits { got: n_qubits }); - } +pub fn parse_pauli_string( + s: &str, + n_qubits: usize, +) -> Result<(Word, Vec), Error> { + check_width::(n_qubits)?; let chars: Vec = s.chars().filter(|c| *c != '_').collect(); if chars.len() != n_qubits { return Err(Error::WrongLength { @@ -79,7 +117,7 @@ pub fn parse_pauli_string(s: &str, n_qubits: usize) -> Result<(Word, Vec), got: chars.len(), }); } - let mut w = Word::new(n_qubits); + let mut w = Word::::new(n_qubits); let mut support = Vec::new(); for (q, c) in chars.into_iter().enumerate() { match c { @@ -105,7 +143,7 @@ pub fn parse_pauli_string(s: &str, n_qubits: usize) -> Result<(Word, Vec), } /// Compute the support (non-identity qubits) of `w`. -pub(crate) fn word_support(w: &Word, out: &mut Vec) { +pub(crate) fn word_support(w: &Word, out: &mut Vec) { out.clear(); for q in 0..w.n_qubits() { if w.xbits[q] || w.zbits[q] { @@ -120,7 +158,7 @@ pub(crate) fn word_support(w: &Word, out: &mut Vec) { /// cached hash once through `FxHasher` and never touches the 32-byte /// payload. #[inline(always)] -pub(crate) fn word_hash(w: &Word) -> u64 { +pub(crate) fn word_hash(w: &Word) -> u64 { use std::hash::{Hash, Hasher}; let mut h = fxhash::FxHasher::default(); w.hash(&mut h); diff --git a/crates/ppvm-lindblad/tests/word_width.rs b/crates/ppvm-lindblad/tests/word_width.rs new file mode 100644 index 000000000..00c5830fe --- /dev/null +++ b/crates/ppvm-lindblad/tests/word_width.rs @@ -0,0 +1,274 @@ +// SPDX-FileCopyrightText: 2026 The PPVM Authors +// SPDX-License-Identifier: Apache-2.0 + +//! Pauli-word width is a const-generic parameter, so a `LindbladSpec` can be +//! instantiated wider than the historical 128-qubit ceiling. +//! +//! The central check is a *padding invariance*: the same physical problem, +//! embedded in words of different widths, must produce bit-identical numbers. +//! A width bug (mask truncation, a stray `W_CHUNKS`, a wrong support index) +//! breaks that immediately. + +use num::Complex; +use ppvm_lindblad::{ + CHUNK_BITS, JumpInput, LindbladSpec, MAX_SUPPORTED_QUBITS, PcStepConfig, Sector, W_CHUNKS, + WIDTHS, Word, chunks_for, codes_from_word, max_qubits, word_from_codes, +}; +use ppvm_pauli_sum::symmetry::{TranslationGroup, canonicalize_pauli_sum_complex}; + +const W128: usize = WIDTHS[0]; +const W256: usize = WIDTHS[1]; +const W512: usize = WIDTHS[2]; + +/// `H = J Σ_{i (Vec<(String, f64)>, Vec) { + let pad = |sites: &[(usize, char)]| -> String { + let mut s = vec!['I'; n_total]; + for &(q, c) in sites { + s[q] = c; + } + s.into_iter().collect() + }; + let mut h = Vec::new(); + for i in 0..n_active - 1 { + h.push((pad(&[(i, 'Z'), (i + 1, 'Z')]), 0.7)); + } + for i in 0..n_active { + h.push((pad(&[(i, 'X')]), 1.3)); + } + let jumps = (0..n_active) + .map(|i| JumpInput { + lincomb: vec![(pad(&[(i, 'Z')]), Complex::new(1.0, 0.0))], + rate: 0.05, + }) + .collect(); + (h, jumps) +} + +/// A decay-type Kossakowski dissipator `A_0 = σ⁻_0`, `A_1 = σ⁻_1` with a +/// non-diagonal pair matrix, on the first two qubits of `n_total`. +fn add_kossakowski(spec: &mut LindbladSpec, n_total: usize) { + let op = |q: usize| -> Vec<(String, Complex)> { + let mut x = vec!['I'; n_total]; + let mut y = vec!['I'; n_total]; + x[q] = 'X'; + y[q] = 'Y'; + vec![ + (x.into_iter().collect(), Complex::new(0.5, 0.0)), + (y.into_iter().collect(), Complex::new(0.0, -0.5)), + ] + }; + let k = vec![ + vec![Complex::new(0.3, 0.0), Complex::new(0.1, 0.05)], + vec![Complex::new(0.1, -0.05), Complex::new(0.2, 0.0)], + ]; + spec.add_kossakowski(&[op(0), op(1)], &k).unwrap(); +} + +/// Evolve `Z_0 Z_1` for a few steps and return the coefficient sum over +/// `{I, Z}`-only strings — i.e. the expectation on the all-`Z = -1` product +/// state, up to the sign convention (identical across widths, which is all +/// this test needs). +fn evolve(n_total: usize, n_active: usize, steps: usize, koss: bool) -> f64 { + let (h, jumps) = model(n_total, n_active); + let mut spec = LindbladSpec::::new(n_total, &h, &jumps).unwrap(); + if koss { + add_kossakowski(&mut spec, n_total); + } + + let mut codes = vec![0u8; n_total]; + codes[0] = 2; // Z + codes[1] = 2; // Z + let mut basis = vec![word_from_codes::(&codes).unwrap()]; + let mut coeffs = vec![1.0f64]; + + let cfg = PcStepConfig { + max_basis: 20_000, + admit_basis: Some(60_000), + drop_tol: 0.0, + tau_add: None, + num_threads: Some(1), + }; + for _ in 0..steps { + spec.pc_step(&mut basis, &mut coeffs, 0.05, &[], &cfg) + .unwrap(); + } + + let mut out = vec![0u8; n_total]; + let mut acc = 0.0; + for (w, c) in basis.iter().zip(&coeffs) { + codes_from_word(w, &mut out); + if out.iter().all(|&b| b == 0 || b == 2) { + let nz = out.iter().filter(|&&b| b == 2).count(); + acc += if nz % 2 == 0 { *c } else { -*c }; + } + } + acc +} + +/// Momentum-sector evolution of `Σ_q X_q` under a translation-invariant ring +/// `H = Σ (X X + Y Y + ½ Z Z) + 0.3 Σ Z` of `n` sites, untruncated. Returns +/// the rep coefficients sorted by word. +fn evolve_orbit(n: usize, k: i32, steps: usize) -> Vec<(Vec, Complex)> { + let word = |ops: &[(usize, char)]| -> String { + let mut s = vec!['I'; n]; + for &(q, c) in ops { + s[q] = c; + } + s.into_iter().collect() + }; + let mut h = Vec::new(); + for i in 0..n { + let j = (i + 1) % n; + h.push((word(&[(i, 'X'), (j, 'X')]), 1.0)); + h.push((word(&[(i, 'Y'), (j, 'Y')]), 1.0)); + h.push((word(&[(i, 'Z'), (j, 'Z')]), 0.5)); + h.push((word(&[(i, 'Z')]), 0.3)); + } + let spec = LindbladSpec::::new(n, &h, &[]).unwrap(); + let group = TranslationGroup::chain_1d(n); + let k_modes = [k]; + + let mut basis: Vec> = (0..n) + .map(|q| { + let mut codes = vec![0u8; n]; + codes[q] = 1; // X + word_from_codes::(&codes).unwrap() + }) + .collect(); + let mut coeffs: Vec> = (0..n) + .map(|q| { + Complex::from_polar( + 1.0, + -2.0 * std::f64::consts::PI * (k as f64) * q as f64 / n as f64, + ) + }) + .collect(); + canonicalize_pauli_sum_complex(&mut basis, &mut coeffs, &group, &k_modes); + + let sector = Sector::new(&group, &k_modes); + let cfg = PcStepConfig { + max_basis: usize::MAX, + admit_basis: None, + drop_tol: 0.0, + tau_add: None, + num_threads: Some(1), + }; + for _ in 0..steps { + spec.pc_step_orbit_rep(&mut basis, &mut coeffs, 0.1, &[], §or, &cfg) + .unwrap(); + } + let mut out: Vec<(Vec, Complex)> = basis + .iter() + .zip(&coeffs) + .map(|(w, c)| { + let mut codes = vec![0u8; n]; + codes_from_word(w, &mut codes); + (codes, *c) + }) + .collect(); + out.sort_by(|a, b| a.0.cmp(&b.0)); + out +} + +#[test] +fn capacity_scales_with_chunk_count() { + assert_eq!(max_qubits::(), 128); + assert_eq!(max_qubits::(), 256); + assert_eq!(max_qubits::(), 512); + assert_eq!(W_CHUNKS, W128); + assert_eq!(W128 * CHUNK_BITS, 128); + assert_eq!(chunks_for(1), Some(W128)); + assert_eq!(chunks_for(128), Some(W128)); + assert_eq!(chunks_for(129), Some(W256)); + assert_eq!(chunks_for(130), Some(W256)); + assert_eq!(chunks_for(512), Some(W512)); + assert_eq!(chunks_for(MAX_SUPPORTED_QUBITS + 1), None); +} + +#[test] +fn wide_words_exceed_the_old_128_qubit_ceiling() { + // The exact case that used to fail with "supports n_qubits ≤ 128". + let (h, jumps) = model(130, 4); + let spec = LindbladSpec::::new(130, &h, &jumps).unwrap(); + assert_eq!(spec.n_qubits(), 130); + + let (h, jumps) = model(512, 4); + let spec = LindbladSpec::::new(512, &h, &jumps).unwrap(); + assert_eq!(spec.n_qubits(), 512); +} + +#[test] +fn too_many_qubits_for_the_width_is_rejected() { + let (h, jumps) = model(200, 4); + let Err(err) = LindbladSpec::::new(200, &h, &jumps) else { + panic!("a 200-qubit spec must not fit 128-qubit words"); + }; + assert_eq!( + err.to_string(), + "LindbladSpec supports n_qubits ≤ 128; got 200" + ); + assert!(LindbladSpec::::new(200, &h, &jumps).is_ok()); + assert!(word_from_codes::(&[0u8; 129]).is_err()); +} + +#[test] +fn padding_into_a_wider_word_changes_nothing() { + // Identical 6-qubit physics, embedded in 64-, 200- and 400-qubit + // registers backed by 128-, 256- and 512-qubit words. + let narrow = evolve::(64, 6, 6, false); + let wide = evolve::(200, 6, 6, false); + let widest = evolve::(400, 6, 6, false); + assert!(narrow.abs() > 1e-6, "test observable is trivially zero"); + assert_eq!(narrow.to_bits(), wide.to_bits(), "{narrow} vs {wide}"); + assert_eq!(narrow.to_bits(), widest.to_bits(), "{narrow} vs {widest}"); + + // The Kossakowski accumulation sums through a hash map, whose iteration + // order follows the word's hash — and a wider word hashes differently. + // Same physics, so agreement to rounding. + let narrow = evolve::(64, 6, 6, true); + for wide in [ + evolve::(200, 6, 6, true), + evolve::(400, 6, 6, true), + ] { + assert!( + (narrow - wide).abs() <= 1e-13 * narrow.abs(), + "{narrow} vs {wide}" + ); + } +} + +#[test] +fn same_width_different_register_size_agrees() { + // Within one width, the spectator qubits must not touch the answer. + assert_eq!( + evolve::(130, 6, 5, false).to_bits(), + evolve::(256, 6, 5, false).to_bits() + ); +} + +#[test] +fn orbit_rep_step_is_width_independent() { + // The momentum-orbit path (canonicalization, masked-shift generators, + // character table) on the same ring stored in 128- and 512-qubit words. + // Untruncated, so both widths hold the same reps; basis *order* follows + // hash-map iteration (width-dependent), so coefficients agree to rounding. + for (n, k, steps) in [(12, 0, 3), (12, 1, 3), (100, 0, 2)] { + let narrow = evolve_orbit::(n, k, steps); + let wide = evolve_orbit::(n, k, steps); + assert!(narrow.len() > 10, "basis did not grow"); + assert_eq!(narrow.len(), wide.len(), "n={n}, k={k}"); + for ((wa, ca), (wb, cb)) in narrow.iter().zip(&wide) { + assert_eq!(wa, wb, "n={n}, k={k}: rep sets differ"); + assert!((ca - cb).norm() <= 1e-12, "n={n}, k={k}: {ca} vs {cb}"); + } + } +} + +#[test] +fn orbit_rep_step_runs_beyond_128_qubits() { + let reps = evolve_orbit::(130, 0, 2); + assert!(reps.len() > 10); +} diff --git a/crates/ppvm-python-native/src/lindblad.rs b/crates/ppvm-python-native/src/lindblad.rs index f7cdd6bf2..796fcb995 100644 --- a/crates/ppvm-python-native/src/lindblad.rs +++ b/crates/ppvm-python-native/src/lindblad.rs @@ -9,12 +9,21 @@ //! only for the Python boundary: decoding the `(N, n_qubits)` numpy uint8 //! arrays into [`ppvm_lindblad::Word`] vectors, and re-encoding outputs //! back into numpy. +//! +//! The core crate is const-generic in the Pauli-word width. [`AnySpec`] +//! holds a spec at the narrowest of the 128/256/512-qubit widths that fits +//! `n_qubits` (picked once, at construction), and every method dispatches +//! on it with [`with_spec!`]. Registers of at most 128 qubits use exactly +//! the historical word layout, so nothing changes for them. use std::collections::HashMap; use num::Complex; use numpy::{Complex64, IntoPyArray, PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2}; -use ppvm_lindblad::{JumpInput, LindbladSpec as CoreSpec, Word, word_from_codes}; +use ppvm_lindblad::{ + JumpInput, LindbladSpec as CoreSpec, MAX_SUPPORTED_QUBITS, WIDTHS, Word, chunks_for, + word_from_codes, +}; use pyo3::{exceptions::PyValueError, prelude::*}; type PyPauliMap<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); @@ -32,8 +41,8 @@ pub(crate) fn map_err(e: ppvm_lindblad::Error) -> PyErr { /// Reject a basis that contains the same Pauli word at two distinct rows. /// Duplicate rows would silently overwrite each other in the generator's /// row-index map and produce an incorrect sparse matrix. -fn assert_basis_unique(basis: &[Word]) -> PyResult<()> { - let mut seen: HashMap<&Word, usize> = HashMap::with_capacity(basis.len()); +fn assert_basis_unique(basis: &[Word]) -> PyResult<()> { + let mut seen: HashMap<&Word, usize> = HashMap::with_capacity(basis.len()); for (i, w) in basis.iter().enumerate() { if let Some(prev) = seen.insert(w, i) { return Err(PyValueError::new_err(format!( @@ -58,20 +67,69 @@ fn to_lincomb(terms: Vec<(String, f64, f64)>) -> Vec<(String, Complex)> { } /// Pack `Vec<(Word, f64)>` into the standard PyO3 return shape. -fn pack_pauli_map<'py>( +fn pack_pauli_map<'py, const C: usize>( py: Python<'py>, - pairs: Vec<(Word, f64)>, + pairs: Vec<(Word, f64)>, n_qubits: usize, ) -> PyResult> { - let (words, coeffs): (Vec, Vec) = pairs.into_iter().unzip(); + let (words, coeffs): (Vec>, Vec) = pairs.into_iter().unzip(); let basis_arr = encode_basis(py, &words, n_qubits)?; Ok((basis_arr, coeffs.into_pyarray(py))) } +/// A [`CoreSpec`] at one of the [`WIDTHS`] (128, 256 or 512 qubits). +enum AnySpec { + W128(CoreSpec<{ WIDTHS[0] }>), + W256(CoreSpec<{ WIDTHS[1] }>), + W512(CoreSpec<{ WIDTHS[2] }>), +} + +/// Run `$body` with `$spec` bound to the concrete-width spec inside +/// `$inner: &AnySpec` and `$C` to its chunk count. +macro_rules! with_spec { + ($inner:expr, $spec:ident, $C:ident => $body:expr) => { + match $inner { + AnySpec::W128($spec) => { + #[allow(dead_code)] + const $C: usize = WIDTHS[0]; + $body + } + AnySpec::W256($spec) => { + #[allow(dead_code)] + const $C: usize = WIDTHS[1]; + $body + } + AnySpec::W512($spec) => { + #[allow(dead_code)] + const $C: usize = WIDTHS[2]; + $body + } + } + }; +} + +/// Kossakowski operators `A_n` (Pauli lincombs) and the pair matrix `K`. +type Kossakowski<'a> = (&'a [Vec<(String, Complex)>], &'a [Vec>]); + +/// Build the core spec at width `C` and attach the optional Kossakowski +/// dissipator. +fn build_spec( + n_qubits: usize, + h: &[(String, f64)], + jumps: &[JumpInput], + koss: Option>, +) -> PyResult> { + let mut spec = CoreSpec::::new(n_qubits, h, jumps).map_err(map_err)?; + if let Some((ops, k)) = koss { + spec.add_kossakowski(ops, k).map_err(map_err)?; + } + Ok(spec) +} + /// PyO3 facade exposing [`ppvm_lindblad::LindbladSpec`] to Python. #[pyclass] pub struct LindbladSpec { - inner: CoreSpec, + inner: AnySpec, } #[pymethods] @@ -117,36 +175,49 @@ impl LindbladSpec { rate, }) .collect(); - let mut inner = CoreSpec::new(n_qubits, &h, &jumps).map_err(map_err)?; - if !kossakowski_ops.is_empty() || !kossakowski_k.is_empty() { - let ops: Vec)>> = - kossakowski_ops.into_iter().map(to_lincomb).collect(); - let k: Vec>> = kossakowski_k - .into_iter() - .map(|row| { - row.into_iter() - .map(|(re, im)| Complex::new(re, im)) - .collect() - }) - .collect(); - inner.add_kossakowski(&ops, &k).map_err(map_err)?; - } + let ops: Vec)>> = + kossakowski_ops.into_iter().map(to_lincomb).collect(); + let k: Vec>> = kossakowski_k + .into_iter() + .map(|row| { + row.into_iter() + .map(|(re, im)| Complex::new(re, im)) + .collect() + }) + .collect(); + let koss = (!ops.is_empty() || !k.is_empty()).then_some((&ops[..], &k[..])); + let inner = match chunks_for(n_qubits) { + Some(c) if c == WIDTHS[0] => AnySpec::W128(build_spec(n_qubits, &h, &jumps, koss)?), + Some(c) if c == WIDTHS[1] => AnySpec::W256(build_spec(n_qubits, &h, &jumps, koss)?), + Some(_) => AnySpec::W512(build_spec(n_qubits, &h, &jumps, koss)?), + None => { + return Err(PyValueError::new_err(format!( + "LindbladSpec supports n_qubits ≤ {MAX_SUPPORTED_QUBITS}; got {n_qubits}" + ))); + } + }; Ok(Self { inner }) } #[getter] fn n_qubits(&self) -> usize { - self.inner.n_qubits() + with_spec!(&self.inner, inner, C => { + inner.n_qubits() + }) } #[getter] fn num_h_terms(&self) -> usize { - self.inner.num_h_terms() + with_spec!(&self.inner, inner, C => { + inner.num_h_terms() + }) } #[getter] fn num_jump_terms(&self) -> usize { - self.inner.num_jump_terms() + with_spec!(&self.inner, inner, C => { + inner.num_jump_terms() + }) } /// Apply `L*` to a single Pauli string `p`. @@ -155,10 +226,12 @@ impl LindbladSpec { py: Python<'py>, p: PyReadonlyArray1<'py, u8>, ) -> PyResult> { - let p_slice = p.as_slice()?; - let p_word = word_from_codes(p_slice).map_err(map_err)?; - let pairs = self.inner.action(&p_word); - pack_pauli_map(py, pairs, self.inner.n_qubits()) + with_spec!(&self.inner, inner, C => { + let p_slice = p.as_slice()?; + let p_word = word_from_codes::(p_slice).map_err(map_err)?; + let pairs = inner.action(&p_word); + pack_pauli_map(py, pairs, inner.n_qubits()) + }) } /// Off-basis component of `L*( Σ_j coeffs[j] · basis[j] )`. @@ -170,22 +243,23 @@ impl LindbladSpec { coeffs: PyReadonlyArray1<'py, f64>, protected: Option>, ) -> PyResult> { - let n_q = self.inner.n_qubits(); - let basis_view = basis.as_array(); - let basis_words = decode_basis(&basis_view, n_q)?; - let coeffs_slice = coeffs.as_slice()?; - check_coeffs_len(coeffs_slice.len(), basis_words.len())?; - let protected_words: Vec = if let Some(ref prot) = protected { - let pv = prot.as_array(); - decode_basis(&pv, n_q)? - } else { - Vec::new() - }; - let pairs = self - .inner - .leakage(&basis_words, coeffs_slice, &protected_words) - .map_err(map_err)?; - pack_pauli_map(py, pairs, n_q) + with_spec!(&self.inner, inner, C => { + let n_q = inner.n_qubits(); + let basis_view = basis.as_array(); + let basis_words = decode_basis::(&basis_view, n_q)?; + let coeffs_slice = coeffs.as_slice()?; + check_coeffs_len(coeffs_slice.len(), basis_words.len())?; + let protected_words: Vec> = if let Some(ref prot) = protected { + let pv = prot.as_array(); + decode_basis::(&pv, n_q)? + } else { + Vec::new() + }; + let pairs = inner + .leakage(&basis_words, coeffs_slice, &protected_words) + .map_err(map_err)?; + pack_pauli_map(py, pairs, n_q) + }) } /// One predictor-corrector adaptive step. @@ -226,36 +300,38 @@ impl LindbladSpec { admit_basis: Option, tau_add: Option, ) -> PyResult> { - let n_q = self.inner.n_qubits(); - let basis_view = basis.as_array(); - let mut basis_words = decode_basis(&basis_view, n_q)?; - assert_basis_unique(&basis_words)?; - let mut coeffs_vec = coeffs.as_slice()?.to_vec(); - check_coeffs_len(coeffs_vec.len(), basis_words.len())?; - let protected_words: Vec = if let Some(ref p) = protected { - decode_basis(&p.as_array(), n_q)? - } else { - Vec::new() - }; - self.inner - .pc_step( - &mut basis_words, - &mut coeffs_vec, - dt, - &protected_words, - &ppvm_lindblad::PcStepConfig { - max_basis, - admit_basis, - drop_tol, - tau_add, - num_threads, - }, - ) - .map_err(map_err)?; + with_spec!(&self.inner, inner, C => { + let n_q = inner.n_qubits(); + let basis_view = basis.as_array(); + let mut basis_words = decode_basis::(&basis_view, n_q)?; + assert_basis_unique(&basis_words)?; + let mut coeffs_vec = coeffs.as_slice()?.to_vec(); + check_coeffs_len(coeffs_vec.len(), basis_words.len())?; + let protected_words: Vec> = if let Some(ref p) = protected { + decode_basis::(&p.as_array(), n_q)? + } else { + Vec::new() + }; + inner + .pc_step( + &mut basis_words, + &mut coeffs_vec, + dt, + &protected_words, + &ppvm_lindblad::PcStepConfig { + max_basis, + admit_basis, + drop_tol, + tau_add, + num_threads, + }, + ) + .map_err(map_err)?; - // Pack output. Basis may have grown; coeffs has the same new length. - let pairs: Vec<(Word, f64)> = basis_words.into_iter().zip(coeffs_vec).collect(); - pack_pauli_map(py, pairs, n_q) + // Pack output. Basis may have grown; coeffs has the same new length. + let pairs: Vec<(Word, f64)> = basis_words.into_iter().zip(coeffs_vec).collect(); + pack_pauli_map(py, pairs, n_q) + }) } /// Same as [`Self::pc_step`] but also returns a dict mapping phase @@ -282,44 +358,45 @@ impl LindbladSpec { admit_basis: Option, tau_add: Option, ) -> PyResult<(PyPauliMap<'py>, Bound<'py, pyo3::types::PyDict>)> { - let n_q = self.inner.n_qubits(); - let basis_view = basis.as_array(); - let mut basis_words = decode_basis(&basis_view, n_q)?; - assert_basis_unique(&basis_words)?; - let mut coeffs_vec = coeffs.as_slice()?.to_vec(); - check_coeffs_len(coeffs_vec.len(), basis_words.len())?; - let protected_words: Vec = if let Some(ref p) = protected { - decode_basis(&p.as_array(), n_q)? - } else { - Vec::new() - }; - let timings = self - .inner - .pc_step_timed( - &mut basis_words, - &mut coeffs_vec, - dt, - &protected_words, - &ppvm_lindblad::PcStepConfig { - max_basis, - admit_basis, - drop_tol, - tau_add, - num_threads, - }, - ) - .map_err(map_err)?; + with_spec!(&self.inner, inner, C => { + let n_q = inner.n_qubits(); + let basis_view = basis.as_array(); + let mut basis_words = decode_basis::(&basis_view, n_q)?; + assert_basis_unique(&basis_words)?; + let mut coeffs_vec = coeffs.as_slice()?.to_vec(); + check_coeffs_len(coeffs_vec.len(), basis_words.len())?; + let protected_words: Vec> = if let Some(ref p) = protected { + decode_basis::(&p.as_array(), n_q)? + } else { + Vec::new() + }; + let timings = inner + .pc_step_timed( + &mut basis_words, + &mut coeffs_vec, + dt, + &protected_words, + &ppvm_lindblad::PcStepConfig { + max_basis, + admit_basis, + drop_tol, + tau_add, + num_threads, + }, + ) + .map_err(map_err)?; - let pairs: Vec<(Word, f64)> = basis_words.into_iter().zip(coeffs_vec).collect(); - let map = pack_pauli_map(py, pairs, n_q)?; - let d = pyo3::types::PyDict::new(py); - d.set_item("leakage1_us", timings.leakage1_us)?; - d.set_item("expand1_us", timings.expand1_us)?; - d.set_item("expm1_us", timings.expm1_us)?; - d.set_item("leakage2_us", timings.leakage2_us)?; - d.set_item("expand2_us", timings.expand2_us)?; - d.set_item("expm2_us", timings.expm2_us)?; - Ok((map, d)) + let pairs: Vec<(Word, f64)> = basis_words.into_iter().zip(coeffs_vec).collect(); + let map = pack_pauli_map(py, pairs, n_q)?; + let d = pyo3::types::PyDict::new(py); + d.set_item("leakage1_us", timings.leakage1_us)?; + d.set_item("expand1_us", timings.expand1_us)?; + d.set_item("expm1_us", timings.expm1_us)?; + d.set_item("leakage2_us", timings.leakage2_us)?; + d.set_item("expand2_us", timings.expand2_us)?; + d.set_item("expm2_us", timings.expm2_us)?; + Ok((map, d)) + }) } /// Per-step orbit-rep predictor-corrector evolution under @@ -370,56 +447,58 @@ impl LindbladSpec { tau_add: Option, num_threads: Option, ) -> PyResult> { - use num::Complex; - use ppvm_lindblad::{Sector, canonicalize_basis_to_rep}; + with_spec!(&self.inner, inner, C => { + use num::Complex; + use ppvm_lindblad::{Sector, canonicalize_basis_to_rep}; - let n_q = self.inner.n_qubits(); - let basis_view = basis.as_array(); - let mut basis_words = decode_basis(&basis_view, n_q)?; - let coeffs_slice = coeffs.as_slice()?; - check_coeffs_len(coeffs_slice.len(), basis_words.len())?; - let mut coeffs_vec: Vec> = coeffs_slice - .iter() - .map(|c| Complex::new(c.re, c.im)) - .collect(); - let protected_words: Vec = if let Some(ref p) = protected { - decode_basis(&p.as_array(), n_q)? - } else { - Vec::new() - }; - let k_slice = momentum.as_slice()?; - check_momentum_len(k_slice.len(), group.core().n_generators())?; - check_group_qubits(n_q, group.core().n_qubits())?; - if canonicalize_first { - canonicalize_basis_to_rep(&mut basis_words, group.core()); - } - // Canonicalization can collapse several input rows onto one rep, - // and the step indexes the basis by Pauli word — so uniqueness is - // checked after the rewrite, not before. - assert_basis_unique(&basis_words)?; - self.inner - .pc_step_orbit_rep( - &mut basis_words, - &mut coeffs_vec, - dt, - &protected_words, - &Sector::new(group.core(), k_slice), - &ppvm_lindblad::PcStepConfig { - max_basis, - admit_basis, - drop_tol, - tau_add, - num_threads, - }, - ) - .map_err(map_err)?; + let n_q = inner.n_qubits(); + let basis_view = basis.as_array(); + let mut basis_words = decode_basis::(&basis_view, n_q)?; + let coeffs_slice = coeffs.as_slice()?; + check_coeffs_len(coeffs_slice.len(), basis_words.len())?; + let mut coeffs_vec: Vec> = coeffs_slice + .iter() + .map(|c| Complex::new(c.re, c.im)) + .collect(); + let protected_words: Vec> = if let Some(ref p) = protected { + decode_basis::(&p.as_array(), n_q)? + } else { + Vec::new() + }; + let k_slice = momentum.as_slice()?; + check_momentum_len(k_slice.len(), group.core().n_generators())?; + check_group_qubits(n_q, group.core().n_qubits())?; + if canonicalize_first { + canonicalize_basis_to_rep(&mut basis_words, group.core()); + } + // Canonicalization can collapse several input rows onto one rep, + // and the step indexes the basis by Pauli word — so uniqueness is + // checked after the rewrite, not before. + assert_basis_unique(&basis_words)?; + inner + .pc_step_orbit_rep( + &mut basis_words, + &mut coeffs_vec, + dt, + &protected_words, + &Sector::new(group.core(), k_slice), + &ppvm_lindblad::PcStepConfig { + max_basis, + admit_basis, + drop_tol, + tau_add, + num_threads, + }, + ) + .map_err(map_err)?; - let out_coeffs: Vec = coeffs_vec - .iter() - .map(|c| Complex64::new(c.re, c.im)) - .collect(); - let basis_arr = encode_basis(py, &basis_words, n_q)?; - Ok((basis_arr, out_coeffs.into_pyarray(py))) + let out_coeffs: Vec = coeffs_vec + .iter() + .map(|c| Complex64::new(c.re, c.im)) + .collect(); + let basis_arr = encode_basis(py, &basis_words, n_q)?; + Ok((basis_arr, out_coeffs.into_pyarray(py))) + }) } /// Sparse generator matrix in COO form: `(rows, cols, vals)`. @@ -428,24 +507,26 @@ impl LindbladSpec { py: Python<'py>, basis: PyReadonlyArray2<'py, u8>, ) -> PyResult> { - let n_q = self.inner.n_qubits(); - let basis_view = basis.as_array(); - let basis_words = decode_basis(&basis_view, n_q)?; - assert_basis_unique(&basis_words)?; - let triplets = self.inner.generator(&basis_words); - let total = triplets.len(); - let mut rows = Vec::with_capacity(total); - let mut cols = Vec::with_capacity(total); - let mut vals = Vec::with_capacity(total); - for (r, c, v) in triplets { - rows.push(r as u64); - cols.push(c as u64); - vals.push(v); - } - Ok(( - rows.into_pyarray(py), - cols.into_pyarray(py), - vals.into_pyarray(py), - )) + with_spec!(&self.inner, inner, C => { + let n_q = inner.n_qubits(); + let basis_view = basis.as_array(); + let basis_words = decode_basis::(&basis_view, n_q)?; + assert_basis_unique(&basis_words)?; + let triplets = inner.generator(&basis_words); + let total = triplets.len(); + let mut rows = Vec::with_capacity(total); + let mut cols = Vec::with_capacity(total); + let mut vals = Vec::with_capacity(total); + for (r, c, v) in triplets { + rows.push(r as u64); + cols.push(c as u64); + vals.push(v); + } + Ok(( + rows.into_pyarray(py), + cols.into_pyarray(py), + vals.into_pyarray(py), + )) + }) } } diff --git a/crates/ppvm-python-native/src/pauli_arr.rs b/crates/ppvm-python-native/src/pauli_arr.rs index 45580d638..c7d08b0e9 100644 --- a/crates/ppvm-python-native/src/pauli_arr.rs +++ b/crates/ppvm-python-native/src/pauli_arr.rs @@ -15,11 +15,42 @@ use pyo3::{exceptions::PyValueError, prelude::*}; use crate::lindblad::map_err; +/// Run `$body` with `$C` bound to the narrowest Pauli-word chunk count +/// that holds `$n` qubits (see [`ppvm_lindblad::chunks_for`]); a +/// `ValueError` above [`ppvm_lindblad::MAX_SUPPORTED_QUBITS`]. +/// +/// Mirrors the width `LindbladSpec` picks, so `n <= 128` keeps the +/// original 128-qubit word layout. +macro_rules! with_width { + ($n:expr, $C:ident => $body:expr) => {{ + let n: usize = $n; + match ppvm_lindblad::chunks_for(n) { + Some(c) if c == ppvm_lindblad::WIDTHS[0] => { + const $C: usize = ppvm_lindblad::WIDTHS[0]; + $body + } + Some(c) if c == ppvm_lindblad::WIDTHS[1] => { + const $C: usize = ppvm_lindblad::WIDTHS[1]; + $body + } + Some(_) => { + const $C: usize = ppvm_lindblad::WIDTHS[2]; + $body + } + None => Err(pyo3::exceptions::PyValueError::new_err(format!( + "Pauli words support n_qubits ≤ {}; got {n}", + ppvm_lindblad::MAX_SUPPORTED_QUBITS + ))), + } + }}; +} +pub(crate) use with_width; + /// Decode a `(N, n_qubits)` uint8 ndarray view into `N` packed [`Word`]s. -pub(crate) fn decode_basis( +pub(crate) fn decode_basis( view: &numpy::ndarray::ArrayView2, n_qubits: usize, -) -> PyResult> { +) -> PyResult>> { let n_basis = view.shape()[0]; let n_cols = view.shape()[1]; if n_cols != n_qubits { @@ -40,9 +71,9 @@ pub(crate) fn decode_basis( } /// Encode packed [`Word`]s back into an `(M, n_qubits)` uint8 array. -pub(crate) fn encode_basis<'py>( +pub(crate) fn encode_basis<'py, const C: usize>( py: Python<'py>, - words: &[Word], + words: &[Word], n_qubits: usize, ) -> PyResult>> { let m = words.len(); diff --git a/crates/ppvm-python-native/src/symmetry.rs b/crates/ppvm-python-native/src/symmetry.rs index cd3dfa32b..4d8129a8a 100644 --- a/crates/ppvm-python-native/src/symmetry.rs +++ b/crates/ppvm-python-native/src/symmetry.rs @@ -17,7 +17,7 @@ use ppvm_pauli_sum::symmetry as core_sym; use pyo3::{exceptions::PyValueError, prelude::*}; use crate::pauli_arr::{ - check_coeffs_len, check_group_width, check_momentum_len, decode_basis, encode_basis, + check_coeffs_len, check_group_width, check_momentum_len, decode_basis, encode_basis, with_width, }; type PyPauliMap<'py> = (Bound<'py, PyArray2>, Bound<'py, PyArray1>); @@ -158,10 +158,13 @@ impl TranslationGroup { self.inner.n_qubits() ))); } - let w = word_from_codes(codes).map_err(|e| PyValueError::new_err(e.to_string()))?; - let canon = self.inner.canonicalize(&w); let mut out = vec![0u8; codes.len()]; - codes_from_word(&canon, &mut out); + with_width!(codes.len(), C => { + let w = word_from_codes::(codes).map_err(|e| PyValueError::new_err(e.to_string()))?; + let canon = self.inner.canonicalize(&w); + codes_from_word(&canon, &mut out); + PyResult::Ok(()) + })?; Ok(out.into_pyarray(py)) } } @@ -195,25 +198,27 @@ pub fn canonicalize_basis_arr_complex<'py>( check_coeffs_len(coeffs_slice.len(), basis_view.shape()[0])?; let k_slice = momentum.as_slice()?; check_momentum_len(k_slice.len(), group.inner.n_generators())?; - let mut basis_words = decode_basis(&basis_view, n_q)?; - let mut coeffs_vec: Vec> = coeffs_slice - .iter() - .map(|c| Complex::new(c.re, c.im)) - .collect(); + with_width!(n_q, C => { + let mut basis_words = decode_basis::(&basis_view, n_q)?; + let mut coeffs_vec: Vec> = coeffs_slice + .iter() + .map(|c| Complex::new(c.re, c.im)) + .collect(); - core_sym::canonicalize_pauli_sum_complex( - &mut basis_words, - &mut coeffs_vec, - &group.inner, - k_slice, - ); + core_sym::canonicalize_pauli_sum_complex( + &mut basis_words, + &mut coeffs_vec, + &group.inner, + k_slice, + ); - let out_coeffs: Vec = coeffs_vec - .iter() - .map(|c| Complex64::new(c.re, c.im)) - .collect(); - let basis_arr = encode_basis(py, &basis_words, n_q)?; - Ok((basis_arr, out_coeffs.into_pyarray(py))) + let out_coeffs: Vec = coeffs_vec + .iter() + .map(|c| Complex64::new(c.re, c.im)) + .collect(); + let basis_arr = encode_basis(py, &basis_words, n_q)?; + Ok((basis_arr, out_coeffs.into_pyarray(py))) + }) } /// Verify that a `(basis_arr, complex_coeffs)` Pauli sum lies in the @@ -238,13 +243,15 @@ pub fn check_momentum_sector_arr<'py>( check_coeffs_len(coeffs_slice.len(), basis_view.shape()[0])?; let k_slice = momentum.as_slice()?; check_momentum_len(k_slice.len(), group.inner.n_generators())?; - let basis_words = decode_basis(&basis_view, n_q)?; - let coeffs_vec: Vec> = coeffs_slice - .iter() - .map(|c| Complex::new(c.re, c.im)) - .collect(); - core_sym::check_momentum_sector(&basis_words, &coeffs_vec, &group.inner, k_slice, tol) - .map_err(|e| PyValueError::new_err(format!("{e}"))) + with_width!(n_q, C => { + let basis_words = decode_basis::(&basis_view, n_q)?; + let coeffs_vec: Vec> = coeffs_slice + .iter() + .map(|c| Complex::new(c.re, c.im)) + .collect(); + core_sym::check_momentum_sector(&basis_words, &coeffs_vec, &group.inner, k_slice, tol) + .map_err(|e| PyValueError::new_err(format!("{e}"))) + }) } /// Merge a `(basis_arr, coeffs)` Pauli sum (the representation used by @@ -272,11 +279,13 @@ pub fn canonicalize_basis_arr<'py>( let coeffs_slice = coeffs.as_slice()?; check_coeffs_len(coeffs_slice.len(), basis_view.shape()[0])?; - let mut basis_words = decode_basis(&basis_view, n_q)?; - let mut coeffs_vec = coeffs_slice.to_vec(); + with_width!(n_q, C => { + let mut basis_words = decode_basis::(&basis_view, n_q)?; + let mut coeffs_vec = coeffs_slice.to_vec(); - core_sym::canonicalize_pauli_sum(&mut basis_words, &mut coeffs_vec, &group.inner); + core_sym::canonicalize_pauli_sum(&mut basis_words, &mut coeffs_vec, &group.inner); - let basis_arr = encode_basis(py, &basis_words, n_q)?; - Ok((basis_arr, coeffs_vec.into_pyarray(py))) + let basis_arr = encode_basis(py, &basis_words, n_q)?; + Ok((basis_arr, coeffs_vec.into_pyarray(py))) + }) } diff --git a/ppvm-python/src/ppvm/lindblad.py b/ppvm-python/src/ppvm/lindblad.py index 8dc3f9455..501eb7d12 100644 --- a/ppvm-python/src/ppvm/lindblad.py +++ b/ppvm-python/src/ppvm/lindblad.py @@ -143,7 +143,10 @@ class Lindbladian: Parameters ---------- n_qubits: - Number of qubits. + Number of qubits, at most 512. The Pauli word width is chosen + automatically: up to 128 qubits uses the original 128-bit layout, + then 256- and 512-bit words. Narrow problems are unaffected by the + wider layouts being available. h_terms: Iterable of ``(pauli_string, coefficient)`` pairs for the Hermitian Hamiltonian ``H = Σ c_i P_i``. Each ``pauli_string`` is diff --git a/ppvm-python/test/lindblad/test_word_width.py b/ppvm-python/test/lindblad/test_word_width.py new file mode 100644 index 000000000..4520bf4e3 --- /dev/null +++ b/ppvm-python/test/lindblad/test_word_width.py @@ -0,0 +1,86 @@ +# SPDX-FileCopyrightText: 2026 The PPVM Authors +# SPDX-License-Identifier: Apache-2.0 +"""Registers wider than 128 qubits (variable Pauli-word width).""" + +import numpy as np +import pytest + +from ppvm import Lindbladian +from ppvm._core import TranslationGroup, canonicalize_basis_arr_complex +from ppvm.lindblad import _basis_to_codes + + +def _word(n, ops): + s = ["I"] * n + for q, p in ops: + s[q] = p + return "".join(s) + + +def _ring(n): + terms = [] + for i in range(n): + j = (i + 1) % n + terms += [ + (_word(n, [(i, "X"), (j, "X")]), 1.0), + (_word(n, [(i, "Y"), (j, "Y")]), 1.0), + (_word(n, [(i, "Z"), (j, "Z")]), 0.5), + (_word(n, [(i, "Z")]), 0.3), + ] + return terms + + +def _orbit_run(n, steps=2): + lind = Lindbladian(n, _ring(n), []) + group = TranslationGroup.chain_1d(n) + mom = np.array([0], dtype=np.int32) + seed = _basis_to_codes([_word(n, [(q, "X")]) for q in range(n)], n) + basis, co = canonicalize_basis_arr_complex(seed, np.ones(n, dtype=np.complex128), group, mom) + for _ in range(steps): + basis, co = lind.pc_step_orbit_rep( + basis, co, 0.1, 10**6, group=group, momentum=mom, drop_tol=0.0 + ) + return basis, co + + +@pytest.mark.parametrize("n", [130, 256, 300, 512]) +def test_wide_lindbladian_constructs_and_steps(n): + # 130 qubits used to fail with "LindbladSpec supports n_qubits ≤ 128". + lind = Lindbladian( + n, + [(_word(n, [(0, "Z"), (n - 1, "Z")]), 0.7), (_word(n, [(n - 1, "X")]), 1.3)], + [(_word(n, [(n - 1, "Z")]), 0.05)], + ) + assert lind.n_qubits == n + basis, _ = lind.pc_step([_word(n, [(n - 1, "Z")])], np.array([1.0]), 0.05, 1000) + assert len(basis) > 1 + assert all(len(b) == n for b in basis) + # The top qubit is live: something acts on it. + assert any(b[n - 1] == "Y" for b in basis) + + +def test_more_than_512_qubits_is_rejected(): + with pytest.raises(ValueError, match="512"): + Lindbladian(513, [(_word(513, [(0, "Z")]), 1.0)], []) + + +def test_orbit_rep_step_beyond_128_qubits(): + basis, co = _orbit_run(130) + assert basis.shape[1] == 130 + assert len(basis) > 5 + assert np.all(np.isfinite(co)) + + +def test_orbit_rep_step_matches_across_the_width_boundary(): + # Rings of 120 and 136 sites use 128- and 256-qubit words. Compare the + # k=0 autocorrelation of M_x, which is ring-size independent here. + def autocorr(n): + basis, co = _orbit_run(n, steps=3) + # The single-X rep, whatever position the canonical form puts it at. + m = np.where(((basis == 1).sum(axis=1) == 1) & ((basis != 0).sum(axis=1) == 1))[0] + assert m.size == 1 + return co[m[0]].real + + # With nearest-neighbour H and 3 short steps the operator front is far + # smaller than either ring, so the value is ring-size independent. + assert autocorr(120) == pytest.approx(autocorr(136), rel=1e-12, abs=1e-14) From 78a4b6db61c9bd3913122860561de1788f23b5c0 Mon Sep 17 00:00:00 2001 From: Alexander Schuckert Date: Mon, 28 Sep 2026 16:44:27 +0300 Subject: [PATCH 6/7] perf(lindblad): store the expm generator cache as blocked flat CSC (#229) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Stacked on #222 (`split/5-kossakowski`). ## Summary `build_mf_cols` / `build_orbit_rep_cols` returned one `Vec` per column, each reserved at the full L* output count although only in-basis entries are kept. The expm generator cache is now a `BlockCsc`: blocks of up to 4096 columns, each an exactly sized `offsets`/`rows`/`vals` triple (structure of arrays, 12 B per real entry instead of 16 B). `CscOp::dot` keeps the same per-thread column partition, so results are bit-identical. Single commit, `crates/ppvm-lindblad/src/mf_expm.rs` only. ## Measurements `pc_step` on the two-leg XY ladder with a local probe (not part of this PR): L = 41 rungs (N = 82), `max_basis = 2^18`, `admit_basis = 3·2^18`, dt = 0.1, 25 steps, 4 threads, `drop_tol = 0`. Apple M4 MacBook Air (fanless, so wall times are only comparable within a rep). Max RSS from `/usr/bin/time -l`; live-heap peak from a counting global allocator. | mimalloc, 3 interleaved reps | wall (s) | max RSS (MiB) | live-heap peak (MiB) | |---|---|---|---| | base (#222) | 53.7 / 53.5 / 62.5 | 484 / 484 / 483 | 441.7 | | this PR | 46.2 / 48.9 / 60.1 | 501 / 491 / 470 | 275.2 | - The live-heap peak drops by 38% and wall time by 4–14%, with each expm call about 20–30% faster. - On macOS + mimalloc the heap saving shows up only weakly in the process footprint: max RSS moves by −3% to +4%, and peak footprint (`/usr/bin/time -l`) drops by about 3–13%. With the macOS system allocator, RSS is not reproducible (768–1237 MiB for the same binary), so the effect cannot be resolved there. - This has not been measured on Linux/glibc, where fragmentation from the per-column allocations of the old layout may matter more. - The final coefficients are bit-identical to the base in every run. Tests: `cargo test -p ppvm-lindblad` passes (8/8) with this change stacked together with the skip-empty-hop2 PR; not yet run on this branch alone (CI will). 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Co-authored-by: David Plankensteiner --- crates/ppvm-lindblad/src/mf_expm.rs | 247 +++++++++++++++++----------- 1 file changed, 154 insertions(+), 93 deletions(-) diff --git a/crates/ppvm-lindblad/src/mf_expm.rs b/crates/ppvm-lindblad/src/mf_expm.rs index caca82c5a..2b99be096 100644 --- a/crates/ppvm-lindblad/src/mf_expm.rs +++ b/crates/ppvm-lindblad/src/mf_expm.rs @@ -35,12 +35,105 @@ use rayon::prelude::*; use std::iter::Sum; use std::ops::{AddAssign, Div, Mul, Sub}; -/// CSC columns of a cached in-basis action: `cols[c]` = `(row, coeff)`. -type Cols = Vec>; /// Per-column `(raw, diag)` for the `μ`/1-norm selection: `raw` bounds /// `Σ_r |M[r,c]|` from above and `diag = M[c,c]`. type PerCol = Vec<(f64, T)>; +/// Scratch buffers for [`LindbladSpec::compute_action_terms`]. +type ActionScratch = (Vec, Vec, FxHashMap, Complex>); + +/// Consecutive CSC columns stored flat: local column `j` holds +/// `rows[offsets[j]..offsets[j + 1]]` and the matching `vals`. +struct CscBlock { + offsets: Vec, + rows: Vec, + vals: Vec, +} + +/// Cached in-basis action in CSC form, stored as blocks of `block` columns. +/// +/// One exactly-sized allocation triple per block replaces one `Vec` per +/// column: at `|basis| ~ 10^6` the per-column layout reserved every `L*` +/// output (in- and out-of-basis) and left ~10^6 small allocations for the +/// system allocator to retain after the expm call. +pub(crate) struct BlockCsc { + blocks: Vec>, + block: usize, + dim: usize, +} + +impl BlockCsc { + /// Visit the columns `range` in order as `(col, rows, vals)`. + fn for_each_col(&self, range: std::ops::Range, mut f: impl FnMut(usize, &[u32], &[T])) { + let mut c = range.start; + while c < range.end { + let b = &self.blocks[c / self.block]; + let base = (c / self.block) * self.block; + let stop = range.end.min(base + b.offsets.len() - 1); + for j in (c - base)..(stop - base) { + let (lo, hi) = (b.offsets[j] as usize, b.offsets[j + 1] as usize); + f(base + j, &b.rows[lo..hi], &b.vals[lo..hi]); + } + c = stop; + } + } +} + +/// Build the [`BlockCsc`] cache and the per-column `(raw, diag)` data for +/// a `dim`-column generator in one parallel pass. `col(c, scratch, rows, +/// vals)` appends the in-basis entries of column `c` to `rows`/`vals` and +/// returns its `(raw, diag)`. +fn build_block_csc( + spec: &LindbladSpec, + dim: usize, + col: F, +) -> (BlockCsc, PerCol) +where + T: Copy + Send + Sync, + F: Fn(usize, &mut ActionScratch, &mut Vec, &mut Vec) -> (f64, T) + Sync, +{ + // ~16 blocks per thread for load balance, but never so small that the + // per-block allocations matter. + let block = dim + .div_ceil(16 * rayon::current_num_threads().max(1)) + .clamp(64, 4096); + let (blocks, per_col): (Vec>, Vec>) = (0..dim.div_ceil(block)) + .into_par_iter() + .map_init( + || { + let scratch: ActionScratch = ( + Vec::with_capacity(spec.n_qubits()), + Vec::with_capacity(128), + FxHashMap::with_capacity_and_hasher(128, FxBuildHasher::default()), + ); + (scratch, Vec::::new(), Vec::::new()) + }, + |(scratch, rows, vals), b| { + let cols = (b * block)..dim.min((b + 1) * block); + rows.clear(); + vals.clear(); + let mut offsets = Vec::with_capacity(cols.len() + 1); + let mut per_col = Vec::with_capacity(cols.len()); + offsets.push(0); + for c in cols { + per_col.push(col(c, scratch, rows, vals)); + offsets.push(u32::try_from(rows.len()).expect("CSC block exceeds u32 entries")); + } + // `to_vec` sizes the stored block exactly; the staging + // buffers are reused for the next block on this thread. + let blk = CscBlock { + offsets, + rows: rows.to_vec(), + vals: vals.to_vec(), + }; + (blk, per_col) + }, + ) + .unzip(); + let per_col = per_col.into_iter().flatten().collect(); + (BlockCsc { blocks, block, dim }, per_col) +} + /// Per-column in-basis action of the real generator `M`, plus the data the /// `(m, s)`/`μ` selection needs — all from ONE action pass over the basis. /// @@ -54,38 +147,24 @@ fn build_mf_cols( spec: &LindbladSpec, basis: &[Word], index: &FxHashMap, u32>, -) -> (Cols, PerCol) { - basis - .par_iter() - .map_init( - || { - ( - Vec::::with_capacity(spec.n_qubits()), - Vec::::with_capacity(128), - FxHashMap::, Complex>::with_capacity_and_hasher( - 128, - FxBuildHasher::default(), - ), - ) - }, - |(s1, s2, lm), p| { - let terms = spec.compute_action_terms(p, s1, s2, lm); - let mut out = Vec::with_capacity(terms.len()); - let mut raw = 0.0; - let mut diag = 0.0; - for (w, c) in terms.iter() { - raw += c.abs(); - if w == p { - diag = *c; - } - if let Some(&row) = index.get(w) { - out.push((row, *c)); - } - } - (out, (raw, diag)) - }, - ) - .unzip() +) -> (BlockCsc, PerCol) { + build_block_csc(spec, basis.len(), |c, (s1, s2, lm), rows, vals| { + let p = &basis[c]; + let terms = spec.compute_action_terms(p, s1, s2, lm); + let mut raw = 0.0; + let mut diag = 0.0; + for (w, v) in terms.iter() { + raw += v.abs(); + if w == p { + diag = *v; + } + if let Some(&row) = index.get(w) { + rows.push(row); + vals.push(*v); + } + } + (raw, diag) + }) } /// Per-column **phase-aware** action of the in-basis-restricted orbit-rep @@ -118,48 +197,33 @@ fn build_orbit_rep_cols( basis: &[Word], index: &FxHashMap, u32>, sector: &Sector<'_>, -) -> (Cols>, PerCol>) { - basis - .par_iter() - .enumerate() - .map_init( - || { - ( - Vec::::with_capacity(spec.n_qubits()), - Vec::::with_capacity(128), - FxHashMap::, Complex>::with_capacity_and_hasher( - 128, - FxBuildHasher::default(), - ), - ) - }, - |(s1, s2, lm), (c, r)| { - // A rep that cannot carry the sector has coefficient zero - // identically, so its column is empty. - let Some(orbit_in) = sector.orbit_size(r) else { - return (Vec::new(), (0.0, Complex::new(0.0, 0.0))); - }; - let terms = spec.compute_action_terms(r, s1, s2, lm); - let mut out = Vec::with_capacity(terms.len()); - let mut raw = 0.0; - let mut diag = Complex::new(0.0, 0.0); - for (q, v) in terms.iter() { - let Some((r_q, phase, orbit_out)) = sector.canonicalize_phase(q) else { - continue; - }; - if let Some(&row) = index.get(&r_q) { - let val = phase * *v * (orbit_in as f64 / orbit_out as f64); - raw += val.norm(); - if row as usize == c { - diag += val; - } - out.push((row, val)); - } +) -> (BlockCsc>, PerCol>) { + build_block_csc(spec, basis.len(), |c, (s1, s2, lm), rows, vals| { + let r = &basis[c]; + // A rep that cannot carry the sector has coefficient zero + // identically, so its column is empty. + let Some(orbit_in) = sector.orbit_size(r) else { + return (0.0, Complex::new(0.0, 0.0)); + }; + let terms = spec.compute_action_terms(r, s1, s2, lm); + let mut raw = 0.0; + let mut diag = Complex::new(0.0, 0.0); + for (q, v) in terms.iter() { + let Some((r_q, phase, orbit_out)) = sector.canonicalize_phase(q) else { + continue; + }; + if let Some(&row) = index.get(&r_q) { + let val = phase * *v * (orbit_in as f64 / orbit_out as f64); + raw += val.norm(); + if row as usize == c { + diag += val; } - (out, (raw, diag)) - }, - ) - .unzip() + rows.push(row); + vals.push(val); + } + } + (raw, diag) + }) } /// Borrowed CSC-style view of an in-basis-restricted generator `M`, backed @@ -172,8 +236,7 @@ fn build_orbit_rep_cols( /// impl for `&T`, so `ExpmOp::from_parts(op, ...)` accepts a `CscOp` by /// value while it keeps borrowing `cols`. pub(crate) struct CscOp<'a, T> { - pub(crate) cols: &'a [Vec<(u32, T)>], - pub(crate) dim: usize, + pub(crate) cols: &'a BlockCsc, } impl LinearOperator for CscOp<'_, T> @@ -188,7 +251,7 @@ where + Sync, { fn dim(&self) -> usize { - self.dim + self.cols.dim } fn parallel_hint(&self) -> bool { @@ -199,7 +262,7 @@ where } fn dot(&self, overwrite: bool, input: &[T], output: &mut [T]) -> Result<(), QuSpinError> { - let n = self.dim; + let n = self.cols.dim; if n == 0 { return Ok(()); } @@ -209,22 +272,20 @@ where // Parallelise over column chunks; each thread accumulates into a dense // local `y` of length `dim`, reading the cached action; the partials // are reduced into `output` sequentially at the end. - let partial_ys: Vec> = self - .cols - .par_chunks(chunk_size) - .enumerate() - .map(|(chunk_idx, chunk)| { - let c_offset = chunk_idx * chunk_size; + let partial_ys: Vec> = (0..n.div_ceil(chunk_size)) + .into_par_iter() + .map(|chunk_idx| { + let cols = (chunk_idx * chunk_size)..n.min((chunk_idx + 1) * chunk_size); let mut y_local = vec![T::zero(); n]; - for (c_local, col) in chunk.iter().enumerate() { - let xc = input[c_offset + c_local]; + self.cols.for_each_col(cols, |c, rows, vals| { + let xc = input[c]; if xc == T::zero() { - continue; + return; } - for &(row, val) in col.iter() { + for (&row, &val) in rows.iter().zip(vals) { y_local[row as usize] += val * xc; } - } + }); y_local }) .collect(); @@ -305,7 +366,7 @@ where /// `select` maps `‖dt·(M−μI)‖₁` to `(m*, s, backward-error tol)`; the two /// call sites differ only in that choice. fn expm_apply_cached( - cols: &Cols, + cols: &BlockCsc, per_col: &PerCol, dt: f64, coeffs: &[T], @@ -323,7 +384,7 @@ where + From + Sum, { - let n = cols.len(); + let n = cols.dim; let trace: T = per_col.iter().map(|(_, d)| *d).sum(); let mu = trace / n as f64; let onenorm = per_col @@ -333,7 +394,7 @@ where let (m_star, s, expm_tol) = select(dt.abs() * onenorm); let mut v = coeffs.to_vec(); - let op = CscOp { cols, dim: n }; + let op = CscOp { cols }; let expm = quspin_expm::ExpmOp::from_parts(op, T::from(dt), mu, s as usize, m_star as usize, expm_tol); expm.apply(ndarray::ArrayViewMut1::from(v.as_mut_slice())) From 2ab9b53abe3e829897bd67306a07edec7f06204e Mon Sep 17 00:00:00 2001 From: Alexander Schuckert Date: Mon, 28 Sep 2026 16:44:36 +0300 Subject: [PATCH 7/7] perf(lindblad): skip the second hop and corrector when hop 2 admits nothing (#230) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Stacked on #222 (`split/5-kossakowski`). ## Summary With `admit_basis` set, once the basis has reached `max_basis` the first leakage admission fills the whole room `admit − |basis|`. The second leakage call then ran with room = 0 and admitted nothing in every capped step, while still evaluating all its candidates, and the corrector repeated the predictor exactly on the same basis and input. `pc_step` and `pc_step_orbit_rep` now skip the second leakage pass when room = 0 and reuse the predicted state whenever the second hop admitted no string. Results are bit-identical; skipped phases report 0 in `PcStepTimings`. Default behaviour is otherwise unchanged: capped steps still carry no second-order admission, as before. Single commit, `crates/ppvm-lindblad/src/step.rs` only. ## Measurements `pc_step` on the two-leg XY ladder with a local probe (not part of this PR): L = 41 rungs (N = 82), `max_basis = 2^18`, `admit_basis = 3·2^18`, dt = 0.1, 25 steps, 4 threads, `drop_tol = 0`. Apple M4 MacBook Air (fanless, so wall times are only comparable within a rep). In the base, 21 of the 25 steps are capped and the second hop (leakage2 + expm2) accounts for 53% of the wall time (34.5 of 65.2 s). | | wall (s), per rep | |---|---| | base (#222), mimalloc | 53.7 / 53.5 / 62.5 | | this PR, mimalloc | 26.6 / 28.7 / 33.3 | | base (#222), system malloc | 65.6 / 79.1 | | this PR, system malloc | 32.4 / 38.4 | - The run is 1.9–2.0× faster end-to-end; memory is unchanged (live-heap peak 441.7 MiB in both). - The final coefficients are bit-identical to the base in every run. Tests: `cargo test -p ppvm-lindblad` passes (8/8) with this change stacked together with the blocked-CSC PR; not yet run on this branch alone (CI will). 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Co-authored-by: David Plankensteiner --- crates/ppvm-lindblad/src/step.rs | 74 ++++++++++++++++++++++---------- 1 file changed, 52 insertions(+), 22 deletions(-) diff --git a/crates/ppvm-lindblad/src/step.rs b/crates/ppvm-lindblad/src/step.rs index e5dce345e..debd9666f 100644 --- a/crates/ppvm-lindblad/src/step.rs +++ b/crates/ppvm-lindblad/src/step.rs @@ -58,6 +58,14 @@ impl LindbladSpec { /// rank cap) per [`PcStepConfig`]. Exact in `dt` within the working /// basis — the only error is basis truncation. /// + /// When the second hop admits no string (in particular when the first + /// hop already filled the admission room `admit_basis − |basis|`, the + /// usual case once the basis has reached `max_basis`), the corrector + /// would repeat the predictor exactly; the second leakage pass (if + /// `room = 0`) and the corrector exponential are then skipped, with + /// bit-identical results. The step's timings report 0 for skipped + /// phases. + /// /// `protected` words are never dropped. All tuning knobs live in `cfg`. pub fn pc_step( &self, @@ -149,23 +157,35 @@ impl LindbladSpec { let coeffs_predict = self.expm_step(basis, dt, coeffs, drop_tol); p.stop(&mut t.expm1_us); - // 3. Second-hop expansion from the predicted state. After leakage2 - // we no longer need `coeffs_predict`. Extend `coeffs` with zeros for - // any newly-added second-hop strings so it remains a valid input - // (pre-step state) for the corrector. - let p = Phase::start(timed); - let leak2 = self.leakage_with_prune(basis, &coeffs_predict, protected, admit, tau_add)?; - p.stop(&mut t.leakage2_us); - drop(coeffs_predict); + // 3. Second-hop expansion from the predicted state. Extend `coeffs` + // with zeros for any newly-added second-hop strings so it remains a + // valid input (pre-step state) for the corrector. Once the basis is + // full, the first hop usually fills the whole admission room; with + // `room = 0` the second hop can admit nothing, so its leakage pass + // is skipped. + let n_predict = basis.len(); + if admit > n_predict { + let p = Phase::start(timed); + let leak2 = + self.leakage_with_prune(basis, &coeffs_predict, protected, admit, tau_add)?; + p.stop(&mut t.leakage2_us); - let p = Phase::start(timed); - add_leakage_capped(basis, coeffs, leak2, admit); - p.stop(&mut t.expand2_us); + let p = Phase::start(timed); + add_leakage_capped(basis, coeffs, leak2, admit); + p.stop(&mut t.expand2_us); + } // 4. Corrector: redo from pre-step state on the doubly-enlarged basis. - let p = Phase::start(timed); - *coeffs = self.expm_step(basis, dt, coeffs, drop_tol); - p.stop(&mut t.expm2_us); + // If the second hop admitted nothing, the corrector would repeat the + // predictor's computation exactly, so the predicted state is kept. + if basis.len() == n_predict { + *coeffs = coeffs_predict; + } else { + drop(coeffs_predict); + let p = Phase::start(timed); + *coeffs = self.expm_step(basis, dt, coeffs, drop_tol); + p.stop(&mut t.expm2_us); + } // 5. Prune basis entries below `drop_tol` (protected words never dropped). prune_basis(basis, coeffs, drop_tol, protected); @@ -250,16 +270,26 @@ impl LindbladSpec { // across every matvec. let coeffs_predict = mf_expm::expm_apply_orbit_rep(self, basis, sector, dt, coeffs); - // 3. Second-hop leakage from the predicted state. - let mut leak2 = self.leakage_orbit_rep(basis, &coeffs_predict, protected, sector, admit)?; - drop(coeffs_predict); - if tau_add > 0.0 { - leak2.retain(|(_, c)| c.norm() > tau_add); + // 3. Second-hop leakage from the predicted state, skipped when the + // first hop left no admission room (see `pc_step_inner`). + let n_predict = basis.len(); + if admit > n_predict { + let mut leak2 = + self.leakage_orbit_rep(basis, &coeffs_predict, protected, sector, admit)?; + if tau_add > 0.0 { + leak2.retain(|(_, c)| c.norm() > tau_add); + } + add_leakage_capped(basis, coeffs, leak2, admit); } - add_leakage_capped(basis, coeffs, leak2, admit); - // 4. Corrector: redo from the pre-step state (the basis grew). - *coeffs = mf_expm::expm_apply_orbit_rep(self, basis, sector, dt, coeffs); + // 4. Corrector: redo from the pre-step state if the basis grew; + // otherwise it would reproduce the predictor exactly. + if basis.len() == n_predict { + *coeffs = coeffs_predict; + } else { + drop(coeffs_predict); + *coeffs = mf_expm::expm_apply_orbit_rep(self, basis, sector, dt, coeffs); + } // 5. Prune by magnitude, then rank-cap to max_basis. prune_basis(basis, coeffs, drop_tol, protected);