mps-xtrap 2.4.1

MPS-based quantum circuit simulator with multilevel Richardson extrapolation, gate fusion, and qubit measurement

mps_xtrap — MPS Quantum Circuit Simulator with Richardson Extrapolation

A production-ready quantum circuit simulator based on Matrix Product States (MPS), featuring multilevel Richardson extrapolation to push expectation value accuracy toward the infinite bond-dimension limit, full GPU acceleration via CuPy, gate fusion to reduce SVD calls, and qubit measurement via projector operators.

222 tests passing across 15 sections — production ready.

---

Licensing

mps_xtrap is released under the Business Source License 1.1 (BUSL-1.1).

• Non-commercial use — free to use under the terms of the BUSL-1.1.
• Commercial use — requires a separate commercial license agreement. Please contact [e-mail] to arrange one.
• Change Date — on 2029-01-01 the Licensed Work converts to the MIT License.

Sponsorships and donations are welcomed at: https://flutterwave.com/pay/6yptuvqab8cu

---

Installation

pip install numpy            # required
pip install cupy-cuda12x     # optional: GPU support (match your CUDA version)

---

Quick Start

from mps_xtrap import Circuit, MPSSimulator, MultiChiRunner, SweepConfig

# --- Simple simulation ---
c = Circuit(4)
c.h(0).cx(0, 1).cx(1, 2).cx(2, 3)   # GHZ state

sim = MPSSimulator(chi=64)
state = sim.run(c)
print(state.expectation_pauli_z(0))   # -> 0.0

# --- GPU simulation ---
sim_gpu = MPSSimulator(chi=64, device='cuda')
state_gpu = sim_gpu.run(c)

# --- With gate fusion ---
sim_fused = MPSSimulator(chi=64, fusion=True)
state = sim_fused.run(c)

# --- Extrapolated simulation (bond dims auto-generated: [16, 32, 64]) ---
runner = MultiChiRunner(SweepConfig(base_chi=16, scaling_factor=2, n_levels=3))
result = runner.run(c, observables={'Z0': ('Z', 0), 'Z1': ('Z', 1)})
print(result.summary())

# --- Correlated-system expectation values (not single-site!) ---
# Two-site correlator <Z0 Z3>, sites need not be adjacent:
print(sim.expectation_value(state, 'ZZ', site=0, site2=3))
# Equivalent, and generalizes to any number of sites:
print(sim.expectation_correlator(state, 'ZZ', sites=[0, 3]))
# N-site correlator <Z0 X2 Y5>:
print(sim.expectation_correlator(state, 'ZXY', sites=[0, 2, 5]))
# Extrapolated correlators work the same way through MultiChiRunner:
result = runner.run(c, observables={'ZZ_03': ('ZZ', 0, 3), 'ZXY_025': ('ZXY', [0, 2, 5])})
print(result.summary())

# --- Qubit measurement (new in 2.2.0) ---
from mps_xtrap import measure_qubit, measure_qubits, sample_counts, full_distribution

mres = measure_qubit(state, site=0)
print(mres.summary())              # P(0)=..., P(1)=...

jres = measure_qubits(state, sites=[0, 1, 2])
print(jres.summary())              # joint probability table

sc = sample_counts(state, sites=[0, 1], shots=1024, seed=42)
print(sc.summary())                # shot histogram

dist = full_distribution(state)    # {bitstring: probability}, n <= 20 only

# --- Extrapolate measurement probabilities (new in 2.2.0) ---
runner = MultiChiRunner(SweepConfig(base_chi=16, scaling_factor=2, n_levels=3))
result = runner.run(
    c,
    observables={'Z0': ('Z', 0)},
    measure={'P0_q0': (0, 'prob0'), 'P1_q0': (0, 'prob1')},
)
print(result.summary())

---

Architecture

