Optimize quantum simulation using a force-gradient integrator

Quantum simulation as proposed by Feynman [1] has shown great potential in many fields due to its powerful computational capabilities and their potential to avoid the fermionic sign problem [2–4] which is an NP hard problem [5]. Consequently, in recent years quantum simulation has developed rapidly and has been applied in many fields [6–9]. With Google's announcement of quantum "supremacy" [10] and the subsequent experiment using photons [11], quantum simulations hold great promise for the future.

One of the biggest problems plaguing quantum simulations is the lack of fidelity. Recently, the fidelity can reach up to $99.9\%$ for one-qubit gate and $94.7\%$ for Clifford gate [12]. The fidelity decays exponentially with the number of gate operations, as a result the required computing resources increase exponentially, and the computational power of a quantum computer will be unable to surpass that of classical computers [13]. This problem has led to widespread researches in quantum error correction [14–21].

The Trotter (or Lie-Trotter-Suzuki) decomposition (TD) is one of the earliest quantum algorithms in quantum simulations [22]. TD belongs to the product formulas, although some post-Trotter methods have been proposed, the product formulas are still found to be competitive especially for the case in which the simulated system has a Lie-algebraic structure, and the error of product formulas has been well studied [23]. A class of optimizations of TD is the high-order TDs [24,25], for example the second-order (or symmetric) TD (STD) has been widely applied [26–30]. Until very recently, TD was still used in simulations [31–34]. The optimization of TD and high-order TDs is very important, since the fidelity decays exponentially, even modest optimization can make a drastic improvement, and may enable some models to be simulated under existing conditions. In addition, the real quantum computers are not yet widespread and many simulations are performed on quantum computer simulators running on a classical computer [35]. The optimization can accelerate these simulations significantly as well.

There is a problem in the field of lattice QCD similar to the improvement on TD and high-order TDs. The STD with two non-commutating terms can correspond to the leapfrog integrator in lattice QCD. However, there are integrators much faster than leapfrog such as the Omelyan integrator [36] and the force-gradient integrators [37–39]. The latter do not belong to the class of high-order TDs. There has been prior work on optimizing the simulations of quantum systems using analytically derived gradients, such as the GRAPE algorithm [40–43]. However, to our knowledge the force-gradient integrators have barely caught the attention of the quantum simulation community so far. The force-gradient integrators can be applied in quantum simulations alongside other optimizations used in product formulas [30,44–46], and with the nested integrator techniques [47]. In this letter, we investigate the feasibility of a force-gradient decomposition (FGD) in quantum simulations.

Considering a system whose Hamiltonian can be written as $H=S+T$, with $[ S,T] \neq 0$, such a system can be simulated using TD, i.e., the $\exp{(\textrm{i}Ht)}$ is approximated by $\exp (\textrm{i}Ht) \approx(\exp (\textrm{i} \tau S) \exp (\textrm{i} \tau T))^{m}$ with $\tau=t/m$. When decomposed to m steps, the total number of exponential operations (denoted as n) required is $n=2m$. Its error can be roughly estimated by the Baker-Campbell-Hausdorff formula. Assume $\tau \lt 1$ and $[ S,T] \sim \mathcal{O}(l)$, the error when evolute to t is $\epsilon \sim \mathcal{O}( t^2 l /2m )$. Similarly, the STD is $\exp (\textrm{i}Ht) \approx (\exp (\textrm{i}\tau S/2) \exp (\textrm{i}\tau T) \exp (\textrm{i}\tau S/2))^m$ with $n=2m+1$ and $\epsilon \sim \mathcal{O}( t^3 l^2 /8m^2 )$. The STD is in the category of higher-order TDs which has the form ${\rm e}^{t_1S}{\rm e}^{t_2T}{\rm e}^{t_3S}{\rm e}^{t_4T}\ldots {\rm e}^{t_MS}$ [24,25], the FGDs are no longer higher-order TDs. A widely used FGD is [39]

