Experiments and Noise Modeling

These examples walk through several of the most common numerical error-correction experiments with the CUDA-Q QEC library: modeling noise at the code-capacity and circuit-level, and running full memory circuit experiments. For the background behind each, see Experiments and Noise Modeling.

Code-Capacity Noise Modeling

This example implements a code-capacity noise experiment: random X/Z errors are applied directly to the data qubits and decoded with a single-error look-up table. See Code-Capacity Noise Modeling for more details.

CUDA-Q QEC Implementation

Here’s how to use CUDA-Q QEC to perform a code capacity noise model experiment in both Python and C++:

import numpy as np
import cudaq_qec as qec

# Get a QEC code
steane = qec.get_code("steane")

# Get the parity check matrix of a code
# Can get the full code, or for CSS codes
# just the X or Z component
Hz = steane.get_parity_z()
print(f"Hz:\n{Hz}")
observable = steane.get_observables_z()
print(f"observable:\n{observable}")

# error probabily
p = 0.1
# Get a decoder
decoder = qec.get_decoder("single_error_lut", Hz)

# Perform a code capacity noise model numerical experiment
nShots = 10
nLogicalErrors = 0
for i in range(nShots):
    print(f"shot: {i}")

    # Generate noisy data
    data = qec.generate_random_bit_flips(Hz.shape[1], p)
    print(f"data: {data}")

    # Calculate which syndromes are flagged.
    syndrome = Hz @ data % 2
    print(f"syndrome: {syndrome}")

    # Decode the syndrome to predict what happened to the data
    results = decoder.decode(syndrome)
    convergence = results.converged
    result = results.result
    data_prediction = np.array(result, dtype=np.uint8)
    print(f"data_prediction: {data_prediction}")

    # See if this prediction flipped the observable
    predicted_observable = observable @ data_prediction % 2
    print(f"predicted_observable: {predicted_observable}")

    # See if the observable was actually flipped
    actual_observable = observable @ data % 2
    print(f"actual_observable: {actual_observable}")
    if (predicted_observable != actual_observable):
        nLogicalErrors += 1

# Count how many shots the decoder failed to correct the errors
print(f"{nLogicalErrors} logical errors in {nShots} shots\n")

# Can also generate syndromes and data from a single line with:
syndromes, data = qec.sample_code_capacity(Hz, nShots, p)
print("From sample function:")
print("syndromes:\n", syndromes)
print("data:\n", data)
// This example shows the primary cudaq::qec types:
// decoder, code
//
// Compile and run with
// nvq++ --target=stim -lcudaq-qec -lcudaq-qec-decoders code_capacity_noise.cpp
// ./a.out

#include <algorithm>
#include <cmath>
#include <random>

#include "cudaq.h"
#include "cudaq/qec/decoder.h"
#include "cudaq/qec/experiments.h"

