Towards adiabatic quantum computing using compressed quantum circuits

We describe tensor network algorithms to optimize quantum circuits for adiabatic quantum computing. To suppress diabatic transitions, we include counterdiabatic driving in the optimization and utilize variational matrix product operators to represent adiabatic gauge potentials. Traditionally, Trotter product formulas are used to turn adiabatic time evolution into quantum circuits and the addition of counterdiabatic driving increases the circuit depth per time step. Instead, we classically optimize a parameterized quantum circuit of fixed depth to simultaneously capture adiabatic time evolution together with counterdiabatic driving over many time steps. The methods are applied to the ground state preparation of quantum Ising chains of sizes N = 7 – 31 with transverse and longitudinal fields. We show that the classically optimized circuits can significantly outperform Trotter product formulas. Furthermore, we discuss how the approach can be used for combinatorial optimization.


I. INTRODUCTION
To make use of current quantum computing technologies, adiabatic time evolution [1][2][3][4] is a promising concept that underpins, for example, the computations of the famous D-Wave device [5] and motivates the popular quantum approximate optimization algorithm [6].In the adiabatic approach, the quantum device prepares a trivial ground state and then realizes time evolution with a specific, time-dependent Hamiltonian, whose initial ground state is the trivial one and whose final ground state encodes the desired quantum computational result.The procedure is successful if the evolution occurs sufficiently slowly and the total evolution time scales ∝ 1/∆ 2 [7] where ∆ is the minimum energy difference between the ground and the first excited state during the evolution.On digital quantum computers, time evolution can be readily realized using Trotter product formulas [8], but long evolution times can lead to deep circuits.One powerful strategy to reduce the adiabatic evolution time is based on the inclusion of additional terms in the Hamiltonian that suppress unwanted transitions during the time evolution [9][10][11][12][13][14][15], which we refer to as counterdiabatic driving [16,17].These additional Hamiltonian terms, however, lead to deeper Trotter circuits per time step on a gate-based quantum device [18][19][20][21].Therefore, adiabatic time evolution over many time steps, with or without counterdiabatic driving, can still be a challenge for current digital quantum computers.
In this article, we extend the tensor network [22][23][24] toolbox for parameterized quantum circuit (PQC) [25][26][27] optimization of Hamiltonian simulation [28][29][30] to tackle adiabatic quantum computing enhanced by counterdiabatic driving.The key ingredients of counterdiabatic driving, that suppress unwanted transitions, are the auxiliary Hamiltonian terms [9][10][11][12][13][14][15] which we refer to as the adiabatic gauge potential (AGP) [16,17,31].We study the suitability of a variational matrix product operator (MPO) [32,33] ansatz to approximate the * conor.mckeever@quantinuum.com AGP.The MPO ansatz naturally fits within our tensor network approach and we also show that it can have advantages over the popular nested commutator (NC) ansatz [34].Following [30], we divide the total evolution time into several chunks of shorter time intervals.For each chunk we optimize a specific PQC, dedicated to that chunk, to capture the counterdiabatic dynamics over the shorter evolution time of the chunk.After the classical optimization, the sequence of PQCs represents the entire counterdiabatic evolution and can be evaluated on a quantum computer.In the context of transverse-field quantum Ising chains with a longitudinal field, we numerically demonstrate that, compared with Trotter circuits, the classically optimized PQCs can improve ground state fidelities by a factor of five and energy accuracies by a factor of three.
Our approach is inspired by several important articles.Firstly, variational quantum algorithms [26,27,35] can be used to simulate quantum dynamics [36][37][38][39] and beautiful proposals exist that make use of counterdiabaticity by including additional gates in the PQC ansatz that come from counterdiabatic driving [21,[40][41][42][43].Because variational quantum algorithms can have a large quantum computational overhead, e.g. when small gradients need to be measured on quantum computers [44][45][46][47][48], we instead perform the entire PQC optimization on classical computers.Furthermore, our PQC ansatz does not necessarily need to contain these additional counterdiabatic gates.Provided that it remains efficient for tensor network methods, the PQC architecture can be freely chosen, e.g. a PQC with a quantum-hardware-efficient structure composed of hardware-native gates can be used.Secondly, classical tensor network algorithms have already successfully simulated the entire evolution corresponding to certain adiabatic protocols [49][50][51][52].In contrast to these simulations, our procedure does not require the classical tensor network optimization to be capable of capturing the entire evolution; instead, only short chunks of the evolution need to be classically simulable.The PQC compression of time evolution for each chunk can be performed such as to exhaust the capabilities of classical computing; then appending all PQCs can create a quantum circuit that is hard to simulate for classical computers but efficient to run on a quantum computer.
The article has the following structure.Section II contains the description of the numerical methods, including the variational optimization of the MPO gauge potential in Sec.II A and the PQCs for the adiabatic evolution in Sec.II B. In Sec.III, we present the results of this study, i.e. the comparison between the nested-commutator and the MPO gauge potentials in Sec.III A as well as the comparison between the Trotter and the classically optimized circuits in Sec.III B. In Sec.III C we compare our methods to other tensor network based techniques.We discuss how to adapt the circuit optimization procedures for other problems, such as combinatorial optimization, in Sec.IV.Technical details are provided in the Appendixes.

II. METHODS
The goal of our approach is to approximate adiabatic evolution, including counterdiabatic terms, as a quantum circuit.The adiabatic dynamics obey a Schrödinger equation where the time-dependent Hamiltonian H(λ(t)) has the form shown in Fig. 1  To proceed, we first describe in Sec.II A how to variationally approximate the AGP at any point along the adiabatic path as an MPO.We then describe in Sec.II B how the full counterdiabatic dynamics can be approximated using a product of classically optimized quantum circuits.

A. Variational matrix product operators for adiabatic gauge potentials
We obtain the gauge potential A(λ) in Fig. 1 (a) using a variational procedure [16,17] based on an MPO ansatz that we numerically optimize to minimize a certain cost function, see Appendix B. The MPO ansatz is characterized by the MPO bond dimension χ that determines the expressive power and computational cost related to the ansatz [32,33].For the optimization, it is useful to write the variational gauge potential operator Ã as a state | Ã⟩ (by making use of the well-known Choi-Jamiołkowski isomorphism).Then the cost function can be conveniently FIG. 1.The proposed method to approximate counterdiabatic dynamics by a product of classically optimized quantum circuits as explained in Sec.II.This example uses M = 3 chunks of L = 2 layer brickwork parameterized quantum circuits (PQCs).
expressed in terms of a system of linear equations cf. Eq. (B16), where 1 denotes the identity operator, H is H(λ) without the gauge potential term, i.e. in our case H = (1 − λ)H i + λH f , and (•) T is the transpose of (•).
Note that the matrix The variational optimization proceeds using tensor network algorithms for linear equations [53,54].These algorithms sweep over the MPO tensors and, for each tensor, solve a specific system of linear equations.Because the matrices of these linear equations can be rank-deficient or ill-conditioned, we add a small regularization parameter η, i.e. the diagonal matrix η1, to them before solving the linear equations.

B. Classical optimization of quantum circuits for counterdiabatic dynamics
For the circuit optimization we use an iterative procedure which generalizes [30] to time-dependent Hamiltonians.The procedure is illustrated in Fig. 1 (b-d).The total evolution over time T is cut into ST slices of shorter evolutions over small time steps τ = 1/S as illustrated in Fig. 1 (b).For the evolution of each slice, we assume λ(t) = λ(t s ) where t s is the midpoint of slice s, so that during every slice evolution the Hamiltonian is time-independent.At each time slice s = 1, 2, . . ., ST , the short-time evolution operator is approximated using a first-order Taylor approximation Here H s , λs and Ãs are equivalent to H, λ and Ã evaluated at time slice s, respectively.The operator ) can be represented efficiently as an MPO using the W I or W II method [55] or, for higher orders of τ , using the approach of [56].The operator Ãs is an MPO representation of the AGP that is calculated using the variational method of Sec.II A. Note that, for Ãs we can alternatively use the NC ansatz [34] or an expansion in terms of a Krylov basis [57].Also, Trotter-based approximations of the small-time evolution operator can be used where appropriate.
In order to approximate the evolution as a circuit, we split the total evolution time T into M chunks and assign the slices to their corresponding chunks, as shown in Fig. 1 (c).For each of the individual chunks, we optimize a PQC U (θ) with variational parameters θ to represent the corresponding evolution as shown in Fig. 1 (d).The PQC corresponding to each chunk m ∈ {1, 2, . . ., M } can be chosen to have any structure amenable to classical optimization.To provide a specific example, we consider the L-layer brickwork circuits U (θ) = U (L) m illustrated in Fig. 1 (e).The PQC optimization then proceeds one slice after another as in [30].More precisely, at slice q of chunk m of the iterative procedure, we search for the parameters F , where ∥•∥ F is the Frobenius norm, U q−1 is the optimized circuit from the previous slice and the Hermitian operator ρ enables us to incorporate knowledge of the initial state.This is equivalent to the minimization of the cost function ) where Re denotes the real part, tr[•] the trace of [•], (•) † adjoint of (•), and U (θ 0 ) = 1.If we have complete knowledge of the initial state |ψ i ⟩, we can set ρ = |ψ i ⟩ ⟨ψ i | for the first chunk (m = 1); we denote the associated cost function as C |ψi⟩ (θ).In the opposite limit, where we have no knowledge of the initial state, the operator ρ is set to the identity ρ = 1 and we denote the associated cost function as C tr (θ).The result of the classical circuit optimization is a set of circuits Parameterized quantum gates Rx(θ) = e −iθX/2 , Rz(θ) = e −iθZ/2 and Uzz(θ) = e −iθZ⊗Z/2 , where θ denotes the variational parameter, for the ansatz in Fig. 1 (e).
which can be evaluated on a quantum computer to prepare the state |ϕ f ⟩ ≈ |ψ f ⟩ as shown in Fig. 1 (e).

