r"""
Feedback-Based Quantum Optimization (FALQON)
============================================
.. meta::
:property="og:description": Solve combinatorial optimization problems without a classical optimizer
:property="og:image": https://pennylane.ai/qml/_images/falqon_thumbnail.png
.. related::
tutorial_qaoa_intro Intro to QAOA
tutorial_qaoa_maxcut QAOA for MaxCut
*Authors: David Wakeham and Jack Ceroni. Posted: 21 May 2021. Last updated: 21 May 2021.*
-----------------------------
While the
`Quantum Approximate Optimization Algorithm (QAOA) `__
is one of the best-known processes for solving combinatorial optimization problems with quantum computers,
it has a major drawback: convergence isn't guaranteed, as the optimization procedure can become "stuck" in local minima.
.. figure:: ../demonstrations/falqon/global_min_graph.png
:align: center
:width: 70%
This tutorial explores the **FALQON** algorithm introduced in `a recent paper by Magann et al. `__
It is similar in spirit to QAOA, but uses iterative feedback steps rather than a global optimization
over parameters, avoiding the use of a classical optimizer.
In this demo, we will implement FALQON to solve the MaxClique problem in graph theory, perform benchmarking, and
combine FALQON with QAOA to create a powerful optimization procedure.
.. note::
If you are not familiar with QAOA, we recommend checking out the
`Intro to QAOA tutorial `__,
since many of the same ideas carry over and will be used throughout this demonstration.
Theory
------
To solve a combinatorial optimization problem with a quantum computer, a typical strategy is to encode
the solution to the problem as the ground state of a *cost Hamiltonian* :math:`H_c`. Then, we use some procedure to drive
an initial state into the ground state of :math:`H_c`.
FALQON falls under this broad scheme.
Consider a quantum system governed by a Hamiltonian of the form :math:`H = H_c + \beta(t) H_d`. These kinds of
Hamiltonians appear often in the theory of `quantum control `__, a
field of inquiry which studies how a quantum system can be driven from one state to another.
The choice of :math:`\beta(t)` corresponds to a "driving strategy", which partially determines how the system evolves
with time.
Suppose our objective is to drive some quantum system
to the ground state of :math:`H_c`. It is a reasonable goal to
construct a quantum control process such that the energy expectation :math:`\langle H_c \rangle_t` decreases with time:
.. math:: \frac{d}{dt} \langle H_c\rangle_t = \frac{d}{dt} \langle \psi(t)|H_c|\psi(t)\rangle = i \beta(t)\langle [H_d, H_c] \rangle_t \leq 0,
where the product rule and
`the Schrödinger equation `__ are used to derive the above formula.
If we pick :math:`\beta(t) = -\langle i[H_d, H_c] \rangle_t,` so that
.. math:: \frac{d}{dt} \langle H_c\rangle_t = -|\langle i[H_d, H_c] \rangle_t|^2 \leq 0,
then :math:`\langle H_c \rangle` is guaranteed to strictly decrease, as desired!
Thus, if we evolve some initial state :math:`|\psi_0\rangle` under the time evolution operator :math:`U` corresponding to :math:`H`,
.. math:: U(T) = \mathcal{T} \exp \Big[ -i \displaystyle\int_{0}^{T} H(t) \ dt \Big] \approx \mathcal{T} \exp \Big[ -i \displaystyle\sum_{k = 0}^{T/\Delta t} H( k \Delta t) \Delta t \Big],
where :math:`\mathcal{T}` is the `time-ordering operator `__ and :math:`\Delta t` is some small time step,
then the energy expectation will strictly decrease, for a large enough value of :math:`T`. This is exactly the procedure used by FALQON to minimize :math:`\langle H_c \rangle`.
In general, implementing a time evolution unitary in a quantum circuit is
difficult, so we use a
`Trotter-Suzuki decomposition `__
to perform approximate time evolution. We then have
.. math:: U(T) \approx \mathcal{T} \exp \Big[ -i \displaystyle\sum_{k = 0}^{T/\Delta t} H( k \Delta t) \Delta t \Big] \approx
e^{-i\beta_n H_d \Delta t} e^{-iH_c \Delta t} \cdots e^{-i\beta_1 H_d \Delta t} e^{-iH_c \Delta t} = U_d(\beta_n) U_c \cdots U_d(\beta_1) U_c,
where :math:`n = T/\Delta t` and :math:`\beta_k = \beta(k\Delta t)`.
For each layer of the time evolution, the value :math:`\beta_k` is required. However,
:math:`\beta_k` is dependent on the state of the system at some time. Recall that
.. math:: \beta(t) = - \langle \psi(t) | i [H_d, H_c] | \psi(t) \rangle.
We let :math:`A(t) := i\langle [H_d, H_c] \rangle_t`. Our strategy is to obtain the values
:math:`\beta_k` recursively, by finding the value of :math:`A(t)` for the **previous time step**. We then set
.. math:: \beta_{k+1} = -A_k = -A(k\Delta t).
This leads to the FALQON algorithm as a recursive process. On step :math:`k`, we perform the following three substeps:
1. Prepare the state :math:`|\psi_k\rangle = U_d(\beta_k) U_c \cdots U_d(\beta_1) U_c|\psi_0\rangle`.
2. Measure the expectation value :math:`A_k = \langle i[H_c, H_d]\rangle_{k \Delta t}`.
3. Set :math:`\beta_{k+1} = -A_k`.
We repeat for all :math:`k` from :math:`1` to :math:`n`, where :math:`n` is a hyperparameter.
At the final step we evaluate :math:`\langle H_c \rangle`,
which gives us an approximation for the ground state of :math:`H_c`.
.. figure:: ../demonstrations/falqon/falqon.png
:align: center
:width: 80%
"""
######################################################################
# Simulating FALQON with PennyLane
# --------------------------------
# To begin, we import the necessary dependencies:
#
import pennylane as qml
from pennylane import numpy as np
from matplotlib import pyplot as plt
from pennylane import qaoa as qaoa
import networkx as nx
######################################################################
# In this demonstration, we will be using FALQON to solve the
# `maximum clique (MaxClique) problem `__: finding the
# largest complete subgraph of some graph :math:`G`. For example, the following graph's maximum
# clique is coloured in red:
#
# .. figure:: ../demonstrations/falqon/max_clique.png
# :align: center
# :width: 90%
#
# We attempt to find the maximum clique of the graph below:
#
edges = [(0, 1), (1, 2), (2, 0), (2, 3), (1, 4)]
graph = nx.Graph(edges)
nx.draw(graph, with_labels=True, node_color="#e377c2")
######################################################################
# We must first encode this combinatorial problem into a cost Hamiltonian :math:`H_c`. This ends up being
#
# .. math:: H_c = 3 \sum_{(i, j) \in E(\bar{G})} (Z_i Z_j - Z_i - Z_j) + \displaystyle\sum_{i \in V(G)} Z_i,
#
# where each qubit is a node in the graph, and the states :math:`|0\rangle` and :math:`|1\rangle`
# represent whether the vertex has been marked as part of the clique, as is the case for `most standard QAOA encoding
# schemes `__.
# Note that :math:`\bar{G}` is the complement of :math:`G`: the graph formed by connecting all nodes that **do not** share
# an edge in :math:`G`.
#
# In addition to defining :math:`H_c`, we also require a driver Hamiltonian :math:`H_d` which does not commute
# with :math:`H_c`. The driver Hamiltonian's role is similar to that of the mixer Hamiltonian in QAOA.
# To keep things simple, we choose a sum over Pauli :math:`X` operations on each qubit:
#
# .. math:: H_d = \displaystyle\sum_{i \in V(G)} X_i.
#
# These Hamiltonians come nicely bundled together in the PennyLane QAOA module:
#
cost_h, driver_h = qaoa.max_clique(graph, constrained=False)
print("Cost Hamiltonian")
print(cost_h)
print("Driver Hamiltonian")
print(driver_h)
######################################################################
# One of the main ingredients in the FALQON algorithm is the operator :math:`i [H_d, H_c]`. In
# the case of MaxClique, we can write down the commutator :math:`[H_d, H_c]` explicitly:
#
# .. math:: [H_d, H_c] = 3 \displaystyle\sum_{k \in V(G)} \displaystyle\sum_{(i, j) \in E(\bar{G})} \big( [X_k, Z_i Z_j] - [X_k, Z_i]
# - [X_k, Z_j] \big) + 3 \displaystyle\sum_{i \in V(G)} \displaystyle\sum_{j \in V(G)} [X_i, Z_j].
#
# There are two distinct commutators that we must calculate, :math:`[X_k, Z_j]` and :math:`[X_k, Z_i Z_j]`.
# This is straightforward as we know exactly what the
# `commutators of the Pauli matrices `__ are.
# We have:
#
# .. math:: [X_k, Z_j] = -2 i \delta_{kj} Y_k \ \ \ \text{and} \ \ \ [X_k, Z_i Z_j] = -2 i \delta_{ik} Y_k Z_j - 2i \delta_{jk} Z_i Y_k,
#
# where :math:`\delta_{kj}` is the `Kronecker delta `__. Therefore it
# follows from substitution into the above equation and multiplication by :math:`i` that:
#
# .. math:: i [H_d, H_c] = 6 \displaystyle\sum_{k \in V(G)} \displaystyle\sum_{(i, j) \in E(\bar{G})} \big( \delta_{ki} Y_k Z_j +
# \delta_{kj} Z_{i} Y_{k} - \delta_{ki} Y_k - \delta_{kj} Y_k \big) + 6 \displaystyle\sum_{i \in V(G)} Y_{i}.
#
# This new operator has quite a few terms! Therefore, we write a short method which computes it for us, and returns
# a :class:`~.pennylane.Hamiltonian` object. Note that this method works for any graph:
#
def build_hamiltonian(graph):
H = qml.Hamiltonian([], [])
# Computes the complement of the graph
graph_c = nx.complement(graph)
for k in graph_c.nodes:
# Adds the terms in the first sum
for edge in graph_c.edges:
i, j = edge
if k == i:
H += 6 * (qml.PauliY(k) @ qml.PauliZ(j) - qml.PauliY(k))
if k == j:
H += 6 * (qml.PauliZ(i) @ qml.PauliY(k) - qml.PauliY(k))
# Adds the terms in the second sum
H += 6 * qml.PauliY(k)
return H
print("MaxClique Commutator")
print(build_hamiltonian(graph))
######################################################################
# We can now build the FALQON algorithm. Our goal is to evolve some initial state under the Hamiltonian :math:`H`,
# with our chosen :math:`\beta(t)`. We first define one layer of the Trotterized time evolution, which is of
# the form :math:`U_d(\beta_k) U_c`. Note that we can use the :class:`~.pennylane.templates.ApproxTimeEvolution` template:
def falqon_layer(beta_k, cost_h, driver_h, delta_t):
qml.ApproxTimeEvolution(cost_h, delta_t, 1)
qml.ApproxTimeEvolution(driver_h, delta_t * beta_k, 1)
######################################################################
# We then define a method which returns a FALQON ansatz corresponding to a particular cost Hamiltonian, driver
# Hamiltonian, and :math:`\Delta t`. This involves multiple repetitions of the "FALQON layer" defined above. The
# initial state of our circuit is an even superposition:
def build_maxclique_ansatz(cost_h, driver_h, delta_t):
def ansatz(beta, **kwargs):
layers = len(beta)
for w in dev.wires:
qml.Hadamard(wires=w)
qml.layer(
falqon_layer,
layers,
beta,
cost_h=cost_h,
driver_h=driver_h,
delta_t=delta_t
)
return ansatz
def expval_circuit(beta, measurement_h):
ansatz = build_maxclique_ansatz(cost_h, driver_h, delta_t)
ansatz(beta)
return qml.expval(measurement_h)
######################################################################
# Finally, we implement the recursive process, where FALQON is able to determine the values
# of :math:`\beta_k`, feeding back into itself as the number of layers increases. This is
# straightforward using the methods defined above:
def max_clique_falqon(graph, n, beta_1, delta_t, dev):
comm_h = build_hamiltonian(graph) # Builds the commutator
cost_h, driver_h = qaoa.max_clique(graph, constrained=False) # Builds H_c and H_d
cost_fn = qml.QNode(expval_circuit, dev) # The ansatz + measurement circuit is executable
beta = [beta_1] # Records each value of beta_k
energies = [] # Records the value of the cost function at each step
for i in range(n):
# Adds a value of beta to the list and evaluates the cost function
beta.append(-1 * cost_fn(beta, measurement_h=comm_h)) # this call measures the expectation of the commuter hamiltonian
energy = cost_fn(beta, measurement_h=cost_h) # this call measures the expectation of the cost hamiltonian
energies.append(energy)
return beta, energies
######################################################################
# Note that we return both the list of :math:`\beta_k` values, as well as the expectation value of the cost Hamiltonian
# for each step.
#
# We can now run FALQON for our MaxClique problem! It is important that we choose :math:`\Delta t` small enough
# such that the approximate time evolution is close enough to the real time evolution, otherwise we the expectation
# value of :math:`H_c` may not strictly decrease. For this demonstration, we set :math:`\Delta t = 0.03`,
# :math:`n = 40`, and :math:`\beta_1 = 0`. These are comparable to the hyperparameters chosen in the original paper.
n = 40
beta_1 = 0.0
delta_t = 0.03
dev = qml.device("default.qubit", wires=graph.nodes) # Creates a device for the simulation
res_beta, res_energies = max_clique_falqon(graph, n, beta_1, delta_t, dev)
######################################################################
# We can then plot the expectation value of the cost Hamiltonian over the
# iterations of the algorithm:
#
plt.plot(range(n+1)[1:], res_energies)
plt.xlabel("Iteration")
plt.ylabel("Cost Function Value")
plt.show()
######################################################################
# The expectation value decreases!
#
# To get a better understanding of the performance of the FALQON algorithm,
# we can create a graph showing the probability of measuring each possible bit string.
# We define the following circuit, feeding in the optimal values of :math:`\beta_k`:
@qml.qnode(dev)
def prob_circuit():
ansatz = build_maxclique_ansatz(cost_h, driver_h, delta_t)
ansatz(res_beta)
return qml.probs(wires=dev.wires)
######################################################################
# Running this circuit gives us the following probability distribution:
#
probs = prob_circuit()
plt.bar(range(2**len(dev.wires)), probs)
plt.xlabel("Bit string")
plt.ylabel("Measurement Probability")
plt.show()
######################################################################
# The bit string occurring with the highest probability is the state :math:`|28\rangle = |11100\rangle`.
# This corresponds to nodes :math:`0`, :math:`1`, and :math:`2`, which is precisely the maximum clique.
# FALQON has solved the MaxClique problem 🤩.
#
graph = nx.Graph(edges)
cmap = ["#00b4d9"]*3 + ["#e377c2"]*2
nx.draw(graph, with_labels=True, node_color=cmap)
######################################################################
# Benchmarking FALQON
# -------------------
#
# After seeing how FALQON works, it is worth considering how well FALQON performs according to a set of benchmarking
# criteria on a batch of graphs. We generate graphs randomly using the
# `Erdos-Renyi model `__, where we start with
# the complete graph on :math:`n` vertices and then keep each edge with probability :math:`p`. We then find the maximum
# cliques on these graphs using the
# `Bron-Kerbosch algorithm `__. To benchmark FALQON,
# the relative error in the estimated minimum energy
#
# .. math:: r_A = \frac{\langle H_C\rangle - \langle H_C\rangle_\text{min}}{|\langle H_C\rangle_\text{min}|}
#
# makes a good figure of merit.
#
# Final results for :math:`r_A`, along with :math:`\beta`, are plotted below,
# with the number of FALQON layers on the horizontal axis. We have averaged over :math:`50` random graphs per node
# size, for sizes :math:`n = 6, 7, 8, 9`, with probability :math:`p = 0.1` of keeping an edge. Running FALQON for
# :math:`40` steps, with :math:`\Delta t = 0.01`, produces:
#
# .. figure:: ../demonstrations/falqon/bench.png
# :align: center
# :width: 60%
#
# The relative error decreases with the number of layers (as we expect from the construction) and graph size (suggesting the errors grows
# more slowly than the minimum energy).
# The exception is :math:`n = 9`, where the step size has become too large
# and the Trotter-Suzuki decomposition breaks down.
# The rate of decrease also slows down. Even though the algorithm will converge to the ground state, it won't always get
# there in a few time steps!
######################################################################
# Seeding QAOA with FALQON (Bird Seed 🦅)
# ---------------------------------------
#
# .. figure:: ../demonstrations/falqon/bird_seed.png
# :align: center
# :width: 90%
#
# Both FALQON and QAOA have unique benefits and drawbacks.
# While FALQON requires no classical optimization and is guaranteed to decrease the cost function
# with each iteration, its circuit depth grows linearly with the number of iterations. The benchmarking data also shows
# how the reduction in cost slows with each layer, and the additional burden of correctly tuning the time step. On the other hand, QAOA
# has a fixed circuit depth, but does require classical optimization, and is therefore subject to all of the drawbacks
# that come with probing a cost landscape for a set of optimal parameters.
#
# QAOA and FALQON also have many similarities, most notably, their circuit structure. Both
# involve alternating layers of time evolution operators corresponding to a cost and a mixer/driver Hamiltonian.
# The FALQON paper raises the idea of combining FALQON and QAOA to yield a new optimization algorithm that
# leverages the benefits of both. In this final section of the tutorial, we will implement this procedure in PennyLane.
#
# Suppose we want to run a QAOA circuit of depth :math:`p`. Our ansatz will be of the form
#
# .. math:: U_{\text{QAOA}} = e^{-i \alpha_p H_m} e^{-i \gamma_p H_c} \cdots e^{-i \alpha_1 H_m} e^{-i \gamma_1 H_c},
#
# for sets of parameters :math:`\{\alpha_k\}` and :math:`\{\gamma_k\}`, which are optimized.
# If we run FALQON for :math:`p` steps, setting :math:`H_d = H_m`, and use the same cost Hamiltonian, we will end up with
# the following ansatz:
#
# .. math:: U_{\text{FALQON}} = e^{-i \Delta t \beta_p H_m} e^{-i \Delta t H_c} \cdots e^{-i \Delta t \beta_1 H_m} e^{-i \Delta t H_c}.
#
# Thus, our strategy is to initialize our QAOA parameters using the :math:`\beta_k` values that FALQON yields.
# More specifically, we set :math:`\alpha_k = \Delta t \beta_k` and :math:`\gamma_k = \Delta t`. We then optimize
# over these parameters. The goal is that these parameters provide QAOA a good place in the parameter space to
# begin its optimization.
#
# Using the code from earlier in the demonstration, we can easily prototype this process. To illustrate the power of
# this new technique, we attempt to solve MaxClique on a slightly more complicated graph:
new_edges = [(0, 1), (1, 2), (2, 0), (2, 3), (1, 4), (4, 5), (5, 2), (0, 6)]
new_graph = nx.Graph(new_edges)
nx.draw(new_graph, with_labels=True, node_color="#e377c2")
######################################################################
# We can now use the PennyLane QAOA module to create a QAOA circuit corresponding to the MaxClique problem. For this
# demonstration, we set the depth to :math:`5`:
depth = 5
dev = qml.device("default.qubit", wires=new_graph.nodes)
# Creates the cost and mixer Hamiltonians
cost_h, mixer_h = qaoa.max_clique(new_graph, constrained=False)
# Creates a layer of QAOA
def qaoa_layer(gamma, beta):
qaoa.cost_layer(gamma, cost_h)
qaoa.mixer_layer(beta, mixer_h)
# Creates the full QAOA circuit as an executable cost function
def qaoa_circuit(params, **kwargs):
for w in dev.wires:
qml.Hadamard(wires=w)
qml.layer(qaoa_layer, depth, params[0], params[1])
@qml.qnode(dev)
def qaoa_expval(params):
qaoa_circuit(params)
return qml.expval(cost_h)
######################################################################
# Now all we have to do is run FALQON for :math:`5` steps to get our initial QAOA parameters.
# We set :math:`\Delta t = 0.02`:
delta_t = 0.02
res, res_energy = max_clique_falqon(new_graph, depth-1, 0.0, delta_t, dev)
params = np.array([[delta_t for k in res], [delta_t * k for k in res]], requires_grad=True)
######################################################################
# Finally, we run our QAOA optimization procedure. We set the number of QAOA executions to :math:`40`:
#
steps = 40
optimizer = qml.GradientDescentOptimizer()
for s in range(steps):
params, cost = optimizer.step_and_cost(qaoa_expval, params)
print("Step {}, Cost = {}".format(s + 1, cost))
######################################################################
# To conclude, we can check how well FALQON/QAOA solved the optimization problem. We
# define a circuit which outputs the probabilities of measuring each bit string, and
# create a bar graph:
@qml.qnode(dev)
def prob_circuit(params):
qaoa_circuit(params)
return qml.probs(wires=dev.wires)
probs = prob_circuit(params)
plt.bar(range(2**len(dev.wires)), probs)
plt.xlabel("Bit string")
plt.ylabel("Measurement Probability")
plt.show()
######################################################################
# The state :math:`|112\rangle = |1110000\rangle` occurs with highest probability.
# This corresponds to nodes :math:`0`, :math:`1`, and :math:`2` of the graph, which is
# the maximum clique! We have successfully combined FALQON and QAOA to solve a combinatorial
# optimization problem 🎉.
#
# References
# ----------
#
# Magann, A. B., Rudinger, K. M., Grace, M. D., & Sarovar, M. (2021). Feedback-based quantum optimization. arXiv preprint `arXiv:2103.08619 `__.