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

Semidefinite programs are optimization methods with a wide array of applications, such as approximating difficult combinatorial problems. One such semidefinite program is the Goemans-Williamson algorithm, a popular integer relaxation technique. We introduce a variational quantum algorithm for the Goemans-Williamson algorithm that uses only $n{+}1$ qubits, a constant number of circuit preparations, and $\text{poly}(n)$ expectation values in order to approximately solve semidefinite programs with up to $N=2^n$ variables and $M \sim O(N)$ constraints. Efficient optimization is achieved by encoding the objective matrix as a properly parameterized unitary conditioned on an auxilary qubit, a technique known as the Hadamard Test. The Hadamard Test enables us to optimize the objective function by estimating only a single expectation value of the ancilla qubit, rather than separately estimating exponentially many expectation values. Similarly, we illustrate that the semidefinite programming constraints can be effectively enforced by implementing a second Hadamard Test, as well as imposing a polynomial number of Pauli string amplitude constraints. We demonstrate the effectiveness of our protocol by devising an efficient quantum implementation of the Goemans-Williamson algorithm for various NP-hard problems, including MaxCut. Our method exceeds the performance of analogous classical methods on a diverse subset of well-studied MaxCut problems from the GSet library.


Introduction
Semidefinite programming (SDP) is a variety of convex programming wherein the objective function is extremized over the set of symmetric positive semidefinite matrices S + [1]. Typically, an N -variable extremization problem is upgraded to an optimization over N vectors of length N , which form the semidefinite matrices of S + . A versatile technique, SDP can be used to approximately solve a variety of problems, including combinatorial optimization problems (e.g., NPhard problems, whose computational complexity grows exponentially in problem size) [2], and is heavily used in fields such as operations research, computer hardware design, and networking [3,4]. In many such cases, semidefinite programs (SDPs) are integer programming relaxations, meaning that the original objective function of integer variables is reformed as a function of continuous vector variables [5], such as the Goemans-Williamson algorithm [6]. This allows the SDP to explore a convex approximation of the problem. Although such solutions are only approximate, SDPs are useful because they can be efficiently solved with a variety of techniques. These include interior-point methods, which typically run in polynomial-time in the number of problem variables N and constraints M [7]. In recent years, more efficient versions of these classical methods have been developed [8,9].
An additional advantage of optimization with SDPs is that many have performance guarantees in the form of approximation ratios. Approximation ratios are a provable worst-case ratio between the value obtained by an approximation algorithm and the ground truth global optimum [10]. In short, SDPs represent an often favorable compromise between computational complexity The weight matrix W is used to generate the unitary U W , which is rotated about the angle α and implemented as a controlled-unitary conditioned on the n+1th-qubit (auxilary qubit). (1b) The population-balancing unitary U P is generated by the diagonal matrix P , which offsets the asymmetric edge weights on certain vertices in proportion to some constant β.
(1c) The Hadamard Test is used to efficiently evaluate the objective function and population balancing constraints. The n-qubit state |ψ⟩ = U V |0⟩ is prepared with a variational quantum circuit U V . Although U V can, in general, be made of any set of parameterized n-qubit quantum gates, in this work, we use the circuits described in Sec. 3.1. The n + 1th (auxilary) qubit is initialized as (|0⟩ − i|1⟩)/ √ 2. Subsequently, the Hadamard Test is carried out: U W or U P is implemented as a controlled-unitary conditioned on the auxilary qubit, which is then measured to compute ⟨σ n+1 ⟩ W,t = Im[⟨ψ|U W |ψ⟩] or ⟨σ n+1 ⟩ P,t = Im[⟨ψ|U P |ψ⟩]. (1d) The M = 2 n SDP amplitude constraints constraints are approximately enforced with only m ∼ poly(n) Pauli string constraints (Eq. 13). These are computed by collecting n-qubit Pauli-z measurements and using marginal statistics to estimate the m expectation values. and solution quality.
Despite the favorable scaling of classical SDPs, they still become intractable for high-dimensional problems. A variety of quantum SDP algorithms (QSDPs) that sample n-qubit Gibbs states to solve SDPs with up to N = 2 n variables and M ∼ O(N ) constraints have been devised [19][20][21][22][23] (Table 1), as have methods for approximately preparing Gibbs states with near-term variational quantum computers [25][26][27]. The former of these algorithms are based on the Arora-Kale method [28] and provide up to a quadratic speedup in N and M . However, they scale significantly poorer in terms of various algorithm parameters, such as accuracy, and are not suitable for nearterm quantum computers. Quantum interiorpoint methods have also been proposed [29,30], in close analogy to the leading family of classical techniques.
Variational methods have long played a role in Table 1: Comparison of common quantum methods for classical optimization. The number of potential variables N and constraints M are given in terms of qubits n. Whether or not the method provides guarantees on its solutions is discussed, as is its suitability for near-term quantum devices, (i.e., fewer than hundreds of qubits with limited error correction [11]). Our Hadmard Test objective function and Approximate Amplitude Constraint Quantum SDP (HTAAC-QSDP, Fig. 1) ensures SDP approximation ratios, is suitable for near-term variational quantum devices, and provides efficient objective function evaluation (via the Hadamard Test) and constraints (via a second Hadamard Test and O(n 2 ) Pauli string constraints).

