Solving Max-Cut Problem with QAOA¶

Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.

Overview¶

In the tutorial on Quantum Approximate Optimization Algorithm, we talked about how to encode a classical combinatorial optimization problem into a quantum optimization problem and slove it with Quantum Approximate Optimization Algorithm [1] (QAOA). In this tutorial, we will take the Max-Cut Problem as an example to further elaborate on QAOA.

Max-Cut Problem¶

The Max-Cut Problem is a common combinatorial optimization problem in graph theory, and it has important applications in statistical physics and circuit design. The maximum cut problem is an NP-hard problem, so there is no efficient algorithm that can solve this problem perfectly.

In graph theory, a graph is represented by a pair of sets G=(V,E)G=(V,E), where the elements in the set VV are the vertices of the graph, and each element in the set EE is a pair of vertices, representing an edge connecting these two vertices. For example, the graph in the figure below is represented by V={0,1,2,3}V={0,1,2,3} and E={(0,1),(1,2),(2,3),(3,0)}E={(0,1),(1,2),(2,3),(3,0)}.

G

Figure 1: A graph with four vertices and four edges

A cut on a graph refers to a partition of the graph's vertex set VV into two disjoint sets. Each cut corresponds to a set of edges, in which the two vertices of each edge are divided into different sets. So we can define the size of this cut as the size of this set of edges, that is, the number of edges being cut. The Max-Cut Problem is to find a cut that maximizes the number of edges being cut. Figure 2 shows a maximum cut of the graph in Figure 1. The size of the maximum cut is 44, which means that all edges in the graph are cut.

Max cut on G

Figure 2: A maximum cut of the graph in Figure 1

Assuming that the input graph G=(V,E)G=(V,E) has n=|V|n=|V| vertices and m=|E|m=|E| edges, we can describe the Max-Cut Problem as a combinatorial optimization problem with nn bits and mm clauses. Each bit corresponds to a vertex vv in the graph GG, and its value zvzv is 00 or 11, corresponding to the vertex belonging to the set S0S0 or S1S1, respectively. Thus, each value zz of these nn bits corresponds to a distinct cut. Each clause corresponds to an edge (u,v)(u,v) in the graph GG. A clause requires that the two vertices connected by its corresponding edge take different values, namely zu≠zvzu≠zv, which means the edge is cut. In other words, when the two vertices connected by the edge are divided into different sets, we say that the clause is satisfied. Therefore, for each edge (u,v)(u,v) in the graph GG, we have

C(u,v)(z)=zu+zv−2zuzv,(1)(1)C(u,v)(z)=zu+zv−2zuzv,

where C(u,v)(z)=1C(u,v)(z)=1 if and only if the edge is cut. Otherwise, the function is equal to 00. The objective function of the entire combinatorial optimization problem is

C(z)=∑(u,v)∈EC(u,v)(z)=∑(u,v)∈Ezu+zv−2zuzv.(2)(2)C(z)=∑(u,v)∈EC(u,v)(z)=∑(u,v)∈Ezu+zv−2zuzv.

Therefore, to solve the maximum cut problem is to find a value zz that maximizes the objective function in equation (2).

Encoding Max-Cut Problem¶

Here we take the Max-Cut Problem as an example to further elaborate on QAOA. In order to transform the Max-Cut Problem into a quantum problem, we need to use nn qubits, where each qubit corresponds to a vertex in the graph GG. A qubit being in a quantum state |0⟩|0⟩ or |1⟩|1⟩ indicates that its corresponding vertex belongs to the set S0S0 or S1S1, respectively. It is worth noting that |0⟩|0⟩ and |1⟩|1⟩ are the two eigenstates of Pauli ZZ gate, and their eigenvalues are respectively 11 and −1−1, namely

Z|0⟩=|0⟩,Z|1⟩=−|1⟩.(3)(4)(3)Z|0⟩=|0⟩,(4)Z|1⟩=−|1⟩.

Therefore, we can use Pauli ZZ gate to construct the Hamiltonian HCHC of the Max-Cut Problem. Because mapping f(x):x→(x+1)/2f(x):x→(x+1)/2 maps −1−1 to 00 and 11 to 11, we can replace zz in equation (2) with (Z+I)/2(Z+I)/2 (II is the identity matrix) to get the Hamiltonian corresponding to the objective function of the original problem:

