Nearly optimal lattice simulation by product formulas

TL;DR

This paper proves that product formulas can simulate lattice Hamiltonians on quantum computers with nearly optimal gate complexity, providing rigorous error bounds and extending previous results.

Contribution

It offers a rigorous proof of the nearly optimal gate complexity for lattice Hamiltonian simulation using product formulas, including error analysis and generalizations.

Findings

Gate complexity is (nt)^{1+o(1)} for lattice Hamiltonian simulation.

Provides error bounds for canonical product formulas including Lie-Trotter-Suzuki.

Extends analysis to time-dependent Hamiltonians and higher dimensions.

Abstract

We consider simulating an $n$-qubit Hamiltonian with nearest-neighbor interactions evolving for time $t$ on a quantum computer. We show that this simulation has gate complexity $(nt)^{1+o(1)}$ using product formulas, a straightforward approach that has been demonstrated by several experimental groups. While it is reasonable to expect this complexity---in particular, this was claimed without rigorous justification by Jordan, Lee, and Preskill---we are not aware of a straightforward proof. Our approach is based on an analysis of the local error structure of product formulas, as introduced by Descombes and Thalhammer and further simplified here. We prove error bounds for canonical product formulas, which include well-known constructions such as the Lie-Trotter-Suzuki formulas. We also develop a local error representation for time-dependent Hamiltonian simulation, and we discuss generalizations to periodic boundary conditions, constant-range interactions, and higher dimensions. Combined with a previous lower bound, our result implies that product formulas can simulate lattice Hamiltonians with nearly optimal gate complexity.

Equations

$${\mathcal{S}}_1(t):$$

$$=e^{-it{H}_{\text{even}}}{e}^{-it{H}_{\text{odd}}}$$

$${\mathcal{S}}_2(t)$$

$${\mathcal{S}}_{2k}(t)$$

$$\Vert {\mathcal{S}}_1(t)-e^{-itH}\Vert$$

$$\Vert {\mathcal{S}}_{2k}(t)-e^{-itH}\Vert$$

$${\mathcal{S}}_1^{\prime }(t)$$

$$\int_{0}^{t}\mathrm{d}\tau_{1}\ e^{-i(t-\tau_{1})H}\big{[}e^{-i\tau_{1}H_{\text{even}}},\ -iH_{\text{odd}}\big{]}e^{-i\tau_{1}H_{\text{odd}}}.$$

$$\displaystyle\left\lVert\mathscr{S}_{1}(t)-e^{-itH}\right\rVert\leq\int_{0}^{t}\mathrm{d}\tau_{1}\int_{0}^{\tau_{1}}\mathrm{d}\tau_{2}\left\lVert\big{[}H_{\text{even}},\ H_{\text{odd}}\big{]}\right\rVert.$$

$$\displaystyle\big{[}H_{\text{even}},\ H_{\text{odd}}\big{]}=\sum_{k=1}^{\frac{n}{2}}\big{[}H_{2k-2,2k-1}+H_{2k,2k+1},\ H_{2k-1,2k}\big{]}.$$

$$\displaystyle\left\lVert\big{[}e^{-itH_{k,k+1}},e^{-itH_{l,l+1}}\big{]}\right\rVert\leq 2t^{2}$$

$$\Vert \prod _{j=1}^{n-1}{e}^{-it{H}_{\sigma (j),\sigma (j)+1}}-e^{-it{H}_{\text{even}}}{e}^{-it{H}_{\text{odd}}}\Vert$$

$$+\Vert e^{-it{H}_{\text{even}}}{e}^{-it{H}_{\text{odd}}}-e^{-itH}\Vert .$$

$${\mathcal{S}}_{2k}(t)-e^{-itH}={\int }_0^t{e}^{-i(t-\tau )H}{\mathcal{S}}_{2k}(\tau )\mathcal{T}(\tau )\mathrm{d}\tau ,$$

$$\mathcal{T}(\tau )=2k{\int }_0^1\mathrm{d}x(1-x{)}^{2k-1}{\mathcal{T}}^{(2k)}(x\tau )\frac{{\tau }^{2k}}{(2k)!}.$$

$$H_{\text{odd}}(t)$$

$$H_{\text{even}}(t)$$