$$\begin{align} &\exp \left\{\textrm{i}m \tau\left(S+T+\mathcal{O}\left(\tau^{4}\right)+\mathcal{O}\left(\tau^{6}\right)\right)\right\}=\nonumber\\ &\left\{\exp \left(\frac{\textrm{i}\tau S}{6}\right) \exp \left(\frac{\textrm{i}\tau T}{2}\right) \exp \left(\frac{2\textrm{i} \tau S}{3}+\frac{\textrm{i}\tau^{3}}{72}[S,[S, T]]\right) \right.\nonumber\\ &\times \left.\exp \left(\frac{\textrm{i}\tau T}{2}\right) \exp \left(\frac{\textrm{i}\tau S}{6}\right)\right\}^{m}, \end{align} \tag{ 1 }$$

with

$$\begin{equation} \begin{split} &\mathcal{O}\left(\tau^{4}\right)=-\frac{\tau^{4}}{155520}\left\{41\left[S,[S,[S,[S, ~T]]]\right]\right.\\ &\left.+36[[S, ~T],[S[S, T]]]+72[[S, ~T],[T,[~S, T]]]\right.\\ &\left.+84[T,[S,[S,[S, ~T]]]]+126[T,[T,[S,[S, ~T]]]]\right.\\ &\left.+54[T,[~T,[T,[S, ~T]]]]\right\}. \end{split} \end{equation} \tag{ 2 }$$

Usually, $[S, [S,[S,T]]]\neq 0$, therefore one needs to further decompose $\exp (2\textrm{i}\tau S/3 +\textrm{i} \tau^3 [S,[S,T]]/72 )$. Both n and depend on how this term is decomposed. Nevertheless, as will be shown, eq. (1) can be usually decomposed as the m-th power of the product of 7 exponents (denoted as 7-stage decomposition). For l-stage decomposition, $n=(l-1)m$.

A general discussion of the error of FGD is beyond our ability, and is model dependent. Instead, we study the optimization brought by the FGD using two specific applications as examples.

To concentrate on the error due to the decomposition, we define $\varepsilon= \sqrt{||\exp(\textrm{i}tH)-M||}/\sqrt{||\exp(\textrm{i}tH)||}$, where $||\ldots||$ denotes the sum of squares of the matrix elements and the matrix M denotes the decomposed matrix, for example $M=(\exp{(\textrm{i}\tau S)}\exp{(\textrm{i}\tau T)})^m$ in the case of TD.

For a small system, the Hamiltonian can be numerically diagonalized, and ε can be calculated on a classical computer. We use two small systems to compare different decompositions. The Omelyan integrator is frequently used in lattice QCD and is also included as an example of high-order TDs, which is denoted as Omelyan decomposition (OD) [36],

$$\begin{align} \exp (\textrm{i}tH) &\approx\!\left\{\exp \left(\textrm{i} \alpha \tau S\right) \exp \left(\frac{\textrm{i}\tau T}{2}\right) \right.\nonumber\\ &\quad\times\! \left.\exp \left(\textrm{i} (1-2\alpha)\tau S\right)\exp \left(\frac{\textrm{i}\tau T}{2}\right)\exp \left(\textrm{i} \alpha \tau S\right)\!\right\}^{m}\!, \end{align} \tag{ 3 }$$

with $\alpha \approx 0.1931833275037836$.

The FGD in eq. (1) is usually a 7-stage decomposition, therefore we also compare FGD with a 7-stage high-order Trotter decomposition (7TD). An optimized 7TD is [38]

$$\begin{align} \exp \left(\textrm{i}tH\right)&\approx \left\{\exp \left(\textrm{i}\frac{\beta}{2} \tau S\right) \exp \left(\textrm{i}\beta \tau T\right) \exp \left(\textrm{i}\frac{1-\beta}{2}\tau S\right)\right.\nonumber\\ &\quad\left.\times \exp \left(\textrm{i}(1-2\beta) \tau T\right) \exp \left(\textrm{i}\frac{1-\beta}{2}\tau S\right)\right.\nonumber\\ &\quad\left.\times \exp \left(\textrm{i}\beta \tau T\right) \exp \left(\textrm{i}\frac{\beta}{2} \tau S\right)\right\}^{m},\nonumber\\ \beta &=\frac{1}{2-\sqrt[3]{2}}. \end{align} \tag{ 4 }$$

Since the number of gate operations is approximately proportional to n for a same model, we use n to quantify the complexity of the decompositions.

