Examples

In this section, we show a basic example on how to use the cuPauliProp API for an end-to-end quantum circuit simulation.

Building code

Assuming your cuPauliProp header files are located in CUPAULIPROP_INCLUDE_DIR and the shared library is located in CUPAULIPROP_LIB_DIR, the example codes can be built via the following commands:

nvcc kicked_ising_example.cpp -I${CUPAULIPROP_INCLUDE_DIR} -L${CUPAULIPROP_LIB_DIR} -lcupauliprop -o kicked_ising_example
nvcc kicked_ising_fused_example.cpp -I${CUPAULIPROP_INCLUDE_DIR} -L${CUPAULIPROP_LIB_DIR} -lcupauliprop -o kicked_ising_fused_example
nvcc kicked_ising_backward_diff_example.cpp -I${CUPAULIPROP_INCLUDE_DIR} -L${CUPAULIPROP_LIB_DIR} -lcupauliprop -o kicked_ising_backward_diff_example
nvcc fused_operators_example.cpp example_utils.cpp -I${CUPAULIPROP_INCLUDE_DIR} -L${CUPAULIPROP_LIB_DIR} -lcupauliprop -o fused_operators_example

When building statically, use the following commands:

nvcc kicked_ising_example.cpp -I${CUPAULIPROP_INCLUDE_DIR} ${CUPAULIPROP_LIB_DIR}/libcupauliprop_static.a -o kicked_ising_example
nvcc kicked_ising_fused_example.cpp -I${CUPAULIPROP_INCLUDE_DIR} ${CUPAULIPROP_LIB_DIR}/libcupauliprop_static.a -o kicked_ising_fused_example
nvcc kicked_ising_backward_diff_example.cpp -I${CUPAULIPROP_INCLUDE_DIR} ${CUPAULIPROP_LIB_DIR}/libcupauliprop_static.a -o kicked_ising_backward_diff_example
nvcc fused_operators_example.cpp example_utils.cpp -I${CUPAULIPROP_INCLUDE_DIR} ${CUPAULIPROP_LIB_DIR}/libcupauliprop_static.a -o fused_operators_example

When using cuQuantum downloaded from the NVIDIA devzone website and extracted in CUQUANTUM_ROOT, the header files are located in ${CUQUANTUM_ROOT}/include and the shared library is located in ${CUQUANTUM_ROOT}/lib. For convenience, we also provide a Makefile and CMakeLists.txt file in the same folder as the example code.

Kicked Ising example

In this example, we demonstrate use of the cuPauliProp API to perform an end-to-end simulation of IBM’s 127 qubit kicked Ising experiment, as presented in Nature volume 618, pages 500–505 (2023). We will specifically simulate a single datum of Figure 4 b, which is the expectation value of observable \(Z_{62}\) after a 20-repetition Trotter circuit. The full demonstration is available at kicked_ising_example.cpp.

Our example code will make use of the below helper macros which handle errors thrown by the CUDA or cuPauliProp APIs.

// ========================================================================
// CUDA and cuPauliProp error handling
// ========================================================================

#define HANDLE_CUPP_ERROR(x)                               \
{                                                          \
  const auto err = x;                                      \
  if (err != CUPAULIPROP_STATUS_SUCCESS)                   \
  {                                                        \
    printf("cuPauliProp error in line %d\n", __LINE__);    \
    fflush(stdout);                                        \
    std::abort();                                          \
  }                                                        \
};


#define HANDLE_CUDA_ERROR(x)                                \
{                                                           \
  const auto err = x;                                       \
  if (err != cudaSuccess)                                   \
  {                                                         \
    const char * error = cudaGetErrorString(err);           \
    printf("CUDA Error: %s in line %d\n", error, __LINE__); \
    fflush(stdout);                                         \
    std::abort();                                           \
  }                                                         \
};

We next decide the memory budget for our simulation. It will turn out that for our chosen combination of simulated circuit, output observable and truncation parameters, only around 60 MiB of total device memory is needed. We will demonstrate how to dedicate all available device memory for simulation however, which this example optionally performs when USE_MAX_VRAM = true.

// ========================================================================
// Memory usage
// ========================================================================

// Each Pauli expansion has two pre-allocated GPU buffers, storing packed
// integers (which encode Pauli strings) and corresponding coefficients.
// As much memory can be dedicated as your hardware allows, while the min-
// imum required is specific and very sensitive to the simulated circuit,
// studied observable, and the chosen truncation hyperparameters.
// Some operations also require additional workspace memory which is also
// ideally pre-allocated, and can be established using the API 'Prepare'
// functions (e.g. cupaulipropPauliExpansionViewPrepareTraceWithZeroState).
// In this demo, we dedicate either (a percentage of) the entirety of GPU
// memory uniformly between the required memory buffers, or instead use a
// fixed hardcoded amount which has been prior tested to be consistent with
// our other simulation parameters (like truncations); these choices are
// toggled via USE_MAX_VRAM below.

// =true to use MAX_VRAM_PERCENT of VRAM, and =false to use fixed memories below
bool USE_MAX_VRAM = false;
double MAX_VRAM_PERCENT = 90; // 0-100%

size_t FIXED_EXPANSION_PAULI_MEM = 16 * (1LLU << 20); // bytes = 16 MiB
size_t FIXED_EXPANSION_COEF_MEM  =  4 * (1LLU << 20); // bytes = 4  MiB
size_t FIXED_WORKSPACE_MEM       = 20 * (1LLU << 20); // bytes = 20 MiB

We here scaffold the quantum circuit we will subsequently instantiate. In short, it is a Trotter circuit encoding the dynamics of a kicked Ising system, with a heavy-hex topology matching the connectivity of IBM’s 127 qubit Eagle device.

// ========================================================================
// Circuit preparation (Trotterised Ising on IBM heavy hex topology)
// ========================================================================

// This demo simulates the circuits experimentally executed by IBM in article
// 'Nature volume 618, pages 500–505 (2023)'. This is a circuit Trotterising
// the evolution operator of a 2D transverse-field Ising model, but where the
// prescribed ZZ rotations have a fixed angle of -pi/2, and where the X angles
// are arbitrarily set/swept; later, we will fix the X angle to be pi/4. The
// Hamiltonian ZZ interactions are confined to a heavy-hex topology, matching
// the connectivity of the IBM Eagle processor 'ibm_kyiv', as we fix below.

const int NUM_CIRCUIT_QUBITS = 127;
const int NUM_ROTATIONS_PER_LAYER = 48;
const int NUM_PAULIS_PER_X_ROTATION = 1;
const int NUM_PAULIS_PER_Z_ROTATION = 2;

const double PI = 3.14159265358979323846;
const double ZZ_ROTATION_ANGLE = - PI / 2.0;

// Indices of ZZ-interacting qubits which undergo the first (red) Trotter round
const int32_t ZZ_QUBITS_RED[NUM_ROTATIONS_PER_LAYER][NUM_PAULIS_PER_Z_ROTATION] = {
  {  2,   1},  { 33,  39}, { 59,  60}, { 66,  67}, { 72,  81}, {118, 119},
  { 21,  20},  { 26,  25}, { 13,  12}, { 31,  32}, { 70,  74}, {122, 123},
  { 96,  97},  { 57,  56}, { 63,  64}, {107, 108}, {103, 104}, { 46,  45},
  { 28,  35},  {  7,   6}, { 79,  78}, {  5,   4}, {109, 114}, { 62,  61},
  { 58,  71},  { 37,  52}, { 76,  77}, {  0,  14}, { 36,  51}, {106, 105},
  { 73,  85},  { 88,  87}, { 68,  55}, {116, 115}, { 94,  95}, {100, 110},
  { 17,  30},  { 92, 102}, { 50,  49}, { 83,  84}, { 48,  47}, { 98,  99},
  {  8,   9},  {121, 120}, { 23,  24}, { 44,  43}, { 22,  15}, { 53,  41}
};

// Indices of ZZ-interacting qubits which undergo the second (blue) Trotter round
const int32_t ZZ_QUBITS_BLUE[NUM_ROTATIONS_PER_LAYER][NUM_PAULIS_PER_Z_ROTATION] = {
  { 53,  60}, {123, 124}, { 21,  22}, { 11,  12}, { 67,  68}, {  2,   3},
  { 66,  65}, {122, 121}, {110, 118}, {  6,   5}, { 94,  90}, { 28,  29},
  { 14,  18}, { 63,  62}, {111, 104}, {100,  99}, { 45,  44}, {  4,  15},
  { 20,  19}, { 57,  58}, { 77,  71}, { 76,  75}, { 26,  27}, { 16,   8},
  { 35,  47}, { 31,  30}, { 48,  49}, { 69,  70}, {125, 126}, { 89,  74},
  { 80,  79}, {116, 117}, {114, 113}, { 10,   9}, {106,  93}, {101, 102},
  { 92,  83}, { 98,  91}, { 82,  81}, { 54,  64}, { 96, 109}, { 85,  84},
  { 87,  86}, {108, 112}, { 34,  24}, { 42,  43}, { 40,  41}, { 39,  38}
};

// Indices of ZZ-interacting qubits which undergo the third (green) Trotter round
const int32_t ZZ_QUBITS_GREEN[NUM_ROTATIONS_PER_LAYER][NUM_PAULIS_PER_Z_ROTATION] = {
  { 10,  11}, { 54,  45}, {111, 122}, { 64,  65}, { 60,  61}, {103, 102},
  { 72,  62}, {  4,   3}, { 33,  20}, { 58,  59}, { 26,  16}, { 28,  27},
  {  8,   7}, {104, 105}, { 73,  66}, { 87,  93}, { 85,  86}, { 55,  49},
  { 68,  69}, { 89,  88}, { 80,  81}, {117, 118}, {101, 100}, {114, 115},
  { 96,  95}, { 29,  30}, {106, 107}, { 83,  82}, { 91,  79}, {  0,   1},
  { 56,  52}, { 90,  75}, {126, 112}, { 36,  32}, { 46,  47}, { 77,  78},
  { 97,  98}, { 17,  12}, {119, 120}, { 22,  23}, { 24,  25}, { 43,  34},
  { 42,  41}, { 40,  39}, { 37,  38}, {125, 124}, { 50,  51}, { 18,  19}
};

We instantiate the circuit as a std::vector of cupaulipropQuantumOperator_t types. In this demonstration, all circuit gates are Pauli rotations of 1 or 2 qubits, chiefly

\[\begin{split}Rx(\theta_X) = \exp\left(-i\theta_X/2 \, \hat{X}\right), \\ Rzz(\theta_Z) = \exp\left(-i\theta_Z/2 \, \hat{Z}\otimes\hat{Z}\right).\end{split}\]

We fix \(\theta_Z\) to be \(-\pi/2\), and will later substitute a value for the X-rotation strength \(\theta_X\).

// ========================================================================
// Circuit construction
// ========================================================================

// Each 'step' of the Trotter circuit alternates a layer of single-qubit X
// rotations on every qubit, then a sequence of two-qubit Z rotations on the
// heavy-hex topology, upon qubit pairs in the red, blue and green lists 
// above. Note that ZZ rotations about -pi/2 are actually Clifford, though
// we still here treat them like a generic Pauli rotation. The functions
// below construct a Trotter circuit with a variable number of steps.


std::vector<cupaulipropQuantumOperator_t> getXRotationLayer(
  cupaulipropHandle_t handle, double xRotationAngle
) {  
  std::vector<cupaulipropQuantumOperator_t> layer(NUM_CIRCUIT_QUBITS);

  const cupaulipropPauliKind_t paulis[NUM_PAULIS_PER_X_ROTATION] = {CUPAULIPROP_PAULI_X};

  for (int32_t i=0; i<NUM_CIRCUIT_QUBITS; i++) {
    HANDLE_CUPP_ERROR(cupaulipropCreatePauliRotationGateOperator(
      handle, xRotationAngle, NUM_PAULIS_PER_X_ROTATION, &i, paulis, &layer[i]));
  }

  return layer;
}


std::vector<cupaulipropQuantumOperator_t> getZZRotationLayer(
  cupaulipropHandle_t handle,
  const int32_t topology[NUM_ROTATIONS_PER_LAYER][NUM_PAULIS_PER_Z_ROTATION]
) {
  std::vector<cupaulipropQuantumOperator_t> layer(NUM_ROTATIONS_PER_LAYER);

  const cupaulipropPauliKind_t paulis[NUM_PAULIS_PER_Z_ROTATION] = {
    CUPAULIPROP_PAULI_Z, CUPAULIPROP_PAULI_Z};

  for (uint32_t i=0; i<NUM_ROTATIONS_PER_LAYER; i++) {
    HANDLE_CUPP_ERROR(cupaulipropCreatePauliRotationGateOperator(
      handle, ZZ_ROTATION_ANGLE, NUM_PAULIS_PER_Z_ROTATION, topology[i], paulis, &layer[i]));
  }

  return layer;
}


std::vector<cupaulipropQuantumOperator_t> getIBMHeavyHexIsingCircuit(
  cupaulipropHandle_t handle, double xRotationAngle, int numTrotterSteps
) {
  std::vector<cupaulipropQuantumOperator_t> circuit;

  for (int n=0; n<numTrotterSteps; n++) {
    auto layerX       = getXRotationLayer (handle, xRotationAngle);
    auto layerRedZZ   = getZZRotationLayer(handle, ZZ_QUBITS_RED);
    auto layerBlueZZ  = getZZRotationLayer(handle, ZZ_QUBITS_BLUE);
    auto layerGreenZZ = getZZRotationLayer(handle, ZZ_QUBITS_GREEN);
    
    circuit.insert(circuit.end(), layerX.begin(),       layerX.end());
    circuit.insert(circuit.end(), layerRedZZ.begin(),   layerRedZZ.end());
    circuit.insert(circuit.end(), layerBlueZZ.begin(),  layerBlueZZ.end());
    circuit.insert(circuit.end(), layerGreenZZ.begin(), layerGreenZZ.end());
  }

  return circuit;
}

This demonstration will simulate the circuit in the Heisenberg picture. As such, our input Pauli expansion will be initialized to an experimentally measured observable operator. We’ll later choose observable \(Z_{62}\), but we here define a function which encodes any Pauli string as a sequence of “packed integers” as accepted by the cuPauliProp API.

// ========================================================================
// Observable preparation
// ========================================================================

// This demo simulates the IBM circuit via back-propagating the measurement
// observable through the adjoint circuit. As such, we encode the measured
// observable into our initial Pauli expansion, in the format recognised by
// cuPauliProp. Pauli strings are represented with "packed integers" wherein
// every bit encodes a Pauli operator upon a corresponding qubit. We maintain
// two masks which separately encode the position of X and Z Pauli operators,
// indicated by a set bit at the qubit index, with a common set bit encoding
// a Y Pauli operator. Simulating more qubits than exist bits in the packed
// integer type (64) requires using multiple packed integers for each X and Z
// mask. We store a Pauli string's constituent X and Z masks contiguously in
// a single array, where the final mask of each per-string is padded with zero
// bits to be an integer multiple of the packed integer size (64 bits).

