Mitigation (Zero-Noise Extrapolation & Density-Matrix Diagnostics)¶
Correcting a quantum measurement result, not a numeric log/trajectory — see Concepts if you're looking for Vector Healing instead.
A real circuit run on noisy hardware gives the wrong answer. Zero-Noise Extrapolation (ZNE) gets closer to the right one without needing a better device: run the same circuit at several deliberately-worsened noise strengths, then extrapolate the trend back to what zero noise would have given.
Step 1. Run the circuit you want to correct¶
import numpy as np
import dense_evolution as de
qasm = 'OPENQASM 2.0; include "qelib1.inc"; qreg q[2]; h q[0]; cx q[0],q[1];'
circuit = de.QASMParser().parse(qasm)
sim = de.DenseSVSimulator(2)
sim.run_circuit_jit(circuit.to_tuples())
sv0 = np.asarray(sim.get_statevector())
round(float(abs(sv0[0]) ** 2), 4)
sv0 is the exact, noiseless Bell state this whole page corrects a noisy version of.
abs(sv0[0]) ** 2 is the population of |00> — exactly 0.5, since a Bell state is an
equal mix of |00> and |11>. Every step below tries to recover this same number from
noisy data.
Step 2. The same circuit, run through a noisy channel at increasing strength¶
from dense_evolution.noise import NoiseModel
scales = (1.0, 2.0, 3.0)
rng = np.random.default_rng(0)
pop00 = [
np.mean([abs(NoiseModel.apply_to_sv(sv0.copy(), 2, "depolarizing", 0.05 * s, rng=rng)[0]) ** 2
for _ in range(20000)])
for s in scales
]
[round(float(x), 4) for x in pop00]
Each entry is |00>'s population averaged over 20000 noisy runs at one noise scale —
1x, 2x, and 3x a base error rate (NoiseModel.apply_to_sv, the same
function Noise introduces). As the scale grows, the noisy population drifts
further from Step 1's ideal 0.5. These three numbers, paired with scales, are exactly
what every extrapolation function below expects as input.
Step 3. Extrapolate back to zero noise¶
from dense_evolution.mitigation import richardson_extrapolate
round(float(richardson_extrapolate(pop00, scales)), 4)
richardson_extrapolate fits an exact curve through the 3 noisy points from Step 2 and
reads off its value at noise scale 0 — recovering Step 1's ideal 0.5 to within rounding
error, from data that never included it. This is the core of every other function on
this page; zero_noise_extrapolation (below) and zne_density_matrix (Step 5) both
call it, or its density-matrix generalization, internally.
Step 4. More noise scales than 3 — polynomial_extrapolate¶
from dense_evolution.mitigation import polynomial_extrapolate
more_scales = (1.0, 2.0, 3.0, 4.0, 5.0)
more_pop00 = [
np.mean([abs(NoiseModel.apply_to_sv(sv0.copy(), 2, "depolarizing", 0.05 * s, rng=rng)[0]) ** 2
for _ in range(20000)])
for s in more_scales
]
round(float(polynomial_extrapolate(more_pop00, more_scales, degree=2)), 4)
richardson_extrapolate needs exactly as many points as it has degrees of freedom, so
its fit becomes numerically unstable with many closely-spaced scales. polynomial_extrapolate
fits a lower-degree polynomial (degree=2 here) by least squares instead — at exactly 3
points it's mathematically identical to Step 3, but with extra points it averages down
noise instead of forcing an increasingly ill-conditioned exact fit through every one of
them.
Step 5. A real experiment gives you a density matrix, not one number¶
import jax.numpy as jnp
from dense_evolution.mitigation import zne_density_matrix
def noisy_rho(p, k, rng):
rho = np.zeros((4, 4), dtype=np.complex128)
for _ in range(k):
s = NoiseModel.apply_to_sv(sv0.copy(), 2, "depolarizing", p, rng=rng)
rho += np.outer(s, s.conj())
return rho / k
rho_at_scales = jnp.stack([noisy_rho(0.05 * s, 200, rng) for s in scales])
corrected = zne_density_matrix(rho_at_scales, scales)
round(float(jnp.trace(corrected).real), 6)
rho_at_scales[i] is the full noisy density matrix at scale scales[i] — the
density-matrix counterpart of Step 2's pop00. Extrapolating a whole matrix the way
Step 4 extrapolates one number does not generally give back a valid density matrix
(negative eigenvalues can appear even though every input matrix was physical);
zne_density_matrix runs polynomial_extrapolate on the whole matrix and then projects
the result onto the nearest true density matrix (project_to_physical, Details below)
so corrected is always physical — trace 1, as shown above — even when the raw
extrapolation wasn't.
Step 6. Grade the correction¶
from dense_evolution.mitigation import uhlmann_fidelity
rho_ideal = jnp.asarray(np.outer(sv0, sv0.conj()), dtype=jnp.complex128)
raw_fidelity = uhlmann_fidelity(rho_at_scales[0], rho_ideal)
corrected_fidelity = uhlmann_fidelity(corrected, rho_ideal)
round(float(raw_fidelity), 4), round(float(corrected_fidelity), 4)
uhlmann_fidelity compares a density matrix against a known ideal state — here, the raw
base-scale noisy result from Step 5 against rho_ideal (built from Step 1's sv0,
never fed into Steps 2-5), versus the same comparison after correction. rho_ideal is
only ever used for this final grading step — feeding a known-ideal state into the
extrapolation or projection steps themselves would be oracle access, not error
mitigation.
Step 7. When the decay is exponential, not polynomial — bounded_exponential_extrapolate¶
from dense_evolution.observables import pauli_expectation
from dense_evolution.mitigation import bounded_exponential_extrapolate
rng = np.random.default_rng(1)
zz_at_scales = [
np.mean([pauli_expectation(np.asarray(
NoiseModel.apply_to_sv(sv0.copy(), 2, "depolarizing", 0.25 * s, rng=rng)), "ZZ")
for _ in range(200)])
for s in scales
]
[round(float(x), 4) for x in zz_at_scales]
round(float(bounded_exponential_extrapolate(zz_at_scales, list(scales))), 4)
At this much higher noise (0.25 * s instead of Step 2's 0.05 * s), sv0's ideal
ZZ expectation of 1.0 decays fast enough that an ordinary unconstrained
a + b*exp(-c*lambda) fit on these same 3 points doesn't just misfire — it fails to
converge at all (scipy.optimize.curve_fit raises RuntimeError: Optimal parameters
not found). bounded_exponential_extrapolate reparametrizes the same exponential
model so the zero-noise value is an explicit, constrained parameter
(Miranskyy, Sorrenti, Thind & Gravel, arXiv:2604.24475) and
recovers the ideal 1.0 exactly. Reach for this instead of Step 3/4's polynomial
fits when the underlying decay is closer to exponential than polynomial — the usual
shape for a single depolarizing-type channel.
Details¶
Healing-adapted zero-noise extrapolation¶
zero_noise_extrapolation(expectation_values, noise_factors) is richardson_extrapolate
by default (Step 3 above), but accepts an optional sigma_at_base_noise: when given
(alongside exactly 3 noise factors — it raises NotImplementedError for any other
count), the 3 Richardson coefficients are perturbed by
dense_evolution.healing's calculate_delta_preemp before renormalizing,
nudging the extrapolation when the measured coherence signal is off an ideal target
instead of trusting the 3 raw points equally.
JSD-informed density-matrix correction¶
jsd_predictive_zne_density_matrix(rho_at_scales, noise_factors) is zne_density_matrix
for exactly 3 equally-spaced scales, with a further nudge based on how much the
noise-scale-to-output-distribution relationship deviates from smooth (measured via
Jensen-Shannon divergence between consecutive scales' measurement distributions,
positive-only — the nudge is rectified to 0 whenever the signal would predict a
different direction, since an earlier unrectified version helped only 5/16 test points
despite the signal itself correlating with success). Validated on 46 real activated
points (photon-loss noise, 6 independent seeds): improves fidelity on 35/46 (76.1%),
mean gain +0.0055. Needs no oracle access to an ideal state, unlike naively reusing
calculate_delta_preemp with an external signal (tried first, found negligible).
project_to_physical¶
Projects a Hermitian, trace-1 matrix onto the nearest true density matrix (Hermitian,
trace 1, positive-semidefinite) in Frobenius distance — the same problem
Smolin, Gambetta & Smith 2012 solve by iterative
eigenvalue clipping, solved here instead as Euclidean projection of the eigenvalues onto
the probability simplex: a fully vectorized, jax.jit-traceable algorithm for the same
convex optimization problem (unique global minimum, so any correct algorithm agrees).
Verified to machine precision (~1e-15) against the paper's own worked example.
uhlmann_fidelity stays finite at degenerate eigenvalues¶
uhlmann_fidelity is differentiable through both arguments, including when rho_A has
(near-)degenerate eigenvalues — common for a near-pure state's noisy density matrix.
JAX's built-in eigh gradient divides by lambda_i - lambda_j and returns NaN there
(a known upstream limitation); uhlmann_fidelity uses a custom JVP rule internally that
masks that term to 0 for near-degenerate pairs instead, verified against finite
differences to exact agreement across non-degenerate, 2-fold, and 3-fold degenerate test
matrices.
jax.jit-compiled entry points¶
Every function above except Step 7's bounded_exponential_extrapolate has a _jit
counterpart (richardson_extrapolate_jit, zero_noise_extrapolation_jit,
polynomial_extrapolate_jit, uhlmann_fidelity_jit, zne_density_matrix_jit) for
callers inside an already-jitted pipeline (e.g.
jax.lax.scan) who don't want a host round-trip per call. Each skips the eager
version's own dtype auto-detection and argument validation — callers pass already-cast
complex128/float64 arrays themselves — and polynomial_extrapolate_jit/
zne_density_matrix_jit require degree as a static argument.
Why zne_density_matrix defaults to degree=2, not exact interpolation¶
More noise-scale points make exact interpolation (richardson_extrapolate) worse
under real statistical noise — Lagrange coefficients grow with point count, so the fit
increasingly forces itself through every noisy sample exactly. Measured directly
(4 qubits, all 5 noise channels, 5 seeds, same total measurement budget): exact
interpolation's mean fidelity gain dropped from +0.148 (3 points) to +0.081 (5 points)
to -0.220 (5 closely-spaced points — actively worse than no correction). A degree-2
least-squares fit on the same extra points instead reduces variance (std 0.062 to
0.035-0.046) at comparable or better mean gain, which is why zne_density_matrix uses
polynomial_extrapolate rather than exact interpolation by default.
zne ¶
Zero-Noise Extrapolation (ZNE)
Standard error-mitigation entry points, named the way the field already
names them (Richardson extrapolation, noise factors, zero-noise
extrapolation -- same vocabulary as e.g. Mitiq's zne API), so callers
and tooling can find "ZNE" without first learning Dense-Evolution's
internal healing vocabulary.
This module composes dense_evolution.healing's existing primitives
(calculate_delta_preemp, ...) -- it does not rename or replace them.
richardson_extrapolate_jit
module-attribute
¶
jax.jit-compiled entry point for richardson_extrapolate. Unlike
polynomial_extrapolate/zne_density_matrix's jitted variants, no
argument needs to be marked static here -- n (the point count) is read
from lambdas.shape[0], itself always static under tracing, not from a
Python degree parameter.
values must already be complex128 or float64 (pick the dtype yourself
before calling -- this skips richardson_extrapolate's np.iscomplexobj
auto-detection, which isn't traceable) and lambdas a float64 array.
Verified to match richardson_extrapolate exactly on real and complex
input; JAX recompiles per distinct input shape/dtype, as usual.
zero_noise_extrapolation_jit
module-attribute
¶
jax.jit-compiled entry point for zero_noise_extrapolation's
healing-adapted branch (the sigma_at_base_noise is not None case --
the plain-Richardson case already has richardson_extrapolate_jit, use
that directly instead). values must already be complex128 or float64
with exactly 3 rows, sigma_at_base_noise a float64 scalar; the
lambdas.shape[0] != 3 validation and dtype auto-detection that
zero_noise_extrapolation does are both skipped here (not traceable) --
callers are responsible for passing exactly-3-row input themselves.
target_sigma_ideal is a plain Python float, fine to leave non-static
since it only ever multiplies/subtracts, no Python branching on its value.
polynomial_extrapolate_jit
module-attribute
¶
polynomial_extrapolate_jit = functools.partial(
jax.jit, static_argnames=("degree",)
)(_polynomial_extrapolate_core)
jax.jit-compiled entry point for polynomial_extrapolate, added for
consistency with every other function in this module (richardson_extrapolate_jit,
zero_noise_extrapolation_jit, uhlmann_fidelity_jit, zne_density_matrix_jit)
-- until now this was the one function whose _core existed (used
internally by zne_density_matrix_jit) but had no standalone public jit
entry point of its own.
values must already be cast to its final dtype (complex128 or float64)
and lambdas a float64 array -- this skips polynomial_extrapolate's
np.iscomplexobj auto-detection, not traceable. degree is static (same
constraint as zne_density_matrix_jit). Verified to match
polynomial_extrapolate exactly.
uhlmann_fidelity_jit
module-attribute
¶
jax.jit-compiled entry point for uhlmann_fidelity. Both rho_A/
rho_B must already be complex128 (this skips uhlmann_fidelity's
own jnp.asarray(..., dtype=jnp.complex128) cast, itself trace-safe, but
kept out of the core to mirror the other _core functions' convention).
Returns a jnp scalar, not a Python float -- call float(...) yourself
if you need one outside a jitted context. Verified to match uhlmann_fidelity
exactly (same underlying math, just not cast to a Python float).
zne_density_matrix_jit
module-attribute
¶
zne_density_matrix_jit = functools.partial(
jax.jit, static_argnames=("degree",)
)(_zne_density_matrix_core)
jax.jit-compiled entry point for zne_density_matrix, for callers
inside a jitted pipeline (e.g. jax.lax.scan in MPSSimulator.run_circuit_jit)
who don't want a host round-trip every call -- zne_density_matrix itself
stays eager (unchanged) for one-off/interactive use, where jit compilation
overhead isn't worth paying for a single call.
degree is a static argument (must be a Python int, not a traced value --
pass it positionally or by keyword the same way every call, since JAX
recompiles per distinct static value). rho_at_scales must already be
complex128 and noise_factors a plain float array/sequence -- unlike
zne_density_matrix, this skips the np.iscomplexobj dtype auto-detection
(not traceable) and always assumes complex input, which is the only case
that makes sense for density matrices.
Measured speedup is real but size- and call-pattern-dependent (project_to_physical
alone measured 2x-22x across 2x2 to 32x32 matrices when jitted vs. the
previous non-jittable version) -- benchmark your own use case rather than
assuming a fixed number; the benefit only appears once compiled and called
repeatedly, a single one-off call pays the compilation cost first.
richardson_extrapolate ¶
Polynomial (Lagrange) Richardson extrapolation to zero noise.
expectation_values[i] is the value measured/simulated at noise scale
noise_factors[i] (e.g. 1x, 2x, 3x folded/scaled noise) -- a scalar,
or itself an array (e.g. a full probability distribution sampled at
that noise scale; extrapolated elementwise). Returns the extrapolated
zero-noise estimate, same shape as one expectation_values[i]. Works
for any number of points and any (not necessarily equally spaced)
noise factors; for the common 3-point case at noise_factors=(1,2,3)
this reduces exactly to the textbook coefficients (3, -3, 1).
expectation_values may be complex (e.g. density matrix entries,
which are complex off-diagonal in general) -- dtype is picked from
the input itself (complex128 if complex, float64 otherwise, matching
this function's previous always-float64 behavior for real input
exactly). Found via a real test case: forcing float64 unconditionally
silently discarded the imaginary part of complex input with no
visible error, only a low-signal ComplexWarning easy to miss --
confirmed directly (richardson_extrapolate([1+2j, 3+4j], ...) used
to return a purely real result, dropping real information).
Passing more noise-scale points makes exact interpolation MORE, not
less, sensitive to shot noise: at noise_factors equally spaced in
[1, 3], the noise-amplification factor kappa (richardson_amplification_factor,
the sum of absolute Lagrange coefficients) is 7 at 3 points, 129 at 5,
2815 at 7 -- verified directly, exact integers, not estimates. Since
the coefficients are fixed constants (independent of the measured
values), i.i.d. shot noise of per-point standard deviation sigma
propagates to an extrapolated-result standard deviation of
sigma * sqrt(sum(coeff_i**2)): at sigma=0.01 this is 0.044 (3
points), 0.67 (5), 13.2 (7) -- verified with this exact formula, not
simulated, and matching a real prior audit's Monte Carlo figures
(0.044 / 0.67 / 13.5) to within sampling noise. Against this,
prog.txt's own audit reports the systematic (noise-free) bias
improving by only about 0.01 over the same 3-to-7-point range on its
real experimental setup -- a bad trade past a handful of points. A
UserWarning fires when kappa exceeds _KAPPA_WARNING_THRESHOLD
(50); see richardson_amplification_factor to check kappa before
extrapolating.
Examples:
>>> from dense_evolution.mitigation import richardson_extrapolate
>>> round(float(richardson_extrapolate([0.90, 0.80, 0.65], [1.0, 2.0, 3.0])), 4)
0.95
Source code in dense_evolution/mitigation/zne.py
richardson_amplification_factor ¶
Noise-amplification factor kappa = sum(|Lagrange coefficients|) for
Richardson extrapolation at the given noise_factors: how much a
unit of i.i.d. shot noise spread evenly across the measured points
gets amplified in the extrapolated zero-noise estimate (each
coefficient can be much larger than 1 and they alternate in sign, so
they do not cancel in the worst case the way their SUM, which is
always exactly 1, might suggest).
Measured at noise_factors equally spaced in [1, 3]: kappa = 7 at 3
points, 129 at 5 points, 2815 at 7 points -- see
richardson_extrapolate's own docstring for the resulting bias/
variance trade-off in these same three cases.
Examples:
>>> from dense_evolution.mitigation import richardson_amplification_factor
>>> round(richardson_amplification_factor([1.0, 2.0, 3.0]), 4)
7.0
Source code in dense_evolution/mitigation/zne.py
zero_noise_extrapolation ¶
zero_noise_extrapolation(
expectation_values,
noise_factors,
sigma_at_base_noise=None,
target_sigma_ideal: float = 10.0,
) -> jnp.ndarray
Zero-Noise Extrapolation -- plain, or healing-adapted when a coherence signal is available.
Without sigma_at_base_noise: standard Richardson ZNE
(richardson_extrapolate).
With sigma_at_base_noise (the measured/simulated coherence sigma at
the base, unscaled noise level): the 3 Richardson coefficients are
perturbed by dense_evolution.healing.calculate_delta_preemp -- the
normalized deviation between the observed sigma and the ideal target
-- then renormalized to sum to 1. This is Dense-Evolution's
"predictive healing" ZNE variant: when the observed coherence is off
the ideal target, the extrapolation is nudged accordingly instead of
trusting the 3 raw noise-scaled points equally.
The healing-adapted path currently only supports exactly 3 noise
factors (the case it has been derived and tested against); passing
sigma_at_base_noise with any other point count raises
NotImplementedError rather than silently generalizing an unverified
formula.
Source code in dense_evolution/mitigation/zne.py
polynomial_extrapolate ¶
Least-squares polynomial extrapolation to zero noise.
Generalizes richardson_extrapolate: fits a degree-degree polynomial
to (noise_factors, expectation_values) by ordinary least squares and
evaluates it at zero. With exactly degree + 1 points the fit is the
unique interpolating polynomial, mathematically identical to
richardson_extrapolate at that point count (verified directly, both
for real and complex input). With MORE than degree + 1 points it
becomes an overdetermined fit -- the extra points average down
statistical noise instead of forcing the polynomial through every
noisy sample exactly, trading a small amount of interpolation bias
for reduced variance.
This matters in practice, not just in theory: adding more noise-scale
points to exact interpolation (richardson_extrapolate) makes
extrapolation WORSE under real statistical noise, because Lagrange
coefficients grow with point count (worse still with closely-spaced
points -- a Runge's-phenomenon-like effect). Measured directly on the
density-matrix healing experiment (experiments/matrix_healing_zne_sweep.py's
setup, n=4 qubits, all 5 noise channels, 5 seeds): exact interpolation's
mean fidelity-delta dropped from +0.148 (3 points) to +0.081 (5 points,
same spacing) to -0.220 (5 points, denser spacing -- actively worse
than not correcting). A degree-2 least-squares fit fed the same extra
points instead REDUCES variance (std 0.062 -> 0.035-0.046) at
comparable or better mean delta, because the extra points are no
longer forced to satisfy an increasingly ill-conditioned exact fit.
This is why zne_density_matrix uses this function (degree=2) instead
of richardson_extrapolate by default.
Raises ValueError if fewer than degree + 1 points are given (the fit
would be underdetermined).
Source code in dense_evolution/mitigation/zne.py
bounded_exponential_extrapolate ¶
bounded_exponential_extrapolate(
expectation_values, noise_factors, bound: float = 1.0
) -> jnp.ndarray
Physically bounded exponential Zero-Noise Extrapolation (Miranskyy, Sorrenti, Thind & Gravel, arXiv:2604.24475, "Improving Zero-Noise Extrapolation via Physically Bounded Models").
Unlike polynomial_extrapolate, this fits an EXPONENTIAL model --
appropriate when the expectation value decays roughly exponentially with
noise strength (the usual case for a depolarizing-type channel), not a
polynomial one -- and is Dense-Evolution's first exponential-family ZNE
model. The zero-noise value is made an explicit model parameter via the
reparametrization
E(lambda) = a + (zeta - a) * exp(-c * lambda), E(0) = zeta
and zeta is constrained to [-bound, bound] during the fit (the
physically valid range for a +-1-eigenvalue Pauli observable, bound=1.0
by default). An ordinary unconstrained a + b*exp(-c*lambda) fit has no
such guarantee and can produce a wildly out-of-range or non-finite
"zero-noise" estimate when the data is noisy or only a few points are
available -- verified directly on this project's own depolarizing-noise
ZZ-expectation setup (real Bell circuit, 30 random noise seeds, 3 noise
scales): the unconstrained fit landed outside [-1, 1] or failed to
converge in 21/30 seeds (mean absolute error over 200000, dominated by
those blow-ups), versus 0/30 for this bounded fit (mean absolute error
0.066) -- see Dense-Evolution-Discovery/scripts/zne_physically_bounded.py
for the full reproduction.
Fit via SciPy's constrained L-BFGS-B, not jax.jit-traceable like the
rest of this module -- the optimization itself runs in plain NumPy, so
(unlike every other function here) there is no _jit variant. Multi-start
(_N_MULTI_STARTS starts: the original single fixed start
[0, values[0], 0.5] first, then _N_MULTI_STARTS - 1 deterministic
random ones, fixed seed _MULTI_START_SEED, keeping the converged fit
with lowest loss) -- keeping the original start first means a case
where it was already the global optimum is unaffected bit-for-bit.
The non-convex 3-parameter fit from a single fixed starting point can
otherwise land in a poor local minimum -- on data generated exactly
from this function's own model, the single-start fit was off by up to
3.9e-3, versus under 1e-5 with multi-start. Raises RuntimeError if every
start fails to converge (result.success), rather than silently
returning an unconverged result.x[1].
Examples:
>>> from dense_evolution.mitigation import bounded_exponential_extrapolate
>>> round(float(bounded_exponential_extrapolate([0.49, 0.1, 0.1], [1.0, 2.0, 3.0])), 4)
1.0
Source code in dense_evolution/mitigation/zne.py
322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 | |
project_to_physical ¶
Project a Hermitian, trace-1 candidate matrix onto the nearest physical density matrix (Hermitian, trace 1, positive-semidefinite) in 2-norm/Frobenius distance -- the same problem Smolin, Gambetta & Smith, "Maximum Likelihood, Minimum Effort" (2012), arXiv:1106.5458, Fig. 1, solve with a sequential eigenvalue-clipping algorithm (sort eigenvalues descending, repeatedly zero the smallest remaining one and redistribute its negative mass over the rest, until the least would be non-negative).
Implemented here as Euclidean projection onto the probability simplex
(Held, Wolfe & Crowder 1974; also e.g. Duchi et al. 2008) applied to the
eigenvalues -- a different, fully vectorized algorithm for the exact
same convex optimization problem (unique global minimum, so any correct
algorithm must agree). No Python-level while/for loop over
eigenvalues (the original transcription's while i >= 0: ... isn't
jax.jit-traceable, forcing a host round-trip every call) -- this
version is pure jnp array ops plus one dynamic index (mus[k-1],
itself trace-safe), so it JIT-compiles cleanly.
Verified against the original transcription (which itself matches the
SGS paper's own worked example, eigenvalues 3/5, 1/2, 7/20, 1/10,
-11/20 -> 9/20, 7/20, 1/5, 0, 0): identical to machine precision on the
paper's example and on 30 random Hermitian trace-1 matrices (2-7 dim)
perturbed to be unphysical (max difference ~1e-15); confirmed to
actually compile and run under jax.jit.
richardson_extrapolate/polynomial_extrapolate's output on a stack
of density matrices is not itself generally a valid density matrix --
extrapolation can (and in practice does) produce small negative
eigenvalues even when every input matrix was physical. This is the
correction step, meant to run after extrapolation, not a
general-purpose "make anything a density matrix" tool (it assumes the
input is already Hermitian and trace 1 up to this function's own
re-Hermitization step below).
Source code in dense_evolution/mitigation/zne.py
uhlmann_fidelity ¶
Uhlmann fidelity F(rho_A, rho_B) = (Tr sqrt(sqrt(rho_A) rho_B sqrt(rho_A)))^2.
Reduces to |zne_density_matrix or any extrapolation step would be using held-out
ground truth to guide the algorithm (oracle access), not a legitimate
error-mitigation technique. Keeping ideal-state comparison to this
function only, rather than plumbing it into the correction functions
at all, makes that boundary structural rather than a convention callers
have to remember.
Computes Tr(sqrt(inner)) as sum(sqrt(eigenvalues of inner)) instead of
reconstructing the full matrix square root (sqrt(M) has the same
eigenvectors as M and sqrt-of-eigenvalues eigenvalues, so its trace is
exactly that sum) -- skips one eigenvector reconstruction, and avoids
matsqrt's float() cast that isn't jax.jit-traceable. Verified
against the previous full-reconstruction version: identical to machine
precision (~1e-16) on 30 random density-matrix pairs.
Differentiable through both arguments, including when rho_A has
(near-)degenerate eigenvalues (e.g. a near-pure state's noisy density
matrix, which typically has several near-zero, near-degenerate
eigenvalues) -- uses _eigh_degenerate_safe internally rather than
jnp.linalg.eigh directly, specifically to keep jax.grad(uhlmann_fidelity, ...)
finite in that case (see _eigh_degenerate_safe's own docstring).
Forward-pass value is bit-identical to jnp.linalg.eigh-based
computation (same underlying eigh call; only the backward rule
differs), verified end-to-end against the previous implementation.
Examples:
>>> import numpy as np
>>> from dense_evolution.mitigation import uhlmann_fidelity
>>> rho = np.array([[1, 0], [0, 0]], dtype=complex)
>>> round(float(uhlmann_fidelity(rho, rho)), 4)
1.0
>>> sigma = np.array([[0.5, 0], [0, 0.5]], dtype=complex)
>>> round(float(uhlmann_fidelity(rho, sigma)), 4)
0.5
Source code in dense_evolution/mitigation/zne.py
global_depolarizing_channel ¶
Global n-qubit depolarizing channel, D_p(rho) = (1-p)rho + (p/dim)I.
Distinct from NoiseModel's 'depolarizing' model, which applies an
independent PER-QUBIT local Kraus channel -- a different physical map
from this GLOBAL channel, which mixes the whole dim-dimensional state
toward the fully mixed state as one unit. Use this one when modeling
e.g. state-prep/measurement (SPAM) error reported as a single joint
depolarizing parameter over the whole register, not per-qubit gate
noise (promoted from a real reproduction of arXiv:2608.16716's own
SPAM model, Dense-Evolution-Discovery Experiment 33).
Source code in dense_evolution/noise/density_matrix_channels.py
amplitude_damping_channel ¶
Single-qubit amplitude-damping channel: E0 @ rho @ E0.conj().T + E1 @ rho @ E1.conj().T, with E0=diag(1, sqrt(1-gamma)) and E1=[[0,sqrt(gamma)],[0,0]] -- population only ever moves |1>->|0>, never the reverse.
Distinct from global_depolarizing_channel (symmetric, mixes toward
the fully-mixed state regardless of which state is |1> or |0>) -- this
one is asymmetric by construction, the real signature of energy-relaxation
(T1) processes and of quasiparticle poisoning (promoted from a real
reproduction of arXiv:2104.05219's measured cosmic-ray-induced error
bursts, Dense-Evolution-Discovery Experiment 34, where this asymmetry is
exactly the mechanism's own reported signature: decay errors only, no
excess excitation errors).
Single-qubit only (rho must be 2x2) -- unlike global_depolarizing_channel,
this is not dimension-generic, since amplitude damping is inherently a
per-qubit process, not a joint-register one.
Source code in dense_evolution/noise/density_matrix_channels.py
cosmic_ray_burst_profile ¶
cosmic_ray_burst_profile(
time_us,
baseline_gamma: float,
ratio_intermediate: float = 2.5,
ratio_peak: float = 3.75,
tau1_us: float = 3.0,
tau2_us: float = 300.0,
tau_decay_ms: float = 25.0,
) -> jnp.ndarray
Time-dependent decay-probability profile for a cosmic-ray/gamma-ray-
induced quasiparticle burst: a two-stage rise (fast to
ratio_intermediatex baseline, slower to ratio_peakx baseline) times
a single-exponential recovery, generalized out of a fixed, paper-number
validation (Dense-Evolution-Discovery Experiment 34, reproducing
arXiv:2104.05219's real measured event on a 26-qubit chip).
Feed the result to continuous_dissipative_evolve alongside
amplitude_damping_channel (or any other single-time-varying-parameter
channel) to inject a realistic burst into any circuit or QEC study,
without re-deriving this shape by hand each time.
The default ratios/timescales are the paper's own real numbers -- see
Experiment 34's docstring for exactly which are paper-fitted (the 25ms
decay) versus chosen to match the paper's two described rise points
(tau1/tau2). All are overridable for a different event severity or
device generation; baseline_gamma is never derived here -- pass
whatever per-slice decay probability corresponds to your own dt/T1
convention (see Experiment 34 for one worked example of that
conversion).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
time_us
|
array_like
|
Time since impact, in microseconds (t=0 is the impact instant). |
required |
baseline_gamma
|
float
|
Undisturbed per-slice decay probability; this profile scales it up, it does not derive it. |
required |
ratio_intermediate
|
float
|
Multiplier on |
2.5
|
ratio_peak
|
float
|
Multiplier on |
2.5
|
tau1_us
|
float
|
Rise timescales for the two saturating-exponential stages (paper defaults 3, 300 -- chosen to match its ~10us/~1ms descriptions, not fitted by the paper itself). |
3.0
|
tau2_us
|
float
|
Rise timescales for the two saturating-exponential stages (paper defaults 3, 300 -- chosen to match its ~10us/~1ms descriptions, not fitted by the paper itself). |
3.0
|
tau_decay_ms
|
float
|
Recovery time constant (paper default 25 -- its own fitted central value, real range 25-30ms across 415 events). |
25.0
|
Returns:
| Type | Description |
|---|---|
ndarray
|
Per-slice decay probability at each entry of |
Source code in dense_evolution/noise/cosmic_ray.py
zne_density_matrix ¶
Zero-Noise Extrapolation for density matrices.
rho_at_scales[i] is a noisy density-matrix estimate (e.g. from a
Monte-Carlo/shot ensemble) at noise scale noise_factors[i].
Extrapolates to zero noise via polynomial_extrapolate (least-squares,
complex-safe, degree=2 by default) and projects the result onto the
nearest physical density matrix (project_to_physical), since the raw
extrapolated matrix is not generally positive-semidefinite even when
every input was.
Uses polynomial_extrapolate rather than exact richardson_extrapolate
because, with exactly 3 noise scales (the original design point), the
two are mathematically identical -- but polynomial_extrapolate stays
well-behaved (reduced variance) when a caller passes MORE than 3 scales,
where exact interpolation instead gets WORSE (see
polynomial_extrapolate's docstring for the measured numbers). This
makes "pass more noise-scale points" a safe thing to try rather than a
trap.
Honest findings, both against a GHZ-state ideal target and
dense_evolution.registry.NoiseModel noise at base_p=0.05, scales
1x/2x/3x unless noted, uhlmann_fidelity against the true ideal state
used only to grade the result -- never as input to any step above:
experiments/matrix_healing_zne.py: 2-qubit Bell state, depolarizing noise, K=200-trajectory estimate per scale, averaged over 4 seeds -- raw fidelity ~0.865, corrected ~0.947 (+0.08), positive on every seed tested.experiments/matrix_healing_zne_sweep.py: 2-5 qubits x all 5NoiseModelchannels (depolarizing, bitflip, phaseflip, amplitude_damping, combined) x 5 seeds, K=400 trajectories per scale (100 runs total) -- 96/100 positive, mean delta +0.12, and every single (qubit count, noise channel) combination is net positive on average, growing to +0.20-0.25 at 5 qubits for depolarizing/bitflip. The 4 remaining negative runs are small (worst -0.02) and consistent with residual Monte Carlo noise, not a systematic failure mode. (These specific numbers were measured with exact 3-point interpolation, which is identical to this function's degree=2 default at 3 points -- unaffected by the switch.)- An earlier draft of this sweep (K=150, 3 seeds) had reported phaseflip/amplitude_damping as "unreliable" -- re-investigated rather than trusted, and confirmed to be a Monte Carlo undersampling artifact (extrapolation coefficients amplify input noise; an undersampled estimate makes the corrected result noisy even when the correction itself is sound), not a real limitation.
- More noise-scale points, SAME total measurement budget (the fair
comparison -- splitting a fixed number of trajectories across more
points, not spending more): 3 points x K=400 (1200 total) vs. 5
points x K=240 (1200 total) vs. 7 points x K=171 (~1200 total),
n=4 qubits, all 5 noise channels, 5 seeds. 5 points matches or
slightly beats the 3-point mean delta (+0.150 vs +0.148) with 19%
lower variance (std 0.050 vs 0.062) -- a real, free improvement at
the same experimental cost, not an artifact of spending more. 7
points trades a little mean (+0.132) for still-lower variance (std
0.043, 30% below baseline) -- a genuine tradeoff point, useful when
reliability matters more than average performance. At exactly 3
points this function is mathematically identical to
richardson_extrapolate(verified to 1e-12) -- there is no free lunch there, the gain only appears once more points are used. - More noise-scale points, FIXED K per point instead (spending more
total measurement, K=400 at every point count): with exact
interpolation this makes things worse, not better (mean delta drops
from +0.148 at 3 points to +0.081 at 5, and to -0.220 at 5
closely-spaced points -- see
polynomial_extrapolate's docstring); with this function's degree=2 default it instead stays comparable in mean with lower variance (std 0.062 -> 0.035-0.046) -- confirms the safety property holds even when not holding budget fixed, on top of the fixed-budget gain above.
Practical implication for callers regardless of degree: correction
quality still depends on rho_at_scales being a reasonably low-noise
estimate to begin with (large enough K, or equivalent) -- an
extrapolation fit through pure noise cannot recover signal that isn't
there. degree trades bias for variance: higher degree fits the true
curve's shape more closely (less bias) but is more sensitive to
per-point noise (more variance); degree=2 was chosen empirically as the
best tested tradeoff, not a theoretical optimum for every regime.
Do not pass a target/ideal density matrix into this function or use one
to pick among candidate corrections -- see uhlmann_fidelity's
docstring for why.
Source code in dense_evolution/mitigation/zne.py
608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 | |
jsd_predictive_zne_density_matrix ¶
Density-matrix ZNE with a Jensen-Shannon-divergence-informed
coefficient nudge, for noise whose scale-to-output-distribution
relationship isn't perfectly smooth (the assumption plain
3-point Richardson extrapolation, which zne_density_matrix
defaults toward at exactly 3 scales, relies on).
Motivated by and validated in Dense-Evolution-Discovery's
scripts/photonic_predictive_zne.py -- prototyped there first for
photon-loss noise (a photonic-relevant channel: photon loss on a
dual-rail-encoded qubit IS this library's amplitude_damping
channel), per this project's cross-repo promotion pattern. Grounded
in real literature: Mills & Mezher, "Mitigating photon loss in
linear optical quantum circuits" (arXiv:2405.02278), find plain
scalar ZNE does not beat postselection for discrete-variable photon
loss -- reproduced directly there (scalar ZNE went unphysical,
fidelity > 1.0, at 14/16 swept points). zne_density_matrix avoids
that failure mode by construction (project_to_physical); this
function asks whether a further, data-driven adaptive correction on
top of it can do even better.
Signal: Jensen-Shannon divergence (_js_divergence, the standard
formula -- see its own docstring for how this differs from
dense_evolution.mps's unrelated "adaptive JSD") between the
measurement-probability distributions (density-matrix diagonals) at
consecutive noise scales. Needs no external calibration or oracle
access to an ideal/target state -- unlike naively reusing
calculate_delta_preemp with an externally-supplied signal (tried
first in the Discovery prototype; found to have a negligible effect
by construction, since that formula's fixed nudge constants 0.01/
0.02 were tuned for a differently-scaled use case elsewhere in this
module).
nonlinearity = (jsd_23 - jsd_12) / (jsd_23 + jsd_12 + eps)
(bounded in [-1, 1]) measures how consistently the JSD grows between
consecutive scales -- near 0 when the noise-scale -> output-
distribution map is locally well-behaved (Richardson's implicit
assumption holding), away from 0 when it isn't. RECTIFIED: the
coefficient nudge is applied only when nonlinearity > 0 -- an
unrectified first version, applying the nudge for both signs,
helped in only 5/16 points on a real run despite the signal itself
being significantly correlated with success (Pearson r=+0.533,
p=0.0334); the fix was clipping to the regime the signal was shown
to work in, not discarding the signal. When nonlinearity <= 0,
this reduces EXACTLY to zne_density_matrix at 3 equally-spaced
scales (verified: max deviation ~1e-8, floating-point noise) --
zero risk in that regime, by construction.
Verified on a real, seed-diverse sample (72 points: 12 photon-loss
rates x 6 independent seeds, K=200 trajectories each) before being
promoted here, not just the small sample that first suggested it:
among 46 points where the mechanism actually activates
(nonlinearity > 0.01), 76.1% (35/46) improve over plain
zne_density_matrix, mean fidelity gain +0.0055, one-sample
t-test against zero p=0.0003 -- and positive in 6/6 independent
seeds tested (not one lucky seed driving the result). The win rate
and effect size were LARGER on the big sample than the small one
that first suggested it, the opposite of the usual small-sample-
regresses-to-null pattern -- checked directly rather than assumed
either way before trusting it.
Only defined for exactly 3 equally-spaced noise factors (1x, 2x,
3x), same restriction as zero_noise_extrapolation's own healing-
adapted branch and for the same reason: the underlying Lagrange
coefficients (3, -3, 1) this nudges are specific to that spacing,
not a general n-point formula.
Source code in dense_evolution/mitigation/zne.py
Standalone density-matrix noise channels¶
Three CPTP channels usable directly on a density matrix, independent of NoiseModel's
per-qubit gate-noise pipeline (circuits.registry). global_depolarizing_channel is
symmetric -- it mixes the whole register toward the fully-mixed state as one unit -- promoted
from a real reproduction of arXiv:2608.16716's SPAM model
(Experiment 33).
amplitude_damping_channel is the opposite kind of asymmetric: population only ever moves
|1>→|0>, never the reverse -- the real signature of T1 decay and of quasiparticle
poisoning, promoted from
Experiment 34's
reproduction of a real cosmic-ray-induced error burst (arXiv:2104.05219). cosmic_ray_burst_profile
is not a channel itself but the time-dependent decay-probability GENERATOR that experiment's
real numbers were extracted into a reusable, parametrized form from -- feed its output straight
to amplitude_damping_channel via continuous_dissipative_evolve. Both are
already covered by the dense_evolution.mitigation.zne API reference above -- this section
is context, not a duplicate listing.
Density-matrix diagnostics¶
Two further density-matrix diagnostics, both originated as Colab proposals with real bugs,
fixed and validated in Dense-Evolution-Discovery
before promotion here: a non-commuting-aware divergence (sandwiched_renyi_divergence) and a
single-qubit non-stabilizerness measure (magic_entropy). Both are validation-only, like
uhlmann_fidelity above -- meant to grade a correction against a known reference state, not to
feed into one.
renyi ¶
Sandwiched Quantum Renyi Divergence for full density-matrix diagnostics (Muller-Lennert, Reeb, Wolf, Wilde, "On quantum Renyi entropies: a new generalization and some applications", arXiv:1306.3142, Definition 1).
D_alpha(rho||sigma) = 1/(alpha-1) * log2 Tr[(sigma^e rho sigma^e)^alpha], e = (1-alpha)/(2*alpha), with the alpha->1 limit reducing to the standard quantum relative entropy and alpha=1/2 reducing to a fidelity-based form.
Originated from a Colab proposal with a real bug in its case_general
branch: tr_inner = jnp.maximum(tr_inner, 1.0) floors the inner trace at
1.0 even when the true value is < 1 (the normal case for non-commuting
rho, sigma), silently forcing every result to log2(1)=0 -- confirmed
directly in the Colab's own printed output (alpha=1.5 gave exactly
0.000000 across an entire rotation sweep). A second, deeper bug survived
the floor-value fix alone: for alpha > 1, a trace below 1 is not a
numerical artifact to clamp away, it is the genuine signature of a
support mismatch (supp(rho) not contained in supp(sigma)), which the
divergence must report as +inf, not a finite (and wrong-signed) number --
verified by hand on two different pure states: Tr[Q^1.5] = 0.6759,
matching the closed-form prediction (|
Fixed and validated in Dense-Evolution-Discovery, Experiment 29 (https://tatopenn-cell.github.io/Dense-Evolution-Discovery/sandwiched_renyi_density_matrix/): against the alpha->1 relative-entropy limit (matches an independent numpy reference to 4 decimal places), the commuting/diagonal case (reduces exactly to the classical Renyi divergence), and the support-violation +inf case (verified at alpha>1, confirmed the branch does not fire spuriously at alpha<1).
Its originally proposed use case -- replacing the JSD-based truncation criterion in dense_evolution.mps's bond-dimension search -- was independently disproven: on the diagonal singular-value spectrum used there, rho and sigma commute, so a non-commuting-aware divergence induces the exact same truncation ordering as JSD (5 benchmark configurations, byte-identical chi_used and truncation error every time) -- nothing for it to add in that setting. Promoted here instead for the genuinely non-commuting full-density-matrix diagnostic use case it WAS validated against: alongside uhlmann_fidelity, tracking a Bell state degraded by amplitude damping, where the two metrics' noise-sensitivity curves visibly diverge from each other.
sandwiched_renyi_divergence_jit
module-attribute
¶
jax.jit-compiled entry point for sandwiched_renyi_divergence. rho/
sigma must already be complex128. Returns a jnp scalar, not a Python
float.
sandwiched_renyi_divergence ¶
Sandwiched quantum Renyi divergence D_alpha(rho||sigma), in bits
(log2). rho, sigma are density matrices of the same dimension;
alpha selects the order (0.5 -> fidelity-based, 1.0 -> standard
relative entropy, both handled as exact closed-form limits rather
than through the general formula's own alpha->0.5/1 numerical
instability).
Zero when rho == sigma at every alpha (verified). For alpha > 1,
returns +inf when supp(rho) is not contained in supp(sigma) --
e.g. two different pure (rank-1) states -- rather than a finite
number; see the module docstring for why this is the mathematically
correct behavior, not an edge-case failure.
KNOWN LIMITATION at exactly alpha=1.0: the same support-violation
check is NOT applied to the alpha=1 (relative-entropy) branch, which
instead clips log(0)-type contributions to 0 rather than diverging --
e.g. D_1(rho||sigma) for two different pure states returns 0.0, not
+inf, even though the true relative entropy diverges there too.
Experiment 29 validated the alpha=1 branch only against full-rank
(depolarized) inputs specifically to sidestep this exactly-singular
case (scipy.linalg.logm itself raises LogmExactlySingularWarning
on singular inputs -- an inherent ill-conditioning of relative
entropy near degenerate support, not unique to this implementation).
Do not rely on alpha=1 to correctly flag a support mismatch; use
alpha slightly above 1 (e.g. 1.001) if that matters for your use case.
Validation-only, like uhlmann_fidelity: meant to grade a correction
against a known reference state, not to feed into one (see
uhlmann_fidelity's docstring for the full "ideal state as oracle"
argument, which applies here identically).
Examples:
>>> import numpy as np
>>> from dense_evolution.mitigation.renyi import sandwiched_renyi_divergence
>>> rho = np.array([[1, 0], [0, 0]], dtype=complex)
>>> round(float(sandwiched_renyi_divergence(rho, rho, alpha=1.5)), 4)
0.0
>>> sigma = np.array([[0.5, 0], [0, 0.5]], dtype=complex)
>>> round(float(sandwiched_renyi_divergence(rho, sigma, alpha=1.5)), 4)
1.0
Source code in dense_evolution/mitigation/renyi.py
magic_entropy ¶
Magic entropy: a single-qubit density-matrix diagnostic built from the 3-fold self-convolution "Key Unitary" construction (Bu, Gu, Jaffe, "Stabilizer testing and magic entropy", arXiv:2306.09292, Definitions 7-8).
Originated from a Colab proposal for a pairwise "Quantum Ruzsa Divergence" (following Bu, Gu, Jaffe, "A convolutional quantum Ruzsa divergence and its applications", arXiv:2401.14385) that turned out to have no valid definition for qubits: that paper's pairwise convolution needs s^2+t^2=1 mod d, which has no solution at d=2. The companion paper above does not patch this with a qubit-specific pairwise formula -- it defines a structurally different, minimum-3-input "Key Unitary" convolution (K quantum registers, K must be ODD, K>=3; there is no K=2 case). For qubits the smallest valid object is therefore the 3-fold SELF-convolution of one state with itself, boxtimes_3(rho,rho,rho), and the entropy of its reduced output register is what the paper calls "magic entropy": zero for stabilizer states, positive for non-stabilizer ("magic") states (the paper's Examples 32/33).
Validated end-to-end in Dense-Evolution-Discovery, Experiment 30 (https://tatopenn-cell.github.io/Dense-Evolution-Discovery/quantum_ruzsa_magic_entropy/): the Key Unitary circuit was checked basis-state by basis-state against the paper's own Lemma 9 combinatorial identity; magic_entropy was checked against all six single-qubit stabilizer states (~0, max observed 4e-11) and the two standard magic states T and H (0.811 bits, matching each other exactly as expected by symmetry); and it was used as a noise diagnostic that is qualitatively distinct from uhlmann_fidelity and the sandwiched Renyi divergence (see renyi.py in this same subpackage) -- under amplitude damping it is non-monotonic (rises then returns to exactly 0, since the p=1 fixed point |0> is a stabilizer state), unlike either of those two, which change monotonically over the same sweep.
Restricted to SINGLE-QUBIT density matrices (2x2) -- the Key Unitary here is built for n=1-qubit registers specifically; a multi-qubit generalization would need a larger Key Unitary circuit (n-qubit registers, Definition 7 in general) not implemented here.
A shadow-measurement-based estimator for this same quantity, using randomized measurement snapshots instead of the exact density matrix, is promoted alongside this module in magic_entropy_shadows.py -- see that module for why it has its own API shape (measurement snapshots in, not a density matrix) rather than a function alongside this one.
magic_entropy_jit
module-attribute
¶
jax.jit-compiled entry point for magic_entropy. rho must already
be complex128. Returns a jnp scalar, not a Python float -- call
float(...) yourself if you need one outside a jitted context.
magic_entropy ¶
Magic entropy of a single-qubit density matrix rho (2x2), in bits
(log2) -- NOTE this differs from
dense_evolution.physics.entropy.von_neumann_entropy's natural-log
(nats) convention; kept as log2 here to match both the source paper's
own convention and the values already published in
Dense-Evolution-Discovery's Experiment 30.
Zero for every single-qubit stabilizer state (|0>, |1>, |+>, |->, |+i>, |-i>), positive for non-stabilizer ("magic") states -- e.g. the T-state and H-state both give 0.811 bits.
Differentiable through jax.grad, including at stabilizer states
where the reduced matrix's eigenvalues are exactly degenerate (e.g.
the fully mixed state I/2 gives eigenvalues [0.5, 0.5]): unlike
uhlmann_fidelity, which needs a custom _eigh_degenerate_safe JVP
rule because it reconstructs eigenVECTORS (ill-defined in a
degenerate eigenspace), this function only ever needs eigenVALUES
(jnp.linalg.eigvalsh, no eigenvectors), whose gradient is
well-defined even at exact degeneracies -- confirmed directly:
jax.grad(magic_entropy) is finite (no NaN) at both the fully mixed
state and a magic state.
Source code in dense_evolution/mitigation/magic_entropy.py
Multi-qubit magic monotone (pure states)¶
A second, unrelated magic quantity: stabilizer_renyi_entropy (Leone, Oliviero, Hamma,
arXiv:2106.12587) is zero for every stabilizer state
and positive otherwise, like magic_entropy above -- but it works on a MULTI-qubit
pure statevector (magic_entropy is single-qubit density-matrices only), and it is a
genuinely different construction (a Walsh-Hadamard transform of Pauli-string overlaps,
not the 3-fold self-convolution "Key Unitary" magic_entropy uses). Promoted from
Dense-Evolution-Discovery's wormhole_magic_entropy.py,
where it tracked how "magic" a wormhole-teleportation state stayed across the protocol,
regardless of how well the teleportation itself worked.
import numpy as np
from dense_evolution.mitigation import stabilizer_renyi_entropy
ghz = np.zeros(8, dtype=complex)
ghz[0] = ghz[-1] = 1.0 / np.sqrt(2.0)
round(stabilizer_renyi_entropy(ghz), 6)
stabilizer_renyi_entropy ¶
Stabilizer Renyi Entropy (SRE): a per-STATE nonstabilizerness ("magic") monotone (Leone, Oliviero, Hamma, "Stabilizer Renyi Entropy", arXiv:2106.12587, Phys. Rev. Lett. 128, 050402 (2022), Eq. 5-8 there, labeled Eq. 14/18 in the paper that motivated promoting this).
NOT the same quantity as dense_evolution.mitigation.magic_entropy
(Bu-Gu-Jaffe's 3-fold self-convolution "Key Unitary" construction,
single-qubit only) or sandwiched_renyi_divergence (Muller-Lennert et
al., a DIVERGENCE between TWO density matrices -- "how different are rho
and sigma", answering a different question entirely; the shared "Renyi"
in both names is a coincidence of both being alpha-generalizations of
entropy applied to different objects, not overlapping math). This SRE is
a genuinely different, MULTI-qubit, SINGLE-state magic monotone: zero for
every stabilizer state, positive otherwise.
M_2(psi) = -log2[ (1/d) * sum_a sum_b |WHTc_a|^4 ], where c_a(x) = conj(psi(x)) * psi(x XOR a) and WHT is the length-d Walsh-Hadamard transform (signmat[b,x] = (-1)^popcount(b AND x)), d = 2**n_qubits.
Verified against known values: every computational-basis/stabilizer state gives exactly 0; a single T state ((cos(pi/8), sin(pi/8)) in the computational basis) gives -log2(0.75) = 0.415037 bits, matching the closed-form derivation of Eq. 5 by hand (not a value copied from elsewhere).
Promoted from Dense-Evolution-Discovery's wormhole_magic_entropy.py (2026-08-29), where it was implemented fresh because the existing magic_entropy is single-qubit only -- this quantity has no dependency on that use case (wormhole teleportation), so it belongs here as a general multi-qubit magic diagnostic.
VECTORIZATION NOTE: the Discovery script's original version had an
explicit Python for a in range(d) loop, each iteration doing its own
(d,d)@(d,) matrix-vector product -- d sequential small matmuls. This
promoted version instead builds the (d,d) matrix of every c_a(x) pair at
once (broadcasting over x and a together) and applies the Walsh-Hadamard
transform as ONE (d,d)@(d,d) matmul -- same O(d^3) total FLOP count, but
expressed as a single large matmul JAX/XLA can execute efficiently
(GPU/TPU-friendly, no per-iteration Python dispatch overhead), matching
this package's JAX-by-default convention. Still O(d^3), not the paper's
own asymptotically-better O(4^n * n) Walsh-Hadamard butterfly algorithm
(n = log2(d) qubits) -- unimplemented here, same as the Discovery
original; fine for d up to a few thousand (n up to ~11-12 qubits), the
sizes this package's exact-statevector backends already target.
stabilizer_renyi_entropy_jit
module-attribute
¶
jax.jit-compiled entry point for stabilizer_renyi_entropy. psi
must already be complex128. Returns a jnp scalar, not a Python float
-- call float(...) yourself if you need one outside a jitted context.
stabilizer_renyi_entropy ¶
Stabilizer Renyi Entropy of a pure state psi (length 2**n_qubits),
in bits (log2) -- the paper's own convention.
Zero for every stabilizer state, positive for non-stabilizer ("magic") states -- e.g. a single T state gives 0.415037 bits.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
psi
|
(array - like, shape(2 ** n_qubits))
|
A normalized pure statevector. |
required |
Returns:
| Type | Description |
|---|---|
float
|
|
Examples:
>>> import numpy as np
>>> psi0 = np.zeros(8, dtype=complex); psi0[0] = 1.0
>>> round(stabilizer_renyi_entropy(psi0), 6) # computational basis state: stabilizer, expect 0
0.0
Source code in dense_evolution/mitigation/stabilizer_renyi_entropy.py
Shadow-based estimation¶
A classical-shadows-based estimator for magic_entropy above, estimating it from randomized
measurement snapshots instead of the exact density matrix. Different API shape from everything
else on this page -- sampling (sample_classical_shadow) and estimation
(magic_entropy_from_shadows) are separate steps, since shadow data can come from this
simulator's own Born-rule oracle sampling or, in principle, real hardware measurement logs
reconstructed the same way. Not jax.jit-compatible (median-of-means uses numpy.median).
magic_entropy_shadows ¶
Classical-shadows-based estimator for magic_entropy (see
magic_entropy.py in this same subpackage) -- estimates the same
single-qubit non-stabilizerness quantity from randomized measurement
snapshots instead of requiring the exact density matrix.
Originated from a Colab proposal for a dense_evolution/circuits/shadows.py
module (following Huang, Kueng, Preskill 2020, "Predicting Many Properties
of a Quantum System from Very Few Measurements") with a real bug in its
purity estimator (a missing transpose in a U-statistic einsum contraction,
silent whenever every snapshot happened to be real-valued). Fixed, then
extended -- using the same multi-copy U-statistic trick the paper says
"readily generalizes to higher order polynomials" -- to estimate
magic_entropy's reduced convolution matrix from shadow snapshots instead
of the exact rho. Matured across Dense-Evolution-Discovery Experiment 31
(three real gaps found and closed there before promotion: the purity
estimator bug, missing median-of-means robustness, and no
sample-complexity guidance):
https://tatopenn-cell.github.io/Dense-Evolution-Discovery/quantum_shadows_magic_entropy/
API SHAPE differs from every other function in this subpackage: sampling
(sample_classical_shadow) and estimation (magic_entropy_from_shadows)
are separate steps, matching how classical shadows work in general -- the
snapshot data can come from this simulator (sample_classical_shadow uses
oracle access to rho's exact Born-rule probabilities, something only a
simulator has) or, in principle, from real hardware measurement outcomes
reconstructed the same way (rho_hat = 3 U^dagger |b><b| U - I per
snapshot, from a recorded basis+outcome).
Restricted to SINGLE-QUBIT density matrices, matching magic_entropy's
own scope.
sample_classical_shadow ¶
Simulates the real single-qubit random-Pauli classical-shadow
measurement protocol against a known rho (2x2, pure or mixed): for
each of n_snapshots independent draws, picks a random Pauli basis
uniformly, samples a computational-basis outcome from the true
Born-rule probability under that basis (this simulator has oracle
access to rho, unlike real hardware), then reconstructs the
classical snapshot rho_hat = 3 U^dagger |b><b| U - I.
Returns an (n_snapshots, 2, 2) complex128 array -- feed this
directly into magic_entropy_from_shadows.
Each individual rho_hat is NOT a valid density matrix on its own
(can have negative eigenvalues) -- only the average over many
snapshots converges to the true rho. Verified directly in
Dense-Evolution-Discovery Experiment 31: the empirical mean over
200,000 snapshots matched the true state to within 0.004.
Source code in dense_evolution/mitigation/magic_entropy_shadows.py
magic_entropy_from_shadows ¶
Estimates magic_entropy(rho) from classical shadow snapshots of
rho (from sample_classical_shadow, or real hardware measurement
data reconstructed the same way) instead of the exact density matrix.
Groups the snapshots into disjoint triples, estimates each entry of
the 3-copy self-convolution's reduced matrix R via median-of-means
over Tr[O_ab . (rho_hat_i (x) rho_hat_j (x) rho_hat_k)] (unbiased,
since each triple's three snapshots are independent unbiased
estimators of rho; real and imaginary parts of each entry are
median-of-means'd separately, the standard practical choice since a
complex median has no single definition), then computes the von
Neumann entropy of the (Hermitized, eigenvalue-clipped, trace-
renormalized) ESTIMATED R classically -- entropy itself is never
shadow-estimated directly, matching how Huang et al. handle their own
Renyi-2 entanglement entropy example.
Not jax.jit-compatible (unlike every other function in this
subpackage): median-of-means uses numpy.median, which has no
equivalent JAX primitive at this scale.
See approx_shadow_std/fit_shadow_sample_complexity for how many
snapshots this needs for a given error tolerance.
Source code in dense_evolution/mitigation/magic_entropy_shadows.py
approx_shadow_std ¶
Rough approximate standard deviation (bits) of
magic_entropy_from_shadows's estimate at a given snapshot count,
from an empirical fit (not a formal theorem) calibrated on a |T>
state in Dense-Evolution-Discovery Experiment 31: 20 independent
trials at each of 4 snapshot counts (3,000-100,000), a log-log linear
regression gave std(n) ~ 11.75 / n^0.546 -- the fitted exponent
(0.546) is close to the ~0.5 ("error shrinks like 1/sqrt(n)") standard
shadow/median-of-means theory predicts.
This is a quick sanity-check fallback, not a guarantee for an
arbitrary state -- call fit_shadow_sample_complexity on YOUR
specific state if you need a real, state-calibrated error bound.
Source code in dense_evolution/mitigation/magic_entropy_shadows.py
fit_shadow_sample_complexity ¶
fit_shadow_sample_complexity(
rho: ndarray,
exact_value: float,
n_snapshots_list,
n_trials: int,
seed_base: int = 0,
)
Empirically measures magic_entropy_from_shadows's standard
deviation across n_trials independent shadow samplings at each
snapshot count in n_snapshots_list, for the SPECIFIC state rho
(rather than trusting approx_shadow_std's built-in T-state-derived
fallback), then fits std(n) ~ C / n^p via log-log linear regression
-- the same method used to derive approx_shadow_std's constants in
the first place (Dense-Evolution-Discovery Experiment 31).
exact_value should be magic_entropy(rho) -- used only to also
report each snapshot count's mean estimation bias alongside the
fitted curve, not part of the fit itself.
Returns (rows, fit_c, fit_p): rows is a list of per-snapshot-count
dicts (n_snapshots, mean_estimate, std_estimate,
mean_abs_error); fit_c/fit_p are the fitted constants for
C / n^p, usable the same way as approx_shadow_std (or pass them to
approx_shadow_std's formula directly: fit_c / n ** fit_p).
Source code in dense_evolution/mitigation/magic_entropy_shadows.py
Classical distribution divergence¶
The classical Kullback-Leibler divergence over probability distributions (Kullback & Leibler,
1951) -- distinct from sandwiched_renyi_divergence above, which operates on density
matrices via matrix logarithms; this operates directly on probability vectors (e.g. a
measurement-outcome distribution jnp.abs(psi) ** 2), no eigendecomposition needed. Additive
to dense_evolution.healing's existing scalar log-ratio signal, not a
replacement for it -- validated in
Dense-Evolution-Discovery, Experiment 32
to be a genuinely different signal on the same states, not a rescaling.
kl_divergence ¶
Classical Kullback-Leibler divergence between probability distributions (Kullback, S. & Leibler, R.A., "On Information and Sufficiency", The Annals of Mathematical Statistics, 22(1), 79-86, 1951).
D_KL(p||q) = sum_x p(x) * log2(p(x)/q(x)), the relative entropy of q from
p, in bits (log2, matching the rest of this subpackage's convention --
sandwiched_renyi_divergence/magic_entropy both use log2).
Checked against the paper's own text directly (Section 2, eq. 2.2-2.3),
not assumed from the textbook formula alone: what this module implements
is what Kullback & Leibler call I(1:2), "the mean information for
discrimination between H1 and H2" -- what the broader literature later
popularized as "the KL divergence". Their OWN word "divergence",
J(1,2) = I(1:2) + I(2:1) (eq. 2.9), names the symmetrized sum of both
directions instead -- deliberately not implemented here, since it would
duplicate the Jensen-Shannon divergence this codebase already uses
(mps.py's bond-dimension search, zne.py's predictive ZNE), which is
bounded and better-behaved at disjoint supports.
Distinct in kind from sandwiched_renyi_divergence(rho, sigma, alpha=1.0),
which reduces to the QUANTUM relative entropy Tr[rho(log rho - log sigma)]
between density MATRICES via matrix logarithms -- this module implements
the simpler classical case directly over probability VECTORS (e.g. a
measurement-outcome distribution |psi|^2, or any other normalized
histogram), with no eigendecomposition needed.
Originated from an honest gap flagged in this subpackage's own healing.py docstring: calculate_vettore_dinamico's core term, log(E_B/E_A), is a single un-weighted log-likelihood ratio between two scalars -- the same elementary quantity this divergence is built from, but not this divergence itself. Built and validated in Dense-Evolution-Discovery, Experiment 32 (https://tatopenn-cell.github.io/Dense-Evolution-Discovery/kullback_leibler_divergence/): against an independent scipy.stats.entropy reference (1e-9 bits across 20 random trials), Gibbs' inequality (D_KL >= 0, 200 random pairs, never negative), a genuine support-violation case (+inf, not a finite wrong number), and a real measurement-distribution application confirming this is not a trivial rescaling of healing.py's existing scalar signal.
Additive, not a replacement for the already-validated healing pipeline.
kl_divergence_jit
module-attribute
¶
jax.jit-compiled entry point for kl_divergence. p/q must already
be float64 arrays. Returns a jnp scalar, not a Python float.
kl_divergence ¶
Classical Kullback-Leibler divergence D_KL(p||q), in bits.
p, q are 1-D real, non-negative probability vectors of the same
length, each summing to 1 (not validated here -- callers pass in a
normalized distribution, e.g. jnp.abs(psi) ** 2 for a statevector's
measurement-outcome probabilities). Not symmetric: D_KL(p||q) !=
D_KL(q||p) in general.
Zero iff p == q (Gibbs' inequality: D_KL(p||q) >= 0 always, with
equality only at p == q). Returns +inf when p has support where q
does not (p(x) > 0, q(x) == 0 for some x) -- the correct value, not an
edge case to avoid; see the module docstring.
Source code in dense_evolution/mitigation/kl_divergence.py
See also: dense_evolution.healing for the predictive-healing primitives
(calculate_delta_preemp) the healing-adapted extrapolation branch is built on, and
NoiseModel for the Kraus-channel noise used to build the noisy ensembles
these functions correct. Full worked example: Density-matrix ZNE healing.