Performant near-term quantum combinatorial optimization
Titus D. Morris Phillip C. Lotshaw

















table/search path=csv/,

pgfplots.groupplots

 






















    

ω̃
ω̂
p̃

∂
α
β̱
γ
Γ
δ̣
Δ
η
ϕ
λ
Λ
κ̨
μ
ν


π
σ
Σ
τ
ϵ
þθ
ζ
χ


    αβ̱γδ̣σΓŁΛω̃ω̂p̃∂𝒜ℬℱℋ̋𝒦𝒩Ø𝒪ℛ𝒲1/21/4a_+a_-Quantum Information Science Section, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USAQuantum Information Science Section, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USAThis manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).
        We present a variational quantum algorithm for solving combinatorial optimization problems with linear-depth circuits.  Our algorithm uses an ansatz composed of Hamiltonian generators designed to control each term in the target combinatorial function, along with parameter updates following a modified version of quantum imaginary time evolution. We evaluate this ansatz in numerical simulations that target solutions to the MAXCUT problem. The state evolution is shown to closely mimic imaginary time evolution, and its optimal-solution convergence is further improved using adaptive transformations of the classical Hamiltonian spectrum, while resources are minimized by pruning optimized gates that are close to the identity.  With these innovations, the algorithm consistently converges to optimal solutions, with interesting highly-entangled dynamics along the way. This performant and resource-minimal approach is a promising candidate for potential quantum computational advantages on near-term quantum computing hardware. 
        Performant near-term quantum combinatorial optimization
            Phillip C. Lotshaw
            April 29, 2024
        =======================================================
        
         Introduction– Combinatorial optimization is widely viewed as a scientific computing challenge where quantum computers may make a significant impact <cit.>. Research on this topic has been motivated by several encouraging results obtained for the quantum approximate optimization algorithm (QAOA) <cit.> as well as its variants <cit.>.  These include theoretical analyses <cit.>, which have shown high QAOA approximation ratios 
        for the SK model <cit.> and for Maximum-cut on large girth graphs <cit.>, as well as a variety of encouraging results on hardware <cit.> and numerical simulations <cit.>, such as a numerical scaling advantage for the low-autocorrelation binary sequence (LABS) problem <cit.>.  Modified QAOA-like ansätze have also determined ways to improve performance by introducing additional parameterized gates <cit.>, including gates with counter-diabatic Hamiltonians <cit.> or based on Hamiltonians selected to maximize the energy convergence rate <cit.>, as well as many others <cit.>. 
        
        Despite encouraging progress there are several drawbacks to current QAOA-based approaches.  The primary limitation is that deep circuits are required for high performance <cit.> while noise at these depths places severe limitations on theoretical and experimental performance <cit.>. Graph instances can even be designed so that any given QAOA depth p cannot outperform classical algorithms in the noiseless case <cit.>. One reason large depths are required is that a single QAOA layer only entangles qubits that are adjacent in a given problem graph; hence the algorithm can require significant numbers of layers to build correlations across all qubits <cit.>, and further layers to remove this entanglement and approach the optimal classical solution <cit.>.  It has been argued that entanglement actually serves as a barrier to QAOA performance due to the prohibitive depths required to remove it <cit.>.  This is concerning—entanglement should be helping our quantum algorithms, not hurting them. There have been approaches to heuristically modify the QAOA ansatz to reduce the depth or improve performance, such as  Refs. <cit.>.  While these and related approaches do give noteworthy empirical improvements, they often do not improve on the theoretical p →∞ guarantee of standard QAOA, and can also introduce their own problems in terms of more involved parameter optimizations, as well as losing QAOA's physical motivation in terms of an approximation to an optimal annealing algorithm <cit.>. It is desirable to find other physically motivated approaches that do not have these shortcomings.
        
        
        In this letter, we present a new algorithm for quantum combinatorial optimization that overcomes several issues of previous approaches such as QAOA. Building on recent work on counter-diabatic inspired ansätze <cit.>, we show that our algorithm is guaranteed to be able to prepare the optimal solution at linear depth given suitable parameters.  We use parameters that are optimized through a systematic procedure, based on variational quantum imaginary time evolution and with an energy-dependent nonlinear objective that isolates the optimal solutions during training.  Numerical simulations show that the final states at different training iterations transition from approximate product states, to volume-law entangled states, and to final states confined to the optimal solution subspace. Our ansatz accesses these highly entangled states during training of short-depth circuits, rather than in deep circuits, as in QAOA. 
        Furthermore, we find that a small number of iterations succeed in identifying optimal solutions to a wide variety of problems we consider. The high-performance, low-depth, and modest training costs distinguish our algorithm as a prime candidate for successful near-term quantum optimization. 
        
         Combinatorial optimization– The goal of a combinatorial optimization problem is to minimize a cost function C( z) with an N-bit argument z = (z_1, …, z_N).  These problems can often be formulated as quadratic unconstrained binary optimization (QUBO) with C( z) = ∑_i<j J_i,j z_i z_j + ∑_i h_i z_i. 
         For quantum computing these problems are mapped to a diagonal Hamiltonian <cit.>
            H_c = ∑_α P_α,    H_c| z⟩ = C( z)| z⟩
         
        
        where P_α are Pauli-Z strings and the optimal solution is the ground state of H_c. Here we focus on the standard benchmarking problem of weighted MAXCUT. Given a graph G=(V,E) with edge weights {w_i,j}, the goal is to bipartition the vertices such that the summed weight of edges with endpoints in different sets is maximized.  The optimal bipartition z_opt minimizes C( z) = ∑_(i,j) ∈ E w_ijz_i z_j with z_i ∈{1,-1} (physics convention, H_c = ∑_(i,j) ∈ E w_ijZ_i Z_j) or C'( z) = -∑_(i,j) ∈ E w_ij(z_i+z_j-2z_iz_j)/2 with z_i ∈{0,1} (computer science convention, H_c = -∑_(i,j) ∈ E w_ij(1-Z_i Z_j)/2). Performance is quantified by the approximation ratio
        
        
            AR = ⟨ C ⟩/C_opt
        
        
        where C_opt is the optimal cost value. The different conventions are offset by a constant ∑_(i,j)∈ E w_ij/2 which affects the approximation ratio.  Following previous work, we use the physics convention for the Sherrington-Kirkpatric model, and the computer science convention for other cases. 
        
        
         Variational Ansatz– A recent series of papers by Xi Chen and collaborators has considered QAOA-like ansatze that incorporate counterdiabatic terms to improve convergence <cit.>, recently culminating in a highly-parameterized ansatz for protein folding problems <cit.> that we build upon here.
        The main ingredients of said ansatz are ZY rotations, which directly rotate between the between the ± 1 sector of ZZ terms of the Hamiltonian.  For MAXCUT, their prescription <cit.> kept only ZY rotations corresponding to ZZ terms that also occur in H_c.  We found that optimization of that ansatz commonly struggled to converge, so we introduce a fully connected ansatz of the form
        
            |Ψ(θ)⟩ = ∏_j = {(r,q) : r<q}∈ |V|exp(iθ_jZ_rY_q)|+⟩^⊗
         
        where j indexes all unique pairs of vertices.  The ordering of the vertices is chosen by summed vertex weight as explained in S1.
        For problems such as MAXCUT consisting of only quadratic ZZ terms, the ground state sector will consist of one or more states locally equivalent (by X flips) to the N-qubit GHZ state.  It can be shown that (<ref>) can reach any such GHZ pair up to inconsequential relative phases, and is overdetermined for such a mapping. 
        Each ZY rotation is implementable on quantum hardware in terms of either two CNOTs or one two-qubit gate of the form exp(-iθ ZZ) as found on certain trapped ion architectures such as Quantinuum devices. The two-qubit gate count for this ansatz is meager: Defining N_p = N(N-1)/2, we need only 2N_p CNOT gates, or half that if using arbitrary ZZ rotations. In comparison, a single layer of QAOA requires a parameterized two-qubit rotation for each ZZ term in H_c, which is less expensive than our ansatze for the least connected graphs and equivalent for fully connected graphs; we will show the cost of our ansatz can be minimized in later sections.
        
         Parameter Optimization– System agnostic, classical parameter optimizers are known to suffer from common pitfalls like local minima <cit.> and vanishing gradients <cit.>.  To improve performance we turn to a novel implementation of quantum imaginary time (QITE).  Prior implementations of QITE relied on McLachlan's principle <cit.>, minimizing 
        
            ||(∂/∂τ+H_c-E_τ)|Ψ(θ⃗)⟩||
        
        with respect to angle variations. Variationally evolving in imaginary time on a quantum computer by Eq. <ref> requires N_p(N_p+1)/2 circuits per iteration, with most relying on controlled ansatz preparation <cit.>.   
        
        We instead derive an alternative of QITE that requires fewer function evaluations and no controlled ansatz preparation.  We notice that it is simpler instead to enforce that the expectation value of each P_α appearing in H_c should evolve according the imaginary time as, 
        
            ∂⟨ P_α⟩/∂τ   = ∑_j 2 ℜ(⟨Ψ(θ⃗)|P_α∂ | Ψ(θ⃗)⟩/∂θ_j)θ̇_̇j̇
               =-⟨Ψ(θ⃗)|{P_α,H_c-E_τ}|Ψ(θ⃗)⟩ .
        
        Performing these evaluations for each P_α in H_c, we arrive at the matrix equation,
        
            G·θ⃗ = D
        
        where 
        
            G_α,j = ℜ(⟨Ψ(θ⃗)|P_α∂ | Ψ(θ⃗)⟩/∂θ_j) 
            
                D_α = -1/2⟨Ψ(θ⃗)|{P_α,H_c-E_τ}|Ψ(θ⃗)⟩ .
        
        For Ising like Hamiltonians encountered in H_c, all P_α commute, so D can be evaluated in a single circuit. Further each column of G, corresponding to a given parameter derivative, requires only two circuits.  To see this note that variations of P_α with respect to a single Pauli rotation [from Eq. (<ref>)] have the form
        
            ⟨ P_α⟩ (θ_j+δ) = a_0+b_0cos(2δ)+2G_α,jsin(2δ)  .
        
        Thus to determine G in QUBO, we need only 2N_p circuit evaluations plus one circuit for energy.  Then equation (<ref>) can be inverted numerically to determine θ̇⃗̇ and parameters can be updated with a simple forward Euler step θ⃗→θ⃗+Δτθ̇⃗̇. Fig. <ref> demonstrates the success of our variational imaginary time (VAR-IT) approach for an 8-vertex fully connected MAXCUT problem with uniform random weights sampled from U(0,1) and the ansatz Eq. (<ref>).
        
         Improving Convergence– For large N graphs, using H_c directly generally imposes harsh requirements on step size as different parameters have rates of change that differ by orders of magnitude.  Additionally, for many families of MAXCUT systems of interest, the spectral gap becomes vanishingly small with large N, leading to slow convergence in imaginary-time-based methods.  To mitigate these difficulties, we leverage the manifestly diagonal form of H_c, introducing the sigmoidal transformation
        
            f(H_c) = (1+exp(-σ_0 (H_c-E_τ)/4σ_τ^2))^-1
        
        where E_τ, σ_τ, and σ_0 are the current energy, current standard deviation, and initial standard deviation in the |+⟩^⊗ N state. This function retains the same optimal cut, but stretches the spectrum of states within a standard deviation of the current energy, while pushing farther states towards energies of zero and one. As the imaginary time evolution decreases E_τ and σ_τ, even dense regions of excited states nearly degenerate with the ground state will approach a gap of 1 in f(H_c).   In the rest of work, Eq. (<ref>) is evaluated using f(H_c) and ⟨ f(H_c)⟩, instead of H_c and E_τ, respectively. 
        
        
        
        Numerical and finite shot considerations makes exactly evaluating derivatives in Eq. (<ref>) a poorly posed task, which manifests most strongly at the end of the imaginary time evolution. We therefor proceed by solving (<ref>) [with f(H_c) in place of H_c] until the lowest energy cut measured achieves half the measurements, discarding the effect of singular values λ_ν < 0.01λ_max in G, where λ_max is the largest singular value. Once this is achieved, we perform Jacobi updates based on a single parameter while holding all the others constant, according to <cit.>; performing this update once for all parameters brings the energy and norm on optimal solutions numerically to one, apart from a single Sherrington-Kirkpatric (SK) instance where VAR-IT converged to a near-degenerate excited state.
        
        
         Numerical Results– We explore three prototypical problem types, 3 regular graphs, random Newman-Watts-Strogatz (NWS) small-world graphs, and the SK spin glass model with more details in Supplemental Information.  These represent linear connectivity, non-trivial clustering type connectivity, and fully connected problem instances, respectively.  In Fig. <ref> we show the approximation ratio error produced by our VAR-IT algorithm and by Eq. (<ref>) using the ADAM optimizer, averaged over one-hundred 16-vertex instances of the three problem types.  VAR-IT always converges to the optimal energy, unlike the ADAM optimization.  In Supplemental Information we show the fraction of instances that obtain <1% norm on  the optimal energy after 100 iterations is increasing rapidly with size for the ADAM optimization of NWS and SK graphs;  it is clear that ADAM will suffer a decreasing, and likely vanishing, probability of obtaining the optimal energy as the system size increases, whereas there is no signature of similar failure in VAR-IT in the sizes explored here.  Additionally, VAR-IT converges in significantly fewer iterations than ADAM.  Also plotted are average approximation ratios for p=8 QAOA, from Ref. <cit.>, which compare unfavorably to both ADAM and VAR-IT. All of the results in Fig. <ref> can be improved by the “relaxed-rounding" postprocessing method of Ref. <cit.>, but still VAR-IT outperforms these leading alternatives, as shown in Supplemental Information.
        
        
        
        
        When the number of ZZ terms in H_c scales linearly in problem size, QAOA requires ∼ N rotation gates while our ansatz requires ∼ N^2 gates. To combat this shortcoming we explore omitting ZY rotations if angles are smaller than a threshold ϵ at any point in the optimization.  In Fig. <ref>, we apply this adaptive procedure to 100 3-regular graphs for ϵ = 0.05,0.10,0.15.  The maximum number of rotations exploited during each graphs solution in terms of equivalent QAOA layers, is averaged and plotted in the left panel.  For 3-regular graphs and a given ϵ, the circuit cost appears to scale linearly at 2-3 layers of QAOA, with a minor increase in the number of iterations until convergence relative to the full ansatz.  In Supplemental Information, iterations required to converge scale linearly in problem size for SK and NWS, similarly as seen for 3 Regular in Fig. <ref>(b).   
        
        
         Entanglement–  
        If our algorithm produces large amounts of entanglement then it suggests the algorithm cannot be simulated at scale and therefor has potential for a quantum advantage, while low amounts of entanglement would seem to preclude quantum advantage due to simulability by tensor networks or related methods <cit.>.
        
        We analyzed the entanglement entropy S produced by our ansatz across 100 instances of the Sherrington-Kirkpatrick model at various N. For each instance we chose a random bipartition into sets of N/2 qubits and computed the entanglement entropy of the bipartition for the state produced at each step in the training, with averages over the instances shown in Fig. <ref>(a) and colored bands denoting the standard error of the mean.  There is a transition from states with low entanglement at small numbers of iterations, where states are close to the initial product state, to highly entangled states midway through training, and finally to states with decreasing entropy as the training approaches optimal solutions; the curves terminate at the average iteration at which the states reach fidelity F=1/2 with the optimal solution subspace, at which point the algorithm switches to Jacobi iterations to isolate the optimal solutions. Fig. <ref>(b) shows the maximum entanglement entropy during these training iterations scales as a volume law S ≈ aN +b. The behavior is similar to what is observed in QAOA <cit.>, except with the crucial difference that we see this transition while  training a shallow circuit, rather than in the dynamics of a deep circuit. Harnessing this entanglement at low depth distinguishes our algorithm as a more practical approach to solving combinatorial problems with a potential advantage.  
        
         Conclusion– We developed an algorithm for quantum combinatorial optimization that makes several important advances over leading alternatives, in terms of optimization procedure and iterations, convergence behaviors, and circuit depth; widespread success of the algorithm is evident in simulations of diverse problems with up to 20 qubits.  Key features of our approach include counterdiabatic Hamiltonian generators capable of preparing optimal solutions, training through QITE with a nonlinear objective that emphasizes the ground state, and resource-efficient generation of volume-law entanglement during training.  We expect the low depth and minimal training requirements of our algorithm to enable quantum computers to solve combinatorial optimization problems with several tens of qubits in the immediate future, with extensions to larger sizes and scaling analyses beyond classical simulability as quantum computers continue to develop.  The low depth, minimal resource requirements, high entanglement, and superior empirical performance make our approach a leading contender for near-term quantum computational advantage in combinatorial optimization.
        
         Acknowledgements– We thank Joseph Wang, Eugene Dumitrescu, Jim Ostrowski, Ryan Bennink, George Siopsis, Rizwanul Alum, Noah Bauer, and Rebekah Herrman for discussions of variational quantum optimization algorithms.  We thank Maxime Dupont for discussion and sharing numerical results from Ref. <cit.>. T.D.M was supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Science Center and devised the ansatz, optimization algorithm, performed the simulations, and contributed to writing. P.C.L. was supported by the DARPA ONISQ program under award W911NF-20-20051 and contributed to design, analysis, writing, and entanglement computations. This research used resources of the Compute and Data Environment for Science (CADES) at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.
        
        
        
        
        
        apsrev4-1
        
        § SUPPLEMENTAL MATERIAL
        
        
         §.§ Ansatz
        
        In pursuit of a solution of the binary optimization problem, we perform an analysis of the QAOA ansatz,  
        
            |θ⃗,γ⃗⟩_p = ∏_j=0^pexp(-iγ_j H_m)exp(-iθ_j H_c)|+⟩^⊗ ,
        
        where H_m is the standard mixing Hamiltonian and H_c is the cost Hamiltonian with an eigenspectrum H_c | z⟩ = C( z)| z⟩ that contains the set of classical solution values 
        
            H_m = -∑_j ∈ VX_j 
            
            H_c = ∑_j,k ∈ Eω_jkZ_j Z_k = ∑_αλ_α P_α .
        
        
        In the limit as p →∞, and with appropriately chosen θ⃗ and γ⃗, the QAOA ansatz yields adiabatic evolution from the ground state |+⟩^⊗ n of H_m to the optimal-solution ground state | z_opt⟩ of H_c. At finite p that are implementable on current quantum computers, QAOA has been shown to be successful at approximating | z_opt⟩ by variationally optimizing the parameters p. Analyzing the QAOA/adiabatic ansatz by combining the two sequential applications of a single layer under its equivalent single exponential is illustrative.
        
        The leading order term after the sum of the driving and mixing Hamiltonians is their single commutator
        
            -γ_j*θ_j/2[H_m,H_c] = -iγ_j*θ_j/2∑_k,l ∈ Eω_kl(Y_k Z_l+Z_k Y_l)  .
        
        Previous work in accelerating convergence within QAOA found these are precisely the terms most important for accelerating convergence <cit.>.  
        
        Is it clear from the following relation
        
            exp(iπ ZY/2)
            {[ |00⟩±|11⟩; |01⟩±|10⟩ ]} = 
            {[ ∓|01⟩∓|10⟩;  |00⟩±|11⟩ ]}
        
        that ZY rotations exactly rotate norm between the ± 1 sector of ZZ terms of the Hamiltonian.  Based on these findings, we employ a parameterized ansatz of the form 
        
            |Ψ(θ)⟩ = ∏_j = {(r,q) : r<q}∈ |V|exp(iθ_jZ_rY_q)|+⟩^⊗
         
        where j indexes all unique pairs of vertices.  
        The ordering is important for robust convergence, so we introduce the quantity 
        
            ρ_j = ∑_k ∈ E|ω_jk|
        
        which represents the relative importance of each vertex in the overall H_c. Thus in this work, we order the ZY rotations in the ansatz by their corresponding size in ρ.  For a vertex set |V| indexed now by their order in ρ, the ansatz exploited in this work has an explicit form
        
            |Ψ(θ)⟩ = exp(iθ_jZ_N-1Y_N)…exp(iθ_jZ_0Y_2)exp(iθ_jZ_0Y_1)|+⟩^⊗ .
        
        
         §.§ Problem types
        
        The different problem types explored in this work are dictated by the ω_ij between vertices, commonly referred to as the adjacency matrix.  All graphs in this work are generated by the Python package NetworkX <cit.>.
        
         Sherrington-Kirkpatrick (SK)– The SK model in this work is defined as the fully connected graph type, with ω_ij chosen between ±1 chosen with equal probability. 
        
         Newman-Watts-Strogatz small-world graphs (NWS)– In order to create models of this type, vertices of a ring graph are connected with its k nearest neighbors.  In this work, k=4.  Afterwards, for each resulting edge (i,j), in the ring with connections to k nearest neighbors, add a new edge (i,l) with probability p, in this work p=0.5, to a randomly-chosen vertex l. Then ω_ij are assigned randomly in the interval (0,1].
        
         Random 3 Regular (3Reg)– This model arises from graphs where each vertex is connected to three other vertices randomly.  ω_ij are all assigned a uniform value of 1 in this case.
        
        
        
         §.§ Relaxed Rounding
        
        In <cit.>, the authors introduced a quantum analogue of the classical relax and round method.  They showed that it can easily be used to augment wavefunction methods like QAOA, and always improves the approximation ratio. It relies on the establishing the correlation matrix χ_ij, defined as 
        
            χ_ij = (δ_ij-1)⟨ Z_iZ_j⟩ .
        
        Then χ can be diagonalized classically to yield its eigenvectors {α_ν∈ℝ^N}.  These can be translated to their closest bitstrings by rounding {z_ν← sign(α_ν) ∈± 1^N}. Then the z_ν which minimizes the C(z) is chosen as the outcome of this quantum,relax, and round procedure.  We apply this procedure at each iteration of our optimization procedure and report the findings in Fig. <ref>.
        
        
        
        
        
         §.§ Further Numerical Results
        
        In numerical results section of this work, we presented the error of the approximation ratio when using our VAR-IT method for optimization of Eq. <ref>, and the ADAM optimizer.  
        
        
        In addition to the VAR-IT method always converging to an optimal solution, Fig. <ref> demonstrates that for the sizes explored in this work, the iterations required to reach convergence for all three problem types explored in this work appear to be scaling linearly in problem size.  Additionally, Fig. <ref> shows that an ansatz optimized with ADAM will yield increasing failures as problem size grows, while there does not appear to be any such signature with the VAR-IT optimization procedure.