$$\left\lVert\mathscr{S}_{\mathcal{T},2k}(t)-\exp_{\mathcal{T}}\bigg{(}{-}i\int_{0}^{t}\mathrm{d}v\ H(v)\bigg{)}\right\rVert=O\big{(}nt^{2k+1}\big{)}$$

$$H_{[1]}=$$

$$H_{[2]}=$$

$$⋮$$

$$H_{[\ell ]}=$$

Full text

A natural class of Hamiltonians that includes many physically reasonable systems is the class of lattice Hamiltonians (6, 11, 12, 13). Lattice Hamiltonians arise in many models of condensed matter physics, including systems of spins (e.g., Ising, XY, and Heisenberg models; Kitaev’s toric code and honeycomb models; etc.), fermions (e.g., the Hubbard model and the $t$-$J$ model), and bosons (e.g., the Bose-Hubbard model). Note that fermion models can be simulated using local interactions among qubits by using a mapping to qubits that preserves locality (14). Digital simulations of quantum field theory also typically involve approximation by a lattice system (6).

For simplicity, we mainly focus on nearest-neighbor lattice systems in one dimension (although we discuss generalizations to other lattice models as well). In this case, $n$ qubits are laid out on a one-dimensional lattice and the Hamiltonian only involves nearest-neighbor interactions. Specifically, a Hamiltonian $H$ is a lattice Hamiltonian if it acts on $n$ qubits and can be decomposed as $H=\sum_{j=1}^{n-1}H_{j,j+1}$, where each $H_{j,j+1}$ is a Hermitian operator that acts nontrivially only on qubits $j$ and $j+1$. We assume that $\max_{j}\left\lVert H_{j,j+1}\right\rVert\leq 1$, for otherwise we evolve under the normalized Hamiltonian $H/\max_{j}\left\lVert H_{j,j+1}\right\rVert$ for time $\max_{j}\left\lVert H_{j,j+1}\right\rVert t$.

Lloyd’s original proposal for an explicit quantum simulation algorithm (15) uses the Lie-Trotter product formula. Subsequent work achieves better asymptotic complexity (16) using higher-order Suzuki formulas (17). We refer to all such formulas as product formulas. The product-formula algorithm is straightforward yet surprisingly efficient for quantum simulation. Indeed, it can conserve certain symmetries of the dynamics (18) and appears to be advantageous for various practical systems (19, 2, 20). Although recent simulation algorithms have better asymptotic complexities (21, 22, 23, 24, 25, 26, 27, 28, 29, 30), the product-formula approach remains a natural choice for experimental simulations (31, 32, 33) due to its simplicity and the fact that it does not require any ancilla qubits. Its study has also illuminated areas beyond quantum computing (34).

One of the main challenges in quantum simulation is to analyze the gate complexity of simulation algorithms. Explicit gate counts are especially desirable for near-term simulation because early quantum computers will only be able to reliably perform a limited number of gates. While existing analysis appears to be tight for recent simulation algorithms, the product-formula bound can be loose by several orders of magnitude (19, 8, 20, 2). This dramatic gap makes it hard to identify the fastest simulation algorithm and to find optimized implementations for near-term applications (19).

Product formulas can simulate a lattice system with fixed accuracy with gate complexity O\big{(}n(nt)^{2}\big{)} in the first-order case and O\big{(}n(nt)^{1+\frac{1}{2k}}\big{)} in the ($2k$)th-order case. However, it is natural to expect a more efficient simulation. Roughly speaking, a system simulates its own evolution for constant time using only constant circuit depth—and hence an extensive number of gates—so one might expect a true simulation complexity of $O(nt)$. Indeed, Jordan, Lee, and Preskill claimed that product formulas can simulate an $n$-qubit lattice system with $(nt)^{1+o(1)}$ gates (6), but they did not provide rigorous justification and it is unclear how to formalize their argument. Subsequent work improves the analysis of the product-formula algorithm using information about commutation among terms in the Hamiltonian (35, 19, 4, 3), the distribution of norms of terms (36), and by randomizing the ordering of terms (37, 38). However, none of these improvements can achieve the claimed gate complexity $(nt)^{1+o(1)}$ for lattice simulation.

