PTSBE end-to-end workflow

PTSBE (Pre-Trajectory Sampling with Batch Execution) is a method for sampling from noisy quantum circuits efficiently. Instead of simulating the full density matrix and sampling once per shot, PTSBE:

  1. Traces the kernel to get the gate sequence and qubit layout.

  2. Extracts noise sites by matching the noise model to the trace (each noisy gate becomes a noise site with a set of Kraus outcomes, e.g. \(I\), \(X\), \(Y\), \(Z\) for depolarization).

  3. Generates trajectories — each trajectory is one possible realization of noise (one Kraus outcome per site). A sampling strategy decides which trajectories to use (e.g. by sampling trajectories proporitional to their likelihood).

  4. Allocates shots across trajectories (e.g. by error likelihood).

  5. Runs batches — for each trajectory, the circuit is run as a noiseless circuit with that trajectory’s outcomes applied; results are collected.

  6. Aggregates all per-trajectory counts into a single SampleResult.

Given sufficient trajectory and shot samples PTSBE will be equivalent to standard trajectory based sampling up to sampling noise. However, it has the additional advantage of being able to batch many shots per trajectory and to control simulation cost via the number of trajectories. This notebook runs the full workflow with a single API call: cudaq.ptsbe.sample().

Set up the environment

[1]:

import sys, os
_root = os.path.join(os.path.expanduser("~/Devel/cudaq-pstbe/vendor/cuda-quantum"), "build")
sys.path.insert(0, os.path.join(_root, "python"))
_build_lib = os.path.join(_root, "lib")
os.environ["LD_LIBRARY_PATH"] = _build_lib + ":" + os.environ.get("LD_LIBRARY_PATH", "")

[2]:

import cudaq

cudaq.set_target("qpp-cpu")
cudaq.set_random_seed(42)

Define the circuit and noise model

Define a kernel and attach a noise model. Each gate you add to the noise model becomes a noise site when that gate appears in the circuit. In this example, for single-qubit gates (e.g. \(H\)) we use a DepolarizationChannel with one qubit and for \(CNOT\) we add the qubit pair [control, target] with Depolarization2 (two-qubit depolarization channel). Here we use a small Bell-style circuit to demonstrate the

[3]:

depol_1q = 0.05
depol_2q = 0.03

@cudaq.kernel
def bell_with_noise():
    q = cudaq.qvector(2)
    h(q[0])
    x.ctrl(q[0], q[1])
    mz(q)

noise = cudaq.NoiseModel()
noise.add_all_qubit_channel("h", cudaq.DepolarizationChannel(depol_1q))
noise.add_all_qubit_channel("x", cudaq.Depolarization2(depol_2q), num_controls=1)

Inline noise with apply_noise

Instead of (or in addition to) attaching noise via a NoiseModel, you can place noise at specific points in the kernel with cudaq.apply_noise. PTSBE traces these calls and includes them as noise sites alongside any model-attached noise.

[4]:

@cudaq.kernel
def bell_inline_noise():
    q = cudaq.qvector(2)
    h(q[0])
    cudaq.apply_noise(cudaq.DepolarizationChannel, depol_1q, q[0])
    x.ctrl(q[0], q[1])
    cudaq.apply_noise(cudaq.Depolarization2, depol_2q, q[0], q[1])
    mz(q)

result_inline = cudaq.ptsbe.sample(
    bell_inline_noise,
    shots_count=10000,
)
print("Inline apply_noise result:")
print(result_inline)

Inline apply_noise result:
{ 00:4873 01:74 10:74 11:4979 }

Run PTSBE sampling

Call cudaq.ptsbe.sample() with the kernel, noise model, and shot count. This example uses OrderedSamplingStrategy plus LOW_WEIGHT_BIAS shot allocation to emphasize likely low-error trajectories with a smaller shot budget.

Optional arguments are:

  • sampling_strategy — how trajectories are chosen:

    • ProbabilisticSamplingStrategy(seed=...) (default): Performs Monte Carlo sampling over trajectories.

    • ExhaustiveSamplingStrategy(): use all possible trajectories (every combination of Kraus outcomes per noise site).

    • OrderedSamplingStrategy(): Select the top-\(k\) trajectories by probability (highest first), up to max_trajectories.

  • shot_allocation — how shots are split across the chosen trajectories:

    • PROPORTIONAL (default): allocate shots in proportion to each sampled trajectory’s weighting.

    • UNIFORM: give each trajectory the same number of shots.

    • LOW_WEIGHT_BIAS: bias more shots toward low-weight (fewer errors) trajectories, optional bias_strength (default 2.0).

    • HIGH_WEIGHT_BIAS: bias more shots toward high-weight trajectories, optional bias_strength (default 2.0).

  • max_trajectories: cap the number of trajectories (useful for large shot counts).

  • return_execution_data: If True, the result includes trace instructions and per-trajectory data (result.ptsbe_execution_data). This is currently an experimental API and subject to change in future releases. See section 5 below.

