Use exact FCI traces for matrix exponentials

crates/core/src/linalg/fci.rs

@@ -29,6 +29,8 @@
 //!   (valid only if `p` is occupied) both carry the sign ``(-1)^m`` where `m` is the number of
 //!   occupied orbitals with index strictly greater than `p`.
 
+use std::ops::Range;
+
 use num_complex::Complex64;
 
 /// The largest number of spatial orbitals supported by the bitmask representation.
@@ -689,13 +691,71 @@ fn swap_if(swap: bool, a: usize, b: usize) -> (usize, usize) {
     if swap { (b, a) } else { (a, b) }
 }
 
+/// Sums `value[k]` over the entries `k` in `range` that sit on the diagonal (`src[k] == dst[k]`).
+///
+/// Every trace contribution in [`CompiledSector::trace`] has this shape -- only the value being
+/// summed and the index range differ -- because a scatter `src -> dst` lands on the matrix diagonal
+/// exactly when its two addresses coincide. Taking the addresses as slices lets this serve both the
+/// whole-array [`ScaledTransitions`] and a single CSR segment of [`PhasedTransitions`]; the `Into`
+/// bound on the summand covers the phase blocks, which store `i8` but accumulate in `f64`.
+fn diagonal_sum<T, V>(range: Range<usize>, src: &[usize], dst: &[usize], value: &[V]) -> T
+where
+    T: std::iter::Sum<T>,
+    V: Copy + Into<T>,
+{
+    range
+        .filter(|&k| src[k] == dst[k])
+        .map(|k| value[k].into())
+        .sum()
+}
+
 impl CompiledSector {
     /// The FCI dimension of the sector this map was compiled for.
     #[inline]
     pub fn dim(&self) -> usize {
         self.dim
     }
 
+    /// The exact trace of the compiled operator on this FCI sector.
+    pub fn trace(&self) -> Complex64 {
+        match &self.kind {
+            CompiledKind::Spinless { entries } => entries
+                .iter()
+                .filter_map(|&(src, dst, weight)| (src == dst).then_some(weight))
+                .sum(),
+            CompiledKind::Spinful(spinful) => {
+                // An identity spin block contributes its own dimension as a multiplicity: every
+                // determinant of the untouched block repeats the other block's diagonal entry.
+                let scaled_trace = |transitions: &ScaledTransitions| {
+                    diagonal_sum::<Complex64, _>(
+                        0..transitions.scale.len(),
+                        &transitions.src,
+                        &transitions.dst,
+                        &transitions.scale,
+                    )
+                };
+                // One mixed term's trace within a single spin block: its CSR segment only.
+                let phased_trace = |block: &PhasedTransitions, term: usize| {
+                    diagonal_sum::<f64, _>(
+                        block.indptr[term]..block.indptr[term + 1],
+                        &block.src,
+                        &block.dst,
+                        &block.phase,
+                    )
+                };
+                let mut trace = spinful.scalar * (spinful.dim_a * spinful.dim_b) as f64;
+                trace += scaled_trace(&spinful.alpha_only) * spinful.dim_b as f64;
+                trace += scaled_trace(&spinful.beta_only) * spinful.dim_a as f64;
+                for (term, &coeff) in spinful.mixed_coeffs.iter().enumerate() {
+                    trace += coeff
+                        * phased_trace(&spinful.mixed_alpha, term)
+                        * phased_trace(&spinful.mixed_beta, term);
+                }
+                trace
+            }
+        }
+    }
+
     /// Applies the compiled operator to a state vector: `out = op @ vec`.
     ///
     /// `vec` must have length [`Self::dim`]. Validated bit-for-bit in the tests against an
@@ -1571,6 +1631,70 @@ mod tests {
         }
     }
 