// The below function accepts a single Pauli string (i.e. a tensor product of
// the given Pauli operators at the specified qubit indices) and returns the
// sequence of packed integers which encode it as per the cuPauliProp API;
// this sequence can be copied directly to the GPU buffer of a Pauli expansion.

std::vector<cupaulipropPackedIntegerType_t> getPauliStringAsPackedIntegers(
  std::vector<cupaulipropPauliKind_t> paulis, 
  std::vector<uint32_t> qubits
) {
  assert(paulis.size() == qubits.size());
  assert(*std::max_element(qubits.begin(), qubits.end()) < NUM_CIRCUIT_QUBITS);

  int32_t numPackedInts;
  HANDLE_CUPP_ERROR(cupaulipropGetNumPackedIntegers(NUM_CIRCUIT_QUBITS, &numPackedInts));

  // A single Pauli string is composed of separate X and Z masks, one after the other
  std::vector<cupaulipropPackedIntegerType_t> out(numPackedInts * 2, 0);
  auto xPtr = &out[0];
  auto zPtr = &out[numPackedInts];

  // Process one input (pauli, qubit) pair at a time
  for (auto i=0; i<qubits.size(); i++) {

    // The qubit corresponds to a specific bit of a specific packed integer
    auto numBitsPerPackedInt = 8 * sizeof(cupaulipropPackedIntegerType_t);
    auto intInd = qubits[i] / numBitsPerPackedInt;
    auto bitInd = qubits[i] % numBitsPerPackedInt;

    // Overwrite a bit of either the X or Z masks (or both when pauli==Y)
    if (paulis[i] == CUPAULIPROP_PAULI_X || paulis[i] == CUPAULIPROP_PAULI_Y)
      xPtr[intInd] = xPtr[intInd] | (1ULL << bitInd);
    if (paulis[i] == CUPAULIPROP_PAULI_Z || paulis[i] == CUPAULIPROP_PAULI_Y)
      zPtr[intInd] = zPtr[intInd] | (1ULL << bitInd);
  }

  return out;
}

With all convenience functions defined, we are ready to begin the main control flow.

// ========================================================================
// Main
// ========================================================================

// Simulation of the IBM utility experiment proceeds as follows. We setup the
// cuPauliProp library, attaching a new stream, then proceed to creating two
// Pauli expansions (since the API is out-of-place, as elaborated upon below).
// One expansion is initialised to the measured observable of the IBM circuit.
// We prepare workspace memory, fix truncation hyperparameters, then create
// the circuit as a list of cuPauliProp operators. We process the circuit in
// reverse, adjointing each operation, applied upon a newly prepared view of
// the input expansion, each time checking our dynamically growing memory costs
// have not exceeded our budgets. Thereafter we compute the overlap between the
// final back-propagated observable and the experimental initial state (the all-
// zero state), producing an estimate of the experimental expectation value.
// Finally, we free all allocated memory like good citizens.


int main(int argc, char** argv) {
  std::cout << "cuPauliProp IBM Heavy-hex Ising Example" << std::endl;
  std::cout << "=======================================" << std::endl << std::endl;

Our first chore is to setup the cuPauliProp library. For simplicity, we choose not to attach any stream, and therefore use the default stream.

  // ========================================================================
  // Library setup
  // ========================================================================

  int deviceId = 0;
  HANDLE_CUDA_ERROR(cudaSetDevice(deviceId));

  // cuPauliProp operations accept a cudaStream_t argument for asynchronous usage
  // We use the default stream (0) in this example
  cudaStream_t stream = 0;
  cupaulipropHandle_t handle;
  HANDLE_CUPP_ERROR(cupaulipropCreate(&handle));

We here choose to use either our pre-decided memory budget, or optionally consult the available device memory and distribute it between the needed data structures. In total, we will reserve memory for two Pauli expansions (each containing a separate Pauli string and coefficient buffer), and a workspace buffer. The needed workspace buffer size is informed by both the capacity of the expansions and the API operations to be later performed. While it can be precomputed using the API “prepare” functions, such as cupaulipropPauliExpansionViewPrepareOperatorApplication(), we here merely reserve as much memory for the workspace as needed by a Pauli expansion, for simplicity.

  // ========================================================================
  // Decide memory usage
  // ========================================================================

  // As outlined in the 'Memory usage' section above, we here either uniformly
  // allocate all (or a high percentage of) available memory between the needed
  // memory buffers, or use the pre-decided fixed values. This demo will create
  // a total of two Pauli expansions (each of which accepts two separate buffers
  // to store Pauli strings and their corresponding coefficients; these have
  // different sizes) and one workspace, hence we arrange for an allocation of
  // five buffers in total.

  size_t expansionPauliMem;
  size_t expansionCoefMem;
  size_t workspaceMem;
  size_t totalUsedMem;

  if (USE_MAX_VRAM) {

    // Find usable device memory
    size_t totalFreeMem, totalGlobalMem;
    HANDLE_CUDA_ERROR(cudaMemGetInfo(&totalFreeMem, &totalGlobalMem));
    size_t totalUsableMem = static_cast<size_t>(totalFreeMem * MAX_VRAM_PERCENT/100);

    // Divide it between the three instances (two expansions, one workspace)
    size_t instanceMem = totalUsableMem / 3;

    // Determine the ideal ratio between an expansion's Pauli and coef buffers
    int32_t numPackedInts;
    HANDLE_CUPP_ERROR(cupaulipropGetNumPackedIntegers(NUM_CIRCUIT_QUBITS, &numPackedInts));
    size_t pauliMemPerTerm = 2 * numPackedInts * sizeof(cupaulipropPackedIntegerType_t);
    size_t coefMemPerTerm = sizeof(double);
    size_t totalMemPerTerm = pauliMemPerTerm + coefMemPerTerm;

    expansionPauliMem = (instanceMem * pauliMemPerTerm) / totalMemPerTerm;
    expansionCoefMem  = (instanceMem * coefMemPerTerm ) / totalMemPerTerm;
    workspaceMem      = instanceMem;

    totalUsedMem = 2*expansionPauliMem + 2*expansionCoefMem + workspaceMem;
    std::cout << "Dedicated memory: " << MAX_VRAM_PERCENT << "% of " << totalFreeMem;
    std::cout << " B free = " << totalUsedMem << " B, divided into..." << std::endl;

  } else {

    // Use pre-decided buffer sizes
    expansionPauliMem = FIXED_EXPANSION_PAULI_MEM;
    expansionCoefMem  = FIXED_EXPANSION_COEF_MEM;
    workspaceMem      = FIXED_WORKSPACE_MEM;

    totalUsedMem = 2*expansionPauliMem + 2*expansionCoefMem + workspaceMem;
    std::cout << "Dedicated memory: " << totalUsedMem << " B = 60 MiB, divided into..." << std::endl;
  }

  std::cout << "  expansion Pauli buffer: " << expansionPauliMem << " B" << std::endl;
  std::cout << "  expansion coef buffer:  " << expansionCoefMem << " B" << std::endl;
  std::cout << "  workspace buffer:       " << workspaceMem << " B\n" << std::endl;

It is our responsibility to allocate and initialise the device memory before passing it to the cuPauliProp Pauli expansion constructor. We create an “input” Pauli expansion initialised to the single Pauli string \(Z_{62}\), and an “output” Pauli expansion which is initially empty. Notice that we indicate to cupaulipropCreatePauliExpansion() that the pre-initialised input strings (of which there is one) happen to be sorted and unique, which cuPauliProp can internally later leverage for potential speedup.

  // ========================================================================
  // Pauli expansion preparation
  // ========================================================================

  // Create buffers for two Pauli expansions, which will serve as 'input' and
  // 'output' to the out-of-place cuPauliProp API. Note that the capacities of
  // these buffers constrain the maximum number of Pauli strings maintained
  // during simulation, and ergo inform the accuracy of the simulation. The
  // sufficient buffer sizes are specific to the simulated system, and we here
  // choose a surprisingly small capacity as admitted by the studied circuit.

  void * d_inExpansionPauliBuffer;
  void * d_outExpansionPauliBuffer;
  void * d_inExpansionCoefBuffer;
  void * d_outExpansionCoefBuffer;
  HANDLE_CUDA_ERROR(cudaMalloc(&d_inExpansionPauliBuffer,  expansionPauliMem));
  HANDLE_CUDA_ERROR(cudaMalloc(&d_inExpansionCoefBuffer,   expansionCoefMem));
  HANDLE_CUDA_ERROR(cudaMalloc(&d_outExpansionPauliBuffer, expansionPauliMem));
  HANDLE_CUDA_ERROR(cudaMalloc(&d_outExpansionCoefBuffer,  expansionCoefMem));

  // Prepare the X and Z masks which encode the experimental observable Z_62,
  // which has a coefficient of unity, as seen in Figure 4. b) of the IBM work.
  std::cout << "Observable: Z_62\n" << std::endl;
  int64_t numObservableTerms = 1;
  double observableCoef = 1.0;
  std::vector<cupaulipropPauliKind_t> observablePaulis = {CUPAULIPROP_PAULI_Z};
  std::vector<uint32_t> observableQubits = {62};

  // Overwrite the 'input' Pauli expansion buffers with the observable data
  auto observablePackedInts = getPauliStringAsPackedIntegers(observablePaulis, observableQubits);
  size_t numObservableBytes = observablePackedInts.size() * sizeof(observablePackedInts[0]);
  HANDLE_CUDA_ERROR(cudaMemcpy(
    d_inExpansionPauliBuffer, observablePackedInts.data(), numObservableBytes, cudaMemcpyHostToDevice));
  HANDLE_CUDA_ERROR(cudaMemcpy(
    d_inExpansionCoefBuffer, &observableCoef, sizeof(observableCoef), cudaMemcpyHostToDevice));

  // Create two Pauli expansions, which will serve as 'input and 'output' to the API.
  // Because we begin from a real observable coefficient, and our circuit is completely
  // positive and trace preserving, it is sufficient to use strictly real coefficients
  // in our expansions, informing dataType below. We indicate that the single prepared
  // term in the input expansion is technically unique, and the terms sorted, which
  // permits cuPauliProp to use automatic optimisations during simulation.

  cupaulipropPauliExpansion_t inExpansion;
  cupaulipropPauliExpansion_t outExpansion;

  cupaulipropSortOrder_t sortOrder = CUPAULIPROP_SORT_ORDER_NONE;
  int32_t hasDuplicates = 0;  // isUnique = !hasDuplicates
  cudaDataType_t dataType = CUDA_R_64F;

  HANDLE_CUPP_ERROR(cupaulipropCreatePauliExpansion( // init to above
    handle, NUM_CIRCUIT_QUBITS,
    d_inExpansionPauliBuffer, expansionPauliMem,
    d_inExpansionCoefBuffer,  expansionCoefMem,
    dataType, numObservableTerms, sortOrder, hasDuplicates, 
    &inExpansion));
  HANDLE_CUPP_ERROR(cupaulipropCreatePauliExpansion( // init to empty
    handle, NUM_CIRCUIT_QUBITS,
    d_outExpansionPauliBuffer, expansionPauliMem,
    d_outExpansionCoefBuffer,  expansionCoefMem,
    dataType, 0, CUPAULIPROP_SORT_ORDER_NONE, 0,
    &outExpansion));

  

We must also ourselves allocate device memory for the workspace.

  // ========================================================================
  // Workspace preparation
  // ========================================================================

  // Some API functions require additional workspace memory which we bind to a
  // workspace descriptor. Ordinarily we use the 'Prepare' functions to precisely
  // bound upfront the needed workspace memory, but in this simple demo, we
  // instead use a workspace memory which we prior know to be sufficient.

  // Create a workspace
  cupaulipropWorkspaceDescriptor_t workspace;
  HANDLE_CUPP_ERROR(cupaulipropCreateWorkspaceDescriptor(handle, &workspace));

  void* d_workspaceBuffer;
  HANDLE_CUDA_ERROR(cudaMalloc(&d_workspaceBuffer, workspaceMem));
  HANDLE_CUPP_ERROR(cupaulipropWorkspaceSetMemory(
    handle, workspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    d_workspaceBuffer, workspaceMem));
  
  // Note that the 'prepare' functions which check the required workspace memory
  // will detach this buffer, requiring we re-call SetMemory above. We can avoid
  // this repeated re-attachment by use of a second, bufferless workspace descri-
  // ptor which we pass to the 'prepare' functions in lieu of this one.

During subsequent Heisenberg evolution of the observable operator \(Z_{62}\), the number of terms in the Pauli expansion will grow exponentially quickly. To mitigate this growth and keep simulation tractable through the entire circuit, we employ truncation strategies. For the simulated system, it is sufficient to perform truncation after every \(10\) gates, wherein we discard all Pauli strings with corresponding coefficient magnitudes smaller or equal to \(10^{-4}\), or which have a Pauli “weight” (the number of non-identity Paulis therein) exceeding \(8\). This strategy keeps total memory costs below 60 MiB (and so within the capacities of our Pauli expansion and workspace buffers), while ultimately producing an expectation value within the error bars of the error-mitigated experimental value.

  // ========================================================================
  // Truncation parameter preparation
  // ========================================================================

  // The Pauli propagation simulation technique has memory and runtime costs which
  // (for generic circuits) grow exponentially with the circuit length, which we
  // curtail through "truncation"; dynamic discarding of Pauli strings in our
  // expansion which are predicted to contribute negligibly to the final output
  // expectation value. This file demonstrates simultaneous usage of two truncation
  // techniques; discarding of Pauli strings with an absolute coefficient less than
  // 0.0001, or a "Pauli weight" (the number of non-identity operators in the string)
  // exceeding eight. For example, given a Pauli expansion containing:
  //    0.1 XYZXYZII + 1E-5 ZIIIIIII + 0.2 XXXXYYYY,
  // our truncation parameters below would see the latter two strings discarded due
  // to coefficient and weight truncation respectively.

  cupaulipropCoefficientTruncationParams_t coefTruncParams;
  coefTruncParams.cutoff = 1E-4;

  cupaulipropPauliWeightTruncationParams_t weightTruncParams;
  weightTruncParams.cutoff = 8;

  const uint32_t numTruncStrats = 2;
  cupaulipropTruncationStrategy_t truncStrats[] = {
    {
      CUPAULIPROP_TRUNCATION_STRATEGY_COEFFICIENT_BASED,
      &coefTruncParams
    },
    {
      CUPAULIPROP_TRUNCATION_STRATEGY_PAULI_WEIGHT_BASED,
      &weightTruncParams
    }
  };

  // It is not necessary to perform truncation after every gate, since the
  // Pauli expansion size may not have grown substantially, and attempting
  // to truncate may incur superfluous memory enumeration costs. In this
  // demo, we choose to truncate only after every tenth applied gate. Note
  // deferring truncation requires additional expansion memory; choosing to
  // truncate after every gate shrinks this demo's costs to 20 MiB total.
  const int numGatesBetweenTruncations = 10;

  std::cout << "Coefficient truncation threshold:  " << coefTruncParams.cutoff << std::endl;
  std::cout << "Pauli weight truncation threshold: " << weightTruncParams.cutoff << std::endl;
  std::cout << "Truncation performed after every:  " << numGatesBetweenTruncations << " gates\n" << std::endl;

With the stage set, we are ready to begin simulation! We fix the angle of all X-rotation gates in our circuit to \(\theta_X = \pi/4\), and use \(20\) Trotter repetitions. Our simulated system therefore corresponds to the antepenultimate experimental datum in Figure 4. b of Nature volume 618, pages 500–505 (2023). Let operator \(\hat{U}\) represent the action of this full, forward experimental circuit. The original experimentalists performed \(\hat{U}\) upon the fiducial zero-state \(|0\rangle^{\otimes 127}\), and estimated the expectation value of observable \(Z_{62}\) through repeated sampling.

In contrast, we will here work in the Heisenberg picture, and simulate the adjoint (ergo reversed) circuit upon the observable operator \(Z_{62}\). For each gate, we create a view of the entire input Pauli expansion; we confirm that the output Pauli expansion has sufficient capacity to store the (worst case) result of applying the gate upon the view; we reattach the workspace buffer which is detached by that query; we apply the gate and overwrite the output Pauli expansion; and finally, we swap the input and output Pauli expansions so that the updated expansion becomes the next input.

  // ========================================================================
  // Back-propagation of the observable through the circuit
  // ========================================================================

  // We now simulate the observable operator being back-propagated through the
  // adjoint circuit, mapping the input expansion (initialised to Z_62) to a
  // final output expansion containing many weighted Pauli strings. We use the
  // heavy-hex fixed-angle Ising circuit with 20 total repetitions, fixing the
  // angle of the X rotation gates to PI/4. Our simulation therefore corresponds
  // to the middle datum of Fig. 4 b) of the IBM manuscript, for which MPS and
  // isoTNS siulation techniques showed the greatest divergence from experiment.

  double xRotationAngle = PI / 4.;
  int numTrotterSteps = 20;
  auto circuit = getIBMHeavyHexIsingCircuit(handle, xRotationAngle, numTrotterSteps);

  std::cout << "Circuit: 127 qubit IBM heavy-hex Ising circuit, with..." << std::endl;
  std::cout << "  Trotter steps: " << numTrotterSteps << std::endl;
  std::cout << "  Total gates:   " << circuit.size() << std::endl;
  std::cout << "  Rx angle:      " << xRotationAngle << " (i.e. PI/4)\n" << std::endl;

  // Constrain that every intermediate output expansion contains unique Pauli
  // strings (forbidding duplicates), but permit the retained strings to be
  // unsorted. This combination gives cuPauliProp the best chance of automatically
  // selecting optimal internal functions and postconditions for the simulation.
  uint32_t adjoint = true;
  sortOrder = CUPAULIPROP_SORT_ORDER_NONE;  // reuse variable from above
  uint32_t keepDuplicates = false;

  std::cout << "Imposed postconditions:" << std::endl;
  if (sortOrder != CUPAULIPROP_SORT_ORDER_NONE) {
    std::cout << "  Pauli strings will be sorted." << std::endl;
  }
  if (!keepDuplicates) {
    std::cout << "  Pauli strings will be unique." << std::endl;
  }
  if (keepDuplicates && sortOrder == CUPAULIPROP_SORT_ORDER_NONE) {
    std::cout << "No postconditions imposed on Pauli strings." << std::endl;
  }
  std::cout << std::endl;

  // Begin timing before any gates are applied
  HANDLE_CUDA_ERROR(cudaStreamSynchronize(stream));
  auto startTime = std::chrono::steady_clock::now();
  int64_t maxNumTerms = 0;

  // Iterate the circuit in reverse to effect the adjoint of the total circuit
  for (int gateInd=circuit.size()-1; gateInd >= 0; --gateInd) {
    cupaulipropQuantumOperator_t gate = circuit[gateInd];

    // Create a view of the current input expansion, selecting all currently
    // contained terms. For very large systems, we may have alternatively
    // chosen a smaller view of the partial state to work around memory limits.
    cupaulipropPauliExpansionView_t inView;
    int64_t numExpansionTerms;
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetNumTerms(handle, inExpansion, &numExpansionTerms));
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetContiguousRange(
      handle, inExpansion, 0, numExpansionTerms, &inView));

    // Track the intermediate expansion size, for our curiousity
    if (numExpansionTerms > maxNumTerms)
      maxNumTerms = numExpansionTerms;

    // Choose whether or not to perform truncations after this gate
    int numPassedTruncStrats = (gateInd % numGatesBetweenTruncations == 0)? numTruncStrats : 0;

    // Check the expansion and workspace memories needed to apply the current gate
    int64_t reqExpansionPauliMem;
    int64_t reqExpansionCoefMem;
    int64_t reqWorkspaceMem;
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewPrepareOperatorApplication(
      handle, inView, gate, sortOrder, keepDuplicates,
      numPassedTruncStrats, numPassedTruncStrats > 0 ? truncStrats : nullptr,
      workspaceMem,
      &reqExpansionPauliMem, &reqExpansionCoefMem, workspace));
    HANDLE_CUPP_ERROR(cupaulipropWorkspaceGetMemorySize(
      handle, workspace,
      CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
      &reqWorkspaceMem));

    // Verify that our existing buffers and workspace have sufficient memory
    if (reqExpansionPauliMem > static_cast<int64_t>(expansionPauliMem) ||
        reqExpansionCoefMem  > static_cast<int64_t>(expansionCoefMem)  ||
        reqWorkspaceMem      > static_cast<int64_t>(workspaceMem))
    {
      std::cout
        << "Insufficient outExpansion capacity and/or workspace buffer size "
        << "to perform operator application. Exiting..."
        << std::endl;
      std::abort();
    }

    // Beware that cupaulipropPauliExpansionViewPrepareOperatorApplication() above
    // detaches the memory buffer from the workspace, which we here re-attach.
    HANDLE_CUPP_ERROR(cupaulipropWorkspaceSetMemory(
      handle, workspace,
      CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
      d_workspaceBuffer, workspaceMem));

    // Apply the gate upon the prepared view of the input expansion, evolving the
    // Pauli strings pointed to within, truncating the result. The input expansion
    // is unchanged while the output expansion is entirely overwritten.
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewComputeOperatorApplication(
      handle, inView, outExpansion, gate,
      adjoint, sortOrder, keepDuplicates,
      numPassedTruncStrats, numPassedTruncStrats > 0 ? truncStrats : nullptr,
      workspace, stream));

    // Free the temporary view since it points to the old input expansion, whereas
    // we will subsequently treat the modified output expansion as the next input
    HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansionView(inView));

    // Treat outExpansion as the input in the next gate application
    std::swap(inExpansion, outExpansion);
  }

  // Restore outExpansion to being the final output for clarity
  std::swap(inExpansion, outExpansion);

