Quantum algorithm and circuit design solving the Poisson equation

We start with the one-dimensional case to introduce the matrix L_h that we will use later in expressing the d-dimensional discretization of the Laplacian, using Kronecker products. We have

$$\begin{eqnarray} & &-\frac{\mathrm{d}^2u(x)}{\mathrm{d}x^2}=f(x), \quad x\in(0,1), \nonumber\\ \\ & &u(0)=u(1)=0, \nonumber \end{eqnarray} \tag{ 3 }$$

where f is a given smooth function and u is the solution we want to compute. We discretize the problem with mesh size h = 1/M and we compute an approximate solution v at M + 1 grid points x_i = ih, $i=0,\ldots ,M$. Let u_i = u(x_i) and f_i = f(x_i), $i=0,\ldots ,M$.

Using finite differences at the grid points to approximate the second derivative (3) becomes

$$\begin{equation} -\frac{\mathrm{d}^2u(x)}{\mathrm{d}x^2}\Bigg|_{x=x_i}=\frac{2u_i-u_{i-1}-u_{i+1}}{h^2}-\xi_i, \end{equation} \tag{ 4 }$$

where ξ_i is the truncation error and can be shown to be $O(h^2\Vert\frac {d^4u}{\mathrm {d}x^4}\Vert_{\infty })$ if f has the fourth derivative uniformly bounded by a constant [31].

Ignoring the truncation error, we solve

$$\begin{equation} h^{-2}(-v_{i-1}+2v_i-v_{i+1})=f_i, \quad 0<i<M. \end{equation} \tag{ 5 }$$

With boundary conditions v₀ = 0 and v_M = 0, we have M − 1 equations and M − 1 unknowns v₁,...,v_M−1:

$$\begin{equation} h^{-2} \, L_h \left( \begin{array}{@{}c@{}} v_1 \\ \vdots \\ \vdots \\ v_{M-1} \end{array} \right) := h^{-2} \left( \begin{array}{@{}cccc@{}} 2 & -1 & \quad & 0 \\ -1 & \ddots & \ddots & \\ \quad & \ddots & \ddots & -1 \\ 0 & \quad & -1 & 2 \end{array} \right) \left( \begin{array}{@{}c@{}} v_1 \\ \vdots \\ \vdots \\ v_{M-1} \end{array} \right) = \left( \begin{array}{@{}c@{}} f_1 \\ \vdots \\ \vdots \\ f_{M-1} \end{array} \right), \end{equation} \tag{ 6 }$$

where L_h is the tridiagonal (M − 1) × (M − 1) matrix above; for the properties of this matrix, including its eigenvalues and eigenvectors see [31, section 6.3].

In two dimensions the Poisson equation is

$$\begin{eqnarray} &&\begin{array}{l} \displaystyle-\frac{\partial^2{u(x,y)}}{\partial{x^2}}-\frac{\partial^2{u(x,y)}}{\partial{y^2}}=f(x,y),\quad (x,y)\in (0,1)^2,\\ \displaystyle u(x,0)=u(0,y)=u(x,1)=u(1,y)=0, \quad x,y\in [0,1]. \end{array} \end{eqnarray} \tag{ 7 }$$

We discretize this equation using a grid with mesh size h = 1/M; see figure 1. Each node is indexed u_j,k, j,k∈{1,2,...,M} (figures 1(a) and (b)). We approximate the second derivatives using

$$\begin{eqnarray*} \frac{\partial^2u}{\partial x^2}(x,y) &\approx&\frac{u(x-h,y)-2u(x,y)+u(x+h,y)}{h^2}, \\ \frac{\partial^2u}{\partial y^2}(x,y) &\approx&\frac{u(x,y-h)-2u(x,y)+u(x,y+h)}{h^2}. \end{eqnarray*} \noindent$$

Zoom In Zoom Out Reset image size

Omitting the truncation error, and denoting the discretized Laplacian by −Δ_h we are led to solve

$$\begin{equation} h^{-2}\left((-v_{j-1,k}+2v_{j,k}-v_{j+1,k})+(-v_{j,k-1}+2v_{j,k}-v_{j,k+1})\right)=f_{j,k}, \end{equation} \noindent \tag{ 8 }$$

where f_j,k = f(jh,kh), j,k = 1,2,...,M − 1 and v_j,k = 0 if j or k∈{0,M} i.e. when we have a point that belongs to the boundary.

Using the fact that the solution is zero at the boundary, we reindex (8) to obtain

$$\begin{equation} h^{-2}(4v_i-v_{i-1}-v_{i+1}-v_{i-M+1}-v_{i+M-1})=f_i, \quad i= 1,2,\ldots,(M-1)^2. \end{equation} \noindent \tag{ 9 }$$

Equivalently, we denote this system by

$$\begin{eqnarray*} &&-\Delta_h \vec v = \vec{f_h}, \end{eqnarray*} \noindent$$

where Δ_h is the discretized Laplacian.

For example, when M = 4, as in figure 1, we have that $\vec {v}=[v_1,\ldots ,v_9]^{\mathrm {T}}$. Furthermore (9) becomes

$$\begin{equation} h^{-2} A \left( \begin{array}{@{}c@{}} v_1 \\ \vdots \\ v_9 \end{array} \right) :=h^{-2} \left( \begin{array}{@{}ccc@{}} B & -I & \\ -I & B & -I \\ & -I & B \end{array} \right) \left( \begin{array}{@{}c@{}} v_1 \\ \vdots \\ v_9 \end{array} \right)= \left( \begin{array}{@{}c@{}} f_1 \\ \vdots \\ f_9 \end{array} \right), \end{equation} \tag{ 10 }$$