HC=∑(u,v)∈EZu+I2+Zv+I2−2⋅Zu+I2Zv+I2=∑(u,v)∈EZu+Zv+2I−(ZuZv+Zu+Zv+I)2=∑(u,v)∈EI−ZuZv2.(5)(6)(7)(5)HC=∑(u,v)∈EZu+I2+Zv+I2−2⋅Zu+I2Zv+I2(6)=∑(u,v)∈EZu+Zv+2I−(ZuZv+Zu+Zv+I)2(7)=∑(u,v)∈EI−ZuZv2.

The expected value of this Hamiltonian for a quantum state |ψ⟩|ψ⟩ is

⟨ψ|HC|ψ⟩=⟨ψ|∑(u,v)∈EI−ZuZv2|ψ⟩=⟨ψ|∑(u,v)∈EI2|ψ⟩−⟨ψ|∑(u,v)∈EZuZv2|ψ⟩=|E|2−12⟨ψ|∑(u,v)∈EZuZv|ψ⟩.(8)(9)(10)(8)⟨ψ|HC|ψ⟩=⟨ψ|∑(u,v)∈EI−ZuZv2|ψ⟩(9)=⟨ψ|∑(u,v)∈EI2|ψ⟩−⟨ψ|∑(u,v)∈EZuZv2|ψ⟩(10)=|E|2−12⟨ψ|∑(u,v)∈EZuZv|ψ⟩.

If we define

HD=−∑(u,v)∈EZuZv,(11)(11)HD=−∑(u,v)∈EZuZv,

then finding the quantum state |ψ⟩|ψ⟩ that maximizes ⟨ψ|HC|ψ⟩⟨ψ|HC|ψ⟩ is equivalent to finding the quantum state |ψ⟩|ψ⟩ such that ⟨ψ|HD|ψ⟩⟨ψ|HD|ψ⟩ is the largest.

Paddle Quantum Implementation¶

Now, let's implement QAOA with Paddle Quantum to solve the Max-Cut Problem. There are many ways to find the parameters →γ,→βγ→,β→. Here we use the gradient descent method in classical machine learning.

To implement QAOA with Paddle Quantum, the first thing to do is to import the required packages. Among them, the networkx package can help us handle graphs conveniently.

In [1]:
from IPython.core.display import HTML
display(HTML("<style>pre { white-space: pre !important; }</style>"))
In [2]:
# Import related modules from Paddle Quantum and PaddlePaddle
import paddle
from paddle_quantum.ansatz import Circuit
from paddle_quantum.qinfo import pauli_str_to_matrix
from paddle_quantum.loss import ExpecVal
from paddle_quantum import Hamiltonian

# Import additional packages needed
import numpy as np
from numpy import pi as PI
import matplotlib.pyplot as plt
import networkx as nx

import warnings
warnings.filterwarnings("ignore")

/usr/local/Caskroom/miniconda/base/envs/pq_new/lib/python3.8/site-packages/paddle/tensor/creation.py:125: DeprecationWarning: `np.object` is a deprecated alias for the builtin `object`. To silence this warning, use `object` by itself. Doing this will not modify any behavior and is safe. 
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  if data.dtype == np.object:
/usr/local/Caskroom/miniconda/base/envs/pq_new/lib/python3.8/site-packages/paddle/tensor/creation.py:125: DeprecationWarning: `np.object` is a deprecated alias for the builtin `object`. To silence this warning, use `object` by itself. Doing this will not modify any behavior and is safe. 
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  if data.dtype == np.object:

Next, we generate the graph GG in the Max-Cut Problem. For the convenience of computation, the vertices here are labeled starting from 00.

In [3]:
# n is the number of vertices in the graph G, which is also the number of qubits
n = 4
G = nx.Graph()
V = range(n)
G.add_nodes_from(V)
E = [(0, 1), (1, 2), (2, 3), (3, 0), (1, 3)]
G.add_edges_from(E)

# Print out the generated graph G
pos = nx.circular_layout(G)
options = {
    "with_labels": True,
    "font_size": 20,
    "font_weight": "bold",
    "font_color": "white",
    "node_size": 2000,
    "width": 2
}
nx.draw_networkx(G, pos, **options)
ax = plt.gca()
ax.margins(0.20)
plt.axis("off")
plt.show()
No description has been provided for this image

Encoding Hamiltonian¶

In Paddle Quantum, a Hamiltonian can be input in the form of list. Here we construct the Hamiltonian HDHD in equation (11).

