Quantum computational phase transition in combinatorial problems

Quantum Approximate Optimization algorithm (QAOA) aims to search for approximate solutions to discrete optimization problems with near-term quantum computers. As there are no algorithmic guarantee possible for QAOA to outperform classical computers, without a proof that bounded-error quantum polynomial time (BQP) ≠ nondeterministic polynomial time (NP), it is necessary to investigate the empirical advantages of QAOA. We identify a computational phase transition of QAOA when solving hard problems such as SAT—random instances are most difficult to train at a critical problem density. We connect the transition to the controllability and the complexity of QAOA circuits. Moreover, we find that the critical problem density in general deviates from the SAT-UNSAT phase transition, where the hardest instances for classical algorithms lies. Then, we show that the high problem density region, which limits QAOA’s performance in hard optimization problems (reachability deficits), is actually a good place to utilize QAOA: its approximation ratio has a much slower decay with the problem density, compared to classical approximate algorithms. Indeed, it is exactly in this region that quantum advantages of QAOA over classical approximate algorithms can be identified.


I. INTRODUCTION
Quantum Approximate Optimization algorithm (QAOA) [1], like all quantum algorithms, aims to utilize quantum hardwares to efficiently solve problems that are hard on classical computers.It is one of the candidates to achieve a quantum supremacy in the noisy intermediate-scale quantum (NISQ) era [2].So far, quantum supremacy has only been realized for random circuit sampling tasks [3,4].For complex but practical problems in the class of nondeterministic-polynomial (NP) time, no quantum advantage has been found, despite trials of QAOA on the Google Sycamore quantum processor [5].To search for quantum supremacy, it is crucial to first understand the difference between what is hard or easy for QAOA and for classical algorithms, which can be best explored via the combinatorial problem of Boolean satisfiability problem (SAT).
In a k-SAT instance, one asks whether multiple clauses, each involving k Boolean variables, can be satisfied simultaneously.Depending on the value of k, the worstcase hardness is drastically different-while 3-SAT is NPcomplete, 2-SAT can be efficiently solved in polynomial time (class P).In addition, the classical empirical hardness of random 3-SAT instances is known to have a computational phase transition [6][7][8][9][10] versus the problem density characterized by the clause-to-variable ratio.When the density is small (large), almost all instances are satisfied (unsatisfied) and easy to solve; while for density approaching the critical point of the SAT-UNSAT phase transition, the 3-SAT problem instances become the hardest to solve.
For quantum algorithms such as QAOA, the above phenomenon is largely unexplored.To begin with, as QAOA always implements SAT problems in their NPhard optimization versions (Max-SAT) [11], it is unclear whether the decision version's NP (k = 3) versus P (k = 2) contrast has any influence on QAOA's performance.It is also unclear how classical empirical hardness of 3-SAT connects to QAOA's performance on its optimization versions.Indeed, Ref. [12] does not find a big difference between QAOA's performance on Max-2-SAT and Max-3-SAT, and only finds QAOA's performance to worsen as the density increases-a phenomenon they call reachability deficits.
In this paper, we reveal a computational phase transition in the trainability of QAOA in solving the positive 1-in-k SAT problem.In terms of trainability characterized by gradient, the typical amplitude of gradient in training SAT problems achieves the minimum at a critical problem density ratio.In general, this quantum critical problem density deviates from the SAT-UNSAT phase transition [6][7][8][9][10], where the instances are hardest classically.We link this gradient transition to the controllability of the quantum systems evolving under the QAOA circuit [13][14][15][16] and the complexity of QAOA circuit [17][18][19][20].In terms of accuracy of the optimization versions of SAT, we find the QAOA's approximation ratio to be robust and decay slowly with the problem density.Moreover, despite the performance decay due to reachability deficits, it is precisely in the large problem density region where a quantum advantage can be identified, when comparing with classical approximate algorithms.In addition, the accuracy in solving Max-2-SAT is higher than that in solving Max-3-SAT, consistent with the P versus NP contrast in the decision version of the problems.Interestingly, for the decision version of the SAT problems, QAOA shows the worst performance at the arXiv:2109.13346v4[quant-ph] 15 Jun 2022 SAT-UNSAT transition, revealing a remnant of classical empirical hardness.Such remnant of classical empirical hardness is also confirmed in quantum adiabatic algorithms (QAA) [21][22][23].

A. Preliminary
Quantum Approximate Optimization algorithm.
To solve an optimization problem, QAOA encodes the cost function into the energy of a problem Hamiltonian H C , defined over spin-1/2 particles (qubits), and then seeks for an approximation of the ground state that encodes the solution to the optimization problem.An nqubit QAOA circuit implements dynamics governed by the problem Hamiltonian H C and a mixing Hamiltonian H B = n i=1 σ x i alternatively in each layer, where σ x i is the Pauli-X operator representing the transverse fields.The output state of a p-layer QAOA is therefore |ψ( γ, β) , where γ = (γ 1 , . . ., γ p ) and β = (β 1 , . . ., β p ) are variational parameters.The initial state is set to be a superposition of all possible spin configurations, |ψ (0, 0) = |+ ⊗n with To solve the problem, variational training is performed over the parameters γ, β to minimize the cost function The variational training terminates when the cost function stops to decrease significantly, and ideally leads to the optimal parameters γ * , β * = argmin γ, β C( γ, β).

