Quantum Approximate Optimization Algorithm

Python

Let’s now see how we can implement the Quantum Approximate Optimization Algorithm (QAOA) to compute the Max-Cut of a rectangular graph by leveraging cudaq:vqe. For more on the QAOA algorithm and the Max-Cut problem, refer to this paper.

import cudaq
from cudaq import spin

from typing import List

import numpy as np

# Here we build up a kernel for QAOA with `p` layers, with each layer
# containing the alternating set of unitaries corresponding to the problem
# and the mixer Hamiltonians. The algorithm leverages the VQE algorithm
# to compute the Max-Cut of a rectangular graph illustrated below.

#       v0  0---------------------0 v1
#           |                     |
#           |                     |
#           |                     |
#           |                     |
#       v3  0---------------------0 v2
# The Max-Cut for this problem is 0101 or 1010.

# The problem Hamiltonian
hamiltonian = 0.5 * spin.z(0) * spin.z(1) + 0.5 * spin.z(1) * spin.z(2) \
       + 0.5 * spin.z(0) * spin.z(3) + 0.5 * spin.z(2) * spin.z(3)

# Problem parameters.
qubit_count: int = 4
layer_count: int = 2
parameter_count: int = 2 * layer_count


@cudaq.kernel
def kernel_qaoa(qubit_count: int, layer_count: int, thetas: List[float]):
    """QAOA ansatz for Max-Cut"""
    qvector = cudaq.qvector(qubit_count)

    # Create superposition
    h(qvector)

    # Loop over the layers
    for layer in range(layer_count):
        # Loop over the qubits
        # Problem unitary
        for qubit in range(qubit_count):
            x.ctrl(qvector[qubit], qvector[(qubit + 1) % qubit_count])
            rz(2.0 * thetas[layer], qvector[(qubit + 1) % qubit_count])
            x.ctrl(qvector[qubit], qvector[(qubit + 1) % qubit_count])

        # Mixer unitary
        for qubit in range(qubit_count):
            rx(2.0 * thetas[layer + layer_count], qvector[qubit])


# Specify the optimizer and its initial parameters. Make it repeatable.
cudaq.set_random_seed(13)
optimizer = cudaq.optimizers.COBYLA()
np.random.seed(13)
optimizer.initial_parameters = np.random.uniform(-np.pi / 8.0, np.pi / 8.0,
                                                 parameter_count)
print("Initial parameters = ", optimizer.initial_parameters)


# Define the objective, return `<state(params) | H | state(params)>`
def objective(parameters):
    return cudaq.observe(kernel_qaoa, hamiltonian, qubit_count, layer_count,
                         parameters).expectation()


# Optimize!
optimal_expectation, optimal_parameters = optimizer.optimize(
    dimensions=parameter_count, function=objective)

# Print the optimized value and its parameters
print("Optimal value = ", optimal_expectation)
print("Optimal parameters = ", optimal_parameters)

# Sample the circuit using the optimized parameters
counts = cudaq.sample(kernel_qaoa, qubit_count, layer_count, optimal_parameters)
print(counts)

C++

Let’s now see how we can implement the Quantum Approximate Optimization Algorithm (QAOA) to compute the Max-Cut of a rectangular graph For more on the QAOA algorithm and the Max-Cut problem, refer to this paper.

#include <cudaq.h>
#include <cudaq/algorithm.h>
#include <cudaq/gradients.h>
#include <cudaq/optimizers.h>
#include <cudaq/spin_op.h>

// Here we build up a CUDA-Q kernel for QAOA with p layers, with each
// layer containing the alternating set of unitaries corresponding to the
// problem and the mixer Hamiltonians. The algorithm leverages the CUDA-Q
// VQE support to compute the Max-Cut of a rectangular graph illustrated below.
//
//        v0  0---------------------0 v1
//            |                     |
//            |                     |
//            |                     |
//            |                     |
//        v3  0---------------------0 v2
// The Max-Cut for this problem is 0101 or 1010.

struct ansatz {
  void operator()(std::vector<double> theta, const int n_qubits,
                  const int n_layers) __qpu__ {
    cudaq::qvector q(n_qubits);

    // Prepare the initial state by superposition
    h(q);

    // Loop over all the layers
    for (int i = 0; i < n_layers; ++i) {
      // Problem Hamiltonian
      for (int j = 0; j < n_qubits; ++j) {

        x<cudaq::ctrl>(q[j], q[(j + 1) % n_qubits]);
        rz(2.0 * theta[i], q[(j + 1) % n_qubits]);
        x<cudaq::ctrl>(q[j], q[(j + 1) % n_qubits]);
      }

      for (int j = 0; j < n_qubits; ++j) {
        // Mixer Hamiltonian
        rx(2.0 * theta[i + n_layers], q[j]);
      }
    }
  }
};

int main() {

  using namespace cudaq::spin;

  cudaq::set_random_seed(13); // set for repeatability

  // Problem Hamiltonian
  const cudaq::spin_op Hp = 0.5 * z(0) * z(1) + 0.5 * z(1) * z(2) +
                            0.5 * z(0) * z(3) + 0.5 * z(2) * z(3);

  // Problem parameters
  const int n_qubits = 4;
  const int n_layers = 2;
  const int n_params = 2 * n_layers;

  // Instantiate the optimizer
  cudaq::optimizers::cobyla optimizer; // gradient-free COBYLA

  // Set initial values for the optimization parameters
  optimizer.initial_parameters = cudaq::random_vector(
      -M_PI / 8.0, M_PI / 8.0, n_params, std::mt19937::default_seed);

  // Call the optimizer
  auto [opt_val, opt_params] =
      cudaq::vqe(ansatz{}, Hp, optimizer, n_params, n_qubits, n_layers);

  // Print the optimized value and the parameters
  printf("Optimal value = %.16lf\n", opt_val);
  printf("Optimal params = (%.16lf, %.16lf, %.16lf, %.16lf) \n", opt_params[0],
         opt_params[1], opt_params[2], opt_params[3]);

  // Sample the circuit using optimized parameters
  auto counts = cudaq::sample(ansatz{}, opt_params, n_qubits, n_layers);

  // Dump the states and the counts
  counts.dump();

  return 0;
}