Improved Time-independent Hamiltonian Simulation

Nhat A. Nghiem nhatanh.nghiemvu@stonybrook.edu Department of Physics and Astronomy, State University of New York at Stony Brook, Stony Brook, NY 11794-3800, USA C. N. Yang Institute for Theoretical Physics, State University of New York at Stony Brook, Stony Brook, NY 11794-3840, USA

Abstract

We describe a simple method for simulating time-independent Hamiltonian $H$ that could be decomposed as $H=\sum_{i=1}^{m}H_{i}$ where each $H_{i}$ can be efficiently simulated. Approaches relying on product formula generally work by splitting the evolution time into segments, and approximate the evolution in each segment by the evolution of composing Hamiltonian $H_{i}$ . This key step incur a constraint, that prohibits a (poly)logarithmic scaling on approximation error. We employ the recently introduced quantum singular value transformation framework to utilize the ability to simulate $H_{i}$ in an alternative way, which then allows us to construct and simulate the main Hamiltonian $H$ with polylogarithmical scaling on the inverse of desired error, which is a major improvement with respect to product formula approaches.

I Introduction

Quantum simulation has stand as one of central application of quantum computer. As envisioned by Feynman feynman2018simulating , some quantum mechanism is necessary to simulate nature, as nature is essentially quantum. Since then, major progress have been made. Some pioneering works include lloyd1996universal , where the author uses Lie-Trotter formula to show that a universal quantum simulator is possible, which means that the ability to simulate some systems can translate directly to the simulation of other systems. In aharonov2003adiabatic , the authors show how a system characterized by a sparse Hamiltonian could be efficiently simulated. In a series of works berry2007efficient ; berry2014high ; childs2021theory , methods based on sparse-access Hamiltonian and product formulas have been constructed. Alternatively, the work in childs2010relationship ; berry2015hamiltonian uses quantum walk technique developed in ambainis2003quantum ; ambainis2007quantum to simulate certain Hamiltonian. Model based on linear combination of unitaries have been constructed in childs2012hamiltonian ; berry2015simulating . Most recently, quantum signal processing was proposed low2017optimal ; low2019hamiltonian as a very fruitful framework for simulating time-independent Hamiltonian. These aforementioned works involve time-independent regime, and in parallel attempts, simulation of time-dependent Hamiltonian have been constructed poulin2011quantum ; berry2020time ; chen2021quantum ; an2022time ; low2018hamiltonian ; kieferova2019simulating .

An important aspect of simulation is the complexity, which measures the number of elementary operations, such as single and two qubits gates, required to achieve desired simulation. For approximation purpose, another important metric is the scaling of resources on the desired error. A common technique to simulate given Hamiltonian used in previous simulation results is Lie-Trotter product formula and its higher order extension due to Suzuki suzuki1991general . The most thorough analysis of scaling of product formula’s approach for simulating Hamiltonian is provided in childs2021theory . We remark that while product formula’s approach allows a straightforward way to simulate a system, given the ability to simulate “its constituents”, the key step in such method is breaking the time interval into multiple segments, and in fact, it induces a fundamental constraints to the complexity, as the number of operation need to be executed needs to be as large as the number of such segments. As we shall show subsequently, it leads to the scaling on error term being at least sublinear.

In this work, motivated by the hurdle of such constraint, we observe an alternative route to remove it. Our method is driven by the recently introduced quantum singular value transformation (QSVT) framework, that allows algebraic quantum operations could be implemented in a simple, effective manner. The essential idea from QSVT is that, if we can build a so-called block encoding of some operator, then we can implement a wide range of function of such operator. As application to quantum simulation, if somehow we can block encode the Hamiltonian $H$ of interest, then we can implement its exponent $\exp(-iHt)$ for some $t$ by making use of the so-called Jacobi-Anger expansion. Building upon such an idea, we would show in this work that the ability to simulate constituting Hamiltonian can be used in a different way, that allows us to construct the block encoding of main Hamiltonian $H$ directly. Hence, it results in a framework with scaling that is (poly)logarithmical in inverse of desired error, which is superpolynomial improvement to the popular product formula approach.