SAT problems.
We will focus on two types of SAT problems, k-SAT problem (k ≥ 2) and the positive 1- uniformly randomly chosen from V .The k elements in each clause can be either positive or negative literal in k-SAT problem with equal probability, while only positive in 1-k-SAT + problems.The conjunctive normal form (CNF) of the SAT instance can be expressed as , where ' ' denotes AND and forces the CNF to be true only when all clauses are satisfied.In a k-SAT problem, each clause is true when at least one element in the clause is true; while in positive 1-in-k SAT, a clause c a {v aj } k j=1 is satisfied if and only if a single variable among {v aj } k j=1 is taken to be true.The (decision version of) SAT problem asks whether F (V ) can be satisfied with an assignment of variables V , while the optimization version-Max-SAT-aims to find an assignment of variables V to minimize the number of clause violations.With the increase of clause-to-variable ratio m/n, it becomes harder to satisfy a random SAT instance, and there exists a phase transition of SAT probability across a critical ratio, m/n = 1 for 2-SAT [24] and m/n = 4.26 for 3-SAT [9], m/n ∼ 0.55 for 1-2-SAT + and m/n ∼ 0.62 for 1-3-SAT + [10], as shown in Fig. 1(a)(b) and Fig. 2

(a)(b).
We study the case of k = 2 and k = 3 for a comparison: while 2-SAT and 1-2-SAT + are in class P and efficiently solvable, 3-SAT and 1-3-SAT + are NP-complete and it takes an exponential amount of time to solve it, e.g., by the well-known algorithm X [25].Despite the contrast in the decision versions, Max-k-SAT and Max-1-k-SAT + are always NP hard, even for k = 2 [26].In addition to the worst case hardness, empirical studies with classical algorithms on different variants of 3-SAT [6][7][8][9][10] show that when m/n is small (large), almost all instances are satisfied (unsatisfied) and easy to solve; while for m/n approaching the critical point of the SAT-UNSAT transition, the SAT problem instances become the hardest to solve.
To solve SAT problems with QAOA, we transform each Boolean variables v i to the spin states of a qubit, with spin-down state |1 (Pauli-Z operator σ z = −1) for true and spin-up state |0 (σ z = 1) for false, and obtain the spin Hamiltonians for 2-SAT and 3-SAT as where A a ,a stands for the literal sign for th element in ath clause with +1(−1) for positive(negative) literal separately and 0 for absence of it in the clause.Similarly, for 1-3-SAT + and 1-2-SAT + as [21][22][23]27] (see Methods) The gate-based implementation of the QAOA for our problem Hamiltonians can be found in Supplementary.Note V B. With the above encoding, an instance is satisfied only if the ground state energy is zero.As QAOA minimizes the cost function, it can be considered as an approximate algorithm for solving Max-1-k-SAT + .By default, via a threshold decision, the solution of the optimization also implies the solution to the decision version.
Our overall goal in this paper is to understand what is hard and what is easy on QAOA, both in terms of trainability and accuracy.

B. Gradient of QAOA training
As a variational circuit, QAOA's cost-function gradients over variables γ, β indicate the shape of cost-   function landscape-larger amplitudes of gradients indicate sharper changes and therefore the problem is easier to train; while small amplitudes of gradients leads to barren plateaus [16,[28][29][30][31] that make the training difficult.As gradients average to zero on random states [28,29], we evaluate the standard deviation (SD) of gradients to characterize their typical amplitudes.
To represent the typical case of training, we evaluate the gradient on random choices of the circuit parameters via a numerical finite-difference.Without loss of generality, we consider the gradient over the first variable γ 1 and denote it as ∂ 1 C [28,29].To enable an easier visualization in Fig. 1(c)(d) and Fig. 2(c)(d), we plot the inverse of the gradient SD, 1/SD ∂ 1 C( γ, β) , so that large values indicate hardness in convergence.We consider different number of layers p in QAOA to obtain a comprehensive picture of it.
For all the problems under study, the inverse gradient SD has a clear peak at a critical clause-to-variable m/n, as shown in Fig. 1(c)(d) and Fig. 2(c)(d).However, this peak is in general different from the classical SAT-UNSAT transition indicated by the dashed line.For the special case of 1-3-SAT + , Fig. 2(c) shows that the peak of the inverse gradient coincides with the SAT-UNSAT transition.A large inverse gradient SD indicates a small gradient in the typical case, and therefore a more barren plateau that makes the training hard at the phase transition.When p is small, the peak disappears; however, at the same time, QAOA fails to provide the accurate solution, making trainability irrelevant.
We notice that the cases of k = 3 have the inverse gradient peaked at a much smaller clause-to-variable density, as a result of the more complex clauses.Overall, the results reveal a transition of the trainability measured by gradient that is different from the classical SAT-UNSAT transition, showing that the empirical hardness for quantum algorithms can be different from classical algorithms.