Transverse Ising model

The transverse Ising model (TIM) is a popular benchmark for quantum simulations. Therefore we also use TIM in 1D and 2D to test the FGD.

The Hamiltonian of TIM can be written as

$$\begin{equation} H=\sum_{\langle ij\rangle} \sigma_{z}(n_i) \sigma_{z}(n_j)+\lambda \sum_i \sigma_{x}(n_i), \end{equation} \tag{ 5 }$$

where $\sigma_{x,z}$ are the Pauli matrices, $\langle ij\rangle$ refers to the nearest-neighbouring pairs. The $\sigma_z(n_i)$ is short for the tensor product of $2\times 2$ matrices on each site, the matrix on n_i is $\sigma_z$, and the others are identity matrices.

Firstly, we consider a 1D transverse Ising chain with only 3 sites and with a periodic boundary condition. The Hamiltonian in this case is a $8\times 8$ matrix acting on 3 qubits. Denoting the indices of the sites as n_{1, 2,3}, $H=T+S$ with

$$\begin{equation} \begin{split} &T=\sigma_{z}(n_{1}) \sigma_{z}(n_{2})+\sigma_{z}(n_{2}) \sigma_{z}(n_{3})+\sigma_{z}(n_{3}) \sigma_{z}(n_{1}),\\ &S=\lambda \sum_{i=1}^3 \sigma_{x}(n_{i}). \end{split} \end{equation} \tag{ 6 }$$

Therefore $[S,[S, T]]=-8 \lambda^{2}(Y-T)$ with

$$\begin{equation} Y=\sigma_{y}(n_{1}) \sigma_{y}(n_{2})+\sigma_{y}(n_{2}) \sigma_{y}(n_{3})+\sigma_{y}(n_{3}) \sigma_{y}(n_{1}).\quad \end{equation} \tag{ 7 }$$

For the $\exp (2\textrm{i}\tau S/3 +\textrm{i} \tau^3 [S,[S,T]]/72 )$ term we use the STD which is found to be optimized among 5-stage decompositions for 3 non-commutating terms [48]. The FGD is then

$$\begin{equation} \begin{split} \exp (\textrm{i}tH) &\approx\left\{\exp \left(\frac{\textrm{i} \tau S}{6}\right) \exp \left(\textrm{i} T\left(\frac{\tau}{2}+\frac{\tau^{3} \lambda^{2}}{18}\right)\right) \right.\\ &\quad\left.\times \exp \left(\frac{\textrm{i} \tau S}{3}\right)\exp \left(-\textrm{i} \frac{\tau^{3} \lambda^{2}}{9} Y\right) \exp \left(\frac{\textrm{i} \tau S}{3}\right) \right.\\ &\quad\left.\times \exp \left(\textrm{i} T\left(\frac{\tau}{2}+\frac{\tau^{3} \lambda^{2}}{18}\right)\right) \exp \left(\frac{\textrm{i} \tau S}{6}\right)\right\}^{m}. \end{split} \end{equation} \tag{ 8 }$$

Similar to the TD, before each exponential evaluation, the qubits should be rotated to $\sigma_{x,y,z}$ representations accordingly. Note that there is a Y operator which is diagonalized in $\sigma_y$ representation which is different from TD and high-order TDs.

In the case of infinite volume, the 1D TIM is self-dual, there is a phase transition at $\lambda_c=1$ [49,50]. Therefore we use $\lambda = 0.5, 1, 1.5$ as examples. When t = 1, we calculate the minimum number of exponential operations needed (denoted as $n_{\min}$) when the error satisfies $\varepsilon \lt 0.1\%$. The results are shown in table 1.

Table 1:. The $n_{\min}$ required to satisfy $\varepsilon \lt 0.1\%$ for the Ising chain when t = 1.

λ		TD	STD	OD	7TD	FGD
0.5	$n_{\min}$	1034	39	33	43	19
	$\varepsilon(\%)$	0.10	0.098	0.086	0.061	0.033
1	$n_{\min}$	1606	55	53	55	25
	$\varepsilon(\%)$	0.10	0.098	0.099	0.088	0.035
