Optimizing counterdiabaticity by variational quantum circuits

1. Introduction

Quantum computing has been of significant interest for a long time now due to the successes in the fields of machine learning [1,2], quantum simulation [3], optimization [4–6] and others. Especially, solving optimization problems with quantum computers as they are expected to outperform classical computers in the noisy intermediate-scale quantum (NISQ) era, see current review [7]. Among various methods [8–11], adiabatic quantum optimization (AQO) [12] is an effective method to tackle these problems. In AQO, the optimization problem is encoded into a problem Hamiltonian such that the ground state corresponds to the solution of the problem. Then, a simple Hamiltonian, whose ground state is easy to prepare, is taken and adiabatically evolved to the problem Hamiltonian. The adiabatic theorem guarantees that the problem Hamiltonian will be in the ground state with a high success probability. The problem with this is that adiabatic processes are slow, and becomes unfeasible. With the digitized adiabatic evolution [13], AQO can be reformulated into a gate model. However, the long time evolution requires lots of gates, which spoils the desired results due to the errors in the noisy quantum devices.

To circumvent this, shortcuts of adiabaticity (STA) was proposed, to implement fast adiabatic-like control protocols [14,15]. Since then, these have become powerful methods with applications in quantum information processing, and more generally quantum computing during the past decade. These methods include counterdiabatic (CD) driving [16,17], also known as transitionless quantum driving [18,19], fast-forward approach [20,21] and invariant-based inverse engineering [22,23]. Particularly, local CD driving and its variations have been used for many-body spin systems, e.g. to find the ground state [24–29]. Recently, the digitized counterdiabaticity for quantum computing has been proposed [30] in the context of quantum gates and circuits. It is found that the inclusion of CD interactions shallows the circuit depth significantly, when compared with previous digitized adiabatic computing [13]. The advantage of this method has been analysed in quantum annealing [31], quantum approximate optimization algorithm (QAOA) [32–35], and variational quantum eigensolver [36], two typical quantum-classical hybrid algorithms, intending to find the ground state of a target Hamiltonian. Apart from this, the implementation of these methods has been also reported in various scenarios, such as factorization [37] and portfolio optimization [38]. Moreover, a wide range of optimization problems has been further tackled using the drastic advantage, illuminating that CD driving serves as a non-stoquastic catalyst in the acceleration of quantum adiabatic algorithms [39].

With the advantages of the counterdiabaticity, a challenging task is to find the efficient CD terms, see [40–42]. Instead of finding CD terms using spectral knowledge, adiabatic gauge potentials [43] have been proposed to achieve the approximate CD terms recently, particularly for many-body spin systems [25,26,29]. The approximate CD terms can be easily found by using the variational principle [25] and Floquet engineering [26], while diagonalizing the Hamiltonian for a many-body spin system is difficult or a time-consuming task to obtain the exact CD term. Along with this, an approximate form of the CD term is pre-selected, and the coefficients are determined by minimizing the action or the nested commutator (NC). However, more high-fidelity control of larger spin systems requires a higher-order NC of CD terms, beyond the two-body interaction [30], which implies the problematic implementation with more CNOT gates. Actually, the optimal (local) CD driving can be achieved by hybridizing the technique of STA through optimal control theory [40,42] and reinforcement learning [41].

In this work, we propose an alternative way to find efficient counterdiabaticity by using a variational quantum circuit, where quantum circuit and classical optimizer are applied to optimize parameters to minimize a cost function. In detail, we first replace the CD coefficients with free optimizable parameters and then use a classical optimizer to optimize them accordingly. The concrete example of Greenberger–Horne–Zeilinger (GHZ) state preparation on the nearest-neighbour Ising model is further considered, see figure 1, to demonstrate that the enhanced performance of optimal CD terms, as compared with digitized adiabatic evolution and QAOA.

• Download figure
• Open in new tab
• Download PowerPoint

