A quantum algorithm for simulating non-sparse Hamiltonians

TL;DR

This paper introduces a quantum algorithm capable of simulating non-sparse Hamiltonians efficiently by leveraging qRAM and quantum walks, achieving poly-logarithmic precision dependence and broad applicability.

Contribution

The paper presents a novel quantum algorithm for simulating general Hamiltonians using qRAM and quantum walks, extending beyond sparse cases with improved efficiency.

Findings

Achieves poly-logarithmic dependence on precision.

Provides a complexity of $O(t ext{sqrt}(N) ext{polylog}(N, t ext{||}H ext{||}, 1/ ext{epsilon}))$.

Applicable as a subroutine for unitary implementation and quantum linear systems.

Abstract

We present a quantum algorithm for simulating the dynamics of Hamiltonians that are not necessarily sparse. Our algorithm is based on the input model where the entries of the Hamiltonian are stored in a data structure in a quantum random access memory (qRAM) which allows for the efficient preparation of states that encode the rows of the Hamiltonian. We use a linear combination of quantum walks to achieve poly-logarithmic dependence on precision. The time complexity of our algorithm, measured in terms of the circuit depth, is $O(t\sqrt{N}\|H\|\,\mathrm{polylog}(N, t\|H\|, 1/\epsilon))$, where $t$ is the evolution time, $N$ is the dimension of the system, and $\epsilon$ is the error in the final state, which we call precision. Our algorithm can be directly applied as a subroutine for unitary implementation and quantum linear systems solvers, achieving $\widetilde{O}(\sqrt{N})$ dependence for both applications.

Figures

Tables

Table 1. Table 1: Comparing our result O ​ ( t ​ d ​ ∥ H ∥ ​ polylog ​ ( t , d , ∥ H ∥ , 1 / ϵ ) ) 𝑂 𝑡 𝑑 delimited-∥∥ 𝐻 polylog 𝑡 𝑑 delimited-∥∥ 𝐻 1 italic-ϵ O(t\sqrt{d}\lVert H\rVert\,\mathrm{polylog}(t,d,\lVert H\rVert,1/\epsilon)) with other quantum and classical algorithms for different models. Since the qRAM model is stronger than the sparse-access model and the classical sampling and query access model, we consider the advantage of our algorithm against others when they are directly applied to the qRAM model.

Model	State-of-the-art	Advantage of our results
Sparse-access with ${\Vert H\Vert }_{\max }$ dependence	$\tilde{O}\left(td{\Vert H\Vert }_{\max }\right)$ [9]	No advantage
Sparse-access with $\parallel H\parallel$ dependence	$O\left({\left(t\sqrt{d}\parallel H\parallel \right)}^{1+o\left(1\right)}/{ϵ}^{o\left(1\right)}\right)$ [30]	Subpolynomial improvement in $t,d$; exponential improvement in $ϵ$
qRAM	$\tilde{O}\left(t\sqrt{d}\parallel H\parallel \right)$ [12]	Same result
Classical sampling and query access	$\mathrm{poly}(t,{\parallel H\parallel }_F,1/ϵ)$ [13]	Polynomial speedup

Equations

$$O_S|i,j⟩|z⟩$$

$$O_H|i,j⟩|z⟩$$

$$O\left(t{\Vert H\Vert }_1{n}^2{\log }^{5/2}(t{\Vert H\Vert }_1/ϵ)\frac{\log (t\Vert H\Vert /ϵ)}{\log \log (t\Vert H\Vert /ϵ)}\right).$$

$$O\left(t\sqrt{N}\Vert H\Vert n^2{\log }^{5/2}(t\sqrt{N}\Vert H\Vert /ϵ)\frac{\log (t\Vert H\Vert /ϵ)}{\log \log (t\Vert H\Vert /ϵ)}\right).$$

$$O\left(t\sqrt{d}\Vert H\Vert n^2{\log }^{5/2}(t\sqrt{d}\Vert H\Vert /ϵ)\frac{\log (t\Vert H\Vert /ϵ)}{\log \log (t\Vert H\Vert /ϵ)}\right).$$

$$H=\left(\begin{array}{cc}0& U\\ U^{†}& 0\end{array}\right),$$

$$|j⟩|0^{\log N}⟩|0⟩\mapsto \frac{1}{\sqrt{{\Vert H\Vert }_1}}|j⟩\sum _{k=0}^{N-1}|k⟩\left(\sqrt{H_{jk}^{\ast }}|0⟩+\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{N}}|1⟩\right),$$

$$U:|\widetilde{{\varphi }_0}⟩|\widetilde{{\varphi }_1}⟩|\tilde{\theta }⟩|0⟩\mapsto |\widetilde{{\varphi }_0}⟩|\widetilde{{\varphi }_1}⟩|\tilde{\theta }⟩\left(e^{i\widetilde{{\varphi }_0}}\cos (\tilde{\theta })|0⟩+e^{i\widetilde{{\varphi }_1}}\sin (\tilde{\theta })|1⟩\right).$$

$$U=\left(\sum _{\widetilde{{\varphi }_0}\in \{0,1{\}}^b}|\widetilde{{\varphi }_0}⟩⟨\widetilde{{\varphi }_0}|\otimes e^{i|0⟩⟨0|\widetilde{{\varphi }_0}}\right)\left(\sum _{\widetilde{{\varphi }_1}\in \{0,1{\}}^b}|\widetilde{{\varphi }_1}⟩⟨\widetilde{{\varphi }_1}|\otimes e^{i|1⟩⟨1|\widetilde{{\varphi }_1}}\right)\left(\sum _{\tilde{\theta }\in \{0,1{\}}^b}|\tilde{\theta }⟩⟨\tilde{\theta }|\otimes e^{-iY\tilde{\theta }}\right),$$

$$\frac{1}{\sqrt{{\Vert H\Vert }_1}}|0⟩|0^{\log N}⟩\left(\sqrt{\sum _{k=0}^{N-1}|H_{jk}^{\ast }|}|0⟩+\sqrt{{\Vert H\Vert }_1-{\sigma }_j}|1⟩\right).$$

$$\frac{1}{\sqrt{{\Vert H\Vert }_1}}|00⟩\left(\sqrt{|c_0|+|c_1|+|c_2|+|c_3|}|0⟩+\sqrt{{\Vert H\Vert }_1-{\sigma }_j}|1⟩\right)=$$