[5]:

shots = 10000
max_traj = 16
strategy = cudaq.ptsbe.OrderedSamplingStrategy()

def print_result(label, sample_result):
    print(f"{label:<36} shots={sample_result.get_total_shots():>6}  counts={sample_result}")

biased_alloc = cudaq.ptsbe.ShotAllocationStrategy(
    cudaq.ptsbe.ShotAllocationType.LOW_WEIGHT_BIAS,
    bias_strength=3.0,
    seed=42,
)
result_biased = cudaq.ptsbe.sample(
    bell_with_noise,
    noise_model=noise,
    shots_count=shots,
    sampling_strategy=strategy,
    shot_allocation=biased_alloc,
    max_trajectories=max_traj,
)

unbiased_alloc = cudaq.ptsbe.ShotAllocationStrategy(
    cudaq.ptsbe.ShotAllocationType.PROPORTIONAL,
    seed=42,
)
result_unbiased = cudaq.ptsbe.sample(
    bell_with_noise,
    noise_model=noise,
    shots_count=shots,
    sampling_strategy=strategy,
    shot_allocation=unbiased_alloc,
    max_trajectories=max_traj,
)

cudaq.set_target("density-matrix-cpu")
result_standard = cudaq.sample(
    bell_with_noise,
    noise_model=noise,
    shots_count=shots,
)
cudaq.set_target("qpp-cpu")

print(f"Comparison with same shots={shots}, max_trajectories={max_traj}")
print_result("PTSBE (ordered + low-weight bias)", result_biased)
print_result("PTSBE (ordered + proportional)", result_unbiased)
print_result("Density-matrix noisy sample", result_standard)

Comparison with same shots=10000, max_trajectories=16
PTSBE (ordered + low-weight bias)    shots= 10000  counts={ 00:5045 01:6 10:4 11:4945 }

PTSBE (ordered + proportional)       shots= 10000  counts={ 00:4877 01:49 10:55 11:5019 }

Density-matrix noisy sample          shots= 10000  counts={ 00:4943 01:83 10:98 11:4876 }

Larger circuit for execution data

To demonstrate PTSBE’s execution data on a larger circuit, we define a 12-qubit GHZ state with depolarization on every gate. This creates 12 noise sites with a large trajectory space.

[6]:

n_qubits = 12

@cudaq.kernel
def ghz(n: int):
    q = cudaq.qvector(n)
    h(q[0])
    for i in range(1, n):
        x.ctrl(q[i - 1], q[i])
    mz(q)

ghz_depol_1q = 0.01
ghz_depol_2q = 0.01

ghz_noise = cudaq.NoiseModel()
ghz_noise.add_all_qubit_channel("h", cudaq.DepolarizationChannel(ghz_depol_1q))
ghz_noise.add_all_qubit_channel("x", cudaq.Depolarization2(ghz_depol_2q), num_controls=1)

Inspecting trajectories with execution data

Pass return_execution_data=True to get the PTSBE execution data via result.ptsbe_execution_data. This reveals why PTSBE is efficient: most probability mass concentrates on a few low-error trajectories, so a small number of circuit simulations captures the bulk of the physics. The highest-probability trajectory (no errors at any noise site) typically receives the vast majority of shots, while multi-error trajectories are exponentially suppressed.

Note: this is an experimental API and may change in future releases.

[7]:

from collections import Counter

exec_shots = 1_000_000
max_traj = 256
result_with_data = cudaq.ptsbe.sample(
    ghz,
    n_qubits,
    noise_model=ghz_noise,
    shots_count=exec_shots,
    max_trajectories=max_traj,
    return_execution_data=True,
)
assert result_with_data.has_execution_data()
data = result_with_data.ptsbe_execution_data

trajs = sorted(data.trajectories, key=lambda t: t.probability, reverse=True)
print(f"Trajectories: {len(trajs)}, Total shots: {exec_shots:,}\n")

# `t.probability` is the per-shot probability mass in the full trajectory distribution.
print("Top 5 trajectories (highest probability):")
cumulative = 0
for rank, t in enumerate(trajs[:5], 1):
    cumulative += t.num_shots
    n_errors = sum(1 for s in t.kraus_selections if s.is_error)
    p_theory = t.probability
    p_empirical = t.num_shots / exec_shots
    print(f"  #{rank}: p_theory={p_theory:.6f}, p_empirical={p_empirical:.6f}, "
          f"shots={t.num_shots:,}, errors={n_errors}, "
          f"cumulative shots={100 * cumulative / exec_shots:.1f}%")

# Lowest-probability trajectory
lowest = trajs[-1]
n_errors_low = sum(1 for s in lowest.kraus_selections if s.is_error)
print(f"\nLowest-probability trajectory:")
print(f"  prob={lowest.probability:.2e}, shots={lowest.num_shots}, errors={n_errors_low}")