1.5	$n_{\min}$	1456	71	69	67	31
	$\varepsilon(\%)$	0.10	0.095	0.094	0.092	0.048

The non-commutativity l grows with λ, therefore the worst case is when $\lambda =1.5$. When $\lambda =1.5,t=1$, the error decays for the decompositions are shown in fig. 1. The results indicates that the simulation of 1D TIM can be significantly optimized by the FGD.

Zoom In Zoom Out Reset image size

Once the FGD of 1D TIM is known, the FGD for higher-dimensional TIM with periodic boundary condition can be obtained because $T=\sum T_{1D}$, and $[S,[S,T]]=\sum [S,[S, T_{1D}]]$.

Take the TIM on a $2\times 3$ lattice with periodic boundary condition as an example, the lattice is depicted in fig. 2, and the Hamiltonian is a $64\times 64$ matrix with

$$\begin{equation} \begin{split} S&=\sum_{i=1}^6 \sigma_x (n_i),\quad T=\sum_{j=1}^5 T_j,\quad Y=\sum_{j=1}^5Y_j,\\ T_1&=\sigma_z (n_1)\sigma_z (n_2)+\sigma_z (n_2)\sigma_z (n_3)+\sigma_z (n_3)\sigma_z (n_1),\\ T_2&=\sigma_z (n_4)\sigma_z (n_5)+\sigma_z (n_5)\sigma_z (n_6)+\sigma_z (n_6)\sigma_z (n_4),\\ T_3&=2\sigma_z (n_1)\sigma_z (n_4),\quad~ T_4=2\sigma_z (n_2)\sigma_z (n_5),\\ T_5&=2\sigma_z (n_3)\sigma_z (n_6),\\ Y_1&=\sigma_y (n_1)\sigma_y (n_2)+\sigma_y (n_2)\sigma_y (n_3)+\sigma_y (n_3)\sigma_y (n_1),\\ Y_2&=\sigma_y (n_4)\sigma_y (n_5)+\sigma_y (n_5)\sigma_y (n_6)+\sigma_y (n_6)\sigma_y (n_4),\\ Y_3&=\sigma_y (n_1)\sigma_y (n_4),\qquad Y_4=\sigma_y (n_2)\sigma_y (n_5),\\ Y_5&=\sigma_y (n_3)\sigma_y (n_6).\\ \end{split} \end{equation} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

It can be verified that $[S,[S,T]]=-8\lambda^2 (Y-T)$ holds. Therefore one can still use the FGD in eq. (8). The critical transverse field $\lambda_c$ for 2D TIM is found to be $2\text{--}4$ [50–54], therefore we choose $\lambda = 1.5, 3, 5$ as examples. When t = 1, $n_{\min}$ required when the error satisfies $\varepsilon \lt 0.1\%$ are listed in table 2. When $\lambda =5,t=1$, the error decays for the decompositions are shown in fig. 3. Similar to the case of 1D TIM, the FGD can significantly reduce the number of gate operations.

Table 2:. The $n_{\min}$ required to satisfy $\varepsilon \lt 0.1\%$ for the Ising chain when t = 1.

λ		TD	STD	OD	7TD	FGD
1.5	$n_{\min}$	3788	121	125	109	37
	$\varepsilon(\%)$	0.10	0.10	0.094	0.088	0.099
3	$n_{\min}$	4042	187	197	157	55
	$\varepsilon(\%)$	0.10	0.098	0.098	0.089	0.091
5	$n_{\min}$	4564	243	257	211	79
	$\varepsilon(\%)$	0.099	0.10	0.10	0.096	0.081

Zoom In Zoom Out Reset image size

Ising gauge model

The $D=d+1$ dimensional Ising gauge model, i.e., $\mathbb{Z}_2$ lattice gauge model [55] at finite temperature can be dual to a d-dimension Ising gauge model in a transverse field at zero temperature [49,56], also known as one of the quantum link models [57]. The Hamiltonian can be written as

$$\begin{equation} H=\sum_{\square} \prod_{l\in \square}\sigma_{z}(l) +k \sum_{l} \sigma_{x}(l), \end{equation} \tag{ 10 }$$

where l is an index of a link, and "□ " denotes a plaquette. We consider the case of a stripe containing two plaquettes with periodic boundary condition on major direction as shown in fig. 4. The Hamiltonian is a $64\times 64$ matrix with

