MPS Simulator

DenseSVSimulator keeps every one of a circuit's 2**n amplitudes in memory — that works for a couple dozen qubits, but beyond that it runs out of RAM. MPSSimulator runs the same circuits at hundreds of qubits by keeping a compact tensor-network representation instead of the full state, at the cost of accuracy on highly-entangled circuits.

When to use MPS: low-entanglement circuits — GHZ chains, shallow local circuits, product-state preparations. For highly-entangled circuits the bond dimension grows exponentially and MPS degrades back toward the same cost as the dense engine. It is a complement to DenseSVSimulator, not a replacement.

---

Step 1. Build a circuit

A GHZ chain: qubit 0 goes into superposition, then every qubit is entangled with the next one. This is exactly the kind of circuit MPS is built for — the entanglement stays low no matter how many qubits the chain grows to.

import dense_evolution as de

n = 50
lines = ["OPENQASM 2.0;", 'include "qelib1.inc";',
         f"qreg q[{n}];", f"creg c[{n}];", "h q[0];"]
lines += [f"cx q[{i}],q[{i+1}];" for i in range(n - 1)]
lines.append("measure q -> c;")

circuit = de.QASMParser().parse("\n".join(lines))

50 qubits is written with a loop instead of by hand, but the result is the same real OpenQASM 2.0 text the parser always takes.

---

Step 2. Run it with MPS instead of the full state

A dense statevector at 50 qubits would require 2**50 complex numbers — roughly 16 PB — far more memory than any machine has. MPSSimulator keeps a bounded bond dimension (max_bond) instead, which is enough for a circuit this simple.

mps = de.MPSSimulator(n_qubits=n, max_bond=16)
mps.run_circuit_jit(circuit.to_tuples())

run_circuit_jit compiles the whole circuit into a single jax.lax.scan-fused, @jax.jit-compiled kernel. Use the eager per-gate methods (apply_gate_1q, apply_gate_2q) instead when memory is the priority over speed — the JIT path pads every gamma/lambda tensor to max_bond for the lifetime of the instance.

---

Step 3. Read out the result

get_top_k_probable_states finds the most likely outcomes via greedy beam search, without ever materialising a 2**50-sized array.

idx, probs = mps.get_top_k_probable_states(k=2)
for i, p in zip(idx, probs):
    print(f"{i:0{n}b}: {p:.4f}")

00000000000000000000000000000000000000000000000000: 0.5000
11111111111111111111111111111111111111111111111111: 0.5000

Only two outcomes are populated — all-zeros and all-ones, each at probability 0.5. That is the GHZ signature; everything else is numerically zero.

Note: beam search recall is not guaranteed for any fixed k. If a known state is missing from the output, increase k (e.g. k=128).

---

Step 4. Check accuracy

print(mps.summary())

MPSSimulator | n=50 | chi_max=16 | chi_used=16 | mem=0.198MB | trunc_err=1.61e-06 | avg_JSD=0.0000 | EE_max=1.000b | budget_violations=0

Field	Meaning
chi_max	Hard cap on bond dimension (max_bond argument)
chi_used	Largest bond dimension actually reached during the run
trunc_err	Cumulative singular-value truncation error
avg_JSD	Mean Jensen-Shannon distance between full and truncated singular-value distributions across all SVD steps
EE_max	Peak entanglement entropy across all bonds (bits)
budget_violations	Times max_bond was hit before jsd_budget could be satisfied — if > 0, raise max_bond

Both trunc_err and avg_JSD are near zero here. budget_violations=0 means the accuracy budget was never hit — the result is reliable at this max_bond.

Tip: enable jax_enable_x64 before importing the package for full float64/complex128 precision. On this circuit it drops trunc_err from 1.61e-06 to 3.75e-15 and reduces chi_used from 16 to 2.

import jax
jax.config.update("jax_enable_x64", True)
import dense_evolution as de

---

Step 5. Check that a result has actually converged with bond dimension

chi_used/avg_JSD/budget_violations above describe a single run at a single max_bond — they can't tell you whether a different max_bond would have given a different answer. bond_convergence runs the same circuit at several bond dimensions and checks whether the observable settles down as max_bond grows.

from dense_evolution.backends.mps import bond_convergence

ops = circuit.to_tuples()
result = bond_convergence(ops, n, observables=["Z" + "I" * (n - 1)], bonds=(2, 4, 8))

print(result.verdicts[0])
print(result.diffs[0])

converged
[0.0, 0.0]

For this GHZ chain the entanglement never grows past bond dimension 2, so <Z0> is exactly identical at max_bond=2, 4, and 8 — both successive differences are exactly zero, and the verdict is converged.