int main() {
  auto steane = cudaq::qec::get_code("steane");
  auto Hz = steane->get_parity_z();
  std::vector<size_t> t_shape = Hz.shape();

  std::cout << "Hz.shape():\n";
  for (size_t elem : t_shape)
    std::cout << elem << " ";
  std::cout << "\n";

  std::cout << "Hz:\n";
  Hz.dump();

  auto Lz = steane->get_observables_x();
  std::cout << "Lz:\n";
  Lz.dump();

  double p = 0.2;
  size_t nShots = 5;
  auto lut_decoder = cudaq::qec::get_decoder("single_error_lut", Hz);

  std::cout << "nShots: " << nShots << "\n";

  // May want a order-2 tensor of syndromes
  // access tensor by stride to write in an entire syndrome
  cudaqx::tensor<uint8_t> syndrome({Hz.shape()[0]});

  int nErrors = 0;
  for (size_t shot = 0; shot < nShots; ++shot) {
    std::cout << "shot: " << shot << "\n";
    auto shot_data = cudaq::qec::generate_random_bit_flips(Hz.shape()[1], p);
    std::cout << "shot data\n";
    shot_data.dump();

    auto observable_z_data = Lz.dot(shot_data);
    observable_z_data = observable_z_data % 2;
    std::cout << "Data Lz state:\n";
    observable_z_data.dump();

    auto syndrome = Hz.dot(shot_data);
    syndrome = syndrome % 2;
    std::cout << "syndrome:\n";
    syndrome.dump();

    auto result = lut_decoder->decode(syndrome);
    cudaqx::tensor<uint8_t> result_tensor;
    // result.result is a std::vector<float_t>, of soft information. We'll
    // convert this to hard information and store as a tensor<uint8_t>.
    cudaq::qec::convert_vec_soft_to_tensor_hard(result.result, result_tensor);
    std::cout << "decode result:\n";
    result_tensor.dump();

    // check observable result
    auto decoded_observable_z = Lz.dot(result_tensor);
    std::cout << "decoded observable:\n";
    decoded_observable_z.dump();

    // check how many observable operators were decoded correctly
    // observable_z_data == decoded_observable_z This maps onto element wise
    // addition (mod 2)
    auto observable_flips = decoded_observable_z + observable_z_data;
    observable_flips = observable_flips % 2;
    std::cout << "Logical errors:\n";
    observable_flips.dump();
    std::cout << "\n";

    // shot counts as a observable error unless all observables are correct
    if (observable_flips.any()) {
      nErrors++;
    }
  }
  std::cout << "Total logical errors: " << nErrors << "\n";

  // Full data gen in function call
  auto [syn, data] = cudaq::qec::sample_code_capacity(Hz, nShots, p);
  std::cout << "Numerical experiment:\n";
  std::cout << "Data:\n";
  data.dump();
  std::cout << "Syn:\n";
  syn.dump();
}

Compile and run with

nvq++ --target=stim -lcudaq-qec -lcudaq-qec-decoders code_capacity_noise.cpp -o code_capacity_noise
./code_capacity_noise

Code Explanation

  1. QEC Code type:
    • CUDA-Q QEC centers around the qec.code type, which contains the data relevant for a given code.

    • In particular, this represents a collection of qubits which represent a single logical qubit.

    • Here we get one of the most well known QEC codes, the Steane code, with the qec.get_code function.

    • We can get the stabilizers from a code with the code.get_stabilizers() function.

    • In this example, we get the parity check matrix of the code. Because the Steane code is a CSS code, we can extract just the Z components of the parity check matrix.

    • Here, we see this matrix has 3 rows and 7 columns, which means there are 7 data qubits (7 possible single bit-flip errors) and 3 Z-stabilizers (parity checks). Note that Z stabilizers check for X type errors.

    • Lastly, we get the logical Z observable for the code. This will allow us to see if the Z observable of our logical qubit has flipped.

  2. Decoder type:
    • A single-error look-up table (LUT) decoder can be acquired with the qec.get_decoder call.

    • Passing in the parity check matrix gives the decoder the required information to associated syndromes with underlying error mechanisms.

    • Once the decode has been constructed, the decoder.decode(syndrome) member function is called, which returns a predicted error given the syndrome.

  3. Noise model:
    • To generate noisy data, we call qec.generate_random_bit_flips(nBits, p) which will return an array of bits, where each bit has probability p to have been flipped into 1, and a 1-p chance to have remained 0.

    • Since we are using the Z parity check matrix H_Z, we want to simulate random X errors on our 7 data qubits.

  4. Logical Errors:
    • Once we have noisy data, we see what the resulting syndromes are by multiplying our noisy data vector with our parity check matrix (mod 2).

    • From this syndrome, we see what errors the decoder predicts occurred in the data.

    • To classify as a logical error, the decoder does not need to exactly identify what happened to the data, but only whether there was a flip in the logical observable.

    • If the decoder guesses this successfully, we have corrected the quantum error. If not, we have incurred a logical error.

  5. Further automation:
    • While this workflow is nice for seeing things step by step, the qec.sample_code_capacity API is provided to generate a batch of noisy data and their corresponding syndromes.