Our desired output quantity is the expectation value of \(Z_{62}\), as given by

\[\langle Z_{62} \rangle = \text{Tr}\left( \hat{U} |0\rangle\langle 0| \hat{U}^\dagger \; \hat{Z}_{62} \right).\]

This is expressible in terms of our final output Pauli expansion \(s_{\text{out}}\) (the Heisenberg-evolved observable) as

\[\langle Z_{62} \rangle = \text{Tr}\left( |0\rangle\langle 0 | \; \hat{U}^\dagger \hat{Z}_{62} \hat{U} \right) = \text{Tr}\left( |0\rangle\langle 0 | \; s_{\text{out}} \right).\]

Therefore we calculate the trace of our output Pauli expansion with the zero state, producing our estimate to the experimental expectation value of \(Z_{62}\).

  // ========================================================================
  // Evaluation of the the expectation value of observable
  // ========================================================================

  // The output expansion is now a proxy for the observable back-propagated
  // through to the front of the circuit (though having discarded components
  // which negligibly influence the subsequent overlap). The expectation value
  // of the IBM experiment is the overlap of the output expansion with the
  // zero state, i.e. Tr(outExpansion * |0><0|), as we now compute.

  // Obtain a view of the full output expansion (we'll free it in 'Clean up')
  cupaulipropPauliExpansionView_t outView;
  int64_t numOutTerms;
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetNumTerms(handle, outExpansion, &numOutTerms));
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetContiguousRange(
    handle, outExpansion, 0, numOutTerms, &outView));

  // Check that the existing workspace memory is sufficient to compute the trace 
  int64_t reqWorkspaceMem;
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewPrepareTraceWithZeroState(
    handle, outView, workspaceMem, workspace));
  HANDLE_CUPP_ERROR(cupaulipropWorkspaceGetMemorySize(
    handle, workspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    &reqWorkspaceMem));
  if (reqWorkspaceMem > static_cast<int64_t>(workspaceMem)) {
    std::cout
      << "Insufficient workspace buffer size to compute the trace. Exiting..."
      << std::endl;
    std::abort();
  }

  // Beware that we must now reattach the buffer to the workspace
  HANDLE_CUPP_ERROR(cupaulipropWorkspaceSetMemory(
    handle, workspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    d_workspaceBuffer, workspaceMem));

  // Compute the trace; the main and final output of this simulation!
  double expecSignificand;
  double expecExponent;
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewComputeTraceWithZeroState(
    handle, outView, &expecSignificand, &expecExponent, workspace, stream));

  // The CUDA stream must be synchronized before accessing the trace output.
  HANDLE_CUDA_ERROR(cudaStreamSynchronize(stream));
  double expec = expecSignificand * std::pow(2.0, expecExponent);

  // End timing after trace is evaluated
  auto endTime = std::chrono::steady_clock::now();
  auto duration = std::chrono::duration_cast<std::chrono::microseconds>(endTime - startTime);
  auto durationSecs = (duration.count() / 1e6);

  std::cout << "Expectation value:       " << expec << std::endl;
  std::cout << "Final number of terms:   " << numOutTerms << std::endl;
  std::cout << "Maximum number of terms: " << maxNumTerms << std::endl;
  std::cout << "Runtime:                 " << durationSecs << " seconds\n" << std::endl;

Finally, we clean up resources and free all device memory reserved during our simulation. This includes any outstanding views, all quantum operators, all expansions and their underlying buffers, and similarly for the workspace. Our final step is to finalise the cuPauliProp library, concluding simulation.

  // ========================================================================
  // Clean up
  // ========================================================================

  HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansionView(outView));
 
  for (auto & gate : circuit) {
    HANDLE_CUPP_ERROR(cupaulipropDestroyOperator(gate));
  }

  HANDLE_CUPP_ERROR(cupaulipropDestroyWorkspaceDescriptor(workspace));
  HANDLE_CUDA_ERROR(cudaFree(d_workspaceBuffer));

  HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansion(inExpansion));
  HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansion(outExpansion));

  HANDLE_CUDA_ERROR(cudaFree(d_inExpansionPauliBuffer));
  HANDLE_CUDA_ERROR(cudaFree(d_outExpansionPauliBuffer));
  HANDLE_CUDA_ERROR(cudaFree(d_inExpansionCoefBuffer));
  HANDLE_CUDA_ERROR(cudaFree(d_outExpansionCoefBuffer));

  HANDLE_CUPP_ERROR(cupaulipropDestroy(handle));

  return EXIT_SUCCESS;
}

Multiprocess Kicked Ising example

kicked_ising_multiprocess_example.cpp distributes the same Kicked Ising simulation across multiple processes, using either MPI or NCCL for cuPauliProp communication. The circuit, observable encoding, and truncation strategy are unchanged from the preceding example, so both examples produce the same simulation result. This small workload will in fact run more slowly with multiple processes than with a single process; its sole purpose is to demonstrate the multiprocess workflow as simply as possible. This walkthrough therefore focuses on configuring and performing distributed execution.

When MPI_ENABLED is defined, the example initializes MPI and assigns a GPU to each rank. If NCCL_ENABLED is also defined, cuPauliProp uses NCCL for communication; MPI still launches the processes and bootstraps the NCCL communicator. CMake builds kicked_ising_multiprocess_example with MPI support when MPI is available. If both MPI and NCCL are available, CMake also builds kicked_ising_multiprocess_nccl_example from the same source, with both macros defined.

After calling cupaulipropCreate(), each rank must call the collective cupaulipropResetDistributedConfiguration() function to configure the handle for multi-process execution. This call must occur before any rank creates handle-owned objects. On every rank, the CUPAULIPROP_COMM_LIB environment variable must identify the selected cuPauliProp communication plugin. The conditional code below shows the setup for both communication providers:

// ========================================================================
// Main
// ========================================================================

