Faster quantum simulation by randomization

Product formulas can be used to simulate Hamiltonian dynamics on a quantum computer by approximating the exponential of a sum of operators by a product of exponentials of the individual summands. This approach is both straightforward and surprisingly efficient. We show that by simply randomizing how the summands are ordered, one can prove stronger bounds on the quality of approximation and thereby give more efficient simulations. Indeed, we show that these bounds can be asymptotically better than previous bounds that exploit commutation between the summands, despite using much less information about the structure of the Hamiltonian. Numerical evidence suggests that our randomized algorithm may be advantageous even for near-term quantum simulation.


Introduction
Simulating quantum dynamics is one of the major potential applications of quantum computers. The apparent intractability of simulating quantum dynamics with a classical computer led Feynman [16] and others to propose the idea of quantum computation. Lloyd gave the first explicit quantum algorithm for simulating the dynamics of local Hamiltonians [22], and later work showed that the more general class of sparse Hamiltonians can also be simulated efficiently [1]. Quantum simulation can be applied to understand the behavior of various physical systems-including quantum chemistry [2,25,32], quantum field theory [20], and many-body physics [27]-and designing new quantum algorithms [9,10,13,15,18].
The main ingredient in Lloyd's algorithm is the Lie product formula, which provides a firstorder approximation to the exponential of a sum as a product of exponentials of the summands. Given Hermitian operators H 1 , . . . , H L (which we refer to as the summands of the Hamiltonian H = L j=1 H j ) and a complex number λ, the Lie product formula approximates the exponentiation in the sense that V (λ) = lim r→∞ S 1 (λ/r) r . Suzuki systematically extended this formula to give a (2k)th-order approximation S 2k , defined recursively by with p k := 1/(4 − 4 1/(2k−1) ) [29]. Again we have V (λ) = lim r→∞ S 2k (λ/r) r , and the approximation obtained with a finite value of r improves as k increases (albeit with a prefactor that grows exponentially in k). We refer to all such formulas as product formulas. Berry et al. gave concrete error bounds for applying the Suzuki product formulas in quantum simulation [4]. Although other simulation algorithms have better asymptotic performance as a function of various parameters [5-8, 17, 23, 24], the product formula algorithm (especially with the first-order formula) is widely used in experimental implementations [3,11,21], due to its simplicity and the fact that it does not require any ancilla qubits. The main challenge in applying product formulas to quantum simulation is to choose the number of segments r to ensure the simulation error is at most some allowed threshold ǫ. To simulate H = L j=1 H j for time t, rigorous error analysis shows that suffices to ensure error at most ǫ for the first-order formula and suffices for (2k)th order [4], where Λ := max j H j . However, numerical simulations suggest that the product formula algorithm can perform significantly better in practice than the best proven error bounds demonstrate [2,14,27,28]. Indeed, recent work suggests that it can even asymptotically outperform more sophisticated simulation algorithms with better proven running times. This dramatic gap between the provable and the actual behavior of product formula simulation suggests that it may be possible to significantly improve their analysis. It is sometimes possible to improve the analysis of product formulas using further information about the form of the Hamiltonian. In particular, the cost of simulation can be reduced when many pairs of summands commute [2,14,22]. However, this approach can only be applied for structured Hamiltonians that contain many commuting summands. Furthermore, the best known bounds of this type give only modest improvement, remaining orders of magnitude away from the empirical performance even in cases where many summands commute [14].
In this paper, we propose a new approach to quantum simulation based on randomization. Previously, Poulin et al. used randomness to give improved simulations of time-dependent Hamiltonians by sampling the Hamiltonian at randomly chosen times [26]. Here we apply randomness in a different way: we randomize the ordering of the operators in the Hamiltonian. Specifically, for any permutation σ ∈ Sym(L) of the L summands, let We show that the (2k)th-order randomized simulation has error where V(−it) and S σ 2k (−it/r) are quantum channels describing the unitary transformation V (−it) and the random unitary S σ 2k (−it/r), respectively, and · ⋄ is the diamond norm (defined in Section 2). We prove this result using a mixing lemma due to Campbell and Hastings [12,19]. Informally, this lemma states that if we can approximate a desired operation as the average over some set of operations, then the overall error depends linearly on the error in the average operation but only quadratically on the error in any individual operation. Standard error bounds for product formulas do not depend on how the summands are ordered, but we show that randomizing the ordering gives a more accurate average evolution. We motivate this approach in Section 2, where we consider the effect of randomizing how the summands are ordered in the simple case of the first-order formula. Assuming Λ := max j H j is constant, our randomized first-order algorithm has gate complexity g rand 1 = O t 1.5 L 2.5 /ǫ 0.5 , improving over g det 1 = O t 2 L 3 /ǫ in the deterministc case. Analyzing the effect of randomization on higher-order formulas is more challenging. For terms of order at most L in the Taylor expansion of a product formula, the majority of the error comes from terms in which no summands are repeated. We call such contributions nondegenerate terms. In Section 3, we give a combinatorial argument to compute nondegenerate terms of the average evolution 1 L! σ∈Sym(L) S σ 2k (λ) in closed form. (In fact, we prove a more general result that applies to the average evolution as a special case.) As a corollary, we show that the nondegenerate terms completely cancel in the randomized product formula. Section 4 presents our main technical result, an upper bound on the error in a randomized higher-order product formula simulation. This bound follows by using the mixing lemma to combine an error bound for the average evolution operator with standard product formula error bounds for the error of the individual terms. Section 5 discusses the overall performance of the resulting algorithm and compares it with deterministic approaches. For the (2k)th-order product formula, assuming Λ := max j H j is constant, our randomized Hamiltonian simulation algorithm has complexity in the deterministic case. Thus our algorithm always improves the dependence on L and sometimes achieves better dependence on t and ǫ as well.
We also show in Section 5 that our bound can outperform a previous bound that takes advantage of the structure of the Hamiltonian. Specifically, we compare our randomized product formula algorithm with the deterministic algorithm using the commutator bound of [14] for a one-dimensional Heisenberg model in a random magnetic field. We find that over a significant range of parameters, the randomized algorithm has better proven performance, despite using less information about the form of the Hamiltonian.
In light of the large gap between proven and empirical performance of product formulas, it is natural to ask whether randomized product formulas still offer an improvement under the best possible error bounds. To address this question, we present numerical comparisons of the deterministic and randomized product formulas in Section 6. In particular, we show that the randomized approach can sometimes outperform the deterministic approach even with respect to their empirical performance.
Finally, we conclude in Section 7 with a brief discussion of the results and some open questions.