where I is the 3 × 3 identity matrix, B is

$$\begin{equation*} \left( \begin{array}{@{}ccc@{}} 4 & -1 & \\ -1 & 4 & -1 \\ & -1 & 4 \end{array} \right). \end{equation*}$$

A is a Hermitian matrix with a particular block structure that is independent of M.

In particular, on a square grid with mesh size h = 1/M we have

$$\begin{equation} -\Delta_h = h^{-2} A \end{equation} \tag{ 11 }$$

and A can be expressed in terms of L_h as follows:

$$\begin{equation} A= \left( \begin{array}{@{}cccccc@{}} L_h+2I & -I & 0 & \cdots & \cdots & 0 \\ -I & L_h+2I & -I & 0 & \cdots & 0 \\ 0 & -I & \ddots & \ddots & 0 & \vdots \\ \vdots & 0 & \ddots & \ddots & -I & 0 \\ \vdots & \vdots & 0 & -I & L_h+2I & -I \\ 0 & 0 & \cdots & 0 & -I & L_h+2I \end{array} \right) \end{equation} \tag{ 12 }$$

and its size is (M − 1)² × (M − 1)² [31].

Recall that L_h is the (M − 1) × (M − 1) matrix shown in (6) and I is the (M − 1) × (M − 1) identity matrix. Moreover, A can be expressed using Kronecker products as follows:

$$\begin{equation} A=L_h\otimes{I}+I\otimes{L_h} . \end{equation} \tag{ 13 }$$

We now consider the problem in d dimensions. Consider the Laplacian

$$\begin{eqnarray*} &&\Delta = \sum_{k=1}^d \frac{\partial^2}{\partial x_k^2}. \end{eqnarray*}$$

We discretize Δ on a grid with mesh size h = 1/M using divided differences.

As before, this leads to a system of linear equations

$$\begin{equation} -\Delta_h\vec{v} = \vec{f_h}. \end{equation} \tag{ 14 }$$

Note that −Δ_h = h⁻²A is a symmetric positive definite matrix and A is given by

$$\begin{equation*} A=\underbrace{L_h\otimes{I}\otimes \cdots \otimes{I}}_{d\,{\rm matrices}}+ I\otimes{L_h}\otimes I\otimes \cdots \otimes{I}+\cdots + I\otimes\cdots \otimes I\otimes{L_h}, \end{equation*}$$

and has size (M − 1)^d × (M − 1)^d. L_h is the (M − 1) × (M − 1) matrix shown in (6) and I is the (M − 1) × (M − 1) identity matrix. See [31] for details.

Observe that the matrix exponential has the form

$$\begin{equation} \mathrm{e}^{\mathrm{i}A \gamma}=\underbrace{\mathrm{e}^{\mathrm{i}L_h\gamma} \otimes \cdots \otimes \mathrm{e}^{\mathrm{i}L_h\gamma}}_{d \,\, \mathrm{matrices}}, \end{equation} \tag{ 15 }$$

for all $\gamma \in \mathbb {R}$, where $\mathrm {i} = \sqrt {-1}$. We will use this fact later to derive the quantum circuit solving the linear system.

We carry out the error analysis to obtain the implementation details. For d = 1 the eigenvalues of the second derivative are

$$\begin{eqnarray*} &&4 M^2 \sin^2(j\pi/(2M)), \quad j=1,\ldots,M-1. \end{eqnarray*}$$

For d > 1, the eigenvalues of −Δ_h are given by sums of the one-dimensional eigenvalues, i.e.

$$\begin{eqnarray*} &&\sum_{k=1}^d \left[4 M^2 \sin^2(j_k\pi/(2M))\right], \quad j_k=1,\ldots,M-1, \quad k=1,\ldots,d. \end{eqnarray*}$$

We consider them in non-decreasing order and denote them by λ_j, $j=1,\ldots , (M-1)^d.$ Then λ₁ = 4dM²sin²(π/(2M)) is the minimum eigenvalue and λ_(M−1)^d = 4dM²sin²(π(M − 1)/(2M)) ⩽ 4dM² is the maximum eigenvalue.

Define E by

$$\begin{equation} \log_2 E=\lceil \log_2 d\rceil + \log_2 (4M^2). \end{equation} \tag{ 16 }$$

Then the eigenvalues are bounded from above by E. Recall that we have already assumed that M is a power of two. Then $E=2^{\lceil {\log _2{d}}\rceil }4M^2\in \mathbb {N}$.

Note that the implementation accuracy of the eigenvalues determines the accuracy of the system solution.

Our algorithm uses approximations $\hat \lambda _j$, such that $|\lambda _j - \hat \lambda _j|\leqslant \frac {17 E}{2^\nu } \leqslant \varepsilon$; see theorem B.2 in appendix B. We use n = log₂ E + ν bits to represent each eigenvalue, of which the most significant log₂ E bits hold each integer part and the remaining bits hold each fractional part. Without loss of generality, we can assume that 2^ν ≫ E. More precisely, we consider an approximation $\hat \Delta _h$ of matrix Δ_h such that the two matrices have the same eigenvectors while their eigenvalues differ by at most ε.

