矩阵乘积态采样 (QFT 电路)

以下代码示例说明了如何为给定的量子线路 (QFT) 定义张量网络状态，然后计算其矩阵乘积态 (MPS) 分解，最后对 MPS 分解状态进行采样。完整的代码可以在 NVIDIA/cuQuantum 仓库中找到 (此处)。

头文件和错误处理

#include <cstdlib>
#include <cstdio>
#include <cassert>
#include <complex>
#include <cmath>
#include <vector>
#include <iostream>

#include <cuda_runtime.h>
#include <cutensornet.h>

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

#define HANDLE_CUTN_ERROR(x) \
{ const auto err = x; \
  if( err != CUTENSORNET_STATUS_SUCCESS ) \
  { printf("cuTensorNet error %s in line %d\n", cutensornetGetErrorString(err), __LINE__); fflush(stdout); std::abort(); } \
};


int main()
{
  static_assert(sizeof(size_t) == sizeof(int64_t), "Please build this sample on a 64-bit architecture!");

  constexpr std::size_t fp64size = sizeof(double);

定义张量网络状态和要生成的所需输出样本数

让我们为具有给定量子比特数的量子线路定义一个张量网络状态，并请求为完整的量子比特寄存器生成给定数量的输出样本。

  // Quantum state configuration
  const int64_t numSamples = 128; // number of samples to generate
  const int32_t numQubits = 16;   // number of qubits
  const std::vector<int64_t> qubitDims(numQubits, 2); // qubit dimensions
  std::cout << "Quantum circuit: " << numQubits << " qubits; " << numSamples << " samples\n";

初始化 cuTensorNet 库句柄

  // Initialize the cuTensorNet library
  HANDLE_CUDA_ERROR(cudaSetDevice(0));
  cutensornetHandle_t cutnHandle;
  HANDLE_CUTN_ERROR(cutensornetCreate(&cutnHandle));
  std::cout << "Initialized cuTensorNet library on GPU 0\n";

在 GPU 内存中定义量子门

  // Define necessary quantum gate tensors in Host memory
  const double invsq2 = 1.0 / std::sqrt(2.0);
  const double pi = 3.14159265358979323846;
  using GateData = std::vector<std::complex<double>>;
  //  CR(k) gate generator
  auto cRGate = [&pi] (int32_t k) {
    const double phi = pi * 2.0 / std::exp2(k);
    const GateData cr {{1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
                       {0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
                       {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0},
                       {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {std::cos(phi), std::sin(phi)}};
    return cr;
  };
  //  Hadamard gate
  const GateData h_gateH {{invsq2, 0.0}, {invsq2, 0.0},
                          {invsq2, 0.0}, {-invsq2, 0.0}};
  //  CR(k) gates (controlled rotations)
  std::vector<GateData> h_gateCR(numQubits);
  for(int32_t k = 0; k < numQubits; ++k) {
    h_gateCR[k] = cRGate(k+1);
  }

  // Copy quantum gates into Device memory
  void *d_gateH {nullptr};
  std::vector<void*> d_gateCR(numQubits, nullptr);
  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateH, 4 * (2 * fp64size)));
  for(int32_t k = 0; k < numQubits; ++k) {
    HANDLE_CUDA_ERROR(cudaMalloc(&(d_gateCR[k]), 16 * (2 * fp64size)));
  }
  std::cout << "Allocated GPU memory for quantum gates\n";
  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateH, h_gateH.data(), 4 * (2 * fp64size), cudaMemcpyHostToDevice));
  for(int32_t k = 0; k < numQubits; ++k) {
    HANDLE_CUDA_ERROR(cudaMemcpy(d_gateCR[k], h_gateCR[k].data(), 16 * (2 * fp64size), cudaMemcpyHostToDevice));
  }
  std::cout << "Copied quantum gates into GPU memory\n";

在 GPU 内存中分配 MPS 张量

