How to define a QAOA cost function

Fire Opal's QAOA Solver's flexible input definition offers flexibility to enable you to solve any type of QAOA problem. One possible input is a NetworkX Graph object, which is automatically converted by Fire Opal to a cost function and mapped to qubits.

Alternatively, the solver accepts a user-defined cost function as input, which makes it possible to define and solve both constrained and unconstrained optimization problems. This guide only covers the use of a cost function as an input. For an example on using an input graph with the QAOA Solver, check out the tutorial on Solving MaxCut using Fire Opal.

This guide provides examples of a few types of optimization problems that can be solved using the QAOA Solver and demonstrates how to construct the cost function as an input. The following problem types are covered:

1. An unconstrained graph problem: MaxCut.
2. A constrained problem, using penalty terms to formulate the constraints: Minimum Vertex Cover (MVC).
3. A generalized unconstrained problem formulation: Quadratic Unconstrained Binary Optimization (QUBO).

The problems explained in ths guide are just a few examples; the QAOA Solver can be used to solve all sorts of optimization problems. For more information on building cost functions for optimization problems, this article provides detailed explanation.

Summary workflow

The following steps describe the general workflow for defining the QAOA problem and solving it using Fire Opal's solve_qaoa function.

1. Define the cost function to minimize

Start by building a graph with a given list of nodes and edges. Using binary variables $n \in {0, 1}$, create a problem-dependent cost function to be minimized.

2. Add penalty terms to the cost function (optional)

If there are constraints to the problem, they can be included as a penalty terms in the cost function.

3. Solve the problem with Fire Opal's QAOA Solver

Solve the problem using Fire Opal's QAOA Solver. The problem (represented by the SymPy cost function or NetworkX graph) and problem type (String) are passed as parameters, and the output Dict is stored in the result variable.

result = fireopal.solve_qaoa(problem, problem_type, ...)

Example: The MaxCut problem

Given some undirected graph $G(V, E)$ with a vertex set $V$ and an edge set $E$, a cut is defined as a partition of the nodes $V$ of a graph into two subsets. The value of a given cut is the sum of weights of all edges that connect two nodes that belongs to different subsets. The MaxCut problem is defined as follows: given some undirected graph $G$, find the cut with the maximal value.

To encode this problem on a quantum computer, a qubit is assigned to each graph node. The following code defines an example problem setup by created a random graph.

import numpy as np

# Given a graph with randomly weighted edges.
# To change the weights, change the seed to any integer.
_rng = np.random.default_rng(seed=99)
node_count = 7
edges = [(0, 1), (1, 3), (1, 2), (3, 4), (4, 5), (4, 6)]
weighted_edges = [(*edge, _rng.random()) for edge in edges]

To model this problem, binary variables can be used such that $n_i = 1$ if vertex $i$ is in one set and $n_j = 0$ if it is in the other set. We can determine if an edge $(i,j)$ is in the cut by evaluating the expression $n_i + n_j − 2n_i n_j$, which is equal to 1 if and only if the endpoints of the edge belong to different subsets. Otherwise $(n_i + n_j − 2n_i n_j)$ is equal to zero and the edge is not in the cut.

Therefore, the objective is to maximize the sum of the weights of edges in the cut, which can be formulated as follows: \begin{equation}\textbf{Maximize}\qquad y = \sum_{(i,j)\in E} \omega_{i,j}(n_i + n_j − 2n_i n_j)\end{equation} where $\omega_{i,j}$ is the weight of edge $(i,j)$.

However, since the Fire Opal's QAOA Solver takes the input cost function as a minimization problem, we can rewrite the objective as follows: \begin{equation}\textbf{Minimize}\qquad y = - \sum_{(i,j)\in E} \omega_{i,j}(n_i + n_j − 2n_i n_j)\end{equation} This is the standard unconstrained binary optimization model for the MaxCut problem.

from sympy import symbols

# Construct the cost function.
variables = symbols([f"n[{i}]" for i in range(node_count)])
cost_function = 0

for i, j, weight in weighted_edges:
    cost_function -= (
        variables[i] + variables[j] - variables[i] * variables[j]
    ) * weight

# Display the cost function.
cost_function

$\displaystyle 0.506030674548394 n[0] n[1] - 0.506030674548394 n[0] + 0.511915959692052 n[1] n[2] + 0.565091631784488 n[1] n[3] - 1.58303826602493 n[1] - 0.511915959692052 n[2] + 0.972186372136213 n[3] n[4] - 1.5372780039207 n[3] + 0.614903141969124 n[4] n[5] + 0.568283497894576 n[4] n[6] - 2.15537301199991 n[4] - 0.614903141969124 n[5] - 0.568283497894576 n[6]$

Example: The Minimum Vertex Cover (MVC) problem

Given an undirected graph $G(V, E)$ with a vertex set $V$ and an edge set $E$, a vertex cover is a subset of the vertices (nodes) such that each edge in the graph is connected to at least one vertex in the subset. The Minimum Vertex Cover problem seeks to find a cover with a minimum number of vertices in the subset.

