Toward the first quantum simulation with quantum speedup

To produce concrete benchmarks, we focus on a specific simulation task. Specifically, we consider a one-dimensional nearest-neighbor Heisenberg model with a random magnetic field in the $z$ direction. This model is described by the Hamiltonianwhere ${\overrightarrow{\sigma }}_j=\left({\sigma }_j^x,{\sigma }_j^y,{\sigma }_j^z\right)$ denotes a vector of Pauli $x$, $y$, and $z$ matrices on qubit $j$. We impose periodic boundary conditions (i.e., ${\overrightarrow{\sigma }}_{n+1}={\overrightarrow{\sigma }}_1$), and $h_j\in \left[-h,h\right]$ is chosen uniformly at random. The parameter $h$ characterizes the strength of the disorder.

This Hamiltonian has been considered in recent studies of self-thermalization and many-body localization (see SI Appendix, section B for more detail). Despite intensive investigation, the details of a transition between thermal and localized phases remain poorly understood. A major challenge is the difficulty of simulating quantum systems with classical computers; indeed, the most extensive numerical study we are aware of was restricted to, at most, 22 spins (18).

Hamiltonian simulation can efficiently access any feature that could be observed experimentally (and more), and there are several proposals for exploring self-thermalization by simulating dynamics (19–21). Since all of these approaches involve only very simple state preparations and measurements, we focus on the cost of simulating dynamics. We consider evolution times comparable to the number of spins, since the system must evolve for this long for self-thermalization to take place (or even for information to propagate across the system, owing to the Lieb–Robinson bound).

Specifically, we produce gate counts for simulations with $h=1$, evolution time $t=n$ (the number of spins in the system), and overall accuracy $ϵ=10^{-3}$. These explicit choices help us to focus on the system-size dependence of quantum simulation algorithms. This is a key consideration for practical applications, yet it has been deemphasized in the literature on sparse Hamiltonian simulation.

There are many distinct quantum algorithms for Hamiltonian simulation, some of which are summarized in Table 1. We implement algorithms based on high-order product formulas (PFs) (13), direct application of the Taylor series (TS) (15), and the recent quantum signal processing method (QSP) (16), all of which are introduced in SI Appendix, section C. We expect these to be among the most efficient approaches to digital quantum simulation. In particular, approaches based on quantum walk (14, 23) appear to incur greater overhead (as discussed in SI Appendix, section D).