Connection to controllability.
To understand the different behaviors of the gradient, we utilize the connection between gradient and controllability measured by the dimension of dynamical Lie algebra (DLA) of QAOA generators, as recently identified in Ref. [16].
As explained in Ref. [13], DLA can be used to test the controllability of the quantum system governed by unitary dynamics.Let us consider an n-qubit system described by a Hilbert space H. Considering an optimal quantum control model described by a unitary is constructed by the repeated and nested commutators of the elements in G.The corresponding dynamical Lie group is therefore obtained by taking the exponential of the DLA e g ≡ {e Generally, for a finite-time evolution governed by Schrödinger's equation, the system is fully controllable when the set of the unitaries obtained during this evolution can cover all unitaries as its elements.This is precisely formulated by the so-called Lie algebra rank condition, which states that the system is fully controllable if and only if dim(g) = 4 n − 1 [32].For quantum systems where the whole Hilbert space is not fully controllable, when the DLA g can be described as the direct sum, i.e. g = j g j , so that the Hilbert space can be written in a form of the direct sum of the subspace H j as H = j H j , dim(g j ) determines the subspace controllability of H j [14][15][16].
With a problem Hamiltonian H C and the mixing Hamiltonian H B as the generators, we can generate the DLA g and provide an estimate of the standard deviation of gradient from the dimension of the DLA dim (g) as where "poly" denotes a polynomial function.Here, for two functions f (x) and g(x), f (x) ∈ Ω(g(x)) means f (x) is bounded below by g(x) asymptotically.
Therefore, we evaluate the DLA dimension numerically to compare with our gradient results.As the numerical evaluation is costly, we are limited to a smaller size of n = 6.Despite the small size, as we see in Fig. 1(g)(h) and Fig. 2(g)(h), the DLA dimension dim (g) essentially has the same behavior versus the clause-to-variable ratio m/n, when compared with the inverse gradient; This manifests a clear connection between the gradient transition and the DLA dimension transition.
QAOA provides a clear physical insight to the concept of trainability of variational quantum algorithms on NISQ device.From Fig. 1(c)-(h) and Fig. 2(c)-(h), the trainability and controllability have a trade off.This can be explained as the following.When the system is more (less) controllable, there are more (less) control protocols available to transform the initial state to the desired final state.Geometrically, these protocols can be described as the the accessible paths characterized by the parameters ( γ, β) from the initial state to the desired state.In this picture, our task is to find the optimal path from all the possible paths.Therefore, from the trainability perspective, it becomes harder (easier) to train ( γ, β) when the system is more (less) controllable.Despite being harder to train, the plurality of the paths also provides more hope to good performance.In addition, one can also connect DLA to controllability via the Quantum Fisher information matrix (QFIM) [30].As the rank of QFIM characterizes the number of independent ways to vary the control parameters to change the generated quantum state, it is intuitive that the dimension of DLA upper bounds the rank of QFIM, which connects the training difficulty and controllability of a quantum model.
An additional insight can be obtained by considering evaluating the DLA dimension at the m n limit for the 1-k-SAT + problems, where the Hamiltonian is symmetric between all qubits (see inset of Fig. 2 (g) and (h)).In this case, we are able to prove an upper bound (see Methods), dim(g) ≤ dim UB ≡ 1 6 n(n 2 + 6n + 11) .We also expect dim(g) to be above n 2 , which is the dimension for a much simpler nearest neighbour Ising model [16].While the upper bound is in general a loose one, it indicates that the gradient in the m n limit is only polynomially small; in contrast, for the hard instances we would expect an exponentially small gradient.This contrast supports the decay of dimension and the increase of gradient when m/n is large.

Barren plateau and complexity.
To further understand the barren plateau phenomena, we study the speed of decay of typical gradient with the number of qubits.To begin with, we pick two values of system size, n = 6 and n = 16 for k-SAT (n = 18 for 1-k-SAT + instead), and evaluate the ratio of the average gradient variance SD(∂ 1 C( γ, β)) for different number of layers p.In Fig. 1 (e)(f) and Fig. 2 (e)(f), we see that right after the peak of the inverse gradients, when the circuit depth p is sufficiently large, the decay ratio saturate to a value independent of the clause-to-variable ratio, indicating the barren plateau.While below the peak of the inverse gradient, the decay ratio of gradient increases gradually with the clause-to-variable ratio.
To confirm an exponential decay of gradient, in Fig. 3, we focus on 1-k-SAT + and plot the SD of gradient versus the layer p for different number of qubits n, while keeping the clause-to-variable ratio m/n to be a constant in each panel.For 1-3-SAT + , when m/n is small in panel (a), the SD of gradient saturates and does not decrease versus p or n, showing no barren plateau; At the critical value in panel (b), we see an exponential decrease of SD versus the number of qubits n at large p, confirming a barren plateau.Above threshold, as shown in panel (c), a barren plateau can still be confirmed, however, with larger gradients than the critical case of panel (b).On the contrary, for 1-2-SAT + as we see in panel (d) and (e), at around the SAT-UNSAT transition we do not see the appearance of a barren plateau.At large m/n in panel (f), the gradient finally starts to show an exponential decay, indicating a barren plateau.
The appearance of barren plateaus are often connected to the complexity of the typical quantum circuit involved [29].The complexity of an ensemble of unitaries can in general be characterized by the closeness to unitary t-design, which reproduces the Haar random expectation values of 2t-point correlators.In this regard, when the quantum circuit forms a 2-design [17,18,33], it has been shown that the variance of the gradient will vanish exponentially with the system size-which leads to a barren plateau of cost function [28,29].
Therefore, we consider the unitary ensemble U p formed by the p-layer QAOA, U QAOA = p =1 e −iβ H B e −iγ H C , with each angle γ ∈ [0, 2π) and each angle β ∈ [0, π) independent and uniform random.To measure the closeness to 2-design, we evaluate the ensemble-averaged infinite-temperature 4-point out-of-time-order correlator (OTOC) where the dimension d = 2 n for an n-qubit system and the average is over the unitary U ∈ E [18].For ensemble E forming a 2-design, we have C OTO (W 1 , W 2 ; E) = −d/(d 2 − 1) saturate to the Haar results [18,33]; while for trivial ensembles, C OTO (W 1 , W 2 ; E) is of order one.Therefore, the decay of OTOC indicates the ensemble being a 2-design.Without loss of generality, we consider the OTOC between single-qubit operators C OTO (σ y 1 , σ y n/2 ; U p ).In Fig. 1 (i)(j) and Fig. 2 (i)(j), we find the OTOC of the QAOA ensemble decays towards the Haar value when clause-to-variable ratio m/n increases to the critical value of the minimum gradient, indicating a transition to 2-design.We also see a difference between the cases of k = 3 versus k = 2-the decay of OTOC for k = 2 is much slower than k = 3.

C. Accuracy of QAOA
In this section, we explore the accuracy of QAOA in solving k-SAT and 1-k-SAT + .To speedup the training, we develop a heuristic pre-optimization initialization strategy (see Supplementary Note V D).To obtain the best accuracy, we perform 10 repetitions on QAOA for each instance to obtain the optimal solution among those results.To benchmark the accuracy of QAOA with the classical algorithm, in the case of Max-k-SAT, we consider the lower bound of state-of-the-art approximation algorithms; In the case of decision versions of k-SAT, we consider success probability of random guess; In the case of 1-k-SAT + , as less results are known about approximation ratios, we reduce the problem to the maximum weighted independent set (MWIS) problem [34,35] and utilize the greedy approximate MWIS algorithms proposed in [36,37] (see Methods).The standard accuracy characterization of approximate algorithms for optimization problem is the approximation ratio [11,36,37].For our case of Max-SAT problems, we define the approximation ratio r ≤ 1 of a solution to be the ratio between the number of clauses satisfied by the solution and the maximum number of clauses that can be satisfied by any solution.As the output state |ψ( γ, β) in QAOA can be in a superposition of multiple solutions, we evaluate the expected approximation ratio via projecting the output state to the computational basis.For Max-k-SAT, a random guess will satisfy on average m(1 − 1/2 k ) number of clauses; For Max 1-k-SAT + , a random assignment will satisfy on average mk/2 k number of clauses.For the instances with most clauses satisfiable, the above corresponds to an approximation ratio of r rand ∼ 1 − 1/2 k and r rand ∼ k/2 k .An exact optimal solution will saturate r = 1 and nontrivial approximate algorithms should have r ∈ [r rand , 1].
In Fig. 4 (a)(b) and Fig. 5 (a)(b), we see that as p increases, QAOA is able to obtain larger approximation ratios.As the clause-to-variable ratio m/n increases, the approximation ratio decays as expected.However, the decay is rather slow and manifest a robustness of QAOA.In the Max-3-SAT case, we see at small p the approximation ratio is already better than the lower bound of r ∼ 0.95 in Ref. [38], similarly, the lower bound of r ≥ 21/22 for the Max-2-SAT [11] case is also overcame  at small depth p.In the case of Max 1-k-SAT + , we consider the approximation ratio of the classical MWIS approximate algorithm for comparison (see Methods).For Max-1-3-SAT + , we identify a clear quantum advantage at around p ∼ 16.For Max-1-2-SAT + , advantages appear even for a shallow depth of p = 8.We want to emphasize that the quantum advantage happens only when the clause-to-variable ratio is large, despite the reachability deficits [12].Indeed, we expect quantum algorithms to be advantageous especially for hard problems, where both classical and quantum algorithms face challenges.Although all Max-SAT problems being considered are NP-hard, we do see some interesting contrast in the performance.In the Max-2-SAT and Max-3-SAT cases, the approximation ratio performance is similar when p is large, consistent with previous results in Ref. [12]; however, for the absolute number of additional violated clauses, Max-2-SAT performs slightly better than Max-3-SAT (see Supplementary Note V D).In the Max-1k-SAT + cases, the accuracy of QAOA is substantially higher for k = 2 than k = 3 with the same number of p layers.We speculate such a contrast in the performance can be caused by the different connectivity and complexity of the Hamiltonian in the problems.
To connect to the empirical hardness transition in classical algorithms, we can also reinterpret each optimization result as a decision of SAT/UNSAT.This can be The probability that state through QAA evolution lies in the ground state done via a threshold decision on the minimized number of UNSAT clauses, e.g., determine an instance as SAT when the expected number of UNSAT clauses is smaller than E th = 0.5 and UNSAT otherwise.To characterize the overall performance, we evaluate the success probability of deciding SAT/UNSAT when solving random instances at a fixed clause-to-variable ratio m/n.The results are shown in Fig. 4 (c)(d) and Fig. 5(c)(d).
The success probability increases with the layer of QAOA p as we expect.For the k = 3 cases in Fig. 4 (c) and Fig. 5(c), there is a valley of low success probability at around the critical point of m/n shown in Fig. 2(a), recovering the same hardness transition identified in empirical studies of classical algorithms [6][7][8][9][10].While for the k = 2 cases in Fig. 4 (d) and Fig. 5(d), despite a similar valley of low success probability at small p, the success probability is almost unity for a circuit depth of p = 16.Similarly, the classical benchmark can be reinterpreted and similar transition versus m/n can be seen.Such a valley at the classical SAT-UNSAT transition indicates a remnant of the classical empirical hardness and is different from the trainability transition identified in Fig. 1 and Fig. 2.
Combining the above, we see that overall QAOA possesses a similar notion of what is hard and easy as classical algorithms, while showing advantage over the classical algorithms being considered in the large problem-density instances.

D. Comparison with quantum adiabatic algorithm
The identified trainability transition in general deviates from the classical computational phase transition, while the performance in solving the SAT problems show a consistent trend with the classical computational phase transition.Such a disparity intrigues us to explore the empirical hardness of SAT instances in other Hamiltonian-based quantum algorithms, such as QAAa popular alternative and also the predecessor of QAOA.
To obtain the solution, QAA prepares the ground state of the problem Hamiltonian in Eqs ( 1) and ( 2) via an adiabatic evolution of the Hamiltonian ⊗n in a superposition of all possible spin configurations.As one tunes the parameter s slowly towards H(1) = H C , the adiabatic theorem guarantees that the state of the system stays in the ground state; therefore, the final state |φ QAA is the ground state of the problem Hamiltonian, which provides the solution to the optimization problem.
From the adiabatic theorem, we can obtain an estimation on the computation time of QAA as T a ∼ 1/(∆E) 2 so that the success probability is close to unity, where ∆E is the minimum gap of the Hamiltonians {H(s), s ∈ [0, 1]} [21,23].In Fig. 6(a)(b) and Fig. 7(a)(b), we evaluate the inverse gap square 1/(∆E) 2 as a measure of the instance hardness with different clause-to-variable ratio using Qutip [39].Note that there exist other more rigorous estimations on the necessary adiabatic evolution time, combining higher-order terms [40,41].However, the inverse gap square as a approximate estimation is sufficient for our purpose.In subplots (a), we identify a computational phase transition for 3-SAT and 1-3-SAT + , where the minimum gap is minimum at about the critical SAT/UNSAT transition, up to some small deviation due to finite size, similar to the decision version of SAT on QAOA in subplots (c) of Fig. 4 and Fig. 5.While for 2-SAT and 1-2-SAT + , we see the minimum gap to be higher than the critical point, qualitatively agreeing with the case of QAOA.
To further confirm the transition, we evaluate the success probability P = D i=1 | ψ i |φ QAA | 2 from the overlap between the final evolved state |φ QAA and the Ddegenerate ground state of the problem Hamiltonian |ψ i .In Fig. 6(c)(d) and Fig. 7(c)(d), at a finite time T a , we see the success probability of QAA decreases before the critical point, while roughly maintaining a constant above the critical point.Such a robustness to the problem density coincides with the slow decay of QAOA's approximation ratio with the problem density.In addition, we also find the success probability of 2-SAT and 1-2-SAT + (subplots (d)) to be much higher than that of 3-SAT and 1-3-SAT + (subplots (c)).

III. DISCUSSION
In this paper, we thoroughly explore the empirical hardness of Hamiltonian-based quantum algorithms in solving SAT problems.In the case of QAOA, we find a trainability phase transition, where the the gradient is minimum at certain critical problem density.Such a phase transition is connected to the controllability and complexity of QAOA circuits.Although the trainability transition in general deviates from the classical SAT-UNSAT transition, in terms of performance, Hamiltonian-based algorithms do show a remnant of the classical SAT-UNSAT transition.Although our results are empirical, we expect analytical results to be challenging, as the classical correspondence of such transition is also empirical due to the complexity of the SAT problems.We also identify quantum advantages of QAOA against several classical greedy approximate algorithms for a relatively small-scale quantum system, potentially realizable in the near-term.Although we have focused on two cases of SAT for convenience, we expect the computational phase transition in QAOA to apply to all combinatorial optimization problems.In particular, as 3-SAT is NP-complete, the clause-to-variable ratio represents a universal characterization of a 'problem density' [12], and the computational phase transition applies to all NPcomplete problems in this regard.

A. Many-body Formulation of the problem Hamiltonian
Here we introduce the many-body formulation of the problem spin Hamiltonian in Eqs. ( 1) and ( 2).For convenience, we first introduce an n × m binary matrix A ij , where A ij = 1(−1) if the variable v i is included in the clause c j as positive(negative) literal and A ij = 0 otherwise; With this matrix in hand, we can express all Hamiltonians in a standard many-body form [35] as for k-SAT problems; and for 1-k-SAT + problems.The notations are introduced as Note that the number of clauses containing v i is equal to |h i |.

B. Dynamical Lie algebra: definition and bounds
Below we give bounds for dim (g) in the m n limit for both 1-2-SAT + and 1-3-SAT + , where all coefficients J ij 's and h i 's approach uniform (see Supplementary Note V A).
For 1-2-SAT + , we have H C,2 + ∝ i<j σ z i σ z j up to a constant.Then, the set of the initial generators for the corresponding DLA For the fully coupled Ising model with transverse fields along x and y axis, the set of the initial generators becomes G x,y ≡ {G 2 , n i=1 σ y i }.
From Ref. [42], the dimension of the corresponding DLA g x,y is where a b ≡ a!/(a − b)!b! is the binomial coefficient.Since the DLAs are generated by the repeated and nested commutators of the generator sets, we must have dim(g Also we know nearest neighbour Ising model has dimension n 2 [16].Therefore, we expect the scaling to be between Ω(n 2 ) and O(n 3 ).
For the 1-3-SAT + , we have + ,H B be the corresponding DLA.Here, because we can write if we start from the initial set of generator , the corresponding DLA of G 3 becomes exactly g x,y .Therefore, we have dim(g H C,3 + ,H B ) ≤ dim(g x,y ), which leads to The lower bound estimation n 2 for the dimension of DLA from nearest neighbour Ising model still holds.

C. Classical approximate algorithms
To solve Max-1-k-SAT + , we transform it to the MWIS problem.Given a 1-k-SAT + instance with n variables and m clauses, one can construct a weighted graph with n vertices {q i } n i=1 , each corresponding to a variable v i and having the weight w(q i ) = |h i | ≥ 0. For every two distinct vertices q i , q j , an edge (q i , q j ) exists if J ij > 0when the corresponding variable v i , v j appear in at least one clause at the same time.
One can verify that the SAT/UNSAT version of 1-k-SAT + problem is reduced to asking whether the weight of maximum independent set is equal to m or not.The reason is simple: an independent set of this graph corresponds to an assignment that does not have more than one true assignment in any clause.To guarantee a solution to the 1-k-SAT + instance, we still need to make sure that all clauses have one true variable.As the total weight of the independent set is equal to how many clauses are satisfied by this assignment, therefore if the total weight is equal to m, all clauses are satisfied.At the same time, the Max-1-k-SAT + can be reduced to solving the MWIS.As a classical benchmark, we can utilize various greedy algorithms for MWIS [36,37] and choose the best performance among them (see Supplementary Note V C).
to the definition J ij = m a=1 A ia A ja , the probability that J ij = J is where we can directly see that the distribution of J is dependent on the number of variables n and number of clauses m.With the distribution of J ij , the mean and variance are Similarly, for 1-2-SAT + , the the probability that A ia A ja = 1 for arbitrary two different i and j is p 2 = 1/ n 2 , and the probability that which is also size-dependent.The mean and variance are As the ratio between standard deviation to the mean of J ij decreases with the clause-to-variable ratio m/n for fixed n, we expect in the limit of large m/n, the coefficients J ij approach uniform for all i, j.The same applies to h i 's for H C,3 + .At the end of the discussion, we want to address the difference between 1-k-SAT + and the wellknown Sherrington-Kirkpatrick (SK) model in spin glass with the Hamiltonian H SK = i<j J ij σ z i σ z j where J ij is independently sampled from the standard normal distribution N (0, 1) [43].

B. Gate-based implementation of QAOA
To implement the Hamiltonian dynamics in QAOA with a quantum circuit, one can decompose the unitary evolution into parallel Pauli-X and Pauli-Z gates as the following for 3-SAT and it is similar for H C,2 by taking K ij = 0 and a factor of 2 in the denominator of exponents.For 1-3-SAT + , the problem Hamiltonian layer is The case of H C,2 + is similar, with all h i 's equaling zero.
The first and second product in Eq. ( 17) and ( 18) correspond to parallel Pauli-Z rotation (RZ) rotation with ZZ interaction (RZZ) gates); and the unique third product in Eq. ( 17) correspond to rotaion with ZZZ interaction (RZZZ gates).Similarly, exp (−iβ k H B ) is implemented by Pauli-X rotation (RX) gates.Numerically, we implement the QAOA with Qulacs [44], a high-performance quantum computing platform for both Python and C++.We employ the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [45][46][47][48], a gradient-based quasi-Newton method implemented in Scipy [49], to find the optimal parameters γ * , β * .The classical optimization stops when either the difference of cost function between steps or gradient norm is smaller than 10 −6 .Our numerical simulations are performed on the Puma HPC from University of Arizona with 50 cores of AMD Zen2 CPU and 250GB of RAM.

C. Details of the classical approximate algorithms
Variants of greedy algorithms are proposed for approximate MWIS problems [36,37].Before applying those algorithms for benchmark, we introduce some notations to avoid any confusion.Given a weighted graph G(Q, E, w), where Q, E, w represent the set of vertices, edges and weights of vertices, we use w(q i ) denotes the weight of vertex q i and w(S) denotes the sum of weight for vertex set S. N (q i ) represents the set of vertices that are adjacent to vertex q i and N + (q i ) = N (q i ) ∪ {q i }.We denote the degree of vertex q in graph G i as d Gi (q).
We briefly summarize the four greedy algorithms that are used for benchmark in this paper, GWMIN, GW-MAX, GWMIN2 [36] and WG [37] algorithms.We also list their corresponding guaranteed lower bounds on maximum weight estimation.To obtain benchmarks, we reduce the 1-k-SAT + instances to MWIS instances with igraph [50] and apply the approximate algorithms introduced above to obtain the approximation ratio results in Fig. 8.As all the four algorithms are variants of greedy algorithms, their performances are similar.The lower bounds of those algorithms in Eq. ( 19) are also plotted in Fig. 8(a)(b) as a reference.The final benchmark presented in the main paper are obtained from the best approximation ratio among the four algorithms for every fixed clause-to-variable ratio m/n separately.

On approximation ratio.
Here we provide a brief summary of known facts of the approximation ratio of problems related to the 1-k-SAT + , as there are not much known results for 1-k-SAT + itself.Note that these results do not carry over to the 1-k-SAT + , as we explain below.
One can reduce a 1-k-SAT + instance to a k-SAT instance.There is a simple polynomial-time algorithm that provides a (1−1/2 k ) approximation ratio for Max-k-SAT (in this paper we always mean exact k variables in each clause).This is also equal to the expected approximation ratio of a random assignment.Ref. [51] shows that the lower bound of approximation ratio for polynomial-time algorithm in solving Max-3-SAT to be r ≥ 0.95.Ref. [11] also shows that it is NP-hard to approximate Max-2-SAT with any approximation ratio above 21/22 0.955.However, it is important to note that the above are worst case results and does not directly apply to 1-k-SAT + due to the reduction.
There are also results on approximate algorithms guaranteeing certain approximate ratios for random instances of Max-k-SAT.However, as instances generated by reducing the random instances of 1-k-SAT + to k-SAT are by no means random, these results also do not apply to random instances of Max-1-k-SAT + .QAOA performance and computing time.In the main text, we have provided the approximation ratio results in Fig. 4 and Fig. 5.Here we also plot the absolute error ∆ (number of additional violated clauses) for k-SAT and 1-k-SAT + in Fig. 9, Fig. 10 We point out that for 1-3-SAT + the ground energy E 0 of H C,3 + could be different from the number of violated clauses, because the energy cost of a violated clause can be different among possible spin configurations; While for 1-2-SAT + there is no such difference.This difference means that despite H C,3 + is a well-accepted Hamitlonian for 1-3-SAT + problem [23,34], there may be other better options.
Initialization strategy of QAOA.
For a p-layer QAOA, the pre-optimization strategy initializes the first p layers  approximation ratio.In particular, as the number of layers p increases, more parameters need to be optimized and the performance of the strategy gets much worse than the pre-optimization strategy, especially for k = 3 with a more complex H C,k .As for the decision version, the performance with the simple ran-dom initialization strategy (shown in Fig. 11(c),(d)) is still worse than the performance of the pre-optimization strategy in Fig. 5(c),(d) of the main paper.In Fig. 11(e)-(h), we plot the steps and computing time for the random strategy; Comparing with Fig. 10(c)-(f), we see the computation cost is at the same order of magnitude with the pre-optimization strategy.

Performance with limitations.
With the existence of barren plateaus at a large depth, finding the best optimal parameters γ * , β * in QAOA consumes a large amount of computing resources for large problem instances.At the same time, the accumulation of errors and noise in near-term devices prohibits a large quantum system to be stable for a long time [2].Therefore, we consider the sub-optimal performance of QAOA under a cutoff T on the optimization steps, shown in Fig. 12.
It turns out that for the n = 10 qubit system p = 16 the performance quickly at only around a hundred steps in most parameter regions.In particular, for the decision version of the problem (Fig. 12 (c)(d)), the easy problems away from the transition only takes a few optimization steps to solve.The only exception is the Max-1-k-SAT + problem at large m/n, where around a thousand steps are necessary.This is due to the optimization versions of the problem being harder at large m/n ratios, recovering the reachability deficits [12].Overall, for a given p, it is efficient to take a limited number of optimization steps in practical implementations to get a balance between accuracy and resource consumption.
a random instance of the SAT problems can be constructed by choosing m clauses C = {c a } m a=1 , each containing k different variables {v aj } k j=1

Figure 1 .
Figure 1.SAT-UNSAT phase transition and trainability of 3-SAT (left column) and 2-SAT.(a)(b) Probability of SAT for different system size n.(c)(d) The mean of 1/SD ∂1C( γ, β) in different p-layer QAOA with n = 16 variables.The inverse is added for easier comparison according to Eq. (3).(e)(f) The ratio of average gradient variance SD(∂1C( γ, β)) for n = 6 variables over n = 16 variables in different p-layer QAOA.Larger ratio indicates barren plateau.(g)(h) The dimension of dynamical Lie algebra dim(g) for generators in QAOA in an n = 6 qubits system.(i)(j) Average of 4-point OTOC for different p-layer QAOA with n = 10 variables.Green horizontal dashed line represents the value given by Haar unitary −2 n /(4 n − 1).Vertical dashed lines indicate the critical SAT-UNSAT transition point.All results are evaluated over 100 instances.

Figure 2 .
Figure 2. SAT-UNSAT phase transition and trainability of 1-3-SAT + (left column) and 1-2-SAT + .(a)(b) Probability of SAT for different system size n.(c)(d) The mean of 1/SD ∂1C( γ, β) in different p-layer QAOA with n = 18 variables.The inverse is added for easier comparison according to Eq. (3).(e)(f) The ratio of average gradient variance SD(∂1C( γ, β)) for n = 6 variables over n = 18 variables in different p-layer QAOA.Larger ratio indicates barren plateau.(g)(h) The dimension of dynamical Lie algebra dim(g) for generators in QAOA in an n = 6 qubits system.The inset log-log plots (g2) and (h2) show dim(g) versus n with a fully symmetric HC .Green and orange curves represent the lower bound estimate dim(g) = n 2 and upper bound of dim UB .(i)(j) Average of 4-point OTOC for different p-layer QAOA with n = 10 variables.Green horizontal dashed line represents the value given by Haar unitary −2 n /(4 n −1).Vertical dashed lines indicate the critical SAT-UNSAT transition point.All results are evaluated over 100 instances.

Figure 3 .
Figure 3. Barren plateau.SD of gradient SD (∂1C (γ, β)) versus the layer of QAOA p for 1-3-SAT + (top) for 1-2-SAT + (bottom) problems with different number of variables n.From left to right we plot at three different ratios m/n = 0.2, 0.8, 2.Due to the finite size n ≤ 18, as shown in Fig.2the transition happens at around m/n = 0.8.

Figure 4 .
Figure 4. Accuracy of QAOA.(a)(b) Approximation ratio r of SAT clauses, (c)(d) Success probability in determining SAT/UNSAT for 3-SAT (left) and 2-SAT (right) versus clause-to-variable ratio m/n with n = 10 variables.Green dashed line represent the lower bound of approximation algorithm r ≥ 0.95 for Max-3-SAT [38] and r ≥ 21/22 for Max-2-SAT [11].The horizontal light green dashed line in (c)(d) represent the success probability of the random guess which are 7/8 and 3/4 for 3-SAT and 2-SAT separately.Vertical black dashed lines in all plots represent critical point of SAT-UNSAT transition.

Figure 5 .
Figure 5. Accuracy of QAOA.(a)(b) Approximation ratio r of SAT clauses, (c)(d) Success probability in determining SAT/UNSAT for 1-3-SAT + (left) (from p = 4 to p = 24) and 1-2-SAT + (right) (from p = 4 to p = 16) versus clause-tovariable ratio m/n with n = 10 variables.Green dots represent the classical approximate results through a reduction to MWIS.The horizontal light green dashed line in (c)(d) represent the success probability of the random guess which are 3/8 and 1/2 for 1-3-SAT + and 1-2-SAT + separately.Vertical black dashed lines in all plots represent critical point of SAT-UNSAT transition.

Figure 6 .
Figure 6.QAA gap size and performance for k-SAT problems.(a)(b)The median of 1/∆E 2 of QAA with n = 10 variables shown by blue circles.The green and purple circles represent the SAT and UNSAT instances separately.(c)(d) The probability that state through QAA evolution lies in the ground state P = D i=1 | ψi|φQAA | 2 .

Figure 7 .
Figure 7. QAA gap size and performance of 1-k-SAT + .(a)(b)The median of 1/∆E 2 of QAA with n = 10 variables shown by blue circles.The green and purple circles represent the SAT and UNSAT instances separately.(c)(d) The probability that the state through QAA evolution lies in ground state P = D i=1 | ψi|φQAA | 2 .
from an ancillary Hamiltonian H B at s = 0 to the problem Hamiltonian H C at s = 1.Here the ancillary Hamiltonian H B = n i=1 |h i |σ x i , where |h i | is the number of times the variable v i appear in the clauses (see Methods).The initial Hamiltonian at H(0) = H B , with a easy to prepare ground state |ψ (0) ∝ (|0 − |1 )

Figure 8 .
Figure 8. Benchmark on approximate ratio r with MWIS approximate algorithms for instances reduced from (a) 1-3-SAT + and (b) 1-2-SAT + .Circles represent approximate ratio r of different algorithms and triangles represent the corresponding lower bound.All results are fixed with n = 10 variables.

Figure 9 .Figure 10 .
Figure 9. (a)(b) Approximation error ∆, and (c)(d) CPU time (in seconds) for a p-layer QAOA versus clause-to-variable ratio m/n to solve 3-SAT (left) and 2-SAT (right) with n = 10 variables.Black dashed line in (a), (b) represents the critical point of SAT-UNSAT transition.

Figure 11 .
Figure 11.(a)(b) Approximation ratio r, (c)(d) success probability, (e)(f) optimization steps and (g)(h) CPU computing time (in seconds) for a p-layer QAOA versus clause-to-variable ratio m/n to solve 1-3-SAT + (left) and 1-2-SAT + (right) with n = 10 variables.All results are obtained with the simple random initialization strategy.
(1 ≤ p < p) by an optimized player QAOA [52], and sample the rest of the parameters {γ k } p k=p +1 , {β k } p k=p +1 randomly uniformly in [0, ].With the initialization, further training gives the optimal parameters.Here we choose = 0.1 to take advantage of the p -layer QAOA results without being trapped in local minima.Now we compare the pre-optimization strategy to a simple random initialization strategy, where all initial parameters {γ k } p k=1 , {β k } p k=1 are randomly sampled uniformly in [0, π].We show the performance on Max-1-k-SAT + in Fig. 11(a)(b).Compared to the preoptimization strategy in Fig. 5(a)(b) of the main paper, the simple random initialization strategy leads to a worse