Circuit-level Noise Modeling

This example runs a circuit-level memory experiment, generating syndromes by executing the stabilizer-measurement circuits under depolarizing noise. See Circuit-level Noise Modeling for more details.

CUDA-Q QEC Implementation

Here’s how to use CUDA-Q QEC to perform a circuit-level noise model experiment in both Python and C++:

import numpy as np
import cudaq
import cudaq_qec as qec

# Get a QEC code
cudaq.set_target("stim")
distance = 5
surface_code = qec.get_code("surface_code", distance=distance)

# Get the Z observables.
Lz = surface_code.get_observables_z()
print(f"Lz:\n{Lz}")

nShots = 1000
nRounds = distance

# Uncomment for repeatability
# cudaq.set_random_seed(13)

# error probability
p = 0.001
noise = cudaq.NoiseModel()
noise.add_all_qubit_channel("x", cudaq.Depolarization2(p), 1)

# prepare logical |0> state, tells the sampler to do z-basis experiment
statePrep = qec.operation.prep0
# our expected measurement in this state is 0
expected_value = 0

# Get the detector error model for this circuit.
dem = qec.z_dem_from_memory_circuit(surface_code, statePrep, nRounds, noise)

# For large runs, set verbose to False to suppress output
verbose = nShots <= 10

# Sample the surface code memory circuit with noise on each cx gate.
syndromes, data = qec.z_sample_memory_circuit(surface_code, statePrep, nShots,
                                              nRounds, noise)

if verbose:
    print("From sample function:\n")
    print("syndromes:\n", syndromes)
    print("data:\n", data)

# Get a decoder
decoder = qec.get_decoder("single_error_lut", dem.detector_error_matrix)
nLogicalErrors = 0

# Logical Mz each shot (use Lx if preparing in X-basis)
logical_measurements = (Lz @ data.transpose()) % 2
# only one logical qubit, so do not need the second axis
logical_measurements = logical_measurements.flatten()
if verbose:
    print("LMz:\n", logical_measurements)

dr = decoder.decode_batch(syndromes)
error_predictions = np.array([e.result for e in dr], dtype=np.uint8)
data_predictions = (dem.observables_flips_matrix @ error_predictions.T) % 2

nLogicalErrorsWithoutDecoding = np.sum(logical_measurements)
nLogicalErrorsWithDecoding = np.sum(data_predictions ^ logical_measurements)
print(
    f'Number of logical errors without decoding (out of {nShots} shots): {nLogicalErrorsWithoutDecoding}'
)
print(
    f'Number of logical errors with decoding    (out of {nShots} shots): {nLogicalErrorsWithDecoding}'
)
// Compile and run with:
// nvq++ --target=stim -lcudaq-qec -lcudaq-qec-decoders circuit_level_noise.cpp
// ./a.out

#include "cudaq.h"
#include "cudaq/qec/decoder.h"
#include "cudaq/qec/experiments.h"
#include "cudaq/qec/noise_model.h"