Main result. Let $H=\sum_{j=1}^{n-1}H_{j,j+1}$ be an $n$-qubit lattice Hamiltonian. We order the terms in the even-odd pattern $H_{1,2},H_{3,4},\ldots,H_{2,3},H_{4,5},\ldots$ obtaining the first-order product formula

[TABLE]

and the ($2k$)th-order product formulas

[TABLE]

with $p_{k}:=1/(4-4^{1/(2k-1)})$. Our main result is an asymptotic upper bound on the product-formula error

[TABLE]

where $\left\lVert\cdot\right\rVert$ denotes the spectral norm.

The above error bound works well only for very small $t$. To simulate for a longer time, we divide the entire evolution into $r$ segments, and within each segment, we simulate using product formulas. To achieve accuracy $\epsilon$, it suffices to choose $r_{1}=O(nt^{2}/\epsilon)$ for the first-order formula and r_{2k}=O\big{(}t(nt/\epsilon)^{\frac{1}{2k}}\big{)} for the ($2k$)th-order formula. Equivalently, we have gate complexity O\big{(}(nt)^{2}\big{)} and O\big{(}(nt)^{1+\frac{1}{2k}}\big{)} for the first- and ($2k$)th-order algorithm, assuming that we simulate with constant accuracy.

For any $\delta>0$, we choose an integer $k$ sufficiently large so that $\frac{1}{2k}\leq\delta$, upper-bounding the gate complexity as O\big{(}(nt)^{1+\frac{1}{2k}}\big{)}=O\big{(}(nt)^{1+\delta}\big{)}. This proves that the product-formula algorithm has asymptotic gate complexity $(nt)^{1+o(1)}$. Combining with the lower bound of $\widetilde{\Omega}(nt)$ established in (12), we have showed that product formulas can simulate a lattice Hamiltonian with nearly optimal gate complexity.

Applications. As an immediate application, our result gives a rigorous proof of the Jordan-Lee-Preskill claim about the complexity of simulating quantum field theory (6). Recent works have analyzed the gate complexity of other quantum field theory simulations (39), including digital simulation of gauge theories (40). The lattice Hamiltonians there have similar locality, so our analysis still applies. We expect our technique can be generalized to speed up the simulation of other systems, such as electronic structure Hamiltonians (9), power-law decaying interactions (41), exponentially decaying interactions (42), and clustered Hamiltonians (43), but we leave a thorough study of such generalizations as a subject for future work (44).

To simulate an $n$-qubit lattice Hamiltonian for time $t$, our algorithm has circuit depth $n^{o(1)}t^{1+o(1)}$. As a side application, our analysis gives a tensor network representation of lattice systems with bond dimension $2^{n^{o(1)}t^{1+o(1)}}$, using the counting argument of (45). This significantly improves a recent construction of (46, Lemma 17) which uses only the first-order Trotter decomposition.

We work primarily with an idealized setting where quantum operations can be performed faithfully. However, in realistic experiments, quantum gates will be imperfectly implemented (47). For such a case, Reference (48) estimates the simulation accuracy as $\frac{\alpha}{r^{2k}}+\beta r$ in diamond-norm distance (49, 50), where $\alpha$ captures the algorithmic error of product formulas and $\beta$ captures gate errors. This leads to an optimal number of segments r=\big{(}\alpha 2k/\beta\big{)}^{\frac{1}{2k+1}}, which can be improved by our result. Specifically, the original analysis in (16) implies \alpha_{\text{orig}}=O\big{(}(nt)^{2k+1}\big{)}. This has been improved by subsequent work (19, 37), although none of these improvements achieves linear scaling in $n$. In contrast, the analysis of this letter gives \alpha_{\text{opt}}=O\big{(}nt^{2k+1}\big{)}, improving the performance as a function of $n$ even in the presence of noise.

Our main goal is to establish the gate complexity of $(nt)^{1+o(1)}$ for the product-formula algorithm. However, our analysis is not only nearly optimal in the asymptotic regime but also appears to be much tighter in practice. For concreteness, we numerically implement and optimize our fourth-order bound, and compare it with previous product-formula analysis, for simulation of a one-dimensional Heisenberg model with a random magnetic field with open boundary conditions (51, Eq.(98)) (see Figure 1). We find that the scaling of our bound matches the empirical performance and the constant prefactor is off by only one order of magnitude, a significant improvement over previous rigorous bounds (19). Further improvements of our bound are possible by optimizing its numerical implementation; we leave a detailed study for future work (44).