$$\frac{1}{\sqrt{{\Vert H\Vert }_1}}\left(\sqrt{|c_0|+|c_1|+|c_2|+|c_3|}|00⟩|0⟩+\sqrt{{\Vert H\Vert }_1-{\sigma }_j}|00⟩|1⟩\right).$$

$$\frac{1}{\sqrt{{\Vert H\Vert }_1}}\left(\left(\sqrt{|c_0|+|c_1|}|00⟩+\sqrt{|c_2|+|c_3|}|10⟩\right)|0⟩+\left(\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{2}}(|00⟩+|10⟩)\right)|1⟩\right).$$

$$\frac{1}{\sqrt{{\Vert H\Vert }_1}}(\sqrt{c_0}|00⟩|0⟩+\sqrt{c_1}|01⟩|0⟩+\sqrt{c_2}|10⟩|0⟩+\sqrt{c_3}|11⟩|0⟩+$$

$$\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{4}}|00⟩|1⟩+\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{4}}|01⟩|1⟩+\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{4}}|10⟩|1⟩+\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{4}}|11⟩|1⟩)=$$

$$\frac{1}{\sqrt{{\Vert H\Vert }_1}}(|00⟩\left(\sqrt{c_0}|0⟩+\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{4}}|1⟩\right)+|01⟩\left(\sqrt{c_1}|0⟩+\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{4}}|1⟩\right)+$$

$$|10⟩\left(\sqrt{c_2}|0⟩+\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{4}}|1⟩\right)+|11⟩\left(\sqrt{c_3}|0⟩+\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{4}}|1⟩\right)).$$

$$T=\sum _{j=0}^{N-1}\sum _{b\in \{0,1\}}(|j⟩⟨j|\otimes |b⟩⟨b|)\otimes |{\phi }_{jb}⟩,$$

$$|{\phi }_{j0}⟩=\frac{1}{\sqrt{{\Vert H\Vert }_1}}\sum _{k=0}^{N-1}|k⟩\left(\sqrt{H_{jk}^{\ast }}|0⟩+\sqrt{\frac{{\Vert H\Vert }_1-{\sigma }_j}{N}}|1⟩\right),$$

$$⟨j|⟨0|T^{†}ST|k⟩|0⟩=\frac{\sqrt{H_{jk}}{\left(\sqrt{H_{jk}^{\ast }}\right)}^{\ast }}{{\Vert H\Vert }_1}=\frac{H_{jk}}{{\Vert H\Vert }_1},$$

$$(I\otimes ⟨0|)T^{†}ST(I\otimes |0⟩)=\frac{H}{{\Vert H\Vert }_1}.$$

$$U=iS(2T{T}^{†}-I).$$

$$U|{\mu }_{\pm }⟩={\mu }_{\pm }|{\mu }_{\pm }⟩,$$

$$|{\mu }_{\pm }⟩=$$

$${\mu }_{\pm }=$$

$$U|{\mu }_{\pm }⟩={\mu }_{\pm }T|\lambda ⟩|0⟩+i\left(1+\frac{2\lambda i}{{\Vert H\Vert }_1}{\mu }_{\pm }\right)ST|\lambda ⟩|0⟩.$$

$$1+\frac{2\lambda i}{{\Vert H\Vert }_1}{\mu }_{\pm }={\mu }_{\pm }^2,$$

$${\mu }_{\pm }=\frac{\lambda i}{{\Vert H\Vert }_1}\pm \sqrt{1-\frac{{\lambda }^2}{{\Vert H\Vert }_1^2}}=\pm e^{\pm i\arcsin (\lambda /{\Vert H\Vert }_1)}.$$

$$\sum _{m=-\infty }^{\infty }{J}_m(z){\mu }_{\pm }^m=\exp \left(\frac{z}{2}\left({\mu }_{\pm }-\frac{1}{{\mu }_{\pm }}\right)\right)=e^{iz\lambda /{\Vert H\Vert }_1},$$

$$V_{\infty }=\sum _{m=-\infty }^{\infty }\frac{J_m(z)}{{\sum }_{j=-\infty }^{\infty }{J}_j(z)}{U}^m=\sum _{m=-\infty }^{\infty }{J}_m(z)U^m=e^{izH/{\Vert H\Vert }_1},$$

$$V_k=\sum _{m=-k}^k\frac{J_m(z)}{{\sum }_{j=-k}^k{J}_j(z)}{U}^m.$$

$$\sum _{m=-k}^k\frac{J_m(z)}{{\sum }_{j=-k}^k{J}_j(z)}{\mu }_{\pm }^m.$$

$$k=O\left(\frac{\log (\Vert H\Vert /({\Vert H\Vert }_1ϵ))}{\log \log (\Vert H\Vert /({\Vert H\Vert }_1ϵ))}\right),$$

Full text

1 Introduction

1.1 Background and main results

Hamiltonian simulation is the problem of simulating the dynamics of quantum systems, which is the original motivation for quantum computers [19, 20]. It has been shown to be BQP-hard and is hence conjectured not to be classically solvable in polynomial time, since such an algorithm would imply an efficient classical algorithm for any problem with an efficient quantum algorithm, including integer factorization [38].

Different input models have been considered in previous quantum algorithms for simulating Hamiltonian evolution. The local Hamiltonian model is specified by the local terms of a Hamiltonian. The sparse-access model for a sparse Hamiltonian $H$ is specified by the following two oracles:

[TABLE]

where $S_{i,j}$ is the $j$-th nonzero entry of the $i$-th row and $\oplus$ denotes the bit-wise XOR. The linear combination of unitaries (LCU) model is specified by a decomposition of a Hamiltonian as a linear combination of unitaries and we are given the coefficients and access to each unitary. Following the first proposal by Lloyd [28] for local Hamiltonians, Aharonov and Ta-Shma gave an efficient algorithm for sparse Hamiltonians [4]. Subsequently, many algorithms have been proposed which improved the runtime [6, 7, 8, 14, 17, 35, 40, 33, 10, 31], mostly in the sparse-access model, and have recently culminated in a few works with optimal dependence on all (or nearly all) parameters for sparse Hamiltonians [9, 32, 30].

While the above-mentioned input models arise naturally in many physics applications and matrix arithmetic applications (i.e., we have access to the local terms of a Hamiltonian or each entry of a Hamiltonian can be efficiently computed), in many machine learning applications, it is more convenient to work with a different input model, namely, the quantum random access memory (qRAM) model, where we assume that the entries of a Hamiltonian are stored in a data structure as in [26] and we have quantum access to the memory. As we receive and process the input data, with little overhead, they can be stored in the data structure, and then the qRAM will facilitate preparing quantum states corresponding to the input data. The use of the qRAM model has been successfully demonstrated in many applications such as quantum principal component analysis [29], quantum support vector machines [36], and quantum recommendation systems [26].