III. RESULTS
In the following, we consider the one-dimensional nearest-neighbour quantum Ising Hamiltonian where J k , g k and h k are parameters of the Hamiltonian at site k and X (Z) is the Pauli X (Z) matrix.
In Sec.III B, we optimize PQCs U (θ) which have a brickwork structure as illustrated in Fig. 1 (e).Each of the M chunks of the brickwork circuits consists of an initial layer of single-qubit blocks followed by L layers of two-qubit blocks arranged in a brickwork pattern.The substructure of each block is shown in Fig. 2.

A. Comparison of variational matrix-product-operator gauge potentials to the nested commutator approach
The NC method [34] approximates the AGP in terms of a series of nested commutators of H and ∂ λ H.The order l of the NC approximation corresponds to the number of NCs retained in the series.To facilitate a direct comparison with [34], we consider an N = 14 qubit quantum Ising Hamiltonian (5) at one point along the adiabatic path, such that J k = 1.0, and h k = g k = 0.3λ for all qubit indices k = 1, 2, . . ., N .The resulting Hamiltonian is To assess the performance of the variational MPO gauge potential Ã, we compute the values of the normalized cost, defined as and the normalized error Results for the normalized cost and normalized error, both of which we seek to minimize during the variational optimization, are illustrated in Fig. 3.We observe that the normalized cost in Fig. 3 (a) decreases with increasing variational MPO bond dimension χ and decreasing regularization strength η.For the normalized error in Fig. 3 (b) we observe that, while the error decreases with increasing bond dimension χ, there is an optimal regularization parameter η which minimizes the error for each χ.Additionally, we find that the variationally optimized gauge potentials outperform the NC gauge potential of order l = 6 in terms of the normalized error.
Figure 4 shows the maximum bond dimension of the NC gauge potential when represented as an MPO.We use two different techniques to construct these MPOs as described in the figure caption of Fig. 4. When considered alongside the results of Fig. 3, we observe that the variationally optimized MPO gauge potentials outperform those constructed using the NC approach while also having a much smaller maximum bond dimension.
FIG. 4. The maximum bond dimensions of matrix product operator (MPO) representations of the nested-commutator (NC) gauge potential of order l for the N = 14 qubit Hamiltonian H(λ) defined in Eq. ( 6).For each NC order l the prefactors α k for k = 1, 2, . . ., l of the NCs are calculated numerically (see [34]).To construct Ã as an MPO, we use numerical multiplication and addition of MPOs representing H and ∂ λ H to evaluate and sum over the NCs.As an alternative method, we simplify the NCs algebraically using QuantumAlgebra.jl[58] into sums of Pauli strings and construct Ã using the so-called AutoMPO method [59] included in the ITensors.jllibrary [60].