The power of randomization
To see how randomness can improve the product formula algorithm, consider a simple Hamiltonian expressed as a sum of two operators, H = H 1 + H 2 . The Taylor expansion of the first-order formula as a function of λ ∈ C is whereas the Taylor series of the ideal exponentiation is Using the triangle inequality, we can bound the spectral-norm error as where Λ := max{ H 1 , H 2 }. Since H 1 and H 2 need not commute, S 1 (λ) approximates V (λ) to first order in λ, as expected. It is clearly impossible to approximate V (λ) to second order using a product of only two exponentials of H 1 and H 2 : any such product can have only one of the products H 1 H 2 and H 2 H 1 in its Taylor expansion, whereas V (λ) contains both of these products in its second-order term. However, we can obtain both products by taking a uniform mixture of S 1 (λ) and Indeed, a simple calculation shows that However, S 1 (−it) + S rev 1 (−it) /2 is not a unitary operation in general. We could in principle implement it using LCU methods [6], but such an approach could have high cost, especially when the Hamiltonian contains many summands. A simpler approach is to apply one of the two operations S 1 (−it) and S rev 1 (−it) chosen uniformly at random, thereby implementing a quantum channel that gives a good approximation to the desired evolution. We now introduce some notation that is useful to analyze the performance of the randomized approach. Let X be a matrix acting on a finite-dimensional Hilbert space H. We write X for its spectral norm (the largest singular value) and X 1 for its trace norm (the sum of its singular values, i.e., its Schatten 1-norm). Let E : X → E(X) be a linear map on the space of matrices on H. The diamond norm of E is where the maximization is taken over all matrices Y on H ⊗ H satisfying Y 1 ≤ 1.
The following mixing lemma bounds how well we can approximate a unitary operation using a random unitary channel. Specifically, the error is linear in the distance between the target unitary and the average of the random unitaries, and only quadratic in the distance between the target unitary and each individual random unitary.
Lemma 1 (Mixing lemma [12,19]). Let V and {U j } be unitary matrices, with associated quantum channels V : ρ → V ρV † and U j : ρ → U j ρU † j , and let {p j } be a collection of positive numbers satisfying j p j = 1. Suppose that Then the average evolution E : To simulate the Hamiltonian H = H 1 + H 2 for time t, we divide the entire evolution into r segments of duration t/r and implement each segment via the random unitary operation (using one bit of randomness per segment). Invoking the mixing lemma with a = O (Λt) 2 /r 2 and b = O (Λt) 3 /r 3 , we find that Since the diamond norm distance between quantum channels accumulates additively under composition [31, p. 178], the error of the entire simulation is Thus the randomized first-order formula is effectively a second-order formula. This approach can easily be extended to the case where the Hamiltonian is a sum of L operators, again effectively making the first-order formula accurate to second order. Keeping track of all the prefactors, we get the following error bound for the randomized first-order formula.
To guarantee that the simulation error is at most ǫ, we upper bound the right-hand side of (20) by ǫ and solve for r. Assuming Λ := max j H j is constant, we find that it suffices to choose r rand 1 = O (tL) 1.5 /ǫ 0.5 , giving a simulation algorithm with gate complexity g rand 1 = O t 1.5 L 2.5 /ǫ 0.5 . In comparison, the gate complexity in the deterministic case is g det 1 = O t 2 L 3 /ǫ . Therefore, the randomized first-order product formula algorithm improves over the deterministic algorithm with respect to all parameters of interest.
It is natural to ask whether a similar randomization strategy can improve higher-order product formulas (as defined in (3)). While it turns out that randomization does not improve the order of the formula, it does result in a significant reduction of the error, and in particular, lowers the dependence on the number of summands in the Hamiltonian. The more complicated structure of higher-order formulas makes this analysis more involved than in the first-order case (in particular, we randomly permute the L summands instead of simply choosing whether or not to reverse them, so we use Θ(L log L) bits of randomness per segment instead of only a single bit). As discussed at the end of Section 1, our proof is based on a randomization lemma (established in the next section) that evaluates the dominant contribution to the Taylor series of the randomized product formula in closed form.
In its Taylor expansion, we call the sum of the form with coefficients α m 1 ...ms ∈ C, the sth-order nondegenerate term. This term contributes Θ(L s ) to the sth-order error, whereas the remaining (degenerate) terms only contribute O(L s−1 ). The following lemma shows how to compute the sth-order nondegenerate term for an arbitrary average evolution.
As an immediate corollary, we compute the sth-order nondegenerate term of the average evolution operator 1 L! σ∈Sym(L) S σ 2k (λ).