In [4]:
# Construct the Hamiltonian H_D in the form of list
H_D_list = []
for (u, v) in E:
    H_D_list.append([-1.0,'z'+str(u) +',z' + str(v)])
print(H_D_list)

[[-1.0, 'z0,z1'], [-1.0, 'z1,z2'], [-1.0, 'z2,z3'], [-1.0, 'z3,z0'], [-1.0, 'z1,z3']]

As you can see, in this example, the Hamiltonian HDHD is

HD=−Z0Z1−Z1Z2−Z2Z3−Z3Z0−Z1Z3.(12)(12)HD=−Z0Z1−Z1Z2−Z2Z3−Z3Z0−Z1Z3.

We can view the matrix form of the Hamiltonian HDHD and get information of its eigenvalues:

In [5]:
# Convert Hamiltonian H_D from list form to matrix form
H_D_matrix = pauli_str_to_matrix(H_D_list, n)
# Take out the elements on the diagonal of H_D
H_D_diag = np.diag(H_D_matrix).real
# Get the maximum eigenvalue of H_D
H_max = np.max(H_D_diag)

print(H_D_diag)
print('H_max:', H_max)

[-5.  1. -1.  1.  1.  3.  1. -1. -1.  1.  3.  1.  1. -1.  1. -5.]
H_max: 3.0

Building the QAOA circuit¶

Earlier we introduced that QAOA needs to apply two unitary transformations UC(γ)UC(γ) and UB(β)UB(β) alternately on the initial state |s⟩=|+⟩⊗n|s⟩=|+⟩⊗n. Here, we use the quantum gates and quantum circuit templates provided in Paddle Quantum to build a quantum circuit to achieve this step. It should be noted that in the Max-Cut Problem, we simplify the problem of maximizing the expected value of the Hamiltonian HCHC to the problem of maximizing the expected value of the Hamiltonian HDHD, so the unitary transformations to be used are UD(γ)UD(γ) and UB(β)UB(β). By alternately placing two circuit modules with adjustable parameters, we are able to build a QAOA circuit

UB(βp)UD(γp)⋯UB(β1)UD(γ1),(13)(13)UB(βp)UD(γp)⋯UB(β1)UD(γ1),

where UD(γ)=e−iγHDUD(γ)=e−iγHD can be constructed with the circuit in the figure below. Another unitary transformation UB(β)UB(β) is equivalent to applying a RxRx gate to each qubit.

U_D circuit

Figure 3: Quantum circuit of unitary transformation eiγZ⊗ZeiγZ⊗Z

Therefore, the quantum circuit that realizes a layer of unitary transformation UB(β)UD(γ)UB(β)UD(γ) is shown in Figure 4.

U_BU_D circuit

Figure 4: Quantum circuit of unitary transformation UB(β)UD(γ)UB(β)UD(γ)

In Paddle Quantum, the default initial state of each qubit is |0⟩|0⟩ (the initial state can be customized by input parameters). We can add a layer of Hadamard gates to change the state of each qubit from |0⟩|0⟩ to |+⟩|+⟩ so that we get the initial state |s⟩=|+⟩⊗n|s⟩=|+⟩⊗n required by QAOA. In Paddle Quantum, we can add a layer of Hadamard gates to the quantum circuit by calling superposition_layer().

In [6]:
def circuit_QAOA(num_qubits, depth, edges, vertices):
    # Initialize the quantum circuit of n qubits
    cir = Circuit(num_qubits)
    # Prepare quantum state |s>
    cir.superposition_layer()
    # Build a circuit with p layers of U_D
    cir.qaoa_layer(edges, vertices, depth)
    
    return cir

After running the constructed QAOA quantum circuit, we obtain the output state

|→γ,→β⟩=UB(βp)UD(γp)⋯UB(β1)UD(γ1)|s⟩.(14)(14)|γ→,β→⟩=UB(βp)UD(γp)⋯UB(β1)UD(γ1)|s⟩.

Calculating the loss function¶

From the output state of the circuit built in the previous step, we can calculate the objective function of the maximum cut problem

Fp(→γ,→β)=⟨→γ,→β|HD|→γ,→β⟩.(15)(15)Fp(γ→,β→)=⟨γ→,β→|HD|γ→,β→⟩.