在这里，我们设置 MPS 张量的形状并为其存储分配 GPU 内存。

  // Define the MPS representation and allocate memory buffers for the MPS tensors
  const int64_t maxExtent = 8; // max bond dimension (level of low-rank MPS approximation)
  std::vector<std::vector<int64_t>> extents;
  std::vector<int64_t*> extentsPtr(numQubits);
  std::vector<void*> d_mpsTensors(numQubits, nullptr);
  for (int32_t i = 0; i < numQubits; ++i) {
    if (i == 0) { // left boundary MPS tensor
      extents.push_back({2, maxExtent});
      HANDLE_CUDA_ERROR(cudaMalloc(&d_mpsTensors[i], 2 * maxExtent * (2 * fp64size)));
    } 
    else if (i == numQubits-1) { // right boundary MPS tensor
      extents.push_back({maxExtent, 2});
      HANDLE_CUDA_ERROR(cudaMalloc(&d_mpsTensors[i], 2 * maxExtent * (2 * fp64size)));
    }
    else { // middle MPS tensors
      extents.push_back({maxExtent, 2, maxExtent});
      HANDLE_CUDA_ERROR(cudaMalloc(&d_mpsTensors[i], 2 * maxExtent * maxExtent * (2 * fp64size)));
    }
    extentsPtr[i] = extents[i].data();
  }
  std::cout << "Allocated GPU memory for MPS tensors\n";

在 GPU 上分配暂存缓冲区

  // Query free memory on Device and allocate a scratch buffer
  std::size_t freeSize {0}, totalSize {0};
  HANDLE_CUDA_ERROR(cudaMemGetInfo(&freeSize, &totalSize));
  const std::size_t scratchSize = (freeSize - (freeSize % 4096)) / 2; // use half of available memory with alignment
  void *d_scratch {nullptr};
  HANDLE_CUDA_ERROR(cudaMalloc(&d_scratch, scratchSize));
  std::cout << "Allocated " << scratchSize << " bytes of scratch memory on GPU: "
            << "[" << d_scratch << ":" << (void*)(((char*)(d_scratch))  + scratchSize) << ")\n";

创建纯张量网络状态

现在让我们为具有给定量子比特数的量子线路创建一个纯张量网络状态。

  // Create the initial quantum state
  cutensornetState_t quantumState;
  HANDLE_CUTN_ERROR(cutensornetCreateState(cutnHandle, CUTENSORNET_STATE_PURITY_PURE, numQubits, qubitDims.data(),
                    CUDA_C_64F, &quantumState));
  std::cout << "Created the initial (vacuum) quantum state\n";

应用量子门

让我们构建没有位反转的 QFT 量子线路，方法是应用相应的量子门。

  // Construct the QFT quantum circuit state (apply quantum gates)
  //  Example of a QFT circuit with 3 qubits (with no bit reversal):
  //  Q0--H--CR2--CR3-------------
  //         |    |
  //  Q1-----o----|----H--CR2-----
  //              |       |
  //  Q2----------o-------o----H--
  int64_t id;
  for(int32_t i = 0; i < numQubits; ++i) {
    HANDLE_CUTN_ERROR(cutensornetStateApplyTensorOperator(cutnHandle, quantumState, 1, std::vector<int32_t>{{i}}.data(),
                      d_gateH, nullptr, 1, 0, 1, &id));
    for(int32_t j = (i+1); j < numQubits; ++j) {
      HANDLE_CUTN_ERROR(cutensornetStateApplyTensorOperator(cutnHandle, quantumState, 2, std::vector<int32_t>{{j,i}}.data(),
                        d_gateCR[j-i], nullptr, 1, 0, 1, &id));
    }
  }
  std::cout << "Applied quantum gates\n";

请求最终量子线路状态的 MPS 分解