The structure of this paper is as follows. In section II, we summarize some preliminaries, including crucial definitions and technical tools of our work. Section III is devoted to our main framework. We beginning in section III by describing the product formula, where we would see the fundamental origin of the constraint we mentioned above. Then, in the next part of the same section III, we outline our method, which is built from the so-called logarithmic of unitary tool developed in gilyen2019quantum . We close and discuss our work in section IV.

II Preliminaries

In this section we recapitulate key recipes that would be central to our construction. We refer the readers to original works gilyen2019quantum for greater details.

Definition 1

Suppose that $A$ is an s-qubit operator, $\alpha,\epsilon\in\mathbb{R}_{+}$ and $a\in\mathbb{N}$ , then we say that the $(s+a)$ -qubit unitary $U$ is an $(\alpha,a,\epsilon)$ -block encoding of $A$ , if

||A-\alpha(\bra{0}^{\otimes a}\otimes\mathbb{I})U(\ket{0}^{\otimes a}\otimes\mathbb{I})||\leq\epsilon

Equivalently, in matrix representation, $U$ is said to be a block encoding of $A/\alpha$ if $U$ has the form

\displaystyle U=\begin{pmatrix}\frac{A}{\alpha}&\cdot\\ \cdot&\cdot\\ \end{pmatrix}.

where $(.)$ in the above matrix refers to irrelevant blocks that could be non-zero.

Lemma 1 (Theorem 2 in rattew2023non )

Given an n-qubit quantum state specified by a state-preparation-unitary $U$ , such that $\ket{\psi}_{n}=U\ket{0}_{n}=\sum^{N-1}_{k=0}\psi_{k}\ket{k}_{n}$ (with $\psi_{k}\in\mathbb{C}$ and $N=2^{n}$ ), we can prepare an exact block-encoding $U_{A}$ of the diagonal matrix $A={\rm diag}(\psi_{0},...,\psi_{N-1})$ with $\mathcal{O}(n)$ circuit depth and a total of $\mathcal{O}(1)$ queries to a controlled- $U$ gate with $n+3$ ancillary qubits.

Lemma 2 (camps2020approximate )

Given the unitary block encoding $\{U_{i}\}_{i=1}^{m}$ of multiple operators $\{M_{i}\}_{i=1}^{m}$ (assumed to be exact encoding), then, there is a procedure that produces the unitary block encoding operator of $\bigotimes_{i=1}^{m}M_{i}$ , which requires a single use of each $\{U_{i}\}_{i=1}^{m}$ and $\mathcal{O}(1)$ SWAP gates.

Lemma 3 (Linear combination of block-encoded matrices)

Given unitary block encoding of multiple operators $\{M_{i}\}_{i=1}^{m}$ and a unitary that prepares the state $\ket{y}=\sum_{i=1}^{m}\sqrt{sign(y_{i})y_{i}/\beta}\ket{i}$ where $\beta=\sum_{i=1}^{m}|y_{i}|$ . Then, there is a procedure that produces a unitary block encoding operator of $\sum_{i=1}^{m}y_{i}M_{i}/\beta$ in complexity $\mathcal{O}(m)$ , using block encoding of each operator $M_{i}$ a single time.

Lemma 4 (Block Encoding of Product of Two Matrices)

Given the unitary block encoding of two matrices $A_{1}$ and $A_{2}$ (assuming to have norm less than 1), then there exists an efficient procedure that constructs a unitary block encoding of $A_{1}A_{2}$ using each block encoding of $A_{1},A_{2}$ one time.

Lemma 5

[gilyen2019quantum Theorem 30] Let $U$ , $\Pi$ , $\widetilde{\Pi}\in{\rm End}(\mathcal{H}_{U})$ be linear operators on $\mathcal{H}_{U}$ such that $U$ is a unitary, and $\Pi$ , $\widetilde{\Pi}$ are orthogonal projectors. Let $\gamma>1$ and $\delta,\epsilon\in(0,\frac{1}{2})$ . Suppose that $\widetilde{\Pi}U\Pi=W\Sigma V^{\dagger}=\sum_{i}\varsigma_{i}\ket{w_{i}}\bra{v_{i}}$ is a singular value decomposition. Then there is an $m=\mathcal{O}\Big{(}\frac{\gamma}{\delta}\log\left(\frac{\gamma}{\epsilon}\right)\Big{)}$ and an efficiently computable $\Phi\in\mathbb{R}^{m}$ such that

