Variational Hamiltonian simulation for translational invariant systems via classical pre-processing

Quantum simulation of real time evolution represents a natural application of quantum computation. In particular, non-equilibrium dynamics is known to lead to highly entangled states [1, 2] and is therefore a strong candidate for quantum advantage, i.e. an application where quantum computers can compute quantities of interest that are not accessible to classical architectures [3]. Given a time-independent Hamiltonian H of a closed, physical system of N qubits, one can define the time evolution operator as the unitary

$$\begin{align} {U_{\text{ex}}} = \exp{\left( -i H \tau \right) }, \end{align} \tag{ 1 }$$

where here and throughout the paper we set $\hbar = 1$. This operator propagates a quantum state $\vert \psi(\tau_0)\rangle$ along a time interval $\tau$ to $\vert \psi(\tau_0 + \tau)\rangle$. The dimension of the underlying Hilbert space grows exponentially in N, which makes simulation of large quantum systems intractable for classical computation, in general. A quantum computer, however, can directly operate on the N qubits with quantum gates. When H contains several non-commuting terms, we are unlikely to have a quantum gate of the form of (1) at hand. Thus, a decomposition into experimentally accessible quantum gates is needed. One such decomposition can be achieved by a Trotter–Suzuki, or product formula, approximation [4–7].

Besides product formula methods, there have been other approaches to implement algorithms for Hamiltonian simulation, such as by quantum signal processing [8], Taylor series truncation [9, 10] and linear combination of unitaries [11]. Recently, the gate complexity in Trotter–Suzuki algorithms has been shown to scale comparably or even better than these so called post-Trotter methods [12].

It has been shown [13] that for a permissible Trotter error $\epsilon$ of a qth order product formula, the gate count G is bounded by

$$\begin{align} G \leqslant N 5^{2q} \frac{(NT)^{1+1/2q}}{\epsilon^{2q}}. \end{align} \tag{ 2 }$$

The Trotter error of qth order with Trotter number m scales as $\epsilon = \mathcal{O}(\frac{T^{2q+1}}{m^{2q}})$. For a nearest-neighbor Hamiltonian, for instance, this bound can be further reduced to scale almost linearly in system size and simulation time in the worst case [12]

$$\begin{align} G = \mathcal{O}\left( (N T)^{1+ \mathcal{O}\left( \frac{1}{q} \right) } \right). \end{align} \tag{ 3 }$$

For the simulation of Hamiltonians with local interactions [14] the gate complexity of higher order Trotter–Suzuki formulas meets a lower bound (for $q \to \infty$ in (2)). This lower bound is realized with product formulas [12] as well as decompositions using Lieb–Robinson bounds [14].

Since quantum simulation on today's noisy intermediate scale quantum (NISQ) devices [15] suffers from gate imperfections, variational methods have been put forward to reduce the total gate count and circuit depth [16–19]. While a Trotter–Suzuki decomposition already yields a good approximation to a Hamiltonian-generated unitary time evolution, it has been shown that alternative methods can improve the total fidelity in specific cases by trading off the improvement of approximation accuracy against the gate count. These alternative methods include randomization [20–25], tensor network techniques [26], variational methods [27–33] and also a non-variational optimization of product formulas if there are only few non-commuting terms in the Hamiltonian [34]. While the optimization potential on the asymptotic behavior of the gate complexity for the simulation of local interaction Hamiltonians is moderate, constant factor improvements (in N) can already result in quantum advantage under NISQ conditions. This is due to the fact that the fidelity for a desired output state decreases exponentially in the gate count for finite fidelity gates.

Apart from the dependence on N and T as in (2), there is a constant factor $5^{2q}$ which grows with the order of the Trotter–Suzuki formula q [13]. This factor makes higher order product formulas unfavorable for NISQ devices. On the other hand, finding an optimal gate sequence with variational methods can improve on that factor by keeping q low and thus represent a good NISQ solution where simulation times and system sizes are moderate such that asymptotic behavior is not dominating the gate number. As product formulas represent a universal, i.e. model-independent, method for Hamiltonian simulation, it is not surprising that one can optimize it for a specific Hamiltonian.

Variational methods make use of this optimization potential in different ways. They might exploit hardware-efficiency for a specific ansatz state [28–30], propose an efficient diagonalization of the Hamiltonian within the available gate-set [31, 32] or focus on low-energy initial states [33].

In this work, we tackle two difficulties that arise in these variational methods. First, employing a quantum optimization routine requires multiple iterations with high sampling cost for each. This is a manifest problem in NISQ algorithms [35–38]. Our method allows to move the optimization to a classical pre-processor for translationally invariant Hamiltonians. Second, most of recent variational methods depend on a specific initial state for which the dynamics are calculated [28–31]. In contrast, a time-universal simulation algorithm needs to be applicable for any initial state. In particular, the same circuit can be applied multiple times on a desired initial state and all the intermediate states. The method we present is state agnostic and therefore not limited to applications on states with bound entanglement as various classical methods [39].

We present an initial state agnostic variational algorithm for the simulation of translational invariant Hamiltonians where the heavy sampling effort of the optimization is avoided via classical pre-processing.

For a translationally invariant Hamiltonian, our generalization scheme allows to obtain the quantum circuit of a large system from the optimal solution of a small, classically solvable system. The upscaling to larger systems can be performed qubit-by-qubit if the variational parameters are chosen to be translationally invariant.

For open boundary conditions, which break translational invariance, we propose a variant of the method which uses two classical optimization routines, one which optimizes the open boundary system and one which yields optimal parameters for a periodic boundary system. The upscaled system then combines the first set of parameters for the qubits near the boundary and the second set for qubits within the bulk. Two-dimensional upscaled open boundary systems potentially require classical optimization of more than one small open boundary system, but otherwise follow the same procedure. This scenario is favored by modern NISQ devices which would require a number of SWAP operations to implement periodic boundaries on two-dimensional architectures. Also, open boundary systems are expected to pose larger challenges for classical simulation methods since no translation symmetry can be exploited in classical algorithms. The process of gluing together optimized systems is discussed in detail in section 4 and illustrated in figure 1.