在这里，我们表达了使用 MPS 分解来分解最终量子线路状态的意图。提供的 MPS 张量的形状（模式范围）是指 MPS 重整化过程中它们的最大尺寸限制。最终 MPS 张量的实际计算形状（模式范围）可能小于其限制。请注意，此处尚未进行任何计算。

  // Specify the target MPS representation (use default strides)
  HANDLE_CUTN_ERROR(cutensornetStateFinalizeMPS(cutnHandle, quantumState, 
                    CUTENSORNET_BOUNDARY_CONDITION_OPEN, extentsPtr.data(), /*strides=*/nullptr ));
  std::cout << "Finalized the form of the MPS representation\n";

配置 MPS 分解过程

在表达了执行最终量子线路状态的 MPS 分解的意图之后，我们还可以通过重置不同的选项（例如，SVD 算法）来配置 MPS 分解过程。

  // Set up the SVD method for bonds truncation (optional)
  cutensornetTensorSVDAlgo_t algo = CUTENSORNET_TENSOR_SVD_ALGO_GESVDJ; 
  HANDLE_CUTN_ERROR(cutensornetStateConfigure(cutnHandle, quantumState, 
                    CUTENSORNET_STATE_MPS_SVD_CONFIG_ALGO, &algo, sizeof(algo)));
  std::cout << "Configured the MPS computation\n";

准备 MPS 分解的计算

让我们创建一个工作区描述符并准备最终量子线路状态的 MPS 分解的计算。

  // Prepare the MPS computation and attach a workspace buffer
  cutensornetWorkspaceDescriptor_t workDesc;
  HANDLE_CUTN_ERROR(cutensornetCreateWorkspaceDescriptor(cutnHandle, &workDesc));
  std::cout << "Created the workspace descriptor\n";
  HANDLE_CUTN_ERROR(cutensornetStatePrepare(cutnHandle, quantumState, scratchSize, workDesc, 0x0));
  int64_t worksize {0};
  HANDLE_CUTN_ERROR(cutensornetWorkspaceGetMemorySize(cutnHandle,
                                                      workDesc,
                                                      CUTENSORNET_WORKSIZE_PREF_RECOMMENDED,
                                                      CUTENSORNET_MEMSPACE_DEVICE,
                                                      CUTENSORNET_WORKSPACE_SCRATCH,
                                                      &worksize));
  std::cout << "Scratch GPU workspace size (bytes) for MPS computation = " << worksize << std::endl;
  if(worksize <= scratchSize) {
    HANDLE_CUTN_ERROR(cutensornetWorkspaceSetMemory(cutnHandle, workDesc, CUTENSORNET_MEMSPACE_DEVICE,
                      CUTENSORNET_WORKSPACE_SCRATCH, d_scratch, worksize));
  }else{
    std::cout << "ERROR: Insufficient workspace size on Device!\n";
    std::abort();
  }
  std::cout << "Set the workspace buffer for MPS computation\n";

计算 MPS 分解

一旦 MPS 分解过程已配置和准备好，让我们计算最终量子线路状态的 MPS 分解。

  // Compute the MPS factorization of the quantum circuit state
  HANDLE_CUTN_ERROR(cutensornetStateCompute(cutnHandle, quantumState, 
                    workDesc, extentsPtr.data(), /*strides=*/nullptr, d_mpsTensors.data(), 0));
  std::cout << "Computed the MPS factorization of the quantum circuit state\n";
  

创建张量网络状态采样器

一旦量子线路状态已经构建并使用 MPS 表示进行分解，让我们为完整的量子比特寄存器（所有量子比特）创建张量网络状态采样器。

  // Create the quantum circuit sampler
  cutensornetStateSampler_t sampler;
  HANDLE_CUTN_ERROR(cutensornetCreateSampler(cutnHandle, quantumState, numQubits, nullptr, &sampler));
  std::cout << "Created the quantum circuit sampler\n";

配置张量网络状态采样器