We use phase estimation with the unitary matrix $\mathrm {e}^{-\mathrm {i}\hat \Delta _h t_0 /E}$ whose eigenvalues are $\mathrm {e}^{2\pi \mathrm {i} \hat \lambda _j t_0/(E 2\pi )}$. Setting t₀ = 2π we obtain the phases $\phi _j = \hat{\lambda}_j /E \in [0,1)$. The initial state of phase estimation is (figure 2)

$$\begin{eqnarray*} &&|{0}\rangle^{\otimes n}|{f_h}\rangle = \sum_{j=1}^{(M-1)^d} \beta_j |{0}\rangle^{\otimes n}|{u_j}\rangle, \end{eqnarray*}$$

where |u_j〉 is the jth eigenvector of −Δ_h and $\beta _j = \left <u_j|f_h\right >$, for j = 1,2,...,(M − 1)^d. Since we are using finite bit approximations of the eigenvalues, we have

$$\begin{eqnarray*} &&\phi_j = \frac{\hat \lambda_j}{E}= \frac{\hat \lambda_j 2^\nu}{2^n}. \end{eqnarray*}$$

Then ϕ_j 2ⁿ is an integer and phase estimation succeeds with probability 1 (see [36, section 5.2, p 221] for details).

The state prior to the application of the inverse Fourier transform in the phase estimation is

$$\begin{equation} \sum_{j=1}^{(M-1)^d} \beta_j \frac{1}{2^{n/2}} \sum_{k=0}^{2^n-1} \mathrm{e}^{2\pi \mathrm{i} \phi_j k} | k\rangle |{u_j}\rangle. \end{equation} \tag{ 17 }$$

After the application of the inverse Fourier transform to the first n qubits we obtain

$$\begin{eqnarray*} &&\sum_{j=1}^{(M-1)^d} \beta_j |{k_j}\rangle|{u_j}\rangle, \end{eqnarray*}$$

where

$$\begin{equation} k_j = 2^n \phi_j = 2^n \hat\lambda_j / E=\hat\lambda_j2^\nu\in\mathbb{N}. \end{equation} \tag{ 18 }$$

Now we need to compute the reciprocals of the eigenvalues. Observe that

$$\begin{eqnarray*} \lambda_1/d &=& 4M^2 \sin^2(\pi/(2M)) = 4 M^2 ( \pi / (2M) + O(M^{-3}))^2 \\ &=& \pi^2 + O(M^{-2}) > 5, \end{eqnarray*}$$

where the last inequality holds trivially for sufficiently large M. This implies $\hat \lambda _j / C_d \geqslant \hat \lambda _1 / C_d \geqslant 4$, where C_d = 2^{⌊log₂ d⌋}, for sufficiently large M. We obtain $k_j = 2^n \hat \lambda _j / E \geqslant \hat \lambda _1 \geqslant 4 C_d$.

Append b qubits initialized to |0〉 on the left (Reg. L in figure 2), to obtain

$$\begin{eqnarray*} &&\sum_{j=1}^{(M-1)^d} \beta_j |{0}\rangle^{\otimes b}|{k_j}\rangle|{u_j}\rangle. \end{eqnarray*}$$

Note that from (18) k_j, $\hat {\lambda }_j$ and $\hat \lambda _j / C_d$ have the same bit representation. The difference between the integer k_j and the other two numbers is the location of the decimal point; it is located after the most significant log₂ E bit in $\hat {\lambda }_j$, and after the most significant log₂(E/C_d) bit in $\hat \lambda _j/C_d$. Therefore, we can use the labels |k_j〉, $|{\hat \lambda _j}\rangle$ and $|{\hat {\lambda }_j / C_d}\rangle$ interchangeably, and write the state above as

$$\begin{eqnarray*} &&\sum_{j=1}^{(M-1)^d} \beta_j |{0}\rangle^{\otimes b}|{\hat \lambda_j/C_d}\rangle|{u_j}\rangle. \end{eqnarray*}$$

Now we need to compute $h_j := h (\hat \lambda _j / C_d) = C_d / \hat \lambda _j$. We do this using Newton iteration. We explain the details in section 4.2. We obtain an approximation $\hat h_j$ such that

$$\begin{equation} \left| \hat h_j - h_j \right| \leqslant \varepsilon_0^2, \end{equation} \tag{ 19 }$$

where $\varepsilon _0 = \min \{\varepsilon , E^{-1}\}$. We store this approximation in the register composed of the leftmost b = 3⌈log₂ ε⁻¹₀⌉ qubits.

This leads to the state

$$\begin{eqnarray*} &&\sum_{j=1}^{(M-1)^d} \beta_j |{\hat h_j}\rangle|{\hat {\lambda}_j/C_d}\rangle|{u_j}\rangle. \end{eqnarray*}$$

We append, on the left, a qubit initialized at |0〉 (Anc. in figure 2). We get

$$\begin{eqnarray*} &&\sum_{j=1}^{(M-1)^d} \beta_j | 0\rangle |{\hat h_j}\rangle|{\hat {\lambda}_j/C_d}\rangle|{u_j}\rangle. \end{eqnarray*}$$

We need to perform the conditional rotation

$$\begin{eqnarray*} &&R | 0\rangle |{\omega}\rangle = \left( \omega | 1\rangle + \sqrt{1 - \omega^2} | 0\rangle \right) |{\omega}\rangle, \quad 0<\omega < 1. \end{eqnarray*}$$

For this, we will approximate the first qubit by

$$\begin{eqnarray*} &&\omega^\prime | 1\rangle + \sqrt{1 - (\omega^\prime)^2} | 0\rangle, \end{eqnarray*}$$