The article is organized as follows. In §2, we explain the CD driving and classical optimizer in detail. In §3, we investigate these methods with the nearest-neighbour Ising model and GHZ state preparation. Section 4 is devoted to the comparison with usual QAOA. Finally, the article ends with a brief conclusion in §5.

2. Counterdiabatic driving and optimization

The use of variational quantum circuits has bloomed a lot since most of them are hybrid quantum-classical algorithms that are considered one of the most crucial applications of the quantum computing paradigm in the NISQ era [44]. It consists of two parts, the quantum part that includes a circuit with variational parameters, and the classical part that includes iterative routines to optimize these parameters. Here, we will apply this method to counterdiabaticity, aiming to minimize a given cost function by tuning the free parameters of the quantum circuit to achieve optimal CD terms.

In AQO, the trivial Hamiltonian ${\hat{H}}_m$ is prepared in its ground state and adiabatically evolved to the problem Hamiltonian ${\hat{H}}_p$. This can be encoded by the following Hamiltonian:

where the time dependence of the system is introduced through the parameter $\lambda (t)$, and for the total time $T$, $\lambda (t)$ should meet the condition: $\lambda (0)=0$ and $\lambda (T)=1$. Typically, the adiabatic evolution takes a long time because of adiabatic criteria. For speeding this process up, we add a CD driving Hamiltonian ${\hat{H}}_{\mathrm{cd}}(t)$, to suppress diabatic transitions, such that where ${\hat{H}}_0(t)$ is the original Hamiltonian (2.1) for the problem and ${\hat{H}}_{\mathrm{cd}}(t)$ is the additional CD Hamiltonian. To elaborate the approximate CD term [25,26], from the variational calculation, we write down the appropriate adiabatic gauge potential ${\hat{A}}_{\lambda }^{\ast }$ in terms of NC method [26], yielding Here, the approximate gauge potential defined is determined by a set of coefficients, like ${\alpha }_1,{\alpha }_2,\dots ,{\alpha }_{\ell }$, and where $\ell$ determines the order of the expansion. The approximate CD coefficient ${\alpha }_{\ell }(t)$ can be finally obtained by minimizing the action $S_{\ell }=Tr[{\hat{G}}_{\ell }^2]$, with the operator ${\hat{G}}_{\ell }$ being defined by ${\hat{G}}_{\ell }={\mathrm{\partial }}_{\lambda }{\hat{H}}_0-i[{\hat{H}}_0,{\hat{A}}_{\lambda }^{(\ell )}]$.

Based on the approximate CD term (NC = 1), we use a variational circuit to find optimal CD parameters by the quantum-classical hybrid algorithm. The circuit can be firstly constructed with gates to digitize the time evolution of Hamiltonian reads

where $n=T/\mathrm{\Delta }t$ is the Trotter step, $T$ is total evolution time and $\mathrm{\Delta }t$ is the interval of every Trotter step. Here, ${\hat{H}}_{\mathrm{cd}}(t)$ can be considered a non-stoquastic catalyst [39]. The unitary $U(0,T)$ generates an output state $|{\psi }_{\mathrm{out}}⟩$ where $|{\psi }_0⟩$ is the $N$ qubit ground state ${|+⟩}^{\otimes N}$.

On the other hand, for the classical optimization routine in quantum-classical hybrid algorithms, a variety of optimization algorithms have been proposed but a gradient-based systematic optimization of parameters is crucial. For that purpose, we use simultaneous perturbation stochastic approximation (SPSA), which is an algorithmic method for optimizing systems with multiple unknown parameters. It is appropriately suited to large-scale population models, adaptive modelling and simulation optimization. The applications of SPSA are widely implemented in many problems, such as neural network training [45], statistical model parameter estimation and fault exposure [46], adaptive control of dynamic systems [47,48], and using QAOA with the SPSA algorithm to solve the Max-cut problem [49].