B. Quantum circuits for counterdiabatic spectral gap traversal
In this section we consider an instance of the quantum Ising Hamiltonian (5) with uniform parameters J k = J, g k = g and h k = h for all qubit indices k = 1, 2, . . ., N ; we refer to this Hamiltonian as H Ising (J, g, h).We consider an adiabatic path that starts with an initial Hamiltonian H GT i = H Ising (0, g * , 0) and ends with a final Hamiltonian H GT f = H Ising (1, g * , 1).The parameter g * is numerically chosen such that the minimum spectral gap appears approximately at λ = 0.5 and the adiabatic evolution realizes a gap traversal (GT) from a paramagnetic to an antiferromagnetic phase.Specifically, for the system sizes N ∈ {7, 15, 23, 31} considered here, we choose g * ∈ {0.48, 0.45, 0.435, 0.43}.The resulting adiabatic Hamiltonian that interpolates between where we choose the path such that the rate of change of the path λ(t) at its beginning and end are zero, so as to smoothly switch on/off the gauge potential at these points.The energy difference between the ground and first excited state of H GT (λ), i.e. the spectral gap, is plotted as a function of λ for system sizes N ∈ {7, 15, 23, 31} in In this example Ã(χ) (λ) is calculated using the method presented in Sec.II A with a bond dimension χ = 4.The entanglement entropy corresponds to that of the matrix-productstate ground state found using standard density-matrix renormalization group (DMRG) algorithms [61].
Fig. 5 (a).We observe that the minimum spectral gap ∆ decreases as the number of qubits increases.
To proceed, we calculate the AGPs for system sizes N ∈ {7, 15, 23, 31}, variational MPO bond dimensions χ ∈ {4, 8} and total evolution times T ∈ {0.1, 0.2, 0.3} using the method of Sec.II A. More precisely, we discretize the total time T into ST slices where, in all cases S = 120, and calculate the AGPs Ã(χ) s for each slice.Since λ(t = 0) = λ(t = T ) = 0, the Hamiltonian with counterdiabatic driving H GT (λ) = H GT (λ) + λ ÃGT (λ) still starts with H GT i and ends at H GT f but differs from H GT (λ) for 0 < λ < 1.The spectral gap of H GT (λ) as a function of λ for N = 23, χ = 4 and T = 0.3 is illustrated in Fig. 5 (b), where we observe that the effect of introducing the gauge potential is to increase the spectral gap as compared to H GT (λ) alone.Additionally we observe in Fig. 5 (c) that the maximum entanglement entropy of the ground state of H GT (λ) along the adiabatic path is increased compared to that of H GT (λ).Here, the measure of entanglement used is the von-Neumann entropy S VN = − j ξ j ln ξ j where ξ j are the Schmidt coefficients resulting from the bipartition of the matrix-product-state ground state at the center of the chain.
Given the choice of interpolating Hamiltonian, the initial eigenstate, i.e. the ground state of H GT (λ = 0) = H GT i , corresponds to the product state |ψ i ⟩ = |−⟩ ⊗N and can easily be incorporated into the cost function C |ψi⟩ (θ) for the first chunk of the circuit optimization procedure.We optimize brickwork circuits of M = 2 chunks and L = 4 layers using the method of Sec.II B. After completing all circuit optimizations we are left with a pair of circuits U m for m ∈ {1, 2} for each parameter set {T, χ, N }.
To compare the performance of the counterdiabatic, classically optimized circuits to more traditional techniques we also construct circuits using a second-order Trotter product formula [62] of the adiabatically interpolating Hamiltonian defined by H GT (λ), i.e., the Trotter circuits do not implement counterdiabatic driving.Importantly, we use R = 8 Trotter layers such that these circuits contain the same number of two-qubit U zz (θ) gates as the classically optimized counterdiabatic circuits of R = M × L = 8 layers.To approximate the adiabatic time dynamics as a Trotter circuit U Trotter , we divide the total adiabatic evolution time T into R = 8 equal time chunks labelled by r = 1, 2, . . ., R. We then construct a second-order Trotter approximation of U Trotter r ≈ e −iT H GT (λ(tr))/R keeping H GT fixed during each chunk, where t r is the time corresponding the the midpoint of the r th chunk.The full Trotterized adiabatic evolution is then given by the time-ordered product In order to quantify the errors due to imperfect adiabatic evolution, we compute the fidelity between the state prepared by each circuit and the corresponding exact Hamiltonian ground states.We also compute the relative error in the energies of these states.These values are defined in terms of the exact ground state of the Hamiltonian H s at the end of slice s, denoted |ψ s ⟩, as well as those of the initial and final Hamiltonians, denoted |ψ i ⟩ and |ψ f ⟩ respectively.We calculate these ground states using standard, numerically exact, density-matrix renormalization group (DMRG) methods.After applying the adiabatic evolution circuits to the initial state, we are left with |ϕ s ⟩ = U s |ψ i ⟩ where U s is the circuit preparing the state at the end of slice s.The target fidelity is defined as and the target energy error is given by the relative error We also consider the instantaneous energy error defined by where we note that E targ s and E inst s are equivalent at the final time slice.
The results are presented in Fig. 6.For the parameters chosen, we observe that classically optimized counterdiabatic circuits outperform the non-counterdiabatic Trotter circuits.In terms of both the target fidelity and energy error, the value obtained is improved using a larger bond dimension for the variationally optimized MPO gauge po-tential.Furthermore, for N = 7 qubits we observe that classically optimized variational circuits outperform the ideal NC dynamics for all orders l ≤ 6. Across all system sizes we observe that classically optimized counterdiabatic circuits outperform the best Trotter result by up to a factor of 5 for the target fidelities and up to a factor of 3 for the ground state energy errors.Note that the energy error at the end of the adiabatic protocol in Fig. 6 appears to be nearly independent of the system size N and, motivated by the similarity with previous results [63], we conjecture that this observation is due to the geometric locality of the terms in the Hamiltonian (5) whose expectation values give the energy error.