with |ω − ω'| ⩽ ε²₁, $\varepsilon _1= \min \{ \varepsilon , 1/(4M^2)\}$. We discuss the cost of implementing this approximation in section 4.3.

The result of approximating the conditional rotation is to obtain $|{\tilde h_j}\rangle$, where $\tilde h_j$ is a q = Θ(log₂ ε⁻¹₁) bit number less than 1 satisfying $|\tilde h_j - \hat h_j|\leqslant \varepsilon _1^2$ and, therefore,

$$\begin{equation} |\tilde h_j - h_j|\leqslant \varepsilon_0^2 + \varepsilon_1^2, \end{equation} \tag{ 20 }$$

for each $j=1,\ldots ,(M-1)^d$.

Ignoring the ancilla qubits needed for implementing the approximation of the conditional rotation, we have the state

$$\begin{eqnarray*} &&\sum_{j=1}^{(M-1)^d} \beta_j \left( \tilde h_j | 1\rangle + \sqrt{1 - \tilde h_j^2} | 0\rangle \right) |{\hat h_j}\rangle|{\hat {\lambda}_j/C_d}\rangle|{u_j}\rangle. \end{eqnarray*}$$

Uncomputing all the qubits except the leftmost gives the state

$$\begin{eqnarray*} &&| \psi\rangle: = \sum_{j=1}^{(M-1)^d} \beta_j \left( \tilde h_j | 1\rangle + \sqrt{1 - \tilde h_j^2} | 0\rangle \right) |{0}\rangle^{\otimes b} |{0}\rangle^{\otimes n} |{u_j}\rangle. \end{eqnarray*}$$

Let P₁ = |1〉〈1|⊗I be the projection acting non-trivially on the first qubit. The system $-\Delta _h \vec v = \skew6\vec f_h$ has solution $\sum _{j=1}^{(M-1)^d} \beta _j\frac {1}{{\lambda }_j} |{u_j}\rangle$. We derive the error as follows:

$$\begin{eqnarray} & & C_d^{-1}\left\| \sum_{j=1}^{(M-1)^d} b_j \frac{C_d}{\lambda_j} |{1}\rangle |{0}\rangle^{\otimes{(b+n)}} |{u_j}\rangle - P_1|{\psi}\rangle \right\| \nonumber\\ & & =C_d^{-1}\left\|\sum_{j=1}^{(M-1)^d} \beta_j \frac{C_d}{\lambda_j} |{u_j}\rangle - \sum_{j=1}^{(M-1)^d} \beta_j \tilde h_j |{u_j}\rangle\right\| \nonumber \\ & & =C_d^{-1}\left\|\sum_{j=1}^{(M-1)^d} \beta_j \frac{C_d}{\lambda_j} |{u_j}\rangle - \sum_{j=1}^{(M-1)^d} \beta_j (\tilde h_j - h_j + h_j)|{u_j}\rangle\right\| \nonumber \\ & &\leqslant\left\|\sum_{j=1}^{(M-1)^d} \beta_j \left( \frac{1}{\lambda_j} - \frac{1}{\hat \lambda_j} \right)|{u_j}\rangle \right\| + \varepsilon_0^2 + \varepsilon_1^2 \leqslant \frac{17 E}{2^\nu} + \varepsilon_0^2 +\varepsilon_1^2, \end{eqnarray} \tag{ 21 }$$

where the second from last inequality is obtained using (20) and the last inequality is due to the fact that

$$\begin{eqnarray*} &&\left|\frac{1}{\lambda} - \frac{1}{\hat\lambda}\right| \leqslant |\lambda - \hat \lambda|, \quad \lambda , \hat \lambda > 1. \end{eqnarray*}$$

Setting ν = ⌈log₂(17E/ε)⌉ gives error ε(1 + o(1)) and the number of matrix exponentials used by the algorithm is O(log₂(E/ε)). Therefore, if we measure the first qubit of the state |ψ〉 and the outcome is 1 the state collapses to a normalized solution of the linear system.

In this part we deal with the computation of the reciprocals of the eigenvalues, which is marked as the 'INV' module in figure 2. For this we use Newton iteration to approximate v⁻¹, v > 1. We perform s iterative steps and obtain the approximation $\hat x_s$. The input and the output of each iterative step are b bit numbers. All of the calculations in each step are performed in fixed precision arithmetic. The initial approximation is $\hat x_0 = 2^{-p}$, 2^p−1 < v ⩽ 2^p. (We use the notation $\hat x_i$ to emphasize that these values have been obtained by truncating a quantity x_i to b bits of accuracy, i = 0,...,s.)

Theorem B.1 of appendix B gives the error of Newton iteration which is

$$\begin{eqnarray*} &&|\hat x_s - v^{-1} | \leqslant \varepsilon_0^2 \leqslant \varepsilon, \end{eqnarray*}$$

where we have $\varepsilon _0=\min \{ \varepsilon , E^{-1}\}$, s = ⌈log₂ log₂(2/ε²₀)⌉ and the number of bits satisfies b ⩾ 2⌈log₂ ε⁻¹₀⌉ + O(log₂ log₂ log₂ ε⁻¹₀).

Therefore, it suffices that the module of the quantum circuit that computes 1/λ_j carries each iterative step with 3⌈log₂ ε⁻¹₀⌉ qubits of accuracy.