In general, there exist many classical optimizers, such as the adaptive moment estimation algorithm, constrained optimization by linear approximations, with their pros and cons. In contrast to others, SPSA uses only the objective function measurements of the objective function. So SPSA is considered to be efficient in high-dimensional problems in terms of providing a good solution for a relatively small number of measurements. More specifically, a differentiable cost function $C(\theta )$ in the SPSA algorithm is considered, and $\theta$ is a L-dimensional vector. The optimization problem can be translated into finding a new $\theta$ at which $\mathrm{\partial }C(\theta )/\mathrm{\partial }\theta =0$. SPSA starts with an initial parameter vector ${\hat{\theta }}_0$, yielding

where ${\hat{g}}_k$ is the estimate of the gradient $g(\theta )=\mathrm{\partial }C(\theta )/\mathrm{\partial }\theta$ at the iterate ${\hat{\theta }}_k$ based on prior measurements of the cost function, and $a_k$ is the size of step, a positive number. It is robust to any noise that may occur when measuring the function $C(\theta )$, therefore the function $y(\theta )=C(\theta )+\epsilon$, where $\epsilon$ denotes perturbation on the output. Then we can get the estimated gradient at each iteration step where $c_k$ is a positive number and each parameter is simultaneously perturbed by either $\pm c_k$. ${\mathrm{\Delta }}_k={({\mathrm{\Delta }}_{k1},{\mathrm{\Delta }}_{k2},\dots ,{\mathrm{\Delta }}_{kL})}^T$ is a perturbation vector, which is generated by using a zero-mean distribution. In the above equations, $a_k$ and $c_k$ reduce over the optimization iterations to converge to a final result, defined as where $A$ is the number of iterations multiply 0.01, and $a$, $c$ are the initial values. Besides, $N_k$ is used to denote the number of iterations. In what follows we shall combine CD driving with the classical optimizer in the Ising model for GHZ state preparation, and compare the results with CD terms obtained with NC ($\ell =1$) and the usual QAOA as well.

3. Nearest-neighbour Ising model and Greenberger–Horne–Zeilinge state preparation

As a heuristic study, we consider the spin 1/2 Ising chain with nearest-neighbour interaction, which is of importance for many-body systems and quantum optimization problems [7,30]. The Hamiltonian of the nearest-neighbour Ising model can be written as

where $N$ is the system size, $J$ is coupling strength between qubits and ${\sigma }_z$ is the Pauli spin 1/2 operator, and the periodic boundary condition ${\sigma }^{N+1}={\sigma }^1$ is assumed. In the ferromagnetic case $(J<0),$ neighbouring vertices prefer to align to the same spin, whereas the opposite happens in the antiferromagnetic setting $(J>0)$. If $J$ is a random number and can either be positive or negative, then we have what is called a spin glass. The Ising antiferromagnet without an external field resembles a graph maximum cut (Max-Cut) problem that can be solved via combinatorial optimization, see [50]. Here, we focus on the preparation of the entangled GHZ state, $|\text{GHZ}⟩=({|0⟩}^{\otimes N}+{|1⟩}^{\otimes N})/\sqrt{2}$, as the ground state of the Hamiltonian (3.1) in the ferromagnetic case of $h_z=0$ and $J=-1$. Our paradigmatic model is different from the transverse field Ising model, in which there exist quantum phase transitions and the finite-rate adiabatic passage across a quantum critical point can be accelerated by using CD driving by mean-field theory [51].

To encode the GHZ state preparation into the adiabatic quantum computing (see equation (2.1)), we write down the time-dependent Hamiltonian,

where the function $\lambda (t)$ is, without loss of generality, assumed as With the help of the NC method (2.4), we can calculate the approximate CD term. By keeping the first order ($\ell$ = 1), we have from which we have the CD terms, implying the two-body spin interactions, consisting of ${\sigma }_y{\sigma }_z$ and ${\sigma }_z{\sigma }_y$ terms. From (2.3), the Hamiltonian of the CD term can be further obtained as