Method
N , M Scaling Solution Guarantee Near-Term Devices Quantum Adiabatic [12][13][14] n If Infinitely Slow Sometimes Quantum Annealing [15,16] n No Yes QAOA [17] n Sometimes Yes Boson Sampling [18] n No Yes QSDPs [19][20][21][22][23] 2 n SDP Approx. Ratios No Variational QSDPs [24] 2 n SDP Approx. Ratios Yes HTAAC-QSDP (Ours) 2 n SDP Approx. Ratios Yes, poly(n) exp. vals./epoch quantum optimization protocols [31] (Table 1), such as adiabatic computation [12][13][14], annealing [15,16], the Quantum Approximate Optimization Algorithm (QAOA) [17], and Boson Sampling [18]. However, only recently have variational QSDPs been proposed [24,32]. Patel et al [24] addresses the same optimization problems as the quantum Arora-Kale and interior-point based methods, but instead uses variational quantum circuits, which are more feasible in the near-term. Like other SDPs, this method offers specific performance guarantees in the form of approximation ratios [10]. While exact methods are efficient for some SDPs, for worst-case problems (e.g., problems with a large number of constraints M , such as MaxCut) they may require the estimation of up to O(2 n ) observables per training epoch. Although it has been demonstrated that solving NP-hard optimization problems on variational quantum devices does not mitigate their exponential overhead, problems such as MaxCut may still retain APX hardness and are upper bounded by the same approximation ratio of classical methods [33]. Likewise, while parameterized quantum circuits can form Haar random 2designs that result in barren plateaus [34], numerous methods of avoiding, mitigating, or and perturbing these systems to effectuate a more trainable space have been developed [35][36][37].
Our Approach: We propose a new variational quantum algorithm for approximately solving QSDPs that uses Hadamard Test objective functions and Approximate Amplitude Constraints (HTAAC-QSDP, Fig. 1).
Theorem 1 HTAAC-QSDP for the Goemans-Williamson algorithm uses n+1 qubits, a constant number of quantum measurements, and O(poly(n)) classical calculations to approximate SDPs with up to N = 2 n variables.
The details of the HTAAC-QSDP implementation must be engineered to each SDP and, in this work, we focus on the Goemans-Williamson algorithm [6], with particular emphasis on its application to MaxCut. In some cases, our method is nearly an exponential reduction in required expectation values, e.g., for high-constraint problems such as MaxCut. As described in Sec. 2, we achieve this, in part, through a unitary objective function encoding with the Hadamard Test [38] ( Fig. 1a), which allows for the extremization of the entire N -dimensional objective by estimating only a single quantum expectation value (Fig. 1c). Our quantum Goemans-Williamson algorimth for MaxCut requires M = N ≤ 2 n amplitude constraints, which we effectively enforce with only 1) a constant number of quantum measurements from a second Hadamard Test (Fig. 1b) and 2) the estimation of a polynomial number of properly selected, commuting Pauli strings (Fig. 1d).
In Sec. 3, we demonstrate the success of the HTAAC-QSDP Goemans-Williamson algorithm (Algorithm 1) by approximating MaxCut [39] for large-scale graphs from the well-studied GSet graph library [40] (Fig. 5). In addition to satisfying the 0.878 MaxCut approximation ratio [6], HTAAC-QSDP achieves cut values that are commensurate with the leading gradient-based classical SDP solvers [41], implying that we reach op-tima that are very near the global optimum of these SDP objective functions.
Finally, in Sec. 4.1 we discuss the feasibility of the Hadamard encoding. For general SDPs, we establish an upper bound (Theorem 2) on the phase α of our Hadamard encoding, such that our technique is a high-quality approximation of the original SDP. The purpose of this upper bound is to demonstrate that tractably large values of the unitary phase α are permissible (i.e., α need not become arbitrarily small) for encoding a wide variety of large-scale problems. We discuss the known difficulty, and thus usefulness, of graph optimization problems under these conditions. Specifically: Theorem 2 Our approximate Hadamard Test objective function U W ∼ iαW (Sec. 2.1) holds for graphs with randomly distributed edges if where e is the number of non-zero edge weights and ξ is the average number of edges per vertex.
We can view the criteria of Theorem 2 in two ways: for SDPs with arbitrarily many variables N , the size of α can be kept reasonably large while the Hadamard Test objective function (see Sec. 2.1) remains valid as long as 1) N does not grow slower than the total number of edges e, or 2) N does not grow slower than the the cube of ξ. Both of these conditions hold for graphs that are not too dense, meaning that they are widely satisfiable. For instance, the majority of interesting and demonstrably difficult graphs for MaxCut are relatively sparse [10,39,40,42,43]. We note that for graphs where edge-density is unevenly distributed, Theorem 2 should hold for the densest region of the graph, i.e., ξ should be the average number of edges per vertex for the most highly connected vertices.

Efficient Quantum Semidefinite Programs
The standard form of an N -variable, Mconstraint SDP is [1,2] minimize X∈S + ⟨W, X⟩ where W is an N × N symmetric matrix that encodes the optimization problem and A µ (b µ ) are N × N symmetric matrices (scalars) that encode the problem constraints. ⟨A, B⟩ denotes the trace inner product In this section, we detail a method of efficient optimization of the above SDP objective and constraints using quantum methods (Fig. 1), specifically by implementing Hadamard Tests and imposing a polynomial number of Pauli constraints. We provide a concrete example in the form of the Goemans-Williamson [6] algorithm for MaxCut [39], as summarized in Algorithm 1.