In this work, we consider the qRAM model for obtaining information about the Hamiltonian that is not necessarily sparse. Quantum access to this data structure allows us to efficiently prepare states that encode the rows of the Hamiltonian. Using the ability to prepare these states in combination with a quantum walk [9], we give the first (to the best of our knowledge) Hamiltonian simulation algorithm in the qRAM model whose time complexity has $\widetilde{O}(\sqrt{N})$ dependence111In this paper, we use $\widetilde{O}(\cdot)$ to hide poly-logarithmic factors. for non-sparse Hamiltonians of dimensionality $N$. As a subroutine of quantum linear systems solver given by [16], our result directly implies a quantum linear systems solver in the qRAM model with square-root dependence on dimension and poly-logarithmic dependence on precision, which exponentially improves the precision dependence of the quantum linear systems solver by [41] with the same input model. Since solving linear systems is a fundamental procedure of many machine learning tasks, our algorithm has extensive potential applications in quantum machine learning.

In [26], a quantum algorithm for recommendation systems was introduced based on an explicit description of a data structure, which resulted in a fast quantum algorithm for estimating singular values for any real matrix. This fast singular value estimation algorithm was used in a quantum algorithm for solving dense linear systems [41]. The data structure in [26, 41] allows us to prepare states that correspond to row vectors. In our algorithm for the Hamiltonian simulation problem, the main hurdle is to efficiently prepare the states which allow for a quantum walk corresponding to $e^{-iH/\lVert H\rVert_{1}}$. These states are quite different from those in [26, 7, 9]: the states required by [9] allow for a quantum walk corresponding to $e^{-iH/(d\left\|H\right\|_{\mathrm{max}})}$, where $d$ is the row-sparsity of $H$, and their states can be prepared with $O(1)$ queries to the sparse-access oracle. However, the states required by our algorithm (See Eq. (7)) is less structural and it is not known how to prepare them with poly-logarithmic cost in the sparse-access model. In the qRAM model, we assume the entries of $H$ are stored in a data structures, which permits state preparation with time complexity (circuit depth) $O(\mathrm{polylog}(N))$. The precise definition of the data structure is presented in Definition 1.

Using the efficient state preparation procedure, we implement the linear combination of quantum walks in order to simulate the time evolution for non-sparse Hamiltonians with only poly-logarithmic dependence on precision. The main result of this work is summarized in the following theorem, which we prove in Sec. 3.

Theorem 1** (Non-sparse Hamiltonian Simulation).**

Let $H\in\mathbb{C}^{N\times N}$ (with $N=2^{n}$) be a Hermitian matrix stored in the data structure as specified in Definition 1. There exists a quantum algorithm for simulating the evolution of $H$ for time $t$ and error $\epsilon$ with time complexity (circuit depth)

[TABLE]

In this work, we use the notation $\left\|\cdot\right\|_{1}$ to denote the induced 1-norm (i.e., maximum absolute row-sum norm), defined as $\left\|H\right\|_{1}=\max_{j}\sum_{k=0}^{N-1}|H_{jk}|$; we use $\lVert\cdot\rVert$ to denote the spectral norm, and use $\left\|\cdot\right\|_{\mathrm{max}}$ to denote the max norm, defined as $\left\|H\right\|_{\mathrm{max}}=\max_{i,j}|H_{jk}|$.

By the fact that $\left\|H\right\|_{1}\leq\sqrt{N}\lVert H\rVert$ (see [15], alternatively, a more generalized version of this relation is shown in Appendix A), we immediately have the following corollary.

Corollary 1**.**

Let $H\in\mathbb{C}^{N\times N}$ (where $N=2^{n}$) be a Hermitian matrix stored in the data structure as specified in Definition. 1. There exists a quantum algorithm for simulating the evolution of $H$ for time $t$ and error $\epsilon$ with time complexity (circuit depth)

[TABLE]

Remarks:

In Theorem 1, the circuit depth scales as $\widetilde{O}(\sqrt{N})$. However, the gate complexity could be as large as $O(N^{2.5}\log^{2}(N))$ in general because of the addressing scheme which allows for quantum access to classical data stored in the data structure as specified in Definition 1. If some structure of $H$ is promised (e.g., the entries of $H$ repeat in some pattern), the addressing scheme could be implemented efficiently. 2. 2.

If $H$ is $d$-sparse (i.e., $H$ has at most $d$ non-zero entries in each row), the time complexity (circuit depth) of our algorithm is

[TABLE]

This follows from Theorem 1 and the fact that $\left\|H\right\|_{1}\leq\sqrt{d}\lVert H\rVert$ (as shown in Appendix A). 3. 3.

The techniques of our algorithm also work for the sparse-access model where we use standard state preparation techniques to prepare the state as in Eq. (7). Compared to the qRAM model, the time complexity of the state preparation for Eq. (7) in the sparse-access model incurs an additional $O(d)$ factor (for computing $\sigma_{j}$). Now, to simulate $H$ for time $t$ in the sparse-access model, the dependence of the time complexity on $t$, $d$, and $\lVert H\rVert$ becomes $O(td^{1.5}\lVert H\rVert)$ as $\lVert H\rVert_{1}\leq\sqrt{d}\lVert H\rVert$. Hence our techniques have no advantage over previous results in the sparse-access model.

1.2 Related work

Hamiltonian simulation with $\lVert H\rVert$ dependence.

For non-sparse Hamiltonians, a suitable model is the black-box model (a special case of the sparse-access model): querying the oracle with an index-pair $\ket{i,j}$ returns the corresponding entry of the Hamiltonian $H$, i.e., $O_{H}$ defined in Eq. (2). With access to a black-box Hamiltonian, simulation with error $\epsilon$ can then be provably performed with query complexity $O((\lVert H\rVert t)^{3/2}N^{3/4}/\sqrt{\epsilon})$ for dense Hamiltonians [7], and it was also empirically observed in [7] that for several classes of Hamiltonians, $O(\sqrt{N}\log(N))$ queries suffice. Whether this $\widetilde{O}(\sqrt{N})$ dependence holds for all Hamiltonians was left as an open problem. After the first version of this work was made public, this open problem was almost resolved by Low [30], where he proposed a quantum algorithm for simulating black-box Hamiltonians with time complexity $O((t\sqrt{N}\lVert H\rVert)^{1+o(1)}/\epsilon^{o(1)})$. Although our input model is stronger than the black-box model, our work distinguishes itself since the two models are still comparable in many quantum machine learning applications and our work gives a better complexity.