C. Quantum circuits for adiabatic ground state preparation
In this section our goal is to compare the capability of our algorithms to adiabatically prepare ground states with the popular tensor network based approach of [64].To achieve high state preparation fidelities, here we utilize that the final state accuracy after adiabatic evolution can be systematically improved simply by increasing the total evolution time T .For large T the rate of change of the adiabatic parameter λ is small such that the counterdiabatic driving is weak and can be ignored.
In the following, we investigate the ability of our methods to find circuits which follow the instantaneous ground state along an adiabatic path defined by the initial Hamiltonian H C i = H Ising (0, −1, 0) and the final Hamiltonian where we choose the path The ground state of the transverse-field Ising Hamiltonian H C f corresponds to the critical point in the thermodynamic limit.
As an ansatz we choose the sequential circuit structure shown graphically in Fig. 7 (a) where the two-qubit FIG. 8. Instantaneous infidelity (16) as a function of time for L = 2 layer classically optimized circuits acting on N = 24 qubits, alongside the ideal adiabatic Trotter evolution.The total adiabatic time is T = 64 and there are S = 8 slices per unit time.The star symbol indicates the target infidelity achieved using a nine-layer sequential circuit by the method of [64].
blocks have the substructure shown in Fig. 7 (b).At each time slice s = 1, 2, . . ., ST , the short-time evolution operator is approximated by a second-order Trotter product formula, which, for the interpolating Hamiltonian (14), can be represented efficiently as an MPO with bond dimension two.
We consider an N = 24 qubit system and optimize the circuit using a single (M = 1) chunk, L = 2 layers of the sequential ansatz and a total adiabatic time T = 64 with S = 8 slices per unit time, i.e. a Trotter time step τ = 0.125.Additionally, we examine the convergence properties of our optimization routine by varying the number of L-BFGS iterations per slice Q ∈ {10, 50, 100}.
Plotted in Fig. 8 is the instantaneous infidelity where |ϕ s ⟩ is the state prepared by the classically optimized circuit (or the ideal adiabatic evolution explained below), and |ψ s ⟩ is the instantaneous ground state calculated using numerically exact DMRG methods.In addition to the instantaneous infidelities associated with the classically optimized circuits, we also plot that of the ideal adiabatic evolution calculated by evolving the initial state along the adiabatic path using the same Trotter circuit in a numerically exact way.We observe that the classically optimized circuit of L = 2 layers maps the initial state to one which closely follows the ideal evolution and the degree to which it does so is dependent on the number of L-BFGS iterations Q performed in each slice.
Using just L = 2 sequential circuit layers, our approach achieves a target state infidelity which is more than one order of magnitude lower than that achieved using a ninelayer sequential circuit by the method of [64], which is indicated by the star symbol at time T = 64 in Fig. 8.