Analysis of the first-order algorithm. The key technique behind our approach is an integral representation of the error $\mathscr{S}(t)-e^{-itH}$ that we develop based on Descombes and Thalhammer’s local error analysis of product formulas (52). In the local error representation, the integrand is expressed as a linear combination of commutators nested with unitary conjugations, where the numbers of summands and nesting layers are both independent of $n$ and $t$. We use this representation to get the correct asymptotic gate count as a function of $n$ and $t$. In contrast, the conventional approach uses the Baker-Campbell-Hausdorff formula or naive Taylor expansion, which requires the manipulation of infinite series and appears to be technically challenging to analyze (53) (51, Section I).

To illustrate the proof idea, we show how to obtain $\left\lVert\mathscr{S}_{1}(t)-e^{-itH}\right\rVert=O(nt^{2})$ for the first-order formula. We differentiate $\mathscr{S}_{1}(t)$ and obtain

[TABLE]

Using the variation-of-constants formula (52) (54, Theorem 4.9) with initial condition $\mathscr{S}_{1}(0)=I$, we find an integral representation of the product-formula error $\mathscr{S}_{1}(t)-e^{-itH}$ as

[TABLE]

We repeat this procedure to analyze the commutator \big{[}e^{-i\tau_{1}H_{\text{even}}},\ -iH_{\text{odd}}\big{]}, obtaining an upper bound on the spectral-norm error

[TABLE]

We expand $H_{\text{odd}}$ and $H_{\text{even}}$ according to their definitions. Fixing an arbitrary term $H_{2k-1,2k}$ in $H_{\text{odd}}$, the commutator \big{[}H_{2l,2l+1},\ H_{2k-1,2k}\big{]} is non-zero only when $l\in\{k-1,k\}$. We thus find that

[TABLE]

Using the triangle inequality, we have $\left\lVert\mathscr{S}_{1}(t)-e^{-itH}\right\rVert=O(nt^{2})$, which proves the claim (3) for the first-order case.

Ordering robustness. Our above bound works when terms of the lattice Hamiltonian are ordered in the even-odd pattern. However, this choice is not necessary: the first-order algorithm has the same asymptotic error bound with respect to any ordering of the lattice terms.

Our analysis relies on an error bound for swapping lattice terms:

[TABLE]

if $|k-l|=1$ and [math] otherwise. In words, we may swap two exponentials $e^{-itH_{k,k+1}}$ and $e^{-itH_{l,l+1}}$ without error unless their supports overlap, in which case the error is $O(t^{2})$.

Let $H=\sum_{j=1}^{n-1}H_{j,j+1}$ be a lattice Hamiltonian. We now simulate it using the first-order formula, but allow terms to be ordered arbitrarily as $\prod_{j=1}^{n-1}e^{-itH_{\sigma(j),\sigma(j)+1}}$, where $\sigma\in S_{n-1}$ is a permutation on the $n-1$ elements $\{1,\ldots,n-1\}$. Then the spectral-norm error $\left\lVert\prod_{j=1}^{n-1}e^{-itH_{\sigma(j),\sigma(j)+1}}-e^{-itH}\right\rVert$ is upper bounded by

[TABLE]

The second term is upper bounded by $O(nt^{2})$. For the first term, we transform $\prod_{j=1}^{n-1}e^{-itH_{\sigma(j),\sigma(j)+1}}$ into $\prod_{k=1}^{\frac{n}{2}-1}e^{-itH_{2k,2k+1}}\prod_{k=1}^{\frac{n}{2}}e^{-itH_{2k-1,2k}}$ by swapping neighboring exponentials. Every time two exponentials are swapped, we use (8) to bound the error. The total number of swaps of exponentials $e^{-itH_{k,k+1}}$ and $e^{-itH_{l,l+1}}$ with $|k-l|=1$ is at most $2n$, incurring error $4nt^{2}=O(nt^{2})$.

We have therefore obtained the same asymptotic error bound for an arbitrary ordering of the Hamiltonian terms. We call this phenomenon ordering robustness. Our analysis shows that the first-order algorithm is ordering-robust. Whether a similar property holds for a general higher-order formula remains an open question.