$$\begin{equation} \begin{split} S&=\prod_{i=1,4,5,3} \sigma_z(l_i)+\prod_{i=2,3,6,4} \sigma_z(l_i),\\ T&=\lambda \sum_{i=1}^6 \sigma_{x}\left(l_{i}\right),\qquad [S,[S,T]]=A+B\\ A&=4\lambda \left(\sum_{i=1,2,5,6}\sigma_{x}\left(l_{i}\right)+2\sigma_{x}\left(l_{3}\right)+2\sigma_{x}\left(l_{4}\right)\right),\\ B&=8\lambda \left(\prod_{i=1,2,5,6} \sigma_z(l_i)\right)\left(\sigma_{x}\left(l_{3}\right)+\sigma_{x}\left(l_{4}\right)\right). \end{split} \end{equation} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

The FGD is

$$\begin{equation} \begin{split} \exp (\textrm{i}tH) &\approx\left\{\exp \left(\frac{\textrm{i} \tau S}{6}\right) \exp \left(\frac{\textrm{i}\tau}{2}T+\frac{\textrm{i}\tau^3 A}{144}\right) \right.\\ &\quad\left.\times \exp \left(\frac{\textrm{i} \tau S}{3}\right) \exp \left(\textrm{i} \frac{\tau^{3} \lambda}{72} B\right) \exp \left(\frac{\textrm{i} \tau S}{3}\right)\right.\\ &\quad\left.\times \exp \left(\frac{\textrm{i}\tau}{2}T+\frac{\textrm{i}\tau^3 A}{144}\right) \exp \left(\frac{\textrm{i} \tau S}{6}\right)\right\}^{m}. \end{split} \end{equation} \tag{ 12 }$$

In the infinite volume, the 2D Ising gauge model in transverse field can be dual to 2D TIM or $D=2+1$ Ising model at finite temperature [49,58]. The $1/k$ corresponds to λ in eq. (5). The $\lambda_c$ for 2D TIM is $2\text{--}4$, indicating a phase transition at $k_c= 0.25\text{--}0.5$, which is also a topological phase transition [59]. Recently, the 2D Ising gauge model is studied by quantum simulation on a $3\times 3$ lattice with periodic boundary condition and k_c is found to be 0.380 [60]. Therefore we use $k =0.1,0.3,1$ as examples. When t = 1, the $n_{\min}$ needed when $\varepsilon \lt 0.1\%$ are shown in table 3. For the case of the largest k, the error decays are shown in fig. 5. Again, the FGD can significantly reduce the number of gate operations.

Table 3:. The $n_{\min}$ required to satisfy $\varepsilon \lt 0.1\%$ for the Ising gauge model when t = 1.

k		TD	STD	OD	7TD	FGD
0.1	$n_{\min}$	376	15	13	19	13
	$\varepsilon(\%)$	0.10	0.092	0.074	0.035	0.036
0.3	$n_{\min}$	1012	29	25	31	19
	$\varepsilon(\%)$	0.10	0.094	0.087	0.067	0.022
1.0	$n_{\min}$	1590	63	53	67	25
	$\varepsilon(\%)$	0.10	0.094	0.10	0.090	0.10

Zoom In Zoom Out Reset image size

Summary of the applications

Quantum simulation is an extremely promising research direction extensively studied recently due to the rapid development of quantum computing technology. Reducing the number of decomposition steps in the product formulas paves a way to practical quantum simulation. While the TD and high-order TDs are widely applied in quantum simulations, we show that significant optimization can be achieved by the force-gradient integrator used in lattice QCD. We use the TIM and the Ising gauge model as examples, the FGD can reduce the number of gate operations up to about a third of those using high-order TDs. In addition, it can be seen that if one wishes to use FGD, the $\exp (2\textrm{i}\tau S/3 +\textrm{i} \tau^3 [S,[S,T]]/72 )$ term needs to be processed, which is sometimes not easy to handle. Nevertheless, FGD shows great advantages in quantum simulations and deserves further studies for various of models.

This work was partially supported by the National Natural Science Foundation of China under Grant No. 12047570.