Computing Tensor Network State Expectation Value

The following code example illustrates how to define a tensor network state, a tensor network operator, and then compute the expectation value of the tensor network operator over the tensor network state. The full code can be found in the NVIDIA/cuQuantum repository (here).

Headers and error handling

 7#include <cstdlib>
 8#include <cstdio>
 9#include <cassert>
10#include <complex>
11#include <vector>
12#include <bitset>
13#include <iostream>
14
15#include <cuda_runtime.h>
16#include <cutensornet.h>
17
18
19#define HANDLE_CUDA_ERROR(x) \
20{ const auto err = x; \
21  if( err != cudaSuccess ) \
22  { printf("CUDA error %s in line %d\n", cudaGetErrorString(err), __LINE__); fflush(stdout); std::abort(); } \
23};
24
25#define HANDLE_CUTN_ERROR(x) \
26{ const auto err = x; \
27  if( err != CUTENSORNET_STATUS_SUCCESS ) \
28  { printf("cuTensorNet error %s in line %d\n", cutensornetGetErrorString(err), __LINE__); fflush(stdout); std::abort(); } \
29};
30
31
32int main(int argc, char **argv)
33{
34  static_assert(sizeof(size_t) == sizeof(int64_t), "Please build this sample on a 64-bit architecture!");
35
36  constexpr std::size_t fp64size = sizeof(double);

Define the tensor network state

Let’s define a tensor network state corresponding to a 16-qubit quantum circuit.

40  // Quantum state configuration
41  constexpr int32_t numQubits = 16; // number of qubits
42  const std::vector<int64_t> qubitDims(numQubits,2); // qubit dimensions
43  std::cout << "Quantum circuit: " << numQubits << " qubits\n";

Initialize the cuTensorNet library handle

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

Define quantum gates on GPU

55  // Define necessary quantum gate tensors in Host memory
56  const double invsq2 = 1.0 / std::sqrt(2.0);
57  //  Hadamard gate
58  const std::vector<std::complex<double>> h_gateH {{invsq2, 0.0},  {invsq2, 0.0},
59                                                   {invsq2, 0.0}, {-invsq2, 0.0}};
60  //  Pauli X gate
61  const std::vector<std::complex<double>> h_gateX {{0.0, 0.0}, {1.0, 0.0},
62                                                   {1.0, 0.0}, {0.0, 0.0}};
63  //  Pauli Y gate
64  const std::vector<std::complex<double>> h_gateY {{0.0, 0.0}, {0.0, -1.0},
65                                                   {0.0, 1.0}, {0.0, 0.0}};
66  //  Pauli Z gate
67  const std::vector<std::complex<double>> h_gateZ {{1.0, 0.0}, {0.0, 0.0},
68                                                   {0.0, 0.0}, {-1.0, 0.0}};
69  //  CX gate
70  const std::vector<std::complex<double>> h_gateCX {{1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
71                                                    {0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
72                                                    {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0},
73                                                    {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}};
74
75  // Copy quantum gates to Device memory
76  void *d_gateH{nullptr}, *d_gateX{nullptr}, *d_gateY{nullptr}, *d_gateZ{nullptr}, *d_gateCX{nullptr};
77  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateH, 4 * (2 * fp64size)));
78  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateX, 4 * (2 * fp64size)));
79  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateY, 4 * (2 * fp64size)));
80  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateZ, 4 * (2 * fp64size)));
81  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateCX, 16 * (2 * fp64size)));
82  std::cout << "Allocated quantum gate memory on GPU\n";
83  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateH, h_gateH.data(), 4 * (2 * fp64size), cudaMemcpyHostToDevice));
84  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateX, h_gateX.data(), 4 * (2 * fp64size), cudaMemcpyHostToDevice));
85  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateY, h_gateY.data(), 4 * (2 * fp64size), cudaMemcpyHostToDevice));
86  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateZ, h_gateZ.data(), 4 * (2 * fp64size), cudaMemcpyHostToDevice));
87  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateCX, h_gateCX.data(), 16 * (2 * fp64size), cudaMemcpyHostToDevice));
88  std::cout << "Copied quantum gates to GPU memory\n";

Allocate the scratch buffer on GPU

92  // Query the free memory on Device
93  std::size_t freeSize{0}, totalSize{0};
94  HANDLE_CUDA_ERROR(cudaMemGetInfo(&freeSize, &totalSize));
95  const std::size_t scratchSize = (freeSize - (freeSize % 4096)) / 2; // use half of available memory with alignment
96  void *d_scratch{nullptr};
97  HANDLE_CUDA_ERROR(cudaMalloc(&d_scratch, scratchSize));
98  std::cout << "Allocated " << scratchSize << " bytes of scratch memory on GPU\n";