with ${\theta }_{\mathrm{cd}}=2\dot{\lambda }(t){\alpha }_1(t)J{h}_0$. By minimizing the action $S_1=Tr[{\hat{G}}_1^2]$, where ${\hat{G}}_1={\mathrm{\partial }}_{\lambda }{\hat{H}}_0-i[{\hat{H}}_0,{\hat{A}}_{\lambda }^{(1)}]$, we finally achieve

It is emphasized again that higher-order NC suggests the higher fidelity for GHZ state preparation for the larger spin systems [30]. However, high-order NC implies the many-body interactions, beyond typical two-body interactions, which makes the physical implementation difficult, for instance, in photonic quantum circuits [8]. Hence, it is quite natural for us to focus on the optimization of CD terms with two-body interaction.

As a consequence, the output state is given by $|{\psi }_{\mathrm{out}}⟩=U(\mathit{\theta })|{\psi }_0⟩$, where the unitary operator $U(\hat{\theta })$ (2.5) is

where $\mathit{\theta }=[{\theta }_x,{\theta }_{zz},{\theta }_{\mathrm{cd}}]$, with ${\theta }_x(\mathrm{\Delta }t)=(1-\lambda (\mathrm{\Delta }t))h_0$, ${\theta }_{zz}(\mathrm{\Delta }t)=\lambda (\mathrm{\Delta }t)J$ and ${\theta }_{\mathrm{cd}}(\mathrm{\Delta }t)$. In principle, one can digitize the CD term (3.5) with the coefficient (3.6) in the gate-based circuit to accelerate the adiabatic preparation of the GHZ state. Instead, we further propose a variational quantum circuit to optimize the digitized CD term for the same proposal, by considering ${\theta }_{\mathrm{cd}}(j\mathrm{\Delta }t)$ as free parameters. To this end, the cost function, $C=1-|⟨{\psi }_{\mathrm{out}}|{\psi }_{\mathrm{tar}}⟩{|}^2$, is minimized for the target state $|{\psi }_{\mathrm{tar}}⟩=|\text{GHZ}⟩$. As we discussed before, the SPSA algorithm is adopted to optimize the variational parameters ${\theta }_{\mathrm{cd}}(j\mathrm{\Delta }t)$, with the boundary conditions ${\theta }_{\mathrm{cd}}(0)=0$ and ${\theta }_{\mathrm{cd}}(T)=0.$ Thus, the previous CD term (3.5) becomes to achieve higher fidelity, where the coefficient ${\tilde{\theta }}_{\mathrm{cd}}$ is inferred from the optimal (site-independent) parameters ${\tilde{\theta }}_{\mathrm{cd}}(j\mathrm{\Delta }t)$, in the digitized quantum circuit, based on the SPSA algorithm.

In figure 1a, we show the form of the quantum circuit, where the parameters are either ${\theta }_{\mathrm{cd}}$ or ${\tilde{\theta }}_{\mathrm{cd}}$ optimized using SPSA. We will implement the results on an ideal quantum simulator based on Qiskit [52], where the number of iterations $N_k=100$, $a=0.15$ and $c=0.1$ are set.

Figure 2 displays the density matrix representation at $t=T$ for the GHZ state preparation in a 5-qubit system with the approximate and optimal CD terms based on NC ($\ell =1$). Obviously, it is exemplified in figure 2 that the results with both types of CD terms adapt to the target GHZ state for $J=-1$. However, figure 3 further shows the optimized CD terms by SPSA have the better performance, where the results with approximate and optimal CD terms based on NC ($\ell =1$) are compared with the normal adiabatic evolution. In detail, without CD driving, the final fidelity of the prepared GHZ states in a 5-qubit system is quite low, due to the adiabatic error when the total time $T=1$ is very short, and does not fulfil the adiabatic criteria. However, the approximate CD terms are not valid at all for all the cases of large spin size $N$ and interaction systems $J$ [30]. With the assistance of the SPSA optimizer, the optimal CD terms improve the fidelity. So our method provides a simple but efficient approach for improving the fidelity of entangled GHZ state preparation by keeping the structure of two-body interaction. Of course, one can also extend this to higher-order NCs, with the many-body interaction involved. Also, other methods including genetic algorithms [53] and reinforcement learning [41] are interesting to incorporate for further exploration with richer structures of CD terms.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