The Hadamard Test as a Unitary Objective
In quantum analogy to the objective function of Eq. 1, we wish to minimize ⟨W, X⟩ over the n-qubit density matrices ρ, which are positive semidefinite by definition. We generate these quantum ensembles ρ from an initial density ma- where U V is a variational quantum circuit. This yields the quantum objective function Throughout most of this work, we consider pure states such that ρ 0 = |0⟩⟨0| and |ψ⟩ = U V |0⟩, although we detail the case of mixed quantum states in Sec. 2.6. In the case of pure states, Eq. 3 yields minimize ⟨ψ|W |ψ⟩.
The Hadamard Test (Fig. 1c) is a quantum computing subroutine for arbitrary n-qubit states |ψ⟩ and n-qubit unitaries U [38]. It allows the real or imaginary component of the 2 n -state inner product ⟨ψ|U |ψ⟩ to be obtained by estimating only a single expectation value ⟨σ z n+1 ⟩, which is the z-axis Pauli spin on the n+1th (auxiliary) qubit. For example, to obtain the imaginary component of ⟨ψ|U |ψ⟩, we prepare the quantum state |ψ⟩⊗ 1 √ 2 (|0⟩−i|1⟩) and apply a controlled-U from the n+1th qubit to |ψ⟩, followed by a Hadamard gate on the n+1th qubit. This produces the state upon which projective measurement yields Rather than estimate the N ≤ 2 n expectation values required to characterize ρ 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αW ) (Fig. 1a). Here, the phase α is a constant scalar. U W is then conditioned on the n+1th (or auxilary) qubit as a controlled-unitary. We then use the Hadamard Test to calculate the objective term in the loss function Pauli constraints compared to max(C SDP ), the best results of classical gradient-based SDPs (specifically, interior points methods) [41]. Our performance on the skewed binary and integer graphs narrowly exceeds that of the classical method (max(C SDP ) < C Q ), while the classical method narrowly outperforms our quantum method for the toroid graphs (max(C SDP ) > C Q ). Overall, the performance of our HTAAC-QSDP implementation and its classical counterpart are commensurate. HTAAC-QSDP exceeds the C Q /C MAX > 0.878 Max-Cut approximation ratio (red dashed line) for all graphs, where C MAX is the true MaxCut of the graph. In this work, we assume C MAX as the largest-known cuts of the GSet graphs, which were obtained from intensive and repeated heuristic searches [44]. The quantum mixed state implementation detailed in Sec. 2.6 is depicted in dark blue. It furnishes a higher rank solution and has marginally improved performance.
The intuition for this objective function is that, for sufficiently small α, Im[U W ] ≈ αW . By restricting ourselves to quantum circuits with realvalued output states, we render the single expectation value ⟨σ z n+1 ⟩ 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, α has a practical upper bound such that Im[U W ] ≈ αW with a reasonably large α, even for arbitrarily large W .