Hamiltonian simulation with $\left\|H\right\|_{\mathrm{max}}$ dependence.

Since the qRAM model is stronger than the black-box model and the sparse-access model, previous quantum algorithms such as [7, 9] can be directly used to simulate Hamiltonians in the qRAM model. For a $d$-sparse Hamiltonian, the circuit depth of the black-box Hamiltonian simulation is given by $\widetilde{O}(td\left\|H\right\|_{\mathrm{max}})$ as shown in [7, 9]. When $H$ is non-sparse, their results imply the scaling $\widetilde{O}(tN\left\|H\right\|_{\mathrm{max}})$. In applications where $\left\|H\right\|_{\mathrm{max}}$ is a measure of cost, our result has no advantage against theirs, as the inequality $\lVert H\rVert\leq\sqrt{N}\left\|H\right\|_{\mathrm{max}}$ implies $\widetilde{O}(t\sqrt{H}\lVert H\rVert)=\widetilde{O}(tN\left\|H\right\|_{\mathrm{max}})$. However, in the case where $\lVert H\rVert$ is a measure of cost, such as solving linear systems [23, 16], and black-box unitary implementation [7], our result has a quadratic improvement in the dimensionality dependence, as the inequality $\left\|H\right\|_{\mathrm{max}}\leq\lVert H\rVert$ implies $\widetilde{O}(tN\left\|H\right\|_{\mathrm{max}})=\widetilde{O}(tN\lVert H\rVert)$.

Hamiltonian simulation in the qRAM model.

Shortly after the first version of this paper was made public, Chakraborty, Gilyén, and Jeffery [12] independently proposed a quantum algorithm for simulating non-sparse Hamiltonians based on the similar qRAM input model which achieved the same time complexity as our method. Their work is based on a very general input model, namely, the block-encoding model, which was originally proposed in [33]. The block-encoding model assumes we are given a unitary \bigl{(}\begin{smallmatrix}H/\alpha&\cdot\\ \cdot&\cdot\end{smallmatrix}\bigr{)} that contains $H/\alpha$ in its upper-left block. Then the evolution $e^{-iHt}$ can be simulated in $\widetilde{O}(\alpha\lVert H\rVert t)$ time. It was shown in [33] that working with the sparse-access model of a $d$-sparse Hamiltonian $H$, a block-encoding of $H$ with $\alpha=d$ can be efficiently implemented, and it hence implies a simulation algorithm with time complexity $\widetilde{O}(d\lVert H\rVert t)$. One of the main results in [12] is that working with the qRAM model of a $d$-sparse Hamiltonian $H$, a block-encoding with $\alpha=\sqrt{d}$ can be efficiently implemented, which yields a simulation algorithm with time complexity $\widetilde{O}(\sqrt{d}\lVert H\rVert t)$. This result has the same complexity as ours, but their techniques are more general. In fact, the techniques of [12] has been generalized in [21] to a quantum framework for implementing singular value transformation of matrices. Another merit of [12] is that they gave detailed analysis for application to a quantum linear systems solvers in the qRAM model.

Quantum-inspired classical algorithms for Hamiltonian simulation.

Recently (after the first version of this paper was made public), several classical algorithms [37, 13] have been proposed for simulating Hamiltonians in a classical input model, namely, the sampling and query access, which is comparable to the qRAM model we use. The algorithm in [37] is efficient in the low-rank and sparse regime, and the algorithm in [13] is efficient when $H$ is low-rank: their time complexity scales as $\mathrm{poly}(t,\lVert H\rVert_{F},1/\epsilon)$, where $\lVert H\rVert_{F}$ is the Frobenius norm of $H$. While these results have ruled out the possibility of exponential speedups of our algorithm in the low-rank regime, noting that the degrees in the polynomials of the time complexity of these classical algorithms are large, our work still has a polynomial speedup over classical algorithm (for low-rank Hamiltonians and dense Hamiltonians).

To summarize this subsection, we provide Table 1 for state-of-the-art algorithms for Hamiltonian simulation in different models and how our techniques compare with theirs in the qRAM model.

1.3 Applications

Unitary implementation.

One immediate application of simulating non-sparse Hamiltonians is the unitary implementation problem: given access to the entries of a unitary $U$, the objective is to construct a quantum circuit to approximate $U$ with precision $\epsilon$. As proposed in [7, 25], unitary implementation can be reduced to Hamiltonian simulation by considering the Hamiltonian

[TABLE]

and the fact that $e^{-iH\pi/2}\ket{1}\ket{\psi}=-i\ket{0}U\ket{\psi}$. If the entries of a unitary matrix can be accessed by a black-box query oracle, it is shown in [7] that this unitary operator can be implemented with $O(N^{2/3}\mathrm{polylog}(N)\mathrm{poly}(1/\epsilon))$ queries to the black box. Assume now that the entries of $U$ are stored in a data structure such as in Definition 1. Using the same reduction to Hamiltonian simulation, our Hamiltonian simulating algorithm implies an implementation of $U$ with time complexity (circuit depth) $O(\sqrt{N}\,\mathrm{polylog}(N,1/\epsilon))$.

Quantum linear systems solver.

Simulating non-sparse Hamiltonians can also be used as a subroutine for solving linear systems of equations for non-sparse matrices. The essence for solving a linear system $A\ket{x}=\ket{b}$ is to apply $A^{-1}$ on $\ket{b}$, assuming here for simplicity that $\ket{b}$ is entirely in the column-space of $A$. When $A$ is $d$-sparse, it is shown in [16] that $A^{-1}$ can be approximated as a linear combination of unitaries of the form $e^{-iAt}$. An efficient quantum algorithm for Hamiltonian simulation such as [9] can then be used as a subroutine, so that this linear system can be solved with gate complexity $O(d\kappa^{2}\mathrm{polylog}(N,\kappa/\epsilon))$. If $A$ is non-sparse, the time complexity of their algorithm scales as $\widetilde{O}(N)$. Based on the data structure in [26], a quantum algorithm for solving linear systems for non-sparse matrices was described in [41], with time complexity (circuit depth) $O(\kappa^{2}\sqrt{N}\mathrm{polylog}(N)/\epsilon)$. If we assume a similar input model as in [26, 41], using our Hamiltonian simulation algorithm, together with the linear combinations of unitaries (LCU) decompositions in [16], we have a quantum algorithm for solving linear systems for non-sparse matrices with time complexity (circuit depth) $O(\kappa^{2}\sqrt{N}\mathrm{polylog}(\kappa/\epsilon))$, which is an exponential improvement in error dependence compared to [41].