Meanwhile, figure 4 illustrates the coefficients of the digitized CD terms, ${\theta }_{\mathrm{cd}}(j\mathrm{\Delta }t)$ and ${\tilde{\theta }}_{\mathrm{cd}}(j\mathrm{\Delta }t)$, in the quantum circuit for the 5-qubit system, where their intensities are compared. As is apparent, the fidelity is improved by the optimized CD terms, using the SPSA algorithm, at the cost of the intensity of CD driving, proportional to the energy. Finally, in figure 5, we check the fidelity of the entangled GHZ state preparation with various interactions $J$, where the rest parameters $J=-1,-0.6,-0.1$ and $h_z=0$ are considered in the Ising model with five qubits. It is concluded that the results of optimal CD terms outperform those of the approximate one, with the same NC ($\ell =1$).

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

4. Comparison with the quantum approximate optimization algorithm

Typically, QAOA is a quantum algorithm that attempts to solve such combinatorial problems [54–56]. QAOA falls under the category of VQAs where $p$-layers of two unitaries $U_m(\beta )$ and $U_p(\gamma )$ are applied iteratively to an initial state $|{\psi }_i⟩={|+⟩}^{\otimes N}$, in the computational basis where $N$ is the system size. Thus the final state is given by

where $U_m(\beta )={\mathrm{e}}^{-\mathrm{i}\beta \sum _{i=1}^N{h}_0{\sigma }_x^i}$ and $U_p(\gamma )={\mathrm{e}}^{-\mathrm{i}\gamma \sum _{i=1}^{N}J{\sigma }_z^i{\sigma }_z^{i+1}}$. $U_m(\beta )$ is known as the mixer unitary and $U_p(\gamma )$ is known as the problem unitary. ($\gamma$, $\beta$) are free optimizable parameters that are tuned by the classical optimizer to minimize the cost function $C^{\mathrm{\prime }}$, which shows infidelity of the output state

We come up with a perspective of QAOA from quantum control by having an eye on equation (4.1). Indeed, QAOA represents the variational quantum control task that alternatively evolves the problem Hamiltonian $H_p$ and the mixing Hamiltonian $H_m$ for the operation time of ${\gamma }_i$ and ${\beta }_i$, respectively. Thus, it shares the same expression in equation (3.2) for our problem, where $\lambda (t)$ is no longer continuous, but a binary value controller $\lambda (t)=\{0,1\}$ instead, reproducing a bang-bang control toward a target state. In this way, we reckon that it is necessary to compare our protocol, a digitized quantum annealing with counterdiabaticity, with QAOA by bounding the operation time, as a critical energetic quantity that affects the performances of bang-bang control and quantum annealing.

In the following, we compare the performance of QAOA ($p=1$ and $p=2$) with the optimized CD evolution in the respect of fidelity with a bounded time of the algorithms. We define the time $T^{\prime }=\sum _{i=1}^p({\gamma }_i+{\beta }_i)$ to compare the performance through preparing the GHZ state for $h_0=-1$ and $J=-1$, where the target state $|{\psi }_{tar}⟩=|\text{GHZ}⟩=({|0⟩}^{\otimes N}+{|1⟩}^{\otimes N})/\sqrt{2}$. Also, we apply the same SPSA as a classical optimizer while performing QAOA, see figure 1b.