To produce concrete circuits, we implement quantum simulation algorithms in a quantum circuit description language called Quipper (24) (see SI Appendix, section E for more details). Wherever possible, we tighten the analysis of algorithm parameters and manually optimize the implementation. We also process all circuits using an automated tool we developed for large-scale quantum circuit optimization (25). Our implementation is available in a public repository (https://github.com/njross/simcount).

We express our circuits over the set of two-qubit cnot gates, single-qubit Clifford gates, and single-qubit $z$ rotations $R_z\left(\theta \right)≔\exp \left(-i{\sigma }^z\theta /2\right)$ for $\theta \in \mathbb{R}$. Such gates can be directly implemented at the physical level with both trapped ions (2) and superconducting circuits (1, 2). In both technologies, two-qubit gates take longer to perform and incur more error than single-qubit gates. Thus, the cnot count is a useful figure of merit for assessing the cost of physical-level circuits on a universal device. We also produce Clifford+$T$ circuits using optimal circuit synthesis (26) so that we can count $T$ gates, which are typically the most expensive gates for fault-tolerant computation.

Our analysis ignores many practical details such as architectural constraints, instead aiming to give a broad overview of potential implementation costs that can be refined for specific systems. When counting qubits, we assume that measured ancillas can be reused later.

Fig. 1 compares gate counts for the PF algorithm (with commutator and empirical error bounds), the TS algorithm, and the QSP algorithm (in both segmented and nonsegmented versions). The TS algorithm uses significantly more qubits than the QSP algorithm (as shown in Fig. 2) while also requiring more gates, so the latter is clearly preferred. In contrast, the QSP algorithm has only slightly greater space requirements than the PF algorithm.

Surprisingly, despite being more involved, the QSP algorithm outperforms the rigorously bounded PF algorithm even for small system sizes. In particular, among the rigorously analyzed algorithms, the segmented QSP algorithm has the best performance, improving over the PF algorithm by about an order of magnitude for cnot count and by almost two orders of magnitude for $T$ count.

Empirical error bounds improve the performance of the PF algorithm by two to three orders of magnitude, making it the preferred approach if rigorous performance guarantees are not required. For the cnot count, the empirical PF algorithm improves over the full QSP algorithm by about an order of magnitude. The advantage in the $T$ count is less significant, but still indicates that the PF algorithm is dominant, especially considering its lower qubit count.

Although we expected that higher-order product formulas would not be advantageous for small system sizes, we find that the fourth- and sixth-order formulas had the best performance for our benchmark system with tens to hundreds of qubits, as shown in Fig. 3. The fourth-order formula with commutator bound gives the best available PF algorithm with a rigorous performance guarantee. Using empirical error bounds, the sixth-order formula outperforms the fourth-order formula for systems of about 30 or more qubits, making the former the method of choice for simulations just beyond the reach of classical computers (again, provided a heuristic error bound can be tolerated). These results suggest that higher-order formulas may be advantageous for other quantum simulations, such as those for quantum chemistry, even though they have not usually been considered (30, 32).

For a system of 50 qubits—which is presumably close to the limits of direct classical simulation for circuits such as ours (33)—the TS algorithm uses 171 qubits and the QSP algorithm uses 67, whereas the PF algorithm uses only 50. [Recent work has demonstrated simulation of 56-qubit computations, but only for circuits of much smaller depth than those considered in our work (34).] At this size, the segmented QSP algorithm is the best rigorously analyzed approach, using about $1.8\times 10^8$ cnot gates (over the set of Clifford+$R_z$ gates) and $2.4\times 10^9\hspace{0.17em}T$ gates (over the set of Clifford+$T$ gates). Using the empirical error bound, the PF algorithm uses about $3\times 10^6$ cnots and $1.8\times 10^8\hspace{0.17em}T$s (over Clifford+$R_z$ and Clifford+$T$, respectively).

For comparison, previous estimates of gate counts for factoring, discrete logarithms, and quantum chemistry simulations are significantly larger, as illustrated in Fig. 4. First, consider factoring a 1,024-bit number, which is beyond the factorization of RSA-768 that was achieved classically in 2009 (36). The best implementation we are aware of uses 3,132 qubits and about $5.7\times 10^9\hspace{0.17em}T$ gates (when realized over the set of Clifford+$T$ gates) (35). (Ref. 35 does not give explicit resource counts; we estimate them as described in SI Appendix, section A.) Quantum algorithms for classically hard instances of the elliptic curve discrete logarithm problem have roughly comparable cost (37). For quantum chemistry, a natural target for a problem just beyond the reach of classical computing is a simulation of FeMoco, the primary cofactor of nitrogenase, an enzyme that catalyzes the nitrogen fixation process. Even for a fairly low-precision simulation, and using nonrigorous estimates of the product formula error, the best implementation we are aware of uses 111 qubits and $1.0\times 10^{14}\hspace{0.17em}T$ gates (30). Thus, it appears that simulation of spin systems is indeed a significantly easier task for near-term practical quantum computation.

For a more detailed discussion of the results, see SI Appendix, section I.

The work described in this paper represents progress toward a possible early application of quantum computers, solving a practical problem that is beyond the reach of classical computation. Our optimized implementations can perform a 50-spin quantum simulation using fewer qubits, and dramatically fewer gates (by over five orders of magnitude in the case of quantum chemistry), than comparable instances of other classically hard problems. Of course, our results only represent upper bounds. While we attempted to optimize the implementation wherever possible, it is likely that further improvements can be found, and it is conceivable that another algorithm (or computational task) may offer better performance. Our work establishes a concrete set of benchmarks that we hope can be improved through future studies.

Demonstrations of digital quantum simulation performed to date (38–41) have been limited in scope, primarily using the first-order formula [except for some limited applications of the second-order formula (38, 39)]. Our results show that higher-order formulas are useful even for simulations of small systems. In the near term, it could be fruitful to demonstrate the utility of these formulas experimentally. Even relatively small experiments might be able to probe the validity of our empirical error bounds.

We have also identified some avenues for future improvement of quantum simulation algorithms. We saw that rigorous error bounds for product formulas are very loose, even with our newly developed commutator bound. This motivates attempting to prove stronger rigorous error bounds for product formulas. Also, the difficulty of computing the angles needed to perform the QSP algorithm prevents us from taking full advantage of the algorithm in practice, so it would be useful to develop a more efficient classical procedure for specifying these angles.

After the initial version of this paper was posted, Haah et al. (42) developed a new quantum algorithm for simulating the dynamics of spatially local Hamiltonians. Their approach uses bounds on the propagation of information in a local system to give a simulation with complexity only linear in the product of the system size and the evolution time (up to log factors), asymptotically improving over all previous methods. While the gate count estimates presented in figure 3 of ref. 42 suggest that, for the example considered in this paper, the PF and QSP methods remain favorable for systems with up to about 100 spins, the techniques developed in ref. 42 are likely to be useful in relatively near-term simulations.

Further reduction of the gate count could be especially significant if it led to a simulation with sufficiently few gates to be performed without invoking fault tolerance. With our current estimate of millions of cnot gates for a superclassical simulation, this is likely out of reach at present. However, further improvement could obviate the need for error correction in a system with highly accurate gates, making an early demonstration of superclassical simulation more accessible.

Finally, our work has considered an idealized system. Practical devices will come with architectural constraints, may use different basic operations than those considered here, may allow parallelization of gates, and will likely require fault tolerance. By incorporating such features, we hope the work begun here will lead to a blueprint for an early practical quantum computation.

We thank Zhexuan Gong, Alexey Gorshkov, Guang Hao Low, Chris Monroe, and Nathan Wiebe for helpful discussions. This work was supported, in part, by the Army Research Office (Multidisciplinary University Research Initiative Award W911NF-16-1-0349), the Canadian Institute for Advanced Research, and the National Science Foundation (Grant 1526380). This material was partially based on work supported by the National Science Foundation during D.M.’s assignment at the foundation. Any opinion, finding, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.