Zoom In Zoom Out Reset image size

The presented methods are numerically exemplified on a one- and two-dimensional transverse field Ising model (TFIM). On small systems for which we can compare the quantum circuit with exact time evolution, our variational approach improves upon the Trotter–Suzuki decomposition in terms of accuracy by up to 3 orders of magnitude in the cost function and at the same time reduces the gate count by a factor of 2. We are then able to suggest quantum simulation circuits for larger systems using the optimal parameters from classical pre-processing, which outperform the Trotter–Suzuki ansatz in a similar way.

Under NISQ conditions, i.e. with imperfect gate fidelities, such improvements turn out to be decisive as to whether an algorithm can be executed with a Trotter–Suzuki ansatz. These scenarios appear when the simulation time τ for fixed circuit depth is large, see figures 2 and 4.

Zoom In Zoom Out Reset image size

In numerical experiments with the one-dimensional TFIM, we found for example that with 50 qubits, a single gate error of $0.1 \%$ and τ = 2.5 our algorithm obtains a NISQ preparation infidelity of ca. $29 \%$, being only a $7 \%$ increase from the minimal infidelity due to gate errors. Compared to approx. $57 \%$ infidelity of the best Trotter–Suzuki ansatz. For this numerical experiment, the optimal parameters for the time evolution on a 6 qubit chain have been used.

We consider a p-local quantum many-body Hamiltonian H which consists of $N \cdot A$ terms

$$\begin{align} H = \sum_{n=1}^N \sum_{a = 1}^A c_a H_{n, a}, \end{align} \tag{ 4 }$$

where $a \in \{1, \ldots, A\}$ counts different interaction types and $n\in \{1, \ldots, N\}$ denotes the first qubit on which $H_{n, a}$ is acting. For a $k-$local term ($k \leqslant p$), for instance $H_{n, a}$ acts on qubits $n, \ldots, n+k-1$ modulo periodic boundary conditions.

The $\{H_{n, a}\}_{(n, a)}$ of the Hamiltonian are chosen such that each $H_{n, a}$ is the generator of a quantum gate that can straight-forwardly be implemented on the considered hardware (depending on the p-locality of H and on the quantum architecture, further decompositions into experimentally accessible gates might be involved).

The unitary time evolution generated by H can be approximated by a first-order Trotter decomposition [4]

$$\begin{align} {U_{\text{ex}}} &= \exp{\left( -i \tau \sum_{n, a} c_a H_{n, a} \right)} = \left[ \prod_{n, a}^\leftarrow \exp{\left( \frac{- i \tau}{m} c_a H_{n, a} \right)} \right]^m + \mathcal{O}\left( \frac{\tau^2}{m} \right), \end{align} \tag{ 5 }$$

where m denotes the Trotter number and the arrow denotes an ordering which starts with the lowest pair (n, a) on the right. The ordering inferred on the tuple (n, a) is not unique and only needs to be fixed when constructing the corresponding circuit. The leading order error of this approximation is given by commutators between the Hamiltonian terms [12]. A qth order Trotter–Suzuki decomposition while using more gates improves the error to leading order $\mathcal{O}(\frac{\tau^{q+1}}{m^{q}})$ [4, 5]. We will stick to first order product formulas as a benchmark in the main text of this paper, since higher order approaches use exponentially more gates in q and hence yield worse total fidelities compared to increasing the Trotter number m. We show in appendix A that also smaller Trotter numbers for higher order Trotter sequences do not save enough gates to yield a better simulation.

Our method is independent of the particular variational quantum circuit ansatz ${U_{\text{var}}}$, which is a gate sequence of parametrized unitaries

$$\begin{align} {U_{\text{var}}}(\vec{\theta}) = U_{1} (\theta_{1}) \cdot U_{2} (\theta_{2}) \cdots . \end{align} \tag{ 6 }$$

We choose the first order Trotter–Suzuki approximation from (5) as a starting point for the circuit ansatz, as this ensures an accuracy that is at least as good as Trotter–Suzuki approximation. In order to exhaust the potential of the ansatz from (5), we will substitute the coefficients $\frac{\tau}{m} c_a$ for each of the Hamiltonians $H_{n, a}$ at each of the m sequences by a variational parameter $\theta_{r, a}$ with $r \in \{1, \ldots, m\}$ labeling the layers and $a \in \{1, \ldots, A \}$ the different gates. The resulting ansatz reads

$$\begin{align} {U_{\text{var}}}(\vec{\theta}) = \prod_{r = 1}^m \left[ \prod_{n, a}^\leftarrow \exp{\left( -\frac{i}{2} \theta_{r, a} H_{n, a} \right)} \right]. \end{align} \tag{ 7 }$$

Similar to the Trotter–Suzuki approximation, the variational ansatz admits an approximation order m which controls the number of parameters. There are thus $A \cdot m$ different parameters to tune. Similar ansätze as in equation (7) are known as variational Hamiltonian ansatz [37] for adiabatic ground state preparation, where they are used with a different cost function than in our case.

The optimization objective is to determine a set of optimal parameters $\theta_{r,a}$ for which the exact time evolution ${U_{\text{ex}}}$ is approximated best by ${U_{\text{var}}}$. The optimal parameters are found by minimizing a cost function measuring the distance of unitaries. We choose a measure of distance induced by the Frobenius norm

$$\begin{align} C(U_\text{var}, \vec{\theta}) &= \frac{1}{2^N} || {U_{\text{var}}}(\vec{\theta}) - {U_{\text{ex}}}\|_F^{\!2} = \frac{1}{2^N} \sum_{a, b=1}^{2^N} \lvert {U_{\text{var}}} {}{^a_{~b}} - {U_{\text{ex}}} {}{^a_{~b}} \rvert^2 \quad \in [0, 1], \end{align} \tag{ 8 }$$

