# Quantum Approximate Optimization Algorithm¶

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)
```

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/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, [&](std::vector<double> params) {
return std::make_tuple(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;
}
```