Corollary 1.
Let {H j } L j=1 be Hermitian operators; let λ ∈ C, k, s ∈ N, and s ≤ L. Then the sth-order nondegenerate term of the average evolution 1 Proof. The fact that S σ 2k (λ) is at least first-order accurate implies that q 1 + · · · + q κ = 1 in (24).
Observe that the sth-order nondegenerate term of V (λ) = exp λ L j=1 H j is also given by (28). Therefore, the sth-order nondegenerate term completely cancels in

Error bounds
In this section we establish our main result, an upper bound on the error of a randomized product formula simulation. To apply the mixing lemma, we need to bound the error of the average evolution. We now present an error bound for an arbitrary fixed-order term in the Taylor expansion of the average evolution operator.
Lemma 3. Let {H j } L j=1 be Hermitian operators; let λ ∈ C and k, s ∈ N. Define the target evolution V (λ) as in (2), and define the permuted (2k)th-order formula S σ 2k (λ) as in (6). Then the sth-order error of the approximation is at most 0 0 ≤ s ≤ 2k, where Λ := max H j .
The proof of this error bound uses the following estimate of a fixed-order degenerate term in the average evolution operator.
be Hermitian operators with Λ := max j H j ; let q 1 , . . . , q κ ∈ R with max k |q k | ≤ 1; and let s ≤ L be a positive integer. Then the norm of the sth-order degenerate term of the ideal evolution operator V (λ) as in (2) is at most and the norm of the sth-order degenerate term of the average evolution operator as in (22) is at most (κΛ|λ|) s s! L s − L(L − 1) · · · (L − s + 1) .
Proof. The sth-order term of V (λ) is We use the following strategy to bound the norms of these terms: (i) bound the norm of a sum of terms by summing the norms of each term; (ii) bound the norm of a product of terms by multiplying the norms of each term; (iii) bound the norm of each summand by Λ; and (iv) replace λ by |λ|.
Applying this strategy, we find that the norm of the sth-order term is at most where the nondegenerate term contributes precisely Taking the difference gives the desired bound (32). According to Lemma 2, the sth-order nondegenerate term of the average evolution is Following the same strategy as for V (λ) and also upper bounding the norm of each q k by 1 as part of step (iv), we find that the norm of this term is at most (κΛ|λ|) s s! L(L − 1) · · · (L − s + 1).
It remains to find an upper bound for the entire sth-order term of the average evolution. To this end, we start with the average evolution (22) and apply the following strategy: (i ′ ) replace each summand of the Hamiltonian by Λ; (ii ′ ) replace each q k by 1 and each λ by |λ|; and (iii ′ ) expand all exponentials into their Taylor series and extract the sth-order term. In other words, we extract the sth-order term of to get (κLΛ|λ|) s s! .
Proof of Lemma 3. We first prove a stronger bound, namely that the sth-order error is at most The first and third cases in this expression are straightforward. The formula S σ 2k is exact for terms with order 0 ≤ s ≤ 2k (this is what it means for the formula to have order 2k), so the error is zero in this case. When s > L, the randomization lemma is not applicable and the error can be bounded as in [14,Proof of Proposition F.3].
To handle the remaining case 2k < s ≤ L, we apply Lemma 4 with κ = 2 · 5 k−1 . The norm of the sth-order degenerate terms can be upper bounded by According to Corollary 1, the sth-order nondegenerate term of (30) cancels, which proves (42) for 2k < s ≤ L.
To finish the proof, we need a unified error expression for order s > 2k. When 2k < s ≤ L, we have with #{·} denoting the size of a set and [L] := {1, . . . , L}, where the inequality follows from the union bound. Therefore, we have If s > L ∈ N, we have s(s − 1) ≥ (L + 1)L ≥ 2L and This completes the proof.
We also use the following standard tail bound on the exponential function [14, Lemma F.2].