where ${U}{^a_{~b}}$ denotes the elements of a matrix representation of the operator U. The additional factor $\frac{1}{2^N}$ eliminates the scaling of $C(\vec \theta)$ in system size N. The exact operator of time evolution ${U_{\text{ex}}}$ serves as a benchmark which is agnostic of an initial state. In a similar fashion, we can consider $C(U_\text{trot}, \theta)$ as a measure of distance between the Trotter–Suzuki ansatz and exact time evolution.

The cost function (8) is a measure of how well time evolution is approximated on average in the full Hilbert space. As we are considering finite-dimensional Hilbert spaces, any error of observable dynamics can be upper bounded by a distance measure

$$\begin{align} | \langle \psi|{U_{\text{ex}}^{\dagger}} \, \mathcal{O} \, {U_{\text{ex}}}| \psi\rangle - \langle \psi| {U_{\text{var}}^{\dagger}}(\theta) \, \mathcal{O} \, {U_{\text{var}}}(\theta)| \psi\rangle | \leqslant 2 ||\mathcal{O}||_\infty \, || {U_{\text{var}}}(\theta) - {U_{\text{ex}}}||_\infty = 2\kappa ||\mathcal{O}||_\infty \, || {U_{\text{var}}}(\theta) - {U_{\text{ex}}}||_F, \end{align} \tag{ 9 }$$

where $||\cdot||_\infty$ denotes the operator norm. The constant κ, in general, depends on the distribution of the error throughout the Hilbert space. For a concentrated error, we have κ = 1 while for the pathological case that only one state is simulated badly while every other state has perfect fidelity, the upper bound gets loose by a factor of $\kappa = 2^n$. For the Trotter algorithm it is known that most of the time, the errors will be concentrated [25] making the cost function (8) a good estimate for the performance of the algorithm.

As an example, the optimization of one- and two-dimensional TFIMs on rectangular lattices is discussed in section 5, but for now let us assume, we have found a set of parameters such that $C(\vec{\theta})$ as defined in (8) is minimal. Note that (8) involves the full time evolution operator which is not accessible on a quantum computer, but has to be calculated on a classical computer. Consequently, there is no quantum advantage, if ${U_{\text{var}}}$ is deployed to a quantum computer. However, a very meaningful, NISQ-friendly application arises, when the found optimal parameters of the trained smaller system are applied to a larger unseen system, which is explained in the following.

Hamiltonians of many-body quantum systems often enjoy translational symmetries, such that larger systems look similar to smaller systems, locally. In particular, $H_{n, a}$ of (4) is the same operator for all n except that it acts on different qubits. Our method exploits this symmetry and generalizes solutions of variational optimizations to larger quantum systems.

To begin with, the upscaling algorithm is explained for fully translational invariant systems. These necessarily have periodic boundaries. Since some modern quantum hardware, however, naturally comes with open boundary conditions for the implementable gates which disrupt translational invariance, we also extend the method by systems, where the translational symmetry is only broken by open boundaries.

4.1. Translationally invariant systems

Translationally symmetric systems can be scaled up by successive application of the gates from the small system to the new qubits. Figure 1(a), for instance, shows a one-dimensional spin chain with periodic boundary conditions. For nearest-neighbor interactions, the interaction between the first two qubits is the same as the interaction between any other neighboring pair of qubits. From this translational invariance, ansätze for larger quantum systems are generated by copying the quantum gates of the small circuit. More precisely, we apply the quantum gates, which acted on qubits N and 1 before, after extension to qubits N and N + 1 as well as 1 and N + 1, respectively. In principle, this qubit-by-qubit upscaling scheme can be applied repetitively to extend an initially few qubit time evolution circuit to a many qubit one. The final extended variational circuit can then be easily inferred from (7) to be

$$\begin{align} {U_{\text{var}}}(\vec{\theta}) = \prod_{r = 1}^m \left[ \prod_{\substack{n=1, \ldots, N+K \\ a=1, \ldots, A}}^\leftarrow \exp{\left( -\frac{i}{2} \theta_{r, a} H_{n, a} \right)} \right]. \end{align} \tag{ 10 }$$

where K denotes the total number of added qubits. Note that (10) encodes more than just tensorating copies of the systems, but restructures gates such that entanglement between the former N qubits and the new K qubits is created (cf. figure 1(a)).

Very importantly, the method itself is insensitive of the precise model as long as it is translational invariant. As before, we will stick to the Trotter–Suzuki inspired variational ansatz and only need to define the Hamiltonian of the extended system $H^{(N+1)}$ and then use the gate sequence from the copied addends in the Hamiltonian as before. To get this, we subtract the old boundary term and insert the new terms for the new qubit

$$\begin{align} H^{(N+1)} = H^{(N)} &- \sum_a \sum_{k=1}^{p-1} H_{N-(k-1) \rightarrow p-k, a} + \sum_a \sum_{k=1}^{p} H_{N-(k-2) \rightarrow p-k, a} \end{align} \tag{ 11 }$$

where we used the notation $H_{n_1 \rightarrow n_2, a}$ to emphasize that this term acts on the qubits $n_1, n_1 + 1, \ldots, n_2$. This procedure is also sketched for p = 2 in figure 1(a). Ultimately, this upscaling can be repeated to end up with the Hamiltonian for N + K qubits.

This process might seem redundant at first. However, if we recall the variational method to improve on trotterized time evolution, we have just derived a generalization scheme which allows us to use the parameters $\vec{\theta}$ found for the N qubit system to suggest a decomposition of quantum gates for the simulation of the N + K qubit system. This is done by, again, breaking the periodic boundary conditions and gluing the gates corresponding to the K new qubits but now with the variational parameters instead of the Trotter parameters $\frac{\tau}{m} c_a$.