In the remainder of this paper, we first define our data structure, then describe the algorithm in detail and finally prove the main results. We then finish with a summary of this work and a discussion of the benefits and possible drawbacks of our algorithm.

2 Data structure and quantum walk

We first give the precise definition of the data structure that stores the entries of the Hamiltonian, and then show how this data structure can be used to prepare states that will allow us to perform fast Hamiltonian simulation even for dense matrices. The data structure was introduced in in [26].

Definition 1** (Data Structure).**

Let $H\in\mathbb{C}^{N\times N}$ be a Hermitian matrix (where $N=2^{n}$), $\left\|H\right\|_{1}$ being the maximum absolute row-sum norm, and $\sigma_{j}:=\sum_{k}\lvert H_{jk}\rvert$. Each entry $H_{jk}$ is represented with $b$ bits of precision. Define $D$ as an array of $N$ binary trees $D_{j}$ for $j\in\{0,\ldots,N-1\}$. Each $D_{j}$ corresponds to the row $H_{j}$, and its organization is specified by the following rules.

The leaf node $k$ of the tree $D_{j}$ stores the value222Note that the conjugation here is necessary. See Eq. (16).* ${H_{jk}^{*}}$ corresponding to the index-entry pair $(j,k,H_{jk})$.* 2. 2.

For the level immediately above the bottom level, i.e., the leaves, and any node level above the leaves, the data stored is determined as follows: suppose the node has two children storing data $a$, and $b$ respectively (note that $a$ and $b$ are complex numbers). Then the entry that is stored in this node is given by $(|a|+|b|)$.

An example of the above data structure is shown Fig. 1. Note that for each binary tree $D_{j}$ in the data structure, the value stored in an internal (non-leaf) node is a real number, while for a leaf node, the value stored is a complex number. The root node of $D_{j}$ stores the value $\sum_{k=0}^{N-1}|H_{jk}^{*}|$ and we can calculate the value $\left\|H\right\|_{1}-\sigma_{j}$ in constant time. In addition, $\left\|H\right\|_{1}$ can be computed by taking the maximum value of the roots of all the binary trees, which can be done during the construction of the data structure.

Using this data structure, we can efficiently perform the mapping described in the following technical lemma for efficient state preparation.

Lemma 1** (State Preparation).**

Let $H\in\mathbb{C}^{N\times N}$ be a Hermitian matrix (where $N=2^{n}$) stored in the data structure as specified in Definition 1. Each entry $H_{jk}$ is represented with $b$ bits of precision. Then the following holds

Let $\left\|H\right\|_{1}$ be defined as above. A quantum computer that has access to the data structure can perform the following mapping for $j\in\{0,\ldots,N-1\}$,

[TABLE]

with time complexity (circuit depth) $O(n^{2}b^{5/2}\log b)$, where $\sigma_{j}=\sum_{k}|H_{jk}|$, and the square-root satisfies \sqrt{H_{jk}}\bigl{(}\sqrt{H_{jk}^{*}}\bigr{)}^{*}=H_{jk}. 2. 2.

The size of the data structure containing all $N^{2}$ complex entries is $O(N^{2}\log^{2}(N))$.

In order to perform the mapping we will need the following Lemma, which enables us to efficiently implement the conditional rotations with complex numbers.

Lemma 2**.**

Let $\theta,\phi_{0},\phi_{1}\in\mathbb{R}$ and let $\widetilde{\theta},\widetilde{\phi_{0}},\widetilde{\phi_{1}}$ be the $b$-bit finite precision representation of $\theta,\phi_{0}$, and $\phi_{1}$, respectively. Then there exists a unitary $U$ that performs the following mapping:

[TABLE]

Moreover, $U$ can be implemented with $O(b)$ 1- and 2-qubit gates.

Proof.

Define $U$ as

[TABLE]

where Y=\bigl{(}\begin{smallmatrix}0&-i\\ i&0\end{smallmatrix}\bigr{)} is the Pauli $Y$ matrix.

To implement the operator $\sum_{\widetilde{\theta}\in\{0,1\}^{b}}\ket{\widetilde{\theta}}\!\bra{\widetilde{\theta}}\otimes e^{-iY\widetilde{\theta}}$, we use one rotation controlled on each qubit of the first register, with the rotation angles halved for each successive bit. The other two factors of $U$ can be implemented in a similar way. Therefore, $U$ can be implemented with $O(b)$ 1- and 2-qubit gates. ∎

Before proving Lemma 1, we first describe the construction and the size of the data structure. Readers may refer to [26] for more details of this data structure.

• •

The data structure is built from $N$ binary trees $D_{i},i\in\{0,\dots,N-1\}$ and we start with an empty tree.

• •

When a new entry $(i,j,H_{ij})$ arrives, we create or update the leaf node $j$ in the tree $D_{i}$, where the adding of the entry takes $O(\log(N))$ time, since the depth of the tree for $H\in\mathbb{C}^{N\times N}$ is at most $\log(N)$. Since the path from the root to the leaf is of length at most $\log(N)$ (under the assumption that $N=2^{n}$), we have furthermore to update at most $\log(N)$ nodes, which can be done in $O(\log(N))$ time if we store an ordered list of the levels in the tree.

• •

The total time for updating the tree with a new entry is given by $\log(N)\times\log(N)=\log^{2}(N)$.

• •

The memory requirements for $k$ entries are given by $O(k\log^{2}(N))$ as for every entry $(j,k,H_{jk})$ at least $\log(N)$ nodes are added and each node requires at most $O(\log(N))$ bits.

Now we are ready to prove Lemma 1.

Proof of Lemma 1.

With this data structure, we can perform the mapping specified in Eq. (7), with the following steps. For each $j$, we start from the root of $D_{j}$. Starting with the initial state $\ket{j}\ket{0^{\log N}}\ket{0}$, first apply the rotation (according to the value stored in the root node and calculating the normalisation in one query) on the last register to obtain the state

[TABLE]

