Quantum Goemans-Williamson Algorithm with the Hadamard Test and Approximate Amplitude Constraints

In quantum analogy to the objective function of Eq. 1, we wish to minimize $\langle W,X\rangle$ over the $n$-qubit density matrices $\rho$, which are positive semidefinite by definition. We generate these quantum ensembles $\rho$ from an initial density matrix $\rho_{0}$ such that $\rho=U_{V}\rho_{0}U_{V}^{\dagger}$, where $U_{V}$ is a variational quantum circuit. This yields the quantum objective function

The Hadamard Test (Fig. 1(c)) is a quantum computing subroutine for arbitrary $n$-qubit states $|\psi\rangle$ and $n$-qubit unitaries $U$ [38]. It allows the real or imaginary component of the $2^{n}$-state inner product $\langle\psi|U|\psi\rangle$ to be obtained by estimating only a single expectation value $\langle\sigma_{n+1}^{z}\rangle$, which is the $z$-axis Pauli spin on the $n{+}1$th (auxiliary) qubit. For example, to obtain the imaginary component of $\langle\psi|U|\psi\rangle$, we prepare the quantum state $|\psi\rangle\otimes\frac{1}{\sqrt{2}}(|0\rangle-i|1\rangle)$ and apply a controlled-$U$ from the $n{+}1$th qubit to $|\psi\rangle$, followed by a Hadamard gate on the $n{+}1$th qubit. This produces the state

Rather than estimate the $N\leq 2^{n}$ expectation values required to characterize $\rho$ and optimize the loss function of Eq. 3, our method encodes the $N$-dimensional objective matrix $W$ as the imaginary part of an $n$-qubit unitary $U_{W}=\exp(i\alpha W)$ (Fig. 1(a)). Here, the phase $\alpha$ is a constant scalar. $U_{W}$ is then conditioned on the $n{+}1$th (or auxilary) qubit as a controlled-unitary. We then use the Hadamard Test to calculate the objective term in the loss function

The intuition for this objective function is that, for sufficiently small $\alpha$, $\text{Im}[U_{W}]\approx\alpha W$. By restricting ourselves to quantum circuits with real-valued output states, we render the single expectation value $\langle\sigma_{n+1}^{z}\rangle_{W}$ proportional to the objective function of Eq. 3, which requires $N$ expectation values to estimate. In Sec. 4.1, we analytically prove that, for many optimization problems, $\alpha$ has a practical upper bound such that $\text{Im}[U_{W}]\approx\alpha W$ with a reasonably large $\alpha$, even for arbitrarily large $W$.

We now illustrate how $\text{Im}[U_{W}]$ can be a close approximation of $\alpha W$, including for optimization problems with an arbitrarily large number of variables $N$. For concreteness, we select the NP-complete problem MaxCut [39], and specifically focus on the corresponding NP-hard optimization problem [45]. This problem is of particular interest due to its favorable $0.878$-approximation ratio with semidefinite programming techniques, notably the Goemans-Williamson algorithm [6], for which we now derive an efficient quantum implementation. The Goemans-Williamson algorithm is also applicable to to numerous other optimization problems, such as MaxSat and Max Directed Cut [6].

For a MaxCut problem with $N$ vertices $v_{i}$, $v_{j}$, let $W$ be the matrix that encodes the up to $N(N-1)/2$ non-zero edge weights in its entries $W_{ij}$. As the vertices lack self-interaction, $W_{ii}=0$. The optimization problem is then defined as

As described by Eq. 3, we can transform the optimization portion of Eq. 9 by substituting the classical positive semidefinite matrix $X$ for the quantum density operator $\rho$. The solution to the SDP is then stored in $|\psi\rangle$, i.e., $v_{i}=\text{sign}(\psi_{i})$ (for more details, see Sec. 2.3). As detailed in Sec. 2.1, the evaluation of this objective function can be optimized by estimating a single expectation value with the Hadamard Test. Likewise, we now introduce a quantum alternative to the constraint $X_{ii}=1$ from Eq. 9. First, note that due to the orthonormality of quantum states, the exact quantum equivalent of Eq. 9 is