noise_instructions = {i: inst for i, inst in enumerate(data.instructions)
                      if inst.type == cudaq.ptsbe.TraceInstructionType.Noise}

def fmt_selection(sel):
    channel = noise_instructions[sel.circuit_location]
    label = "error" if sel.is_error else "no-error"
    return (f"    site {sel.circuit_location} [{channel.name} on q{channel.targets}]: "
            f"K{sel.kraus_operator_index} ({label})")

highest = trajs[0]
print(f"\nHighest-probability trajectory (id={highest.trajectory_id}):")
print(f"  prob={highest.probability:.6f}, shots={highest.num_shots:,}")
print(f"  Kraus selections:")
for sel in highest.kraus_selections:
    print(fmt_selection(sel))

print(f"\nLowest-probability trajectory (id={lowest.trajectory_id}):")
print(f"  prob={lowest.probability:.2e}, shots={lowest.num_shots}")
print(f"  Kraus selections:")
for sel in lowest.kraus_selections:
    print(fmt_selection(sel))

# Error count histogram
error_counts = Counter(
    sum(1 for s in t.kraus_selections if s.is_error) for t in trajs
)
print("\nTrajectories grouped by error count:")
for n_err in sorted(error_counts):
    n_traj = error_counts[n_err]
    total = sum(t.num_shots for t in trajs
                if sum(1 for s in t.kraus_selections if s.is_error) == n_err)
    print(f"  {n_err} errors: {n_traj} trajectories, "
          f"{total:,} shots ({100 * total / exec_shots:.1f}%)")

Trajectories: 147, Total shots: 1,000,000
[2026-03-12 00:13:55.316] [warning] [PTSBESampleResult.cpp:20] PTSBE execution data API is experimental and may change in a future release.

Top 5 trajectories (highest probability):
  #1: p_theory=0.886385, p_empirical=0.888297, shots=888,297, errors=0, cumulative shots=88.8%
  #2: p_theory=0.002984, p_empirical=0.002378, shots=2,378, errors=1, cumulative shots=89.1%
  #3: p_theory=0.002984, p_empirical=0.002333, shots=2,333, errors=1, cumulative shots=89.3%
  #4: p_theory=0.002984, p_empirical=0.004269, shots=4,269, errors=1, cumulative shots=89.7%
  #5: p_theory=0.000597, p_empirical=0.000434, shots=434, errors=1, cumulative shots=89.8%

Lowest-probability trajectory:
  prob=2.71e-10, shots=359, errors=3

Highest-probability trajectory (id=0):
  prob=0.886385, shots=888,297
  Kraus selections:
    site 1 [depolarization_channel on q[0]]: K0 (no-error)
    site 3 [depolarization2 on q[1, 0]]: K0 (no-error)
    site 5 [depolarization2 on q[2, 1]]: K0 (no-error)
    site 7 [depolarization2 on q[3, 2]]: K0 (no-error)
    site 9 [depolarization2 on q[4, 3]]: K0 (no-error)
    site 11 [depolarization2 on q[5, 4]]: K0 (no-error)
    site 13 [depolarization2 on q[6, 5]]: K0 (no-error)
    site 15 [depolarization2 on q[7, 6]]: K0 (no-error)
    site 17 [depolarization2 on q[8, 7]]: K0 (no-error)
    site 19 [depolarization2 on q[9, 8]]: K0 (no-error)
    site 21 [depolarization2 on q[10, 9]]: K0 (no-error)
    site 23 [depolarization2 on q[11, 10]]: K0 (no-error)

Lowest-probability trajectory (id=108):
  prob=2.71e-10, shots=359
  Kraus selections:
    site 1 [depolarization_channel on q[0]]: K0 (no-error)
    site 3 [depolarization2 on q[1, 0]]: K13 (error)
    site 5 [depolarization2 on q[2, 1]]: K12 (error)
    site 7 [depolarization2 on q[3, 2]]: K0 (no-error)
    site 9 [depolarization2 on q[4, 3]]: K0 (no-error)
    site 11 [depolarization2 on q[5, 4]]: K7 (error)
    site 13 [depolarization2 on q[6, 5]]: K0 (no-error)
    site 15 [depolarization2 on q[7, 6]]: K0 (no-error)
    site 17 [depolarization2 on q[8, 7]]: K0 (no-error)
    site 19 [depolarization2 on q[9, 8]]: K0 (no-error)
    site 21 [depolarization2 on q[10, 9]]: K0 (no-error)
    site 23 [depolarization2 on q[11, 10]]: K0 (no-error)

Trajectories grouped by error count:
  0 errors: 1 trajectories, 888,297 shots (88.8%)
  1 errors: 126 trajectories, 103,948 shots (10.4%)
  2 errors: 19 trajectories, 7,396 shots (0.7%)
  3 errors: 1 trajectories, 359 shots (0.0%)