Quantum Goemans-Williamson Algorithm
We now illustrate how Im[U W ] can be a close approximation of αW , including for optimization problems with an arbitrarily large number of variables N . For concreteness, we select the NPcomplete 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, no-tably 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 which can be mapped to the classical SDP with M = N constraints 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 ρ. The solution to the SDP is then stored in |ψ⟩, i.e., v i = sign(ψ 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 This rescaling changes neither the effectiveness nor the guarantees of the semidefinite program, because the salient feature of the constraint is that all of the quantum states have the same amplitude magnitude |ψ i |, such that all of the vertices are of equal magnitude and none are disproportionately favored. The solutions ρ and X differ only by a constant factor, such that ρ = X/2 n . This yields the quantum MaxCut SDP minimize ⟨W, ρ⟩ 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. 1a) As W is real, the odd powers of l in Eq. 12 are the imaginary components. The condition that Note that, when this condition holds, ⟨σ z n+1 ⟩ W approximates ⟨W ⟩ with only vanishing error and a rescaling by α. In Sec. 4.1, we prove Theorem 2, demonstrating that this condition is achievable with a tractable α (i.e., α larger than some fixed finite value that is constant in problem size N ) for a wide variety of graphs.
Next, we note that enforcing the M = N ≤ 2 n amplitude constraints ρ ii = 2 −n , i ≤ N would require the estimation of all z-axis Pauli strings of order k ≤ n (all Pauli strings with k ≤ n Pauli-z operators) of the state |ψ⟩. 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. 1d). For example, consider the set of m = n(n − 1)/2 + n ∼ n 2 /2 Pauli strings of length k ≤ 2 where σ z a are the Pauli z-operators on the ath qubit. We can use these equalities as partial constraints for the n-qubit output state |ψ⟩. This set of m ∼ n 2 /2 constraints approximates the same restrictions as the set of M = N constraints of Eq. 11 by limiting quantum correlations, as these indicate subsets of states with unequal populations. That is, each such z-axis Pauli string is the difference of the total populations of two equal partitions of state space. To gain intuition, let us consider a two-qubit state |ψ⟩ = [ψ 00 , ψ 01 , ψ 10 , ψ 11 ] T and the k = 1 Pauli string σ z 0 ⊗ I. Using the constraints for Eq. 13, we enforce the equality |ψ 00 | 2 + |ψ 01 | 2 = |ψ 10 | 2 + |ψ 11 | 2 Table 2: MaxCut statistics for all 800-vertex graphs studied by the leading gradient-based classical SDP (interior points) method [41]. The highest known MaxCut values (C MAX , found by intensive heursitics [44]) are greater the highest results obtained by the classical method (max(C SDP )), but the approximation ratio max(C SDP )/C MAX > 0.878 is satisfied. The largest cut values of our HTAAC-QSDP Goemans-Williamson algorithm (max(C Q )) are comparable with max(C SDP ), as are the average results (mean(C Q )). which, while not fully enforcing the constraints of Eq. 11 (e.g., not enforcing that all states have equal populations), does enforce equal populations among a subset of states. This results in a lighter-weight and more flexible set of constraints.
We again emphasize that these k ≤ 2 constraints only approximately enforce the SDP constraint ρ 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 |ψ⟩ with all of the Pauli strings of length k ≤ n. However, as there exist n choose k z-axis Pauli strings of order k, this requires estimating n k=1 ( n k ) = 2 n − 1 different expectation values and greatly decreases the efficiency of the algorithm. Sec. 3 details that, in practice for 800vertex graphs, competitive results are obtained using only k ≤ 2 constraint terms ( Fig. 2 and Table 2), and optimization performance is largely saturated with terms k ≤ 4 ( Fig. 4 and Table 3).
In order to explicitly see how the m ∼ n 2 /2 constraints of Eq. 13 largely enforce the con- Figure 3: The cut value estimated by quantum observables ⟨C Q ⟩ obtained by HTAAC-QSDP vs the true SDP rounded cut value C Q for the G11 (toroid), G14 (binary), and G20 (integer) graphs. As with classical SDP methods, low loss function values are correlated with high cut values. The strong correlation between the observable estimated cut and the rounded true cut value not only illustrates the convergence of the HTAAC-QSDP Goemans-Williamson algorithm despite its approximate nature, it also demonstrates its ability to extract cut values of many variables from few measurements. straint ρ 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 |ψ⟩ = 1 r,s,p=0 ψ rsp |ψ rsp ⟩ and its constraints from Eq. 13 are Combined with the normalization constraint ⟨ψ|ψ⟩ = 1, the above system of equations nearly guarantees that Eq. 10 is fulfilled. However, it still permits a small subset of states that do not satisfy Eq. 10 due to three-qubit correlations, e.g., States with higher-order correlations such as |ψ * ⟩, which neither satisfy Eq. 10 nor are disallowed by Eq. 13, can be avoided by adding higher-order Pauli string constraints. For the above n = 3 example, we would add the k = n = 3 constraint ⟨σ z 1 σ z 2 σ z 3 , ρ⟩ = 0, which would disallow |ψ * ⟩ as ⟨σ z 1 σ z 2 σ z 3 , ρ * ⟩ = 1. 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 Figure 4: The effect of including higher-order HTAAC-QSDP Pauli string amplitude constraints in MaxCut optimization on the G11 (toroid), G14 (binary), and G20 (integer) graphs [40]. (Left) the performance of HTAAC-QSDP increases as higher-order Pauli strings are used to constrain state amplitudes. The algorithm's performance saturates with k ≈ 4, indicating that the benefits saturate with less than a polynomial number of Pauli string constraints (k ≤ n). As illustrated by this work (e.g., Fig. 2 and Table 2), k = 2 (m ≈ n 2 /2) is often sufficient for competitive SDP optimization. (Right, solid lines) the variance of state magnitude σ ρ = var(ρ ii ) = var(|ψ i | 2 ) vs the order k of Pauli strings constraints. As k increases, σ ρ decreases considerably, although this effect is largely saturated by k ≈ 4. (Right, dashed line) in the absence of competing dynamics (i.e., ⟨σ n+1 ⟩ W and ⟨σ n+1 ⟩ P ), the Pauli string constraints are fully enforced such that 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. 1b), which is implemented on |ψ⟩ via a second Hadamard Test (Sec. 2.1, Fig. 1c) and adds the loss function term ⟨σ n+1 ⟩ P . Specifically, U P = exp(iβP ) where P is some diagonal operator of edge weights P ii = −(P max − j |ω ij |), where β is an adjustable hyperparameter and P max = max i ( j |ω 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 ⟨σ n+1 ⟩ 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 each time step t by preparing a quantum state ρ t on a variational quantum computer. The scalar λ is the penalty hyperparameter. While for simplicity we have chosen a single, time-constant λ for all constraints, in principle each constraint j could be parameterized with a distinct λ j , each of which could also vary as a function of t. The number of quantum circuit preparations required to optimize our HTAAC-QSDP protocol is constant with respect to the number of qubits n (and thus to the maximum number of vertices N = 2 n ), as ⟨U W , ρ t ⟩ and ⟨U P , ρ t ⟩ each require only the Pauli-z measurement ⟨σ z n+1 ⟩ on the auxilary qubit, and the m ∼ poly(n) amplitude constraint terms can be calculated from a single set of n-qubit measurements on the state |ψ⟩. The classical complexity of each training step scales as just m ∼ poly(n), as one marginal expectation value is calculated from the |ψ⟩ measurements for each of the m constraints.

Retrieving SDP Solutions
At the end of our protocol, the SDP solution is encoded into |ψ⟩. Like in other QSDP protocols, |ψ⟩ may either be used to extract the full Nvariable 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 where W sum = j<i W ij . We note that this latter sum can be computed on a quantum device as where H is the Hadamard gate, as H ⊗n |0⟩ is the positive-phase and equal superposition of all states. Thus, W sum can be estimated with a Hadamard test where U W is applied to the input state, which we here denote ⟨σ n+1 ⟩ ⊗H W . Likewise, note that the second summation in Eq. 16 is given by in the case that the constraints are well-enforced. This allows us to estimate the cut value with ⟨σ n+1 ⟩ W , the Hadamard test using the variationally prepared |ψ⟩. Thus, the cost function can be estimated directly as If the full solution |ψ⟩ is desired, then full realspace tomography of |ψ⟩ must be conducted by calculating the N marginal distributions of all k ≤ 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 |ψ⟩ is determined, it suffices to assign the partition of each vertex as v i = sign(ψ i ), or the sign of the state component ψ 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 † T , where T is an upper diagonal matrix. The sign of each vertex v i is then designated as where t i are the column vectors of T and g is a length-N vector of normally distributed random variables g j ∼ N (0, 1). We define the quantum parallel by noting that as ρ = |ψ⟩⟨ψ|, its Cholesky decomposition is simply the 2 n × 2 n matrix that has the first row ⟨ψ| and and all other entries equal to zero. In this decomposition, Eq. 20 reduces to As MaxCut has Z 2 symmetry, the cut values are symmetric under inversion, or flipping the sign of all vertices. This makes the sign of the normally distributed g 0 irrelevant to the graph partitioning. Without loss of generality, we can therefore set g 0 = 1 and classify each vertex as v i = sign(ψ i ).