Then a sequence of conditional rotations is applied on each qubit of the second register to obtain the state as in Eq. (7). At level $\ell$ of the binary tree $D_{j}$, a query to the data structure is made to load the data $c$ (stored in the node) into a register in superposition, the rotation to perform is proportional to \bigl{(}\sqrt{c},\sqrt{(\left\|H\right\|_{1}-\sigma_{j})/2^{\ell}}\bigr{)} (assuming at the root, $\ell=0$, and for the leaves, $\ell=\log N$). Then the rotation angles will be determined by calculating the square root and trigonometric functions on the output of the query: this can be implemented with $O(b^{5/2})$ 1- and 2-qubit gates using simple techniques based on Taylor series and long multiplication as in [9], where the error is smaller than that caused by truncating to $b$ bits. Then the conditional rotation is applied by the circuit described in Lemma 2, and the cost for the conditional rotation is $O(b)$. There are $n=\log(N)$ levels, so the cost excluding the implementation of the oracle is $O(nb^{5/2})$. To obtain quantum access to the classical data structure, a quantum addressing scheme is required. One addressing scheme described in [22] can be used. Although the circuit size of this addressing scheme is $\widetilde{O}(N)$ for each $D_{j}$, its circuit depth is $O(n)$. Therefore, the time complexity (circuit depth) for preparing the state in Eq. (7) is $O(n^{2}b^{5/2}\log n)$.

We use the following rules to determine the sign of the square-root of a complex number: if $H_{jk}$ is not a negative real number, we write $H_{jk}=re^{i\varphi}$ (for $r\geq 0$ and $-\pi\leq\varphi\leq\pi$) and take $\sqrt{H_{jk}^{*}}=\sqrt{r}e^{-i\varphi/2}$; when $H_{jk}$ is a negative real number, we take $\sqrt{H_{jk}^{*}}=\mathrm{sign}(j-k)i\sqrt{|H_{jk}|}$ to avoid the sign ambiguity. With this recipe, we have \sqrt{H_{jk}}\bigl{(}\sqrt{H_{jk}^{*}}\bigr{)}^{*}=H_{jk}. ∎

In the following, we demonstration a state preparation procedure based on the data structure in Fig. 1. In this example, the initial state (omitting the first register) is $\ket{00}\ket{0}$. Let $\sigma_{j}=|c_{0}|+|c_{1}|+|c_{2}|+|c_{3}|$. Apply the first rotation, we obtain the state

[TABLE]

Then, apply a rotation on the first qubit of the first register conditioned on the last register, we obtain the state

[TABLE]

Next, apply a rotation on the second qubit of the first register conditioned on the first qubit of the first register and last register, we obtain the desired state:

[TABLE]

Based on the data structure specified in Definition 1 and the efficient state preparation in Lemma 1, we construct a quantum walk operator for $H$ as follows. First define the isometry $T$ as

[TABLE]

with $\ket{\varphi_{j1}}=\ket{0}\ket{1}$ and

[TABLE]

where $\sigma_{j}=\sum_{k=0}^{N-1}|H_{jk}|$. Let $S$ be the swap operator that maps $\ket{j_{0}}\ket{b_{0}}\ket{j_{1}}\ket{b_{1}}$ to $\ket{j_{1}}\ket{b_{1}}\ket{j_{0}}\ket{b_{0}}$, for all $j_{0},j_{1}\in\{0,\ldots,N-1\}$ and $b_{0},b_{1}\in\{0,1\}$. We observe that

[TABLE]

where the second equality is ensured by the choice of the square-root as in the proof of Lemma 1. This implies that

[TABLE]

The quantum walk operator $U$ is defined as

[TABLE]

A more general characterization of the eigenvalues of quantum walks is presented in [39]. Here we give a specific proof on the relationship between the eigenvalues of $U$ and $H$ as follows.

Lemma 3**.**

Let the unitary operator $U$ be defined as in Eq. (18), and let $\lambda$ be an eigenvalue of $H$ with eigenstate $\ket{\lambda}$. It holds that

[TABLE]

where

[TABLE]

Proof.

By the fact that $T^{{\dagger}}T=I$ and $(I\otimes\bra{0})T^{{\dagger}}ST(I\otimes\ket{0})=H/\left\|H\right\|_{1}$, and $(I\otimes\bra{1})T^{{\dagger}}ST(I\otimes\ket{0})=0$, we have

[TABLE]

In order for this state being an eigenstate, it must hold that

[TABLE]

and the solution is

[TABLE]

∎

3 Linear combination of unitaries and Hamiltonian simulation

To see how to convert the quantum walk operator $U$ to Hamiltonian simulation, we first consider the generating function for the Bessel function, denoted by $J_{m}(\cdot)$. According to [2, (9.1.41)], we have

[TABLE]

where the second equality follows from Eq. (21) and the fact that $\sin(x)=(e^{ix}-e^{-ix})/2i$. This leads to the following linear combination of unitaries:

[TABLE]

where the second equality follows from the fact that $\sum_{j={-\infty}}^{\infty}J_{j}(z)=1$.

Now we consider an approximation to $e^{-izH/\left\|H\right\|_{1}}$ in terms of Eq. (26):

[TABLE]

Here the coefficients are normalized by $\sum_{j=-k}^{k}J_{j}(z)$ so that they sum to 1. This will minimize the approximation error (see the proof of Lemma 4 in Appendix B, and the normalization trick was originated in [9]). The eigenvalues of $V_{k}$ are

[TABLE]

Note that each eigenvalue of $V_{k}$ does not depend on $\pm$ as $J_{-m}(z)=(-1)^{m}J_{m}(z)$.

To bound the error in this approximation, we have the technical lemma, and the proof is shown in Appendix B.

Lemma 4**.**

Let $V_{k}$ and $V_{\infty}$ be defined as above. There exists a positive integer $k$ satisfying $k\geq|z|$ and

[TABLE]

such that

[TABLE]

In the following, we provide technical lemmas for implementing linear combination of unitaries. Suppose we are given the implementations of unitaries $U_{0}$, $U_{1}$, …, $U_{m-1}$, and coefficients $\alpha_{0},\alpha_{1},\ldots,\alpha_{m-1}$. Then the unitary

[TABLE]

can be implemented probabilistically by the technique called linear combination of unitaries (LCU) [27]. Provided $\sum_{j=0}^{m-1}|\alpha_{j}|\leq 2$, $V$ can be implemented with success probability $1/4$. To achieve this, we define the multiplexed-$U$ operation, which is denoted by $\text{multi-}U$, as

[TABLE]

We summarize the probabilistic implementation of $V$ in the following lemma, whose proof is shown in Appendix B.

Lemma 5**.**

Let $\mathrm{multi}$-$U$ be defined as above. If $\sum_{j=0}^{m-1}|\alpha_{j}|\leq 2$, then there exists a quantum circuit that maps $\ket{0}\ket{0}\ket{\psi}$ to the state

