The Lie-Trotter product formula is commonly used in quantum algorithms for Hamiltonian simulation, where one can approximate the Hamiltonian evolution of H = A + B in terms of exponentials of A and B when these operators do not commute. For short time, a standard approximation is the second-order symmetric formula S2(t) = e−iAt/2e−iBte−iAt/2, which satisfies e−iHt = S2(t) + 𝒪(t3). More generally, the error in an order k formula is 𝒪(tk+1). Longer times are simulated by breaking the evolution into many repetitions of a short time. The number of repetitions needed is reduced with the order, motivating the search for higher-order product formulae. A systematic method to produce arbitrarily high-order formulae is the fractal method [1,2,3,4], which has found applications in Hamiltonian simulation [5]. The first explicit use of product formulae for quantum simulation was given in [6], applying it for quantum evolution under local Hamiltonians. Later work considered the broader class of sparse Hamiltonians [7] and higher orders [5], as well as methods beyond product formulae [8,9,10,11].
Recent work has shown that despite its simplicity, the Lie-Trotter product formula can compete with other Hamiltonian simulation algorithms due to the low error that it achieves in practice [12]. Methods based on linear combinations of unitaries [13,14] or quantum signal processing [15] give complexity logarithmic in the inverse error, but the error is not required to be extremely small, meaning those methods do not provide as large an advantage as might be expected. Product formula error bounds can be further improved by considering particular physical systems [16,17,18] or leveraging randomisation [19,20,21]. Moreover, Trotter formulae are expected to be relevant for both noisy intermediate-scale quantum (NISQ) computation and fault-tolerant computation. It is then of great importance to seek efficient implementations of product formulae as it can have a great impact on the efficiency of Hamiltonian simulation algorithms in practice.
The downside of the fractal methods to generate higher-order formulae is that the number of exponential operators to implement it grows very rapidly. Fractal product formulae are usually assumed in quantum computing, but they can be greatly improved upon. An alternative method by Yoshida [3] can be used to obtain product formulae with a smaller number of exponentials. Similar to fractal formulae, they are given as a product of S2 for different time intervals, but in contrast to the fractal approach there is not an explicit analytic form for the higher-order formulae. Instead one needs to derive and solve a complicated set of simultaneous non-linear polynomial equations.
In classical computing, product formulae are equivalent to symplectic integrators for numerical solution of differential equations. Since the initial work of Yoshida, there have been many results on better integrators of many different orders [22,23,24,25,26,27,28,29,30,31,32]. See Ref. [33] for a comprehensive review of product formulae/integrators in the literature. In particular, Refs. [24] and [25] give what appear to be the optimal 8th order product formulae with 15 and 17 stages (the number of S2 in the product formula). Reference [29] gives improved 8th order product formulae with 19 and 21 stages. For 10th order, Ref. [29] provides solutions with 31, 33, and 35 stages. These solutions from 2005 in Ref. [29] appear to be the best in prior work for 8th and 10th order (see Ref. [33]).
Another approach is that of processed product formulae [31,34,35,36,37,38,39,40]. Instead of the product formula being symmetric, it is of the form PΣP−1 for kernel Σ and processor P . The P and P−1 cancel when using the product of many of these for evolution over a long time, so the cost is dominated by the number of stages in the kernel. Since the kernel has fewer restrictions on it, it can have fewer stages than a normal symmetric product formula. According to Ref. [33], the best 8th order product formula with processing is from 2006 [31].
The objective of this work is to find improved product formulae and compare their performance in a consistent manner. We show that the performance of product formulae is better quantified by the error in the eigenvalues, rather than the spectral-norm error as usually considered in prior work. This is because it is the eigenvalue error that dominates the error for evolution at longer times. We also derive a method to consistently compare the performance of product formulae with different orders and numbers of exponentials.
We use these methods to compare between the performance of our product formulae, as well as to compare prior product formulae in the literature. By our numerical search, we found hundreds of thousands of solutions at 8th order, including many that outperform those previously reported in the literature. Our best 8th order product formula without processing improves over the best prior work by a factor of about 100, and our best 8th order processed product formula improves over the best prior processed product formulae by a factor of nearly 300. Our best product formula without processing has slightly lower error, but our product formula with processing provides better performance due to being shorter. For 10th order, one solution from Ref. [29] provides extremely high accuracy, a factor of nearly 300 times better than any others in prior work. We provide a detailed numerical comparison of product formulae in Section 5.
When comparing product formulae of different orders, it is better to use higher-order formulae for larger values of T /ε, where T is the total evolution time and ε is the required precision. For smaller T /ε it is better to use lower-order formulae, but as it is increased there are threshold values beyond which it is optimal to use higher-order formulae. We derive methods for determining these thresholds and show that our 8th order product formula is best for T /ε from around 107 to 1016, which includes the range of typical values for quantum chemistry applications. This means that the best 8th order product formula we have found in this work will be best suited to real applications.
The upper threshold means T /ε ≳ 1016 would be needed for 10th order product formulae to outperform our 8th order product formula, which is unrealistically large. Although the 10th order solution from Ref. [29] provides high accuracy, it is not competitive for quantum simulation due to its greater length. It is more than twice as long as our processed 8th order formula, so for a similar number of exponentials the 8th order product formula can use a time step half as long, resulting in much better performance for realistic values of the parameters. The lower threshold means T /ε ≲ 107 is needed for 6th order to outperform 8th order. Accurate estimation of this threshold requires testing of the product formulae for larger time steps. It is possible to reduce this threshold by adjusting the 8th order product formula to provide better performance for larger time steps.
The organisation of this work is as follows. In Section 2 we give a pedagogical explanation of the product formulae and the method of solution. First, we define product formulae and introduce the fractal methods for generating higher-order product formulae in Section 2.1, then in Section 2.2 we give a summary of Yoshida's method, then we explain processors for product formulae in Section 2.3, more general types of product formulae in Section 2.4, and the Taylor series method in Section 2.5. In Sections 3.2 and 3.3 we present the results regarding new product formulae. In Section 4 we explain the method for comparing product formulae, and then in Section 5 we give the comparison of the product formulae based on numerical testing. We finish in Section 6 by discussing what the results mean for choosing appropriate product formulae for quantum simulation, as well as searching for better-performing product formulae.
In this section, we review the background and methods for product formulae. Readers who are familiar with the methods may wish to skip to Section 3 where we present new results. See also Ref. [33] for a more complete recent review of this topic. We begin by defining product formulae and the Baker-Campbell-Hausdorff formula, then we introduce fractal methods and Yoshida's method to obtain higher-order formulae, and describe the processed product formulae and the Taylor series method.
It is well known that, for any non-commuting operators X and Y ,
Let X and Y be any operators such that ‖X‖ + ‖Y ‖ < ln 2. We have for an operator Z that exp(X) exp(Y) = exp(Z), with
The standard second-order symmetric product formula is as given in the definition below.
Let X and Y be non-commuting operators and let t be a real variable. Define
The operators X and Y used in the definition of S2 should always be clear from context. More generally, when there is a sum of Xj, the product formula is
Let X and Y be any operators such that ‖X‖ + ‖Y ‖ < ln 2 and let t be a real variable. Define Z such that S2(t) = exp(Z). Then there is a sequence αℓ consisting of linear combinations of ℓ-term commutators in X and Y such that
Reference [3] also shows that even terms are zero for more general symmetric product formulae. The first three non-zero αℓ terms from above are
Let X and Y be any non-commuting operators. Given a natural number n, a product formula of order n is a pair
Hence S2 is a length-3 product formula. The number of exponentials used in the product formula is a crucial measure of its efficiency, and in quantum simulation it is proportional to the required number of gates.
Here we describe two fractal methods to obtain higher-order product formulae. Starting from the symmetrised product formula in Eq. (3), fractal methods produce product formulae at all even orders. The first fractal (or triple jump) technique to generate a product formula of order k = 2κ is [1,2,3,42,43]
The history of this approach is discussed in detail in a review article by Forest [44]. This form for a fourth-order integrator was originally derived by Ruth about 1984, but not published until 1990 [42]. That work specifically considered a Hamiltonian that is a sum of two parts. The iterative form was given by Creutz and Gocksch in 1989 [1], but citing unpublished work by Campostrini and Rossi for the fourth-order case. That work was published in 1990 [43], and cited Creutz and Gocksch for the iterative form. Suzuki [2] and Yoshida [3] independently derived the iterative form in 1990.
The iterative form given by Creutz and Gocksch in 1989 [1] was the more general expression
For the first fractal method, the total number of exponentials for a given order 2κ = 4, 6, 8, ... in the product formula S2κ is given by
Here, we give a pedagogical explanation of the method to derive product formulae.
Rather than using Eq. (11) and Eq. (13), Yoshida uses the general procedure [3]
The method is to expand Eq. (16) using the BCH formula from Theorem 1 recursively as follows. First, note that from Corollary 1,
We are for the moment considering sums of two terms X + Y . Define
The case m = 0 is just the symmetric BCH formula, so it is clear that Eq. (21) holds with
To derive a 6th order formula, the lower-order terms in the exponential in Eq. (21) must be zero (see also Eq. (5.16) of [3]), which gives the four conditions
The product formulae obtained through the Yoshida method also work for exponentials of sums of more operators. The S2 product formula can again be used as the building block for the product formula, and we can write
Another technique to obtain higher-order product formulae is that of processing [31,34,35,36,37,38,39]. In this technique a product formula Sk of order k is generated by the composition of a kernel Σ conjugated by a processor operator P as
Typically, one chooses a kernel Σ as a product formula that has a smaller order than k, but conjugating by processors gives an order k integrator. Some of the best-performing kernels are given in Ref. [31], and recently new kernels were reported in [40]. The advantage of this method is that the number of non-linear equations required to be solved to find new product formulae is reduced. One could use a shorter minimum-length product formula, but best performance is provided by increasing the length of a product formula beyond the minimum. The more the number of free parameters exceeds the number of equations, the more freedom there is to find product formulae with low error. Reducing the number of equations by solving for the kernel enables reduced error for a given length.
To be more specific, Ref. [31] gives a set of conditions for the kernels and processors in Table 4 of that work. Reference [31] uses the notation Yj rather than αj, so S2(t) = etY1+t3Y3⋯. It then uses the notation Ej,i for multicommutator expressions, with Ej,1 = Yj ≡ αj and for example
To obtain the conditions for a 6th order kernel from Ref. [31], one should use the conditions for fj,i up to q = 5 in Table 4 of that work, which give f1,1 = 1, f3,1 = 0, f5,1 = 0, which in the above notation are equivalent to
For 4th order there are the same number of conditions, so the kernel is a valid product formula of that order. For the processors, Ref. [31] uses the notation pj,i for the coefficients of Ej,i in the expansion of P . These need to be determined from the BCH formula rather than symmetric BCH, and instead of only odd-order terms there need to be only even-order terms. Some symmetries mean that low-order odd terms automatically cancel, but higher-order ones need to be made zero by the choice of P . Table 4 in Ref. [31] gives the conditions on the even-order terms. The first processor condition depending on fj,i is for q = 5, which means that for 6th order and higher the processors are dependent on the kernel. The 4th order kernel is already 4th order, and the processors will yield another 4th order product formula, but can be chosen to reduce the error (which arises from higher-order terms). See Appendix B for more detail on the method of determining processors.
The form of product formulae constructed from products of S2 can be generalised in a number of ways. First, one may choose different weightings in the two parts of Eq. (4), so one has
In the case where the Hamiltonian is a sum of only two parts, then a product formula may be written in the form
In the case where the Hamiltonian is a sum of two terms, the product of S(t1, t2) gives a form as in Eq. (34), and conversely the product formula as in Eq. (34) can be expressed as a product of S(t1, t2). In Section 5 we give results for testing some product formulae of this type, and find that they give the best results for 4th and 6th order. This method is mostly considered only up to 6th order because for higher orders the number of order conditions becomes significantly larger (see Table I of Ref. [31]).
Another class of product formulae is those constructed from products of fourth-order product formulae S4. Some product formulae have been constructed in that way; for example, see Ref. [31]. We have also tested those (using the fourth-order fractal product formulae), but do not include them in the comparison in Section 5. We found that, although they reduced the error, the much larger number of exponentials in the product meant that they were not competitive. Part of the reason is that the best gains are obtained when there are significantly more free parameters to adjust than the number of equations to satisfy. This is the reason to use product formulae with more than the minimum length. Constructing a product formula from S4 increases the length without the added advantage of more free parameters to adjust.
Another option is to construct product formulae from a different second-order product formula. For example, a more accurate second-order method is provided in Ref. [24]. In our testing, we find that yields slightly more accurate product formulae, but with twice the length, making them uncompetitive with product formulae constructed from the standard second-order symmetric product formula S2.
Another solution method is based on computing the Taylor expansion of both the exact exponential and its product formula approximation with given wj. This method is not normally used because it is much more computationally intensive than the symmetric BCH approach. Here, we describe how to effectively code this method for numerical solution. The Taylor expansion is performed on both sides up to the desired order of approximation for the integrator. By imposing equality on terms of the same order one can obtain a series of equations for wj which can be solved. For higher orders, a large number of simultaneous equations are obtained, so we need an automated way of generating them.
To make precise the problem we are solving, denote as 𝒯k[⋅] the map giving the Taylor expansion in t around 0, truncated at order k, so
The product formula obtained from the solution procedure needs to be independent of the choice of X and Y , so we need to match the coefficients for each sequence of products of X and Y . Because X and Y are non-commuting, we need to record coefficients for each ordered sequence. To do this in an automated way we construct a data structure to store the coefficients.
Given operators of the form
To construct a bit string, we map each X to 0 and each Y to 1, then place a 1 on the left to flag the length of the string, as shown in the second line of Eq. (39). Then, to obtain the operator products, simply remove the leading 1 and then map 0 to X and 1 to Y . The empty string corresponds to the identity I. Now, to store the coefficients in a sum of products of X and Y , convert each product to a binary integer as above, then store the coefficient in the corresponding location in a vector. This approach enables products of Taylor expansions of exponentials to be rapidly calculated.
By way of illustration, consider the polynomial 10I + 3X + 2Y + 2X2 + Y X. This operator would be stored in an array as [10, 3, 2, 2, 0, 1, 0, ⋯]. In this way, any polynomial of X and Y can be efficiently stored in a vector. We denote this vector storing the coefficients of operators of order up to k as veck[p(X, Y)], where k denotes that the vector will only store the coefficients of the corresponding operators up to order k (so a vector of length 2k+1 − 1) and p(X, Y) is the polynomial in terms of X and Y.
First we describe the numerical methods, then summarise solutions we found at 8th and 10th order.
The size of the product formula we consider as in Eq. (16) is governed by the value of m, so the number of S2 in the product is M = 2m + 1. The number of free parameters is m, and it needs to be a least a large as the number of independent equations in order for there to be a solution. Yoshida [3] finds that m = 7 is minimal for 8th order, and m = 15 is minimal for 10th order [29]. Choosing the minimal m typically yields product formulae with large error, and for that reason we focus on larger values of m.
For the numerical solution of the simultaneous non-linear equations, we found that Matlab's fsolve (using the Levenberg-Marquardt algorithm) provided rapid results using the vector of errors. The method was to repeatedly solve from random starting vectors
The Taylor series approach could be used to solve for product formulae, but is more computationally intensive because it produces an error vector that is exponentially large in the order. However, it is very effective at further improving the performance of solutions that have already been found. If a solution is of 8th order (for example), then we can consider the error using the Taylor expansion up to and including the 9th order. This will give non-zero error for the 8th order product formula, so minimising the error gives a new product formula. That can be used as a starting point to again solve for an 8th order formula, and often it is found that the product formula adjusted in that way has reduced error.
Following these methods, we have searched for new product formulae and found thousands of product formulae of 8th order. Once we have obtained the potential solutions, we generate random Hamiltonians and compute the product formula errors for a chosen short time. The error for 8th order product formulae is given by δ(t) ∝ t9, but the constant will depend on the measure of error used. One could use the spectral-norm error for a single time step, but a more appropriate measure of the performance of the product formulae is the accuracy of the eigenvalues. In the following we use χ for the constant factor based on spectral-norm error and ζ for the constant based on eigenvalue error, with χ = δ/t9 or ζ = δ/t9.
For each product formula, we compute an average constant factor error; this average corresponds to the geometric mean of the constant factors computed for each random Hamiltonian. See Section 5 for more details of the numerical methods, and the comprehensive comparison of product formulae. This method of comparing the performance of product formulae through the estimation of the constant factor in the error has been used before (see for example [32]) and is considered a good approximation to the performance of the product formula in practice.
The case m = 7 yields product formulae of minimal length because the number of equations is equal to the number of unknowns. We have found over 600 solutions, and the search now finds almost only repeated solutions and very few new solutions. This indicates that we have found nearly all solutions, but it is also possible that there are many more solutions with large values of
This search for m = 7 (15 stages) verified that the best-performing product formula corresponds to that previously reported in Refs. [24] and [25]. Similarly, a numerical search with m = 8 (17 stages) yielded the best-performing product formula close to solutions s17odr8a and s17odr8b given in Ref. [25]. It is not exactly the same, but the small difference is likely due to the system of equations being underdetermined. The extensive nature of our search indicates that these product formulae reported in Refs. [24,25] are optimal for 8th order product formulae with m = 7 and m = 8.
To provide improved performance over prior work, we extended the search to m = 10 (21 stages), and found over 100,000 solutions. This size of 8th order product formula was previously studied in Ref. [29], but we find solutions with significantly lower error. Our best solutions for spectral-norm and eigenvalue error are given in Table 1. The solution selected for minimum spectral-norm error has a constant of χ = 5.8 × 10−8. Its constant factor for eigenvalue error ζ = 7.0 × 10−9 is significantly smaller, but the best performance is obtained by selecting for minimum eigenvalue error. The second solution given in Table 1 has the significantly smaller constant factor in the error ζ = 5.4 × 10−10.
Our best-performing 8th order solutions when setting m = 10.
| Best 8th order for spectral-norm error | Best 8th order for eigenvalue error | |
|---|---|---|
| w1 | 0.59358060400850625863514059265224 | 0.10467636532245895252340732579853 |
| w2 | −0.46916012347004197296293264921328 | −0.57896999331780988041471955125778 |
| w3 | 0.2743566425898467907228242878146 | 0.57503350160061785946141563279891 |
| w4 | 0.17193879484656773059919074965377 | 0.12231011868707029786561397542663 |
| w5 | 0.23439874482541384415430578747541 | 0.27793149999039524816733903301747 |
| w6 | −0.48616424480326193899617759997914 | −0.37349605088056728482635987352576 |
| w7 | 0.49617367388114660354871757044906 | 0.11575566589480463220616543972403 |
| w8 | −0.32660218948439130114501815323814 | 0.1464645610975800618712569230326 |
| w9 | 0.23271679349369857679445410270557 | −0.39443578322284085764474498594073 |
| w10 | 0.098249557414708533273471906180643 | 0.44370228726021218923197141183196 |
These constant factors all improve over those provided by product formulae in prior work. Our solution with smallest eigenvalue error is well over 100 times more accurate than the m = 10 product formula in Ref. [29], and about 100 times more accurate than the m = 9 product formula in that work. See Section 5 for the detailed numerical testing of the product formulae.
A further improvement in product formulae can be obtained by using processing. We have performed a search for high-accuracy 8th order product formulae using this method with m = 8, and found one particular example that provides excellent performance; see Table 2. This kernel provides a constant of ζ = 8.1 × 10−10, only slightly larger than the best length m = 10 product formula, while providing better performance due to the shorter length. This error is nearly 300 times smaller than provided by 8th order processed product formulae in prior work [31]. Therefore, we have provided a sequence of three solutions that all improve over prior work, with successively better and better performance.
Our best-performing 8th order kernel and processor when setting m = 8. The value of γ10 for the processor is chosen such that the sum is zero.
| Best 8th order kernel | Processor for kernel | ||
|---|---|---|---|
| w1 | 0.21784176681731006074681969186513 | −0.44324901019570126590495430949294 | γ1 |
| w2 | 0.1947017706053903224022456342907 | 0.25459857192003772850622377066944 | γ2 |
| w3 | 0.18372413281145589944261642180363 | −0.73862036266779261573694538099739 | γ3 |
| w4 | −0.37307499512657736825709230652023 | −0.00024139614958652134370419495289618 | γ4 |
| w5 | 0.15757644257569146373033662060461 | 0.73873460354125365739379753874964 | γ5 |
| w6 | −0.33342207567391682979227850551172 | −0.20285971152536085519251666906017 | γ6 |
| w7 | 0.51788649682987924281787142226803 | 0.44989521689676869571827637424046 | γ7 |
| w8 | 0.21456475499897766986381219621761 | 0.29538398007876871184026747505657 | γ8 |
| −0.3364996155865700091428329802017 | γ9 | ||
To provide an accurate result in terms of the spectral norm (rather than just the eigenvalues) it is necessary to solve for the processor. For the kernel there are many solutions for the processor, but we list the one with the best performance of the solutions we have found. This one was found by solving for many (hundreds of thousands) of processors, picking the one with the best performance, and further optimising. Moreover, using the processor the constant for the spectral-norm error is χ = 5.3 × 10−8. That is even better than the best length m = 10 (nonprocessed) formula above.
We have also used our solution procedure to find new 10th order product formulae. For Yoshida's method applied to 10th order, there are 15 independent equations to be solved (see Appendix A). We performed searches for solutions with m = 15 (the minimal number) to m = 18. Again the larger values of m give the flexibility to adjust the solution to reduce the error. As in Section 3.2, we compare the performance of product formulae of 10th order by computing the constant factors χ and ζ for random Hamiltonians.
Out of our solutions, the best eigenvalue performance is given by a m = 17 solution (see Table 3), which yields ζ = 6.1 × 10−11. This formula is slightly outperformed by one product formula in Ref. [29], but for the purpose of quantum computing we will show that 8th order provides better performance than 10th order. For this reason we did not perform further searches for processed 10th order product formulae or proceed to any higher orders.
The 10th order solution with lowest eigenvalue error in our numerical search.
| 10th order solution with m = 17 | |
|---|---|
| w1 | −0.28371232689144296279654621726493 |
| w2 | 0.046779504778147381605331000278223 |
| w3 | 0.36845892382797770619657504217539 |
| w4 | 0.19186204094674514739760408197461 |
| w5 | −0.53123134392680669702873064192428 |
| w6 | −0.0081253242720827266680816105600661 |
| w7 | −0.16389450414378567860032917538393 |
| w8 | 0.18514766119291405032528647881 |
| w9 | 0.5383584694754681989174668806505 |
| w10 | −0.30583981835573485697292316732177 |
| w11 | 0.43199935609523301289295473774488 |
| w12 | 0.1510502301631786853020124612813 |
| w13 | −0.35051099204829676098801520498121 |
| w14 | 0.1032971125844291674511513007661 |
| w15 | 0.15043936943817152697371946806229 |
| w16 | 0.12118469498650736511410491586846 |
| w17 | 0.10437742779547826358296681557444 |
In the search for 10th order product formulae, unlike in the case of 8th order, we find that almost all new solutions found are different from those found before. That indicates there is an extremely large number of solutions, and we have only found a very small proportion of them. There is still the potential for better-performing solutions to be found in this order, but we expect that the longer length will make them unsuitable for quantum computing.
A fair comparison between product formulae of different lengths may be made in the following way. An order k product formula for time t will have an error δ = ζtk+1 (using ζ for eigenvalue error, or δ = χtk+1 for the case of spectral-norm error). Let T be the total evolution time for a product formula of order k, and ε be the maximum allowable error. Subdivide the evolution time T into r subintervals, so t = T /r is the length of each time subinterval. We thus have ζ(T /r)k+1 ≈ ε/r, which gives
As another motivation for this form, consider the case where the task is to estimate an eigenvalue of the Hamiltonian. Then one would use an ancilla register to control the time of evolution, and use phase estimation. If we are controlling between the forward and reverse time evolution as per Ref. [45], then the (root-mean square) error in estimating the phase of e−iHT /r is given by
If we wish to compare product formulae of different order, then we need to take account of the values of T and ε. Assume we have two product formulae of order k1 and k2, with corresponding constants ζ1, ζ2. Given T and ε, when the two product formulae use the same number of exponentials we have M1r1 = M2r2, thus
A limitation of this analysis is that it is assumed the time step (ε/(ζT))1/k is small enough for the scaling law for the error to hold. To adjust the threshold for cases where the time step is large, we can consider a more general functional dependence of the error on the time interval as f(t). Then we require for time T and error ε that
Here, we give further explanation of why it is important to consider eigenvalue error. If the goal is to estimate eigenvalues of the Hamiltonian (as is often the case in quantum chemistry), then the measure of the error should be that in the eigenvalues. That is, how close are the eigenvalues of the product formula to those of the exact Hamiltonian evolution. In the case that the goal is to accurately reproduce the final state after the evolution, then an appropriate measure of the error should be the spectral norm of the difference of the unitary operators. This error accounts for both error in the eigenvalues and basis.
That error will upper bound the 2-norm error in the generated quantum state, but using the triangle inequality to bound the error for long evolution times overestimates the error. This is because the error beyond that in the eigenvalues cancels when using a product of many short time intervals. Let us denote the exact evolution operator for short time U , and the approximate evolution operator provided by the product formula as
It is convenient to write
Then, for the case of a large number of time steps r, we can bound the overall error as
In practice we find that in many cases the error in the basis scales as 𝒪(tk), resulting in the expected contribution to the spectral-norm error 𝒪(tk+1). However, in other cases we find scaling of the error in the basis as 𝒪(tk+1), resulting in a higher-order contribution to the spectral-norm error. Let us give the error in the eigenvalues as ζtk+1, and the error in the basis as μtk+ν where ν ∈ [0, 1] to account for the various scalings found numerically.
Then the choice of r to make the error in the eigenvalues smaller than ε is
In what follows, we report on the comparison between product formulae of 4th, 6th, 8th, and 10th order based on the threshold provided above. We analyse those found in this work, and those from prior work including processed product formulae. When giving the number of stages of a processed product formula, we refer to the number of stages in the kernel (not P) because those are the ones that would be repeated for a simulation over an extended time.
We list 4th order product formulae in Table 4, 6th order in Table 5, 8th order in Table 6, and 10th order in Table 7, and give the constant factors in the error scaling χ, ζ (spectral norm or eigenvalue error). These product formulae include some constructed with even M rather than M = 2m + 1 as for our product formulae. We generated 1024 random Hamiltonians for each product formula, and then estimated the constants χ and ζ by taking the geometric mean. Each random Hamiltonian is generated as a sum of two random Hamiltonians H = A + B, where A and B are random Hermitian matrices of dimension 64 × 64 and norm 1.
The list of 4th order product formulae. The column labelled “processing” indicates whether the formula uses processor. The label is the name we will refer to the formula by. The constant factor in the error is denoted χ for spectral-norm error and ζ for eigenvalue error, and the corresponding quantities Mχ1/k and Mζ1/k are given. The column labelled S2 indicates whether the product formula is a product of S2. The best result is highlighted in green, and is from 2002 [47]. The 2006 result PPBCM4m6 is slightly higher error, but it is not evident to 2 significant figures.
| label | M | S2 | processing | reference | χ | Mχ1/k | ζ | Mζ1/k |
|---|---|---|---|---|---|---|---|---|
| S4m1 | 3 | Y | N | first fractal product [2] | 4.9 × 10−2 | 1.41 | 2.3 × 10−2 | 1.17 |
| S4m2 | 5 | Y | N | second fractal product [4] | 3.0 × 10−3 | 1.17 | 3.3 × 10−4 | 0.67 |
| O4M5 | 5 | N | N | Ostmeyer Eq. (40) of [46] | 3.0 × 10−4 | 0.66 | 7.7 × 10−5 | 0.47 |
| BM4M6 | 6 | N | N | S6 Table 2 of [47] | 1.6 × 10−4 | 0.67 | 2.4 × 10−5 | 0.42 |
| PPBCM4m6 | 6 | N | Y | P64 Table 5 of [31] | 5.7 × 10−5 | 0.52 | 2.4 × 10−5 | 0.42 |
| BCE4m3 | 6 | Y | Y |
| 1.6 × 10−2 | 2.15 | 1.5 × 10−3 | 1.17 |
| BCE4m4 | 8 | Y | Y |
| 3.4 × 10−4 | 1.09 | 9.8 × 10−5 | 0.80 |
| BCE4m5 | 10 | Y | Y |
| 6.7 × 10−5 | 0.90 | 2.1 × 10−5 | 0.68 |
| BCE4m6 | 12 | Y | Y |
| 2.6 × 10−5 | 0.85 | 7.0 × 10−6 | 0.62 |
| BCE4m7 | 14 | Y | Y |
| 1.3 × 10−5 | 0.85 | 3.0 × 10−6 | 0.58 |
| BCE4m8 | 16 | Y | Y |
| 7.9 × 10−5 | 0.85 | 1.5 × 10−6 | 0.56 |
| BCE4m9 | 18 | Y | Y |
| 5.3 × 10−6 | 0.86 | 8.6 × 10−7 | 0.55 |
The list of 6th order product formulae with their average errors. The column labelled S2 indicates whether it is a product of S2. The best result is highlighted in green, and is from 2006 [31].
| label | M | S2 | processing | reference | χ | Mχ1/k | ζ | Mζ1/k |
|---|---|---|---|---|---|---|---|---|
| S6m1 | 9 | Y | N | first fractal product [2] | 4.5 × 10−2 | 5.36 | 2.3 × 10−2 | 4.81 |
| S6m2 | 25 | Y | N | second fractal product [4] | 1.3 × 10−5 | 3.83 | 2.0 × 10−7 | 1.91 |
| Y6m3a | 7 | Y | N | Solution A in Table 1 of [3] | 2.0 × 10−3 | 2.49 | 1.2 × 10−3 | 2.28 |
| KL6s9a | 9 | Y | N | s9odr6a in Appendix A of [25] | 2.9 × 10−4 | 2.31 | 1.7 × 10−4 | 2.12 |
| KL6s9b | 9 | Y | N | s9odr6b in Appendix A of [25] | 2.9 × 10−4 | 2.31 | 1.7 × 10−4 | 2.11 |
| SS6s11 | 11 | Y | N | Section 4.2 of [29] | 4.2 × 10−5 | 2.05 | 1.4 × 10−5 | 1.71 |
| SS6s13 | 13 | Y | N | Section 4.2 of [29] | 2.1 × 10−5 | 2.16 | 3.4 × 10−6 | 1.59 |
| BM6M10 | 10 | N | N | S10 Table 2 of [47] | 6.0 × 10−6 | 1.35 | 1.4 × 10−6 | 1.06 |
| PPBCM6m9 | 9 | N | Y | P96 Table 5 of [31] | 5.0 × 10−6 | 1.18 | 4.3 × 10−7 | 0.78 |
| PPBCM6m5 | 11 | Y | Y | P116 in Table 6 of [31] | 1.9 × 10−6 | 1.22 | 9.0 × 10−7 | 1.08 |
| PPBCM6m6 | 13 | Y | Y | P136 in Table 6 of [31] | 4.8 × 10−7 | 1.15 | 2.5 × 10−7 | 1.03 |
| BCE6m5 | 10 | Y | Y |
| 5.8 × 10−3 | 4.24 | 3.2 × 10−3 | 3.84 |
| BCE6m6 | 12 | Y | Y |
| 7.7 × 10−5 | 2.47 | 2.0 × 10−5 | 1.98 |
| BCE6m7 | 14 | Y | Y |
| 1.6 × 10−5 | 2.22 | 2.4 × 10−6 | 1.62 |
| BCE6m8 | 16 | Y | Y |
| 4.6 × 10−6 | 2.06 | 3.2 × 10−7 | 1.32 |
| BCE6m9 | 18 | Y | Y |
| 7.3 × 10−6 | 2.51 | 1.3 × 10−7 | 1.28 |
| BCE6m10 | 20 | Y | Y |
| 3.2 × 10−6 | 2.42 | 9.9 × 10−9 | 0.93 |
| BCE6m11 | 22 | Y | Y |
| 7.9 × 10−6 | 3.11 | 1.4 × 10−8 | 1.08 |
The list of 8th order product formulae with their average errors. This table shows our new results for product formulae that improve over those in the prior literature. The result that provides the most efficient simulations is highlighted in (dark) green, and is our processed product formula. That provides the best performance due to the shorter length, but the lowest error is provided by our product formula highlighted in (dark) blue. The lowest error results from prior work are highlighted in light blue (for the non-processed case), and light green (for the processed case). The processed product formula PP8s13 (P138) has larger error than the highlighted PP8s19, but better performance due to its shorter length.
| label | M | processing | reference | χ | Mχ1/k | ζ | Mζ1/k | |
|---|---|---|---|---|---|---|---|---|
| S8m1 | 27 | N | first fractal product [2] | 5.6 × 10−2 | 18.8 | 1.6 × 10−2 | 16.2 | |
| S8m2 | 125 | N | second fractal product [4] | 6.7 × 10−9 | 11.9 | 5.2 × 10−13 | 3.64 | |
| Y8m7d | 15 | N | Solution D in Table 2 of [3] | 1.1 × 10−3 | 6.41 | 1.3 × 10−4 | 4.89 | |
| MC8s15 | 15 | N | Table 2 of [24] | 6.5 × 10−6 | 3.37 | 2.0 × 10−6 | 2.90 | |
| MC8s17 | 17 | N | Table 2 of [24] | 7.9 × 10−7 | 2.94 | 3.3 × 10−7 | 2.63 | |
| KL8s17a | 17 | N | s17odr8a in [25] | 6.1 × 10−7 | 2.84 | 2.7 × 10−7 | 2.56 | |
| KL8s17b | 17 | N | s17odr8b in [25] | 5.9 × 10−7 | 2.83 | 2.5 × 10−7 | 2.54 | |
| SS8s19 | 19 | N | Section 4.3 of [29] | 1.8 × 10−7 | 2.72 | 5.3 × 10−8 | 2.34 | |
| SS8s21 | 21 | N | Section 4.3 of [29] | 3.4 × 10−7 | 3.27 | 7.5 × 10−8 | 2.70 | |
| PP8s13 | 13 | Y | P138 in Table 6 of [31] | 1.2 × 10−6 | 2.37 | 6.5 × 10−7 | 2.19 | |
| PP8s19 | 19 | Y | P198 in Table 6 of [31] | NA | NA | 2.4 × 10−7 | 2.83 | |
| Y8m10 | 21 | N | Table 1 (our new result) | 5.8 × 10−8 | 2.61 | 7.0 × 10−9 | 2.01 | |
| Y8m10b | 21 | N | Table 1 (our new result) | 6.3 × 10−7 | 3.53 | 5.4 × 10−10 | 1.46 | |
| YP8m8 | 17 | Y | Table 2 (our new result) | 5.3 × 10−8 | 2.09 | 8.1 × 10−10 | 1.24 | |
The list of 10th order product formulae with their average errors. The best result is highlighted in green, and is the product formula from 2005 [29]. That product formula provides exceptional performance for 10th order, being about 280 times better than all others from prior work, and half the error of our solution (bottom line).
| label | M | processing | reference | χ | Mχ1/k | ζ | Mζ1/k | |
|---|---|---|---|---|---|---|---|---|
| S10m1 | 81 | N | first fractal product [2] | 9.0 × 10−2 | 63.7 | 2.7 × 10−3 | 44.8 | |
| S10m2 | 625 | N | second fractal product [4] | 4.1 × 10−13 | 36.0 | 6.1 × 10−19 | 9.44 | |
| KL10s31a | 31 | N | s31odr10a in Appendix A of [25] | 8.7 × 10−6 | 9.67 | 5.8 × 10−6 | 9.28 | |
| KL10s31b | 31 | N | s31odr10b in Appendix A of [25] | 8.4 × 10−5 | 12.1 | 4.1 × 10−5 | 11.3 | |
| Tsi10s33 | 33 | N | Table II of [26] | 3.3 × 10−6 | 9.33 | 6.4 × 10−7 | 7.93 | |
| SS10s31 | 31 | N | Section 4.4 of [29] | 4.2 × 10−8 | 5.66 | 2.4 × 10−8 | 5.36 | |
| SS10s33 | 33 | N | Section 4.4 of [29] | 1.4 × 10−8 | 5.42 | 8.8 × 10−9 | 5.17 | |
| SS10s35 | 35 | N | Section 4.4 of [29] | 1.0 × 10−9 | 4.41 | 3.1 × 10−11 | 3.11 | |
| Alberdi31 | 31 | N | Appendix A of [48] | 1.5 × 10−7 | 6.44 | 1.0 × 10−7 | 6.20 | |
| Alberdi33 | 33 | N | Appendix A of [48] | 8.4 × 10−8 | 6.47 | 5.5 × 10−8 | 6.20 | |
| Alberdi35 | 35 | N | Appendix A of [48] | 1.5 × 10−8 | 5.79 | 9.5 × 10−9 | 5.52 | |
| PP10s19 | 19 | Y | P1910 in Table 6 of [31] | NA | NA | 6.2 × 10−6 | 5.73 | |
| PP10s23 | 23 | Y | P2310 in Table 6 of [31] | 1.4 × 10−6 | 5.96 | 1.6 × 10−8 | 3.82 | |
| Y10m17 | 35 | N | Table 3 | 1.9 × 10−8 | 5.91 | 6.1 × 10−11 | 3.33 | |
The size of the time step used for the estimation is t = e−5/2 ≈ 0.08, which we found to be small enough for accurate estimation of the constants. An exception is the 10th order product formula PP10s23 in the second-last line of Table 7. That was evaluated for t = e−3/2 ≈ 0.2, rather than t ≈ 0.08 for the other product formulae. Estimates for t ≈ 0.08 are significantly larger, possibly because the coefficients for that product formula were not given to high precision. We provide the estimates of χ, ζ for larger t to evaluate it as fairly as possible.
To compare among product formulae of the same order, we compare Mζ1/k for eigenvalue error, and the best results in each table are highlighted. We also provide values of Mχ1/k for comparison, but these are not used to select the best product formulae. For 8th order it can be seen that our new processed product formula provides the best performance out of all formulae tested. Our two new product formulae of length m = 10 (M = 21) also give better performance than any prior product formulae tested, including processed product formulae. Hence we have given a sequence of three new solutions here, all of which provide improved performance over prior work. There are also other results in our database of over 100,000 solutions that outperform prior work, but we have just given a selection of the best here.
In terms of eigenvalue error, our solution Y8m10b is about 100 times more accurate than the lowest-error prior result, SS8s19. Our processed product formula YP8m8 is about a factor of 300 times more accurate than the lowest-error prior result for processed product formulae, PP8S19, which is longer. These surprisingly large improvements demonstrate that the solution space had not been fully explored in prior work. Each increase in m greatly expands the solution space to explore, and it is likely that there are further improvements still to be found. The error for our processed product formula YP8m8 has slightly larger error than our non-processed formula Y8m10b, but its shorter length means that Mζ1/k is smaller, thus providing the best performance out of all product formulae tested.
These improvements are consistent between the samples of the Hamiltonians. In comparing our formula Y8m10b to the best from prior work SS8s19, our smallest improvement (out of 1,024 samples) was about a factor of 35. A histogram of the improvement for our new product formula is shown in Figure 1. We also find that the improvements are consistent between sizes of the matrices. We have tested sizes 6, 16, and 32 as well as 64 (which is shown in the tables), and found significant improvements at all sizes tested.
A histogram of the ratio of the error of the most accurate prior product formula SS8s19 to our new solution, in the case without processing. The histogram is over the 1024 samples of 64 × 64 matrices. The average improvement is about a factor of 100, with the lowest being about 35.
A recent review [33] also provides a comparative analysis of many of the product formulae in the literature, and for 8th order their recommended methods are SS8s21 in the case without processing, and PP8s13 for the case with processing. In our testing we find that SS8s19 is slightly better than SS8s21; this difference may be because we are testing Hamiltonians here, whereas Ref. [33] considers general matrices (which are non-Hermitian). PP8s19 provides lower error than PP8s13, but PP8s13 provides better performance due to the shorter length (and correspondingly smaller value of Mζ1/k). The three new results we have listed in Table 6 all outperform those methods recommended by Ref. [33].
As a further test, we have calculated the error for simulating the transverse-field Ising model for an 8-qubit chain. We used a time step of 1, with the two components of the Hamiltonian normalised to unit norm, so it is
In the case of 10th order, we find that the best performance is provided by SS10s35, from Sofroniou and Spaletta [29] (which is recommended by Ref. [33]). The performance of that product formula is extraordinarily good, with error nearly a factor of 280 times less than any others we have tested from prior work. Our solution Y10m17 provides intermediate performance, having error twice as large as SS10s35, but a factor of 140 times better than any others. For example, it outperforms the processed product formula PP10s23, which was also recommended by Ref. [33]. The fact that there is a very exceptional solution at one particular length for 10th order suggests that there are better solutions still to be found, motivating further searches in future work.
For 4th order the best result for eigenvalue error is given by the formula BM4M6 of Blanes and Moan [47], which is not in the form of a product of S2. The processed product formula PPBCM4m6 also provides excellent performance, with a marginally larger error in our testing. A recently proposed formula O4M5 of Ostmeyer [46] also gives good performance. In the 6th order case we find that the best performance is given by PPBCM6m9, which is a processed product formula not constructed as a product of S2. The recently proposed processed product formula BCE6m10 also gives excellent performance with low error, but is much longer, resulting in a larger value of Mζ1/k.
In all cases the fractal products provide poor performance. The second fractal product yields very low error but very large M, equal to 125 for 8th order and 625 for 10th order. This is far beyond that for any other product formulae considered here, making this product formula uncompetitive despite the small error. The first fractal product yields large error as well as large M, making its performance by far the worst out of any product formula tested. We also tested various combinations of the two formulae (by using different formulae in the iteration to increase the order) or using fractal formulae to increase the order of lower-order formulae. None of these provided competitive performance, so are not included in the tables.
To consider the threshold to use a higher-order product formula over a lower-order one, we first consider the asymptotic expression using Mζ1/k for eigenvalue error. For 4th order the value of this quantity is 0.42 (BM4M6) and for 6th order is 0.78 (PPBCM6m9), which gives a threshold of 1,700. The value for the best 8th order is 1.24 (our result YP8m8), which gives a threshold of around 6000 when compared with 4th order and 68,000 when compared with 6th order.
At 10th order the best value is 3.11 (SS10s35), which when compared with 8th order gives a threshold of 9.3 × 1015. This is far greater than can be expected for realistic simulations. To see the reason for this, if we use Eq. (46) for the approximate number of exponentials needed for phase estimation (where the threshold is for 1/εH), it would be more than 3 × 1018. In realistic implementations of quantum chemistry it takes about 104 Toffoli gates per exponential, so this would correspond to a Toffoli complexity over 1022. For comparison, in the analysis of Shor's algorithm in Ref. [49], 3 × 109 Toffolis were estimated to take 8 hours. Even with improvements in gate speed, a simulation sufficiently large for 10th order product formulae to provide an advantage would require a computation time of many millions of years.
This threshold between 8th and 10th order can be expected to be accurate (at least for the estimate of the order of magnitude) because it corresponds to small time steps below 0.3. On the other hand, the other thresholds would correspond to large time steps, so would be inaccurate because the asymptotic expressions break down. The calculated threshold between 6th and 8th order would correspond to a time step larger than 3, which is far too large for the scaling of the error to be accurate.
In this case we find numerically that the threshold between these 6th and 8th order formulae (PPBCM6m9 and YP8m8) is increased dramatically, to about 3.2 × 107, corresponding to a time step of about 1.22. If the target ε is 0.01 for example, that would correspond to about 260,000 time steps. It is possible to optimise our kernel for the larger time step size in order to provide better performance. This optimisation provides a kernel (see Table 8) such that the threshold is reduced about a factor of twelve to 2.8 × 106. This now corresponds to 15,000 time steps with a time step size of about 1.82. This is fewer time steps than expected for simulations of classically intractable systems. Thus the optimisation provides another order of magnitude where 8th order product formulae are optimal, albeit with a product formula tailored to the time step size.
The 8th order kernel tailored for large time step size t = 1.82 to improve the threshold.
| kernel for large time step | |
|---|---|
| w1 | 0.1777372900430394 |
| w2 | 0.2862580532195395 |
| w3 | 0.1701306063199336 |
| w4 | −0.3746748008394162 |
| w5 | 0.1485267804844835 |
| w6 | −0.3773225725485588 |
| w7 | 0.5395886879620081 |
| w8 | 0.2210419534887659 |
One must also be careful in considering thresholds based on eigenvalue error because there needs to be a large number of time steps such that the spectral-norm error in the complete evolution time is dominated by the eigenvalue error. There is a factor of about 65 between the spectral-norm error and the eigenvalue error for the 8th order product formula YP8m8. For these thresholds the number of time steps is orders of magnitude larger than this, which indicates that the eigenvalue error is an appropriate measure.
These threshold calculations are for Hamiltonians that are sums of two terms, each of which is normalised. In the case where the Hamiltonian is not normalised, the threshold should be scaled by the norm. In the simple case where the Hamiltonian is a sum of two terms with equal norms ‖A‖ = ‖B‖, then the threshold will be for ‖A‖T/ɛ rather than T/ɛ. More generally the Hamiltonian may be a sum of any number of terms with different norms, which would change the threshold and make it unclear what quantity the threshold should be considered in terms of. There is also the possibility that the threshold can be changed for Hamiltonians that have structure to them instead of being random, for example those for quantum chemistry.
Fermionic Hamiltonians encountered in quantum chemistry often have the form
Let
To fairly compare product formulae as applied to quantum chemistry, we can define ω so that
To analyse the non-asymptotic regime, we rewrite the error to be proportional to
To see why that is the case, note that Lemma 2 of Ref. [50] describes the error at order k + 1 (p + 1 in the notation of that work) in terms of a sum over norms of multicommutators of k + 1 terms in the Hamiltonian. It is easily seen that the error at each of the higher orders will also be given by similar multicommutator forms. The method of proof of Theorem 4 in Ref. [50] is based on simplifying the multicommutator expressions from Lemma 2, and proceeds in exactly the same way for the higher orders.
Therefore, the sum over errors of all orders would give a total error of the form
Then following the same procedure to determine the threshold as before, we will have a total error given by
The threshold can be calculated for Hamiltonians with a given value of ‖τ‖1 + ‖ν‖1,[η], then for some new Hamiltonian with different values of the norms, the threshold should still be for the same value of t1(‖τ‖1 + ‖ν‖1,[η]). This implies that the right-hand side of Eq. (69) should be the same for the threshold with the new Hamiltonian. As a result, we expect that the threshold should be in terms of the left-hand side of Eq. (69), which is the same expression as found in the asymptotic regime.
In our numerical testing, the coefficients τpq, νpq were chosen uniformly from the interval [−1, 1]. We computed ω for a selection of the product formulae with the best performance for random Hamiltonians. In Table 9, Table 10, and Table 11 we give the computed result for d = 6 orbitals, assuming half-filling of the orbitals. Our numerics indicate that the error is roughly proportional to the bound in Eq. (61), independent of the number of orbitals, though the constant factors are small. In Appendix C we give the constant factors for the case d = 4, showing that the computed constant does not change much with a different d.
Comparison of constant factors ω for a selection of the lowest-error best product formulae for 6th order. We generate 1, 000 random Hamiltonians with d = 6 orbitals as in Eq. (60) and compute the average ω.
| label | M | processing | reference | ω | Mω1/k |
|---|---|---|---|---|---|
| PPBCM6m9 | 9 | Y | P96 Table 5 of [31] | 2.8 × 10−9 | 0.33 |
| PPBCM6m6 | 13 | Y | P136 in Table 6 of [31] | 1.2 × 10−9 | 0.42 |
| BCE6m10 | 20 | Y | 3.4 × 10−11 | 0.36 | |
Comparison of constant factors ω for a selection of the lowest-error product formulae for 8th order. We generate 1, 000 random Hamiltonians with d = 6 orbitals as in Eq. (60) and compute the average ω.
| label | M | processing | reference | ω | Mω1/k |
|---|---|---|---|---|---|
| SS8s19 | 19 | N | Section 4.3 of [29] | 3.5 × 10−11 | 0.94 |
| PP8s13 | 13 | Y | P138 in Table 6 of [31] | 4.4 × 10−10 | 0.88 |
| Y8m10 | 21 | N | Table 1 (our new result) | 8.5 × 10−12 | 0.87 |
| Y8m10b | 21 | N | Table 1 (our new result) | 1.3 × 10−12 | 0.68 |
| YP8m8 | 17 | Y | Table 2 (our new result) | 1.7 × 10−12 | 0.57 |
Comparison of constant factors ω for a selection of the lowest-error product formulae for 10th order. We generate 1, 000 random Hamiltonians with d = 6 orbitals as in Eq. (60) and compute the average ω.
| label | M | processing | reference | ω | Mω1/k |
|---|---|---|---|---|---|
| SS10s35 | 35 | N | Section 4.4 of [29] | 3.0 × 10−15 | 1.24 |
| PP10s23 | 23 | Y | P2310 in Table 6 of [31] | 2.3 × 10−12 | 1.58 |
| Y10m17 | 35 | N | Table 3 | 2.6 × 10−14 | 1.53 |
From the results in these tables, we see that the relative performance for the product formulae in unchanged from the results with the random matrices. This is a useful independent verification that the assessment of the performance of the product formulae is relatively independent of the class of matrices. In particular, we find that for 8th order our three new results all outperform the best product formulae in prior work (where we have compared to the two prior formulae with best performance). In the case of 10th order the product formula SS10s35 of Ref. [29] still provides exceptional performance.
For calculating the thresholds based on the asymptotic formula in Eq. (63), we should compare our best 8th order formula YP8m8 to SS10s35 in the 10th order case, and PPBCM6m9 in the 6th order case. Then we obtain a threshold of 3 × 1013 for 8th to 10th order, and 5 × 105 for 6th to 8th order, indicating about 9 orders of magnitude for T /ε where our 8th order product formula is optimal.
Similar to the case for general matrices, we expect that the threshold for 8th to 10th order is sufficiently large that the time steps are small and the asymptotic formulae are accurate, but the threshold between 6th and 8th order may be unreliable. In that case we solve Eq. (68) for t1 with t2 = t1M2/M2, and find the threshold using Eq. (69). We expect that the 6th order product formula PPBCM6m9 can give better performance for larger time steps, and indeed we find that the threshold is further increased to 1.7 × 107. By using our 8th order kernel optimised for large time step in Table 8 we can improve the threshold between PPBCM6m6 and our 8th order formula in the non-asymptotic case. We then obtain a threshold of around 3 × 106, which is an improvement by a factor of 5 times.
Next we compare these thresholds to the parameters for specific systems. From Ref. [50] the norms can be expected to scale as
For comparison, we also consider the Alpha + Hydrogen system from Ref. [51]. The parameters η, N , and Ω are given in Table 1 of that work as 218, 533, and 2420, respectively. Some of the most challenging values of T , ε are given in Figure 4 of that work, with for example T = 40 and ε = 10−4. With these parameters the LHS of Eq. (63) is approximately 8.6 × 109, which is below the threshold (between 8th and 10th order) by a few orders of magnitude. The smallest values of T /ε from Ref. [51] are obtained with T = 10 and ε = 0.01 (see Table IV). That results in a value for the LHS of Eq. (63) of 2.2 × 107, which is above the threshold between 6th and 8th order, so our 8th order product formula would still be optimal. This indicates that the new 8th order product formula found here should be optimal for the full range of parameter values used in Ref. [51]. Nonetheless, for smaller values of T /ε it may be advantageous to further customise the product formula to provide better performance at larger time steps, as described above to improve the threshold between 6th order and 8th order.
Note also that the fact that the threshold here is different than for random Hamiltonians means that the threshold will depend on the class of Hamiltonians. Selecting a different type of Hamiltonians can change the threshold. Similarly, it will be expected that choosing a particular class of fermionic Hamiltonians (such as those arising from plane waves) rather than random fermionic Hamiltonians can give a different threshold. Nevertheless, changing the class of Hamiltonians still resulted in an extremely large threshold for 10th order product formulae to be optimal. This indicates that the 8th order formula will be optimal for realistic simulations regardless of the class of Hamiltonians.
We have found 8th order product formulae with significantly lower error than in prior work. For the case of processed product formulae, our best solution is about 300 times more accurate than that in prior work, and in the case of product formulae without processing our solution is about 100 times more accurate. The size of the improvement is estimated based on testing with random 64 × 64 matrices. We have also tested on a range of different matrix sizes, fermionic Hamiltonians of different sizes, and a range of sizes of the transverse-field Ising model and find consistent improvements with well over an order of magnitude in all cases.
The size of the improvement is surprising, especially considering the best prior results are from 2005 (without processing) and 2006 (with processing). To achieve this improvement we have used distinct techniques from prior work. The processed product formulae from Ref. [31] set multiple coefficients equal to one another. Specifically, in the P198 formula, six of the coefficients are identical, following a heuristic suggested by McLachlan [27]. However, our results suggest that using distinct values for each coefficient yields better outcomes. The advantage of taking longer product formulae with distinct coefficients is that there are significantly more parameters than constraints, giving more freedom for optimisation of the product formula. The solutions which we find give the best performance have no clear pattern to the coefficients, indicating that assuming these patterns is too restrictive.
The way that product formulae were evaluated in Ref. [31] was based on the ℓp norms of the coefficients for p = 1 and p = k + 1. Those norms are both smaller for P198 given in Ref. [31] than our lower-error solution. Similarly, these norms are both smaller for the solutions from Ref. [29] than our non-processed product formula. Also, Ref. [48] found 10th order integrators with 1-norm lower than those in Ref. [29], but in our testing their error is larger. In all these cases the solutions with lower error have larger norms. The ℓk+1 norm corresponds to just one of the 𝒪(tk+1) error terms, and a product formula that minimises that term while giving large values for the others will not perform well. We find that the strategy of selecting product formulae with low error for examples can yield much better results than attempting to minimise norms.
In this work we compare the eigenvalue error, rather than the spectral-norm error as in prior work because it is the eigenvalue error in a single step which dominates the error over longer time evolutions, with basis error cancelling. In many applications of quantum computing the goal is to estimate eigenvalues, so it is solely the eigenvalue error which should be considered. For most product formulae tested, the spectral norm of the error (for a single step) is substantially larger than the eigenvalue error, meaning that it can be an unrealistic measure of accuracy. This also suggests that techniques for reducing error based on randomisation of the term order [21] or randomised corrections [53] may provide less improvement than anticipated because if successive steps are not identical then the basis error will not cancel.
We have shown how to fairly compare the performance of product formulae with different numbers of factors and different orders. The optimal order of integrator to use depends primarily on the ratio T /ε, the ratio of the total evolution time to the allowed error. This is because, as T /ε becomes larger, the error allowed for each individual time step becomes smaller, and the time steps must become shorter to make the error sufficiently small. That results in a larger overall number of time steps, and so a larger complexity. The higher-order integrators better reduce the error with shorter time steps, and so for large T /ε that improvement more than compensates for the larger number of exponentials for the higher-order product formulae. For eigenvalue estimation there is an equivalent threshold in 1/εH , in terms of the allowed error in the eigenvalue.
The threshold for 10th order product formulae to outperform 8th order product formulae is very large, so 10th order would only be optimal for simulations well beyond the scale considered for quantum algorithms. Similarly, it would not be useful to consider 12th or higher order product formulae for the purpose of Hamiltonian simulation, as the threshold would be even larger. For large quantum simulations that are still reasonable to perform on a quantum computer, such as those considered in Refs. [51] or [52], the best performance is provided by our 8th order product formulae (out of the extensive set of product formulae tested). For more moderately sized simulations 6th order will be optimal, and in our testing we find that the product formula S96 of Ref. [31] provides excellent performance. For smaller simulations, as in NISQ experiments, 4th order product formulae would be ideal, with S6 of Ref. [47] and P64 of [31] providing the best performance.
Care is required in determining the thresholds between the orders. Using the formula based on the asymptotic scaling can yield estimates of the T /ε threshold corresponding to large time steps where the asymptotic scaling no longer holds. When properly accounting for the error for a given time step, the T /ε threshold can change by many orders of magnitude. This can also be a problem more generally in estimating the complexity of an algorithm. Estimating the number of time steps based on the asymptotic expression may correspond to unrealistically large time steps, so the asymptotic formula is inaccurate. We find that significant performance gains can be obtained by fine-tuning a product formula for larger time steps.
A potential topic for future work is further customising product formulae for large time steps. We have customised the 8th order product formula to improve the threshold, but it would also be possible to customise the 6th and 4th order product formulae. It would also be possible to customise the product formulae for larger time steps with the fermionic Hamiltonians. This approach could be used to provide better performance for small-scale simulations suitable for NISQ devices. It could also potentially extend the parameter range where 6th order product formulae provide best performance to larger simulations on fault tolerant quantum computers.