\left(\bra{+}\otimes\widetilde{\Pi}_{\leq\frac{1-\delta}{\gamma}}\right)U_{\Phi}\left(\ket{+}\otimes\Pi_{\leq\frac{1-\delta}{\gamma}}\right)=\sum_{i\colon\varsigma_{i}\leq\frac{1-\delta}{\gamma}}\tilde{\varsigma}_{i}\ket{w_{i}}\bra{v_{i}},\text{ where }\Big{|}\!\Big{|}\frac{\tilde{\varsigma}_{i}}{\gamma\varsigma_{i}}-1\Big{|}\!\Big{|}\leq\epsilon.

(1)

Moreover, $U_{\Phi}$ can be implemented using a single ancilla qubit with $m$ uses of $U$ and $U^{\dagger}$ , $m$ uses of C_ΠNOT and $m$ uses of C ${}_{\widetilde{\Pi}}$ NOT gates and $m$ single qubit gates. Here,

•

C_ΠNOT $:=X\otimes\Pi+I\otimes(I-\Pi)$ and a similar definition for C ${}_{\widetilde{\Pi}}$ NOT; see Definition 2 in gilyen2019quantum ,
•

$U_{\Phi}$ : alternating phase modulation sequence; see Definition 15 in gilyen2019quantum ,
•

$\Pi_{\leq\delta}$ , $\widetilde{\Pi}_{\leq\delta}$ : singular value threshold projectors; see Definition 24 in gilyen2019quantum .

III Main Framework

Review of Product Formula

Let the Hamiltonian of interest $H$ be $H=\sum_{i=1}^{m}H_{i}$ where each $H_{i}$ is a time-independent and that it can be efficiently simulated. For example, it can be done by using sparse access method berry2007efficient , quantum walk childs2010relationship ; berry2015hamiltonian , quantum signal processing low2017optimal , etc. The basic idea of product formula is to split the time evolution of $H$ into simpler evolutions for (sufficiently) small time steps. The most general theory has been developed in childs2021theory . Using recursively high-order Suzuki’s technique suzuki1991general , we have the following childs2021high ; berry2007efficient :

Lemma 6

Let $H=\sum_{i=1}^{m}H_{i}$ . Define $S_{2}(\lambda)=\prod_{i=1}^{m}e^{H_{i}\lambda/2}\prod_{i^{\prime}=m}^{1}e^{H_{i}^{\prime}\lambda/2}$ , and recursion relation:

S_{2k}(\lambda)=[S_{2k-2}(p_{k}\lambda]^{2}S_{2k-2}((1-4p_{k})\lambda)[S_{2k-2}(p_{k}\lambda]^{2}

with $p_{k}=(4-4^{1/(2k-1)})^{-1}$ for $k>1$ . Then we have that:

\displaystyle||\exp(-iHt)-S_{2k}(-it)||\leq 2\Gamma^{p+1}\frac{t^{p+1}}{p+1}\alpha^{(p)}

(2)

where $\Gamma=2\times 5^{k-1}$ and

\alpha^{(p)}=\sum_{i,j,...,2k+1=1}^{m}||[H_{2k+1}[H_{2k},...[H_{2},H_{1}]]]||

Furthermore, dividing the time interval $[0,t]$ into $r$ segments, we have:

\displaystyle||\exp(-iHt)-S_{2k}(-it/r)^{r}||\leq r||\exp(-iH\frac{t}{r}-S_{2k}(\frac{-it}{r})||\leq 2\Gamma^{2k+1}\frac{t^{p+1}}{r^{2k}(2k+1)}\alpha^{(2k)}

(3)

The above result suggests that, in order to have a total error to be $\delta$ , we need:

\displaystyle 2\Gamma^{2k+1}\frac{t^{p+1}}{r^{2k}(2k+1)}\alpha^{(2k)}=\delta

(4)

which implies that:

\displaystyle r=\mathcal{O}(\frac{1}{\delta^{1/2k}})

(5)

Thus, the complexity of simulating $H$ is (at best) sublinear in $1/\delta$ , as we need to use the simulation of each composing Hamiltonian $H_{i}$ at least $\mathcal{O}(r)$ times. Hence, the product formula itself possess a fundamental constraint on the (inverse of) error scaling, which seems inprohibitive. Now we proceed to describe a method for improving such scaling to polylogarithmic in $1/\delta$ .

Improved Framework

Let $\exp(-iH_{i}t)$ be the evolution operator being simulated at time $t$ . As it is a unitary operator, it block encodes itself (see further gilyen2019quantum and appendix 1 for reference). Our framework is based on the following result, which is Corollary 71 from gilyen2019quantum :

Lemma 7

Suppose that $U=exp(-iH)$ , where $H$ is a Hamiltonian of (spectral) norm $||H||\leq 1/2$ . Let $\epsilon\in(0,1/2]$ then we can implement a $(2/\pi,2,\epsilon)$ -block encoding of $H$ with $\mathcal{O}(\log(1/\epsilon)$ uses of controlled-U and its inverse, using $\mathcal{O}(\log(1/\epsilon)$ two-qubits gates and using a single ancilla qubit.

The above result, combined with the ability to simulate $\exp(-iH_{i})$ (basically simulating $H_{i}$ up to $t=1$ ), allows us to obtain the block encoding of $\pi H_{i}/2$ . Let $T_{i}$ be the complexity (could be query, or the gate complexity, i.e, the number of queries or gates required) of obtaining $\exp(-iH_{i})$ . Then the total complexity of obtaining the block encoding of $\pi H_{i}/2$ is $\mathcal{O}(T_{i}\log\frac{1}{\epsilon})$ (with extra 2 more qubits).

Given the block encodings of $\pi H_{i}/2$ (for $i=1,2,...,m$ ), the technique developed in Lemma 52 of gilyen2019quantum (see lemma 3 of appendix) allows us to use the block encoding of each $\pi H_{i}/2m$ to construct the block encoding of:

\displaystyle\frac{1}{m}\sum_{i=1}^{m}\frac{\pi}{2}H_{i}=\frac{\pi H}{2m}

(6)

Then again, one can use quantum singular value transformation technique gilyen2019quantum to simulate $\exp(-i\frac{\pi H}{2m}t)$ . More thoroughly, the central result of quantum singular value transformation framework is:

Lemma 8

[gilyen2019quantum Theorem 56] Suppose that $U$ is an $(\alpha,a,\epsilon)$ -encoding of a Hermitian matrix $A$ . (See Definition 43 of gilyen2019quantum for the definition.) If $P\in\mathbb{R}[x]$ is a degree- $d$ polynomial satisfying that

•

for all $x\in[-1,1]$ : $|P(x)|\leq\frac{1}{2}$ .

Then, there is a quantum circuit $\tilde{U}$ , which is an $(1,a+2,4d\sqrt{\frac{\epsilon}{\alpha}})$ -encoding of $P(A/\alpha)$ , and consists of $d$ applications of $U$ and $U^{\dagger}$ gates, a single application of controlled- $U$ and $\mathcal{O}((a+1)d)$ other one- and two-qubit gates.

To apply the above lemma, we make use of Jacobi-Anger expansion low2017optimal ; low2019hamiltonian ; childs2017lecture to approximate $\exp(-ixt)$ (for some real $x$ ) by some polynomials, i.e.,

\displaystyle\exp(-ixt)\approx J_{0}(-t)+2\sum_{k=1}^{K}i^{k}J_{k}(-t)T_{k}(x)

(7)

where $J_{k}$ is Bessel function and $T_{k}$ is Chebyshev polynomial. To obtain good approximation, for example, to make the approximation $\mathcal{O}(\epsilon)$ , the value of $K$ to be $\mathcal{O}(t+\log(\frac{1}{\epsilon})/\log(e+\ln(\frac{1}{\epsilon})/t)$ . Hence, these polynomials could be implemented using QSVT result gilyen2019quantum in a simple manner, as summarized in the above lemma.

Thus, it is straightforward to obtain the block encoding of $\exp(-i\frac{\pi H}{2m}t)$ for some fixed $t$ . One can then note that, by scaling $t\longrightarrow 2t/\pi$ , and repeating $m$ applications:

\displaystyle\prod_{1}^{m}\exp(-i\frac{\pi H}{2m}\frac{2t}{\pi})=\exp(-iHt)

(8)

Now we analyze the error resulting from approximating the block encoding of individual $H_{i}$ . Remind that using lemma 7 resulted in the block encoding of $\epsilon$ -approximated operator to $\pi H_{i}/2$ . Hence, the error accumulated when constructing combination $\frac{1}{m}\sum_{i=1}^{m}\frac{\pi H_{i}}{2}$ is $m\epsilon$ . Hence, according to lemma 8, the error resulted in for the approximation of block encoding of $\exp(-i\pi Ht/2m)$ is $\mathcal{O}(4(t+\log(\frac{1}{\epsilon}))\sqrt{m\epsilon})$ . The error further accumulated in the products: $\prod_{1}^{m}\exp(-i\frac{\pi H}{2m}\frac{2t}{\pi})$ , led to the error of approximating $\exp(-iHt)$ is $\mathcal{O}(4(t+\log(\frac{1}{\epsilon}))\sqrt{m\epsilon}m)$ . To make this error $\delta$ , we need to set:

	$\displaystyle 4(t+\log\frac{1}{\epsilon})\sqrt{m\epsilon}m=\delta$		(9)
	$\displaystyle\longrightarrow\epsilon=\frac{\delta^{2}}{16(t+\log\frac{1}{\epsilon})^{2}m^{3}}$		(10)

Given such $\epsilon$ , the total complexity for simulating $\exp(-iHt)$ up to error $\delta$ is then:

\mathcal{O}(m^{2}(\sum_{i=1}^{m}T_{i})\log(\frac{16(t+\log(1/\epsilon))^{2}m^{3}}{\delta^{2}})(t+\log(\log(\frac{1}{\delta^{2}}))))

which is logarithmical dependence on inverse of error. We remark that $T_{i}$ is the complexity of simulating $H_{i}$ for time $t=1$ , for which optimal result have been founded low2017optimal , which is

\mathcal{O}(d_{i}||H_{i}||+\log\frac{1}{\epsilon})

where $d_{i}$ is the sparsity of $H_{i}$ and $||H_{i}||$ is the spectral norm. As we have assumed the ability to simulate such constituent Hamiltonian, defining $d_{\max}=\max_{i}\{d_{i}\}$ and $||H||_{\max}=\max_{i}\{||H_{i}||\}$ , then we have that $\sum_{i}T_{i}=\mathcal{O}(m(d_{\max}||H||_{\max}+\log(1/\delta))$ . Hence, the total complexity is:

\mathcal{O}(m^{3}(d_{\max}||H||_{\max}+\log\frac{1}{\delta})\log(\frac{16(t+\log(1/\delta))^{2}m^{3}}{\delta^{2}})(t+\log(\log(\frac{1}{\delta^{2}}))))

Comparing to the standard product formula approach, our framework achieves polylogarithmcal dependence on inverse of error tolerance, which is superpolynomial improvement.

IV Conclusion

We have provided a simple method for constructing evolution operator of some Hamiltonian that could be written as summation of efficiently simulable Hamiltonian. Our method is motivated by the key constraint possessed by the standard product formula, e.g., Lie-Trotter-Suzuki formula, which prohibits logarithmical scaling on inverse of error tolerance. Thanks to the recently introduced quantum singular value transformation framework, we are able to use the simulation ability of composing Hamiltonian in an indirect manner, to construct the desired Hamiltonian directly. Then it can be efficiently simulated again, using technique from quantum singular value transformation framework.

Acknowledgement

We acknowledge support from Center for Distributed Quantum Processing, Stony Brook University.

References

(1) Richard P Feynman. Simulating physics with computers. In Feynman and computation, pages 133–153. CRC Press, 2018.
(2) Seth Lloyd. Universal quantum simulators. Science, 273(5278):1073–1078, 1996.
(3) Dorit Aharonov and Amnon Ta-Shma. Adiabatic quantum state generation and statistical zero knowledge. In Proceedings of the thirty-fifth annual ACM symposium on Theory of computing, pages 20–29, 2003.
(4) Dominic W Berry, Graeme Ahokas, Richard Cleve, and Barry C Sanders. Efficient quantum algorithms for simulating sparse hamiltonians. Communications in Mathematical Physics, 270(2):359–371, 2007.
(5) Dominic W Berry. High-order quantum algorithm for solving linear differential equations. Journal of Physics A: Mathematical and Theoretical, 47(10):105301, 2014.
(6) Andrew M Childs, Yuan Su, Minh C Tran, Nathan Wiebe, and Shuchen Zhu. Theory of trotter error with commutator scaling. Physical Review X, 11(1):011020, 2021.
(7) Andrew M Childs. On the relationship between continuous-and discrete-time quantum walk. Communications in Mathematical Physics, 294(2):581–603, 2010.
(8) Dominic W Berry, Andrew M Childs, and Robin Kothari. Hamiltonian simulation with nearly optimal dependence on all parameters. In 2015 IEEE 56th annual symposium on foundations of computer science, pages 792–809. IEEE, 2015.
(9) Andris Ambainis. Quantum walks and their algorithmic applications. International Journal of Quantum Information, 1(04):507–518, 2003.
(10) Andris Ambainis. Quantum walk algorithm for element distinctness. SIAM Journal on Computing, 37(1):210–239, 2007.
(11) Andrew M Childs and Nathan Wiebe. Hamiltonian simulation using linear combinations of unitary operations. arXiv preprint arXiv:1202.5822, 2012.
(12) Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D Somma. Simulating hamiltonian dynamics with a truncated taylor series. Physical review letters, 114(9):090502, 2015.
(13) Guang Hao Low and Isaac L Chuang. Optimal hamiltonian simulation by quantum signal processing. Physical review letters, 118(1):010501, 2017.
(14) Guang Hao Low and Isaac L Chuang. Hamiltonian simulation by qubitization. Quantum, 3:163, 2019.
(15) David Poulin, Angie Qarry, Rolando Somma, and Frank Verstraete. Quantum simulation of time-dependent hamiltonians¡? format?¿ and the convenient illusion of hilbert space. Physical review letters, 106(17):170501, 2011.
(16) Dominic W Berry, Andrew M Childs, Yuan Su, Xin Wang, and Nathan Wiebe. Time-dependent hamiltonian simulation with $l1$ -norm scaling. Quantum, 4:254, 2020.
(17) Yi-Hsiang Chen, Amir Kalev, and Itay Hen. Quantum algorithm for time-dependent hamiltonian simulation by permutation expansion. PRX Quantum, 2(3):030342, 2021.
(18) Dong An, Di Fang, and Lin Lin. Time-dependent hamiltonian simulation of highly oscillatory dynamics and superconvergence for schrödinger equation. Quantum, 6:690, 2022.
(19) Guang Hao Low and Nathan Wiebe. Hamiltonian simulation in the interaction picture. arXiv preprint arXiv:1805.00675, 2018.
(20) Mária Kieferová, Artur Scherer, and Dominic W Berry. Simulating the dynamics of time-dependent hamiltonians with a truncated dyson series. Physical Review A, 99(4):042314, 2019.
(21) Masuo Suzuki. General theory of fractal path integrals with applications to many-body theories and statistical physics. Journal of mathematical physics, 32(2):400–407, 1991.
(22) András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe. Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics. In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204, 2019.
(23) Arthur G Rattew and Patrick Rebentrost. Non-linear transformations of quantum amplitudes: Exponential improvement, generalization, and applications. arXiv preprint arXiv:2309.09839, 2023.
(24) Daan Camps and Roel Van Beeumen. Approximate quantum circuit synthesis using block encodings. Physical Review A, 102(5):052411, 2020.
(25) Andrew M Childs, Jin-Peng Liu, and Aaron Ostrander. High-precision quantum algorithms for partial differential equations. Quantum, 5:574, 2021.
(26) Andrew M Childs. Lecture notes on quantum algorithms. Lecture notes at University of Maryland, 2017.