An Exclusive-Sum-of-Products Pipeline for QAOA

TL;DR

This paper introduces an ESOP-based method for encoding constraints in QAOA, improving approximation ratios for combinatorial problems like maximum independent set compared to traditional penalization techniques.

Contribution

It presents a novel ESOP encoding approach for constraints in QAOA, demonstrating improved performance on maximum independent set problems.

Findings

ESOP constraint formulations outperform standard penalization in approximation ratios.

Up to 30.3% improvement in approximation ratios with ESOP encoding.

Higher approximation ratios achieved after one QAOA layer on 64% of tested graphs.

Abstract

The quantum approximate optimization algorithm is commonly used to solve combinatorial optimization problems. While unconstrained problems map naturally into the algorithm, incorporating constraints typically requires penalizing constraint violations in the objective function. In this work, we propose an alternative approach that encodes constraints as Boolean expressions in exclusive-sum-of-products (ESOP) form before penalization. We test this method on the maximum independent set problem using graphs with 3 to 20 vertices and find that ESOP constraint formulations achieve higher approximation ratios than standard constraint penalization methods, with percent increases of up to 30.3%. Furthermore, ESOP constraint formulations result in higher approximation ratios than standard QAOA penalization approaches after one layer of the algorithm on approximately 64% of the tested graphs.

Tables

Table 1. Table 1: The table from [ 15 ] that shows how to write a Boolean expression f ​ ( x ) f(x) as a Hamiltonian H f H_{f} .

Equations

$$|\gamma ,\beta ⟩=U(B,{\beta }_p)U(C,{\gamma }_p)\cdots U(B,{\beta }_1)U(C,{\gamma }_1)|s⟩.$$

$${\gamma }^{\ast },{\beta }^{\ast }=\underset{\gamma ,\beta }{\mathrm{argmin}}⟨C⟩=\underset{\gamma ,\beta }{\mathrm{argmin}}⟨\gamma ,\beta |C|\gamma ,\beta ⟩.$$

$$\frac{⟨C⟩-C_{\max }}{C_{\min }-C_{\max }}$$

$$\underset{x\in \{0,1{\}}^{|V|}}{\max }$$

$$x_i+x_j\le 1,\forall (i,j)\in E.$$

$$\min -\sum _{i\in V}{x}_i+J\sum _{ij\in E}{x}_i{x}_j,$$

$$a\vee b$$

$$H_{f\oplus g}=H_f+H_g-2H_f{H}_g.$$

$$H_{f\oplus g}$$

$$=\frac{1}{4}(I-Z_i)(I-Z_j)+\frac{1}{4}(I+Z_j)(I-Z_k)-\frac{1}{8}(I-Z_i)(I-Z_j)(I+Z_j)(I-Z_k)$$

$$=\frac{1}{4}(I-Z_i)(I-Z_j)+\frac{1}{4}(I+Z_j)(I-Z_k)$$

$$=H_f+H_g$$