mps_xtrap/
├── core/
│   └── mps.py              # MPS tensor train, SVD truncation, GPU/CPU backend
├── gates/
│   └── __init__.py         # Full gate library (H, CNOT, Rx, ZZ, ...)
├── circuits/
│   └── __init__.py         # Circuit builder API + MPSSimulator engine
├── extrapolation/
│   └── __init__.py         # Richardson extrapolation engine + MultiChiRunner
├── fusion/
│   └── __init__.py         # Gate fusion (GateFuser, fuse)      <- added 2.1.0
├── measurement/
│   └── __init__.py         # Projector measurement + sampling    <- added 2.2.0
├── tests/
│   ├── test_all.py         # Core unit tests (40 tests)
│   ├── test_fusion.py      # Fusion unit tests
│   ├── test_measurement.py # Measurement unit tests (59 tests)   <- added 2.2.0
│   └── test_production.py  # Production suite (123 tests, 14 sections)
├── examples/
│   └── examples.py         # Runnable examples
└── cli.py                  # Command-line interface

---

Qubit Measurement

Theory

A projective measurement on qubit s for outcome k in {0, 1} uses the projector:

P_k = |k><k|

P_0 = [[1, 0],   (projects onto |0>)
       [0, 0]]

P_1 = [[0, 0],   (projects onto |1>)
       [0, 1]]

The Born-rule probability of outcome k at site s is:

prob(k, s) = <psi| P_k^(s) |psi>

This is a single-site expectation value of the projector — computed efficiently using the MPS canonicalization machinery already in the simulator. Post-measurement collapse applies P_k as a gate to a copy of the MPS and renormalizes the tensor at site s by 1/sqrt(prob(k,s)).

Single-qubit measurement

from mps_xtrap import measure_qubit

res = measure_qubit(state, site=2)
print(res.prob0, res.prob1)      # Born-rule probabilities (always sum to 1)

# With collapse: returns the normalized post-measurement MPS for each outcome
res = measure_qubit(state, site=2, collapse=True)
state_if_0 = res.post_state_0   # MPS conditioned on observing |0>
state_if_1 = res.post_state_1   # MPS conditioned on observing |1> (None if prob=0)

Multi-qubit joint measurement

from mps_xtrap import measure_qubits

# Joint probability via chain rule: P(b0,b1,...) = P(b0)*P(b1|b0)*...
res = measure_qubits(state, sites=[0, 1, 2])
print(res.probabilities)            # {(0,0,0): p, (0,0,1): p, ...}
print(res.marginals[0].prob0)       # marginal P(0) at site 0
print(res.most_likely_bitstring())  # e.g. (0, 0, 1)
print(res.summary())                # formatted table

Simulated shot sampling

from mps_xtrap import sample_counts

sc = sample_counts(state, sites=[0, 1, 2, 3], shots=1024, seed=42)
print(sc.counts)            # {'0000': 512, '1111': 512, ...}
print(sc.probabilities())   # normalize to [0, 1]
print(sc.summary())

Full probability distribution (small systems)

from mps_xtrap import full_distribution

dist = full_distribution(state)   # {'0000': 0.5, '1111': 0.5, ...}
# Only feasible for n <= 20 qubits

Projector constants

from mps_xtrap import projector, P0, P1

P0 = projector(0)   # [[1,0],[0,0]]
P1 = projector(1)   # [[0,0],[0,1]]

# P0 and P1 can be passed directly to MPS.expectation_single()
prob0 = float(state.expectation_single(P0, site=2).real)

---

Extrapolating Measurement Probabilities

Since probabilities are projector expectation values, they obey the same chi^{-alpha} power law as Pauli observables and can be Richardson-extrapolated.

MultiChiRunner.run() accepts a measure= keyword argument:

runner = MultiChiRunner(SweepConfig(base_chi=16, scaling_factor=2, n_levels=3))
result = runner.run(
    circuit,
    observables={'Z0': ('Z', 0)},
    measure={
        'P0_q1': (1, 'prob0'),   # extrapolate P(outcome=0) at site 1
        'P1_q1': (1, 'prob1'),   # extrapolate P(outcome=1) at site 1
    },
)
print(result.summary())
p0 = result.observables['P0_q1'].extrapolated
p0_raw = result.observables['P0_q1'].raw_values   # one value per chi level

Both observables and measure are optional keyword arguments and can be used together or independently.

---

Gate Fusion

Gate fusion merges consecutive gates on the same qubit(s) into a single combined gate. Accuracy is unchanged; only SVD call count is reduced.

from mps_xtrap import MPSSimulator, GateFuser
from mps_xtrap.fusion import fuse