It becomes apparent that the gluing method is identical to Trotterization on a coarser level. This can be understood best, when looking at the reverse process using the Trotter–Suzuki parameters $\frac{\tau}{m} c_a$. Given the Hamiltonian $H^{(N+K)}$ of the large system, $\exp(-i\tau H^{(N+K)})$ is split into a product of N + K identical gates with the only difference of acting on different qubits as described in (11). Two subsequent Trotterizations, first into $H^{(N)}$ and $H^{(K)}$ and afterwards into 2-qubit gates yield the same result as a Trotterization directly into 2-qubit gates. As a result, if the variational parameters are reset to the Trotter–Suzuki parameters, there is no difference in a Trotter–Suzuki decomposition of $\exp(-i \tau H^{(N+K)})$ and gluing K qubits to the system $H^{(N)}$.

Since the leading error of Trotter approximations is given by commutators between Hamiltonian terms, this error scales linearly in the qubit number for local Hamiltonians [12]. Due to the similar structure of our ansatz in equation (10), one may expect comparable scaling for our variational method. This suggests that the upscaling procedure we employ only causes a degrading of approximation accuracy that scales linearly in the system size. This expectation is confirmed by our numerical experiments presented below.

4.2. Open boundary conditions

On some quantum hardware, such as superconducting circuit devices, the implementation of periodic boundary conditions on two-dimensional lattice models requires additional SWAP gates. Hence, the simulation of physical systems with open boundary conditions is better suited for such devices, particularly under NISQ conditions where each additional operation exponentially decreases the success probability. To adapt the upscaling scheme to systems with open boundary conditions, additionally to a small system with periodic boundaries, a small system of size N_o with open boundaries is variationally optimized. The corresponding Hamiltonian has the same form as (4) but without terms reaching beyond the periodic boundary. An intuition for why this should work is that if the depth of the circuit is short enough, the bulk qubits of the open system are affected by approximately the same operations as for a periodic boundary system. This is why we are allowed to glue an arbitrary periodic system into the middle of an open boundary system.

Time evolution can be approximated by the operators ${U_{\text{trot}}} {}_{, o}$ and ${U_{\text{var}}} {}_{, o}(\vec \theta)$ in an analogue way as before. The upscaled system is constructed from both the small periodic and non-periodic one, by cutting the small open boundary system of size N_o in half and inserting an arbitrary number K of periodic boundary qubits in the middle (see figure 1(b)). The gluing of the gates for each time step works as described before. The final circuit for time evolution of a single time step (m = 1) for the glued system with $N_o + K$ qubits with open boundaries reads

$$\begin{align} {U_{\text{var}, o}^{(1, N_o + K)}} \approx \, & {U_{\text{var}, o}^{(1, \frac{N_o}{2})} }(\vec \theta_{o}) \, {U_{\text{var}}}(\vec \theta) \, {U_{\text{var}, o}^{(K + \frac{N_o}{2} + 1, K+N_o)}} (\vec \theta_{o}) \end{align} \tag{ 12 }$$

where we used ${U_{\text{var}}}(\vec \theta)$ from (10) for the time evolution of the bulk. ${U_{\text{var}, o}^{(1, X)}}$ denotes a circuit simulating the open boundary dynamics on qubits 1 to X. Furthermore, we split the optimized parameter set θ_o of the open boundary N_o qubit simulation in half in order to use them for the time evolution of the left (or respectively right) boundary. For the sake of simplicity, we assume N_o to be even. This way, the three unitaries in (12) correspond to the three segments of figure 1(b) that are being glued. The boundary parameters θ_o and the bulk parameters θ result from two separate optimizations. We are thus making an approximation neglecting open boundary effects on the bulk qubits. We will see in numerical experiments that for a one- and two-dimensional TFIM it suffices to optimize only few dedicated boundary gates separate from the periodic bulk.

To numerically test the performance of the variational gluing method in comparison to Trotter–Suzuki approaches, the TFIM in one- and two-dimensional is studied. For this, we simulate the exact time evolution ${U_{\text{ex}}}$ generated by the Hamiltonian

$$\begin{align} H_\mathrm{TFIM} = J_z \sum_{\langle i, j \rangle} Z_i Z_j + h_x \sum_i X_i, \end{align} \tag{ 13 }$$

with the Pauli operators $X, Z$ and the coupling constants J_z and h_x . We denoted a Pauli Z operator acting on the $i\text{th}$ qubit by Z_i (and analogously for X). Furthermore, the notation $\langle i, j \rangle$ only counts nearest neighbor terms and we implicitly understood periodic boundary conditions. The transverse field Ising Hamiltonian relates to our previous notation via $H_{n, a = 1} = Z_n Z_{n+1}, H_{n, a = 2} = X_n$ and $c_1 = J_z, c_2 = h_x$. To demonstrate our method on another model, we show additional results for an XY model with a transverse field in appendix C.

The cost values for the Trotter–Suzuki ${U_{\text{trot}}}$ and variational sequence ${U_{\text{var}}}$ are compared for different system sizes, as well as for multiple repetitions of the circuit (i.e. multiple time steps). To emphasize that near term NISQ hardware can benefit from our classical pre-training algorithm for time evolution, simulation results of two-point spin correlations including gate infidelities are presented.

We simulated and optimized systems with a total number of up to 12 qubits in the discussed fashion. For the study of larger systems limitations of memory and time resources forced us to approximate the cost function (8) by sampling over a set of random input vectors $\mathcal{V}$

$$\begin{align} C_\mathrm{var}(\vec \theta) &\approx \frac{1}{|\mathcal{V}|} \sum_{\psi \in \mathcal{V}} \left( 1 - |\langle \psi| {U_{\text{ex}}^\dagger} {U_{\text{var}}}(\vec{\}theta) | \psi\rangle| \right), \end{align} \tag{ 14 }$$

where the absolute value $|\cdot|$ ensures the equality of the two operators up to a multiplicative global phase. This way, we are able to present results for system sizes up to 24 qubits. For the circuit simulations we used the Python library Cirq [40] and for optimization we used the Adam gradient descent method [41]. In all the optimizations, we have used a maximal number of training steps of 10⁴, a learning rate of 10⁻⁴ and the hyperparameters $\beta_1 = 0.9, \beta_2 = 0.99$ to control the decay rates of the gradient and squared gradient average.