int main(int argc, char** argv) {

  // To be overwritten below
  int rank = 0;
  int numRanks = 1;

  // Both distributed modes use MPI to establish the process topology and map
  // each rank to a GPU. In the NCCL build, MPI performs only launch, bootstrap,
  // and reporting duties; cuPauliProp's data communication uses NCCL instead.
  // Without MPI, the initialized defaults describe one process on device 0.
#if defined(MPI_ENABLED)
  MPI_Init(&argc, &argv);
  HANDLE_MPI_ERROR(MPI_Comm_set_errhandler(MPI_COMM_WORLD, MPI_ERRORS_RETURN));
  HANDLE_MPI_ERROR(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
  HANDLE_MPI_ERROR(MPI_Comm_size(MPI_COMM_WORLD, &numRanks));

  MPI_Comm localComm;
  HANDLE_MPI_ERROR(MPI_Comm_split_type(
    MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, rank, MPI_INFO_NULL, &localComm));
  int localRank = 0;
  int localSize = 1;
  HANDLE_MPI_ERROR(MPI_Comm_rank(localComm, &localRank));
  HANDLE_MPI_ERROR(MPI_Comm_size(localComm, &localSize));
  HANDLE_MPI_ERROR(MPI_Comm_free(&localComm));

  // Oversubscribed nodes map ranks round-robin onto the visible devices.
  int numDevices = 0;
  HANDLE_CUDA_ERROR(cudaGetDeviceCount(&numDevices));
  if (numDevices == 0) {
    fprintf(stderr, "No visible CUDA devices\n");
    abortAllProcesses();
  }
  HANDLE_CUDA_ERROR(cudaSetDevice(localRank % numDevices));
#else
  HANDLE_CUDA_ERROR(cudaSetDevice(0));
#endif

  if (rank == 0) {
    std::cout << "cuPauliProp IBM Heavy-hex Ising Example" << std::endl;
    std::cout << "========================================================"
              << std::endl << std::endl;
  }

  cudaStream_t stream = 0;
  cupaulipropHandle_t handle;
  HANDLE_CUPP_ERROR(cupaulipropCreate(&handle));

  // Select the provider at runtime from the plugin filename. Both
  // communicator variables exist for the whole run because teardown must
  // free whichever one was created after the handle is destroyed.
#if defined(MPI_ENABLED)
  const bool useNccl = commLibRequestsNccl();
  MPI_Comm mpiCommunicator = MPI_COMM_NULL;
#endif
#if defined(NCCL_ENABLED)
  ncclComm_t ncclCommunicator = nullptr;
#endif

#if defined(NCCL_ENABLED)
  if (useNccl) {
    // NCCL communicators require an out-of-band bootstrap: rank 0 creates a
    // unique ID, MPI broadcasts it, and every rank constructs its NCCL
    // endpoint. The resulting opaque ncclComm_t is passed to the NCCL
    // provider.
    ncclUniqueId ncclId;
    if (rank == 0) {
      HANDLE_NCCL_ERROR(ncclGetUniqueId(&ncclId));
    }
    HANDLE_MPI_ERROR(
      MPI_Bcast(&ncclId, sizeof(ncclId), MPI_BYTE, 0, MPI_COMM_WORLD));

    HANDLE_NCCL_ERROR(
      ncclCommInitRank(&ncclCommunicator, numRanks, ncclId, rank));
    HANDLE_CUPP_ERROR(cupaulipropResetDistributedConfiguration(
      handle,
      CUPAULIPROP_DISTRIBUTED_PROVIDER_NCCL,
      &ncclCommunicator,
      sizeof(ncclCommunicator)));
  }
#endif
#if defined(MPI_ENABLED)
  if (!useNccl) {
    // The MPI provider can consume an MPI communicator directly. Duplicate
    // MPI_COMM_WORLD so the example owns the communicator given to cuPauliProp
    // and can preserve it until the handle has been destroyed.
    HANDLE_MPI_ERROR(MPI_Comm_dup(MPI_COMM_WORLD, &mpiCommunicator));
    HANDLE_CUPP_ERROR(cupaulipropResetDistributedConfiguration(
      handle,
      CUPAULIPROP_DISTRIBUTED_PROVIDER_MPI,
      &mpiCommunicator,
      sizeof(mpiCommunicator)));
  }
#endif
#if defined(MPI_ENABLED) && !defined(NCCL_ENABLED)
  if (useNccl) {
    if (rank == 0) {
      fprintf(stderr,
        "CUPAULIPROP_COMM_LIB names an NCCL plugin, but this binary was "
        "built without NCCL support; point it at the MPI plugin or rebuild "
        "with NCCL.\n");
    }
    abortAllProcesses();
  }
#endif

  if (rank == 0) {
    // Report the selected cuPauliProp communication path.
#if defined(MPI_ENABLED)
    if (useNccl) {
      std::cout << "Execution mode: distributed, NCCL backend (MPI bootstrap)" << std::endl;
    } else {
      std::cout << "Execution mode: distributed, MPI backend" << std::endl;
    }
#else
    std::cout << "Execution mode: single-process" << std::endl;
#endif
    std::cout << "Number of ranks: " << numRanks << "\n" << std::endl;
  }

Memory remains process-local. Each process allocates buffers on its assigned GPU for the local expansion data and workspace. These buffers hold the process’s share (a “segment”) of two distributed expansions. Consequently, the sizes reported by the prepare APIs and the memory budget below apply only to the calling process, not to the aggregate memory across all processes.

The fixed configuration reserves 20 MiB for each expansion (16 MiB for Pauli strings and 4 MiB for coefficients) and 128 MiB for workspace, totaling 168 MiB per process. When USE_MAX_VRAM is enabled, each expansion receives approximately 20% of usable memory and workspace receives the remaining approximately 60%. The larger workspace share accommodates overlapping communication staging and deduplication scratch used by distributed operator application.

  // ========================================================================
  // Decide memory usage
  // ========================================================================

  // Process-local GPU memory sizes, to be allocated by every process.
  size_t expansionPauliMem;
  size_t expansionCoefMem;
  size_t workspaceMem;
  size_t totalUsedMem;

  if (USE_MAX_VRAM) {
    size_t totalFreeMem, totalGlobalMem;
    HANDLE_CUDA_ERROR(cudaMemGetInfo(&totalFreeMem, &totalGlobalMem));
    size_t totalUsableMem = static_cast<size_t>(totalFreeMem * MAX_VRAM_PERCENT / 100);

    int32_t numPackedInts;
    HANDLE_CUPP_ERROR(cupaulipropGetNumPackedIntegers(NUM_CIRCUIT_QUBITS, &numPackedInts));
    size_t pauliMemPerTerm = 2 * numPackedInts * sizeof(cupaulipropPackedIntegerType_t);
    size_t coefMemPerTerm = sizeof(double);
    size_t totalMemPerTerm = pauliMemPerTerm + coefMemPerTerm;

    // Reserve approximately 20% for each expansion and 60% for workspace.
    // Distributed redistribution and deduplication require overlapping scratch
    // buffers, so they need a larger share than either local expansion.
    size_t perExpansionMem = totalUsableMem / 5;
    size_t expansionTermCapacity = perExpansionMem / totalMemPerTerm;
    expansionPauliMem = expansionTermCapacity * pauliMemPerTerm;
    expansionCoefMem  = expansionTermCapacity * coefMemPerTerm;
    workspaceMem =
      totalUsableMem - 2 * expansionPauliMem - 2 * expansionCoefMem;

    totalUsedMem = 2 * expansionPauliMem + 2 * expansionCoefMem + workspaceMem;
    if (rank == 0) {
      std::cout << "Dedicated memory per process: " << MAX_VRAM_PERCENT << "% of "
                << totalFreeMem << " B free = " << totalUsedMem << " B = "
                << totalUsedMem / (1LLU << 20) << " MiB, divided into..." << std::endl;
    }
  } else {
    expansionPauliMem = FIXED_EXPANSION_PAULI_MEM;
    expansionCoefMem  = FIXED_EXPANSION_COEF_MEM;
    workspaceMem      = FIXED_WORKSPACE_MEM;

    totalUsedMem = 2 * expansionPauliMem + 2 * expansionCoefMem + workspaceMem;
    if (rank == 0) {
      std::cout << "Dedicated memory per process: " << totalUsedMem << " B = "
                << totalUsedMem / (1LLU << 20) << " MiB, divided into..." << std::endl;
    }
  }

  if (rank == 0) {
    const size_t aggregateUsedMem = totalUsedMem * static_cast<size_t>(numRanks);
    std::cout << "  expansion Pauli buffer: " << expansionPauliMem << " B" << std::endl;
    std::cout << "  expansion coef buffer:  " << expansionCoefMem << " B" << std::endl;
    std::cout << "  workspace buffer:       " << workspaceMem << " B" << std::endl;
    std::cout << "Number of processes: " << numRanks << std::endl;
    std::cout << "Aggregate dedicated memory: " << aggregateUsedMem << " B = "
              << aggregateUsedMem / (1LLU << 20) << " MiB\n" << std::endl;
  }

In a distributed run, cupaulipropCreatePauliExpansion() consumes only the calling process’s local storage and initialized terms. Its buffer pointers, buffer sizes, and numLocalTerms argument are process-local, so their values may differ across ranks. The union of all process-local segments forms the complete Pauli expansion. In this example, only rank 0 initializes its segment with the single \(Z_{62}\) term; every other rank starts with an empty segment.

Other creation arguments describe properties of the full, global, distributed expansion. For example, setting hasDuplicates to zero asserts that no duplicate Pauli strings exist within or across segments. Every process must provide the same value. Somewhat distinctly, the sortOrder argument describes the ordering within each local segment and must also be consistent across processes; however, global sortedness between all ranks is not necessary.

  // ========================================================================
  // Pauli expansion preparation
  // ========================================================================

  void* d_inExpansionPauliBuffer;
  void* d_outExpansionPauliBuffer;
  void* d_inExpansionCoefBuffer;
  void* d_outExpansionCoefBuffer;
  HANDLE_CUDA_ERROR(cudaMalloc(&d_inExpansionPauliBuffer,  expansionPauliMem));
  HANDLE_CUDA_ERROR(cudaMalloc(&d_inExpansionCoefBuffer,   expansionCoefMem));
  HANDLE_CUDA_ERROR(cudaMalloc(&d_outExpansionPauliBuffer, expansionPauliMem));
  HANDLE_CUDA_ERROR(cudaMalloc(&d_outExpansionCoefBuffer,  expansionCoefMem));

  if (rank == 0) {
    std::cout << "Observable: Z_62\n" << std::endl;
  }
  std::vector<cupaulipropPauliKind_t> observablePaulis = {CUPAULIPROP_PAULI_Z};
  std::vector<uint32_t> observableQubits = {62};
  double observableCoef = 1.0;

  int64_t numLocalTerms = 0;
  if (rank == 0) {
    numLocalTerms = 1;
    auto observablePackedInts =
      getPauliStringAsPackedIntegers(observablePaulis, observableQubits);
    size_t numObservableBytes =
      observablePackedInts.size() * sizeof(observablePackedInts[0]);
    HANDLE_CUDA_ERROR(cudaMemcpy(
      d_inExpansionPauliBuffer, observablePackedInts.data(),
      numObservableBytes, cudaMemcpyHostToDevice));
    HANDLE_CUDA_ERROR(cudaMemcpy(
      d_inExpansionCoefBuffer, &observableCoef,
      sizeof(observableCoef), cudaMemcpyHostToDevice));
  }

  cupaulipropPauliExpansion_t inExpansion;
  cupaulipropPauliExpansion_t outExpansion;
  cupaulipropSortOrder_t sortOrder = CUPAULIPROP_SORT_ORDER_NONE;
  int32_t hasDuplicates = 0;
  cudaDataType_t dataType = CUDA_R_64F;

  // Each call creates only this rank's local segment of the distributed
  // expansion. The buffers, their capacities, and the local term count may
  // differ between ranks; together, the local segments form the full expansion.
  HANDLE_CUPP_ERROR(cupaulipropCreatePauliExpansion(
    handle, NUM_CIRCUIT_QUBITS,
    d_inExpansionPauliBuffer, expansionPauliMem,
    d_inExpansionCoefBuffer, expansionCoefMem,
    dataType, numLocalTerms, sortOrder, hasDuplicates,
    &inExpansion));
  HANDLE_CUPP_ERROR(cupaulipropCreatePauliExpansion(
    handle, NUM_CIRCUIT_QUBITS,
    d_outExpansionPauliBuffer, expansionPauliMem,
    d_outExpansionCoefBuffer, expansionCoefMem,
    dataType, 0, CUPAULIPROP_SORT_ORDER_NONE, 0,
    &outExpansion));

The workspace descriptor and its attached allocation are also local to each process. Prepare calls detach the allocation on each rank, so each rank must reattach its local buffer before computation.

  // ========================================================================
  // Workspace preparation
  // ========================================================================

  cupaulipropWorkspaceDescriptor_t workspace;
  HANDLE_CUPP_ERROR(cupaulipropCreateWorkspaceDescriptor(handle, &workspace));

  void* d_workspaceBuffer;
  HANDLE_CUDA_ERROR(cudaMalloc(&d_workspaceBuffer, workspaceMem));
  HANDLE_CUPP_ERROR(cupaulipropWorkspaceSetMemory(
    handle, workspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    d_workspaceBuffer, workspaceMem));

During distributed propagation, cupaulipropPauliExpansionGetNumTerms() returns the term count for only the calling process. cupaulipropPauliExpansionGetContiguousRange() likewise indexes only that process’s local segment. Together, the full-range local views represent the complete distributed expansion.

In a distributed run, the operator prepare and compute calls are collective. Every process must specify the same operator, requested postconditions, and truncation strategies. However, the output-buffer and workspace sizes returned by the prepare call are process-local. Prepare cannot in general know exactly how many resulting terms each process will own. It therefore includes statistical slack for potential inter-process imbalances. Consequently, cupaulipropPauliExpansionViewComputeOperatorApplication() can still return CUPAULIPROP_STATUS_INSUFFICIENT_OUT_EXPANSION in the extremely unlikely event that the actual inter-process imbalance exceeds this prepared capacity, even when the reported memory sizes were provided. This failure status is propagated to every participating process.

The keepDuplicates argument has global meaning: setting it to false, as shown below, requires Pauli strings to be unique both within and across processes. The truncation strategies are passed directly to the operator prepare and compute calls on the iterations where truncation is requested. When deduplication and truncation are both requested, cuPauliProp always completes global deduplication before applying the truncation strategies. At the end of the compute call, the calling process will contain a possibly distinct segment of the resulting output expansion, in outExpansion.

  // ========================================================================
  // Back-propagation of the observable through the circuit
  // ========================================================================

  double xRotationAngle = PI / 4.;
  int numTrotterSteps = 20;
  auto circuit = getIBMHeavyHexIsingCircuit(handle, xRotationAngle, numTrotterSteps);

  if (rank == 0) {
    std::cout << "Circuit: 127 qubit IBM heavy-hex Ising circuit, with..." << std::endl;
    std::cout << "  Trotter steps: " << numTrotterSteps << std::endl;
    std::cout << "  Total gates:   " << circuit.size() << std::endl;
    std::cout << "  Rx angle:      " << xRotationAngle << " (i.e. PI/4)\n" << std::endl;
  }

  uint32_t adjoint = true;
  sortOrder = CUPAULIPROP_SORT_ORDER_NONE;

  // Global deduplication, requested with keepDuplicates=false, can require
  // redistributing Pauli terms between processes. This can significantly
  // increase communication costs and therefore simulation runtime.
  int32_t keepDuplicates = false;

  if (rank == 0) {
    std::cout << "Imposed postconditions:" << std::endl;
    if (sortOrder != CUPAULIPROP_SORT_ORDER_NONE) {
      std::cout << "  Pauli strings will be sorted." << std::endl;
    }
    if (!keepDuplicates) {
      std::cout << "  Pauli strings will be unique." << std::endl;
    }
    if (keepDuplicates && sortOrder == CUPAULIPROP_SORT_ORDER_NONE) {
      std::cout << "No postconditions imposed on Pauli strings." << std::endl;
    }
    std::cout << std::endl;
  }

  // Exclude preceding asynchronous CUDA work and align all ranks before timing.
  HANDLE_CUDA_ERROR(cudaStreamSynchronize(stream));
#if defined(MPI_ENABLED)
  HANDLE_MPI_ERROR(MPI_Barrier(MPI_COMM_WORLD));
#endif
  auto startTime = std::chrono::steady_clock::now();
  int64_t maxNumLocalTerms = 0;

  for (int gateInd = circuit.size() - 1; gateInd >= 0; --gateInd) {
    cupaulipropQuantumOperator_t gate = circuit[gateInd];

    cupaulipropPauliExpansionView_t inView;
    int64_t currentNumLocalTerms;
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetNumTerms(
      handle, inExpansion, &currentNumLocalTerms));
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetContiguousRange(
      handle, inExpansion, 0, currentNumLocalTerms, &inView));

    if (currentNumLocalTerms > maxNumLocalTerms)
      maxNumLocalTerms = currentNumLocalTerms;

    bool applyTruncation = (gateInd % numGatesBetweenTruncations == 0);
    int32_t appliedNumTruncStrats = applyTruncation ? numTruncStrats : 0;
    const cupaulipropTruncationStrategy_t* appliedTruncStrats =
      applyTruncation ? truncStrats : nullptr;
      
    // In multi-process execution, Prepare returns process-local sizing with
    // statistical slack for potential inter-process term imbalances. Unlike
    // in single-process execution, passing workspaces and output expansions
    // of the sizes prescribed by Prepare to the subsequent Compute function
    // does not guarantee it will succeed, though the probability of failure
    // is very small, and vanishes for meaningfully large systems.
    int64_t reqExpansionPauliMem;
    int64_t reqExpansionCoefMem;
    int64_t reqWorkspaceMem;
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewPrepareOperatorApplication(
      handle, inView, gate, sortOrder, keepDuplicates,
      appliedNumTruncStrats, appliedTruncStrats,
      workspaceMem,
      &reqExpansionPauliMem, &reqExpansionCoefMem, workspace));
    HANDLE_CUPP_ERROR(cupaulipropWorkspaceGetMemorySize(
      handle, workspace,
      CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
      &reqWorkspaceMem));

    if (reqExpansionPauliMem > static_cast<int64_t>(expansionPauliMem) ||
        reqExpansionCoefMem  > static_cast<int64_t>(expansionCoefMem)  ||
        reqWorkspaceMem      > static_cast<int64_t>(workspaceMem))
    {
      std::cout
        << "Insufficient outExpansion capacity and/or workspace buffer size "
        << "to perform operator application. Exiting..."
        << std::endl;
      abortAllProcesses();
    }

    HANDLE_CUPP_ERROR(cupaulipropWorkspaceSetMemory(
      handle, workspace,
      CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
      d_workspaceBuffer, workspaceMem));

    const auto operatorStatus =
      cupaulipropPauliExpansionViewComputeOperatorApplication(
        handle, inView, outExpansion, gate,
        adjoint, sortOrder, keepDuplicates,
        appliedNumTruncStrats, appliedTruncStrats,
        workspace, stream);
    if (operatorStatus == CUPAULIPROP_STATUS_INSUFFICIENT_OUT_EXPANSION
        && numRanks > 1 && rank == 0) {
      std::cerr
        << "The actual inter-process term distribution exceeded the "
           "probabilistic capacity reported by PrepareOperatorApplication. "
           "This outcome is extremely unlikely."
        << std::endl;
    }
    HANDLE_CUPP_ERROR(operatorStatus);

    HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansionView(inView));
    std::swap(inExpansion, outExpansion);
  }

  std::swap(inExpansion, outExpansion);

  int64_t numOutTerms;
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetNumTerms(
    handle, outExpansion, &numOutTerms));
  if (numOutTerms > maxNumLocalTerms)
    maxNumLocalTerms = numOutTerms;

  // GetNumTerms reports a rank-local count. Use MPI to find the largest local
  // shard seen in either distributed mode; NCCL builds can also use this reduction
  // because MPI remains active for reporting. In single-process mode the local
  // maximum already is the global maximum.
