docs/examples.md

Examples

Three runnable, end-to-end examples, each lifted from real experiments/tests already in the repository rather than written fresh for this page:

• Density-matrix ZNE healing — from experiments/matrix_healing_zne.py, the exact script zne_density_matrix's docstring reports measured numbers from.
• MPS for low-entanglement circuits — from test_mps.py::test_run_circuit_jit_ghz_chain.
• Differentiable VQE — from test_autodiff.py::TestCircuitToEnergyFn.

Density-matrix ZNE healing

Zero-noise extrapolation on full density matrices: run a circuit at several scaled noise levels, extrapolate back to the zero-noise limit, then project the (generally unphysical) extrapolated result back onto the physical cone of valid density matrices (Smolin-Gambetta-Smith 2012).

The ideal state (rho_ideal) is built once here purely to grade the result at the end. It is never fed into the noise ensemble, the extrapolation, or the physical-projection step — feeding it back in anywhere but the final fidelity check would be oracle access, not error mitigation.

import numpy as np
import jax.numpy as jnp
import dense_evolution as de
from dense_evolution.registry import NoiseModel
from dense_evolution.mitigation import uhlmann_fidelity, zne_density_matrix

N_QUBITS = 2
BASE_P = 0.05
SCALES = (1.0, 2.0, 3.0)
K_TRAJECTORIES = 200


def bell_state_sv():
    sim = de.DenseSVSimulator(N_QUBITS)
    sim.run_circuit([("h", 0), ("cx", 0, 1)])
    return np.asarray(sim.get_statevector())


def noisy_density_matrix(ideal_sv, p, k, rng):
    dim = len(ideal_sv)
    rho = np.zeros((dim, dim), dtype=np.complex128)
    for _ in range(k):
        sv_noisy = NoiseModel.apply_to_sv(ideal_sv.copy(), N_QUBITS, 'depolarizing', p, rng=rng)
        rho += np.outer(sv_noisy, sv_noisy.conj())
    rho /= k
    return jnp.asarray(rho, dtype=jnp.complex128)


ideal_sv = bell_state_sv()
rho_ideal = jnp.asarray(np.outer(ideal_sv, ideal_sv.conj()), dtype=jnp.complex128)

rng = np.random.default_rng(0)
rho_at_scales = jnp.stack([
    noisy_density_matrix(ideal_sv, BASE_P * scale, K_TRAJECTORIES, rng)
    for scale in SCALES
])

raw_fidelity = uhlmann_fidelity(rho_at_scales[0], rho_ideal)          # base-scale, uncorrected
corrected = zne_density_matrix(rho_at_scales, SCALES)                 # Richardson + physical projection
corrected_fidelity = uhlmann_fidelity(corrected, rho_ideal)           # grading only

print(f"raw fidelity:       {raw_fidelity:.4f}")
print(f"corrected fidelity: {corrected_fidelity:.4f}")

Measured on this exact script (5 seeds, 400 trajectories each, 2–5 qubits, 5 noise channels): 96/100 runs improve fidelity, mean delta +0.12. See the Changelog and dense_evolution.mitigation docstrings for the full sweep and the honest negative results (predictive-healing coefficient perturbation was tested and rejected for the matrix case — negligible effect even amplified 100x).

Run it yourself:

python experiments/matrix_healing_zne.py

MPS for low-entanglement circuits

MPSSimulator keeps a bounded bond dimension (max_bond), trading exactness for circuits that stay low-entanglement (GHZ chains, shallow QAOA layers, most NISQ ansätze) at qubit counts a dense statevector could never hold. run_circuit_jit fuses the whole circuit into a single jax.lax.scan.

import numpy as np
from dense_evolution import MPSSimulator

n = 6
ops = [["h", 0]] + [["cx", q, q + 1] for q in range(n - 1)]  # GHZ chain

mps = MPSSimulator(n_qubits=n, max_bond=8)
mps.run_circuit_jit(ops)

prob = np.abs(np.asarray(mps.contract_to_statevector())) ** 2
print(prob[0], prob[2 ** n - 1])  # both ~0.5, everything else ~0 -- the GHZ signature

See dense_evolution.mps for the truncated-SVD mechanics and the jax.lax.scan fusion this relies on.

Differentiable VQE

circuit_to_energy_fn turns a parsed circuit into a (theta, hamiltonian) -> (energy, statevector) function that is differentiable end-to-end with jax.grad/jax.value_and_grad — no manual parameter-shift rule needed. This runs a short Adam loop against a random diagonal Hamiltonian.

import numpy as np
import jax
import jax.numpy as jnp
import dense_evolution as de

VQE_QASM = (
    'OPENQASM 2.0; include "qelib1.inc"; qreg q[2]; creg c[2]; '
    'ry(0.5) q[0]; rx(0.5) q[1]; cx q[0],q[1]; rz(0.2) q[1]; cx q[0],q[1]; '
    'ry(0.5) q[0]; rx(0.5) q[1]; measure q -> c;'
)


def random_hamiltonian(n_qubits, seed=7):
    rng = np.random.default_rng(seed)
    values = np.sort(rng.uniform(-2.5, 2.5, 2 ** n_qubits))
    return jnp.diag(jnp.array(values, dtype=jnp.float64))


circ = de.QASMParser().parse(VQE_QASM)
energy_fn, n_params = de.circuit_to_energy_fn(circ, circ.n_qubits)
h_matrix = random_hamiltonian(circ.n_qubits)

rng = np.random.default_rng(3)
theta = rng.uniform(-np.pi, np.pi, n_params)

energy_and_grad = jax.jit(jax.value_and_grad(energy_fn, argnums=0, has_aux=True))
m, v = np.zeros(n_params), np.zeros(n_params)
lr, beta1, beta2 = 0.1, 0.9, 0.999

for epoch in range(1, 31):
    (energy, _), grad = energy_and_grad(jnp.asarray(theta), h_matrix)
    grad = np.asarray(grad)
    m = beta1 * m + (1 - beta1) * grad
    v = beta2 * v + (1 - beta2) * (grad ** 2)
    m_hat, v_hat = m / (1 - beta1 ** epoch), v / (1 - beta2 ** epoch)
    theta = theta - lr * m_hat / (np.sqrt(v_hat) + 1e-8)

print(f"final energy:  {float(energy):.4f}")
print(f"ground state:  {float(jnp.min(jnp.diag(h_matrix))):.4f}")

circuit_to_energy_fn also accepts an optional noise= (dense_evolution.NoiseSpec, a JAX pytree) to trace noisy VQE runs natively — composable with jax.jit, jax.grad, and jax.vmap over noise-key batches. See dense_evolution.autodiff.