int main() {
  // Choose a QEC code
  auto steane = cudaq::qec::get_code("steane");

  // Access the parity check matrix
  auto H = steane->get_parity();
  std::cout << "H:\n";
  H.dump();

  // Access the logical observables
  auto observables = steane->get_pauli_observables_matrix();
  auto Lz = steane->get_observables_z();

  // Data qubits the logical Z observable is supported on
  std::cout << "Lz:\n";
  Lz.dump();

  // Observables are stacked as Z over X for mat-vec multiplication
  std::cout << "Obs:\n";
  observables.dump();

  // How many shots to run the experiment
  int nShots = 3;
  // For each shot, how many rounds of stabilizer measurements
  int nRounds = 4;

  // can set seed for reproducibility
  // cudaq::set_random_seed(1337);
  cudaq::noise_model noise;

  // Add a depolarization noise channel after each cx gate
  noise.add_all_qubit_channel("x", cudaq::depolarization2(/*probability*/ 0.01),
                              /*numControls*/ 1);

  // Perform a noisy z-basis memory circuit experiment
  auto [syndromes, data] = cudaq::qec::sample_memory_circuit(
      *steane, cudaq::qec::operation::prep0, nShots, nRounds, noise);

  // With noise, many syndromes will flip each QEC cycle, these are the
  // syndrome differences from the previous cycle.
  std::cout << "syndromes:\n";
  syndromes.dump();

  // With noise, Lz will sometimes be flipped
  std::cout << "data:\n";
  data.dump();

  // Use z-measurements on data qubits to determine the logical mz
  // In an x-basis experiment, use Lx.
  auto logical_mz = Lz.dot(data.transpose()) % 2;
  std::cout << "logical_mz each shot:\n";
  logical_mz.dump();

  // Select a decoder
  auto decoder = cudaq::qec::get_decoder("single_error_lut", H);

  // Initialize a pauli_frame to track the logical errors
  cudaqx::tensor<uint8_t> pauli_frame({observables.shape()[0]});

  // Start a loop to count the number of logical errors
  size_t numLerrors = 0;
  for (size_t shot = 0; shot < nShots; ++shot) {
    std::cout << "shot: " << shot << "\n";

    for (size_t round = 0; round < nRounds; ++round) {
      std::cout << "round: " << round << "\n";

      // Access one row of the syndrome tensor
      size_t count = shot * nRounds + round;
      size_t stride = syndromes.shape()[1];
      cudaqx::tensor<uint8_t> syndrome({stride});
      syndrome.borrow(syndromes.data() + stride * count);
      std::cout << "syndrome:\n";
      syndrome.dump();

      // Decode the syndrome
      auto result = decoder->decode(syndrome);
      cudaqx::tensor<uint8_t> result_tensor;
      cudaq::qec::convert_vec_soft_to_tensor_hard(result.result, result_tensor);
      std::cout << "decode result:\n";
      result_tensor.dump();

      // See if the decoded result anti-commutes with observables
      auto decoded_observables = observables.dot(result_tensor);
      std::cout << "decoded observable:\n";
      decoded_observables.dump();

      // update from previous stabilizer round
      pauli_frame = (pauli_frame + decoded_observables) % 2;
      std::cout << "pauli frame:\n";
      pauli_frame.dump();
    }

    // prep0 means we expected to measure out 0.
    uint8_t expected_mz = 0;
    // Apply the pauli frame correction to our logical measurement
    uint8_t corrected_mz = (logical_mz.at({0, shot}) + pauli_frame.at({0})) % 2;

    // Check if Logical_mz + pauli_frame_X = 0?
    std::cout << "Corrected readout: " << +corrected_mz << "\n";
    std::cout << "Expected readout: " << +expected_mz << "\n";
    if (corrected_mz != expected_mz)
      numLerrors++;
    std::cout << "\n";
  }

  std::cout << "numLogicalErrors: " << numLerrors << "\n";
}

Compile and run with

nvq++ --target=stim -lcudaq-qec -lcudaq-qec-decoders circuit_level_noise.cpp -o circuit_level_noise
./circuit_level_noise