[TABLE]

where $(\ket{0}\!\bra{0}\otimes\ket{0}\!\bra{0}\otimes I)\ket{\Phi^{\bot}}=0$. Moreover, this quantum circuit uses $O(1)$ applications of $\mathrm{multi}$-$U$ and $O(m)$ 1- and 2-qubit gates.

Let $W$ be the quantum circuit in Lemma 5, and let $P$ be the projector defined as $P=\ket{0}\!\bra{0}\otimes\ket{0}\!\bra{0}\otimes I$. We have

[TABLE]

If $\sum_{j=0}^{m-1}\alpha_{j}U_{j}$ is a unitary operator, one application of the oblivious amplitude amplification operator $-W(1-2P)W^{{\dagger}}(1-2P)W$ implements $\sum_{j=0}^{m-1}U_{j}$ with certainty [8]. However, in our application, the unitary operator $\widetilde{W}$ implements an approximation of $V_{\infty}$ in the sense that

[TABLE]

with $\lVert V_{k}-V_{\infty}\rVert\leq\epsilon$. The following lemma shows that the error caused by the oblivious amplitude amplification is bounded by $O(\epsilon)$, and the proof is given in Appendix B.

Lemma 6**.**

Let the projector $P$ be defined as above. If a unitary operator $\widetilde{W}$ satisfies $P\widetilde{W}\ket{0}\ket{0}\ket{\psi}=\frac{1}{2}\ket{0}\ket{0}\widetilde{V}\ket{\psi}$ where $\lVert\widetilde{V}-V\rVert\leq\epsilon$. Then $\lVert-\widetilde{W}(I-2P)\widetilde{W}^{{\dagger}}(I-2P)\widetilde{W}\ket{0}\ket{0}\ket{\psi}-\ket{0}\ket{0}V\ket{\psi}\rVert=O(\epsilon)$.

Now we are ready to prove Theorem 1.

Proof of Theorem 1.

The proof we outline here follows closely the proof given in [9]. The intuition of this algorithm is to divide the simulation into $O(t\left\|H\right\|_{1})$ segments, with each segment simulating $e^{-iH/2}$. To implement each segment, we use the LCU technique to implement $V_{k}$ defined in Eq. (27), with coefficients $\alpha_{m}=J_{m}(z)/\sum_{j=-k}^{k}J_{j}(z)$. When $z=-1/2$, we have $\sum_{j=-k}^{k}|\alpha_{j}|<2$. Actually, this holds for all $|z|\leq 1/2$ because

[TABLE]

where the first inequality follows from the fact that $\sum_{j=-k}^{k}J_{j}(z)\geq 1-4|z/2|^{k+1}/(k+1)!$ (see [9]), the second inequality follows from the fact that $|J_{m}(z)|\leq|z/2|^{|m|}/|m|!$ (see [2, (9.1.5)]), and the third inequality uses the assumption that $|z|\leq 1/2$. Now, Lemmas 5 and 6 can be applied. By Eq. (25) and using Lemma 4, set

[TABLE]

and we obtain a segment that simulates $e^{-iH/(2\left\|H\right\|_{1})}$ with error bounded by $O(\epsilon^{\prime})$. Repeat the segment $O(t\left\|H\right\|_{1})$ times with error $\epsilon^{\prime}=\epsilon/(t\left\|H\right\|_{1})$, and we obtain a simulation of $e^{-iHt}$ with error bounded by $\epsilon$. It suffices to take

[TABLE]

By Lemma 5, each segment can be implemented by $O(1)$ application of $\text{multi-}U$ and $O(k)$ 1- and 2-qubit gates, as well as the cost for computing the coefficients $\alpha_{m}$ for $m\in\{-k,\ldots,k\}$. The cost for each $\text{multi-}U$ is $k$ times the cost for implementing the quantum walk $U$. By Lemma 1, the state in Eq. (15) can be prepared with time complexity (circuit depth) $O(n^{2}b^{5/2})$, where $b$ is the number of bit of precision. To achieve the overall error bound $\epsilon$, we choose $b=O(\log(t\left\|H\right\|_{1}/\epsilon))$. Hence the time complexity for the state preparation is $O(n^{2}\log^{5/2}(t\left\|H\right\|_{1}/\epsilon))$, which is also the time complexity for applying the quantum walk $U$. Therefore, the time complexity for one segment is

[TABLE]

Considering $O(t\left\|H\right\|_{1})$ segments, the time complexity is as claimed. ∎