There is an ambiguity in the choice of the ordering of the gate sequences which correspond to $Z_i Z_j$ and X_i . This leads to different operators when doing Trotterization. Here, we fixed the ordering for the sake of an efficient implementation to

$$\begin{align} {U_{\text{var}}}(\vec{\}theta) &= \prod_{r = 1}^m \Bigg[ \prod_{\substack{n = 1 \\ n \text{odd}}}^{N} e^{-i \theta_{r, 1} Z_{n} Z_{n+1}} \prod_{\substack{n = 1 \\ n \text{even}}}^{N} e^{- i \theta_{r, 1} Z_{n} Z_{n+1}} \prod_{n = 1}^N e^{-i \theta_{r, 2} X_n } \Bigg]. \end{align} \tag{ 15 }$$

Each of these gates can be implemented directly on current superconducting hardware via tunable couplers [42–46]. In this numeric example, four cases are studied. These are one- and two-dimensional systems with either periodic or open boundary conditions. In the case of two-dimensional periodic lattices, the gluing can be performed in either of the two dimensions by adding a row or column of qubits analogously to the discussion in section 4. Demonstration of scaling up in system size on two-dimensional systems with two open boundaries is computationally prohibitive, as the smallest fully open, two-dimensional and non-trivially glued system would be of size $5\times7$ (cf. appendix B). We therefore employ open boundaries in only one direction, while retaining periodic boundary conditions for the orthogonal direction. Optimal parameters can be classically found but no measure of accuracy can be evaluated any longer as we exceed our limit of numerical feasibility with a total size of 35 qubits.

5.1. Single time step

To begin with, our variational circuit is compared to Trotter–Suzuki approaches with equal or higher gate count and depth, where higher gate count is usually traded against a better approximation. Afterwards, the optimal parameters are reused in larger systems and also compared to the Trotter–Suzuki benchmarks. Figures 2(a) and (d) show that the variational ansatz is able to beat Trotterizations with a higher Trotter number along all numerically feasible system sizes (while being less interesting as it requires more gates, higher order Trotter–Suzuki approximations can also be beaten, see appendix A). Notably, the cost function values for upscaled systems hardly differ from the systems on which the optimization is performed. For the one-dimensional system, optimization was performed on a 6-qubit chain. An upscaling to 24 qubits only raises the cost value by a factor of 3. Analogously, for a two-dimensional system, the optimal parameters found for a $3 \times 3$ grid are extended to a $4\times 6$ grid by increasing the cost value by around the same factor. The cost values for the two-dimensional models are generally higher since each qubit has more interaction partners.

A linear fit in system size gives an estimate of how the improvement develops for system sizes which can no longer be classically simulated. In all studied cases, we observe an improvement factor in cost of the variational method compared to Trotter–Suzuki of more than 1000 in the one-dimensional and more than 100 in the two-dimensional case while using only half of the number of gates. The linear fits indicate that these improvement factors are constant over several larger system sizes.

NISQ gates have finite infidelities, leading to large compound errors when gate counts are high. It is therefore important to be as gate efficient as possible when running Hamiltonian simulation on NISQ hardware. To take into account gate infidelities, we quantify the total infidelity of a simulated NISQ experiment by

$$\begin{align} 1 - \mathscr{F}_\mathrm{NISQ} = 1 - \langle \mathscr{F}_\mathrm{approx.}\rangle \mathscr{F}_\mathrm{gate}, \end{align} \tag{ 16 }$$

where $\langle\mathscr{F}_\mathrm{approx}\rangle$ denotes the averaged fidelity of the output state (variational or Trotter) when the hardware runs without error and $\mathscr{F}_\mathrm{gates}$ is the average fidelity of the implementations of the employed gates. $\langle\mathscr{F}_\mathrm{approx}\rangle$ is inferred from the cost function used for optimization, where (8) takes an average over the whole Hilbert space while (14) only takes an average over a randomly chosen subset. For the gate fidelity $\mathscr{F}_\mathrm{gate}$, we consider a simple error model assuming the same fidelity p_g for each individual gate. The circuit fidelity is given as the product of the individual gate fidelities for every gate in the circuit, i.e. $\mathscr{F}_\mathrm{gate} = p_g^G$.

The total number of gates G, in general, is determined by the Trotter number m. In particular, for the Trotter–Suzuki decomposition of a TFIM with periodic boundaries, depth D and gate count G of a single time step read

$$\begin{align} D = (d + 1) \, m, \qquad \qquad G = (d+1) \, m N - bdm = (N-b)D + bm, \end{align} \tag{ 17 }$$

where d denotes the spatial dimension and b is 0 for periodic or 1 for open boundary conditions. We assume an optimistic gate fidelity of $p_g = 99.9\%$ [47], so that our results should remain relevant for future years.

The improvement in $\mathscr{F}_\mathrm{var}$ and hence in $\mathscr{F}_\mathrm{NISQ}$ strongly depends on which single time step size τ is chosen. For the Trotter–Suzuki approximation, a small ratio of τ and m is favored to keep the error small (see (4)). This can be done by increasing m and paying the price of using more gates. Figures 2(b) and (e) show that for small τ the gate infidelity dominates the error, whereas for large τ, the algorithmic error takes over. For the one-dimensional case, Trotter–Suzuki error becomes important at around τ = 1.5 for m = 3 and at τ = 2.0 for m = 6. The variational algorithm has the latest turning point, where the algorithmic error becomes dominant, around τ = 2.5. This turning point is earlier (around τ ≈ 1.5) in the two-dimensional model.