We also numerically compare the first-order algorithm with the even-odd ordering and the ordering of (19). Although they have the same asymptotic error bound, in practice the even-odd ordering has smaller exponent and constant prefactor. Details can be found in (51, Section VI).

Analysis of higher-order algorithms. Analyzing higher-order product formulas is more challenging. To this end, we represent them in a canonical form, which is easy to manipulate and encompasses well-known constructions such as the Lie-Trotter-Suzuki formulas $\mathscr{S}_{1}(t)$, $\mathscr{S}_{2k}(t)$ as special cases. We then use the variation-of-constants formula to write

[TABLE]

where \mathscr{T}(t)=\mathscr{S}_{2k}^{\dagger}(t)\big{[}\frac{\mathrm{d}}{\mathrm{d}t}\mathscr{S}_{2k}(t)-(-iH)\mathscr{S}_{2k}(t)\big{]}. As a ($2k$)th-order formula, $\mathscr{S}_{2k}(t)$ satisfies an order condition $\mathscr{S}_{2k}(t)=e^{-itH}+O(t^{2k+1})$, which further implies by Taylor’s theorem

[TABLE]

Canonical product formulas and their order conditions are further discussed in (51, Section II).

A direct expansion of $\mathscr{T}(t)$ gives the correct $t$-dependence $O(t^{2k+1})$ of the product-formula error, but the scaling in $n$ is incorrect. Instead, we seek an alternative expression for the integrand that consists of a linear combination of commutators nested with unitary conjugations, where the number of summands and nested layers are both independent of $n$ and $t$. Such an expression is referred to as a local error representation in (52). However, the result of (52) depends on auxiliary functions whose recursive structure is hard to unravel. Instead, we develop a simpler representation of the local error structure (51, Section III).

In our representation, the operator $\mathscr{T}(\tau)$ can be written as a linear combination of operator-valued functions of the form $e^{i\tau X_{1}}\cdots e^{i\tau X_{l}}Ye^{-i\tau X_{l}}\cdots e^{-i\tau X_{1}}$, where operators $X_{j},Y\in\{H_{\text{even}},H_{\text{odd}}\}$. As such, its higher-order derivatives consist of unitary conjugations and commutators. When a commutator is composed, we imitate (7) to show that the support of the operator is expanded by at most a constant factor. When a unitary conjugation is composed, we decompose the unitary operators and cancel exponentials with non-overlapping supports. Throughout this procedure, we only introduce $O(n)$ error in the innermost layer, proving the claim in (3) for the higher-order cases. This error analysis is discussed in more details in (51, Section IV).

Generalized lattice Hamiltonians. We have so far focused on time-independent one-dimensional systems with nearest-neighbor interactions and open boundary conditions. However, our analysis can be easily adapted to handle time-dependent Hamiltonians, periodic boundary conditions, constant-range interactions, and higher-dimensional systems, again with nearly optimal gate complexity.

When the Hamiltonian $H(t)$ is time-dependent, the problem of quantum simulation becomes more difficult (55). Then there no longer exists a closed-form solution to the Schrödinger equation. Furthermore, some quantum simulation algorithms (22, 25) that behave well in the time-independent case fail to handle the time-dependent Hamiltonian simulation. Nevertheless, we show that product formulas can simulate time-dependent lattice Hamiltonians with nearly optimal gate complexity. We group the terms in the even-odd pattern

[TABLE]

and simulate using the time-dependent Lie-Trotter-Suzuki formulas $\mathscr{S}_{\mathcal{T},2k}(t)$ (55). We show that

[TABLE]

where $\exp_{\mathcal{T}}$ denotes the time-ordered matrix exponential. Similar to the time-independent case, we find that the total gate complexity is O\big{(}(nt)^{1+\frac{1}{2k}}\big{)}. See (51, Section V) for detailed discussions.

