Scaling adaptive quantum simulation algorithms via operator pool tiling

Adaptive variational quantum simulation algorithms use information from the quantum computer to dynamically create optimal trial wavefunctions for a given problem Hamiltonian. A key ingredient in these algorithms is a predefined operator pool from which trial wavefunctions are constructed. Finding suitable pools is critical for the efficiency of the algorithm as the problem size increases. Here, we present a technique called operator pool tiling that facilitates the construction of problem-tailored pools for arbitrarily large problem instances. By first performing an ADAPT-VQE calculation on a smaller instance of the problem using a large, but computationally inefficient operator pool, we extract the most relevant operators and use them to design more efficient pools for larger instances. Given that many problems, such as those arising in condensed matter physics, have a naturally repeating (lattice) structure, we expect the pool tiling method to be widely applicable. We demonstrate the technique here on strongly correlated quantum spin models in one and two dimensions, finding that the resulting state preparation circuits are significantly shorter compared to existing methods.

Variational quantum eigensolvers (VQEs) comprise an important and much-studied class of algorithms for nearterm quantum computers [1][2][3][4].As examples of hybrid quantum-classical algorithms, they divide the computational task of optimizing an objective function f (θ) (for instance, the ground state energy of a physical system) between a quantum and a classical processor, in which the former efficiently queries the objective function, while the latter determines updated guesses for the variational parameters θ.In a VQE, the objective function takes the form of an expectation value f (θ) = H θ of an observable H in the state |ψ(θ) .The latter is prepared by a parameterized quantum circuit, and is subsequently measured to determine the relevant expectation value.
While the choice of H is typically fixed by the problem at hand, the success of the algorithm depends strongly on the form of the circuit that prepares |ψ(θ) (also known as the "ansatz") and the optimization method used.Many different ansatze have been proposed, including hardware-efficient [5] and physics-or chemistryinspired [1,[6][7][8] ones.In general, a good ansatz ought to be expressive enough to closely approximate the solution for suitable values of θ, while maintaining a low circuit depth, for compatibility with the limited coherence times of noisy quantum devices.
An important development in this regard was the introduction of the ADAPT-VQE method [9], which put forward the idea of adaptive VQEs.This approach iteratively constructs an ansatz where θ (n) = (θ 1 , . . ., θ n ) are the variational parameters at step n in the algorithm.In adaptive VQEs, the operators A n are selected from a fixed, predefined "operator pool" according to a selection criterion.In the ADAPT-VQE appraoch, the operator chosen at each step is the one that maximizes the gradient of the objective function, i.e., the A k for which ∂ H /∂θ n+1 | θn+1=0 is maximized, where H is the objective function operator.Intuitively, this criterion seeks to identify the operators A k which are most effective for decreasing the objective function at a given step.The ADAPT-VQE algorithm was first applied to finding the ground state energies of small molecules, where it obtained high accuracy with shallower circuit depths compared to other commonly used ansatze.The approach was subsequently generalized to a number of other applications, including optimization problems [10,11], real and imaginary time evolution [12,13], the preparation of excited states [14], and strongly-correlated lattice models [15].The adaptive philosophy towards ansatz construction has also found other expressions, as in evolutionary algorithms and related methods [16,17].
As one might expect, the success of ADAPT is crucially dependent on the choice of operator pool.The first application of the method used pools consisting of fermionic operators that correspond to single and double excitations of electrons.It was then discovered that qubit-based pools, consisting of individual Pauli strings, can also reach the ground states of molecular systems but with considerably reduced counts of two-qubit gates [18].Subsequent work [19] demonstrated that even more savings can be achieved by symmetrizing the pool of Ref. [18].
In general, the optimal choice of pool for a given problem is subject to a tradeoff: on the one hand, pools with more operators are more likely to lead to convergence (and more quickly), due to the greater flexibility of operator selection at each step in the algorithm.On the other hand, since |∂ H /∂θ| is calculated for each operator in the pool at every step, this can lead to considerable measurement overhead when estimating the expectation val-arXiv:2206.14215v1[quant-ph] 28 Jun 2022 ues on a quantum computer.In seeking to reduce this, it was proven that "minimally complete" pools can be constructed that contain only 2L − 2 operators (where L is the number of qubits) [18].These pools have the property that all states in the Hilbert space can be obtained from an ansatz using only the operators of the given pool.However, due to the local nature of the ansatz updates, it is not guaranteed that the sequence of operators required can be found by the algorithm.One method for alleviating some of these issues is to explicitly encode information about symmetries of the problem into the choice of operator pool [20].
Despite these advances, the space of possible choices for operator pools remains staggeringly large, growing exponentially in the number of qubits.At present, identifying ultimately successful pools is somewhat of an art, often based on a trial-and-error approach in which various candidates are investigated and compared, slowing the pace of discovery.
In this work, we present a general systematic procedure for constructing pools that are both small in size and also facilitate rapid convergence to the ground state.This is achieved by tiling successful pools from small problem instances to construct pools that are effective for larger problems.We then illustrate the method on the example of the spin-1/2 XXZ model, a prototypical correlated spin model.We find that our approach leads to shorter state preparation ansatze compared to previous methods.Finally, we discuss other variations of the underlying idea and their connection to many-body expansions used in quantum chemistry.
ADAPT instances that use larger pools are expected to converge more often and with fewer iterations, owing to the greater number of states that can be accessed at each step.Among qubit-based pools, the largest possible one is the "full Pauli pool", consisting of all Pauli strings on L qubits, excluding the identity.As this pool includes all other qubit pools as proper subsets, it performs at least as well as each of them, in the sense of making the locally optimal choice.However, the exponential scaling of the size, 4 L − 1, clearly prevents its application beyond small problem instances.
Many problems, such as those arising in condensed matter physics, have a simple repeated (lattice) structure that relates instances defined for different numbers of qubits, up to finite-size effects.Exploiting the lattice structure of these systems on near-term quantum devices has attracted increasing attention [21][22][23].Example problems include spin models (as used in the field of quantum magnetism), as well as fermionic ones like the Hubbard model, well-known from the study of hightemperature superconductivity.The lattice structure of these models suggests the following approach to designing effective pools: (1) Start with a small instance of the particular model or problem of interest on L 1 qubits and run ADAPT-VQE with a large, highly expressive pool that can quickly converge to the desired solution.(2) Take the operators chosen by ADAPT during the execution of the algorithm as the basis for a new pool defined on a larger instance of the same problem on L 2 > L 1 qubits.This requires taking the tensor product of each operator with the single-qubit identity L 2 − L 1 times, such that operators in the new pool have the correct dimension.Since the products can be applied either on the left or the right, a single operator of the original pool yields L 2 − L 1 + 1 operators for the new one, corresponding to "tiling" the operator across the system in the onedimensional case (with straightforward generalizations to higher dimensions).Importantly, while the original pool size can be exponential in the number of qubits L 1 (as with the full Pauli pool), the number of operators in the new pool scales linearly with the system size L 2 .
We now demonstrate the viability of this strategy through numerical simulations on a concrete example.We consider both one-dimensional and two-dimensional spin-1/2 XXZ models.The first case consists of a chain of L sites with open boundary conditions, whose Hamiltonian is given by where X i , Y i , and Z i are Pauli operators acting on site i.We set J xy = 1 throughout this work.For the first part of the tiling method, we solve the L = 4 case of the XXZ chain using the full Pauli pool.ADAPT converges in five steps for this simple problem, choosing the operators: YXII, IIYX, ZYXI, XZYI, and YXZI in sequential order (here the site indices have been suppressed and the identity is written explicitly).Convergence of the algorithm is defined to occur when the maximum gradient for a given step drops below a threshold ε.We note that in obtaining this result, we used a modified operator selection compared to the original version of ADAPT: If multiple operators have the same maximal gradient at a given step, we select the one that has the lowest Pauli weight (i.e. the number of non-identity Pauli operators in the string).This is useful both for reducing the number of two-qubit gates (which is better suited for noisy devices) and for identifying particular motifs that can be tiled when going to a larger system.Indeed, given the simple structure of the set of operators chosen by ADAPT, we first generalize the starting L = 4 pool to include other related Pauli strings, before tiling these to produce pools for larger systems.Thus we start from a "Z-decorated XY pool" consisting of the operators: {YX, XY, ZYX, ZXY, YZX, XZY, YXZ, XYZ}.We then tile these operators to create pools for systems with L > 4 and run ADAPT.Fig. 1(a) shows the relative error in the ground state energy compared to the exact solution for this model, for system sizes 4 ≤ L ≤ 16 and several values of the interaction strength J z .In each case, the error is less than 0.04%, indicating good convergence across a range of J z that includes both the paramagnetic and ordered regimes.In Fig. 1(b) we present the number of ADAPT iterations required to find the ground state as a function of L. We note that this value is also equal to the number of optimization parameters in the final ansatz.The behavior over the simulated range appears to be close to linear, such that a simple extrapolation would suggest that 10 3 -10 4 parameters are needed for a system with L ≈ 100 sites.To compare with another popular VQE ansatz, we consider Ref. [24], which studied the XXZ chain for J z = 1 and closed boundary conditions using the Hamiltonian variational ansatz (HVA).For a system with L = 12 sites, 136 variational parameters were required to converge to the ground state.In contrast, the present simulations with ADAPT use only 72 parameters in the final ansatz for the L = 12 open chain, illustrating the ability of ADAPT with the tiling pool to obtain highly compact representations of the ground state.
We now turn to the details of how ADAPT-VQE approaches convergence in the present example.The energy and the absolute value of the maximum gradient at each iteration of the algorithm are shown in Fig. 2(a),(b) respectively.Three distinct regimes are apparent in both quantities.The energy decreases sharply during the first L/2 iterations, while the maximum gradient remains at its initial value for L/2 − 1 iterations.This is followed by a second region in which the energy decreases roughly linearly with iteration number, while the gradient shows a plateau with oscillations between several values.The final regime shows a slower decay both in the energy and the maximum gradient.One can obtain some physical insight into this behavior by inspecting the operators chosen by ADAPT and the corresponding optimal parameter values.The first L/2 operators in the ansatz are of the form e −iθ k X 2k−1 Y 2k with k = 1, . . ., L/2 and θ k ≈ −π/4.Acting on the initial Néel state, this sequence produces a superposition of basis states in which pairs of up and down spins have been swapped.The terms in this superposition all have the same magnitude, but with a nontrivial relative sign structure between them.The other two regimes in the ADAPT convergence process involve the Z-decorated three-local operators, and possess different optimal parameter values which do not form an easily discernible pattern.Nevertheless, it is clear that the operations performed during the middle regime are more effective for reducing the energy overall, while the final regime involves finer tuning of the wave function.
For the two-dimensional XXZ model we study the case of open boundary conditions and nearest-neighbor interactions on rectangular lattices of N = L x L y sites: where i, j denotes nearest neighbors.We first consider the two-leg ladder with L y = 2 and variable L x .In generating the ADAPT operator pool, it is natural to consider tiles consisting of 2×2 blocks, given the two-dimensional nature of the problem.Inspired by the success of the Zdecorated XY pool in the one-dimensional case, we consider two types of 2×2 blocks: (1) those having X and Y on two of the four sites, the other two sites being I; (2) those having X, Y, and Z on three of the four sites, the remaining site being I.In each case, we use all permutations of the locations of the nontrivial operators on the block.These blocks are then tiled across the ladder (in overlapping fashion) to create the operator pool.The procedure is illustrated in Fig. 3 In Fig. 4 we show the convergence properties of the two-leg XXZ ladder (L y = 2) as a function of L x .In each case, the algorithm achieves a ground state energy close to the exact value [Fig.4(a)], with a number of ADAPT steps that is roughly linear in the system length.This behavior is consistent with that of the XXZ chain, as expected from the quasi-1D nature of the model.
A more challenging problem to consider is the 2D XXZ model with L x = L y .In Fig. 5 we present the en- ergy error [Fig.5(a)] and number of steps to convergence [Fig.5(b)] for this case.While the ground state is accurately approximated here as well, the computational difficulty associated with the L 2 x growth in the total number of sites limits our calculations to only three values (L x = 2, 3, 4).Although a quantitative extrapolation to large systems would not be warranted, it is clear that the required number of ADAPT steps grows more quickly with system size than in the 1D case.Notably, this nonlinear growth remains when the number of steps is plotted as a function of L 2 x , suggesting that the computational difficulty is at least partially due to the intrinsic connectivity of the 2D lattice, and not just the total number of sites.In any case, it appears that the 2D XXZ model presents a significant challenge for VQE methods in the quantum advantage regime, due to the large number of optimization steps needed.Nevertheless, the rapid progress in variational quantum algorithms lends hope that such difficulties may be surmounted.
Having demonstrated the power of the "tiling" approach to designing operator pools for ADAPT, we turn to possible generalizations of the method.While the above examples focused on 1D and 2D systems, it is clear how to generalize to higher dimensions (though likely to be computationally challenging).Generalizations to other popular models in condensed matter physics, such as the Hubbard model, are also straightforward.While such constructions are quite natural for lattice systems with short-range interactions, it is less clear that a simple tiling procedure will work for more complex geome- tries and/or long-range interactions (as found in quantum chemistry applications, for instance).In this case, one can consider the more general approach of starting from the full problem, selecting out various subsets of the degrees of freedom, and solving these reduced problems using ADAPT with a large operator pool.The resulting operators chosen from each subproblem can then be collected to generate a pool for executing ADAPT on the original problem.This method bears a close resemblance to many-body expansion techniques used in quantum chemistry.We leave further development of these ideas for future work.
To conclude, we have proposed and demonstrated a systematic operator tiling method for designing small but accurate operator pools for use with the ADAPT-VQE algorithm.The ability of these pools to produce highly compact ansatzes is a significant advantage for using ADAPT on noisy quantum devices.We have illustrated the approach on the open XXZ spin chain and found that it performs favorably compared to the popular Hamiltonian variational ansatz, requiring a significantly reduced number of variational parameters to reach convergence.Confirming the usefulness of the operator tiling method for other models in strongly correlated physics and chemistry is an interesting direction for future work.
Numerical calculations of the exact ground state energies were performed with the QuSpin Python library [25].This work was supported by the Department of Energy.E.B. and N.J.M. acknowledge Award No. DE-SC0019199, and S.E.E.acknowledges the DOE Office of Science, National Quantum Information Science Research Centers, Co-design Center for Quantum Advantage (C2QA), contract number DE-SC0012704.

FIG. 1 .
FIG. 1. (a)Relative error in the ground state energy versus system size L for the XXZ chain, using the Z-decorated XY pool.E is the energy found by ADAPT after convergence, while E0 is the exact ground state energy.(b) Number of ADAPT steps required to converge the algorithm versus system size L, using the Z-decorated XY pool.The initial state is the Néel state | ↑↓↑↓ • • • , and we set ε = 0.01.Variational parameters are optimized using the LBFGS algorithm, with θn+1 = 0 and {θj}j=1,...,n set to the optimal values from the previous ADAPT iteration.

FIG. 2 .FIG. 3 .FIG. 4 .
FIG. 2. (a) Energy and (b) absolute value of the maximum gradient as a function of ADAPT iteration number for the XXZ chain, using the Z-decorated XY pool.Simulation parameters are the same as for Fig. 1 with Jz = 1.

FIG. 5 .
FIG. 5. (a) Relative error in the ground state energy and (b) number of ADAPT steps required for convergence for the 2D XXZ model with a square geometry, using the pool described in Fig. 3. Simulation parameters are Ly = 2, Jz = 1.