可选地，我们可以通过设置张量网络收缩路径查找器要使用的超样本数量来配置张量网络状态采样器。

  // Configure the quantum circuit sampler
  const int32_t numHyperSamples = 8; // desired number of hyper samples used in the tensor network contraction path finder
  HANDLE_CUTN_ERROR(cutensornetSamplerConfigure(cutnHandle, sampler,
                    CUTENSORNET_SAMPLER_OPT_NUM_HYPER_SAMPLES, &numHyperSamples, sizeof(numHyperSamples)));

准备张量网络状态采样器

现在让我们准备张量网络状态采样器。

  // Prepare the quantum circuit sampler
  HANDLE_CUTN_ERROR(cutensornetSamplerPrepare(cutnHandle, sampler, scratchSize, workDesc, 0x0));
  std::cout << "Prepared the quantum circuit state sampler\n";

设置工作区

现在我们可以设置所需的工作区缓冲区。

  // Attach the workspace buffer
  HANDLE_CUTN_ERROR(cutensornetWorkspaceGetMemorySize(cutnHandle,
                                                      workDesc,
                                                      CUTENSORNET_WORKSIZE_PREF_RECOMMENDED,
                                                      CUTENSORNET_MEMSPACE_DEVICE,
                                                      CUTENSORNET_WORKSPACE_SCRATCH,
                                                      &worksize));
  std::cout << "Scratch GPU workspace size (bytes) for MPS Sampling = " << worksize << std::endl;
  assert(worksize > 0);
  if(worksize <= scratchSize) {
    HANDLE_CUTN_ERROR(cutensornetWorkspaceSetMemory(cutnHandle, workDesc, CUTENSORNET_MEMSPACE_DEVICE,
                      CUTENSORNET_WORKSPACE_SCRATCH, d_scratch, worksize));
  }else{
    std::cout << "ERROR: Insufficient workspace size on Device!\n";
    std::abort();
  }
  std::cout << "Set the workspace buffer\n";

执行最终量子线路状态的采样

一旦一切都设置好，我们执行量子线路状态的采样并打印输出样本。

  // Sample the quantum circuit state
  std::vector<int64_t> samples(numQubits * numSamples); // samples[SampleId][QubitId] reside in Host memory
  HANDLE_CUTN_ERROR(cutensornetSamplerSample(cutnHandle, sampler, numSamples, workDesc, samples.data(), 0));
  std::cout << "Performed quantum circuit state sampling\n";
  std::cout << "Bit-string samples:\n";
  for(int64_t i = 0; i < numSamples; ++i) {
    for(int64_t j = 0; j < numQubits; ++j) std::cout << " " << samples[i * numQubits + j];
    std::cout << std::endl;
  }

释放资源

  // Destroy the workspace descriptor
  HANDLE_CUTN_ERROR(cutensornetDestroyWorkspaceDescriptor(workDesc));
  std::cout << "Destroyed the workspace descriptor\n";

  // Destroy the quantum circuit sampler
  HANDLE_CUTN_ERROR(cutensornetDestroySampler(sampler));
  std::cout << "Destroyed the quantum circuit state sampler\n";

  // Destroy the quantum circuit state
  HANDLE_CUTN_ERROR(cutensornetDestroyState(quantumState));
  std::cout << "Destroyed the quantum circuit state\n";

  // Free GPU buffers
  for (int32_t i = 0; i < numQubits; i++) {
    HANDLE_CUDA_ERROR(cudaFree(d_mpsTensors[i]));
  }
  HANDLE_CUDA_ERROR(cudaFree(d_scratch));
  for(auto ptr: d_gateCR) HANDLE_CUDA_ERROR(cudaFree(ptr));
  HANDLE_CUDA_ERROR(cudaFree(d_gateH));
  std::cout << "Freed memory on GPU\n";

  // Finalize the cuTensorNet library
  HANDLE_CUTN_ERROR(cutensornetDestroy(cutnHandle));
  std::cout << "Finalized the cuTensorNet library\n";

  return 0;
}