#if defined(MPI_ENABLED)
  int64_t maxNumLocalTermsAcrossRanks = 0;
  HANDLE_MPI_ERROR(MPI_Allreduce(
    &maxNumLocalTerms, &maxNumLocalTermsAcrossRanks,
    1, MPI_INT64_T, MPI_MAX, MPI_COMM_WORLD));
  maxNumLocalTerms = maxNumLocalTermsAcrossRanks;
#endif

In a distributed run, trace computation is collective. cuPauliProp combines the local contributions and writes the same global significand and exponent on every process. The expectation value therefore requires no explicit MPI reduction. Term-count queries remain local, so the example uses MPI_Allreduce to sum the final local counts for display.

  // ========================================================================
  // Evaluation of the expectation value of the observable
  // ========================================================================

  cupaulipropPauliExpansionView_t outView;
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetContiguousRange(
    handle, outExpansion, 0, numOutTerms, &outView));

  // Sum rank-local output sizes for display in either distributed mode. This
  // explicit MPI reduction is only for the statistic: cuPauliProp's trace call
  // below already combines distributed trace contributions through the selected
  // MPI or NCCL provider. No reduction is needed in single-process mode.
#if defined(MPI_ENABLED)
  int64_t totalNumOutTerms = 0;
  HANDLE_MPI_ERROR(MPI_Allreduce(
    &numOutTerms, &totalNumOutTerms, 1, MPI_INT64_T, MPI_SUM, MPI_COMM_WORLD));
  numOutTerms = totalNumOutTerms;
#endif

  int64_t reqWorkspaceMem;
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewPrepareTraceWithZeroState(
    handle, outView, workspaceMem, workspace));
  HANDLE_CUPP_ERROR(cupaulipropWorkspaceGetMemorySize(
    handle, workspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    &reqWorkspaceMem));
  if (reqWorkspaceMem > static_cast<int64_t>(workspaceMem)) {
    std::cout
      << "Insufficient workspace buffer size to compute the trace. Exiting..."
      << std::endl;
    abortAllProcesses();
  }

  HANDLE_CUPP_ERROR(cupaulipropWorkspaceSetMemory(
    handle, workspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    d_workspaceBuffer, workspaceMem));

  double expecSignificand;
  double expecExponent;
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewComputeTraceWithZeroState(
    handle, outView, &expecSignificand, &expecExponent, workspace, stream));

  // Complete local GPU work and align all ranks before ending timing.
  HANDLE_CUDA_ERROR(cudaStreamSynchronize(stream));
#if defined(MPI_ENABLED)
  HANDLE_MPI_ERROR(MPI_Barrier(MPI_COMM_WORLD));
#endif
  auto endTime = std::chrono::steady_clock::now();
  double expec = expecSignificand * std::pow(2.0, expecExponent);

  auto duration =
    std::chrono::duration_cast<std::chrono::microseconds>(endTime - startTime);
  auto durationSecs = duration.count() / 1e6;

  if (rank == 0) {
    std::cout << "Expectation value:                  " << expec << std::endl;
    std::cout << "Final global number of terms:       " << numOutTerms << std::endl;
    std::cout << "Maximum local terms on any rank:    " << maxNumLocalTerms << std::endl;
    std::cout << "Runtime:                            "
              << durationSecs << " seconds\n" << std::endl;
  }

Finally, each process destroys its local views, expansion handles, workspace, device buffers, and library handle. The MPI-enabled paths keep the MPI or NCCL communicator alive until after cupaulipropDestroy(). They then destroy the selected communicator and finalize MPI.

  // ========================================================================
  // Clean up
  // ========================================================================

  HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansionView(outView));

  for (auto& gate : circuit) {
    HANDLE_CUPP_ERROR(cupaulipropDestroyOperator(gate));
  }

  HANDLE_CUPP_ERROR(cupaulipropDestroyWorkspaceDescriptor(workspace));
  HANDLE_CUDA_ERROR(cudaFree(d_workspaceBuffer));

  HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansion(inExpansion));
  HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansion(outExpansion));

  HANDLE_CUDA_ERROR(cudaFree(d_inExpansionPauliBuffer));
  HANDLE_CUDA_ERROR(cudaFree(d_outExpansionPauliBuffer));
  HANDLE_CUDA_ERROR(cudaFree(d_inExpansionCoefBuffer));
  HANDLE_CUDA_ERROR(cudaFree(d_outExpansionCoefBuffer));

  HANDLE_CUPP_ERROR(cupaulipropDestroy(handle));

  // The provider communicator had to outlive the cuPauliProp handle. Destroy
  // whichever concrete communicator was created during setup, then finalize
  // the MPI process environment used by both distributed modes. The
  // single-process path owns neither resource and therefore has no
  // corresponding cleanup.
#if defined(NCCL_ENABLED)
  if (useNccl) {
    HANDLE_NCCL_ERROR(ncclCommDestroy(ncclCommunicator));
  }
#endif
#if defined(MPI_ENABLED)
  if (mpiCommunicator != MPI_COMM_NULL) {
    HANDLE_MPI_ERROR(MPI_Comm_free(&mpiCommunicator));
  }
  MPI_Finalize();
#endif

  return EXIT_SUCCESS;
}

Kicked Ising fused Clifford example

This companion example simulates the same 127-qubit circuit and observable as kicked_ising_example.cpp, but exploits the fact that every fixed \(Rzz(-\pi/2)\) rotation is Clifford. During Heisenberg back-propagation, the adjoint of each rotation is represented by a \(CX\)-\(S\)-\(CX\) sequence. The three heavy-hex interaction layers are concatenated and applied with one fused Clifford call per Trotter step, while the non-Clifford \(Rx(\pi/4)\) gates are still applied individually.

The example uses the fused Prepare API to verify its fixed output and workspace allocations before computation. Truncation is skipped for the non-branching Clifford sequence and performed only after every tenth \(Rx\) rotation. Consequently, its truncation schedule differs from the original example and can produce a slightly different approximate expectation value. The full source is available in kicked_ising_fused_example.cpp.

Kicked Ising backward differentiation example

This example demonstrates the cuPauliProp facilities for calculating gradients of expectation values with respect to circuit parameters via reverse-mode automatic differentiation (AD). For simplicity, we employ the same circuit as in the above example, parameterised by the \(X\) and \(ZZ\) rotation angles which we label \(\theta_X\) and \(\theta_Z\) respectively. Where the above example computed the expectation value \(\langle Z_{62} \rangle\) for \((\theta_X,\theta_Z) = (\pi/4, -\pi/2)\), we will here compute the gradient of this value at the same coordinate. That is, we seek

\[\nabla_{\theta_X, \theta_Z} \langle Z_{62} \rangle = ( \frac{\partial}{\partial \theta_X} \langle Z_{62} \rangle, \frac{\partial}{\partial \theta_Z} \langle Z_{62} \rangle ).\]

We will skip all code which is similar to the above example, and focus instead on the additional or particular steps involved in use of the AD API. The full demonstration is available in kicked_ising_backward_diff_example.cpp.

Our first algorithmic step is to simulate the non-differentiated circuit and obtain the output expectation value, in an identical fashion to the example above. We evaluate the circuit in the Heisenberg picture, beginning from the \(Z_{62}\) observable and applying the adjoint of the full circuit. Somewhat confusingly, this is the “forward pass” of AD, which just so happens to be a backward enumeration of the experimental circuit. Below, inExpansion is initialised to observable operator \(Z_{62}\).

  // ========================================================================
  // Forward pass: back-propagate observable through adjoint circuit
  // ========================================================================
  auto startTime = std::chrono::high_resolution_clock::now();
  int64_t maxNumTerms = 0;

  for (int i = static_cast<int>(circuit.size()) - 1; i >= 0; --i) {
    int64_t inTerms = 0, reqXZ = 0, reqCoef = 0, reqWs = 0;
    cupaulipropPauliExpansionView_t inView{};
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetNumTerms(handle, inExpansion, &inTerms));
    if (inTerms > maxNumTerms) maxNumTerms = inTerms;
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetContiguousRange(handle, inExpansion, 0, inTerms, &inView));
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewPrepareOperatorApplication(
        handle, inView, circuit[i].op, sortOrder, 0, 2, truncStrats, WORKSPACE_MEM, &reqXZ, &reqCoef, workspace));
    HANDLE_CUPP_ERROR(cupaulipropWorkspaceGetMemorySize(
        handle, workspace, CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH, &reqWs));
    if (reqXZ   > static_cast<int64_t>(EXPANSION_PAULI_MEM) ||
        reqCoef > static_cast<int64_t>(EXPANSION_COEF_MEM)  ||
        reqWs   > static_cast<int64_t>(WORKSPACE_MEM))
    {
      std::cout
          << "Insufficient outExpansion capacity and/or workspace buffer size "
          << "to perform operator application. Exiting..."
          << std::endl;
      std::abort();
    }
    reattachWorkspace(handle, workspace, dWorkspace);
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewComputeOperatorApplication(
        handle, inView, outExpansion, circuit[i].op, 1, sortOrder, 0, 2, truncStrats, workspace, stream));
    HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansionView(inView));
    std::swap(inExpansion, outExpansion);
  }

We then obtain the expectation value \(\langle Z_{62} \rangle\) at \((\theta_X,\theta_Z) = (\pi/4, -\pi/2)\). Notice below that the value is given to us as a separate significand and exponent, which we will here denote as \(s\) and \(p\) respectively, such that \(\langle Z_{62} \rangle = s \times 2^p\).

  // ========================================================================
  // Trace evaluation (expectation value)
  // ========================================================================
  int64_t finalTerms = 0;
  cupaulipropPauliExpansionView_t finalView{};
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetNumTerms(handle, inExpansion, &finalTerms));
  if (finalTerms > maxNumTerms) maxNumTerms = finalTerms;
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetContiguousRange(handle, inExpansion, 0, finalTerms, &finalView));
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewPrepareTraceWithZeroState(handle, finalView, WORKSPACE_MEM, workspace));
  reattachWorkspace(handle, workspace, dWorkspace);
  double s = 0.0, p = 0.0;
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewComputeTraceWithZeroState(handle, finalView, &s, &p, workspace, stream));
  const double expec = s * std::exp2(p);

We now proceed to evaluating the gradient of this quantity. While above we obtained the result of applying the full circuit, computing the gradient will make use of each of the intermediate states reached by evolving the observable through only a partial circuit. We did obtain these intermediate states in-turn but did not store them, which would necessitate as many Pauli expansions as operators in the circuit. Instead, we will leverage that the kicked Ising circuit is unitary and therefore invertible, enabling us to re-obtain intermediate states by applying the inverse of the last applied operator. This trick is an instance of “tape replay” and means that in addition to our input and output expansions used above, we need only maintain two additional Pauli expansions during gradient evaluation, independent of the circuit length. We prepare these initially empty Pauli expansions just like those used above.

  cupaulipropPauliExpansion_t inExpansion{}, outExpansion{}, cot0Expansion{}, cot1Expansion{};
  const auto sortOrder = CUPAULIPROP_SORT_ORDER_NONE;
  HANDLE_CUPP_ERROR(cupaulipropCreatePauliExpansion(
      handle, NUM_CIRCUIT_QUBITS, dInXZ, EXPANSION_PAULI_MEM, dInCoef, EXPANSION_COEF_MEM,
      CUDA_R_64F, 1, sortOrder, 0, &inExpansion));
  HANDLE_CUPP_ERROR(cupaulipropCreatePauliExpansion(
      handle, NUM_CIRCUIT_QUBITS, dOutXZ, EXPANSION_PAULI_MEM, dOutCoef, EXPANSION_COEF_MEM,
      CUDA_R_64F, 0, sortOrder, 0, &outExpansion));
  HANDLE_CUPP_ERROR(cupaulipropCreatePauliExpansion(
      handle, NUM_CIRCUIT_QUBITS, dCot0XZ, EXPANSION_PAULI_MEM, dCot0Coef, EXPANSION_COEF_MEM,
      CUDA_R_64F, 0, sortOrder, 0, &cot0Expansion));
  HANDLE_CUPP_ERROR(cupaulipropCreatePauliExpansion(
      handle, NUM_CIRCUIT_QUBITS, dCot1XZ, EXPANSION_PAULI_MEM, dCot1Coef, EXPANSION_COEF_MEM,
      CUDA_R_64F, 0, sortOrder, 0, &cot1Expansion));