We also consider lattice Hamiltonians with periodic boundary conditions $H=\sum_{j=1}^{n-1}H_{j,j+1}+H_{1,n}$, where $H_{j,k}$ represents a local term that acts nontrivially only on qubits $j$ and $k$. To simulate such a system, we decompose $H$ as $H=H_{\text{bndry}}+H_{\text{even}}+H_{\text{odd}}$, where $H_{\text{bndry}}=H_{1,n}$. Correspondingly, we also use a canonical product formula with three exponentials per stage. With a similar analysis for the open boundary condition, we find that the product-formula error is $O(nt^{2k+1})$ as expected.

A generalization of this approach allows us to simulate a $D$-dimensional lattice Hamiltonian with nearly optimal gate complexity. We use a $2D$-coloring of the edges of the lattice to decompose the Hamiltonian into $2D$ terms, each of which is a sum of commuting terms. We also extend the definition of canonical product formulas to allow for $2D$ exponentials per stage. An analysis of the local error structure shows that this algorithm has gate complexity $O((L^{D}t)^{1+\frac{1}{2k}}/\epsilon^{\frac{1}{2k}})=O((nt)^{1+\frac{1}{2k}}/\epsilon^{\frac{1}{2k}})$, where $n$ is the total number of lattice sites and $L=n^{\frac{1}{D}}$ is the linear size of the lattice.

Finally, our algorithm can also simulate lattice Hamiltonians with constant-range interactions $H=\sum_{j=1}^{n-\ell+1}H_{j,\ldots,j+\ell-1}$. To achieve nearly-optimal gate complexity, we classify the Hamiltonian terms into the $\ell$ groups

[TABLE]

and use a product formula with $\ell$ elementary exponentials per stage.

Discussion. The product-formula algorithm is arguably the simplest approach to quantum simulation. We have showed that this approach can simulate lattice Hamiltonians with nearly optimal gate complexity. Our algorithm invokes product formulas by ordering terms in an even-odd pattern, which is conceptually easy to understand and straightforward to implement. Beyond the one-dimensional time-independent system with nearest-neighbor interactions and open boundary conditions, our analysis is also applicable to periodic boundary conditions, constant-range interactions, higher dimensions, and the time-dependent case, all with nearly optimal gate complexity. Our result also gives product-formula bounds that are much tighter in practice.

Recently, Haah, Hastings, Kothari and Low (HHKL) proposed another nearly optimal algorithm for lattice simulation (12). Instead of analyzing the product-formula approach, they develop a new approach motivated by the Lieb-Robinson bound (56, 57), which quantifies how fast information can propagate in a system with local interactions. HHKL decomposes the entire evolution into blocks, where each block involves forward and backward evolution on a small region. Using product formulas within each block, their approach gives an ancilla-free algorithm for lattice simulation with asymptotic gate complexity $(nt)^{1+o(1)}$. However, this results in a much larger constant prefactor in practice than the pure product-formula algorithm analyzed here (51, Section VI).

The near optimality of HHKL depends essentially on the use of a Lieb-Robinson bound. As noted in (12), it may be difficult to apply this idea to Hamiltonians whose interactions are described by general graphs. Our approach directly exploits the commutation of lattice terms without the help of Lieb-Robinson bounds, which we expect could illuminate the simulation of other physical systems (9, 41, 42, 43).

Our local error analysis represents the product-formula error as an integral of a linear combination of commutators nested with unitary conjugations. Similar techniques have been used to establish the Lieb-Robinson bound and to study computational complexity aspects of many-body physics (56, 57, 41, 58). We leave it as an avenue for future work to explore whether these techniques could find more applications in the study of locality in quantum systems.

Acknowledgments

Acknowledgements.

Y.S. thanks Nathan Wiebe, Guang Hao Low, Minh Cong Tran, Stephen Jordan, Jeongwan Haah, Robin Kothari, Su-Kuan Chu, Rolando Somma, Xin Wang, Alexey Gorshkov, James R. Garrison, Anurag Anshu, Xiaodi Wu, Leonard Wossnig, Salini Karuvade, Xiao Yuan, Simon Benjamin, Shuchen Zhu, Scott Lawrence, and Zohreh Davoudi for helpful discussions. We thank anonymous referees for their comments. This work was supported in part by the Army Research Office (MURI award W911NF-16-1-0349), the Canadian Institute for Advanced Research, the National Science Foundation (grants 1526380 and 1813814), and the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Quantum Algorithms Teams and Quantum Testbed Pathfinder programs.