Create a pure tensor network state

Now let’s create a pure tensor network state for a 16-qubit quantum circuit.

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

Apply quantum gates

Let’s construct the GHZ quantum circuit by applying the corresponding quantum gates.

110  // Construct the final quantum circuit state (apply quantum gates) for the GHZ circuit
111  int64_t id;
112  HANDLE_CUTN_ERROR(cutensornetStateApplyTensor(cutnHandle, quantumState, 1, std::vector<int32_t>{{0}}.data(),
113                    d_gateH, nullptr, 1, 0, 1, &id));
114  for(int32_t i = 1; i < numQubits; ++i) {
115    HANDLE_CUTN_ERROR(cutensornetStateApplyTensor(cutnHandle, quantumState, 2, std::vector<int32_t>{{i-1,i}}.data(),
116                      d_gateCX, nullptr, 1, 0, 1, &id));
117  }
118  std::cout << "Applied quantum gates\n";

Construct a tensor network operator

Now let’s create an empty tensor network operator for 16-qubit states and then append three components to it, where each component is a direct product of Pauli matrices scaled by some complex coefficient (like in the Jordan-Wigner representation).

122  // Create an empty tensor network operator
123  cutensornetNetworkOperator_t hamiltonian;
124  HANDLE_CUTN_ERROR(cutensornetCreateNetworkOperator(cutnHandle, numQubits, qubitDims.data(), CUDA_C_64F, &hamiltonian));
125  // Append component (0.5 * Z1 * Z2) to the tensor network operator
126  {
127    const int32_t numModes[] = {1, 1}; // Z1 acts on 1 mode, Z2 acts on 1 mode
128    const int32_t modesZ1[] = {1}; // state modes Z1 acts on
129    const int32_t modesZ2[] = {2}; // state modes Z2 acts on
130    const int32_t * stateModes[] = {modesZ1, modesZ2}; // state modes (Z1 * Z2) acts on
131    const void * gateData[] = {d_gateZ, d_gateZ}; // GPU pointers to gate data
132    HANDLE_CUTN_ERROR(cutensornetNetworkOperatorAppendProduct(cutnHandle, hamiltonian, cuDoubleComplex{0.5,0.0},
133                      2, numModes, stateModes, NULL, gateData, &id));
134  }
135  // Append component (0.25 * Y3) to the tensor network operator
136  {
137    const int32_t numModes[] = {1}; // Y3 acts on 1 mode
138    const int32_t modesY3[] = {3}; // state modes Y3 acts on
139    const int32_t * stateModes[] = {modesY3}; // state modes (Y3) acts on
140    const void * gateData[] = {d_gateY}; // GPU pointers to gate data
141    HANDLE_CUTN_ERROR(cutensornetNetworkOperatorAppendProduct(cutnHandle, hamiltonian, cuDoubleComplex{0.25,0.0},
142                      1, numModes, stateModes, NULL, gateData, &id));
143  }
144  // Append component (0.13 * Y0 X2 Z3) to the tensor network operator
145  {
146    const int32_t numModes[] = {1, 1, 1}; // Y0 acts on 1 mode, X2 acts on 1 mode, Z3 acts on 1 mode
147    const int32_t modesY0[] = {0}; // state modes Y0 acts on
148    const int32_t modesX2[] = {2}; // state modes X2 acts on
149    const int32_t modesZ3[] = {3}; // state modes Z3 acts on
150    const int32_t * stateModes[] = {modesY0, modesX2, modesZ3}; // state modes (Y0 * X2 * Z3) acts on
151    const void * gateData[] = {d_gateY, d_gateX, d_gateZ}; // GPU pointers to gate data
152    HANDLE_CUTN_ERROR(cutensornetNetworkOperatorAppendProduct(cutnHandle, hamiltonian, cuDoubleComplex{0.13,0.0},
153                      3, numModes, stateModes, NULL, gateData, &id));
154  }
155  std::cout << "Constructed a tensor network operator: (0.5 * Z1 * Z2) + (0.25 * Y3) + (0.13 * Y0 * X2 * Z3)" << std::endl;

Create the expectation value object

Once the quantum circuit and the tensor network operator have been constructed, let’s create the expectation value object that will compute the expectation value of the specified tensor network operator over the final state of the specified quantum circuit.

159  // Specify the quantum circuit expectation value
160  cutensornetStateExpectation_t expectation;
161  HANDLE_CUTN_ERROR(cutensornetCreateExpectation(cutnHandle, quantumState, hamiltonian, &expectation));
162  std::cout << "Created the specified quantum circuit expectation value\n";

Configure the expectation value calculation

Now we can configure the expectation value object by setting the number of hyper-samples to be used by the tensor network contraction path finder.