# At construction
sim = MPSSimulator(chi=64, fusion=True)

# Fine-grained control
fuser = GateFuser(fuse_single=True, fuse_two=True, max_window=4)
sim = MPSSimulator(chi=64, fusion=fuser)

# Replace after construction
sim.fusion = GateFuser(fuse_single=False, fuse_two=True)

# Fuse a circuit directly
fc = fuse(c)
report = GateFuser().fusion_report(c)
# {'original_count': 20, 'fused_count': 11, 'saved': 9, ...}

---

GPU Support

sim = MPSSimulator(chi=128, device='cuda')
state = sim.run(circuit)

state_gpu = state.to('cuda')
state_cpu = state_gpu.to('cpu')

Automatic CPU fallback if CuPy is missing or no GPU is detected. Most beneficial at large bond dimensions (chi >= 64).

---

Richardson Extrapolation

MPS truncation error follows a power law in chi:

<O>(chi) ~ <O>(inf) + a1/chi^alpha + a2/chi^(2*alpha) + ...

Richardson extrapolation cancels successive error orders. The same applies to Born-rule probabilities, since they are projector expectation values.

from mps_xtrap import MultiChiRunner, SweepConfig

runner = MultiChiRunner(
    SweepConfig(base_chi=16, scaling_factor=2, n_levels=3, alpha=None),
    device='cpu',
    verbose=True,
)
result = runner.run(
    circuit,
    observables={'Z0': ('Z', 0)},
    measure={'P0_q0': (0, 'prob0'), 'P1_q0': (0, 'prob1')},
)
r = result.observables['Z0']
print(r.extrapolated, r.uncertainty, r.is_reliable, r.alpha)

Reliability diagnostics are built in: non-monotone convergence, corrections not decreasing, large uncertainty vs signal, and implausible alpha values are all flagged automatically.

---

Gate Library

Single-qubit: I, X, Y, Z, H, S, T, Sdg, Tdg, Rx(theta), Ry(theta), Rz(theta), P(phi), U(theta,phi,lambda)

Two-qubit: CNOT/CX, CZ, SWAP, iSWAP, XX(theta), YY(theta), ZZ(theta), CRz(theta), CP(phi)

---

Circuit Builder API

import numpy as np
from mps_xtrap import Circuit, MPSSimulator

c = Circuit(6)
c.h(0).cx(0, 1).rz(np.pi/4, 2).zz(0.5, 3, 4).swap(4, 5)

state = MPSSimulator(chi=64).run(c)
print(state.expectation_pauli_z(0))
print(state.bond_dimensions())
print(state.total_truncation_error())

---

Running the Tests

python -m unittest mps_xtrap.tests.test_all          # 40 core tests
python -m unittest mps_xtrap.tests.test_fusion       # fusion tests
python -m unittest mps_xtrap.tests.test_measurement  # 59 measurement tests
python -m unittest mps_xtrap.tests.test_production   # 123 production tests

---

API Reference

MPSSimulator(chi, svd_threshold=1e-14, device='cpu', fusion=False)

• .run(circuit) -> MPS
• .run_from(circuit, state) -> MPS
• .expectation_value(state, observable, site, site2=None) -> float Single-site by default (e.g. 'Z'). Pass site2 for a two-site correlator (e.g. observable='ZZ', site, site2); the two sites need not be adjacent.
• .expectation_correlator(state, observable, sites) -> float N-site Pauli-string correlator for any number of sites, adjacent or not, e.g. expectation_correlator(state, 'ZXZ', [0, 2, 5]). len(observable) must match len(sites).
• .fusion — read/replace the active fuser

GateFuser(fuse_single=True, fuse_two=True, max_window=None)

• __call__(circuit_or_instructions) -> fused Circuit or list
• .fusion_report(circuit) -> dict

fuse(circuit, ...) -> Circuit

MPS