After the forward pass already seen, which computed the output expectation value \(\langle Z_{62} \rangle = s 2^p\), we are ready to proceed with backward differentiation. This involves computing the backward derivative of each API operation involved in evaluating \(\langle Z_{62} \rangle\), the final of which (and ergo our first to now process) was the trace function; that which output \(s\) and \(p\). So we “seed” our AD cotangent vector with the gradient of \(\langle Z_{62} \rangle\) (our “cost function”) with respect to these symbols.

\[\nabla_{s,p} \langle Z_{62} \rangle = ( \frac{\partial}{\partial s} \langle Z_{62} \rangle, \frac{\partial}{\partial p} \langle Z_{62} \rangle ) = ( 2^p, s 2^p \log_e(2) ).\]

And as always, it is prudent to use the prepare functions to ensure that our cotangent expansion has sufficient memory to store the output, which requires subsequently re-attaching the workspace buffer.

  // ========================================================================
  // Backward: seed cotangent from trace
  // ========================================================================
  const double cotS = std::exp2(p);
  const double cotP = s * std::exp2(p) * std::log(2.0);
  int64_t reqCotXZ = 0, reqCotCoef = 0;
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewPrepareTraceWithZeroStateBackwardDiff(
      handle, finalView, WORKSPACE_MEM, &reqCotXZ, &reqCotCoef, workspace));
  if (reqCotXZ   > static_cast<int64_t>(EXPANSION_PAULI_MEM) ||
      reqCotCoef > static_cast<int64_t>(EXPANSION_COEF_MEM))
  {
    std::cout
        << "Insufficient cot0Expansion capacity to differentiate the trace. Exiting..."
        << std::endl;
    std::abort();
  }

This modified cot0Expansion, our cotangent expansion. We now enumerate backward through the operators applied upon our original Pauli expansion, which happens to be the forward order of the experimental circuit, merely because we originally worked in the Heisenberg picture. Each processed operator will contribute to either the \(\theta_X\) or \(\theta_Z\) component of the gradient, which we store separately.

  reattachWorkspace(handle, workspace, dWorkspace);
  HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewComputeTraceWithZeroStateBackwardDiff(
      handle, finalView, &cotS, &cotP, cot0Expansion, workspace, stream));
  HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansionView(finalView));

  // ========================================================================
  // Backward pass: reverse-mode AD with tape replay
  // ========================================================================
  double gradX = 0.0, gradZZ = 0.0;
  auto bwdStartTime = std::chrono::high_resolution_clock::now();

  for (size_t i = 0; i < circuit.size(); ++i) {
    int64_t inTerms = 0;
    cupaulipropPauliExpansionView_t inView{};
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetNumTerms(handle, inExpansion, &inTerms));
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetContiguousRange(handle, inExpansion, 0, inTerms, &inView));

For each operator, we apply its inverse (which for unitary operators is merely the adjoint) in order to restore the intermediate expansion from before that operator was applied (this is the tape replay).

    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewPrepareOperatorApplication(
        handle, inView, circuit[i].op, sortOrder, 0, 2, truncStrats, WORKSPACE_MEM, &reqCotXZ, &reqCotCoef, workspace));
    reattachWorkspace(handle, workspace, dWorkspace);
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewComputeOperatorApplication(
        handle, inView, outExpansion, circuit[i].op, 0, sortOrder, 0, 2, truncStrats, workspace, stream));
    HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansionView(inView));
    std::swap(inExpansion, outExpansion);

We pass this intermediate expansion, along with our cotangent expansion, to the backward derivative function of the operator, obtaining a term of a derivative of \(\langle Z_{62} \rangle\).

    int64_t cotTerms = 0, reqXZ = 0, reqCoef = 0, reqWs = 0;
    cupaulipropPauliExpansionView_t viewIn{}, cotOutView{};
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetNumTerms(handle, inExpansion, &inTerms));
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetContiguousRange(handle, inExpansion, 0, inTerms, &viewIn));
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetNumTerms(handle, cot0Expansion, &cotTerms));
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionGetContiguousRange(handle, cot0Expansion, 0, cotTerms, &cotOutView));
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewPrepareOperatorApplicationBackwardDiff(
        handle, viewIn, cotOutView, circuit[i].op, sortOrder, 0, 2, truncStrats, WORKSPACE_MEM, &reqXZ, &reqCoef, workspace));
    HANDLE_CUPP_ERROR(cupaulipropWorkspaceGetMemorySize(
        handle, workspace, CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH, &reqWs));
    if (reqXZ   > static_cast<int64_t>(EXPANSION_PAULI_MEM) ||
        reqCoef > static_cast<int64_t>(EXPANSION_COEF_MEM)  ||
        reqWs   > static_cast<int64_t>(WORKSPACE_MEM))
    {
      std::cout
          << "Insufficient cot1Expansion capacity and/or workspace buffer size "
          << "to differentiate operator application. Exiting..."
          << std::endl;
      std::abort();
    }
    reattachWorkspace(handle, workspace, dWorkspace);

    double gateGrad = 0.0;
    HANDLE_CUPP_ERROR(cupaulipropQuantumOperatorAttachCotangentBuffer(
        handle, circuit[i].op, &gateGrad, sizeof(double), CUDA_R_64F, CUPAULIPROP_MEMSPACE_HOST));
    HANDLE_CUPP_ERROR(cupaulipropPauliExpansionViewComputeOperatorApplicationBackwardDiff(
        handle, viewIn, cotOutView, cot1Expansion, circuit[i].op, 1, sortOrder, 0, 2, truncStrats, workspace, stream));
    HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansionView(viewIn));
    HANDLE_CUPP_ERROR(cupaulipropDestroyPauliExpansionView(cotOutView));

If all our circuit operators were uniquely parameterised (e.g. \(\hat{U} = \hat{U}_2(\theta_2)\hat{U}_1(\theta_1)\)), this quantity would be an element of the gradient, corresponding to e.g.

\[\frac{\partial \hat{U}}{\partial \theta_2} = \frac{\partial\hat{U}_2}{\partial \theta_2} \hat{U}_1(\theta_1).\]

But because our \(\theta_X\) and \(\theta_Z\) parameters are repeated between gates (akin to \(\hat{U} = \hat{U}_2(\theta)\hat{U}_1(\theta)\)), the product-rule prescribes

\[\frac{\partial \hat{U}}{\partial \theta} = \frac{\partial\hat{U}_2}{\partial \theta} \hat{U}_1(\theta) + \hat{U}_2(\theta) \frac{\partial\hat{U}_1}{\partial \theta},\]

and so the output quantity (which is linear in \(\hat{U}\)) is one of the above summands, contributing to either the \(\theta_X\) or \(\theta_Z\) gradient element.

    if (circuit[i].isX) gradX += gateGrad; else gradZZ += gateGrad;
    std::swap(cot0Expansion, cot1Expansion);
  }

After fully enumerating the circuit operators in this way, gradX and gradZZ contain the desired gradient.

  auto endTime = std::chrono::high_resolution_clock::now();
  auto bwdDuration = std::chrono::duration_cast<std::chrono::microseconds>(endTime - bwdStartTime);
  auto totalDuration = std::chrono::duration_cast<std::chrono::microseconds>(endTime - startTime);
  double bwdSecs = bwdDuration.count() / 1e6;
  double totalSecs = totalDuration.count() / 1e6;

  // ========================================================================
  // Results
  // ========================================================================
  std::cout << "Backward pass completed in " << bwdSecs << " seconds" << std::endl;
  std::cout << std::endl;
  std::cout << "Expectation value:       " << expec << std::endl;
  std::cout << "d<Z_62>/d(x_angle):      " << gradX << std::endl;
  std::cout << "d<Z_62>/d(zz_angle):     " << gradZZ << std::endl;
  std::cout << "Final number of terms:   " << finalTerms << std::endl;
  std::cout << "Maximum number of terms: " << maxNumTerms << std::endl;
  std::cout << "Total runtime:           " << totalSecs << " seconds" << std::endl;
  std::cout << std::endl;

Fused operators example

In this example, we demonstrate use of cuPauliProp’s fused API, to efficiently apply a sequence of operators. This interface relaxes some of the memory requirements of the single-operator API seen in the above example, and is generally higher performing, though involves more nuanced usage. This example will comparatively benchmark the fused and single-operator APIs upon an arbitrary, random Pauli expansion. Note that presently, only Clifford gates and Pauli noise channels are supported by the fused API, and showcased herein. Note this example makes use of some utility functions defined in example_utils.cpp. The full demonstration is available at fused_operators_example.cpp.

This example will make use of the below adjustable hyperparameters.

// ========================================================================
// HYPER-PARAMETERS (can be changed)
// ========================================================================


constexpr int32_t NUM_QUBITS = 700;

constexpr int64_t TARGET_NUM_TERMS = 10'000;

constexpr int32_t NUM_CLIFFORD_GATES = 300;

constexpr int32_t NUM_PAULI_NOISE_CHANNELS = 300;

Parameters NUM_QUBITS and TARGET_NUM_TERMS inform the size of the example’s utilised Pauli expansions, which may (due to their random preparation, detailed later) actually contain twice as many terms as TARGET_NUM_TERMS, but no fewer. These parameters determine the device memory used by the example.

Parameters NUM_CLIFFORD_GATES and NUM_PAULI_NOISE_CHANNELS are the respective number of Clifford gates and Pauli noise channels in the sequences that will be passed to the new fused API, and also applied in-turn via the existing single-operator API. They inform the fused API speedup, and also the required device capabilities.

For simplicity, this example will fix the below hyperparameters, which must not be changed here without other logical revisions to the example.

// ========================================================================
// CONSTANTS (are strictly fixed)
// ========================================================================


// Our Pauli expansions maintain real, double-precision coefficients
constexpr auto COEF_TYPE = CUDA_R_64F;

// We will only apply operators in the Schrodinger picture
constexpr int32_t ADJOINT = 0;

// We disable any sorting or de-duplication of pauli expansion terms, since they
// do not affect the interaction with, nor performance of, the fused-operator API
constexpr auto SORT_ORDER = CUPAULIPROP_SORT_ORDER_NONE;
constexpr int DUPLICATES = 1;

// We disable truncation for simplicity, since although accepted by the
// fused-operator API, it has only a trivial/additive performance overhead
constexpr int32_t NUM_TRUNC_STRATS = 0;
constexpr cupaulipropTruncationStrategy_t* TRUNC_STRATS = nullptr;

We also define the below helper functions, which calculate the memory needed to prepare expansions informed by the hyperparameters, or wrap cuPauliProp API calls for brevity.

// ========================================================================
// SMALL HELPER FUNCTIONS
// ========================================================================


auto getPauliStrAndCoefSizes() {

  // Size of Pauli string depends upon the number of qubits
  int32_t numPackedIntsPerXZMask;
  HANDLE_CUPP_ERROR( cupaulipropGetNumPackedIntegers(NUM_QUBITS, &numPackedIntsPerXZMask) );
  size_t numBytesPerPauliStr = 2 * numPackedIntsPerXZMask * sizeof(cupaulipropPackedIntegerType_t);

  // Size of coefficient depends upon type
  size_t numBytesPerCoef = sizeof(double); // COEF_TYPE = CUDA_R_64F

  return std::array{numBytesPerPauliStr, numBytesPerCoef};
}


auto getNewExpansionBufferSizes() {

  auto [numBytesPerPauliStr, numBytesPerCoef] = getPauliStrAndCoefSizes();
  size_t numBytesPerTerm = numBytesPerPauliStr + numBytesPerCoef;

  // Allocate enough for TWICE our desired expansion size, to handle if we over-populate
  auto numTerms = 2 * TARGET_NUM_TERMS;
  size_t pauliBufferSize = numBytesPerPauliStr * numTerms;
  size_t coefBufferSize  = numBytesPerCoef     * numTerms;
  return std::array{ pauliBufferSize, coefBufferSize };
}


auto getExistingExpansionBufferSizes(cupaulipropHandle_t handle, cupaulipropPauliExpansion_t expansion) {

  int64_t xzBitsBufferSize;
  int64_t coefBufferSize;

  void* dummyBuffer;
  int64_t dummyNumTerms;
  cupaulipropMemspace_t dummyLocatin;

  HANDLE_CUPP_ERROR( cupaulipropPauliExpansionGetStorageBuffer(
    handle, expansion, 
    &dummyBuffer, &xzBitsBufferSize, 
    &dummyBuffer, &coefBufferSize, 
    &dummyNumTerms, &dummyLocatin) );

  return std::array{xzBitsBufferSize, coefBufferSize};
}


int64_t getExistingExpansionTermCapacity(cupaulipropHandle_t handle, cupaulipropPauliExpansion_t expansion) {

  auto [xzBitsBufferSize, coefBufferSize] = getExistingExpansionBufferSizes(handle, expansion);
  auto [numBytesPerPauliStr, numBytesPerCoef] = getPauliStrAndCoefSizes();

  return std::min(
    xzBitsBufferSize / static_cast<int64_t>(numBytesPerPauliStr),
    coefBufferSize   / static_cast<int64_t>(numBytesPerCoef));
}


int64_t getPopulatedExpansionMemorySize(int64_t numTerms) {

  int32_t numPackedIntsPerXZMask;
  HANDLE_CUPP_ERROR( cupaulipropGetNumPackedIntegers(NUM_QUBITS, &numPackedIntsPerXZMask) );

  int64_t numBytesPerPauliStr =
    2LL * numPackedIntsPerXZMask * static_cast<int64_t>(sizeof(cupaulipropPackedIntegerType_t));

  int64_t numBytesPerCoef = sizeof(double); // COEF_TYPE = CUDA_R_64F
  return numTerms * (numBytesPerPauliStr + numBytesPerCoef);
}


auto createViewOfAllTerms(cupaulipropHandle_t handle, cupaulipropPauliExpansion_t expansion) {

  int64_t numTerms;
  HANDLE_CUPP_ERROR( cupaulipropPauliExpansionGetNumTerms(handle, expansion, &numTerms) );

  cupaulipropPauliExpansionView_t inView;
  HANDLE_CUPP_ERROR( cupaulipropPauliExpansionGetContiguousRange(
    handle, expansion, 0, numTerms, &inView) );

  return inView;
}

We next define functions applySingleOperator and applyFusedOperators, which will showcase the respective workflows of applying one or applying many operators.

// ========================================================================
// INDIVIDUAL AND FUSED OPERATOR APPLICATION
// ========================================================================

The applySingleOperator function uses the existing single-operator API to deterministically (i.e. without any expected failure) modify a Pauli expansion under the action of a single operator. It follows the standard paradigm witnessed in the above Kicked Ising example, whereby we prepare a view of the input expansion, call Prepare to validate our memory buffers are sufficient for the worst-case expansion growth, re-attach the workspace buffer, deterministically apply the operator, then destroy the view. The workspace and expansion buffers are prior allocated, and merely passed here.