The quantum circuit computing the initial approximation $\hat x_0$, of the Newton iteration is given in figure 3. The second register holds |v〉 and is n qubits long, of which the first log₂(E/C_d) qubits represent the integer part of v and the remaining ones its fractional part. The first register is b qubits long. Recall that $\hat {\lambda }_j / C_d \geqslant 4$. So input values below 4 do not correspond to meaningful eigenvalue estimates and we do not need to compute their reciprocals altogether; they can be ignored. Hence the circuit implements the unitary transformation |0〉^⊗b|v〉 → |0〉^⊗b|v〉, if the first log₂(E/C_d) − 2 bits of v are all zero. Otherwise, it implements the initial approximation $\hat x_0$ through the transformation $|{0}\rangle ^{\otimes b} | v\rangle \to |{\hat x_0}\rangle | v\rangle$.

Zoom In Zoom Out Reset image size

Each iteration step x_i+1 = −vx²_i + 2x_i is implemented using a quantum circuit of the form shown in figure 4 that computes $|{\hat x_{i}}\rangle | v\rangle \to |{\hat x_{i+1}}\rangle | v\rangle$. This involves quantum circuits for addition and multiplication which have been studied in the literature [37].

Zoom In Zoom Out Reset image size

The register holding |v〉 is n qubits long and the register holding the $|{\hat x_i}\rangle$ and $|{\hat x_{i+1}}\rangle$ is b qubits long. Note that internally the modules performing the iteration steps may use more than b qubits, say, double precision, so that the addition and multiplication operations required in the iteration are carried out exactly and then return the most significant b qubits of the result. The total number of qubits required for the implementation of each of these modules is O(log ε⁻¹₀) and the total number of gates is a low degree polynomial in log ε⁻¹₀.

We now consider the implementation of the controlled rotation

$$\begin{eqnarray*} &&R | 0\rangle |{\omega}\rangle = \left( \omega | 1\rangle + \sqrt{1 - \omega^2} | 0\rangle \right) |{\omega}\rangle, \quad 0<\omega < 1. \end{eqnarray*}$$

Assume for a moment that we have obtained |θ〉, a q qubit state, corresponding to an angle θ such that sin θ approximates ω. Then we can use controlled rotations R_y about the y-axis to implement R. We consider the binary representation of θ and have

$$\begin{eqnarray*} &&\theta =\theta_1 \ldots \theta_q = \sum_{j=1}^q \theta_j 2^{-j}, \quad \theta_j \in \{0,1\}. \end{eqnarray*}$$

Then

$$\begin{eqnarray*} R_y(2\theta) &=& \mathrm{e}^{-\mathrm{i} \theta Y} = \left( \begin{array}{@{}cc@{}} \sqrt{1-\sin^2\theta} & -\sin\theta \\ \sin \theta & \sqrt{1-\sin^2\theta} \end{array} \right) \\ &=& \prod_{j=1}^q \mathrm{e}^{-\mathrm{i} Y \theta_j / 2^j} = \prod_{j=1}^q R_y^{\theta_j} \left( 2^{1-j}\right) , \end{eqnarray*}$$

where Y is the Pauli Y operator and θ∈[0,π/2]. The detailed circuit is shown in figure 5.

Zoom In Zoom Out Reset image size

We now turn to the algorithm that calculates |θ〉 from |ω〉. Since ω corresponds to the reciprocal of an approximate eigenvalue of the discretized Laplacian, we know that sin⁻¹(ω) belongs to the first quadrant and sin⁻¹(ω) = Ω(1/M²). Therefore, we can find an angle θ such that |sin(θ) − ω| ⩽ ε²₁, $\varepsilon _1=\min \{ 1/(4M^2),\varepsilon \}$, using bisection and an approximation of the sine function.

In appendix B we show the error in approximating the sine function using fixed precision arithmetic. In section 5 we show the details of the resulting quantum algorithm computing the approximation to the sine function. These results, with a minor adjustment in the number of bits needed can be used here. We will not deal with the details of the quantum algorithm for the sine function in this section since we present them in section 5 that deals with the simulation of Poisson's matrix. We will only describe the steps of the algorithm and its cost.

An evaluation at the midpoint of an interval yields a value that satisfies either the condition of step 3, or that of step 4, or |s_θ − ω| ⩽ ε²₁/2. If at any time both the conditions of steps 3 and 4 are false then θ will not change its value until the end. Then, at the end, we have |sin(θ) − ω| ⩽ |sin(θ) − s_θ| + |s_θ − ω| ⩽ ε²₁, since the error in computing the sine is ε²₁/2. On the other hand, if θ is updated until the very end of the algorithm the final value of theta also satisfies |sin(θ) − w| ⩽ ε²₁, because in the final interval we have |sin(θ) − ω| ⩽ |θ − sin⁻¹(ω)| ⩽ ε²₁.

In a way similar to that of propositions 1 and 2 of appendix B we carry out the steps of the algorithm in q bit fixed precision arithmetic, $q=\max \{ 2\nu + 9, 13 + \nu + 2\log _2 M\}$ and sufficiently large ν to satisfy the accuracy requirements. (The last expression for q is slightly different form that in proposition 2 because it accounts for the fact that in the case we are dealing with here the angle is Ω(1/M²).) This gives us an approximation to the sine with error 2^−(ν−1). We set

$$\begin{eqnarray*} &&\nu = \lceil \log_2 \varepsilon_1^{-2} \rceil+ 1. \end{eqnarray*} \noindent$$