• .expectation_pauli_z/x/y(site) -> float — single-site only
• .expectation_single(op, site) -> complex — single-site only, arbitrary 2x2 operator
• .expectation_two(op, site_i) -> complex — two-site operator on the adjacent pair (site_i, site_i + 1)
• .expectation_two_site(op_i, op_j, site_i, site_j) -> complex — two-site correlator <op_i(site_i) op_j(site_j)> for arbitrary, possibly non-adjacent sites
• .expectation_multi_site(ops, sites) -> complex — the general, correlated-system building block: <op_0(sites[0]) op_1(sites[1]) ...> for any number of single-site operators at arbitrary, possibly non-adjacent sites, computed in a single O(n · chi^3) sweep. expectation_two_site is a thin wrapper around this.
• .to_statevector() -> ndarray (n <= 20)
• .bond_dimensions() -> list, .total_truncation_error() -> float
• .to(device) -> MPS, .copy() -> MPS

Single-site vs. correlated systems: expectation_pauli_z/x/y and expectation_single only ever act on one site and say nothing about correlations between qubits. To get expectation values of the correlated system — e.g. <Z0 Z3>, <X1 Y4 Z7> — use expectation_two_site / expectation_multi_site on MPS directly, or the higher-level MPSSimulator.expectation_value(..., site2=...) / MPSSimulator.expectation_correlator(...) wrappers above. These are also what MultiChiRunner.run(..., observables=...) uses under the hood — see below.

MultiChiRunner(config, device='cpu', verbose=True)

• .run(circuit, observables=None, measure=None) -> MultiObservableResult
  • observables: {label: spec}, where spec is one of:
    • (obs_type, site) — single-site Pauli expectation value, e.g. {'Z0': ('Z', 0)}
    • (obs_type, site, site2) — two-site correlator of the correlated system, sites need not be adjacent, e.g. {'ZZ_03': ('ZZ', 0, 3)}
    • (obs_type, [site0, site1, ...]) — N-site correlator over any number of (possibly non-adjacent) sites, e.g. {'ZXZ_025': ('ZXZ', [0, 2, 5])}, where len(obs_type) must match the number of sites
  • measure: {label: (site, 'prob0'|'prob1')} — Born-rule probabilities

SweepConfig(base_chi=64, scaling_factor=2.0, n_levels=3, alpha=None)

• .bond_dims — auto-generated chi list
• .effective_chi()

ExtrapolationResult

• .extrapolated, .raw_values, .uncertainty, .alpha, .scaling_factor
• .is_reliable, .reliability_notes
• .summary(), .improvement_factor()

MultiObservableResult

• .observables — {name: ExtrapolationResult}
• .bond_dims, .summary()

projector(k) -> ndarray

Returns |k><k| for k in {0, 1}.

P0, P1

Module-level projector constants.

MeasurementEngine(state)

• .measure_qubit(site, collapse=False) -> SingleMeasurementResult
• .measure_qubits(sites) -> MultiQubitMeasurementResult
• .joint_probability(sites, outcomes) -> float
• .sample_counts(sites, shots, seed) -> SampledCounts
• .full_distribution() -> dict

measure_qubit(state, site, collapse=False) -> SingleMeasurementResult

• .prob0, .prob1, .post_state_0, .post_state_1
• .most_likely() -> int, .summary() -> str

measure_qubits(state, sites) -> MultiQubitMeasurementResult

• .probabilities, .marginals, .most_likely_bitstring(), .summary()

sample_counts(state, sites, shots=1024, seed=None) -> SampledCounts

• .counts, .probabilities(), .shots, .summary()

full_distribution(state) -> dict

Full {bitstring: probability} for n <= 20.

show_license() -> None

Prints the mps_xtrap licensing notice.

---

CLI

mps_xtrap simulate    --circuit ghz   --n 10 --chi 64
mps_xtrap extrapolate --circuit ising --n 8  --chi-base 16 --levels 3
mps_xtrap benchmark   --circuit ising --n 8  --chi-start 8 --chi-levels 4
mps_xtrap info

---

Limitations

• Non-adjacent two-qubit gates use SWAP chains (increases depth)
• to_statevector() and full_distribution() only for n <= 20
• measure_qubits() enumerates 2^m bitstrings — only for m <= 20; use sample_counts for large m
• Richardson extrapolation assumes power-law error decay; may not hold for highly entangled circuits
• GPU requires NVIDIA hardware and matching CuPy

---

Dependencies

• Python 3.8+
• NumPy >= 1.21
• CuPy (optional, for GPU)

---