Error Convergence
We now demonstrate that our method has an approximately O(1/t) error convergence rate. Assuming that our loss function 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.
for all parameter inputs x and y. L will be L-Lipschitz continuous if each of its components are independently so. Starting with the objective function ⟨σ z n+1 ⟩ W = ⟨ψ|U W |ψ⟩, note that where V k is the generator of the unitary matrix parameterized by θ k . Eq. 23 is now composed of two terms, each comprised entirely of normal matrices and vectors, save for perhaps the Hermitian generators V k , with some extremal eigenvalue a k . Then, |∂⟨σ z n+1 ⟩ W /∂θ k | ≤ 2a k and |∇ k ⟨σ z n+1 ⟩ W − ∇ k ⟨σ z n+1 ⟩ W | ≤ 4a k . This proof doubles for the constraint terms, save that we replace the unitary matrix U W with some Pauli string observable. 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.

Extensions to Other SDPs
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.

Max and Min Bisection
As one example, consider the Min/Max Bisection problems [51]. Min/Max Bisection are particularly relevant to very-large-scale integration (VLSI) for integrated circuit design [52], a vital application area for large-scale SDPs.
The SDPs for estimating the Max Bisection problem has the standard form: The first constraint is equivalent to requiring that half of the variables of X be partitioned equally, hence the term "bisection". In analogy with Eq. 11, Eq. 24 can be written as minimize ⟨W, ρ⟩ subject to i,j ρ ij ≤ −N/2 and ρ ii = 2 −n , ∀i ≤ N.

(25)
The second of these two constraints can be enforced by the Pauli strings constraints of Eq. 13. For large N and assuming no systematic correlations between the ordering of the vertices, the first constraint can be ensured by adding any single Pauli string constraint where O x is any Pauli string of σ x operators. To see how Eq. 26 enforces the first constraint of Eq. 11, consider that any operator O x induces a bit-flip on a subset of qubits, such that each state ψ i is mapped to another state ψ i ′ . This means that ⟨O x ⟩ = ⟨ψ|O x |ψ⟩ is the sum of 2 n−1 Table 3: MaxCut statistics for the 800-vertex graphs G11 (toroid), G14 (skew binary), and G20 (skew integer) for different orders of Pauli string constraints k. We compare the best cut value max(C SDP ) produced by the leading classical method [41] compared to those produced by HTAAC-QSDP, with each entry providing the ratio max(C Q )/max(C SDP ) (mean(C Q )/max(C SDP )). With relatively few Pauli string constraints (k = 4), our method exceeds the performance of classical methods on all graphs studied. products 2ψ * i ψ i ′ , where for each i, |ψ i | ≈ 2 −n/2 , as enforced by the Pauli-Z constraints of Eq. 13. If the probability that a random state ψ i of |ψ⟩ is positive is p, then in the limit of large N and uncorrelated vertex assignment The above yields ⟨O x ⟩ = 0 iff p = 1/2, which would correspond to the equal partitioning of the vertices required by the Bisection problems. In the case of correlated vertex encodings, the average of several Pauli-X strings ⟨O x ⟩ can be considered. We note that Eq. 26 can be modified to enforce any partition ratio by solving for ⟨O x ⟩ (Eq. 27) with the desired p.