In table 1, we have compared the fidelities in the same time with the optimized CD driving, which shows the fidelities of QAOA are lower than the results from optimized CD terms at $T^{\mathrm{\prime }}=T=1$. We have to emphasize that the shorter time is of significance to prevent the state evolution from decoherence, and gate errors. As expected, the fidelities decrease as the system size increases. Compared with hybrid quantum-classical algorithms like QAOA, the optimal CD works better with bounded time, where the parameters for SPSA are the same, $N_k=100$, $a=0.05$ and $c=0.05$. In fact, QAOA can be also accelerated by adding the CD terms, $\alpha {\hat{H}}_{\mathrm{cd}}$, as a richer ansatz in a variational circuit [32,33], which is beyond the scope of our paper.

Last but not least, we briefly analyse the errors in our protocols, resulting in the state preparation’s infidelity. The adiabatic error, i.e. excitation induced by the non-adiabatic process, is not cancelled perfectly, even if we have introduced optimized CD terms. It can be reduced by prolonging the operation time, or equivalently scaling up the mixing and problem Hamiltonian. The Trotter error is also introduced by the first-order Suzuki–Trotter decomposition [57] in our digital quantum computing paradigm. We realize that the Trotter error increases since either two of the mixing Hamiltonian, problem Hamiltonian, and the CD terms do not commute. However, the Trotter error is still on the scale of $O(\mathrm{\Delta }{t}^2)$ by extending the Baker–Campbell–Hausdorff formula to the case of three components as

where $A=-i{H}_m\mathrm{\Delta }t$, $B=-i{H}_p\mathrm{\Delta }t$ and $C=-i{H}_{\mathrm{cd}}\mathrm{\Delta }t$. Actually, the detailed error analysis can be worked out separately using Trotter formulae [58], with/without CD terms involved.

5. Conclusion

In summary, we have implemented a variational quantum circuit to find optimal coefficients of CD driving, in order to speed up digitized adiabatic quantum computing. By suppressing non-adiabatic transitions but keeping high-fidelity, we have exemplified this method for the GHZ entangled state preparation in the nearest-neighbour Ising model and within a short time. The results have demonstrated that our method is superior to the approximate CD term, based on NC, by minimizing the action. With the first-order NC ($\ell =1$), the coefficients of two-body interaction have been optimized by SPSA, a classical optimizer with variational quantum circuits. In addition, our method is still better than the usual QAOA in terms of fidelity with bounded time, with reminiscences of the comparison between shortcuts to adiabaticity and optimal control theory with respect to robustness and operation time, etc. [59]. Of course, by using a richer ansatz of CD terms (e.g. with N-body interaction), one can also further increase the fidelity in quantum circuits, but the physical implementation becomes tough with the state-of-the-art architecture of NISQ devices. Moreover, other techniques of optimization, such as reinforcement learning [32], greedy and genetic algorithms [53] are worthy of further exploration and comparison. Finally, it would be interesting to apply this method for various complex systems [7], e.g. the Ising antiferromagnet model and transverse-field Ising model, solving the combinatorial optimization problems, and addressing the fundamental and relevant issues on the circuit complexity, gate error mitigation, Trotter and adiabatic errors.

Data accessibility

This article has no additional data.

Authors' contributions

D.S.: formal analysis, investigation, writing—original draft, writing—review and editing; P.C.: formal analysis, investigation, writing—original draft; Z.-H.X.: investigation, writing—original draft, writing—review and editing; X.C.: conceptualization, investigation, methodology, project administration, supervision, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

The authors declare no competing interests.

Funding

This work was supported by NSFC (grant no. 12075145), STCSM (grants no. 2019SHZDZX01-ZX04), EU FET Open Grant EPIQUS (899368), QUANTEK project (KK-2021/00070), the Basque Government through Grant No. IT1470-22, the project grant PID2021-126273NB-I00 funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe” and “ERDF Invest in your Future” and the Ramón y Cajal program (RYC-2017-22482).