Lemma 5.
For any x ∈ C and κ ∈ N, we have We now establish the main theorem, which upper bounds the error of a higher-order randomized product formula.
Proof. We first prove that To this end, note that the sth-order error of On the other hand, Lemma 3 implies that the sth-order error of V (λ) − 1 so again Lemma 5 gives Equation (51) now follows from Lemma 1 by setting To simulate the evolution for time t, we divide it into r segments. The error within each segment is obtained from (51) by setting λ = −it/r. Then subadditivity of the diamond norm distance gives which completes the proof.

Algorithm performance and comparisons
We now analyze the complexity of our randomized product formula algorithm. Assume that k ∈ N is fixed, Λ = O(1) is constant, and r > tL. By Theorem 2, the asymptotic error of the (2k)th-order randomized product formula is To guarantee that the simulation error is at most ǫ, we upper bound the right-hand side of (58) by ǫ and solve for r. We find that it suffices to use While this bound quantifies the simulation error in terms of the spectral-norm distance, it can easily be adapted to the diamond-norm distance using either Lemma 1 or [8, Lemma 7]. This translation introduces only constant-factor overhead, so we have Therefore, the number of segments that suffice to ensure error at most ǫ satisfies giving an algorithm with elementary gates. Comparing to (60), we see that the randomized product formula strictly improves the complexity as a function of L. Indeed, the (2k)th-order randomized approach either provides an improvement with respect to all parameters of interest over the (2k)th order deterministic approach (if the first term of (60) obtains the maximum), or has better dependence on the number of terms in the Hamiltonian than any deterministic formula (if the second term dominates). We can also compare our result to the commutator bound of [14], which depends on the specific structure of the Hamiltonian. For concreteness, we consider a one-dimensional nearest-neighbor Heisenberg model with a random magnetic field, as studied in [14]. Specifically, let with periodic boundary conditions (i.e., σ n+1 = σ 1 ), and h j ∈ [−h, h] chosen uniformly at random, where σ j = (σ x j , σ y j , σ z j ) denotes a vector of Pauli x, y, and z matrices on qubit j. The (2k)th-order deterministic formula with the commutator bound has error at most [14,Eq. (146)] where we have again used Lemma 1 (or [8,Lemma 7]) to relate the spectral-norm distance to the diamond-norm distance. To guarantee that the simulation error is at most ǫ, it suffices to choose elementary gates. Comparing to the corresponding bound (60) for randomized product formulas, we see that the only difference is that the exponent 1/(2k + 1) for the commutator bound becomes 1/(4k + 1) in the randomized case. Thus the randomized approach can provide a slightly faster algorithm despite using less information about the structure of the Hamiltonian. More specifically, the relationship between t and L determines whether the randomized approach offers an improvement. If t = Ω(L 2k ), then the second term of (68) achieves the maximum, and both approaches have asymptotic complexity O tL 2 t ǫ 1 2k . However, if t = o(L 2k ), then the randomized formula is advantageous.