As graph weights are real-valued and symmetric (i.e., $W_{ij}=W_{ji}$), $W$ is Hermitian. We can thus use it as the generator of $U_{W}$ such that (Fig. 1(a))

Next, we note that enforcing the $M=N\leq 2^{n}$ amplitude constraints $\rho_{ii}=2^{-n},\hskip 2.84544pti\leq N$ would require the estimation of all $z$-axis Pauli strings of order $k\leq n$ (all Pauli strings with $k\leq n$ Pauli-$z$ operators) of the state $|\psi\rangle$. This would be a total of $N-1$ expectation values. As an alternative to this large overhead, HTAAC-QSDP proposes the use of Approximate Amplitude Constraints (Fig. 1(d)). For example, consider the set of $m=n(n-1)/2+n\sim n^{2}/2$ Pauli strings of length $k\leq 2$

We again emphasize that these $k\leq 2$ constraints only approximately enforce the SDP constraint $\rho_{ii}=2^{-n}$. Fully satisfying this constraint would require restricting Pauli-$z$ correlations between any subset of the $n$ qubits, such that no states of unequal amplitude magnitudes are permitted. Eq. 10 can be fully satisfied if we constrain $|\psi\rangle$ with all of the Pauli strings of length $k\leq n$. However, as there exist $n$ choose $k$ $z$-axis Pauli strings of order $k$, this requires estimating $\sum_{k=1}^{n}\left(\begin{smallmatrix}n\\ k\end{smallmatrix}\right)=2^{n}-1$ different expectation values and greatly decreases the efficiency of the algorithm. Sec. 3 details that, in practice for $800$-vertex graphs, competitive results are obtained using only $k\leq 2$ constraint terms (Fig. 2 and Table 2), and optimization performance is largely saturated with terms $k\leq 4$ (Fig. 4 and Table 3).

In order to explicitly see how the $m\sim n^{2}/2$ constraints of Eq. 13 largely enforce the constraint $\rho_{ii}=1/2^{n}$, let us take the example of a three-qubit state ($n=3$), which can encode up to eight vertices ($N=2^{n}=8$) using HTAAC-QSDP. Any real-valued $n=3$ state can be written generically as

Eq. 10 can also be systematically undermined by the unequal distribution of graph edges among the quantum states. For instance, the asymmetrically distributed edge-weights in skewed graphs (Fig. 5 left and Sec. 3). With such graphs, the minimization of the loss function can lead to outsized state populations for quantum states that encode high-degree (high edge-weight) vertices. Moreover, as the number of classical variables will not generally be a power of two, there will often be quantum states that are not encoded with a classical variable. For example and as detailed in Sec. 3, we use $n=10$ qubits ($2^{n}=1024$ states) to optimize the $800$-vertex ($N=800$) GSet graphs, such that the states $801$ to $1024$ are absent from the objective function. In such cases, the minimization of the loss function can lead to outsized state populations of quantum states that are present in the optimization function. In principle, these imbalances can be addressed by increasing the magnitude of the Pauli string amplitude constraints, but this is known to cause poor objective function convergence [46].

To redress this systematic skew, we add a population-balancing unitary $U_{P}$ (Fig. 1(b)), which is implemented on $|\psi\rangle$ via a second Hadamard Test (Sec. 2.1, Fig. 1(c)) and adds the loss function term $\langle\sigma_{n+1}\rangle_{P}$. Specifically, $U_{P}=\exp(i\beta P)$ where $P$ is some diagonal operator of edge weights $P_{ii}=-(P_{max}-\sum_{j}|\omega_{ij}|)$, where $\beta$ is an adjustable hyperparameter and $P_{max}=\text{max}_{i}(\sum_{j}|\omega_{ij}|)$ is the maximum magnitude of edge weights for any given vertex. $U_{P}$ works to balance the state populations by premiating the occupation of states that are lesser represented by or absent from the objective function $\langle\sigma_{n+1}\rangle_{W}$.