If we choose the single time step at this turning point, we find that the variational sequence increases the lower infidelity bound due to gate infidelities only marginally, while the Trotter–Suzuki decompositions introduce a larger additional error. For example, with 50 qubits and periodic boundaries in one dimension the NISQ infidelity $1-\mathscr{F}_\mathrm{NISQ} \approx 29 \%$ is increased by only $7 \%$ due to the variational approximation infidelity, whereas $1 - \mathscr{F}_\mathrm{NISQ} \approx 57 \%$ for the best (m = 6) Trotter–Suzuki ansatz. Both, a higher gate count and approximation infidelity lead to this deterioration. The standard deviation σ of the linear $\mathscr{F}_\mathrm{approx.}$ extrapolations in figures 2(c) and (f) are small compared to $\mathscr{F}_\mathrm{approx.}$: $\sigma = 0.19 \% \cdot \mathscr{F}_\mathrm{approx.}$ (1D) and $\sigma = 0.03 \% \cdot \mathscr{F}_\mathrm{approx.}$ (2D). Therefore, the fits in figure 2 yield a faithful extrapolation in system size.

5.2. Multiple time steps

Since we determine ${U_{\text{var}}}$ via a state agnostic optimization, it can be applied repeatedly to the time evolved state. In this way larger simulation times T can be covered. Figure 3 show this for a 6- (3(a)–(c)) and a 15-qubit system (3(d)–(f)) both using the parameters optimized on 6 qubits. The variational sequence is optimized for a single time step, and improvements in the cost function decrease for later times. Despite this, the variational algorithm remains better than the Trotter–Suzuki benchmark even after one hundred repetitions for most values of the Hamiltonian coefficients J_z and h_x . However, when one of these coefficients is small the Trotter error decreases, which can lead to Trotter–Suzuki overtaking the variational algorithm in certain regions of parameter space. Indeed, the Trotter–Suzuki approximation becomes exact in the limit that either J_z or h_x goes to zero. In extreme cases, the accuracy of the Trotter–Suzuki sequence exceeds that of the variational algorithm after 50 repetitions. However, as we shall show later, at 50 repetitions the total error is dominated by gate infidelities rather than approximation error.

Zoom In Zoom Out Reset image size

We obtain a similar behavior of the cost function for the case of d = 2 and open boundary conditions. This indicates that the open boundary variation of our method can be applied without significant change of performance. In this numerical demonstration, our variational algorithm exhibits even a better approximation fidelity for half-open boundaries than for periodic boundaries, and hence a fully translational invariant system. This further stresses the feasibility of our method for NISQ-inspired open boundary quantum systems.

5.3. Precision of time-evolved observables

When computing the quantum dynamics of a system, one is typically interested in the evolution of physical observables. From (9), it is clear that the cost values only yield an upper bound for the error in predicting observables. Although the variational algorithm has a smaller cost value, Trotter–Suzuki could thus, in principle, describe the dynamics of a specific observable for a specific initial state with greater accuracy. The approximation error in the expectation value of a time-evolved observable hence depends on the precise observable and on the chosen initial state. In this section, we show results for exemplary observable dynamics and an initial state.

Note that the gluing procedure is not sensitive to entanglement of the input state, since the time evolution operator is agnostic of the initial state. States of high entanglement (or strongly entangling dynamics) are particularly interesting, as the Trotter error of localized physics has been found to be independent of simulation time [48].

In figure 4, we illustrate the relative error on two-point correlation of two neighboring spins on a 15 qubit system ($3 \times 5$ for d = 2) far away from the boundary, $\mathcal{O} = \langle \psi(t)|Z_8 Z_9|\psi(t)\rangle$ ($\mathcal{O} = \langle \psi(t)|Z_{(2,3)} Z_{(2,4)}|\psi(t)\rangle$ for d = 2). We find that for the initial state $\vert \psi(0)\rangle = \vert +\rangle^{\otimes N}$, the two-point correlation $\mathcal{O}$ predicted by ${U_{\text{var}}}$ lies closer to the exact dynamics than any of the Trotter–Suzuki benchmarks.

Taking a product state as initial state is particularly interesting, as it is accessible in experiment, yet also yields an entangled state after a few repetitions. Although the m = 3 variational circuit and the m = 6 Trotter–Suzuki benchmark achieve similar approximation error, the latter contains twice as many gates and thus suffers from a squared gate fidelity.

These imperfections are illustrated in figures 4(c) and (d), where we show the overall infidelity of (16), for up to 20 repetitions of the single time step τ = 1.0. Besides the variational sequence beating both Trotter–Suzuki benchmarks, the overall infidelity is minimized even more when optimizing for a larger time step 2τ using m = 6 layers. Despite a significant noise floor set by the gate infidelities, we achieve an improvement of infidelity by a factor between 3 (d = 2) and 4 (d = 1) for a single repetition. Since the best Trotter–Suzuki decomposition starts with an infidelity of $44.5 \%$, this shows that the simulation time $2\tau = 5.0$ cannot yield satisfactory results in the sketched NISQ scenario, while the variational sequence does (infidelity of $18.8\%$). Since also twice as many parameters are involved, it is expected that the variational sequence optimized on longer times performs better. This comes at the cost of a more difficult optimization task in the classical pre-processing step. The trade-off between more time spent once in classical pre-processing against a permanently better quantum circuit is very favorable on NISQ devices. Generalizing this to longer simulation times $k\tau \gt 5.0$ and more parameters m > 6, however, does not further improve the long time performance and thus is omitted in figure 4.

Zoom In Zoom Out Reset image size

All the numerical experiments support the maxim to keep the number of repetitions low while using as large time steps τ as possible. To do this, an increase of the Trotter number for the variational single time step m might be worth evaluating (cf. figure 4). For small τ (as τ = 0.3 in figures 2 and 3, for instance), raising the Trotter number m does not significantly improve the cost values. Therefore, the data points for m > 3 are omitted.

Also remarkably, the (half-)open boundary systems yield even smaller infidelity values than their periodic counterparts. Besides supporting the statement that open boundaries can be glued to get NISQ-friendly systems, this phenomenon can be explained, since in open boundary systems fewer gates contribute to both approximation and gate error.