Thus ν and q are both Θ(log₂ ε₁).

The algorithm for the sine function is based on an approximation of the exponential function using repeated squaring. Each square requires O(q²) quantum operations and O(q) qubits. This is repeated ν times before the approximation to the sine is obtained. Thus the cost of one bisection step requires O(νq²) quantum operations and O(νq) qubits. So, in terms of ε₁, the total cost of bisection is proportional to (log₂ ε⁻¹₁)⁴ quantum operations and (log₂ ε⁻¹₁)³ qubits.

We start with the implementation of e^ih−2L_hγ, γ = 2π2^t/E, E = 4M² when d = 1 and t = 0,1,...,n − 1, where n is the number of qubits in register C; see (17). The form of L_h is shown in (6) and is positive definite. It is a Toeplitz matrix and it is known that this type of matrix can be diagonalized via the sine transform S [42]. We have L_h = SΛS, where Λ is an (M − 1) × (M − 1) diagonal matrix containing the eigenvalues 4sin²(jπ/(2M)), j = 1,...,M − 1, of L_h and S = {S_i,j}_{i,j=1,2,...,M−1} is the sine transform where $S_{i,j}=\sqrt {\frac {2}{M}}$sin $(\frac {{\pi }ij}{M})$, i,j = 1, $\ldots$, M − 1. Thus

$$\begin{equation} \mathrm{e}^{\mathrm{i}h^{-2} L_h \gamma } = S \,\mathrm{e}^{\mathrm{i} h^{-2} \Lambda \gamma} S. \end{equation} \noindent \tag{ 23 }$$

The relationship between the sine and cosine transforms and the Fourier transform can be found in [40, theorem 3.10].

In particular, using the notation in [40], we have

$$\begin{equation} T_{M}^\dagger{F_{2M}}T_{M}=C_{M+1}\oplus({-\mathrm{i}S_{M-1}}) = \left(\begin{array}{@{}cc@{}} C_{M+1} & 0 \\ 0 & -\mathrm{i} S_{M-1} \end{array} \right), \end{equation} \noindent \tag{ 24 }$$

where C_M+1, S_M−1 denote the cosine and sine transforms, and the subscripts M − 1 and M + 1 emphasize the size of the respective matrix. F_2M is the 2M × 2M matrix of the Fourier transform. The matrix T_M has size 2M × 2M and is given by (25)

$$\begin{equation} T_{M}= \left( \begin{array}{@{}cccccccc@{}} 1 & & & & & & & \\ & \displaystyle\frac{1}{\sqrt{2}} & & & & \displaystyle\frac{\mathrm{i}}{\sqrt{2}} & & \\ & & \ddots & & & & \ddots & \\ & & & \displaystyle\frac{1}{\sqrt{2}} & & & & \displaystyle\frac{\mathrm{i}}{\sqrt{2}} \\ & & & & 1 & & & \\ & & & \displaystyle\frac{1}{\sqrt{2}} & & & & -\displaystyle\frac{\mathrm{i}}{\sqrt{2}} \\ & & \vdots & & & & \vdots & \\ & \displaystyle\frac{1}{\sqrt{2}} & & & & -\displaystyle\frac{\mathrm{i}}{\sqrt{2}} & & \end{array} \right). \end{equation} \noindent \tag{ 25 }$$

The quantum circuits for implementing the unitary transformation T_M is discussed in [38]. The action of T_M can be described by [38]

$$\begin{equation} \begin{array}{@{}l@{}} T_M|0x\rangle = \displaystyle\frac{1}{\sqrt{2}}|0x\rangle+\frac{1}{\sqrt{2}}|1x'\rangle, \\ T_M|1x\rangle = \displaystyle\frac{\mathrm{i}}{\sqrt{2}}|0x\rangle-\frac{\mathrm{i}}{\sqrt{2}}|1x'\rangle, \end{array} \end{equation} \tag{ 26 }$$

where i² = −1, x is an n-bit number ranging 1 ⩽ x < 2ⁿ and x' denotes its complement of 2 i.e. x' = 2ⁿ − x. The basic idea of implementing T_M is to separate its operation into an operator D, which ignores the complement of 2 in T_M, and a controlled permutation π, which transforms the state |bx〉 to |bx'〉 only if b is 1. Therefore the action of D and π can be written as

$$\begin{equation} \begin{array}{@{}l@{}} \displaystyle D|0x\rangle = \frac{1}{\sqrt{2}}|0x\rangle+\frac{1}{\sqrt{2}}|1x\rangle, \\ \displaystyle D|1x\rangle = \frac{\mathrm{i}}{\sqrt{2}}|0x\rangle-\frac{\mathrm{i}}{\sqrt{2}}|1x\rangle, \\ \displaystyle \pi|0x\rangle = |0x\rangle, \\ \displaystyle \pi|1x\rangle = |1x'\rangle. \end{array} \end{equation} \tag{ 27 }$$

Clearly, T_M = Dπ and the overall circuit for implementing operation T_M is shown in figure 8.

By (24) the sine transform S can be implemented by cascading the quantum circuits in figure 8 with the circuit for the Fourier transform [36]. An ancilla bit is added to register b. It is kept in the state |1〉 in order to select the lower-right block

$$\begin{equation} { \left( \begin{array}{@{}cc@{}} a & 0 \\ 0 & -\mathrm{i}S_{M-1} \end{array} \right) } \end{equation} \tag{ 28 }$$