Combining both the efficient Hadamard Test objective function and the Approximate Amplitude Constraints, we can use simple gradient descent-based penalty methods [46] to find the solution. Specifically, we minimize the HTAAC-QSDP loss function for the Goemans-Williamson algorithm

At the end of our protocol, the SDP solution is encoded into $|\psi\rangle$. Like in other QSDP protocols, $|\psi\rangle$ may either be used to extract the full $N$-variable solution or for less computationally intensive analysis (i.e., to characterize the features of the solution or as an input state for further quantum protocols). For many SDPs, such as MaxCut, a good approximation of the solution (here, cut value) can also be extracted from the expectation values comprising the cost function.

To extract this approximate cut value, we note that Eq. 8 can be rewritten as

If the full solution $|\psi\rangle$ is desired, then full real-space tomography of $|\psi\rangle$ must be conducted by calculating the $N$ marginal distributions of all $k\leq n$ Pauli strings along the $z$ and $x$-axes [47], although partial and approximate methods could make valuable contributions [48, 49]. We now show that once $|\psi\rangle$ is determined, it suffices to assign the partition of each vertex as $v_{i}=\text{sign}(\psi_{i})$, or the sign of the state component $\psi_{i}$.

In classical semidefinite programming algorithms, such as the Goemans-Williamson algorithm [6], the optimal solution $X^{*}$ is factorized by Cholesky decomposition into the product $X^{*}=T^{\dagger}T$, where $T$ is an upper diagonal matrix. The sign of each vertex $v_{i}$ is then designated as

We define the quantum parallel by noting that as $\rho=|\psi\rangle\langle\psi|$, its Cholesky decomposition is simply the $2^{n}\times 2^{n}$ matrix that has the first row $\langle\psi|$ and and all other entries equal to zero. In this decomposition, Eq. 20 reduces to

We now demonstrate that our method has an approximately $O(1/t)$ error convergence rate. Assuming that our loss function $\mathcal{L}$ is a convex function, which is generally only locally and/or approximately satisfied, it has been established that it will converge as $O(1/t)$ for $t$ iterations if it is $L$-Lipschitz continuous [50].

While parameterized quantum circuits are known to be $L$-Lipschitz continuous [24, 34], we here, in the name or thoroughness, explicitly demonstrate this for our particular loss function. Specifically, a function $f(x)$ is $L$-Lipschitz continuous if

Starting with the objective function $\langle\sigma_{n+1}^{z}\rangle_{W}=\langle\psi|U_{W}|\psi\rangle$, note that

Fig. 6 (left) demonstrates this approximate $O(1/t)$ convergence for the cut value of G11. The approximate nature of this convergence stems from both the discontinuous rounding process and the nonconvexities of the optimization space.

As explained above, the Goemans-Williamson algorithm [6] can be applied to numerous other optimization algorithms, such as MaxSat and Max Directed Cut [6]. Moreover, HTAAC-QSDP can be adapted to accommodate the constraints of various other SDPs. The use of our method is particularly advantageous when the constraints of a problem are amenable to being expressed through a tractable number of Pauli strings, or when these constraints can be approximately enforced by such a set of strings. The precise mapping of constraints to limited sets of Pauli strings depends on the problem at hand and may require some engineering. We here provide a few such examples.

We now overview the implementation of the HTAAC-QSDP Goemans-Williamson algorithm with mixed quantum states, such as might occur on a noisy quantum device or by interacting the utilized qubits with a set of unmeasured qubits. While the formalism for measuring Pauli strings on such systems is well known, we demonstrate how the required Hadamard Test remains viable. We start with $\rho_{R}$, some mixed state equivalent of the variational state $|\psi\rangle$. Upon application of the controlled-$U_{W}$ (without loss of generality for $U_{P}$) conditioned on the $n{+}1$th qubit, we obtain the density matrix