Empirical performance
While randomization provides a useful theoretical handle for establishing better provable bounds, those bounds may still be far from tight. As described in Section 1, our original motivation for considering randomization was the observation that product formulas appear to perform dramatically better in practice than the best available proven bounds would suggest. To investigate the empirical behavior of product formulas, we numerically evaluate their performance for simulations of the Heisenberg model (65) with t = n and h = 1, targeting error ǫ = 10 −3 , as previously considered in [14]. We collect data for the first-, fourth-, and sixth-order formulas as the latter two orders have the best performance in practice for small n and the first-order formula offers a qualitatively better theoretical improvement.
For the deterministic formula, we order the operators of the Hamiltonian in the same way as [14], namely We compute the error in terms of the spectral-norm distance and convert it to the diamond-norm distance using Lemma 7 of [8] (i.e., we multiply by 2). To analyze the randomized formula, we would like to numerically evaluate the diamond-norm distances and While the diamond norm can be computed using a semidefinite program [30], direct computation is prohibitive as the channel contains (L!) r Kraus operators. Instead, we use Lemma 1 to estimate the error. We randomly choose the ordering of the summands in each of the r segments, exponentiate each individual operator, and construct a unitary operator by concatenating the exponentials according to the given product formula. We use this procedure to obtain a Monte Carlo estimate of the average error using only three samples, as the standard deviations are already negligibly small (about 10 −5 ). We then invoke Lemma 1 to bound the diamond-norm error in (71). To the extent that the bound of Lemma 1 is loose, we expect the empirical performance to be better in practice. Using five randomly generated instances for each value of n, we apply binary search to determine the smallest number of segments r that suffices to give error at most 10 −3 . Figure 1 shows the resulting data for the first-, fourth-, and sixth-order formulas, which are well-approximated by power laws. Fitting the data, we estimate that segments should suffice to give error at most 10 −3 . We thus observe that the empirical complexity of the randomized algorithm is still significantly better than the provable performance We see that the randomized bound offers significantly better empirical performance at first order, consistent with the observation that randomization improves the order of approximation in this case. The fourth-order formula slightly improves both the exponent and the constant factor. While this improvement is small, it is nevertheless notable since it involves only a minor change to the algorithm. At sixth order we see negligible improvement. Since the proven bounds give less improvement with each successive order, it is perhaps not surprising to see that the empirical performance shows similar behavior.
To illustrate the effect of using different formulas and different error bounds to simulate larger systems, Figure 2 compares the cost of simulating our model system for sizes up to n = 100 with deterministic and randomized formulas of orders 4 and 6, using both proven error bounds and the above empirical estimates. (We omit the first-order formula since it is not competitive even at such small sizes.) We give rigorous bounds for deterministic formulas using the minimized bound of [14], and for fourth order we also show the result of using the commutator bound. We see that randomization gives a significant improvement over the deterministic formula using the minimized bound, although the commutator bound outperforms the randomized bound at the system sizes shown here. For sufficiently large n, the randomized bound gives lower complexity, but this requires a fairly large n since the difference in exponents is small and the commutator bound achieves a favorable constant prefactor. Empirical estimates of the error improve the performance by several orders of magnitude, with randomization giving a small advantage for the fourth-order formula as indicated above. However, for systems of size larger than about n = 25, the sixth-order bound prevails, and in this case randomization no longer offers a significant advantage.  Figure 2: Comparison of the total number of elementary exponentials for product formula simulations of the Heisenberg model using deterministic and randomized product formulas of fourth and sixth order with both rigorous and empirical error bounds. Note that since the empirical performance of deterministic and randomized sixth-order product formulas is almost the same, the latter data points are obscured by the former.

Discussion
We have shown that randomization can be used to establish better performance for quantum simulation algorithms based on product formulas. By simply randomizing how the summands in the Hamiltonian are ordered, we introduce terms in the average evolution that could not appear in any deterministic product formula approximation of the same order, and thereby give a more efficient algorithm. Indeed, this approach can outperform the commutator bound even though that method uses more information about the structure of the Hamiltonian. A randomized product formula simulation algorithm is not much more complicated than the corresponding deterministic formula, using only O(L log L) bits of randomness per segment and no ancilla qubits. Furthermore, we showed that randomization can even offer improved empirical performance in some cases.
While randomization has allowed us to make some progress on the challenge of proving better bounds on the performance of product formulas, our strengthened bounds remain far from the apparent empirical performance. We expect that other ideas will be required to establish tight bounds. Although our bounds have better asymptotic n-dependence than the previous commutator bound, they only offer an improvement if the system is sufficiently large. It could be fruitful to establish bounds for randomized product formulas that take advantage of the structure of the Hamiltonian, perhaps offering better performance both asymptotically and for small system sizes. More generally, it may also be of interest to investigate other scenarios in which random choices can be used to improve the analysis of quantum simulation and other quantum algorithms.