In this work, we demonstrated how classical pre-processing can be leveraged to determine a variational representation of time evolution operator of a translation symmetric quantum system that reliably exceeds comparable Trotterization in approximation accuracy. The main idea for our variational algorithm consists of classically optimizing a variational time evolution ansatz circuit on a small tractable system and then reusing these optimal parameters for larger, classically intractable systems by gluing gate sequences together.

These results are particularly interesting due to their applicability to NISQ devices. Most promisingly, it can avoid any sampling cost and limitation associated with variational optimization on quantum hardware. Further, the significant increase in approximation fidelity compared to Trotter–Suzuki with the same gate count for practically relevant system sizes can be taken advantage of to realize the same approximation fidelity at a notably reduced gate count. We found that, for one- and two-dimensional example problems, the observed advantages persist and are in fact enhanced for systems with open boundaries, despite the loss of translational symmetry on the edges. Such open boundary systems are important as they are easier to realize on NISQ hardware.

Our results open up a highly versatile method for variational quantum simulation. Since we have shown that it is sometimes sufficient to maximize the simulation quality locally, possible limitations of this gluing scheme are worth investigating further. We also expect to gain further insights into the scope and performance of our method from investigations of the relations between the range of correlations in simulated quantum states and the size of the system considered in classical optimization. Since the optimization is however state agnostic, depending only on the Hamiltonian of the system, it is expected to perform well in simulating the time evolution of highly-entangled states with long-ranged correlations.

Another way of optimizing quantum simulation algorithms can be achieved by a quantum optimization. The cost function of (14) can be measured by comparing initial and output state of a circuit which first performs the variational ansatz circuit which propagates forward in time. Afterwards the state is propagated back to its initial state by a high order Trotter–Suzuki approximation. This way, the variational circuit can be optimized to invert the Trotter–Suzuki decomposition with fewer gates and later scaled up in system size. Further investigation is needed to determine which system sizes are suitable for optimization.

Randomization techniques have attracted much attention recently and have been shown to compress Trotter circuits [21, 24]. It remains to be investigated whether randomized ansätze can be used for variational algorithms like the one presented in this work.

The data that support the findings of this study are available upon reasonable request from the authors.

We thank Petr Zapletal for valuable feedback and discussions and Naeimeh Mohseni for feedback on the manuscript. This work received funding from the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 828826 'Quromorphic.', from the German Federal Ministry of Education and Research via the funding program quantum technologies—from basic research to the market under contract number 13N15577 'MANIQU'. It is also part of the Munich Quantum Valley, which is supported by the Bavarian state government with funds from the Hightech Agenda Bayern Plus. This work was partially supported by NHR@FAU, which is funded by the State of Bavaria and by the Federal Ministry of Education and Research. TE acknowledges funding from the International Max Planck Research School Physics of Light.

No new data were created or analysed in this study.

A more accurate decomposition as in (5) can be achieved by higher order product formulas. A second order Trotter–Suzuki decomposition is achieved by only reordering the operators in a symmetric way. For a single Trotter step (m = 1) this reads

$$\begin{align} {U_{\text{ex}}} &= \left[ \prod_{n, a}^\leftarrow \exp{\left( \frac{- i \tau}{2} c_a H_{n, a} \right)} \prod_{n, a}^\rightarrow \exp{\left( \frac{- i \tau}{2} c_a H_{n, a} \right)} \right] + \mathcal{O}\left( \tau^3 \right) =: S_2(\tau), \end{align} \tag{ A.1 }$$

Furthermore, one can achieve arbitrarily accurate approximations of an exponential $e^{-i\tau H}$ by applying the qth order Trotter–Suzuki formula which is defined recursively by [4]

$$\begin{align} S_{q}(\tau) = [S_{q - 2} (\nu_q \tau)]^2 S_{q - 2} ((1 - 4 \nu_q) \tau)[S_{q - 2} (\nu_q \tau)]^2, \quad \text{for } q > 2 \end{align} \tag{ A.2 }$$

with $\nu_q = \left( 4 - 4^{\frac{1}{q - 1}} \right)^{-1}$. In addition to (A.2), one can also pick a Trotter number m > 1 to raise the accuracy further. Employing a higher m is achieved by substituting $\tau \mapsto \frac{\tau}{m}$ in (A.2), but repeating the sequence m times, thus using more gates by a factor of m.

We will stick to the $m \in \{1, 2\}$ case here and show that the use of higher order formulas as a benchmark for our method is practically not relevant for NISQ-focused variational optimization, since the number of required gates, which scales with 5^q [13], will dominate the error for finite gate infidelities. We provide results for $q \in \{2, 4, 6, 8\}$.

The overall circuit fidelity is given by the product of the imperfect gate fidelity $\mathscr{F}_\mathrm{gate}$ and the fidelity due to the approximate algorithm

$$\begin{align} \mathscr{F}_\mathrm{var}(\vert \psi\rangle) = \langle \psi |U_\text{ex}^\dagger U_\text{var} | \psi\rangle \end{align} \tag{ A.3 }$$

for an arbitrary initial state $\vert \psi\rangle$. Analogously, one can define the algorithmic fidelity $\mathscr{F}_\mathrm{trot}$ of a Trotter–Suzuki algorithm described by the product formula in (A.2). We combine the notations $\mathscr{F}_\mathrm{var}$ and $\mathscr{F}_\mathrm{trot}$ into $\mathscr{F}_\mathrm{approx.}$. Figure A1 shows the total infidelity $1 -\langle\mathscr{F}_\mathrm{approx.}\rangle\cdot\mathscr{F}_\mathrm{gate}$ for up to 100 time steps. One can see that while second order Trotter–Suzuki would admit fewer gates than the variational algorithm (m = 3), it suffers from greater algorithmic errors. Higher order Trotter–Suzuki algorithms are already dominated by gate infidelities after a single time step. The variational algorithm hence outruns all of the benchmark algorithms in the presented simulation problem. For this reason, no higher order q Trotter–Suzuki formulas are discussed in the main text. Other product formulas [49, 50] behave very similar to Trotter in this respect. As Trotter exhaust theoretical complexity lower bounds for higher orders [12], we only discuss higher order Trotter formulas here.