$${\ell }_{\sigma (u)}^i=\{\begin{array}{ll}{x}_u& \sigma (u)=0,\\ \overline{x_u}& \sigma (u)=1\end{array}.$$

$$H_{{⨁}_i{f}_i}=\sum _i{H}_{f_i}.$$

$$v_{\mathrm{IS}}(\overrightarrow{x_G})$$

$$\underset{ij\in E}{\overset{|E|}{\bigwedge }}\overline{x_i\wedge x_j},$$

$$v_{\mathrm{MIS}}(\overrightarrow{x_G})=\overline{\underset{ij\in E}{\overset{|E|}{\bigvee }}{x}_i\wedge x_j}.$$

$$\overline{v_{\mathrm{IS}}(\overrightarrow{x_G})}=\underset{ij\in E}{\bigvee }{x}_i\wedge x_j.$$

$$\overline{v_{\mathrm{IS}}(\overrightarrow{x_G})}$$

$$=e_0\vee e_2\vee \cdots \vee e_{|E|-1}.$$

$$\overline{v_{\mathrm{MIS}}(\overrightarrow{x_G})}$$

$$=(e_0\wedge \overline{e_1}\wedge \overline{e_2})\oplus (e_1\wedge \overline{e_2})\oplus (e_2)\vee e_3$$

$$\vee \cdots \vee e_{|E|-1}$$

$$=(e_0\wedge \overline{e_1}\wedge \overline{e_2}\wedge \overline{e_3})\oplus (e_1\wedge \overline{e_2}\wedge \overline{e_3})\oplus (e_2\wedge \overline{e_3})$$

$$\vee \cdots \vee e_{|E|-1}.$$

$$\overline{v_{\mathrm{IS}}(\overrightarrow{x_G})}$$

$$\oplus e_1\wedge \overline{e_2}\wedge \overline{e_3}\wedge \overline{e_4}\wedge \cdots \wedge \overline{e_{|E|-1}}$$

$$\oplus e_2\wedge \overline{e_3}\wedge \overline{e_4}\wedge \cdots \wedge \overline{e_{|E|-1}}$$

$$\oplus e_3\wedge \overline{e_4}\wedge \cdots \wedge \overline{e_{|E|-1}}$$

$$\oplus \cdots \oplus e_{|E|-2}\wedge \overline{e_{|E|-1}}\oplus e_{|E|-1}$$

$$\overline{v_{\mathrm{IS}}(\overrightarrow{x_G})}$$

$$=\left[\underset{k=0}{\overset{|E|-2}{⨁}}(x_{i_k}\wedge x_{j_k})\underset{r=k+1}{\overset{|E|-1}{\bigwedge }}\overline{(x_{i_r}\wedge x_{j_r})}\right]\oplus (x_{i_{|E|-1}}\wedge x_{j_{|E|-1}})$$

$$=\left[\underset{k=0}{\overset{|E|-2}{⨁}}(x_{i_k}\wedge x_{j_k})\underset{r=k+1}{\overset{|E|-1}{\bigwedge }}(\overline{x_{i_r}}\vee \overline{x_{j_r}})\right]\oplus (x_{i_{|E|-1}}\wedge x_{j_{|E|-1}})$$

$$=\left[\underset{k=0}{\overset{|E|-2}{⨁}}(x_{i_k}\wedge x_{j_k})\underset{r=k+1}{\overset{|E|-1}{\bigwedge }}\left[(\overline{x_{i_r}}\wedge x_{j_r})\oplus \overline{x_{j_r}}\right]\right]\oplus (x_{i_{|E|-1}}\wedge x_{j_{|E|-1}}).$$

$$c_i=\underset{j=1}{\overset{k}{\bigwedge }}{x}_j\to H_{c_i}=\frac{1}{2^k}\prod _j(-1{)}^{j+1}(I-Z_j).$$

$$H_c:=H_{{\oplus }_i{c}_i}=\sum _i{H}_{c_i}.$$

$$H_C=-H_{\mathrm{MAX}}-2|V|H_c=-H_{\mathrm{MAX}}-2|V|\sum _i{H}_{c_i}.$$

$$\overline{v_{\mathrm{IS}}(\overrightarrow{x_{P_4}})}$$

Full text

1 Introduction

The quantum approximate optimization algorithm (QAOA) [12] is a variational quantum algorithm that has the potential to provide a speed-up over classical algorithms for solving combinatorial optimization problems [26, 39, 42]. Recent QAOA research has resulted in numerous variations of the algorithm that seek to either decrease the circuit depth or problem size [32, 48], allow for more expressible circuits [17, 44, 46, 43], or reduce the space of feasible solutions [13, 45, 41].

While unconstrained combinatorial optimization problems such as MaxCut can easily be mapped into QAOA [28], solving constrained combinatorial optimization problems, such as the Maximum Independent Set (MIS) problem, with the QAOA is more challenging. One of the main methods for handling constraints is to penalize them in the objective function, turning the problem into an unconstrained one that can be naturally mapped to QAOA [18, 21]. These approaches can require ancillary qubits and can involve squaring sums of constraints, leading to products of several variables. This penalization can lead to complex cost Hamiltonians whose optimization landscape is difficult to navigate due to barren plateaus and local minima [19, 20, 3]. Another method of handling constraints is to encode the constraints directly into the initial state of the QAOA circuit, reducing the search space to only the feasible states [6, 44, 13, 45, 41, 14]. The drawback to this method is the difficulty in creating an efficient (in circuit depth or number of two-qubit gates) state preparation circuit to encode the feasible states. For single constraints, such as linear equality constraints [7, 44], the state preparation can be efficient, however it is not clear if this is easy for general sets of constraints, specifically constraints that share large amounts of variables. Other approaches for constrained optimization for QAOA try to avoid penalization altogether [4, 10].

The foundational work by Hadfield et al. [16] first established systematic methods for encoding Boolean constraints into phase oracles and Hamiltonians for quantum algorithms. Building on this foundation, recent approaches such as Boolean-Hamiltonian-Transformation QAOA (BHT-QAOA) [2] have extended these techniques to generate QAOA cost Hamiltonians by marking bit-strings that do not satisfy Boolean constraint expressions. In contrast to other penalization methods, BHT-QAOA formulates each constraint as an exclusive-sum-of-product (ESOP) Boolean expression. Interestingly, ESOP expressions have been used in quantum circuit compilation algorithms to reduce T gate count [25, 38] and in quantum walk assisted QAOA (QW-QAOA) to solve constrained combinatorial optimization problems [23], though the QW-QAOA approach updates the mixing Hamiltonian to search only over feasible solutions.

In this work, we first introduce background on QAOA, MIS, and ESOPs in Sec. 2. Then, we describe how to use ESOP expressions to encode constraints into QAOA, using MIS as a test problem, in Sec. 3. In Sec. 4, we compare the ESOP constraint approach to standard QAOA penalty approaches to solve MIS on over 3000 graphs with between 3 and 20 vertices and find that encoding constraints using the ESOP formulation results in up to a $30.3\%$ percent increase in approximation ratio over the standard QAOA constraint encoding approach. Finally, we conclude with a discussion in Sec. 5.

2 Background

In this section, we provide background information for QAOA, the MIS problem, SOPs, and ESOPs.

2.1 QAOA

QAOA is a variational quantum algorithm in which two unitaries, a cost unitary $U(C,\gamma)=e^{-iC\gamma}$ and a mixing unitary $U(B,\beta)=e^{-iB\beta}$, are applied to an initial state ${\left|s\right\rangle}$. The cost Hamiltonian, $C$, is a diagonal matrix that encodes the combinatorial optimization problem $f(x)$ being solved. The construction of $C$ for the MIS problem will be described in more detail later. The mixer Hamiltonian $B$ allows for probability amplitude to be transferred between the computational basis states. In order for QAOA to converge to the optimal solution to the combinatorial optimization problem, the initial state ${\left|s\right\rangle}$ should be an eigenstate of $B$. $B$ is often chosen to be the transverse field mixer, $B=\sum_{i}\sigma_{i}^{x}$, and ${\left|s\right\rangle}={\left|+\right\rangle}^{\otimes n}$. The final state after $p$ iterations of the algorithm is

[TABLE]

The real-valued parameters $\beta$ and $\gamma$, often called angles, are selected to minimize the expected value of the cost Hamiltonian $C$,

[TABLE]

Classical optimization routines such as BFGS [34] and COBYLA [33] are commonly used to find the angles.

A common measure of success for QAOA is the approximation ratio (AR), which for constrained optimization, is defined as

[TABLE]

where $C_{\mathrm{max}}$ ($C_{\mathrm{min}}$) is the largest (smallest) eigenvalue of $C$. The AR is a real number between 0 and 1, with values closer to 1 indicating better solutions to the combinatorial optimization problem.

2.2 MIS

The goal of the MIS problem on a simple graph $G=(V,E)$ is to find the largest subset of vertices $I\subset V$ such that for all $i,j\in I$, $ij\notin E$. For an example, see Fig. 1. We can formulate the MIS problem as an integer program

[TABLE]

The MIS problem is commonly written as a quadratic unconstrained binary optimization (QUBO) problem of the form [47, 31, 11, 8]

[TABLE]

which is mapped to $H_{C}$ using the substitution $x_{i}\rightarrow\frac{I-Z_{i}}{2}$ [15]. Note that this formulation recasts MIS as a minimization problem by negating the objective function. Throughout this work, we shall refer to this encoding as the “standard QAOA constraint encoding" method.

2.3 SOPs and ESOPs

Boolean expressions are commonly written in a SOP formulation [1, 5, 22], in which products of Boolean variables (ANDs) are added together (ORs). Any Boolean expression can be converted to an equivalent ESOP, which is a Boolean expression consisting solely of exclusive sums of products of literals that can be negated. One can convert a Boolean expression to an ESOP expression using the identity

[TABLE]

ESOP formulations of Boolean expressions have applications in machine learning [30] and Toffoli gate synthesis [37]. While converting Boolean expressions to ESOP expressions is hard, there exist parallel algorithms for converting an arbitrary Boolean expression to an ESOP formulation [29]. Furthermore, there exist efficient algorithms for minimizing ESOP expressions [40]. The author of [15] introduced rules for converting Boolean expressions into sums and products of Pauli-$Z$ operators, which are summarized in Table 1.

Another identity within [15] that will be commonly used in this work is that for two Boolean expressions $f$ and $g$ represented by Hamiltonians $H_{f}$ and $H_{g}$ respectively, then

[TABLE]

Note in particular that for a Boolean expression of the form $(x_{i}\wedge x_{j})\oplus(\overline{x_{j}}\wedge x_{k})$, where $f=(x_{i}\wedge x_{j})$ and $g=(\overline{x_{j}}\wedge x_{k})$, then $H_{f}=\frac{1}{4}(I-Z_{i})(I-Z_{j})$ and $H_{g}=\frac{1}{4}(I+Z_{j})(I-Z_{k})$. Substituting these into Eq. 1 yields

[TABLE]

as $(I-Z_{j})(I+Z_{j})=0$. Thus, one can eliminate higher-order terms from the Hamiltonian formulation of $f\oplus g$ if there exists a variable of the form $x_{j}$ in one expression and $\overline{x_{j}}$ in the other. In general, suppose we have a set of Boolean expressions $\{f_{i}\}_{i}$ where $f_{i}=\bigwedge_{u\in V_{i}}\ell_{\sigma(u)}(x_{u})^{i}$, $V_{i}\subset V$, $\sigma(u)\in\{0,1\}$, and

[TABLE]

If $|\{u\in V:\exists i,j\text{ s.t. }\ell_{\sigma(u)}^{i}(x_{u})\neq\ell_{\sigma(u)}^{j}(x_{u})\}|>1$, that is, if there is at least one vertex for which $x_{u}$ is in $f_{i}$ and $\overline{x_{u}}$ is in $f_{j}$, then

[TABLE]

3 ESOP pipeline

The following sections describe the framework for generating the QAOA cost Hamiltonian, $H_{C}$, for the MIS problem. In this framework, we first convert each constraint in the constrained optimization problem to a Boolean expression. Then, we rewrite the Boolean expression into an ESOP formulation. The ESOP formulation is then multiplied by a large penalty parameter and added to the objective function of the optimization problem. We shall refer to this as the "ESOP pipeline".The objective function is then converted to a cost Hamiltonian, $H_{C}$, via the rules from [15]. Once $H_{C}$ has been generated, we then use it in QAOA with the standard Pauli-$X$ mixer. To automate and test the procedure, we developed a Python library (found at [9]) which takes a Networkx graph as input and produces a QAOA circuit for the MIS problem on said graph. These circuits are then simulated using the Qiskit Aer simulator [35] in Qiskit version 1.3.2 [36]. Classical optimization of angles $\gamma$ and $\beta$ is employed using COBYLA [33].

Note that the methods in this section apply to any constrained optimization problem, however we use MIS as an illustrative example. The rest of this section will apply the ESOP pipeline to MIS problems.

3.1 MIS Boolean encoding

The first step of the ESOP pipeline for a given graph $G$ is to generate a Boolean function $v_{\mathrm{IS}}(\overrightarrow{x_{G}}):\{0,1\}^{n}\rightarrow\{0,1\}$ satisfying

[TABLE]

To do this, we will require a Boolean function that compares incident vertices. The independence constraints for the MIS problem can be written as

[TABLE]

since if both $x_{i}$ and $x_{j}$ are 1, the expression is equal to [math], signaling that both $i$ and $j$ cannot be included in $\overrightarrow{x_{G}}$. Using De’Morgans law, this is equivalent to

[TABLE]

This formulation of the constraints can be refactored into an ESOP using algebraic techniques and then added to the objective function.

3.2 Converting Boolean constraints to ESOP

By inverting both sides of Eq. (3), notice that

[TABLE]

For simplicity, we rewrite each $x_{i}\wedge x_{j}$ as some $e_{k}$ for $0\leq k\leq|E|-1$. Then for all edges $e_{k}\in E$,

[TABLE]

Applying the identity $a\vee b=(a\wedge\neg{b})\oplus b$ three times yields

[TABLE]

It is clear that continuing this process until there are no more OR operations will yield an equivalent expression of $v_{\mathrm{IS}}(\overrightarrow{x_{G}})$, written only using the XOR, AND, NOT, operations:

[TABLE]

This expression can be written more succinctly as

[TABLE]

The above process transforms the individual products into even longer terms of conjunctions and negated conjunctions that can be simplified with De’Morgan’s law. See Sec. 3.4 for an example of the ESOP derived for the path on four vertices.

3.3 Transformation from Boolean expression to a Hamiltonian

Once the constraints have been converted into the ESOP formulation $\bigoplus_{i}c_{i}$, where $c_{i}=\bigwedge_{j}x_{j}$ with some $x_{j}$ being negated, we will convert each term into a Hamiltonian using Table 1. Specifically, each variable $x_{j}\rightarrow\frac{1}{2}(I-Z_{j})$ and each $\overline{x_{j}}\to\frac{1}{2}(I+Z_{j})$. Where we diverge from Table 1 is how we handle $\bigwedge_{j=1}^{k}x_{j}$. Similar to [2], we will multiply every second variable in the term $c_{i}$ by $-1$, resulting in

[TABLE]

Since there exist $x_{j_{r}}$ and $\overline{x_{j_{r}}}$ terms in the ESOP formulation, we can use Eq. 2 to write the constraint penalty Hamiltonian as

[TABLE]

With the MIS objective Hamiltonian $H_{\mathrm{MAX}}=\sum_{j\in V}\frac{1}{2}(I-Z_{j})$, the full cost Hamiltonian is given by

[TABLE]

Note that normally when minimizing an objective function, one adds the Hamiltonians representing constraints to the objective. However, since we derived $\overline{v_{\mathrm{IS}}(\overrightarrow{x_{G}})}$ above, we subtract it from the objective instead of adding it. Also note that we use $2|V|$ as a penalty in our simulations, but any sufficiently large penalty parameter suffices.

3.4 Example: $P_{4}$

Consider the path graph on four vertices, denoted $P_{4}$, in Fig. 1. We will apply the ESOP pipeline to this graph to find a cost Hamiltonian that encodes MIS on this problem. Then, we simulate QAOA solving this problem using the Qiskit Aer simulator. It is clear that there are three solutions to this problem: $\{x_{1},x_{2}\}$, $\{x_{2},x_{3}\}$, and $\{x_{3},x_{4}\}$. Furthermore, the objective of this problem is $\max\sum_{i=1}^{4}x_{i}$, or equivalently, $\min\;-\sum_{i=1}^{4}x_{i}$.

The ESOP that encodes the MIS constraints is

[TABLE]

We shall refer to $(\overline{x_{1}}\wedge x_{2}\wedge{x_{4}})$ as $c_{1}$, $(x_{1}\wedge\overline{x_{3}}\wedge{x_{4}})$ as $c_{2}$, and $(x_{1}\wedge x_{3})$ as $c_{3}$. Note that by Eq. 2, $H_{c_{1}\oplus c_{2}\oplus c_{3}}=H_{c_{1}}+H_{c_{2}}+H_{c_{3}}$. Thus,

[TABLE]

where

[TABLE]

Fig. 2 shows the counts when simulating QAOA with this cost Hamiltonian. Note that the optimal solutions of $\{x_{1},x_{2}\}$, $\{x_{2},x_{3}\}$, and $\{x_{3},x_{4}\}$ correspond to bitstrings $1100$, $0110$, and $0011$ respectively.

4 Results

We implement the ESOP pipeline for penalizing constraints with QAOA on over 3000 graphs: all 3-vertex (2 instances), 4-vertex (6 instances), 5-vertex (21 instances), and 6-vertex graphs graphs (112 instances), 700 randomly sampled 7-vertex graphs, 500 randomly generated graphs for all vertex counts with 8 vertices and between 11 and 14 vertices, and 50 randomly generated graphs for all vertex counts between 15 and 20 vertices. All tested graphs are connected and nonisomorphic. Note that graphs with less than 8 vertices were generated using the g6 files from [24]. We calculate the approximation ratio for each instance and then compare the results to the standard QAOA penalization technique described in Sec. 2.2.

The approximation ratios for tested graphs with between 3 and 14 vertices is found in Fig. 3, while the approximation ratios for graphs with between 15 and 20 vertices is found in Fig. 4. Note that when $p=1$, the ESOP constraint encoding approach results in higher average approximation ratios across each graph size except $n=3$, which is a dataset consisting of only two graphs. At $p=2$ and $p=3$, the ESOP constraint encoding method results in higher approximation ratios than QAOA for all size graphs except $n=7$. The average approximation ratios of both approaches are listed in Table 2, as well as the percent change between the standard MIS constraint encoding for QAOA and the ESOP constraint encoding. When $p=1$, the ESOP constraint formulations result in higher approximation ratios than standard QAOA penalization approaches on approximately 64% of the tested graphs.

5 Discussion

In this work, we develop an ESOP pipeline for encoding constrained optimization problems into QAOA. In this pipeline, we first write all constraints as Boolean expressions, and then write the Boolean expressions into an ESOP formulation, in contrast to other QAOA constraint encoding methods. Then, we multiply these terms by a large penalty parameter and add it to the objective function. Our experimental results demonstrate that rewriting constraints as ESOP statements and penalizing them in the QAOA cost Hamiltonian results in significantly higher approximation ratios than standard constraint penalization methods. Out of 54 test configurations spanning different graph sizes (3-20 vertices) and QAOA layer depths ($p=1,2,3$), ESOP-QAOA achieves higher average approximation ratios in 51 cases. The method shows particularly large approximation ratio improvements on smaller graphs, with notable increases of $30.3\%$ for $4$-vertex graphs at $p=2$, $28.7\%$ for $8$-vertex graphs at $p=1$, and $27.2\%$ for $12$-vertex graphs at $p=1$. Even for larger problem instances (15-20 vertices), the ESOP constraint encoding method shows positive increases in approximation ratios, suggesting scalability.

This work highlights the importance of problem encoding in the quantum approximate optimization algorithm, as the ESOP formulation of constraints implements the same logic as the usual Boolean implementation. The ESOP formulation appears to work well for MIS because of the existence of $x_{i}$ and $\overline{x_{i}}$ terms in the Boolean expression, which allows for constraint Hamiltonian simplifications. Future interesting research avenues include developing problem-based ESOP simplification methods, similar to the higher order terms canceling in MIS, for other constrained optimization problems. Other future avenues of research include identifying problem classes where ESOP constraints guarantee simpler implementations with shorter circuit depth and fewer ancilla qubits than standard penalty methods, and developing parameter transfer techniques [27] to further enhance optimization efficiency. Finally, the ESOP constraint formulations result in higher approximation ratios than standard constraint encoding methods on approximately 64% of the tested graphs. Characterizing these graphs would allow one to select the best encoding method for a particular MIS instance.

Code and Data Availability

The code and data for this research can be found at https://github.com/mxttbrunet/Quantum-Walk-Project.

{credits}

5.0.1 Acknowledgements

M.B., S.S., and R.H. acknowledge NSF CNS-2244512. A.W. and R.H. acknowledge NSF CCF 2210063.

5.0.2 \discintname

The authors have no competing interests to declare that are relevant to the content of this article.