+    fn reference_trace(
+        norb: u32,
+        n_alpha: u32,
+        n_beta: Option<u32>,
+        terms: &[(Complex64, Vec<bool>, Vec<u32>)],
+    ) -> Complex64 {
+        let table = BinomialTable::new(norb);
+        let dim = table.num_strings(norb, n_alpha)
+            * n_beta.map_or(1, |n_beta| table.num_strings(norb, n_beta));
+        let mut trace = Complex64::new(0.0, 0.0);
+        for diagonal in 0..dim {
+            let mut basis = vec![Complex64::new(0.0, 0.0); dim];
+            basis[diagonal] = Complex64::new(1.0, 0.0);
+            trace += reference_matvec(norb, n_alpha, n_beta, terms, &basis)[diagonal];
+        }
+        trace
+    }
+
+    #[test]
+    fn compiled_trace_matches_reference() {
+        let spinless_terms = vec![
+            (Complex64::new(1.2, -0.1), vec![], vec![]),
+            (Complex64::new(0.5, 0.2), vec![true, false], vec![0, 0]),
+            (Complex64::new(0.7, 0.0), vec![true, false], vec![0, 1]),
+            (
+                Complex64::new(1.3, -0.4),
+                vec![true, false, true, false],
+                vec![0, 0, 2, 2],
+            ),
+        ];
+        let spinless = SpinlessSector::new(4, 2)
+            .compile(term_views(&spinless_terms))
+            .unwrap();
+        let expected = reference_trace(4, 2, None, &spinless_terms);
+        assert!((spinless.trace() - expected).norm() < 1e-12);
+
+        let norb = 3;
+        let spinful_terms = vec![
+            (Complex64::new(0.8, 0.1), vec![], vec![]),
+            (Complex64::new(0.4, 0.0), vec![true, false], vec![0, 0]),
+            (
+                Complex64::new(-0.2, 0.3),
+                vec![true, false],
+                vec![norb, norb],
+            ),
+            (
+                Complex64::new(1.1, -0.2),
+                vec![true, false, true, false],
+                vec![0, 0, norb, norb],
+            ),
+            (
+                Complex64::new(0.6, 0.0),
+                vec![true, true, false, false],
+                vec![0, norb + 2, norb, 2],
+            ),
+        ];
+        let spinful = SpinfulSector::new(norb, 2, 1)
+            .unwrap()
+            .compile(term_views(&spinful_terms))
+            .unwrap();
+        let expected = reference_trace(norb, 2, Some(1), &spinful_terms);
+        assert!((spinful.trace() - expected).norm() < 1e-12);
+    }
+
     #[test]
     fn spinless_matvec_matches_reference() {
         // Cover a range of orbital counts / fillings and several kinds of terms.

crates/pyext/src/linalg/fci.rs

@@ -42,6 +42,7 @@ type Matvec = Box<dyn Fn(&[Complex64]) -> Result<Vec<Complex64>, FciMatvecError>
 #[pyclass(module = "qiskit_fermions.linalg.fci", name = "FciLinearOperator")]
 pub struct FciLinearOperator {
     dim: usize,
+    trace: Complex64,
     matvec: Matvec,
     rmatvec: Matvec,
 }
@@ -52,9 +53,10 @@ impl FciLinearOperator {
     /// `dim` is the FCI sector dimension (the length of the state vectors this operator acts on) and
     /// is exposed to Python as the square :attr:`shape` `(dim, dim)`. `matvec` applies the operator;
     /// `rmatvec` applies its adjoint (`A.H @ v`), which SciPy's `expm_multiply` requires.
-    pub fn new(dim: usize, matvec: Matvec, rmatvec: Matvec) -> Self {
+    pub fn new(dim: usize, trace: Complex64, matvec: Matvec, rmatvec: Matvec) -> Self {
         Self {
             dim,
+            trace,
             matvec,
             rmatvec,
         }
@@ -93,6 +95,12 @@ impl FciLinearOperator {
         numpy.getattr("dtype")?.call1(("complex128",))
     }
 
+    /// The exact trace of the operator on this FCI sector.
+    #[getter]
+    fn trace(&self) -> Complex64 {
+        self.trace
+    }
+
     /// Applies the operator to a state vector: returns ``op @ vec``.
     ///
     /// Args:

crates/pyext/src/operators/fermion_operator.rs

@@ -1277,10 +1277,11 @@ impl PyFermionOperator {
                     .compile_fci_spinless(&sector)
                     .map_err(crate::value_err)?,
             );
+            let trace = compiled.trace();
             let (compiled_mv, compiled_rmv) = (Arc::clone(&compiled), compiled);
             let matvec = Box::new(move |vec: &[Complex64]| compiled_mv.apply(vec));
             let rmatvec = Box::new(move |vec: &[Complex64]| compiled_rmv.apply_conj(vec));
-            Ok(FciLinearOperator::new(dim, matvec, rmatvec))
+            Ok(FciLinearOperator::new(dim, trace, matvec, rmatvec))
         } else {
             let (n_alpha, n_beta) = nelec.extract::<(u32, u32)>()?;
             let sector = SpinfulSector::new(norb, n_alpha, n_beta).map_err(crate::value_err)?;
@@ -1290,10 +1291,11 @@ impl PyFermionOperator {
                     .compile_fci_spinful(&sector)
                     .map_err(crate::value_err)?,
             );
+            let trace = compiled.trace();
             let (compiled_mv, compiled_rmv) = (Arc::clone(&compiled), compiled);
             let matvec = Box::new(move |vec: &[Complex64]| compiled_mv.apply(vec));
             let rmatvec = Box::new(move |vec: &[Complex64]| compiled_rmv.apply_conj(vec));
-            Ok(FciLinearOperator::new(dim, matvec, rmatvec))
+            Ok(FciLinearOperator::new(dim, trace, matvec, rmatvec))
         }
     }
 }

python/qiskit_fermions/circuit/library/_expm_multiply.py

@@ -0,0 +1,66 @@
+# This code is a Qiskit project.
+#
+# (C) Copyright IBM 2026.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at https://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Shared matrix-exponential helpers for fermionic circuit simulation."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, cast
+
+import numpy as np
+import scipy.sparse.linalg
+
+from qiskit_fermions.protocols.linear_operator import scipy_linear_operator_from_kernel
+
+if TYPE_CHECKING:
+    from qiskit_fermions.protocols.linear_operator import _SupportsFciLinearOperator
+
+
+def _expm_multiply_fci(
+    operator: _SupportsFciLinearOperator,
+    vec: np.ndarray,
+    norb: int,
+    nelec: int | tuple[int, int],
+    scale: complex = 1.0,
+) -> np.ndarray:
+    r"""Applies ``exp(scale * operator)`` to ``vec`` on the ``(norb, nelec)`` FCI sector.
+
+    This goes through the internal ``_fci_linear_operator_`` carrier rather than the public
+    :func:`~qiskit_fermions.linalg.linear_operator` helper, because the native kernel it returns
+    exposes the operator's exact fixed-sector ``trace`` alongside its matrix-vector action. SciPy
+    uses that trace to precondition the exponential -- it factors out ``exp(traceA / n)`` -- which is
+    not a correctness input but is a large win in both speed and accuracy for an operator with a
+    sizeable trace. A SciPy ``LinearOperator`` has nowhere to carry the value (and drops attributes
+    when scaled), so the trace is read off the kernel here instead of being transported through the
+    public wrapper.
+
+    ``scale`` multiplies the operator, so it must multiply the trace too: ``trace(c * A)`` is
+    ``c * trace(A)``. It is applied to the trace *before* SciPy sees it, since ``scale * linop``
+    produces a composed operator that no longer carries the kernel's metadata.
+
+    Args:
+        operator: the operator to exponentiate, exposing the internal FCI kernel carrier.
+        vec: the state vector to apply the exponential to.
+        norb: the number of spatial orbitals.
+        nelec: the electron count -- an integer for a spinless sector, or an ``(n_alpha, n_beta)``
+            pair for a spinful one.
+        scale: a scalar multiplying the operator inside the exponential.
+
+    Returns:
+        The vector ``exp(scale * operator) @ vec``.
+    """
+    kernel = operator._fci_linear_operator_(norb, nelec)
+    linop = scipy_linear_operator_from_kernel(kernel)
+    return cast(
+        np.ndarray,
+        scipy.sparse.linalg.expm_multiply(scale * linop, vec, traceA=scale * kernel.trace),
+    )

python/qiskit_fermions/circuit/library/evolution.py

@@ -16,12 +16,10 @@
 
 from typing import TYPE_CHECKING
 
-import scipy.sparse.linalg
-
-from qiskit_fermions.linalg import linear_operator
 from qiskit_fermions.operators import OperatorTrait
 
 from .. import FermionicGate
+from ._expm_multiply import _expm_multiply_fci
 
 if TYPE_CHECKING:
     import numpy as np
@@ -113,8 +111,8 @@ def _apply_unitary_placed_(
         """Applies ``exp(-i * time * operator)`` after relabeling the operator to global modes.
 
         The operator is relabeled onto its global modes and turned into a ``scipy`` ``LinearOperator``
-        via :func:`.linear_operator` (backed by a native FCI matrix-vector kernel), then applied to the
-        vector via ``scipy.sparse.linalg.expm_multiply``. This mirrors ffsim's own ``_apply_unitary_``
+        backed by a native FCI matrix-vector kernel, then applied to the vector via
+        ``scipy.sparse.linalg.expm_multiply``. This mirrors ffsim's own ``_apply_unitary_``
         implementations (e.g. for its UCCSD operators).
 
         Args:
@@ -178,9 +176,4 @@ def _apply_unitary_placed_(
                 + (" of each spin species" if not isinstance(nelec, int) else "")
                 + f" (norb={norb}, nelec={nelec})."
             )
-        linop = linear_operator(operator, norb, nelec)
-        # ``traceA`` is only a balancing hint for scipy (it factors out ``exp(traceA / n)`` to
-        # improve conditioning), not a correctness input: an inexact value costs at most some
-        # numerical conditioning. Passing 0.0 avoids scipy estimating the trace itself and mirrors
-        # ffsim's own ``_apply_unitary_`` implementations.
-        return scipy.sparse.linalg.expm_multiply(-1j * self.params[0] * linop, vec, traceA=0.0)
+        return _expm_multiply_fci(operator, vec, norb, nelec, scale=-1j * self.params[0])

python/qiskit_fermions/circuit/library/orbital_rotation.py

@@ -15,15 +15,14 @@
 from __future__ import annotations
 
 import numbers
-from typing import cast
 
 import numpy as np
 import scipy.linalg
-import scipy.sparse.linalg
 
 from qiskit_fermions.utils.optionals import HAS_FFSIM
 
 from .. import FermionicGate
+from ._expm_multiply import _expm_multiply_fci
 
 
 class OrbitalRotation(FermionicGate):
@@ -120,10 +119,9 @@ def _apply_unitary_placed_(
         - **General path**: otherwise (i.e. when ``ffsim`` is unavailable) the rotation is applied as
           the evolution :math:`\exp(G)` under its generator
           :math:`G = \sum_{ij} \log(U)_{ij} a^\dagger_i a_j`, where :math:`U` is the embedded matrix.
-          :math:`G` is turned into a ``scipy`` ``LinearOperator`` via
-          :func:`.linear_operator` (backed by the native FCI matrix-vector kernel) and applied via
-          :func:`scipy.sparse.linalg.expm_multiply`. This mirrors
-          :meth:`.Evolution._apply_unitary_placed_`.
+          :math:`G` is turned into a ``scipy`` ``LinearOperator`` backed by the native FCI
+          matrix-vector kernel and applied via :func:`scipy.sparse.linalg.expm_multiply`. This
+          mirrors :meth:`.Evolution._apply_unitary_placed_`.
 
         Args:
             vec: the state vector to act on.
@@ -203,7 +201,6 @@ def _apply_via_generator(
         ``logm`` preserves that block-diagonal structure, so the generator ``G`` conserves the
         ``(norb, nelec)`` sector and no amplitude is dropped.
         """
-        from qiskit_fermions.linalg import linear_operator
         from qiskit_fermions.operators import FermionOperator
 
         if copy:
@@ -224,10 +221,4 @@ def _apply_via_generator(
         }
         generator = FermionOperator.from_dict(terms)  # type: ignore[arg-type]
 
-        # ``_linear_operator_`` is monkeypatched onto ``FermionOperator`` at import time (in
-        # ``qiskit_fermions.operators``); the stubs type ``from_dict``'s result as the compiled
-        # ``_lib`` type, which does not carry the patched method, so mypy still needs the ignore.
-        linop = linear_operator(generator, norb, nelec)  # type: ignore[arg-type]
-        # ``traceA=0.0`` mirrors Evolution._apply_unitary_placed_: it is only a scipy conditioning
-        # hint (it factors out ``exp(traceA / n)``), not a correctness input.
-        return cast(np.ndarray, scipy.sparse.linalg.expm_multiply(linop, vec, traceA=0.0))
+        return _expm_multiply_fci(generator, vec, norb, nelec)