To maximize the objective function is equivalent to minimizing −Fp−Fp. Therefore, we define L(→γ,→β)=−Fp(→γ,→β)L(γ→,β→)=−Fp(γ→,β→) as the loss function, that is, the function to be minimized. Then, we use a classical optimization algorithm to find the optimal parameters →γ,→βγ→,β→. The following code shows the loss function constructed with Paddle Quantum:

In [7]:
# construct the loss function
loss_func = ExpecVal(Hamiltonian(H_D_list))

Training quantum neural network¶

After defining the quantum neural network for QAOA, we use the gradient descent method to update the parameters in the network to maximize the expected value in equation (15).

In [8]:
depth = 4      # Number of layers in the quantum circuit
ITR = 120  # Number of training iterations
LR = 0.1   # Learning rate of the optimization method based on gradient descent
SEED = 1024 # Set global RNG seed

Here, we optimize the loss function defined above in PaddlePaddle.

In [9]:
paddle.seed(SEED)

cir = circuit_QAOA(n, depth, E, V)
# Use Adam optimizer
opt = paddle.optimizer.Adam(learning_rate=LR, parameters=cir.parameters())

for itr in range(1, ITR + 1):
    state = cir()
    # Calculate the gradient and optimize
    loss = -loss_func(state)
    loss.backward()
    opt.minimize(loss)
    opt.clear_grad()
    if itr % 10 == 0:
        print("iter:", itr, "  loss:", "%.4f" % loss.numpy())

iter: 10   loss: -2.5212
iter: 20   loss: -2.7688
iter: 30   loss: -2.9486
iter: 40   loss: -2.9832
iter: 50   loss: -2.9907
iter: 60   loss: -2.9969
iter: 70   loss: -2.9990
iter: 80   loss: -2.9997
iter: 90   loss: -2.9999
iter: 100   loss: -3.0000
iter: 110   loss: -3.0000
iter: 120   loss: -3.0000

Decoding the quantum solution¶

After obtaining the minimum value of the loss function and the corresponding set of parameters →γ∗,→β∗γ→∗,β→∗, our task has not been completed. In order to obtain an approximate solution to the Max-Cut Problem, it is necessary to decode the solution to the classical optimization problem from the quantum state |→γ∗,→β∗⟩|γ→∗,β→∗⟩ output by QAOA. Physically, to decode a quantum state, we need to measure it and then calculate the probability distribution of the measurement results:

p(z)=|⟨z|→γ∗,→β∗⟩|2.(16)(16)p(z)=|⟨z|γ→∗,β→∗⟩|2.

Usually, the greater the probability of a certain bit string, the greater the probability that it corresponds to an optimal solution of the Max-Cut problem.

Paddle Quantum provides a function to view the probability distribution of the measurement results of the state output by the QAOA quantum circuit:

In [10]:
state = cir()
# Repeat the simulated measurement of the circuit output state 1024 times
prob_measure = state.measure(plot=True)
No description has been provided for this image

After measurement, we can find the bit string with the highest probability of occurrence. Let the vertices whose bit values are 00 in the bit string belong to the set S0S0 and the vertices whose bit values are 11 belong to the set S1S1. The set of edges between these two vertex sets is a possible maximum cut of the graph.

The following code selects the bit string with the greatest chance of appearing in the measurement result, then maps it back to the classic solution, and draws the corresponding maximum cut:

  • The red vertex belongs to the set S0S0,
  • The blue vertex belongs to the set S1S1,
  • The dashed line indicates the edge being cut.
In [11]:
# Find the most frequent bit string in the measurement results
cut_bitstring = max(prob_measure, key=prob_measure.get)
print("The bit string form of the cut found:", cut_bitstring)

# Draw the cut corresponding to the bit string obtained above on the graph
node_cut = ["blue" if cut_bitstring[v] == "1" else "red" for v in V]

edge_cut = []
for u in range(n):
    for v in range(u + 1, n):
        if (u, v) in E or (v, u) in E:
            if cut_bitstring[u] == cut_bitstring[v]:
                edge_cut.append("solid")
            else:
                edge_cut.append("dashed")
nx.draw(G, pos, node_color=node_cut, style=edge_cut, **options)
ax = plt.gca()
ax.margins(0.20)
plt.axis("off")
plt.show()

The bit string form of the cut found: 0101
No description has been provided for this image

As you can see, in this example, QAOA has found a maximum cut on the graph.

References¶

[1] Farhi, E., Goldstone, J. & Gutmann, S. A Quantum Approximate Optimization Algorithm. arXiv:1411.4028 (2014).