Digitized-Counterdiabatic Quantum Optimization

We propose digitized-counterdiabatic quantum optimization (DCQO) to achieve polynomial enhancement over adiabatic quantum optimization for the general Ising spin-glass model, which includes the whole class of combinatorial optimization problems. This is accomplished via the digitization of adiabatic quantum algorithms that are catalysed by the addition of non-stoquastic counterdiabatic terms. The latter are suitably chosen, not only for escaping classical simulability, but also for speeding up the performance. Finding the ground state of a general Ising spin-glass Hamiltonian is used to illustrate that the inclusion of k-local non-stoquastic counterdiabatic terms can always outperform the traditional adiabatic quantum optimization with stoquastic Hamiltonians. In particular, we show that a polynomial enhancement in the ground-state success probability can be achieved for a finite-time evolution, even with the simplest 2-local counterdiabatic terms. Furthermore, the considered digitization process, within the gate-based quantum computing paradigm, provides the flexibility to introduce arbitrary non-stoquastic interactions. Along these lines, using our proposed paradigm on current NISQ computers, quantum speed-up may be reached to find approximate solutions for NP-complete and NP-hard optimization problems. We expect DCQO to become a fast-lane paradigm towards quantum advantage in the NISQ era.

Introduction.-Many important optimization problems in science and industry can be formulated as solving combinatorial optimization problems [1]. In general, these problems are known to be NP-hard, so that no classical or quantum computers are expected to solve them efficiently. However, there is a hope that quantum computers might give some polynomial speed-up, which helps reduce the resources and hence the cost of solving many practical problems of interest. Especially, adiabatic quantum optimization (AQO) algorithms are developed to tackle such problems [2][3][4]. Here, the solution to the optimization problem is encoded in the ground state of an Ising spin-glass Hamiltonian [5]. The adiabatic theorem states that the system will remain in the instantaneous ground state if the evolution from the ground state of an initial Hamiltonian to the final Hamiltonian is sufficiently slow enough. The corresponding time-dependent Hamiltonian is given by where σ z and σ x are the Pauli operators. And, λ(t) ∈ [0, 1] is a scheduling function represents the interpolation from the initial Hamiltonian to the final Hamiltonian. In Eq. (1), the first term corresponds to the Ising spin-glass Hamiltonian with all-to-all interactions, finding its ground state in the worst case scenario is NP-hard [6]. The second term with transverse field represents the initial Hamiltonian, corresponding to quantum fluctuation.
The Hamiltonian H ad (λ) has off-diagonal matrix elements that are real and non-positive in the computational basis, in other words, stoquastic and quantum Monte Carlo simulations can tackle such problems without facing any sign problem. It is believed that adiabatic quantum optimization or quantum annealing with stoquastic Hamiltonian might not give significant enhancement over classical algorithms, but some counterexamples have recently been found [7][8][9]. The adiabatic quantum computation with non-stoquastic Hamiltonians is known to be universal. However, the role of non-stoquastic catalysts to speed-up adiabatic quantum optimization problems is an unresolved problem. There are some problems where the non-stoquastic catalysts are advantageous [10][11][12][13][14][15][16][17], and others are known to worsen the performance compared to their stoquastic counterparts [18]. The main reason for this ambiguity is that the non-stoquastic terms are chosen randomly in all the previous works.
The counterdiabatic (CD) technique, borrowing from shortcuts to adiabaticity (STA) [19,20], was introduced to speed up adiabatic evolutions by adding Hamiltonian terms to suppress the non-adiabatic transitions [21][22][23][24]. Recently, a number of developments have shown the advantage of CD techniques in digitized-adiabatic quantum computing [25][26][27][28], quantum annealing [29][30][31][32][33], and quantum approximate optimization algorithms (QAOA) [34][35][36]. In this manuscript, we propose digitized-counterdiabatic quantum optimization (DCQO) as a novel paradigm to solve the general class of combinatorial optimization problems with quantum speed-up. We show that the suitably designed CD terms appearing during the nonadiabatic evolution act as non-stoquastic catalysts. This provides us with a guaranteed enhancement over traditional adiabatic quantum optimization with stoquastic Hamiltonians, while solving the long-range Ising spin-glass problem. Moreover, we consider local approximate CD terms that can be obtained without knowing any prior detail of the Hamiltonian spectra [37][38][39]. We remark that we derive a general analytical expression for the CD coefficients, so that one does not have to calculate them explicitly for each case. Finally, we follow the gate-model approach for the digitized-adiabatic quantum evolution [25,40], allowing the digital implementation of ar-arXiv:2201.00790v1 [quant-ph] 3 Jan 2022 bitrary non-stoquastic CD terms in current noisy intermediatescale quantum (NISQ) computers. Consequently, the proposed DCQO paradigm will be useful to solve the general class of combinatorial optimization problems with quantum speed-up in the NISQ era.
Counterdiabatic driving as a non-stoquastic catalyst.-The concept of STA was originally proposed last decade [19], and was found to have wide applications in many fields, ranging from quantum physics, classical physics to stochastic physics [20]. Among all techniques of STA, CD driving, also called transitionless driving, provided the possibility of tailoring the Hamiltonian of quantum many-body systems [37,38] to speed up the adiabatic process. Here, the source adiabatic Hamiltonian is added by a CD term that takes the form Here, H i and H f are the initial and the problem Hamiltonians connected by λ(t), satisfying λ(0) = 0 and λ(T ) = 1, and H cd is the CD term with scheduling f (λ), which vanishes at the beginning and end of the protocol. In principle, introducing the exact CD term help to follow the instantaneous ground state of the original Hamiltonian H ad (λ) at all times during the evolution. This motivates us to define the CD term as H cd (λ) =λA λ . We can notice that, when the evolution is very fast, H cd increases dramatically, and in the adiabatic limit, i.e., |λ| → 0, the CD term vanishes. Here, | is known as adiabatic gauge potential corresponding to the non-adiabatic transitions between the eigen states |m(λ) [41], which satisfies the condition [i∂ λ H ad − [A λ , H ad ] , H ad ] = 0. Obtaining the exact A λ for a many-body system is challenging. Its implementation is not optimal because, in many cases, A λ can contain exponentially many terms with non-local many-body interactions. An alternative approach is to consider the approximate form of the adiabatic gauge potential [37,39] that can be obtained from a nested commutator (NC) [38] A (l) Here, α k (t) is the CD coefficient, obtained by minimizing the operator distance between the exact gauge potential and the approximate gauge potential. This is similar to minimizing the action S = Tr[G 2 λ ], where the Hilbert-Schmidt operator For the real-valued Hamiltonian in Eq. (1), the adiabatic gauge potential is always imaginary, so the terms appearing in the NC expansion always contain an odd number of Pauli-y terms. So, by restricting to lower-order terms obtained from Eq. (3) and postselecting only one-spin and two-spin interaction terms, one can construct the general 2-local CD term as where the CD coefficients α i , β i j , and γ i j is obtained by variational minimization [37,42]. We can notice that, by con-struction, the approximate CD terms have imaginary numbers as the off-diagonal matrix elements, which makes it nonstoquastic. So, the NC ansatz in Eq. (3) serves as a prescription for choosing the non-stoquastic catalyst. However, not all the terms in the NC ansatz give significant enhancement. Therefore, preselecting the correct operators can help reduce the cost of implementation. Besides, by introducing more control parameters, the use of machine learning techniques and quantum variational algorithms for obtaining optimal values might further enhance the performance [34,[42][43][44][45].
We start with a simple local adiabatic gauge potential, which takes the formÃ λ = Even though the off-diagonal matrix elements forÃ λ contains complex values, with a simple change of basis, one can make it stoquastic. Also, from the first order expansion of the NC ansatz we obtain the 2-local CD term as Here, the CD coefficient α In the last term of Eq. (6), we have the following additional constraints: i = k or j = l, and equivalently i = l or j = k. In Eq. (5), we only have a single variational parameter α 1 (t).
For better performance, one can also consider the CD term in Eq. (4), which requires optimizing N 2 parameters considering β i j = β ji , and γ i j = γ ji . For convenience, we use the shorthand notation, Y =λ N i β i (t)σ y i , and Y|ZY = H (1) cd (λ). Digitized counterdiabatic driving.-For adiabatic quantum optimization problems, the CD terms are non-stoquastic, making it challenging to realize such Hamiltonian on existing quantum annealers [14,46]. Also, for obtaining higher success probability in many-body systems, it is essential to consider higher-order k-local CD terms. This is a challenging task for analog quantum computers and quantum annealers. Along with that, the lack of flexibility in realizing arbitrary interactions is one of the main drawbacks of analog quantum computers and quantum annealers. To overcome all these problems, we used digitized-adiabatic quantum computing techniques [25,40], which provides the flexibility to introduce arbitrary multi-qubit and non-stoquastic interactions. The model is also consistent with error correction [47], and error mitigation techniques are being developed for NISQ computers [48]. Here, we fix the total evolution time T = 1 and the number of trotter steps to 20, i.e. δt = 0.05. The interaction strengths and the local fields are chosen randomly from a Gaussian distribution for 1000 random instances. The blue curve corresponds to the evolution without the CD term, the red curve is with the local single qubit CD term, and the green curve is for the 2-local CD term in Eq. (5). In the latter case, a polynomial enhancement in the success probability can be observed.
For the time-dependent Hamiltonian in Eq. (2), the evolved state is given by |ψ and T is the time ordering operator. The total Hamiltonian can be decomposed into sum of local terms, i.e., H(λ) = j γ j (t)H j . We discretize the total time T into many small intervals of size δt. Using the first order Trotter-Suzuki formula, we obtained the approximate time evolution operator given by where M corresponds to the number of trotter steps. For better approximation one can also consider the recently proposed commutator product formulas [49]. For the evolution of a time-dependent Hamiltonian, the trotter step size δt should be less than the fluctuation time scale of the Hamiltonian [50], i.e., δt ∂H/∂t −1 . However, recently, it has been shown that the digitized-adiabatic evolution is robust again discretization error [51]. This loosens the restrictions on δt. The cost associated with the digitized-counterdiabatic evolution is given by Cost = T max λ H(λ(t)) , which corresponds to the total gate count, while the number of trotter steps decide the circuit depth. For the fully-connected Ising spin-glass problem, total N(N − 1)/2 entangling operations are needed for each trotter step, and the inclusion of the 2-local CD terms increases it by a constant factor. Recently, it has been shown that implementing such fully connected Ising spin-glass problems on current quantum annealers, and also on parity quantum computers, results in huge time overhead because of the embedding schemes [52]. In this regard, it is argued that gate-model quantum computing on a 2-D grid has an advantage compared to other architectures. In this sense, trapped-ion systems with all-to-all connectivity would be an ideal choice but they are not strictly necessary.
In our simulation, we fix the total time T = 1, and the step size δt = 0.05. The scheduling function is chosen as λ(t) = sin 2 π 2 sin 2 πt 2T , so that CD terms vanish at the beginning and the end of the evolution. Each matrix exponential term in Eq. (7) is implemented using standard two-qubit CNOT gates and single-qubit rotations. For all the cases, the number of shots is chosen between 10000-100000 depending on the system size. By measuring the qubits in the computational basis, we obtain the success probability given by P s = ψ g |ψ f (λ = 1) 2 . Here, |ψ g is the actual ground state, and |ψ f (λ = 1) is the time evolved state at t = T . For the Hamiltonian in Eq. (1), the coupling terms J i j and the local fields h i are chosen randomly from a continuous Gaussian distribution with unit variance and zero mean. For estimating the fraction of instances where the inclusion of the CD term gives an enhancement, we define a metric called enhancement ratio, given by R enh = L cd /L. Here, L cd is the number of instances with enhanced performance by including the CD term, and L is the total number of instances, where we set L = 1000. In order to quantifying the improvement in the success probability, we define a quantity called success probability enhancement, given by P enh = P cd s /P ad s . Here, P cd s and P ad s are the success probability with and without the CD terms, respectively.
In Fig. 1, the average ground-state success probabilities for the naive stoquastic Hamiltonian with only transverse field in Eq. (1), including the local CD term Y and 2-local CD term obtained from NC ansatz Y|ZY, are depicted for system size up to 18 qubits. We see that the success probability decreases rapidly with increasing system size for the naive nonadiabatic approach, the inclusion of the 2-local CD term H (1) CD gives a polynomial enhancement, and the local single-spin CD gives a constant enhancement. In Fig. 2, the histogram shows the ground-state success probability distribution for evolution with and without including the CD term for 1000 random instances. The top panel of Fig. 3 shows the average success probability enhancement (P Avg enh ) by including the CD terms Y and Y|ZY. As the system size grows, P Avg enh increases polynomially for the CD term obtained from the first order NC ansatz. With the single spin CD term Y, we obtained an enhancement by a factor of ≈ 3, irrespective of the system size. The bottom panel of Fig. 3 depicts the fraction of instances where the inclusion of the CD term gives enhancement. For the local CD Y, we obtained R enh ≈ 75.6%, while the CD term Y|ZY gives the enhancement ratio R enh ≈ 100%, indicating that the 2-local CD term gives a guaranteed enhancement for all the random instances. In general, the inclusion of CD terms does not help reduce the minimum gap, and its enhancement is explained by the fact that the additional terms help suppress the matrix elements responsible for the excitations between the eigenstates. However, to analyse the effect of non-stoquastic CD terms on the minimum gap ∆ min , we study the instanta-  The average success probability enhancement (P Avg enh ) and probability enhancement ratio (R enh ) as a function of system size for the Ising spin-glass problem is depicted. For the 2-local CD term from the first order NC ansatz, P Avg enh increases with the system size. Whereas for the local CD term (Y), a constant enhancement by a factor of 3 is observed. In the bottom panel, we see that, the 2-local CD term always give enhancement for all the 1000 random instances, whereas the local CD has an average enhancement ratio 0.756. neous energy spectrum as a function of time. Surprisingly, we noticed that the inclusion of approximate CD terms obtained from the NC ansatz increases the minimum energy gap between the ground state and the first excited state during the evolution. This increased gap helps to reduce the excitations, resulting in increased success probability. In Fig. 4, the energy gap between the ground state and the first excited state (|E 1 − E 0 |) as a function time is plotted for four randomly chosen instances with system size N = 10.
Discussion and conclusion.-We have answered a longdebated problem in adiabatic quantum optimization, i.e., the speed-up role of non-stoquastic catalysts. We showed that cleverly chosen non-stoquastic counterdiabatic Hamiltonians achieve enhanced performance compared to traditional stoquastic adiabatic methods. We considered the general Ising- spin glass Hamiltonian with all-to-all connectivity to show that a polynomial enhancement in the ground-state success probability can be obtained, even with 2-local non-stoquastic CD terms stemming from the NC ansatz. As an outlook, considering higher-order k-local CD terms may further enhance the already observed quantum speed up of DCQO paradigm.
In this work, we also provided a general analytical expression for scheduling these CD terms, whose calculation does not require any prior knowledge of the Hamiltonian spectra or the structure of the eigenstates. In conclusion, we proved that the proposed digitized counterdiabatic quantum optimization (DCQO) paradigm involving suitable non-stoquastic CD terms is superior to the traditional adiabatic quantum optimization (AQO) using stoquastic Hamiltonians. In this sense, DCQO might help to achieve quantum advantage for obtain-ing approximate solutions to combinatorial optimization problems on noisy intermediate-scale quantum (NISQ) computers from hundreds to few thousand qubits.