Code Explanation

  1. QEC Code and Decoder types:
    • As in the code capacity example, our central objects are the qec.code and qec.decoder types.

  2. Clifford simulation backend:
    • As the size of QEC circuits can grow quite large, Clifford simulation is often the best tool for these simulations.

    • cudaq.set_target("stim") selects the highly performant Stim simulator as the simulation backend.

  3. Noise model:
    • To add noisy gates we use the cudaq.NoiseModel type.

    • CUDA-Q supports the generation of arbitrary noise channels. Here we use a cudaq.Depolarization2 channel to add a depolarization channel.

    • This is added to the CX gate by adding it to the X gate with 1 control.

    • This noisy gate is added to every qubit via the noise.add_all_qubit_channel function.

  4. Getting circuit-level noisy data:
    • The qec.code is the first input parameter here, as the code’s stabilizer_round determines the circuits executed.

    • Each memory circuit runs for an input number of nRounds, which specifies how many stabilizer_round kernels are run.

    • After nRounds the data qubits are measured and the run is over.

    • This is performed nShots number of times.

    • During a shot, each stabilizer round’s syndrome is xor’d against the preceding syndrome, so that we can track a sparser flow of data showing which round each parity check was violated.

    • The first round returns the syndrome as is, as there is nothing preceding to xor against.

  5. Data qubit measurements:
    • The data qubits are only read out after the end of each shot, so there are nShots worth of data readouts.

    • The basis of the data qubit measurements depends on the state preparation used.

    • Z-basis readout when preparing the logical |0> or logical |1> state with the qec.operation.prep0 or qec.operation.prep1 kernels.

    • X-basis readout when preparing the logical |+> or logical |-> state with the qec.operation.prepp or qec.operation.prepm kernels.

  6. Logical Errors:
    • From here, the decoding procedure is again similar to the code capacity case, except that we use a Pauli frame to track errors that happen each QEC cycle.

    • The final values of the Pauli frame tell us how our logical state flipped during the experiment, and what needs to be done to correct it.

    • We compare our known initial state (corrected by the Pauli frame), against our measured data qubits to determine if a logical error occurred.

The CUDA-Q QEC library thus provides a platform for numerical QEC experiments. The qec.code can be used to analyze a variety of QEC codes (both library or user provided), with a variety of decoders (both library or user provided). The CUDA-Q QEC library also provides tools to speed up the automation of generating noisy data and syndromes.

Memory Circuit Experiments

The sample_memory_circuit API runs a memory circuit experiment end to end – preparing a logical state, running rounds of stabilizer measurement under noise, and measuring the data qubits. See Memory Circuit Experiments for more details.

Function Variants

import cudaq
import cudaq_qec as qec

# Use the stim backend for performance in QEC settings
cudaq.set_target("stim")

# Get a code instance
code = qec.get_code("steane")

# Basic memory circuit with |0⟩ state
syndromes, measurements = qec.sample_memory_circuit(
    code,           # QEC code instance
    numShots=1000,  # Number of circuit executions
    numRounds=1     # Number of stabilizer rounds
)

# Memory circuit with custom initial state
syndromes, measurements = qec.sample_memory_circuit(
    code,                     # QEC code instance
    op=qec.operation.prep1,   # Initial state
    numShots=1000,            # Number of shots
    numRounds=1               # Number of rounds
)

# Memory circuit with noise model
noise = cudaq.NoiseModel()
# Configure noise
noise.add_all_qubit_channel("x", cudaq.Depolarization2(0.01), 1)
syndromes, measurements = qec.sample_memory_circuit(
    code,             # QEC code instance
    numShots=1000,    # Number of shots
    numRounds=1,      # Number of rounds
    noise=noise       # Noise model
)
// Basic memory circuit with |0⟩ state
auto [syndromes, measurements] = qec::sample_memory_circuit(
    code,       // QEC code instance
    numShots,   // Number of circuit executions
    numRounds   // Number of stabilizer rounds
);

// Memory circuit with custom initial state
auto [syndromes, measurements] = qec::sample_memory_circuit(
    code,               // QEC code instance
    operation::prep1,   // Initial state preparation
    numShots,           // Number of circuit executions
    numRounds           // Number of stabilizer rounds
);

// Memory circuit with noise model
auto noise_model = cudaq::noise_model();
noise_model.add_channel(...);  // Configure noise
auto [syndromes, measurements] = qec::sample_memory_circuit(
    code,         // QEC code instance
    numShots,     // Number of circuit executions
    numRounds,    // Number of stabilizer rounds
    noise_model   // Noise model to apply
);