void applySingleOperator(
  cupaulipropHandle_t handle,
  cudaStream_t stream,
  cupaulipropPauliExpansion_t inExpansion,
  cupaulipropPauliExpansion_t outExpansion,
  cupaulipropQuantumOperator_t quantumOp,
  cupaulipropWorkspaceDescriptor_t workspace,
  void* d_workspaceBuffer,
  size_t workspaceBufferSize)
{
  // Create a view of all terms in the expansion, in order for all to be processed
  auto inView = createViewOfAllTerms(handle, inExpansion);

  // Call Prepare() to obtain the REQUIRED expansion and workspace sizes; if
  // these are not satisfied, the corresponding Compute() call will error
  int64_t requiredPauliSize;
  int64_t requiredCoefSize;
  HANDLE_CUPP_ERROR( cupaulipropPauliExpansionViewPrepareOperatorApplication(
    handle, inView, quantumOp,
    SORT_ORDER, DUPLICATES, NUM_TRUNC_STRATS, TRUNC_STRATS,
    workspaceBufferSize, &requiredPauliSize, &requiredCoefSize, workspace) );

  // The required workspace size is attached to the descriptor
  int64_t requiredWorkspaceSize;
  HANDLE_CUPP_ERROR( cupaulipropWorkspaceGetMemorySize(
    handle, workspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    &requiredWorkspaceSize) );

  // Ensure sufficient memory, guaranteeing the Compute() call will succeed
  auto [pauliBufferSize, coefBufferSize] = getExistingExpansionBufferSizes(handle, outExpansion);
  if (requiredPauliSize     > static_cast<int64_t>(pauliBufferSize) ||
      requiredCoefSize      > static_cast<int64_t>(coefBufferSize)  ||
      requiredWorkspaceSize > static_cast<int64_t>(workspaceBufferSize))
  {
    std::cout
      << "Insufficient outExpansion capacity and/or workspace buffer size "
      << "to perform single-operator application. Exiting..."
      << std::endl;
    std::abort();
  }

  // Re-attach the workspace memory buffer, which is detached by the Prepare() call
  HANDLE_CUPP_ERROR( cupaulipropWorkspaceSetMemory(
    handle, workspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    d_workspaceBuffer, workspaceBufferSize) );

  // Apply the individual operator, which always exceeds (except due to a bug!)
  HANDLE_CUPP_ERROR( cupaulipropPauliExpansionViewComputeOperatorApplication(
    handle, inView, outExpansion, quantumOp,
    ADJOINT, SORT_ORDER, DUPLICATES,
    NUM_TRUNC_STRATS, TRUNC_STRATS, workspace, stream) );

  // Destroy the view, which has been invalidated by the above Compute() call
  HANDLE_CUPP_ERROR( cupaulipropDestroyPauliExpansionView(inView) );
}

Worst-case expansion growth refers to the scenario where every term in the input Pauli expansion is modified by the operator to become a sum of distinct terms, with the maximum number of summands possible of the operator. This number is sometimes referred to as the “branching factor” of the operator, and is equal to the maximum number of non-zero terms in any column of the operator’s Pauli transfer matrix. For example, the branching factor of Pauli rotation gates and amplitude damping channels is 2, while that of Pauli noise channels and Clifford gates is 1. This means applying a Pauli rotation gate via the single-operator API above requires providing an output expansion (and workspace) with sufficient capacity to store (or process) at least twice as many terms as exist in the input expansion. Such is the output of cupaulipropPauliExpansionViewPrepareOperatorApplication(), although in this demonstration, the input operators happen to have branching factors of 1.

The fused API necessarily relaxes these requirements, for reasons we elaborate upon below. Its expected usage, demonstrated in applyFusedOperators below, follows a similar pattern to applySingleOperator but has some notable nuances.

Firstly, it accepts a std::vector of quantum operators, each of which has a corresponding adjoint flag (which in this example, are all fixed to ADJOINT=0). And secondly, it involves the simultaneous management of three workspaces:

void applyFusedOperators(
  cupaulipropHandle_t handle,
  cudaStream_t stream,
  cupaulipropPauliExpansion_t inExpansion,
  cupaulipropPauliExpansion_t outExpansion,
  const std::vector<cupaulipropQuantumOperator_t>& quantumOps,
  void* d_workspaceBuffer,
  size_t workspaceBufferSize,
  cupaulipropWorkspaceDescriptor_t minWorkspace,
  cupaulipropWorkspaceDescriptor_t maxWorkspace,
  cupaulipropWorkspaceDescriptor_t avgWorkspace)
{
  // Create a view of all terms in the expansion, in order for all to be processed
  auto inView = createViewOfAllTerms(handle, inExpansion);
  const std::vector<int32_t> adjoints(quantumOps.size(), ADJOINT);

This is because, unlike the deterministic single-operator API which demanded upfront memory buffers sufficient for worst-case growth, the fused API accepts smaller buffers. If these smaller buffers turn out insufficient during computation, the API function will fail on the fly. Applying a sequence of quantum operators upon a Pauli expansion in this manner suggests three meaningful thresholds of memory requirements.

• Best case: every term in the expansion is mapped to another single term, by each operator in the sequence. The output expansion needs only capacity to store the input expansion. Ignoring the unimportant possibility of a reduction in the number of terms (via duplication, or by a non-physical operator mapping some terms to have zero coefficients), this scenario has the minimum memory requirements. Providing less memory than this minimum is invalid, and will guarantee that fused operator application fails.

• Worst case: every operator exhibits the worst-case scenario of single-operator application, whereby every term in the expansion is mapped to as many terms as the operator’s branching factor. After n operators each with branching factor b, the output Pauli expansion has a factor pow(b,n) more terms than the input expansion. When b>1 (such as for Pauli rotations and amplitude damping channels), this exponentially-growing factor is quickly intractable; 30 Pauli rotations can grow the expansion by a factor of one billion! Providing an output expansion with sufficient capacity for this scenario (along with corresponding workspace) will guarantee that the fused operator application succeeds, but is often infeasible, and gratuitously expensive when the actual operator growth is expected to be much less.

• Average case: every operator exhibits the expected growth of the expansion, when treating every input term (even those produced by each intermediate operator) as uniformly random in Hilbert space. This is not a simple interpolation between the best and worst case scenarios; the expected growth is instead informed by the Pauli transfer matrix of the fused operators. Though still often growing exponentially with the number of fused operators, the average case memory requirements are significantly smaller than the worst case requirements. For example, n Pauli rotations has an average growth factor of pow(3/2, n), while n amplitude damping channels has pow(5/4, n), despite their equal branching factors (of 2).

The fused API Prepare function, chiefly cupaulipropPauliExpansionViewPrepareOperatorFusedApplication, provides the output expansion and corresponding workspace capacities of these scenarios.

  // Call Prepare() to obtain the expansion capacities and workspaces necessary to
  // handle the best-case, worst-case and average-case Compute() scenarios. These
  // correspond to different magnitudes of expansion growth under the operators.
  int64_t minOutCapacity;
  int64_t maxOutCapacity; // -1 if overflowed
  int64_t avgOutCapacity; // -1 if overflowed
  cupaulipropStatus_t prepareStatus = cupaulipropPauliExpansionViewPrepareOperatorFusedApplication(
    handle, inView, quantumOps.size(), quantumOps.data(), adjoints.data(),
    NUM_TRUNC_STRATS, TRUNC_STRATS,
    static_cast<int64_t>(workspaceBufferSize),
    &minOutCapacity, minWorkspace,
    &avgOutCapacity, avgWorkspace,
    &maxOutCapacity, maxWorkspace);

We note that in this demonstration, only Clifford gates and Pauli noise channels will be simulated, which are both non-branching. As such, the minimum, maximum and average memories are all equivalent. We however showcase the general usage pattern anticipated for when the fused API accepts branching operators, such as Pauli rotation gates and amplitude damping channels.

The above fused-operator Prepare call can fail, returning a non-success status code, in all the ways possible of the single-operator Prepare call. However, it can also fail in some additional ways:

  // Process the additional, "recoverable" failure modes of the fused-operator API,
  // which are not exhibited by the single-operator API.
  switch (prepareStatus)
  {
    case CUPAULIPROP_STATUS_SUCCESS:
      break;
    case CUPAULIPROP_STATUS_INSUFFICIENT_DEVICE_PROPERTY:
      std::cout
        << "Too many fused operators (" << quantumOps.size() << ") were given, exceeding "
        << "a hardware limitation of the user's device. One could dynamically shrink "
        << "the number passed, but in this simple example, we treat this as an error."
        << std::endl;
      std::abort();
    default:
      HANDLE_CUPP_ERROR(prepareStatus);

    // Below are not thrown by the Prepare() call
    // case CUPAULIPROP_STATUS_INSUFFICIENT_WORKSPACE:
    // case CUPAULIPROP_STATUS_INSUFFICIENT_OUT_EXPANSION:
  }

As we will handle below, the average and maximum memory sizes can also incidentally overflow, because intractably many operators were passed. In that scenario, those outputs are set to -1.

We next process the six outputs of the Prepare call, comparing them to our given expansions and workspace, which will be allocated to a fixed size prior to calling applyFusedOperators.

  // Learn how many terms our pre-allocated outExpansion can store
  int64_t outTermCapacity = getExistingExpansionTermCapacity(handle, outExpansion);

  // Learn the workspace buffer sizes advised by Prepare() above
  int64_t minWorkspaceSize;
  int64_t avgWorkspaceSize;
  int64_t maxWorkspaceSize;
  HANDLE_CUPP_ERROR( cupaulipropWorkspaceGetMemorySize(
    handle, minWorkspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    &minWorkspaceSize) );
  HANDLE_CUPP_ERROR( cupaulipropWorkspaceGetMemorySize(
    handle, avgWorkspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    &avgWorkspaceSize) );
  HANDLE_CUPP_ERROR( cupaulipropWorkspaceGetMemorySize(
    handle, maxWorkspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    &maxWorkspaceSize) );

  // Compare our pre-allocated expansion and workspace buffers to the min-memory scenario,
  // for which the outputs cannot overflow without a returned error code. Note that we
  // already know the workspace is sufficient, since the above Prepare() call would 
  // otherwise return error CUPAULIPROP_STATUS_INSUFFICIENT_WORKSPACE.
  if (minOutCapacity > outTermCapacity || minWorkspaceSize > workspaceBufferSize)
  {
    std::cout
      << "The best-case memory sizes exceed our outExpansion and/or workspace buffer sizes. "
      << "This is invalid and will trigger an error in the Compute() call. Exiting..."
      << std::endl;
    std::abort();
  }

  // Compare our pre-allocated expansion and workspace buffers to the max-memory scenario
  if (maxOutCapacity == -1 || maxWorkspaceSize == -1)
  {
    std::cout
      << "  (The worst-case memory sizes overflowed and are reported as -1. "
      << "This example will continue, with a non-zero likelihood of failure.)"
      << std::endl;
  } 
  else if (maxOutCapacity > outTermCapacity || maxWorkspaceSize > workspaceBufferSize)
  {
    std::cout
      << "  (The worst-case memory sizes exceed our outExpansion and/or workspace buffer sizes. "
      << "This example will continue, with a non-zero likelihood of failure.)"
      << std::endl;
  }
  else
  {
    std::cout
      << "  (The worst-case memory sizes are smaller than our pre-allocated buffers. "
      << "As such, the fused API is guaranteed to succeed.)"
      << std::endl;
  }

  // Compare our pre-allocated expansion and workspace buffers to the average-memory scenario
  if (avgOutCapacity == -1 || avgWorkspaceSize == -1)
  {
    std::cout
      << "  (The average-case memory estimate overflowed and is reported as -1. "
      << "This example will continue, with a considerable likelihood of failure.)"
      << std::endl;
  } 
  else if (avgOutCapacity > outTermCapacity || avgWorkspaceSize > workspaceBufferSize)
  {
    std::cout
      << "  (The average-case memory sizes exceed our outExpansion and/or workspace buffer sizes. "
      << "This example will continue, with a possibly large likelihood of failure.)"
      << std::endl;
  } 
  else
  {
    // When the buffers exceed the average case, there's nothing significant to report!  
  }

We remind again that because this example only features non-branching gates, which cause the minimum, maximum and average memory scenarios to agree, the above pathways are not actually exercised, and our prior allocated buffers are guaranteed to be sufficient.

  // Attach our pre-allocated workspace buffer to a workspace; we arbitrarily choose avgWorkspace
  HANDLE_CUPP_ERROR( cupaulipropWorkspaceSetMemory(
    handle, avgWorkspace,
    CUPAULIPROP_MEMSPACE_DEVICE, CUPAULIPROP_WORKSPACE_SCRATCH,
    d_workspaceBuffer, workspaceBufferSize) );

We are finally ready to apply the sequence of operators, performed in a single call to cupaulipropPauliExpansionViewComputeOperatorFusedApplication. Because application can in principle fail on the fly, there are additional statuses to process.

  // Attempt to apply all operators at once
  cupaulipropStatus_t computeStatus = cupaulipropPauliExpansionViewComputeOperatorFusedApplication(
    handle, inView, outExpansion, quantumOps.size(), quantumOps.data(), adjoints.data(),
    NUM_TRUNC_STRATS, TRUNC_STRATS, avgWorkspace, stream);

  switch (computeStatus) {
    case CUPAULIPROP_STATUS_SUCCESS:
      break;
    case CUPAULIPROP_STATUS_INSUFFICIENT_WORKSPACE:
    case CUPAULIPROP_STATUS_INSUFFICIENT_OUT_EXPANSION:
      std::cout
        << "Fused application failed on-the-fly, due to the provided expansion and/or "
        << "workspace being insufficient to contain the expansion growth under the given "
        << "operators. One could dynamically shrink the number of passed operators, or "
        << "enlarge the memory capacities, but in this simple example, we opt to exit."
        << std::endl;
      std::abort();
    default:
      HANDLE_CUPP_ERROR(computeStatus);

    // Below are made impossible by our prior Prepare() call:
    // case CUPAULIPROP_FAILURE_NUM_OPERATORS_EXCEEDS_DEVICE_LIMITATION:
  }

  HANDLE_CUPP_ERROR( cupaulipropDestroyPauliExpansionView(inView) );

While we simply exit after failure to keep this example simple, it is straightforward to reduce the number of operators in the sequence and retry, re-calling the Prepare function and Compute functions, ping-ponging the input and output expansions.

That concludes the core demonstration of this example. The remaining code includes comparative benchmarking.

// ========================================================================
// BENCHMARKING
// ========================================================================

Invoking the single-operator API involves creating a simple workspace, applying each operator in-turn by ping-ponging between the input and output expansions. Beware that inExpansion below is therefore also mutated, while outExpansion will contain the final state.

auto timeSingleOperatorSequence(
  cupaulipropHandle_t handle,
  cudaStream_t stream,
  cupaulipropPauliExpansion_t inExpansion,
  cupaulipropPauliExpansion_t outExpansion,
  const std::vector<cupaulipropQuantumOperator_t>& quantumOps,
  void* d_workspaceBuffer,
  size_t workspaceBufferSize)
{
  cupaulipropWorkspaceDescriptor_t workspace;
  HANDLE_CUPP_ERROR( cupaulipropCreateWorkspaceDescriptor(handle, &workspace) );

  // Warm up before timing (leaves inExpansion unchanged)
  applySingleOperator(
    handle, stream,
    inExpansion, outExpansion, quantumOps[0],
    workspace, d_workspaceBuffer, workspaceBufferSize);
  HANDLE_CUDA_ERROR( cudaStreamSynchronize(stream) );

  const auto start = std::chrono::steady_clock::now();

  // Apply every operator in-turn
  for (auto quantumOp : quantumOps) {

    applySingleOperator(
      handle, stream,
      inExpansion, outExpansion, quantumOp,
      workspace, d_workspaceBuffer, workspaceBufferSize);

    std::swap(inExpansion, outExpansion);
  }

  HANDLE_CUDA_ERROR( cudaStreamSynchronize(stream) );

  const auto end = std::chrono::steady_clock::now();

  // Destroy workspace outside of timing
  HANDLE_CUPP_ERROR( cupaulipropDestroyWorkspaceDescriptor(workspace) );

  return end - start;
}

Invoking the fused-operator API requires preparing three workspace descriptors, as seen above. The inExpansion passed to the below benchmarking function is not mutated.

auto timeFusedOperatorSequence(
  cupaulipropHandle_t handle,
  cudaStream_t stream,
  cupaulipropPauliExpansion_t inExpansion,
  cupaulipropPauliExpansion_t outExpansion,
  const std::vector<cupaulipropQuantumOperator_t>& quantumOps,
  void* d_workspaceBuffer,
  size_t workspaceBufferSize)
{
  cupaulipropWorkspaceDescriptor_t minWorkspace;
  cupaulipropWorkspaceDescriptor_t maxWorkspace;
  cupaulipropWorkspaceDescriptor_t avgWorkspace;
  HANDLE_CUPP_ERROR( cupaulipropCreateWorkspaceDescriptor(handle, &minWorkspace) );
  HANDLE_CUPP_ERROR( cupaulipropCreateWorkspaceDescriptor(handle, &maxWorkspace) );
  HANDLE_CUPP_ERROR( cupaulipropCreateWorkspaceDescriptor(handle, &avgWorkspace) );

  // Warm up before timing (leaves inExpansion unchanged)
  applyFusedOperators(
    handle, stream,
    inExpansion, outExpansion, quantumOps,
    d_workspaceBuffer, workspaceBufferSize,
    minWorkspace, maxWorkspace, avgWorkspace);
  HANDLE_CUDA_ERROR( cudaStreamSynchronize(stream) );

  const auto start = std::chrono::steady_clock::now();

  applyFusedOperators(
    handle, stream,
    inExpansion, outExpansion, quantumOps,
    d_workspaceBuffer, workspaceBufferSize,
    minWorkspace, maxWorkspace, avgWorkspace);
  HANDLE_CUDA_ERROR( cudaStreamSynchronize(stream) );

  const auto end = std::chrono::steady_clock::now();

  HANDLE_CUPP_ERROR( cupaulipropDestroyWorkspaceDescriptor(minWorkspace) );
  HANDLE_CUPP_ERROR( cupaulipropDestroyWorkspaceDescriptor(maxWorkspace) );
  HANDLE_CUPP_ERROR( cupaulipropDestroyWorkspaceDescriptor(avgWorkspace) );

  return end - start;
}

It is straightforward to comparatively benchmark use of the single-operator and fused APIs simulating a sequence of quantum operators.

void benchmarkOperatorSequence(
  cupaulipropHandle_t handle,
  cudaStream_t stream,
  cupaulipropPauliExpansion_t inExpansion,
  cupaulipropPauliExpansion_t outExpansion,
  const std::vector<cupaulipropQuantumOperator_t>& quantumOps,
  void* d_workspaceBuffer, size_t workspaceBufferSize)
{
  // Time the fused API first, since it does not modify inExpansion
  auto fusedDuration = timeFusedOperatorSequence(
    handle, stream,
    inExpansion, outExpansion, quantumOps,
    d_workspaceBuffer, workspaceBufferSize);
  
  // Time the single-operator API second, because it involves ping-ponging
  // between inExpansion and outExpansion, mutating both! 
  auto sequentialDuration = timeSingleOperatorSequence(
    handle, stream,
    inExpansion, outExpansion, quantumOps,
    d_workspaceBuffer, workspaceBufferSize);

  // Report the measured speedup of the fused API over the single-operator API.
  const auto fusedMs      = std::chrono::duration<double, std::milli>(fusedDuration     ).count();
  const auto sequentialMs = std::chrono::duration<double, std::milli>(sequentialDuration).count();
  std::cout << std::endl;
  std::cout << "Single-operator API: " << sequentialMs << " ms" << std::endl;
  std::cout << "Fused-operator API:  " << fusedMs      << " ms" << std::endl;
  std::cout << "Speedup:             " << sequentialMs / fusedMs << "x\n" << std::endl;
}

Our final chore is to invoke the above benchmarking from main.

// ========================================================================
// MAIN
// ========================================================================


int main()
{
  std::cout << "cuPauliProp fused operators example" << std::endl;
  std::cout << "===================================" << std::endl << std::endl;

This involves preparing the cuPauliProp library, creating the input and output expansions and their corresponding memory buffers, and initialising the input expansion; to an arbitrary single X operator on qubit 0.

  /*
   * PREPARE LIBRARY AND MEMORY BUFFERS
   */

  HANDLE_CUDA_ERROR( cudaSetDevice(0) );
  cudaStream_t stream = 0;

  cupaulipropHandle_t handle{};
  HANDLE_CUPP_ERROR( cupaulipropCreate(&handle) );

  // Prepare buffers for two Pauli expansions
  const auto [pauliBufferSize, coefBufferSize] = getNewExpansionBufferSizes();
  void* d_inPauliBuffer;
  void* d_inCoefBuffer;
  void* d_outPauliBuffer;
  void* d_outCoefBuffer;
  HANDLE_CUDA_ERROR( cudaMalloc(&d_inPauliBuffer,   pauliBufferSize) );
  HANDLE_CUDA_ERROR( cudaMalloc(&d_inCoefBuffer,    coefBufferSize ) );
  HANDLE_CUDA_ERROR( cudaMalloc(&d_outPauliBuffer,  pauliBufferSize) );
  HANDLE_CUDA_ERROR( cudaMalloc(&d_outCoefBuffer,   coefBufferSize ) );

  // Initialise the "in" Pauli expansion to a single term; X on qubit 0
  const cupaulipropPackedIntegerType_t initX = 1;
  const double initCoef = 1.0;
  HANDLE_CUDA_ERROR( cudaMemset(d_inPauliBuffer, 0, pauliBufferSize) );
  HANDLE_CUDA_ERROR( cudaMemcpy(d_inPauliBuffer, &initX,    sizeof(initX),    cudaMemcpyHostToDevice) );
  HANDLE_CUDA_ERROR( cudaMemcpy(d_inCoefBuffer,  &initCoef, sizeof(initCoef), cudaMemcpyHostToDevice) );

  // Create two Pauli expansions, to ping-pong as input and output to the cuPP API
  cupaulipropPauliExpansion_t inExpansion;
  cupaulipropPauliExpansion_t outExpansion;
  HANDLE_CUPP_ERROR( cupaulipropCreatePauliExpansion(
    handle, NUM_QUBITS,
    d_inPauliBuffer, pauliBufferSize,
    d_inCoefBuffer,  coefBufferSize, COEF_TYPE,
    /*numTerms=*/1,
    SORT_ORDER, DUPLICATES,
    &inExpansion) );
  HANDLE_CUPP_ERROR( cupaulipropCreatePauliExpansion(
    handle, NUM_QUBITS,
    d_outPauliBuffer, pauliBufferSize,
    d_outCoefBuffer,  coefBufferSize, COEF_TYPE,
    /*numTerms=*/0,
    SORT_ORDER, DUPLICATES,
    &outExpansion) );

  // Allocate workspace buffer; descriptors are created where needed.
  void* d_workspaceBuffer;
  size_t workspaceBufferSize = 2 * (pauliBufferSize + coefBufferSize);
  HANDLE_CUDA_ERROR( cudaMalloc(&d_workspaceBuffer, workspaceBufferSize) ); // attached in Prepare()

Because this example will only showcase non-branching operators (Clifford gates and Pauli noise channels), the input expansion will not grow from the current single term by application of the operators. It is therefore necessary to prepare a larger expansion; as large as the hyperparameter TARGET_NUM_TERMS. We do this by making use of the growExpansionViaRandomRotations() function defined in example_utils.cpp, populating inExpansion with numInputTerms-many arbitrary terms.

  /*
   * PREPARE A RANDOM EXPANSION
   */

  // Grow inExpansion (not outExpansion) to size >= TARGET_NUM_TERMS via random Pauli rotations
  growExpansionViaRandomRotations(
    handle, stream,
    inExpansion, outExpansion,
    NUM_QUBITS, TARGET_NUM_TERMS,
    pauliBufferSize, coefBufferSize,
    d_workspaceBuffer, workspaceBufferSize);

  int64_t numInputTerms;
  HANDLE_CUPP_ERROR( cupaulipropPauliExpansionGetNumTerms(handle, inExpansion, &numInputTerms) );
  const int64_t inputExpansionBytes = getPopulatedExpansionMemorySize(numInputTerms);
  const double inputExpansionMiB = static_cast<double>(inputExpansionBytes) / (1024.0 * 1024.0);
  std::cout << "Input expansion size: " << numInputTerms
            << " terms (" << NUM_QUBITS << " qubits)" << std::endl;
  std::cout << "Input expansion memory: " << inputExpansionBytes
            << " bytes (" << std::fixed << std::setprecision(2) << inputExpansionMiB
            << " MiB)" << std::endl << std::endl;

We first comparatively benchmark the single-operator and fused APIs with sequences of Pauli noise channels. We use an arbitrary list of one and two-qubit channels, with random target qubits and error probabilities, obtained from createRandomPauliNoiseChannel(), defined in example_utils.cpp.

  /*
   * BENCHMARK RANDOM PAULI NOISE CHANNELS
   */

  std::cout << "Applying " << NUM_PAULI_NOISE_CHANNELS << " Pauli noise channels...\n" << std::endl;

  // Prepare a list of random Pauli channels (1 and 2 qubit)
  std::mt19937_64 rng(5678);
  std::vector<cupaulipropQuantumOperator_t> pauliNoiseChannels;
  for (int32_t i = 0; i < NUM_PAULI_NOISE_CHANNELS; ++i)
    pauliNoiseChannels.push_back(createRandomPauliNoiseChannel(handle, rng, NUM_QUBITS));

  benchmarkOperatorSequence(
    handle, stream,
    inExpansion, outExpansion, pauliNoiseChannels,
    d_workspaceBuffer, workspaceBufferSize);

Note that Pauli noise channels do not mutate Pauli strings, and instead only change the expansion term coefficients, which themselves do not affect the runtime of subsequent operations performed in this example. As such, there is no need to restore the original state in inExpansion before proceeding.

We then comparatively benchmark simulation of random Clifford gates, again using utilities defined in example_utils.cpp.

  /*
   * BENCHMARK RANDOM CLIFFORD GATES
   */

   std::cout << "Applying " << NUM_CLIFFORD_GATES << " random Clifford gates...\n" << std::endl;

  // Prepare a list of random Clifford gates (1 and 2 qubit, and random kinds)
  std::vector<cupaulipropQuantumOperator_t> cliffordGates;
  for (int32_t i = 0; i < NUM_CLIFFORD_GATES; ++i)
    cliffordGates.push_back(createRandomCliffordGate(handle, rng, NUM_QUBITS));

  // Comparative benchmarking penalises fused API with print statements - no big deal!
  benchmarkOperatorSequence(
    handle, stream,
    inExpansion, outExpansion, cliffordGates,
    d_workspaceBuffer, workspaceBufferSize);

With benchmarking complete, our final chore is to clean up.

  /*
   * CLEAN UP
   */
  
  for (auto gate : cliffordGates)
    HANDLE_CUPP_ERROR( cupaulipropDestroyOperator(gate) );
  for (auto channel : pauliNoiseChannels)
    HANDLE_CUPP_ERROR( cupaulipropDestroyOperator(channel) );
  HANDLE_CUPP_ERROR( cupaulipropDestroyPauliExpansion(inExpansion) );
  HANDLE_CUPP_ERROR( cupaulipropDestroyPauliExpansion(outExpansion) );
  HANDLE_CUPP_ERROR( cupaulipropDestroy(handle) );
  HANDLE_CUDA_ERROR( cudaFree(d_workspaceBuffer) );
  HANDLE_CUDA_ERROR( cudaFree(d_inPauliBuffer) );
  HANDLE_CUDA_ERROR( cudaFree(d_inCoefBuffer) );
  HANDLE_CUDA_ERROR( cudaFree(d_outPauliBuffer) );
  HANDLE_CUDA_ERROR( cudaFree(d_outCoefBuffer) );

  return EXIT_SUCCESS;
}

With the hyperparameters shown in the first code section, and using an NVIDIA A100 GPU, we anticipate an output such as:

cuPauliProp fused operators example
===================================

Input expansion size: 11402 terms (700 qubits)
Input expansion memory: 2097968 bytes (2.00 MiB)

Applying 300 Pauli noise channels...

  (The worst-case memory sizes are smaller than our pre-allocated buffers. As such, the fused API is guaranteed to succeed.)
  (The worst-case memory sizes are smaller than our pre-allocated buffers. As such, the fused API is guaranteed to succeed.)

Single-operator API: 20.09 ms
Fused-operator API:  1.21 ms
Speedup:             16.57x

Applying 300 random Clifford gates...

  (The worst-case memory sizes are smaller than our pre-allocated buffers. As such, the fused API is guaranteed to succeed.)
  (The worst-case memory sizes are smaller than our pre-allocated buffers. As such, the fused API is guaranteed to succeed.)

Single-operator API: 17.33 ms
Fused-operator API:  1.04 ms
Speedup:             16.60x

Useful tips

• For debugging, one can set the environment variable CUPAULIPROP_LOG_LEVEL=n. The level n = 0, 1, …, 5 corresponds to the logger level as described in the table below. The environment variable CUPAULIPROP_LOG_FILE=<filepath> can be used to redirect the log output to a custom file at <filepath> instead of stdout.

Level	Summary	Long Description
0	Off	Logging is disabled (default)
1	Errors	Only errors will be logged
2	Performance Trace	API calls that launch CUDA kernels will log their parameters and important information
3	Performance Hints	Hints that can potentially improve the application’s performance
4	Heuristics Trace	Provides general information about the library execution, may contain details about heuristic status
5	API Trace	API calls will log their parameter and important information