MaxSat
MaxSat problems are another branch of optimization tasks with constraints that focus on equally weighted vertices. In Max k-Sat problems, the number of logical boolean strings of length k are maximized over a given set of boolean variables v i [53]. For example, Max 2-Sat is given by [ where a ij and b ij are the coefficients of the problem. To convert this problem into an SDP, we optimize the objective function where W − ij = a ij − b ij . Note that Eq. 29 is equal to Eq. 9 and that the number of satisfied boolean strings can be extracted from an expectation value like Eq. 19, save that it is now paired

Correlation Matrix Calculation
Correlation matrices are key to a wide array of statistical applications and can be estimated with limited information using SDPs [54,55]. In particular, autocorrelation dictates that correlation matrices X have unit diagonals (X ii = 1, ∀i ≤ N ), much like MaxCut, Min/Max Bisection, and MaxSat, and can thus be addressed with the rescaled quantum version (ρ ii = 2 −n , Eq. 10) and approximated by the z-axis Pauli string constraints. Meanwhile, the extremization of certain correlations (e.g., maximize/minimize X ij ) and the enforcement of inequality and equality constraints (e.g., 0.2 < X ij < 0.4 or X ij = 0.4) can be enforced by additional constraints with either Pauli strings or the tomography of a select few state components.

HTAAC-QSDP with Mixed Quantum States
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 ρ R , some mixed state equivalent of the variational state |ψ⟩. Upon application of the controlled-U W (without loss of generality for U P ) conditioned on the n+1th qubit, we obtain the density matrix Upon applying a Hadamard gate on the n+1th qubit and measuring it along the z-axis, we obtain (31) which is proportional to Eq. 3 with the same coefficient α as prescribed by Eq. 7. This enables not only optimization, but also evaluation of the cut value estimation ⟨C Q ⟩ (Eq. 19), which is an accurate representation of the true cut value (Fig.  3). Alternatively, the element-wise rounded cut value calculated using Eq. 21 can still be utilized as an approximation to the traditional Goemans-Williamson rounding scheme. Although it would not result in the typical inner-product rounding, the sign of each vertex could still defined as the relative sign between each ψ i and ψ 0 . 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 higherrank 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.

Simulations and Results
The viability of our HTAAC-QSDP Goemans-Williamson algorithm is displayed in Fig. 2 and Table 2. We compare C Q , the cut values obtained by HTAAC-QSDP, to max(C SDP ), the best results obtained by the leading gradient-based classical method [41]. We study all of the 800vertex MaxCut problems explored in [41] (Table 2) in order to make an extensive comparison with the leading classical gradient-based interior point method. These graphs represent a broad sampling from the well-studied GSet graph library [40]. Graphs G11, G12, and G13 have vertices that are connected to nearest-neighbors on a toroid structure and ±1 weights (Fig 5,  right), while G14 and G15 (G20 and G21) have binary weights 0 and 1 (integer weights ±1) and randomly distributed edge density that is highly skewed towards the lower numbered vertices (Fig 5, left).
HTAAC-QSDP with k ≤ 2-Pauli term constraints exceeds the performance of its classical counterpart on skewed binary and skewed integer graphs, and falls narrowly short of classical performance on toroid graphs ( Fig. 2 and Table 2). The slight differences between the HTAAC-QSDP solution quality and those of the classical solver are typical of comparing different SDP solvers, which often differ slightly in their answers due to different numerical factors, including sparsity tolerance, rounding conventions (especially in the context of degenerate SDP solutions), and other hyperparameters [56]. All trajectories converge above the 0.878-approximation ratio C Q /C max (dashed red line) guaranteed by classical semidefinite programming, where C max is the highest known cut of each graph found by intensive, multi-shot, classical heuristics [44]. As SDPs are approximations of the optimization problem, the extremization of the loss function and the figure of merit (here, cut value) are highly correlated, particularly for well-enforced constraints. Fig. 3 demonstrates the strong correlation between the cut values estimated by quantum observables ⟨C Q ⟩ (Eq. 19) and the fully rounded SDP result C Q , indicating that the HTAAC-QSDP Goemans-Williamson algorithm measurement of few quantum observables closely approximates the rounded and composited cut values of all variables.
The addition of Pauli string amplitude constraints with k > 2 can better enforce Eq. 10, leading to higher-quality solutions to the SDPs. Fig. 4 and Table 3 demonstrate that increasing k produces moderately higher C Q values for the 800-vertex graphs, until the performance increases saturate k ≈ 4. Moreover, we note that at k values ≈ 4, HTAAC-QSDP outperforms the analogous classical algorithm for all graph types, indicating a well-conditioned solution. Likewise, the population variance (solid lines) σ ρ = var(|ψ i | 2 ) decreases substantially un- Figure 5: The structure of the G11 and G14/G20 GSet graphs, where non-zero edges between two vertices are marked as blue dots. Left) the edges of the toroid graphs (G11, G12, G13) follow a fixed connectivity, with edges extending between neighboring vertices on a torus structure. Right) the skewed graphs have connectivity that is drawn from a random distribution, with edges extending between arbitrary vertices (G14 and G15 binary weights 1 and 0, G20 and G21 integer weights ±1). The degrees of each vertex are disproportionately biased towards the vertices of lower index, with edge density decaying as vertex number increases. We compare with all 800-vertex graphs considered in the leading classical analog [41].
Note that, in the absence of the competing objective function (⟨σ z n+1 ⟩ W ) dynamics, all Pauliz correlations become restricted as k → n. This results in the complete constraint |ψ i | = N −1/2 , ∀i, such that σ ρ → 0 (black dashed line in Fig. 4, left). We can understand this behavior in the context of other penalty methods [57]. In particular, consider the penalty methods formulation q(t, λ) = f (x t ) + λ t g(x t ), where in our case x t are the time-dependent optimized parameters, f (x t ) = ⟨σ z n+1 ⟩ W,t , and g(x t ) represents the constraints. In the limits t → ∞ and λ t → ∞, it is known that x t → x, where x is the fully enforced solution of the hard constraints in some neighborhood of x 0 [58]. That is, the constraint-only dynamics (black dashed line of Fig. 4) represents the solution quality with respect to constraints in the λ t → ∞ limit.
The HTAAC-QSDP Goemans-Williamson algorithm was also tested against the G81 graph, a 20,000 vertex graph with a similar structure to G11-G13, and which is the largest MaxCut graph available that has been benchmarked against classical SDP methods. The average ratio between the HTAAC-QSDP obtained cut value for G81 and that of the classical counterpart was within 1% of the 800-vertex toroid graphs tested (ratio of 0.93). While the requisite circuit-depth for some quantum objective functions is known to grow linearly with the size of the circuit's Hilbert space [59], this has not been shown for objective functions of the form of Eq. 15. This is in agreement with our observation that, while Hilbert space required to solve the G81 graphs is 32 times larger than that for G11-G21, it's circuit-depth is < 8 times larger. This reduction in requisite circuitdepth could also be due to the influence of constraints, which may effectively limit the degrees of freedom of the quantum circuit, a phenomenon that has been referred to as a "Hamiltonian informed" model [59].
That being said, even if the worst-case exponential bounds for quantum optimization were saturated, that would not preclude them from having lesser overhead than classical SDPs. We note that classically solving an NP-hard problem of 2 n variables is exponentially hard with respect to the 2 n variables (that is, it scales as O(2 2 n )) and that, while these problems can be classically approximated using polynomial time and memory with respect to the number of variables, this still yields a classical complexity ∼ polynomial(2 n ). For instance, the leading classical techniques, interior point methods [56], have a memory cost that scales as O(N 4 ) ∼ O(2 4n ) and a computational cost that scales as O(N 2 ) ∼ O(2 2n ), making them in fact less efficient than the exponential bounds put on some worst-case quantum objective functions of O(2 n /2) = O(2 n−1 ) ∼ O(2 n ) memory cost and O(2 n−1 ) ∼ O(2 n ) computational cost, a quartic and quadratic reduction, respectively.
Even among the lighter-weight yet often less powerful first-order methods, the classical overhead is considerable. For example, the popular branch of projection-based methods has arithmetic and memory complexity of O(2 3n ) and O(2 2n ), respectively, for such 2 n variable problems [60]. Alternative algorithms, such as the Arora-Kale algorithm, can have time complexity as efficient asÕ(M ) [61], although for problems such as MaxCut this still reduces toÕ(2 n )). On a similar topic of scale, we note that G81 was optimized with constraints of orders k ≤ 4, which results in a similar ratio of approximate to exact constraints as that of smaller experiments. Figure 6: (Left) The convergence of the cut value for the G11 graph. The cut value approaches its maximum trajectory value as ∼ 1/t, with some fluctuations due to both the discontinuity of the rounding process and the nonconvexities of the optimization space. The error convergence is markedly faster than 1/ √ t. (Right) Although the general Hilbert space of quantum operators requires 4 n operators to span (black), in practice, many SDPs of interest, such as toroid graphs, can be encoded with a polynomial number of Pauli strings (gray). Moreover, many of these Pauli terms are relatively insignificant to the problem structure, such that the approximate matrix W approx can be encoded with little error using only a fraction of the already reduced operator set (red).