from the unitary operation T^†_MF_2MT_M (24), $a\in \mathbb {C}$. Considering the state |f_h〉, that corresponds to the right-hand side of (6), and for $b_i = \left <i|f_h\right >$ we have

$$\begin{equation} { ({0},\underbrace{b_1,b_2,\ldots ,b_{M-1}}_{ \begin{array}{@{}c@{}}\scriptsize \textrm{ values on the} \\ \scriptsize \textrm{$(M-1)$ nodes} \end{array} }) }, \end{equation} \tag{ 29 }$$

then the element a in equation (28) has no effect, and the circuit in figure 6 is equivalent to applying (S_M−1 e^2πiΛ2t/E S_M−1) onto the (M − 1) elements of |f_h〉. This is also equivalent to simulating the Hamiltonian e^{2πih−2Λ2t/E} with the state |f_h〉 stored in register b.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We implement ${\mathrm { e}}^{2\pi \mathrm {i} h^{-2} \hat \Lambda 2^t / E}$ where $\hat \Lambda = \{\hat {\lambda }_j\}_{j = 1,\ldots ,M-1}$ is a diagonal matrix approximating Λ = {λ_j}_{j=1,...,M−1}.

We obtain each $\hat {\lambda }_j$, j = 1,...,M − 1 by the following algorithm. The general idea is to approximate sin x = ℑ(e^ix) = ℑ((e^ix/r)^r) with W^r where W = 1 − ix/r + x²/r² is the Taylor expansion of e^ix/r up to the second order term. W^r is computed efficiently in fixed point arithmetic using repeated squaring. The detailed steps are the following.

Eigenvalue simulation algorithm (ESA)

• (i)  
  Let r = 2^ν+7 where ν is positive integer which is related to the accuracy of the result. The inputs and the outputs of the modules below are $s = \max \{2\nu +9, 11 + \nu + \log _2 M\}$ bit numbers. Internally the modules may carry out calculations in higher precision O(s), but the results are returned using s bits. This value of s follows from the error estimates in proposition B.2.
• (ii)  
  We perform the transformation

  $$\begin{eqnarray*} &&| j\rangle |{0}\rangle^{\otimes s}\rightarrow | j\rangle \underbrace{|{\hat y_j = \hat x_j / r}\rangle}_{s \,\, \rm qubits}, \end{eqnarray*}$$

  where $\hat x_j$ is the s bit truncation of $x_j = \frac {\pi j}{2M}$. Note that y_j = x_j/r∈(0,1) and $\hat y_j$ is the s bit truncation of y_j. Recall that r ⩾ 2 and 2M are powers of 2. Calculations are to be performed in fixed precision arithmetic, so division does not actually need to be performed. All one needs to do is multiply j by π with O(s) bits of accuracy, keeping in track the position of the decimal point and then take the most significant s bits of the result.
• (iii)  
  We compute the real and imaginary parts of the complex number $\hat W_1$ by truncating, if necessary, the respective parts of $\hat W_0 = 1 - \hat y^2 + \mathrm {i} \hat y$ to s bits of accuracy; see (B.7) in proposition B.1. This is expressed by the transformation

  $$\begin{eqnarray*} &&|{\hat y_j}\rangle |{0}\rangle^{\otimes s} |{0}\rangle^{\otimes s} \rightarrow |{\hat y_j}\rangle |{\Re (\hat W_1)}\rangle |{\Im (\hat W_1)}\rangle. \end{eqnarray*}$$

  Note that since $|{\hat y_j}\rangle$ is s qubits long, $\hat W_0$ can be computed exactly using double precision and ancilla qubits and the final result can be returned in s qubits.Complex numbers are implemented using two registers, holding the real and imaginary parts. Complex arithmetic is performed by computing the real and imaginary parts of the result.
• (iv)  
  We compute $\hat W_r$ approximating $\hat W_1^r$ using repeated squaring. Each step of this procedure is accomplished by the transformation

  $$\begin{eqnarray*} &&|{\Re (\hat W_{2^j})}\rangle|{\Im (\hat W_{2^j})}\rangle |{0}\rangle^{\otimes s} |{0}\rangle^{\otimes s} \rightarrow |{\Re (\hat W_{2^j})}\rangle|{\Im (\hat W_{2^j})}\rangle |{\Re (\hat W_{2^{j+1}})}\rangle|{\Im (\hat W_{2^{j+1}})}\rangle, \end{eqnarray*}$$

  which describes the steps in (B.7). The registers holding real and imaginary parts of the numbers are s qubits long.
• (v)  
  $\Im (\hat W_r)$ approximates sin(πj/(2M)) with error 2^−(ν−1). Hence $\Im ^2(\hat W_r)$ approximates the sin²(πj/(2M)). We compute the square of $\Im (\hat W_r)$ exactly and multiply it by 4M² (this involves only shifting). We keep the most significant ν + log₂(4M²) bits of the result, which we denote by ℓ_j. This means that the log₂(4M²) bits of the binary string representing ℓ_j compose the integer part and the last ν bits compose the fractional parts of the approximation to λ_j. Then

  $$\begin{eqnarray*} &&|{\lambda}_j - \ell_j| \leqslant 17 \times 2^{-\nu} M^2. \end{eqnarray*}$$

  For details of the error estimate see proposition B.2. When d = 1, n (the number of qubits in register C) and ν are related by n = ν + log₂(4M²). Moreover, in the one-dimensional case $\hat {\lambda }_j = \ell _j$.