The following code defines and sets up a random MVC problem.

# Given a graph with randomly weighted nodes.
# To change the weights, change the seed to any integer.
_rng = np.random.default_rng(seed=99)
node_count = 7
edges = [(0, 3), (0, 4), (1, 3), (2, 3), (2, 4), (2, 5), (4, 5), (4, 6), (5, 6)]
weighted_nodes = [(node, _rng.random()) for node in range(node_count)]

A standard optimization model for weighted MVC can be formulated as follows. First, a penalty must be added for any case where an edge is not connected to a vertex in the subset. Therefore, let $n_i = 1$ if vertex $i$ is in the cover (i.e., in the subset) and $n_i = 0$ otherwise. Second, the goal is to minimize the total number of vertices in the subset, which can be represented by the following function: \begin{equation}\textbf{Minimize}\qquad y = \sum_{i\in V} \omega_i n_i\end{equation}

# Construct the cost function.
variables = symbols([f"n[{i}]" for i in range(node_count)])
cost_function = 0

for i, weight in weighted_nodes:
    cost_function += variables[i] * weight

# Display the cost function.
cost_function

$\displaystyle 0.286786722869132 n[0] + 0.554511453248897 n[1] + 0.467523525550902 n[2] + 0.610058008019437 n[3] + 0.930442497679666 n[4] + 0.245885383640433 n[5] + 0.309438338643525 n[6]$

Now every edge in the graph should include at least one end point from the cover, which can be expressed as the inequality: \begin{equation}n_i + n_j \ge 1 \texttt{ for all } (i,j)\in E\end{equation}

Any case where an edge is not connected to the vertex of cover must be penalized. This can be represented in the cost function by adding a penalty of the form $P(1-n_i-n_j+n_i n_j)$ where $P$ is a positive penalty constant. Thus, an unconstrained alternative to the constrained inequality for weighted MVC is: \begin{equation}\textbf{Minimize}\qquad y = \sum_{i\in V}\omega_i n_i + P(\sum_{(i,j)\in E}(1 - n_i - n_j + n_i n_j))\end{equation}

# Add penalty term.
penalty_constant = 2
for i, j in edges:
    cost_function += penalty_constant * (
        1 - variables[i] - variables[j] + variables[i] * variables[j]
    )

# Display the cost function.
cost_function

$\displaystyle 2 n[0] n[3] + 2 n[0] n[4] - 3.71321327713087 n[0] + 2 n[1] n[3] - 1.4454885467511 n[1] + 2 n[2] n[3] + 2 n[2] n[4] + 2 n[2] n[5] - 5.5324764744491 n[2] - 5.38994199198056 n[3] + 2 n[4] n[5] + 2 n[4] n[6] - 7.06955750232033 n[4] + 2 n[5] n[6] - 5.75411461635957 n[5] - 3.69056166135648 n[6] + 18$

Example: The Quadratic Unconstrained Binary Optimization (QUBO)

All the cost functions constructed prior have a common feature: they all contain weighted sums on linear ($n_i$) or quadratic ($n_in_j$) terms of binary variables. Such cost functions can be written as a bilinear form of a vector $\vec{N}^T = (n_1,...,n_{|G|})$ and a symmetric matrix $Q$: \begin{equation}C = \vec{N}^TQ\vec{N} = \sum\limits_{i\in G} Q_{ii}n_i + 2\sum\limits_{i<j\in G} Q_{ij}n_in_j\end{equation} where $n_i^2 = n_i$. Hence, any symmetric matrix $Q$ defines a valid quadratic unconstrained binary optimization that may include penalty terms.

# Given any symmetric matrix Q, we can formulate the cost function as the input of QAOA solver.
Q = np.array(
    [
        [2, 1, -1, 0, 0],
        [1, 2, 0, -1, 0],
        [-1, 0, 3, -1, -1],
        [0, -1, -1, 3, -1],
        [0, 0, -1, -1, 2],
    ]
)

# Construct the cost function.
node_count, _ = Q.shape
variables = symbols([f"n[{i}]" for i in range(node_count)])
cost_function = 0

for i in range(node_count):
    for j in range(node_count):
        if i == j:
            cost_function += variables[i] * Q[i, j]
        if i < j:
            cost_function += variables[i] * variables[j] * 2 * Q[i, j]

# Display the cost function.
cost_function

$\displaystyle 2 n[0] n[1] - 2 n[0] n[2] + 2 n[0] - 2 n[1] n[3] + 2 n[1] - 2 n[2] n[3] - 2 n[2] n[4] + 3 n[2] - 2 n[3] n[4] + 3 n[3] + 2 n[4]$

from fireopal import print_package_versions

print_package_versions()

| Package   | Version |
| --------- | ------- |
| Python    | 3.11.3  |
| networkx  | 2.8.8   |
| numpy     | 1.23.5  |
| sympy     | 1.12    |
| fire-opal | 5.3.2   |