Simulation Details
All simulations are done using a one-dimensional ring qubit connectivity, such that each qubit has two neighbors and the nth qubit neighbors the 1st 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 1 , while a modified version of cuStateVec [64] is used for G81. Gradient descent was conducted an ADAM [65] optimizer, with learning rate η = 0.01 (η = 0.005) for graphs G11-G21 (G81), as well as hyperparameters β 1 = 0.9, and β 2 = 0.999.
The evolution angle α was set as α = 0.01 for all graphs. The values of β used in this work were β = 1/1.2 for the toroid graphs and β = 1/3 for the skew binary and skew integer graphs. We did not employ a unitary U P for the G81 graph. β values should be chosen such that β < 1, as diagonal entries are always satisfiable (i.e., some population can always be placed on the state, lowering the loss function), in con-trast to edge cuts, which are not (i.e., not every edge can be cut with any given partition for general graphs). β values can be tuned on the device in the real time by monitoring the Pauli string constraints and choosing a β that leads to largely satisfied Pauli constraints with relatively small coefficients λ, 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 λ ∝ α/m, to keep the total influence of m constraint terms in proportion to the objective term ⟨σ z n+1 ⟩ W ≈ αW . Specifically, for the 800vertex graphs, we choose λ = 100α/m for the toroid and skew binary graphs and λ = 50α/m for the skew integer graphs. For the G81 graph, λ = 2000α/m.

Theoretical Analysis of Hadamard Test Unitaries
This subsection addresses the implementation of the Hadamard tests found in this work, including the approximation of W by U W with a finite phase α, the construction of W for difficult problems, and the implementation of the prescribed Hadamard Tests using ZX-calculus.

Finite Phase α for Unitary Objective Function
In this subsection, we derive Theorem ??, which we here restate for completeness: where e is the number of non-zero edge weights and ξ is the average number of edges per vertex.
As discussed above, Theorem 2 can be understood in two ways: that α 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 ξ. We again note that Theorem 2 should hold for the densest graph region if the edge density is assymetrically distributed, i.e., ξ 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 [10,39,40,42,43]. Many classes of graphs for which MaxCut is NP-hard satisfy Theorem 2 with tractably large α, 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 ξ ≥ 3 in the limit of large N [67] 2 . In accordance with Theorem 2, constant values of ξ actually permit α to grow as N 1/2 , while for constant α ξ 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., ξ = 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 ∼ U [0,b] , where 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.
Proof sketch of Theorem 2: • The Hadamard Test encoding is a good approximation when U W ∝ iαW .
• This is satisfied when α 3 3! |W 3 | ij ≪ α 1! |W | ij for typical edges between vertices i,j. 2 Other 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 ≥ 1/N and thus ξ ≥ 1/2 [68].
• The mean of the non-zero elements in W is That is, the additive error between the Hadamard encoding terms and the matrix elements W ij scales as b 3 d 3 /8. .

Construction of W
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 threedimensional (non-planar), topological structure [71]. While encoding difficult problems, these data sets can often be represented in just 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 approx and permitting a small amount of error |W approx −W |/|W | ≤ 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.
Likewise, we can use similarly few terms to construct Erdös-Rényi random graphs, in which edges between any two vertices are equally likely and occur with probability p [70]. As Pauli strings are a spanning set, these same statistics are replicated when such strings are chosen randomly. Moreover, we note that each Pauli string adds O(2 n ) edges, such that the graph is rapidly populated.

Construction of Controlled Unitaries
The construction of unitary rotations U W and U P follows naturally from ZX calculus [72]. Specifically, Pauli Gadgets can be used to generalize unitary rotations from the qubit to which they are applied to multiple qubits through the use of O(n) CNOT gates, one on either size of each qubit and its rotated counterparts in a conjugated ladder scheme [73]. Moreover, rotations along distinct Pauli axes are achieved by conjugating these qubits with π/2 rotation gates along said axis. The gates are selected to match the terms of W or P and α is the phase applied.
As this method generalizes to all rotations, it can also be paired with the controlled gates required for the Hadamard Test. Specifically, the Pauli rotation gate applied to the auxiliaryadjacent qubit is fashioned as a controlled rotation, with the control conditioned on the auxiliary qubit. Moreover, the small values of α used in this work make the addition of multiple Pauli terms by Trotterization favorable, as the error of this technique is bounded by α 2 /2 times the spectral norm [74].