A highly-entangled circuit tells a different story: on a 40-qubit, 4-layer brickwall circuit, bonds=(4, 8, 32) gives successive <Z0> differences of about 4.7e-2 then 1.2e-2 — shrinking, but nowhere near a reasonable tol, so the verdict is not_converged. bond_convergence needs at least 3 bond values to make that call at all: two values alone give a single difference, with no way to tell whether it is closing in on tol or has already stalled (see the function's own docstring for the exact numbers). If even the largest bond dimension in bonds is still hitting its own cap, the verdict is undecidable instead of converged or not_converged — no tolerance can be certified from data where the truncation never had room to breathe.

---

Details

Troubleshooting

Problem	Likely cause	Fix
budget_violations > 0 or avg_JSD is high	max_bond is too small: the bond-dimension search hits the cap before jsd_budget is satisfied	Increase max_bond (e.g. 16 → 64). Raising jsd_budget only loosens the tolerance — it does not improve accuracy.
contract_to_statevector raises MemoryError	n > 24 is a hard cutoff, not a guideline	Use get_probabilities_sampled or get_top_k_probable_states instead — neither materialises a (2**n,) array.
get_top_k_probable_states misses a known state	Greedy beam search recall is not guaranteed for a fixed k	Increase k (e.g. 32 → 128).
Simulation is slow or memory is high after run_circuit_jit	Tensors are padded to max_bond for the instance lifetime	Reduce max_bond, or use the eager per-gate path for short, low-entanglement circuits.

Internal: bucketed SVD dispatch

run_circuit_jit's 2-qubit SVD step no longer always runs at a fixed max_bond-padded size -- it dispatches to the smallest provably-sufficient bucket size for the real bond dimension at that cut, via jax.lax.switch, inside the same single compiled kernel. No API or behavior change; measured 68.80x-73.96x faster on CPU and ~1.41x faster on GPU (measured correctly through this same API, not a standalone reimplementation -- an earlier 2.74x GPU claim here was wrong, see the correction below) on a real N=50 TFIM Trotter circuit. Validated first in Dense-Evolution-Discovery's bucketed-SVD experiment and its GPU timing correction, including a real bug (an under-sized bucket could silently drop genuine Schmidt weight already on the bond between the two gated qubits) found and fixed before promotion here.

Optional: gate blocking (fuse_gates=True)

run_circuit_jit(ops, fuse_gates=True) fuses consecutive gates acting on the same (or a growing) qubit pair into one matrix on the host before compiling -- exact, no approximation -- cutting the number of scan steps. Measured ~2x faster than the bucketed dispatch alone on GPU (~2.77x-2.87x total over the original fixed-size SVD), at the cost of coarser truncation_errors/entanglement_entropy/bond-history bookkeeping (one entry per fused step instead of per original gate) -- default is False, unchanged behavior. Validated in Dense-Evolution-Discovery's gate-blocking redesign, including against non-adjacent-gate (SWAP-chain) and CCX circuits.

Performance

DenseSVSimulator and Chunk both require 2**n × 16 bytes for the statevector alone — around 17 GB at 30 qubits, before any computation. MPSSimulator scales as O(n × max_bond²) instead, so a low-entanglement circuit like the GHZ chain above runs at 50+ qubits on an ordinary machine where the dense engines would raise MemoryError.

Exact wall-clock numbers depend heavily on the machine, circuit depth, and entanglement structure, so none are quoted here as a general result. Run mps.summary() to see the numbers for your own machine and circuit, the same way Step 4 above does.

---

mps

MPSSimulator - Matrix Product State statevector simulator, JAX-backed.

Ported from the "TurboQuant TUREQ MPSSimulator v8.2 MatryoshkaFlash" prototype (private research notebook, never published as part of the dense-evolution package). Two real bugs were found and fixed by independent verification against DenseSVSimulator before this module existed in its current form:

1. The original applied Lloyd-Max quantization to the SVD singular values on every truncation ("PolarQuantizer"). Measured a real ~0.5% Total Variation Distance error against DenseSVSimulator on an 8-qubit entangling test circuit, with ZERO bond-dimension savings to show for it. Dropped entirely -- this module keeps only the plain adaptive SVD truncation (JSD-budget-driven bond dimension, the author's own stopping criterion -- standard SVD truncation, non-standard stopping metric).

2. get_top_k_probable_states (originally "_extract_top_k_paths") picked a single "best" bond index via argmax at each step instead of correctly summing over the bond dimension. Measured 0/8 correct states against the exact contraction on the same test circuit, values off by ~30x. Fixed by propagating the true partial-contraction vector through each bond (matches exactly, to machine precision, on every state it finds) -- but note it's a genuine greedy beam search, not an exact top-k finder: recall of the true top states grows with beam width k but isn't guaranteed complete for any fixed k.

Originally ported in plain numpy (matching the prototype), then converted to jax.numpy so the core tensor contractions (einsum, SVD) run on the same backend as the rest of dense_evolution instead of a second, inconsistent numerics stack. Re-verified against DenseSVSimulator after the conversion -- see test_mps.py.

Uses whatever jax_enable_x64 precision is currently active in the process (does not toggle it itself) -- same convention as DenseSVSimulator/Chunk, which rely on the caller (dashboard_core.py's run_simulation) to set precision, since jax_enable_x64 is a process-wide flag and toggling it locally would leak to unrelated code running later in the same process.

For circuits with LOW entanglement (product states, GHZ/Bell-like chains, shallow local circuits), the bond dimension stays small regardless of qubit count, so this scales to hundreds of qubits where DenseSVSimulator (or Chunk) cannot -- see get_probabilities_sampled and get_top_k_probable_states, neither of which ever materializes a (2**n,)-shaped array. For HIGHLY entangled circuits the bond dimension grows and this degrades back toward the same exponential cost DenseSVSimulator has -- it is not a universal replacement, it is complementary.

MPSSimulator

MPSSimulator(
    n_qubits: int,
    max_bond: int = 64,
    svd_cutoff: Optional[float] = None,
    jsd_budget: float = 1e-05,
    use_float32: Optional[bool] = None,
)

Matrix Product State simulator with adaptive SVD-truncated bond dimension (JSD-budget driven), no lossy post-truncation quantization. JAX-backed core (einsum, SVD).

Parameters:

Name	Type	Description	Default
n_qubits	int		required
max_bond	int		64
svd_cutoff	Optional[float]	outright. None (default) resolves to a value appropriate for the dtype this instance actually runs at: ~1e-12 for complex128, ~1e-6 (roughly 10x float32's own machine epsilon) for complex64 -- a fixed 1e-12 in complex64 sits below that dtype's noise floor, so numerical noise gets counted as real Schmidt weight and chi never shrinks below max_bond regardless of jsd_budget. An explicitly passed value always wins verbatim, never rescaled.	None
jsd_budget	float	full and truncated singular-value distributions at each cut; chi is grown by 1 until satisfied or max_bond is hit.	1e-05
use_float32	bool or None -- None (default) follows the process-wide	jax_enable_x64 flag, same convention as this module always used (see module docstring). True forces complex64 (and the complex64- appropriate svd_cutoff default) even if x64 is enabled. False requests complex128 -- since complex128 arrays don't exist in JAX at all unless the process-wide flag is on, this calls the same lazy ensure_x64() DenseSVSimulator uses, mirroring its own use_float32=False handling. That call is a no-op if precision was already pinned via set_precision(), same deference DenseSVSimulator gives it; the flag is re-read afterward rather than assumed, so dtype/eps stay consistent with whatever precision is really active even in that pinned-False edge case.	None

Source code in dense_evolution/backends/mps.py

def __init__(
    self,
    n_qubits: int,
    max_bond: int = 64,
    svd_cutoff: Optional[float] = None,
    jsd_budget: float = 1e-5,
    use_float32: Optional[bool] = None,
):
    self.n = n_qubits
    self.chi = max_bond
    if use_float32 is None:
        x64_active = jax.config.jax_enable_x64
    elif use_float32:
        x64_active = False
    else:
        ensure_x64()
        x64_active = jax.config.jax_enable_x64
    dtype = jnp.complex128 if x64_active else jnp.complex64
    self.eps = svd_cutoff if svd_cutoff is not None else (1e-12 if x64_active else 1e-6)
    self.jsd_budget = jsd_budget

    self.gammas: List[jnp.ndarray] = []
    self.lambdas: List[jnp.ndarray] = [jnp.ones(1)] * (n_qubits + 1)
    # Real (non-max_bond-padded) bond dimension at every cut, kept in
    # sync by BOTH the eager path (apply_gate_2q, below) and
    # run_circuit_jit -- needed by the bucketed-SVD dispatch in
    # _build_mps_runner to pick a provably-sufficient bucket size
    # without ever inferring it from zero-counting (see that
    # function's own docstring for why that would be unreliable).
    self._real_chi: np.ndarray = np.ones(n_qubits + 1, dtype=np.int64)

    self.truncation_errors: List[float] = []
    self.jsd_per_bond: List[float] = []
    self.entanglement_entropy = np.zeros(max(n_qubits - 1, 0))
    self._bond_history: List[int] = []
    # Counts truncations where max_bond was hit before jsd_budget could
    # be satisfied -- the while loop below exits silently in that case,
    # and avg_JSD (a mean over all steps) can look deceptively low even
    # when the final contracted state is badly wrong (verified: TVD
    # ~0.97 against DenseSVSimulator on an 8-qubit/15-layer entangling
    # circuit with max_bond=2, while avg_JSD read 0.0534).
    self.budget_violations: int = 0

    # Cached compiled closure for run_circuit_jit -- built lazily on
    # first use (self.n/self.chi/self.eps/self.jsd_budget are fixed
    # for this instance's lifetime), never rebuilt per call. Same
    # caching pattern as Chunk.__init__'s self._multi_chunk_runner.
    self._mps_runner = None
    self._fused_mps_runner = None

    for _ in range(n_qubits):
        g = jnp.zeros((1, 2, 1), dtype=dtype)
        g = g.at[0, 0, 0].set(1.0)
        self.gammas.append(g)

apply_gate_1q

apply_gate_1q(gate: ndarray, qubit: int) -> None

O(chi^2) -- updates only Gamma[qubit].

Source code in dense_evolution/backends/mps.py

def apply_gate_1q(self, gate: jnp.ndarray, qubit: int) -> None:
    """O(chi^2) -- updates only Gamma[qubit]."""
    gate = jnp.asarray(gate)
    self.gammas[qubit] = jnp.einsum("ij,ljr->lir", gate, self.gammas[qubit])

apply_gate_2q

apply_gate_2q(gate_2q: ndarray, q1: int, q2: int) -> None

2-qubit gate with adaptive SVD truncation. O(chi^3).

Vidal's full two-site update (prog.txt P0 fix): theta is built from BOTH outer Lambdas (Lambda[q1], Lambda[q2+1]) as well as the middle one, not just the middle one -- so its singular values are the true global Schmidt coefficients at this cut, not an artifact of the local 2-site reduced state. New Gamma tensors are recovered by dividing the outer Lambdas back out (regularized: entries at or below svd_cutoff map to a zero inverse instead of blowing up -- exactly the padded/zero entries in the JIT path's fixed-size arrays, and the trivial size-1 boundary Lambda everywhere else).

Source code in dense_evolution/backends/mps.py

def apply_gate_2q(self, gate_2q: jnp.ndarray, q1: int, q2: int) -> None:
    """2-qubit gate with adaptive SVD truncation. O(chi^3).

    Vidal's full two-site update (prog.txt P0 fix): theta is built from
    BOTH outer Lambdas (Lambda[q1], Lambda[q2+1]) as well as the middle
    one, not just the middle one -- so its singular values are the true
    global Schmidt coefficients at this cut, not an artifact of the
    local 2-site reduced state. New Gamma tensors are recovered by
    dividing the outer Lambdas back out (regularized: entries at or
    below svd_cutoff map to a zero inverse instead of blowing up --
    exactly the padded/zero entries in the JIT path's fixed-size
    arrays, and the trivial size-1 boundary Lambda everywhere else)."""
    gate_2q = jnp.asarray(gate_2q)
    if abs(q1 - q2) != 1:
        self._apply_nonlocal_2q(gate_2q, q1, q2)
        return
    if q1 > q2:
        q1, q2 = q2, q1
        gate_2q = jnp.transpose(gate_2q, (1, 0, 3, 2))

    lam_l = self.lambdas[q1]
    g1 = self.gammas[q1]
    lam_m = self.lambdas[q2]
    g2 = self.gammas[q2]
    lam_r = self.lambdas[q2 + 1]

    theta = jnp.einsum("l,lik,k,kjr,r->lijr", lam_l, g1, lam_m, g2, lam_r)
    chiL, d1, d2, chiR = theta.shape

    theta_new = jnp.einsum("abcd,ecdf->eabf", gate_2q, theta)
    theta_mat = theta_new.reshape(chiL * d1, d2 * chiR)

    U_t, S_t, Vh_t, trunc_err, jsd_val = self._svd_truncate(theta_mat)
    chi_new = len(S_t)

    lam_l_inv = jnp.where(lam_l > self.eps, 1.0 / lam_l, 0.0)
    lam_r_inv = jnp.where(lam_r > self.eps, 1.0 / lam_r, 0.0)

    new_g1 = jnp.einsum("l,lir->lir", lam_l_inv, U_t.reshape(chiL, d1, chi_new))
    new_g2 = jnp.einsum("ljr,r->ljr", Vh_t.reshape(chi_new, d2, chiR), lam_r_inv)

    p_dist = S_t ** 2  # already unit-norm (see _svd_truncate)
    mask = p_dist > 1e-20
    ee = float(-jnp.sum(jnp.where(mask, p_dist * jnp.log2(jnp.where(mask, p_dist, 1.0)), 0.0)))
    if q1 < len(self.entanglement_entropy):
        self.entanglement_entropy[q1] = ee

    self.lambdas[q2] = S_t
    self.gammas[q1] = new_g1
    self.gammas[q2] = new_g2
    self._real_chi[q2] = chi_new

    self._bond_history.append(chi_new)
    self.jsd_per_bond.append(jsd_val)

apply_ccx

apply_ccx(c1: int, c2: int, tgt: int) -> None

Toffoli via standard T-gate decomposition (all 1q/2q gates).

Source code in dense_evolution/backends/mps.py

def apply_ccx(self, c1: int, c2: int, tgt: int) -> None:
    """Toffoli via standard T-gate decomposition (all 1q/2q gates)."""
    inv2 = 1.0 / np.sqrt(2.0)
    h = inv2 * jnp.array([[1, 1], [1, -1]], dtype=complex)
    t = jnp.array([[1, 0], [0, jnp.exp(1j * jnp.pi / 4)]], dtype=complex)
    tdg = jnp.array([[1, 0], [0, jnp.exp(-1j * jnp.pi / 4)]], dtype=complex)

    self.apply_gate_1q(h, tgt)
    self.apply_cx(c2, tgt)
    self.apply_gate_1q(tdg, tgt)
    self.apply_cx(c1, tgt)
    self.apply_gate_1q(t, tgt)
    self.apply_cx(c2, tgt)
    self.apply_gate_1q(tdg, tgt)
    self.apply_cx(c1, tgt)
    self.apply_gate_1q(t, c2)
    self.apply_gate_1q(t, tgt)
    self.apply_gate_1q(h, tgt)
    self.apply_cx(c1, c2)
    self.apply_gate_1q(t, c1)
    self.apply_gate_1q(tdg, c2)
    self.apply_cx(c1, c2)

get_probabilities_sampled

get_probabilities_sampled(
    n_samples: int = 100000, seed: Optional[int] = None
) -> dict

Returns a {bitstring: empirical_probability} dict from n_samples sequential draws -- the only entry point safe for n_qubits > 24.

Source code in dense_evolution/backends/mps.py

def get_probabilities_sampled(
    self, n_samples: int = 100_000, seed: Optional[int] = None
) -> dict:
    """Returns a {bitstring: empirical_probability} dict from n_samples
    sequential draws -- the only entry point safe for n_qubits > 24."""
    from collections import Counter

    rng = np.random.default_rng(seed)
    counts: Counter = Counter()
    for _ in range(n_samples):
        bits = self._sample_bitstring(rng)
        counts["".join(map(str, bits))] += 1
    return {bitstr: c / n_samples for bitstr, c in counts.items()}

get_top_k_probable_states

get_top_k_probable_states(
    k: int = 128,
) -> Tuple[np.ndarray, np.ndarray]

Greedy beam search (beam width k) for approximately-most-probable basis states, without ever contracting to a full statevector.

Returns (indices, probabilities): indices are computational-basis integers, probabilities are exact for the states found (not approximated), sorted descending. Recall of the TRUE top states improves with k but is not guaranteed for any fixed k -- see the module docstring.

Source code in dense_evolution/backends/mps.py

def get_top_k_probable_states(self, k: int = 128) -> Tuple[np.ndarray, np.ndarray]:
    """Greedy beam search (beam width k) for approximately-most-probable
    basis states, without ever contracting to a full statevector.

    Returns (indices, probabilities): indices are computational-basis
    integers, probabilities are exact for the states found (not
    approximated), sorted descending. Recall of the TRUE top states
    improves with k but is not guaranteed for any fixed k -- see the
    module docstring."""
    paths: List[Tuple[int, jnp.ndarray]] = [(0, jnp.array([1.0 + 0.0j]))]
    for i in range(self.n):
        candidates = []
        gamma = self.gammas[i]
        lam = self.lambdas[i + 1] if (i + 1) < len(self.lambdas) else jnp.ones(gamma.shape[2])
        for idx_p, vec_p in paths:
            for bit in (0, 1):
                new_vec = jnp.einsum("l,lr->r", vec_p, gamma[:, bit, :]) * lam
                weight = float(jnp.sum(jnp.abs(new_vec) ** 2))
                candidates.append(((idx_p << 1) | bit, new_vec, weight))
        # heapq.nlargest instead of a full sort-then-slice (prog.txt
        # point 5e): only the top k by weight are ever used below, and
        # this avoids materializing/sorting the full candidates list
        # when len(candidates) >> k. Final probability order is
        # re-derived from scratch at the end of this function anyway
        # (`order = np.argsort(-probabilities)`), so which of two
        # equal-weight candidates heapq happens to prefer over sort's
        # stable order has no effect on the result.
        paths = [(idx, vec) for idx, vec, _ in heapq.nlargest(k, candidates, key=lambda c: c[2])]

    indices = np.array([p[0] for p in paths])
    amplitudes = np.array([
        complex(vec[0]) if len(vec) == 1 else complex(jnp.sum(vec))
        for _, vec in paths
    ])
    probabilities = np.abs(amplitudes) ** 2
    order = np.argsort(-probabilities)
    return indices[order], probabilities[order]

run_circuit_jit

run_circuit_jit(
    ops: List, fuse_gates: bool = False
) -> None

Runs an entire circuit through a single jax.lax.scan-fused, @jax.jit-compiled kernel instead of one eager Python call per gate -- the eager path (apply_gate_1q/apply_gate_2q/_apply_nonlocal_2q, all still available and unchanged) has zero @jax.jit anywhere and pays a host-device sync on every 2-qubit gate's bond-dimension search; measured 88.9s vs Qiskit Aer's 0.64s on a 60-qubit stress circuit -- see README changelog for the real before/after number this method produces on that same circuit.

Trade-off, explicit and intentional (not hidden): every gamma/ lambda is kept at a fixed max_bond-padded size for the rest of this instance's lifetime after this call. Structurally correct either way (zero-padding is mathematically transparent to every other method here -- contract_to_statevector, get_top_k_probable_ states, etc. all still work correctly on the padded arrays, verified), just not memory-minimal for genuinely low-entanglement circuits, which is this module's whole point for very large qubit counts. Use the eager methods directly instead when memory, not speed, is the priority -- this is an addition, not a replacement.

ops: same convention as DenseSVSimulator.run_circuit_jit_beast_mode -- list of (name, *args) tuples/lists. Unlike that method, SWAP is never decomposed into 3xCX (kept as one real gate, see _compile_mps_ops's docstring for why that matters here).

fuse_gates: opt-in, default False. When True, consecutive gates acting on the same (or a growing) qubit pair are fused into one matrix on the host before compiling (exact -- matrix multiplication, no approximation), cutting the number of scan steps and measurably faster on GPU (~2x on top of the bucketed SVD dispatch alone, see Dense-Evolution-Discovery's mps_gate_blocking_redesign_v2 experiment for the full validation, including verification against non-adjacent-gate and CCX circuits). The trade-off: self._bond_history/jsd_per_bond/ truncation_errors/entanglement_entropy get one entry per FUSED step instead of per original gate -- real diagnostics, just coarser-grained. Defaults to False so existing behavior and per-gate bookkeeping granularity are unchanged unless requested.

Source code in dense_evolution/backends/mps.py

def run_circuit_jit(self, ops: List, fuse_gates: bool = False) -> None:
    """Runs an entire circuit through a single jax.lax.scan-fused,
    @jax.jit-compiled kernel instead of one eager Python call per gate
    -- the eager path (apply_gate_1q/apply_gate_2q/_apply_nonlocal_2q,
    all still available and unchanged) has zero @jax.jit anywhere and
    pays a host-device sync on every 2-qubit gate's bond-dimension
    search; measured 88.9s vs Qiskit Aer's 0.64s on a 60-qubit stress
    circuit -- see README changelog for the real before/after number
    this method produces on that same circuit.

    Trade-off, explicit and intentional (not hidden): every gamma/
    lambda is kept at a fixed max_bond-padded size for the rest of
    this instance's lifetime after this call. Structurally correct
    either way (zero-padding is mathematically transparent to every
    other method here -- contract_to_statevector, get_top_k_probable_
    states, etc. all still work correctly on the padded arrays,
    verified), just not memory-minimal for genuinely low-entanglement
    circuits, which is this module's whole point for very large qubit
    counts. Use the eager methods directly instead when memory, not
    speed, is the priority -- this is an addition, not a replacement.

    ops: same convention as DenseSVSimulator.run_circuit_jit_beast_mode
    -- list of (name, *args) tuples/lists. Unlike that method, SWAP is
    never decomposed into 3xCX (kept as one real gate, see
    _compile_mps_ops's docstring for why that matters here).

    fuse_gates: opt-in, default False. When True, consecutive gates
    acting on the same (or a growing) qubit pair are fused into one
    matrix on the host before compiling (exact -- matrix
    multiplication, no approximation), cutting the number of scan
    steps and measurably faster on GPU (~2x on top of the bucketed
    SVD dispatch alone, see Dense-Evolution-Discovery's
    mps_gate_blocking_redesign_v2 experiment for the full validation,
    including verification against non-adjacent-gate and CCX
    circuits). The trade-off: self._bond_history/jsd_per_bond/
    truncation_errors/entanglement_entropy get one entry per FUSED
    step instead of per original gate -- real diagnostics, just
    coarser-grained. Defaults to False so existing behavior and
    per-gate bookkeeping granularity are unchanged unless requested.
    """
    dtype = self.gammas[0].dtype
    lambda_dtype = self.lambdas[0].dtype

    if fuse_gates:
        compiled_rows = _compile_mps_ops(ops, self.n)
        fused = _fuse_compiled_rows(compiled_rows, dtype) if compiled_rows else []

        if self._fused_mps_runner is None:
            self._fused_mps_runner = _build_fused_mps_runner(self.n, self.chi, self.eps, self.jsd_budget)

        gammas_padded = _pad_all_gammas(tuple(self.gammas), self.chi, dtype)
        lambdas_padded = _pad_all_lambdas(tuple(self.lambdas), self.chi, lambda_dtype)
        real_chi_initial = jnp.asarray(self._real_chi, dtype=jnp.int32)

        if fused:
            xs = _fused_entries_to_arrays(fused, dtype)
            final_gammas, final_lambdas, final_real_chi, diag = self._fused_mps_runner(
                gammas_padded, lambdas_padded, real_chi_initial, xs)

            self.gammas = [final_gammas[i] for i in range(self.n)]
            self.lambdas = [final_lambdas[i] for i in range(self.n + 1)]
            self._real_chi = np.asarray(final_real_chi)

            q1_ids = np.asarray([entry[1] for entry in fused])
            is_2q_mask = np.asarray([entry[0] == '2q' for entry in fused])
            self._record_diag_bookkeeping(diag, q1_ids, is_2q_mask)
        return

    compiled_rows = _compile_mps_ops(ops, self.n)
    ops_dtype = _real_dtype_for(dtype)

    if compiled_rows:
        ops_array = jnp.array(compiled_rows, dtype=ops_dtype)
    else:
        ops_array = jnp.zeros((0, 5), dtype=ops_dtype)

    if self._mps_runner is None:
        self._mps_runner = _build_mps_runner(self.n, self.chi, self.eps, self.jsd_budget)

    gammas_padded = jnp.stack([_pad_gamma(g, self.chi).astype(dtype) for g in self.gammas])
    lambdas_padded = jnp.stack([_pad_lambda(l, self.chi).astype(lambda_dtype) for l in self.lambdas])
    real_chi_initial = jnp.asarray(self._real_chi, dtype=jnp.int32)

    final_gammas, final_lambdas, final_real_chi, diag = self._mps_runner(
        gammas_padded, lambdas_padded, real_chi_initial, ops_array)

    self.gammas = [final_gammas[i] for i in range(self.n)]
    self.lambdas = [final_lambdas[i] for i in range(self.n + 1)]
    self._real_chi = np.asarray(final_real_chi)

    if compiled_rows:
        g_ids = np.asarray([row[0] for row in compiled_rows])
        q1_ids = np.asarray([int(row[1]) for row in compiled_rows])
        is_2q_mask = g_ids >= 20
        self._record_diag_bookkeeping(diag, q1_ids, is_2q_mask)

mps_pauli_expectation

mps_pauli_expectation(
    mps: MPSSimulator, pauli_terms
) -> complex

/ for a single Pauli string P, contracted directly against the MPS (Gamma/Lambda tensors).

pauli_terms accepts the same three forms as physics.observables.pauli_expectation (a string, e.g. 'XIZ'; a dict {qubit: 'X'|'Y'|'Z'}; or an iterable of (qubit, pauli) pairs) -- reuses that module's own _normalize_terms so both functions agree on parsing by construction, not by parallel reimplementation. See _mps_transfer_sweep for why the division is needed.

Source code in dense_evolution/backends/mps.py

def mps_pauli_expectation(mps: "MPSSimulator", pauli_terms) -> complex:
    """<psi|P|psi> / <psi|psi> for a single Pauli string P, contracted
    directly against the MPS (Gamma/Lambda tensors).

    pauli_terms accepts the same three forms as
    physics.observables.pauli_expectation (a string, e.g. 'XIZ'; a dict
    {qubit: 'X'|'Y'|'Z'}; or an iterable of (qubit, pauli) pairs) -- reuses
    that module's own `_normalize_terms` so both functions agree on
    parsing by construction, not by parallel reimplementation. See
    _mps_transfer_sweep for why the division is needed.
    """
    assignment = _normalize_terms(pauli_terms, mps.n)
    raw, norm_sq = _mps_transfer_sweep(mps, assignment, need_norm=True)
    return raw / norm_sq

mps_pauli_sum_expectation

mps_pauli_sum_expectation(
    mps: MPSSimulator, terms
) -> complex

sum_i coeff_i * / -- same terms format as physics.observables.pauli_sum_expectation: an iterable of (coeff, pauli_terms) pairs. does not depend on which Pauli string is being measured, so it is computed once via its own sweep and applied to the whole sum, instead of once per term.

Source code in dense_evolution/backends/mps.py

def mps_pauli_sum_expectation(mps: "MPSSimulator", terms) -> complex:
    """sum_i coeff_i * <psi|P_i|psi> / <psi|psi> -- same terms format as
    physics.observables.pauli_sum_expectation: an iterable of
    (coeff, pauli_terms) pairs. <psi|psi> does not depend on which Pauli
    string is being measured, so it is computed once via its own sweep
    and applied to the whole sum, instead of once per term."""
    terms = list(terms)
    if not terms:
        return 0j
    _, norm_sq = _mps_transfer_sweep(mps, {}, need_norm=True)
    raw_sum = sum(
        coeff * _mps_transfer_sweep(mps, _normalize_terms(pauli_terms, mps.n), need_norm=False)[0]
        for coeff, pauli_terms in terms
    )
    return raw_sum / norm_sq

bond_convergence

bond_convergence(
    ops: List,
    n_qubits: int,
    observables: list,
    bonds: List[int],
    tol: float = 0.001,
    **mps_kwargs,
) -> BondConvergenceResult

Runs the same circuit at every value in bonds (increasing) and checks whether the reported observables have actually converged with respect to bond dimension, instead of trusting a single run's own internal diagnostics.

Requires len(bonds) >= 3. Two bonds give exactly one discrepancy, which is a single number with no way to tell whether it is still shrinking toward tol or has already stalled -- measured on a 40-qubit, 4-layer brickwall circuit, chi=4->8->32 gave || discrepancies of ~4.7e-2 then ~1.2e-2 (chi_used never hit its own cap, so this is a real not_converged, not an artifact of running out of bond dimension): a two-bond check (chi=4 vs 8) would see only the first number and have no basis to call it anything, while three bonds show a trend that is decreasing but still two orders of magnitude above any reasonable tol.

A verdict of "converged" additionally requires the successive discrepancies to be monotonically non-increasing, not just that the last one is below tol -- a single small discrepancy proves nothing about the trend on its own, which is the same failure mode as the two-bond case above, one level up. (Ties count as non-increasing: an exactly-converged observable, e.g. a GHZ chain whose bond dimension never needs to grow, produces identical values -- and therefore zero discrepancies -- at every bond, which must count as converged.)

avg_jsd and budget_violations (from the underlying MPSSimulator runs) are reported per bond for context only, never used to decide the verdict -- a low average JSD is computed per truncation step and says nothing about whether the specific observable being tracked has settled down as max_bond grows.

If max_bond_used() at the highest bond still equals that bond's cap, the truncation never had headroom below max_bond at any cut, so no tolerance can be certified from this data: every observable's verdict becomes "undecidable" regardless of its own discrepancies.

Source code in dense_evolution/backends/mps.py

def bond_convergence(
    ops: List, n_qubits: int, observables: list, bonds: List[int],
    tol: float = 1e-3, **mps_kwargs,
) -> BondConvergenceResult:
    """Runs the same circuit at every value in `bonds` (increasing) and
    checks whether the reported observables have actually converged with
    respect to bond dimension, instead of trusting a single run's own
    internal diagnostics.

    Requires len(bonds) >= 3. Two bonds give exactly one discrepancy,
    which is a single number with no way to tell whether it is still
    shrinking toward `tol` or has already stalled -- measured on a
    40-qubit, 4-layer brickwall circuit, chi=4->8->32 gave |<Z0>|
    discrepancies of ~4.7e-2 then ~1.2e-2 (chi_used never hit its own
    cap, so this is a real not_converged, not an artifact of running out
    of bond dimension): a two-bond check (chi=4 vs 8) would see only the
    first number and have no basis to call it anything, while three
    bonds show a trend that is decreasing but still two orders of
    magnitude above any reasonable `tol`.

    A verdict of "converged" additionally requires the successive
    discrepancies to be monotonically non-increasing, not just that the
    last one is below `tol` -- a single small discrepancy proves nothing
    about the trend on its own, which is the same failure mode as the
    two-bond case above, one level up. (Ties count as non-increasing: an
    exactly-converged observable, e.g. a GHZ chain whose bond dimension
    never needs to grow, produces identical values -- and therefore
    zero discrepancies -- at every bond, which must count as converged.)

    avg_jsd and budget_violations (from the underlying MPSSimulator runs)
    are reported per bond for context only, never used to decide the
    verdict -- a low average JSD is computed per truncation step and says
    nothing about whether the specific observable being tracked has
    settled down as `max_bond` grows.

    If max_bond_used() at the highest bond still equals that bond's cap,
    the truncation never had headroom below max_bond at any cut, so no
    tolerance can be certified from this data: every observable's verdict
    becomes "undecidable" regardless of its own discrepancies.
    """
    if len(bonds) < 3:
        raise ValueError(f"bond_convergence needs at least 3 bonds to detect a trend, got {len(bonds)}")

    chi_used, avg_jsd, budget_violations = [], [], []
    values = [[] for _ in observables]
    for bond in bonds:
        mps = MPSSimulator(n_qubits=n_qubits, max_bond=bond, **mps_kwargs)
        mps.run_circuit_jit(ops)
        chi_used.append(mps.max_bond_used())
        avg_jsd.append(mps.avg_jsd())
        budget_violations.append(mps.budget_violations)
        for obs_idx, obs in enumerate(observables):
            values[obs_idx].append(mps_pauli_expectation(mps, obs))

    diffs = [
        [abs(vals[i + 1] - vals[i]) for i in range(len(vals) - 1)]
        for vals in values
    ]

    undecidable = chi_used[-1] >= bonds[-1]
    verdicts = []
    for d in diffs:
        if undecidable:
            verdicts.append("undecidable")
        elif all(d[i + 1] <= d[i] for i in range(len(d) - 1)) and d[-1] < tol:
            verdicts.append("converged")
        else:
            verdicts.append("not_converged")

    return BondConvergenceResult(
        bonds=list(bonds), chi_used=chi_used, avg_jsd=avg_jsd,
        budget_violations=budget_violations, values=values, diffs=diffs,
        verdicts=verdicts,
    )

See also: DenseSVSimulator for exact statevector simulation when entanglement is too high for a bounded bond dimension, and Chunk for anti-OOM dense simulation at large qubit counts without bond-dimension truncation.