Return Values

The functions return a tuple containing:

  1. Syndrome Measurements (tensor<uint8_t>):

    • Shape: (num_shots, num_detectors)

    • Columns follow the layout [ B  S  S    S  B ], where:

      • B (boundary block) = numAncZ = code.get_num_z_stabilizers() for Z-basis preparations (prep0/prep1), or numAncX = code.get_num_x_stabilizers() for X-basis preparations (prepp/prepm)

      • S (inter-round block) = numAncZ + numAncX detectors per round transition (num_rounds - 1 blocks total)

      • Total: num_detectors = 2*B + (num_rounds - 1)*S

    • Values are 0 or 1 representing measurement outcomes

  2. Data Measurements (tensor<uint8_t>):

    • Shape: (num_shots, block_size)

    • Contains final data qubit measurements

    • Used to verify logical state preservation

Example Usage

Example of running a memory experiment:

import cudaq
import cudaq_qec as qec

# Use the stim backend for performance in QEC settings
cudaq.set_target("stim")

# Create code and decoder
code = qec.get_code('steane')
decoder = qec.get_decoder('single_error_lut',
                          code.get_parity())

# Configure noise
noise = cudaq.NoiseModel()
noise.add_all_qubit_channel("x", cudaq.Depolarization2(0.01), 1)

# Run memory experiment
syndromes, measurements = qec.sample_memory_circuit(
    code,
    op=qec.operation.prep0,
    numShots=1000,
    numRounds=10,
    noise=noise
)

# Analyze results
for shot in range(1000):
    # Get syndrome for this shot
    syndrome = syndromes[shot].tolist()

    # Decode syndrome
    result = decoder.decode(syndrome)
    if result.converged:
        # Process correction
        pass
// Compile and run with:
// nvq++ --target=stim -lcudaq-qec -lcudaq-qec-decoders example.cpp
// ./a.out

#include "cudaq.h"
#include "cudaq/qec/decoder.h"
#include "cudaq/qec/experiments.h"
#include "cudaq/qec/noise_model.h"

int main(){
  // Create a Steane code instance
  auto code = cudaq::qec::get_code("steane");

  // Configure noise model
  cudaq::noise_model noise;
  noise.add_all_qubit_channel("x", cudaq::depolarization2(0.1),
                      /*num_controls=*/1);

  // Run memory experiment
  auto [syndromes, data] = cudaq::qec::sample_memory_circuit(
      *code,                          // Code instance
      cudaq::qec::operation::prep0,   // Prepare |0⟩ state
      1000,                           // 1000 shots
      1,                              // 1 rounds
      noise                           // Apply noise
  );

  // Analyze results
  auto decoder = cudaq::qec::get_decoder("single_error_lut", code->get_parity());
  for (std::size_t shot = 0; shot < 1000; shot++) {
    // Get syndrome for this shot
    std::vector<cudaq::qec::float_t> syndrome(syndromes.shape()[1]);
    for (std::size_t i = 0; i < syndrome.size(); i++)
      syndrome[i] = syndromes.at({shot, i});

    // Decode syndrome
    auto results = decoder->decode(syndrome);
    // Process correction
    // ...
  }
}

Additional Noise Models

noise = cudaq.NoiseModel()

# Add multiple error channels
noise.add_all_qubit_channel('h', cudaq.BitFlipChannel(0.001))

# Specify two qubit errors
noise.add_all_qubit_channel("x", cudaq.Depolarization2(p), 1)
cudaq::noise_model noise;

// Add multiple error channels
noise.add_all_qubit_channel(
    "x", cudaq::bit_flip_channel(/*probability*/ 0.01));

// Specify two qubit errors
noise.add_all_qubit_channel(
    "x", cudaq::depolarization2(/*probability*/ 0.01),
    /*numControls*/ 1);