IV. DISCUSSION
In this article, we numerically demonstrate that adiabatic evolution with counterdiabatic driving for specific, one-dimensional quantum many-body systems can be accurately represented by appending shallow onedimensional quantum circuits that are optimized using standard tensor network techniques.We note that the specific quantum many-body ground states considered throughout this article can alternatively be calculated on classical computers with high accuracies using DMRG and the resulting matrix product states can directly be transformed into quantum circuits [64,[66][67][68][69][70].We have provided one example in Sec.III C that shows how our methods can outperform this approach and therefore constitute a valuable addition to the arsenal of tensor network techniques available for this purpose.We anticipate that our adiabatic quantum computing procedure will outperform fully classical approaches in the preparation of highly entangled states for which MPS fail, e.g.large two-dimensional quantum many-body systems [71] or quantum circuits representing adiabatic evolution for combinatorial optimization [72,73].
Our approach can be readily applied to twodimensional quantum many-body systems by making use of the corresponding two-dimensional tensor-network algorithms [22][23][24].
The quantum circuit optimization procedures can also be used to adiabatically solve combinatorial optimization problems, which have historically always been an important application for adiabatic quantum computing [1][2][3][4].For one-dimensional classical systems, it makes sense to use the same approach as for the quantum systems and optimize shallow one-dimensional quantum circuits to realize the adiabatic protocol.Figure 9 shows results for 10 random instances of classical nearest-neighbour Ising chains.We see that, similar to the quantum counterpart, classically optimized circuits with counterdiabatic driving outperform non-counterdiabatic Trotter circuits with the same number of two-qubit gates.We emphasize, however, that one-dimensional classical spin systems can be efficiently solved using classical computers.
For the more interesting, hard classical problems with nonlocal interactions, including quadratic unconstrained binary optimization [79,80], we propose to use a variational circuit architecture consisting of three parts.The first part is a circuit structure that can be dealt with efficiently using classical techniques, for instance a shallow brickwork circuit or a tensor-network-inspired quantum circuit as considered, e.g., in [81][82][83].The second part is composed of commuting many-qubit gates, e.g. one layer of the quantum approximate optimization algorithm (QAOA) [6] where, in contrast to [6], each gate contributes an independent variational parameter as in multi-angle QAOA [84].The third part consists of one For the initial Hamiltonian we choose g k = 0.5 ∀k and J k = h k = 0 ∀k.For the final Hamiltonian, couplings J k and local fields h k are chosen randomly from the sets J k ∈ ±{0.2, 0.4, 0.6, 0.8, 1.0}, h k ∈ ±{0.0, 0.2, 0.4, 0.6, 0.8} and g k = 0 ∀k.For each of the 10 instances, the target energy error (12) of the classically optimized counterdiabatic circuits along the linear adiabatic path, λ(t) = t/T with T = 0.2, from λ = 0 to λ = 1 is indicated by the dashed green lines.Additionally, for each instance we plot the energy error achieved at the end of the adiabatic path, i.e. at λ = 1 only, by a secondorder Trotter circuit following the non-counterdiabatic linear adiabatic path for a range of total times T , up to a maximum total time T = 20.In all cases, the minimum energy achieved is indicated by a pentagon (circle) for Trotter (classically optimized counterdiabatic) and is labelled with its instance number alongside the numerical value of the error.layer of arbitrary single-qubit gates.A quantum circuit consisting of these three parts can be efficiently optimized using classical tensor network algorithms for socalled weighted graph states [85][86][87][88][89].We present an example circuit architecture of this type in Fig. 10.For the classical optimization of such circuit structures to be efficient, however, it needs to be based on expectation values of few-qubit operators [85][86][87][88][89]. Therefore, in the context of these circuit architectures, we propose to replace the cost function ( 4) by the corresponding one derived from McLachlan's variational principle [90], see Eqs. (11) and (12) in [37].Additionally, to represent the AGP, we propose to replace the variational MPO by a variational few-body ansatz whose nonlocality can be systematically increased [16,20,34,91].An interesting possibility is to go beyond previous proposals and study a variational ansatz defined in terms of generic tensors similar to the one in [92] which has proven to be a suc- cessful approach in a different context.
Note added.A few days before our article appeared on the arXiv, a related article was published there [99] (see final journal version [100]) in which the authors also propose to use variational MPOs for AGPs.While their numerical procedures for the MPO approximation of the AGP are similar to ours, they use them in a different context, namely to improve classical tensor network optimization and scan phase diagrams.In contrast, in our article we apply the variational MPO approximation of the AGP to the quantum circuit optimization for adiabatic quantum computing.In this way, the two articles complement each other and we recommend reading both of them.

FIG. 3 .
FIG. 3. (a)Normalized cost C( Ã)(7) and (b) normalized error E( Ã) (8) as a function of the strength of regularization used in the variational matrix product operator (MPO) method, see Sec.II A, for different bond dimensions.The dashed horizontal line indicates the error achieved by the nested commutator (NC) method of order l = 6 as reported in[34].

FIG. 5 .
FIG. 5. (a)The spectral gap of H GT (λ) defined in Eq. (9) as a function of the adiabatic parameter λ for a range of qubit counts N ∈ {7, 15, 23, 31}.(b) The spectral gap and (c) ground-state entanglement entropy for the N = 23 qubit Hamiltonian with and without counterdiabatic driving (CD).In this example Ã(χ) (λ) is calculated using the method presented in Sec.II A with a bond dimension χ = 4.The entanglement entropy corresponds to that of the matrix-productstate ground state found using standard density-matrix renormalization group (DMRG) algorithms[61].

FIG. 6 .
FIG.6.Target fidelity(11) (left column) and energy errors(12) and (13) (right column) achieved using either a secondorder Trotterization of the standard adiabatic path or classically optimized counterdiabatic circuits for the quantum Ising gap traversal problem defined in Sec.III B. For each system size N , all circuits considered contain the same number of two-qubit gates R(N − 1) where R = M × L = 8.The solid black line indicates the target fidelity and energy error achieved by the Trotter circuit at the end of a standard adiabatic sweep, i.e. at λ(T ) = 1, for total evolution times T up to a maximum T = 10.The best values achieved by the Trotter circuits are indicated by a star symbol.The dashed blue and red lines indicate the target fidelities or instantaneous energy errors achieved during the adiabatic evolution using the classically optimized counterdiabatic circuits where the counterdiabatic gauge potentials are calculated as MPOs of bond dimension χ.The maximum target fidelities and minimum energy errors are achieved at the end of the adiabatic path and are indicated by blue circles (red pentagons) for χ = 4 (χ = 8).Additionally, the solid green lines indicate the target fidelities and energy errors achieved by ideal adiabatic evolution aided by an order l nested commutator (NC) gauge potential.

FIG. 7 .
FIG. 7. (a) A sequential circuit ansatz for N = 7 qubits and L = 2 layers.(b) Substructure of the two-qubit blocks of the sequential circuit where three CNOT gates and fifteen onequbit rotations parameterize an arbitrary two-qubit unitary operator [65].

FIG. 9 .
FIG.9.Energy error achieved by circuits of R = M × L = 8 layers for random instances, id ∈ {1, 2, . . ., 10}, of the nearest-neighbour Ising Hamiltonian(5) with N = 8 qubits.For the initial Hamiltonian we choose g k = 0.5 ∀k and J k = h k = 0 ∀k.For the final Hamiltonian, couplings J k and local fields h k are chosen randomly from the sets J k ∈ ±{0.2, 0.4, 0.6, 0.8, 1.0}, h k ∈ ±{0.0, 0.2, 0.4, 0.6, 0.8} and g k = 0 ∀k.For each of the 10 instances, the target energy error(12) of the classically optimized counterdiabatic circuits along the linear adiabatic path, λ(t) = t/T with T = 0.2, from λ = 0 to λ = 1 is indicated by the dashed green lines.Additionally, for each instance we plot the energy error achieved at the end of the adiabatic path, i.e. at λ = 1 only, by a secondorder Trotter circuit following the non-counterdiabatic linear adiabatic path for a range of total times T , up to a maximum total time T = 20.In all cases, the minimum energy achieved is indicated by a pentagon (circle) for Trotter (classically optimized counterdiabatic) and is labelled with its instance number alongside the numerical value of the error.

FIG. 10 .
FIG. 10.Part I is a brickwork circuit of L = 1 layer.The component inside the dashed box needs to be included only in the first chunk optimization of Fig. 1.Part II consists of commuting two-qubit gates acting on all possible pairs of qubits and part III contains generic single-qubit gates.The gate parameterization is defined in Fig. 2.