Zoom In Zoom Out Reset image size

We construct the smallest non-trivially glued fully open boundary system in two dimensions. It turns out that, although optimal parameters to use on fully open two-dimensional systems can be found, no measure of accuracy can be numerically evaluated. In figure B2 every colored square corresponds to a qubit and thick lines represent open boundaries. As every qubit contributes to the Hamiltonian with a number of A interactions, the squares can be thought of as a set of $A \cdot m$ optimal parameters and the figure depicts how parameters are copied when scaling up. The blue squares all share the same parameters for each gate.

Zoom In Zoom Out Reset image size

Optimized open boundary systems have to be cut in half to insert bulk qubits in between (cf. section 4). In two dimensions however, this means the system needs to be cut along both dimensions. Then first, a half-open system is inserted to faithfully keep the open boundary upright (orange blocks in figure B2). Finally, completely periodic qubits can be filled in the very bulk of the system (blue blocks), but also the green blocks have added at least one column of qubits to keep everything rectangular. The additional green column copies parameters from the inner most qubit of one half of the boundary (for one column it does not matter whether it copies the parameter from the left or the right). This way, the inner most qubit of each half of the fully open system has to be considered equivalent to a qubit from a half-open system. If we would not choose at least 4 qubits along this open boundary, at least one boundary qubit would coincide with a bulk qubit which is contradictory.

If there is a periodic boundary along one dimension nothing has to be cut in half, however, the optimization should treat every qubit such that it has two distinct neighbors along the periodic dimension. This can only be achieved with at least three qubits, as otherwise left and right neighbors would coincide. As a result, the smallest periodic boundary block is of size $3\times 3$, the smallest half-open block is $3 \times 4$ qubits (3 along the periodic boundary and 4 along the open boundary) and the smallest fully open block is of size $4 \times 4$.

Thus, the smallest glued open boundary system in two dimensions is of size $5 \times 7$ consisting of three optimized blocks, a half-open one of size $3 \times (2+2)$, a fully open block of size $5 \times (2+2)$ (the optimized block is of size $(2+2) \times (2+2)$) and a row of 3 qubits in the middle, see figure B2. The optimization of the single elements can all be performed on a classical computer, but the cost function for a $5 \times 7$ qubit system can no longer be evaluated, as we exceed our limit of numerical feasibility. The results for a half-open system of size $3 \times (4+K)$ can be evaluated and hence represent a numerical example of gluing open boundaries to two-dimensional spin lattice models.

Next to the study of the TFIM, in this appendix the XY model with a transverse field (TFXY) is discussed using similar quality measures as in section 5. Coming from the TFIM, we just add a YY interaction to neighboring qubits and get an XY model with an external field. The Hamiltonian reads

$$\begin{align} H_\mathrm{TFXY} = J_y \sum_{\langle i, j \rangle} Y_i Y_j + J_z \sum_{\langle i, j \rangle} Z_i Z_j + h_x \sum_i X_i. \end{align} \tag{ C.1 }$$

Since this model includes more non-commuting terms, we expect the Trotter error to be larger, making it another interesting model for NISQ simulation. This claim is supported by figure C3 which compares cost values of variational and Trotter–Suzuki sequence scaled in system size and in simulation time, as well as errors of observable dynamics.

Zoom In Zoom Out Reset image size

Figures C3(a) and (b) show that the first order Trotter–Suzuki approximation has a larger cost value compared to the TFXY model with the same interaction coefficients as in section 5 except adding a YY interaction with $J_y = 0.5$. As before, the improvement from optimization of a single time step is larger in one-dimensional than two-dimensional systems and also the scaling in system size, resembles the behavior of the TFIM well. The extrapolation of the cost value in systems size keeps up an improvement in cost value of a factor of 1000 (100 for d = 2) compared to Trotter–Suzuki with Trotter number m = 3. For a one-dimensional system with 50 qubits, for instance, the Trotter–Suzuki decomposition cost value is extrapolated to around 0.01 for m = 3 and $5 \cdot 10^{-3}$ for m = 6. The variational sequence has a cost value of 10⁻⁵. For a two-dimensional grid of 50 qubits, in total, the cost values are of the order of 10⁻¹ (Trotter m = 3), 10⁻² (Trotter m = 6) and 10⁻³ (variational m = 3).

The scaling in simulation time is shown in figures C3(c)–(f). The general trend of the plots again resembles the case of the TFIM which indicates a certain model independence of the variational method. The slight increase in approximation error in the XY model however also reduces the number of circuit repetitions which can be performed until the variational sequence falls back to accuracies of Trotter–Suzuki. For an optimization on m = 3 layers, this happens after around 40 repetitions in the one-dimensional, and 20 repetitions in the two-dimensional case. To push this limit to longer times, an optimization on double the time step $2\tau = 0.6$ and double the layers m = 6 can be performed. The improvement gained for larger time steps is traded off against a more complex classical optimization algorithm, here (cf. figure 4 and discussion).

The simulation of observable dynamics shall complete our excursion to the XY model. In figures C3(e) and (f), the relative error of a two-point correlation along 50 repetitions is shown. The initial state $\psi = \vert +\rangle^{\otimes 15}$ and observable $\langle Z_8 Z_9\rangle$ are the same as in the analogue discussion for the TFIM. What is different is the single time step τ = 0.3 and that we also compare the variational sequence optimized on m = 6 layers. Although the cost values of the two different variational sequences show a clear hierarchy for multiple time steps, there is no clear winner identifiable when it comes to the approximation of the observable dynamics. For the two-dimensional system, the m = 6 Trotter–Suzuki sequence sometimes yields better approximations to the chosen observable, whereas it is beaten in all repetitions in the one-dimensional case. This emphasizes the volatility of observable errors, although also here an advantageous behavior of the variational algorithms is visible. Note that still m = 6 Trotter–Suzuki uses double the number of gates than the other sequences which needs to be taken into account in the case of imperfect gate fidelities.