This mixed state formalism of the HTAAC-QSDP Goemans-Williamson algorithm works not only in principle, but also in practice. Fig. 2 displays the best cut value of a mixed state formalism with otherwise equal parameters (dark blue). The mixed state was generated by adding an unmeasured qubit to the randomly parameterized quantum circuit, which was then traced over prior to minimization and cut classification. The higher rank states moderately improved the performance on most graphs, with a mean higher-rank SDP value of $0.96$ (G11), $1.01$ (G14), and $0.98$ (G20), compared to $0.94$, $1.00$, and $0.98$ in the rank-1 case.

All simulations are done using a one-dimensional ring qubit connectivity, such that each qubit has two neighbors and the $n$th qubit neighbors the $1$st qubit. The circuit ansatz of the simulations for graphs G11-G21 (G81) is $120$ ($900$) repetitions of two variationally parameterized $y$-axis rotations interleaved with CNOT gates, alternating between odd-even and even-odd qubit control. The TensorLy-Quantum simulator [62, 63] is used for graphs G11-G21¹¹1Source code at https://github.com/tensorly/quantum/blob/main/doc/source/examples/htaacqsdp.py., while a modified version of cuStateVec [64] is used for G81. Gradient descent was conducted an ADAM [65] optimizer, with learning rate $\eta=0.01$ ($\eta=0.005$) for graphs G11-G21 (G81), as well as hyperparameters $\beta_{1}=0.9$, and $\beta_{2}=0.999$.

The evolution angle $\alpha$ was set as $\alpha=0.01$ for all graphs. The values of $\beta$ used in this work were $\beta=1/1.2$ for the toroid graphs and $\beta=1/3$ for the skew binary and skew integer graphs. We did not employ a unitary $U_{P}$ for the G81 graph. $\beta$ values should be chosen such that $\beta<1$, as diagonal entries are always satisfiable (i.e., some population can always be placed on the state, lowering the loss function), in contrast to edge cuts, which are not (i.e., not every edge can be cut with any given partition for general graphs). $\beta$ values can be tuned on the device in the real time by monitoring the Pauli string constraints and choosing a $\beta$ that leads to largely satisfied Pauli constraints with relatively small coefficients $\lambda$, such that the convergence of the algorithm is not hindered by large constraints that outweigh the objective function or lead to unstable convergence. In this work, we set $\lambda\propto\alpha/m$, to keep the total influence of $m$ constraint terms in proportion to the objective term $\langle\sigma_{n+1}^{z}\rangle_{W}\approx\alpha W$. Specifically, for the $800$-vertex graphs, we choose $\lambda=100\alpha/m$ for the toroid and skew binary graphs and $\lambda=50\alpha/m$ for the skew integer graphs. For the G81 graph, $\lambda=2000\alpha/m$.

In this subsection, we derive Theorem LABEL:theorem.general_claims, which we here restate for completeness:

As discussed above, Theorem 2 can be understood in two ways: that $\alpha$ satisfies the approximation of Eq. 32 while remaining tractably large for SDPs of arbitrary $N$, as long as $N$ does not 1) grow slower than the total number of edges $e$, or 2) grow slower than the the cube of $\xi$. We again note that Theorem 2 should hold for the densest graph region if the edge density is assymetrically distributed, i.e., $\xi$ should be the average number of edges for the densest vertices. As the conditions of Theorem 2 hold for graphs that are not too dense, they are widely satisfiable as the majority of interesting and demonstrably difficult graphs for MaxCut are relatively sparse [39, 10, 40, 42, 43]. Many classes of graphs for which MaxCut is NP-hard satisfy Theorem 2 with tractably large $\alpha$, even for arbitrarily large $N$.