• (vi)  
  Let k_j be the binary string representing ℓ_j. For a fixed t, we implement the transformation

  $$\begin{equation} \underbrace{|{k_j}\rangle}_{n \,\, \rm qubits} |{0}\rangle^{\otimes n} \rightarrow |{k_j}\rangle \underbrace{|{k_j 2^t}\rangle}_{n \,\, \rm qubits}. \end{equation} \tag{ 30 }$$

  This is accomplished using CNOTs with the circuit shown in figure 9, since t ⩽ n the total number of quantum operations and qubits required to implement the circuit for all the values of t is O(n²).
• (vii)  
  Finally, we use phase kickback (see e.g. [43]) to obtain e^2πiϕ_j2t from the |k_j2^t〉 state where ϕ_j is the phase corresponding to the eigenvalue ℓ_j that approximates λ_j; see (18).

Zoom In Zoom Out Reset image size

To implement e^−iΔ_hγ, γ = 2π2^t/E, E defined in (16) and t = 0,...,n − 1 we use

$$\begin{equation} \mathrm{e}^{-\mathrm{i}\Delta_h \gamma}=\underbrace{\mathrm{e}^{\mathrm{i} h^{-2} L_h \gamma} \otimes \cdots \otimes \mathrm{e}^{\mathrm{i} h^{-2} L_h \gamma}}_{d \,\, \mathrm{matrices}}. \end{equation} \tag{ 31 }$$

Therefore the quantum circuit implementing e^−iΔ_hγ in d dimensions is obtained by the replication and parallel application of the circuit simulating e^ih−2L_hγ. For example, when d = 2 we have the circuit in figure 6. The register B of figure 2 contains dm qubits, m = log₂M and its initial state is assumed to have the form

$$\begin{equation} { \underbrace{(0,\ldots ,0}_{M^d-(M-1)^d},\underbrace{b_1,b_2,\ldots ,b_{(M-1)^d})}_{ \begin{array}{@{}c@{}} {\scriptsize \textrm{ values on the nodes of}} \\ \scriptsize \textrm{$(M-1)^{(\times{d})}$ grid} \end{array} }, } \end{equation} \tag{ 32 }$$

where $b_i = \left <i|f_h\right >$. This way we select the S_M−1 block in T^†_MF_2MT_M in (24) in each circuit for e^ih−2L_hγ. Recall that |f_h〉 corresponds to the right-hand side of (14).

The eigenvalues in the d-dimensional case are given as sums of the one-dimensional eigenvalues. We do not need to form the sums explicitly for the simulation of −Δ_h; they are computed by the tensor products. The difference between the d-dimensional and one-dimensional cases is that the register C in figure 2 has ⌈log₂ d⌉ additional qubits; i.e n = ⌈log₂ d⌉ + log₂ 4M² + ν. Accordingly, we generate the one-dimensional approximations to the eigenvalues using steps 1–5 of the eigenvalue estimation algorithm of the previous section. Then we append ⌈log₂ d⌉ qubits initialized to |0〉^{⊗⌈log₂ d⌉} to the left of the register holding the |ℓ_j〉 and carry out the remaining two steps, 6 and 7, with n = ⌈log₂ d⌉ + log₂ 4M² + ν. The error in the approximate eigenvalues is equal to 17M²d/2^ν; see theorem B.2.

Simulating the sine and cosine transforms (24) requires O(m²), m = log₂ M quantum operations and O(m) qubits [38]. The diagonal eigenvalue matrix of the one-dimensional case (23) is simulated by ESA. Its steps 1–3 and 5 require O(s²) quantum operations and O(s) qubits. In step 4 repeated squaring is performed ν + 7 times. Each repetition or step of the procedure requires O(s²) quantum operations and O(s) qubits. The total cost of step 4 is proportional to ν·O(s²) quantum operations and ν·O(s) qubits, accounting for any ancilla qubits used in repeated squaring. Step 6 requires O(n + t) quantum operations and qubits for fixed t. Step 7 requires O(n²) quantum operations, due to the Fourier transform, and O(n) qubits.

Using theorem B.2, and requiring error ε in the approximation of the eigenvalues, we have

$$\begin{eqnarray*} &&\frac{17E}{2^\nu}\leqslant \varepsilon, \end{eqnarray*}$$

$$\begin{eqnarray*} &&\nu = \left\lceil\log_2 \frac{17 E}{\varepsilon} \right\rceil, \end{eqnarray*}$$

i.e. ν = Θ(log₂ d + m + log₂ ε⁻¹). We also have n = Θ(ν) and s = Θ(n).

We derive the simulation cost taking the following facts into account.

• Steps 1–5 deal with the approximation of the eigenvalues. These computations are not repeated for every $t=0,\ldots ,n-1$. The total cost of these steps is O(n³) quantum operations and O(n²) qubits.
• The total cost of step 6, resulting from all the values of t, is O(n²) quantum operations and qubits.
• The total cost of step 7, that applies phase kickback for all values of t, does not exceed O(n³) quantum operations and O(n²) qubits.

Therefore the total cost to simulate e^ih−2L_hγ, γ = 2π2^t/E, for all $t=0,\ldots ,n-1$, is O(n³) quantum operations and O(n²) qubits. From (22) we conclude that the cost to simulate Poisson's matrix for the d-dimensional problem is d·O(n³) quantum operations and d·O(n²) qubits.

Finally, we remark that the dominant component of the cost is the one resulting from the approximation of the eigenvalues (i.e. the cost of steps 1–5).