Conclusion
The efficient optimization of very large-scale SDPs on variational quantum devices has to the potential to revolutionize their use in operations, computer architecture, and networking applications. In this manuscript, we have introduced HTAAC-QSDP, which uses n + 1 qubits to approximate SDPs of up to N = 2 n variables and M ∼ O(N ) constraints by taking only a constant number of quantum measurements and a polynomial number of classical calculations per epoch. As we approximately encode the SDP objective function into a unitary operator, the Hadamard Test can be used to optimize arbitrarily large SDPs by estimating a constant number of expectation values. Likewise, we demonstrate that the constraints of many SDPs can also be efficiently enforced with approximate amplitude constraints.
Devising a quantum implementation of the Goemans-Williamson algorithm, we approxi-mately enforce the M = 2 n constraints with a population-balancing Hadamard Test and the estimation of as few as m ∼ n 2 /2 Pauli string expectation values. We demonstrate our method on a wide array of graphs from the GSet library [40], approaching and often exceeding the performance of the leading gradient-based classical SDP solver on all graphs [41]. Finally, we note that by increasing the order k of our Pauli string constraints, we improve the accuracy of our results, exceeding the classical performance on all graphs while still estimating only polynomiallymany expectation values.
The approximate amplitude constraints of HTAAC-QSDP make it particularly helpful for problems with a large number of constraints M . The benefits of using the Hadamard Test objective function depend on the original optimization problem. The optimization matrix of many NPhard problems can be encoded with controlledunitaries of polynomially-many Pauli terms, such that the Hadamard Test would be efficient to implement. While such cases could instead be optimized by directly estimating polynomially-many different non-commuting expectation values, use of the Hadamard Test circumvents the need to prepare an ensemble of output states for each Pauli term, eliminating these extra circuit preparations. Conversely, optimizing worst-case objective functions with the Hadamard Test would require controlled-unitaries with up to O(2 n ) Pauli terms. While exact implementation of these problems with HTAAC-QSDP on purely gatebased quantum computers would be inefficient, such objective functions could be engineered as the natural time-evolution of quantum devices with rich interactions (e.g., quantum simulators [75][76][77] in their future iterations), or by approximate means.
Due to the immense importance of SDPs in scientific and industrial optimization, as well as the ongoing efforts to generate effective quantum SDP methods that are often limited by poor scaling in key parameters such as accuracy and problem size, our work provides a variational alternative with tractable overhead. In particular, the largest SDPs solved via classical methods, which required over 500 teraFLOPs on nearly ten-thousand CPUs and GPUs [78], could be addressed by our method with just ∼ 20 qubits.
In future work, the techniques of this manuscript can be extended to additional families of SDPs. For instance, SDPs that extremize operator eigenvalues are a natural application for quantum circuits [79]. Similarly, variational quantum linear algebra techniques [80] can potentially be adapted to enforce the more general constraints ⟨A µ , X⟩ = b µ , ∀µ ≤ M of Eq. 1. In many cases, more general constraints are likewise satisfiable with the Pauli string constraints, as suggested in this work. For instance, when the number of requisite constraints M is much smaller than the number of variables N , or, as is the case with our quantum implementation of the Goemans-Williamson algorithm, by enforcing a relatively small subset of the constraints. (32) for typical edges between vertices i,j. By induction, the criterion in Eq. 32 also guarantees that odd (imaginary) powers >3 will likewise be smaller than the first order term, and are thus also negligible. While this condition can always be satisfied with an arbitrarily small α, we in practice require that α maintain some finite size to avoid unitary rotations with vanishingly small gate times τ ∝ α and imaginary components ⟨σ n+1 ⟩ W ∝ α. We now demonstrate that this criteria can be met for a wide array of graphs with NP-complete MaxCut optimization complexity.
First, we consider Erdös-Rényi random graphs [70] with elements W ij ∼ U [0,b] , which are uniformly distributed on the interval [0, b]. The graphs are said to have edge density d, which is the fraction of non-zero edges e over total possible edges E = N (N − 1)/2. Typical elements (W 3 ) ij are the sum of ∼ N 2 terms W ij W jk W kl , with expectation value W ij W jk W kl = b 3 d 3 /8, such that the matrix elements of W 3 have the expectation value (W 3 ) ij = N 2 b 3 d 3 /8. As the mean of the non-zero elements in W is W ij = b/2, the criterion of Eq. 32 becomes We can rewrite this criterion in terms the number of non-zero edges e by noting that graph density d scales as d = e/E, where E = N (N − 1)/2 ≈ N 2 /2 is the number of non-zero edges possible for an N vertex graph. Likewise, the average number of edges per vertex is then ξ = e/N , and Eq. 33 can be rewritten as For graphs where edge density d is not uniformly distributed, the above conditions should hold for the most densely connected vertices of the graph.
We briefly illustrate how our approximation holds for a few other classes of graphs. For instance, graphs with elements W ij ∼ U [−b,b] drawn from uniform distributions with both positive and negative components generally require α ranges that are even more permissible (i.e., can be even larger) than those of the positive case, with the criterion of Eq. 33 serving as a small lower-bound.
Similar proofs of implementability can also be done for graphs with normally distributed weights W ij ∼ N (µ, σ 2 ) of mean µ and variance σ 2 . For the case µ ̸ ≪ σ, (W 3 ) ij = N 2 µ 3 d 3 and α need only satisfy which requires the same permissive scaling between N and d (e or ξ) as the condition Eq. 32 (Eq. 33) for positive uniform distributions. Likewise, for normal distributions where σ ≫ µ, Eq. 35 with µ → σ would be a large upper bound.