As an example, we consider non-planar graphs, for which optimization problems like MaxCut are typically NP-complete. While planar graphs can be solved in polynomial time [66], a graph is guaranteeably non-planar when $e>3N-6$, which reduces to $\xi\geq 3$ in the limit of large $N$ [67]²²2Other families of easy graphs are even more restrictive, such as graphs that lack a giant component. In the limit of large $N$, these graphs only occur in more than a negligible fraction of all possible graphs when $d\geq 1/N$ and thus $\xi\geq 1/2$ [68].. In accordance with Theorem 2, constant values of $\xi$ actually permit $\alpha$ to grow as $N^{1/2}$, while for constant $\alpha$ $\xi$ can grow as $N^{1/3}$, such that our approximation is valid for a wide variety of large-scale non-planar graphs. Indeed, most standard benchmarking graph sets have a small average number of edges per vertex, e.g., $\xi=3$ [40, 42, 43], as sparse edge-density is common among graphs with real-world applications. In fact, solving MaxCut with many classes of dense graphs (i.e., graphs with nearly all non-zero edges) is provably less challenging, and therefore less interesting, than with their relatively sparse counterparts [69].

We here sketch a brief proof of Theorem 2 for Erdös–Rényi random graphs [70] with edge weights $W_{ij}\sim\mathcal{U}_{[0,b]}$, where $\mathcal{U}_{[0,b]}$ is the uniform distribution on the interval $[0,b]$. The edge density of a graph is described as $d=e/E$, where $e$ is the number of non-zero edges $e$ and $E=N(N-1)/2$ is the number of total possible edges. We provide a detailed proof of this and other graph types in the Appendix A.

• •

  The Hadamard Test encoding is a good approximation when $U_{W}\propto i\alpha W$.

• •

  This is satisfied when $\frac{\alpha^{3}}{3!}|W^{3}|_{ij}\ll\frac{\alpha}{1!}|W|_{ij}$ for typical edges between vertices $i$,$j$.

• •

  The mean of the non-zero elements in $W$ is $\overline{W_{ij}}=b/2$ ³³3The mean value of all elements of $W$ is $\overline{W_{ij}}^{\prime}=db/2$, however the relevant comparison is between the elements of $W^{3}$ and the non-zero elements of $W$..

• •

  Elements $(W^{3})_{ij}$ are the sum of $\sim N^{2}$ terms $W_{ij}W_{jk}W_{kl}$, with expectation value $\overline{W_{ij}W_{jk}W_{kl}}=b^{3}d^{3}/8$. That is, the additive error between the Hadamard encoding terms and the matrix elements $W_{ij}$ scales as $b^{3}d^{3}/8$.

• •

  $\frac{\alpha^{3}}{3!}\overline{(W^{3})_{ij}}\ll\frac{\alpha}{1!}\overline{W_{ij}}\rightarrow\alpha^{2}\ll 24/(N^{2}b^{2}d^{3})$.

• •

  Substituting $d=2e/N(N-1)\approx 2e/N^{2}$ and $\xi=e/N$, we obtain Theorem 2.

While some optimization problems of $O(2^{n})$ variables may only be represented by graphs $W$ of $O(4^{n})$ distinct Pauli strings, we here illustrate that there are many interesting (indeed, NP-hard optimization problems) for which this is not the case. In particular, we focus on the MaxCut problem and discuss toroid and Erdos Renyi random graphs.

Toroid graphs, or graphs that can be embedded on a toroid such that none of the edges connecting vertices cross, have a regular, yet still three-dimensional (non-planar), topological structure [71]. While encoding difficult problems, these data sets can often be represented in just $\text{poly}(n)$ Pauli strings, as is the case with the 8-100 tourus family to which G11 pertains (Fig. 6, right, gray). What is more, the number of Pauli strings can be reduced even further by instead constructing an approximate operator $W_{\text{approx}}$ and permitting a small amount of error $\overline{|W_{\text{approx}}{-}W|}/\overline{|W|}\leq 0.015$ (Fig. 6, right, red). The population balancing graphs $P$ are a similar subset of structured graphs, whose diagonality renders them expressible with Pauli strings of only $z$-axis and identity gates.