Note that the coefficients $\alpha_{-k},\ldots,\alpha_{k}$ (for $k$ defined in Eq. (37)) in Lemma 5 can be classically computed using the methods in [11, 34], and the cost is $O(k)$ times the number of bits of precision, which is $O(\log(t\lVert H\rVert/\epsilon)$. This is no larger than the quantum time complexity.

4 Discussion

We presented a quantum Hamiltonian simulation algorithm whose time complexity has $\widetilde{O}(\sqrt{N})$ dependence for non-sparse Hamiltonians with access to a seemingly more powerful input model. Our technique for Hamiltonian simulation combines ideas from linear combination of quantum walks and an efficient memory model which prepares a special states to provide the improved performance. Our algorithm can be directly applied as a subroutine for the unitary implementation problem and for a quantum linear systems solver. Especially for the latter application, many quantum machine learning algorithms are based on solving linear systems in the same input model, and our algorithm implies exponentially improved error dependence over [41] for this application.

We note that the data structure in qRAM may require a large overhead in practice, if the data were not already stored in it. Furthermore it might be hard to implement such a data structure physically due to the exponential amount of quantum resources [1, 3, 18]. However, if the Hamiltonian is highly structured (i.e., the entries repeat in some pattern), the memory model could be efficiently implemented. Another potential point of criticism is the required error rate of such a device, as some computations will require an error rate per gate of $O(1/\mathrm{poly}(N))$ to retain a feasible error rate for applications [5]. Whereas, not all computations might need such low error rates [5] and hence the feasibility of our algorithm as a subroutine in an explicit algorithm must be further validated experimentally. To summarize, our algorithm inherits many problems of Grover’s search algorithm for unordered classical data, as well as many qRAM-based quantum machine learning algorithms.

However, if the above mentioned caveats can be overcome, our algorithm still supplies a polynomial speedup over known quantum algorithms for Hamiltonian simulation and furthermore allows for a polynomial speedup in comparison with the best known classical algorithms for several practical problems, such as linear regression or linear systems.

5 Acknowledgement

We thank Richard Cleve and Simone Severini for the discussion and comments on this project. We also thank anonymous reviewers for their valuable suggestions and comments on this paper. CW acknowledges financial support by a David R. Cheriton Graduate Scholarship. LW acknowledges financial support by the Royal Society through a Research Fellow Enhancement Award.

Appendix A The relation between $\left\|H\right\|_{1}$ and $\lVert H\rVert$

We prove the following proposition.

Proposition 2**.**

If $A\in\mathbb{C}^{N\times N}$ has at most $d$ non-zero entries in any row, it holds that $\left\|A\right\|_{1}\leq\sqrt{d}\lVert A\rVert$.

Proof.

First observe that $\lVert A\rVert^{2}\leq\sum_{i}\lambda_{i}(A^{{\dagger}}A)=\mathrm{Tr}\left(A^{{\dagger}}A\right)=\lVert A\rVert^{2}_{F}$, and furthermore we have that $\sum_{ij}|a_{ij}|^{2}\leq d\max_{j\in[N]}\sum_{i}|a_{ij}|^{2}$. From this we have that $\lVert A\rVert\leq\sqrt{d}\lVert A\rVert_{1}$ for a $d$-sparse $A$. By [24, Theorem 5.6.18], we have that $\lVert A\rVert_{1}\leq C_{M}(1,*)\lVert A\rVert$ for $C_{M}(1,*)=\max_{A\neq 0}\frac{\lVert A\rVert}{\lVert A\rVert_{1}}$ and using the above we have $C_{M}(1,*)\leq\sqrt{d}$. Therefore we find that $\lVert A\rVert_{1}\leq\sqrt{d}\lVert A\rVert$ as desired. ∎

It immediately follows that $\left\|A\right\|_{1}\leq\sqrt{N}\lVert A\rVert$ for dense $A$, i.e., $d=N$.

Appendix B Proofs of technical lemmas

Proof of Lemma 4.

The proof outlined here follows closely the proof of Lemma 8 in [9]. Recalling the definition of $V_{k}$ and $V_{\infty}$, we define the weights in $V_{k}$ by

[TABLE]

where $C_{k}=\sum_{l=-k}^{k}J_{l}(z)$. The normalization here is chosen so that $\sum_{m}a_{m}=1$ which will give the best result [9].

Since

[TABLE]

observe that we have two error sources. The first one comes from the truncation of the series, and the second one comes from the different renormalization of the terms which introduces an error in the first $\lvert m\rvert\leq k$ terms in the sum. We therefore start by bounding the normalization factor $C_{k}$. For the Bessel-functions for all $m$ it holds that $\lvert J_{m}(z)\rvert\leq\frac{1}{\lvert m\lvert!}\left\lvert\frac{z}{2}\right\rvert^{|m|}$, since $J_{-m}(z)=(-1)^{m}J_{m}(z)$ [2, (9.1.5)]. For $\lvert m\rvert\leq k$ we can hence find the following bound on the truncated part

[TABLE]

Since $\sum_{m}J_{m}(z)=1$, based on the normalization, we hence find that

[TABLE]

which is a lower bound on the normalization factor $C_{k}$. Since $a_{m}=\frac{J_{m}(z)}{C_{k}}$, the correction is small, which implies that

[TABLE]

and we have a multiplicative error based on the renormalization.

Next we want to bound the error in the truncation before we join the two error sources.

From Eq. (25) we know that

[TABLE]

by the normalization of $\sum_{m}J_{m}(z)$. From this we can see that we can hence obtain a bound on the truncated $J_{m}(z)$ as follows.

[TABLE]

Therefore we can upper bound the left-hand side in terms of the exact value of $V_{\infty}$, i.e. $e^{iz\lambda/\Lambda}$ if we can bound the right-most term in Eq. (46). Using furthermore the bound in Eq. (44) we obtain

[TABLE]

which reduced with $\left\lvert 2^{iz\lambda/\Lambda}-1\right\rvert\leq\lvert z\lambda/\Lambda\rvert$ and $|z|\leq k$ to

[TABLE]

We can then obtain the desired bound $\lVert V_{\infty}-V_{k}\rVert$ by reordering the above equation, and using that $\sum_{m=-k}^{k}a_{m}(z)=1$ such that we have

[TABLE]

We hence only need to bound the right-hand side.

For $\mu_{+}$ we can use that $\lvert\mu_{+}^{m}-1\rvert\leq 2\lvert m\lambda/\Lambda\rvert=:2\lvert m\nu\rvert$ and obtain the bound $2\frac{\lvert\nu\rvert}{k!}\left\lvert\frac{z}{2}\right\rvert^{k+1}$ [9]. For the $\mu_{-}$ case we need to refine the analysis and will show that the bound remains the same. Let $\nu:=\lambda/\Lambda$ as above. First observe that $J_{m}(z)\mu_{-}^{m}+J_{-m}(z)\mu_{i}^{-m}=J_{-m}(z)\mu_{+}^{-m}+J_{m}(z)\mu_{+}^{m}$, and it follows that

[TABLE]

Therefore we only need to treat the $\mu_{+}$ case.

[TABLE]

Using this bound we hence obtain from Eq. (49),

[TABLE]

In order for the above equation being upper-bounded by $\epsilon$, it suffices to choose some $k$ that is upper bounded as claimed. ∎

Proof of Lemma 5.

Let $s=\sum_{j=0}^{m-1}|\alpha_{j}|$. We first define the unitary operator $B$ to prepare the coefficients:

[TABLE]

Define the unitary operator $W$ as $W=(B^{{\dagger}}\otimes I)(I\otimes\mbox{multi-}U)(B\otimes I)$. We claim that $W$ performs the desired mapping, as

[TABLE]

where $\ket{\Phi^{\bot}}$ is a state satisfying $(\ket{0}\!\bra{0}\otimes\ket{0}\!\bra{0}\otimes I)\ket{\Phi^{\bot}}=0$, and $\gamma$ is some normalization factor.

The number of applications of $\mbox{multi-}U$ is constant, as in the definition of $W$. To implement the unitary operator $B$, $O(m)$ 1- and 2-qubit gates suffice. ∎

Proof of Lemma 6.

We have

[TABLE]

Because $\lVert\widetilde{V}-V\rVert\leq\epsilon$ and $V$ is a unitary operator, we have $\lVert\widetilde{V}^{{\dagger}}\widetilde{V}-I\rVert=O(\epsilon)$. Therefore, we have

[TABLE]

Thus

[TABLE]

∎