166  // Configure the computation of the specified quantum circuit expectation value
167  const int32_t numHyperSamples = 8; // desired number of hyper samples used in the tensor network contraction path finder
168  HANDLE_CUTN_ERROR(cutensornetExpectationConfigure(cutnHandle, expectation,
169                    CUTENSORNET_EXPECTATION_OPT_NUM_HYPER_SAMPLES, &numHyperSamples, sizeof(numHyperSamples)));

Prepare the expectation value calculation

Let’s create a workspace descriptor and prepare the computation of the desired expectation value.

173  // Prepare the specified quantum circuit expectation value for computation
174  cutensornetWorkspaceDescriptor_t workDesc;
175  HANDLE_CUTN_ERROR(cutensornetCreateWorkspaceDescriptor(cutnHandle, &workDesc));
176  std::cout << "Created the workspace descriptor\n";
177  HANDLE_CUTN_ERROR(cutensornetExpectationPrepare(cutnHandle, expectation, scratchSize, workDesc, 0x0));
178  std::cout << "Prepared the specified quantum circuit expectation value\n";

Set up the workspace

Now we can set up the required workspace buffer.

182  // Attach the workspace buffer
183  int64_t worksize {0};
184  HANDLE_CUTN_ERROR(cutensornetWorkspaceGetMemorySize(cutnHandle,
185                                                      workDesc,
186                                                      CUTENSORNET_WORKSIZE_PREF_RECOMMENDED,
187                                                      CUTENSORNET_MEMSPACE_DEVICE,
188                                                      CUTENSORNET_WORKSPACE_SCRATCH,
189                                                      &worksize));
190  std::cout << "Required scratch GPU workspace size (bytes) = " << worksize << std::endl;
191  if(worksize <= scratchSize) {
192    HANDLE_CUTN_ERROR(cutensornetWorkspaceSetMemory(cutnHandle, workDesc, CUTENSORNET_MEMSPACE_DEVICE,
193                      CUTENSORNET_WORKSPACE_SCRATCH, d_scratch, worksize));
194  }else{
195    std::cout << "ERROR: Insufficient workspace size on Device!\n";
196    std::abort();
197  }
198  std::cout << "Set the workspace buffer\n";

Compute the requested expectation value

Once everything has been set up, we compute the requested expectation value and print it. Note that the returned expectation value is not normalized. The 2-norm of the tensor network state is returned as a separate argument.

202  // Compute the specified quantum circuit expectation value
203  std::complex<double> expectVal{0.0,0.0}, stateNorm{0.0,0.0};
204  HANDLE_CUTN_ERROR(cutensornetExpectationCompute(cutnHandle, expectation, workDesc,
205                    static_cast<void*>(&expectVal), static_cast<void*>(&stateNorm), 0x0));
206  std::cout << "Computed the specified quantum circuit expectation value\n";
207  std::cout << "Expectation value = (" << expectVal.real() << ", " << expectVal.imag() << ")\n";
208  std::cout << "State 2-norm = (" << stateNorm.real() << ", " << stateNorm.imag() << ")\n";

Free resources

212  // Destroy the workspace descriptor
213  HANDLE_CUTN_ERROR(cutensornetDestroyWorkspaceDescriptor(workDesc));
214  std::cout << "Destroyed the workspace descriptor\n";
215
216  // Destroy the quantum circuit expectation value
217  HANDLE_CUTN_ERROR(cutensornetDestroyExpectation(expectation));
218  std::cout << "Destroyed the quantum circuit state expectation value\n";
219
220  // Destroy the tensor network operator
221  HANDLE_CUTN_ERROR(cutensornetDestroyNetworkOperator(hamiltonian));
222  std::cout << "Destroyed the tensor network operator\n";
223
224  // Destroy the quantum circuit state
225  HANDLE_CUTN_ERROR(cutensornetDestroyState(quantumState));
226  std::cout << "Destroyed the quantum circuit state\n";
227
228  HANDLE_CUDA_ERROR(cudaFree(d_scratch));
229  HANDLE_CUDA_ERROR(cudaFree(d_gateCX));
230  HANDLE_CUDA_ERROR(cudaFree(d_gateZ));
231  HANDLE_CUDA_ERROR(cudaFree(d_gateY));
232  HANDLE_CUDA_ERROR(cudaFree(d_gateX));
233  HANDLE_CUDA_ERROR(cudaFree(d_gateH));
234  std::cout << "Freed memory on GPU\n";
235
236  // Finalize the cuTensorNet library
237  HANDLE_CUTN_ERROR(cutensornetDestroy(cutnHandle));
238  std::cout << "Finalized the cuTensorNet library\n";
239
240  return 0;
241}