Direct Application of the Phase Estimation Algorithm to Find the Eigenvalues of the Hamiltonians

TL;DR

This paper introduces a method to directly apply the phase estimation algorithm to Hamiltonians without needing matrix exponential, using an ancilla framework and successive applications to efficiently estimate eigenvalues.

Contribution

The authors propose a novel approach that allows direct use of Hamiltonians in phase estimation, avoiding the costly matrix exponential computation and analyzing its complexity.

Findings

Efficient eigenvalue estimation for Hamiltonians composed of simple unitaries.

Reduction in quantum operations needed for phase estimation.

Successful simulation of hydrogen molecule ground state energy.

Abstract

The eigenvalue of a Hamiltonian, $\mathcal{H}$, can be estimated through the phase estimation algorithm given the matrix exponential of the Hamiltonian, $exp(-i\mathcal{H})$. The difficulty of this exponentiation impedes the applications of the phase estimation algorithm particularly when $\mathcal{H}$ is composed of non-commuting terms. In this paper, we present a method to use the Hamiltonian matrix directly in the phase estimation algorithm by using an ancilla based framework: In this framework, we also show how to find the power of the Hamiltonian matrix-which is necessary in the phase estimation algorithm-through the successive applications. This may eliminate the necessity of matrix exponential for the phase estimation algorithm and therefore provide an efficient way to estimate the eigenvalues of particular Hamiltonians. The classical and quantum algorithmic complexities of the framework are analyzed for the Hamiltonians which can be written as a sum of simple unitary matrices and shown that a Hamiltonian of order $2^n$ written as a sum of $L$ number of simple terms can be used in the phase estimation algorithm with $(n+1+logL)$ number of qubits and $O(2^anL)$ number of quantum operations, where $a$ is the number of iterations in the phase estimation. In addition, we use the Hamiltonian of the hydrogen molecule as an example system and present the simulation results for finding its ground state energy.

Tables

Table 1. Table 1: The values of the coefficients in Eq.( 37 ).

Coefficient	Value	Coefficient	Value
${\beta }_1$	-0.8126	${\beta }_9$	0.1659
${\beta }_2$	0.1712	${\beta }_{10}$	0.1205
${\beta }_3$	0.1712	${\beta }_{11}$	0.1743
${\beta }_4$	-0.2228	${\beta }_{12}$	-0.0453
${\beta }_5$	-0.2228	${\beta }_{13}$	0.0453
${\beta }_6$	0.1686	${\beta }_{14}$	0.0453
${\beta }_7$	0.1205	${\beta }_{15}$	-0.0453
${\beta }_8$	0.1659	$\kappa$	20.117

Equations

$$\tilde{\mathcal{H}}=\left(\mathbf{I}-\frac{i\mathcal{H}}{\kappa }\right)$$

$${\tilde{\lambda }}_j=1-i\frac{{\lambda }_j}{\kappa }=\left|1+i\frac{{\lambda }_j}{\kappa }\right|e^{i\frac{{\lambda }_j}{\kappa }}.$$

$${\mathbf{U}}^{(1)}=\left(\begin{array}{cc}\tilde{\mathcal{H}}& \bullet \\ \bullet & \bullet \end{array}\right),$$

$${\mathbf{U}}^{(1)}|0⟩|\phi ⟩=\sqrt{p_0}|0⟩\tilde{\mathcal{H}}|\phi ⟩+\sqrt{1-p_0}|\phi ⟩.$$

$${\tilde{\mathbf{U}}}^{(1)}\approx \left(\begin{array}{cc}\tilde{\mathcal{H}}& 0\\ 0& \bullet \end{array}\right),$$

$$\underset{k\text{ times}}{\underbrace{{\tilde{\mathbf{U}}}^{(1)}\dots {\tilde{\mathbf{U}}}^{(1)}}}\approx \left(\begin{array}{cc}{\tilde{\mathcal{H}}}^k& 0\\ 0& \bullet \end{array}\right).$$

$${\mathbf{I}}^{\otimes 1}\otimes {\mathbf{U}}^{(1)}=\left(\begin{array}{cccc}\tilde{\mathcal{H}}& \bullet & 0& 0\\ \bullet & \bullet & 0& 0\\ 0& 0& \tilde{\mathcal{H}}& \bullet \\ 0& 0& \bullet & \bullet \end{array}\right).$$

$$\left({\mathbf{I}}^{\otimes 1}\otimes {\mathbf{U}}^{(1)}\right)\left(\begin{array}{cccc}\mathbf{I}& 0& 0& 0\\ 0& 0& \mathbf{I}& 0\\ 0& \mathbf{I}& 0& 0\\ 0& 0& 0& \mathbf{I}\end{array}\right)\left({\mathbf{I}}^{\otimes 1}\otimes {\mathbf{U}}^{(1)}\right)=\left(\begin{array}{cccc}{\tilde{\mathcal{H}}}^2& \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \end{array}\right),$$

$$|{\psi }_0⟩=|0⟩|0⟩|{\phi }_j⟩,$$

$$|{\psi }_1⟩=\frac{1}{\sqrt{2}}\left(|0⟩|0⟩|{\phi }_j⟩+|1⟩|0⟩|{\phi }_j⟩\right).$$

$$|{\psi }_2⟩=\frac{1}{\sqrt{2}}\left(|0⟩|0⟩|{\phi }_j⟩+|1⟩{\mathbf{U}}^{(2^{\mathbf{k}})}|0⟩|{\phi }_j⟩\right).$$

$$\begin{array}{rl}|{\psi }_2⟩=& \frac{1}{\sqrt{2}}\left(|0⟩|0⟩|{\phi }_j⟩\right)\\ & +\frac{|1⟩}{\sqrt{2}}\left(|0⟩\hat{p}{e}^{i2\pi {\varphi }_j{2}^k}|{\phi }_j⟩+\sqrt{1-{\hat{p}}^2}|\Phi ⟩\right),\end{array}$$

$$|{\psi }_3⟩=\frac{1}{\sqrt{2}}\left(|0⟩|0⟩|{\phi }_j⟩+\hat{p}|1⟩|0⟩e^{i2\pi {\varphi }_j{2}^k}|{\phi }_j⟩\right).$$

$$|{\psi }_3⟩=\frac{1}{\sqrt{2}}\left(|0⟩+\hat{p}|1⟩e^{i2\pi (0.x_k)}\right)|0⟩|{\phi }_j⟩,$$

$$\begin{array}{rl}|{\psi }_{final}⟩=& \frac{1}{2}\left(1+\hat{p}{e}^{i2\pi (0.x_k)}\right)|0⟩|0⟩|{\phi }_j⟩\\ & +\frac{1}{2}\left(1-\hat{p}{e}^{i2\pi (0.x_k)}\right)|1⟩|0⟩|{\phi }_j⟩\end{array}$$

$$\frac{1}{2}\left(\left(1+\hat{p}\right)|0⟩+\left(1-\hat{p}\right)|1⟩\right)|0⟩|{\phi }_j⟩$$

$$\frac{1}{2}\left(\left(1-\hat{p}\right)|0⟩+\left(1+\hat{p}\right)|1⟩\right)|0⟩|{\phi }_j⟩$$

$$\mathbf{A}|0⟩=\sum _{x\in X_{good}}{\alpha }_x|\mathbf{x}⟩|\mathbf{\Phi }⟩+\sum _{x\in X_{bad}}{\alpha }_x|\mathbf{x}⟩|\mathbf{\Phi }⟩,$$

$$\mathbf{Q}=\mathbf{A}{{\mathbf{U}}_0}^{\perp }{\mathbf{A}}^{†}{\mathbf{U}}_{\mathbf{f}}.$$

$$\mathbf{Q}={\mathbf{U}}^{(1)}\left({{\mathbf{U}}_0}^{\perp }\otimes {\mathbf{I}}^{\otimes n}\right){{\mathbf{U}}^{(1)}}^{†}\left({{\mathbf{U}}_0}^{\perp }\otimes {\mathbf{I}}^{\otimes n}\right),$$

$$\Vert \tilde{\mathcal{H}}{\tilde{\mathcal{H}}}^{†}-\mathbf{I}\Vert =\Vert {\tilde{\mathcal{H}}}^2/{\kappa }^2\Vert .$$

$$\mathcal{H}=\sum _{l=1}^L{\alpha }_l{\mathbf{H}}_{\mathbf{l}},$$

$${\mathbf{U}}_{\mathbf{r}}\approx \sum _{k=0}^K\frac{1}{k!}(\frac{-i\mathcal{H}t}{r}{)}^k.$$

$${\mathbf{U}}_{\mathbf{r}}\approx \sum _{k=0}^K\frac{(-it/r{)}^k}{k!}\sum _{l1,\dots ,l_k=1}^L{\alpha }_{l_1}{\mathbf{H}}_{{\mathbf{l}}_1}\dots {\alpha }_{l_k}{\mathbf{H}}_{{\mathbf{l}}_{\mathbf{k}}}=\sum _{j=0}^M{\beta }_j{\mathbf{V}}_{\mathbf{j}},$$

$$select(\mathbf{V})|j⟩|\psi ⟩=|j⟩{\mathbf{V}}_{\mathbf{j}}|\psi ⟩.$$

$$\mathbf{B}|0⟩=\sqrt{p_0}\sum _{j=0}^M\sqrt{{\beta }_j}|j⟩;$$

$$\mathbf{U}=({\mathbf{B}}^{†}\otimes \mathbf{I})select(\mathbf{V})(\mathbf{B}\otimes \mathbf{I})=\left(\begin{array}{cc}{\mathbf{U}}_{\mathbf{r}}& \bullet \\ \bullet & \bullet \end{array}\right),$$

$$\tilde{\mathcal{H}}=\left(\mathbf{I}-\frac{i\mathcal{H}}{\kappa }\right)=\frac{1}{\kappa }\sum _{l=0}^L-i{\alpha }_l{\mathbf{H}}_{\mathbf{l}}=\sum _{l=0}^L{\beta }_l{\mathbf{V}}_{\mathbf{l}},$$

$${\tilde{\mathcal{H}}}^k=\frac{1}{{\kappa }^k}\sum _{l1,\dots ,l_k=0}^L(-i{)}^k{\alpha }_{l1}{\mathbf{H}}_{{\mathbf{l}}_1}\dots {\alpha }_{lk}{\mathbf{H}}_{{\mathbf{l}}_{\mathbf{k}}}.$$

$${\mathbf{U}}^{(k)}=\left(\begin{array}{cc}{\tilde{\mathcal{H}}}^k& \bullet \\ \bullet & \bullet \end{array}\right).$$

Full text

I Introduction

With the recent efforts such as digitalizing adiabatic quantum computers Barends et al. (2016), a digital quantum simulation of the real-time dynamic of a lattice gauge theoryMartinez et al. (2016), a simulation of the Hubbard modelSalfi et al. (2016), and the race to build universal quantum computers among the big companies Castelvecchi (2017), quantum computers become closer to be used in unsolved real-world applications.

As Feynman suggestedFeynman (1982), one of the breakthroughs of quantum computers is expected to be in the simulation of quantum systems, a known difficult problem for classical computers because of the exponential growth of the system size with the number of qubits (e.g., see Ref.Troyer and Wiese (2005)). In the standard formalism of quantum mechanics, the Hamiltonian of a quantum system is considered as a Hermitian matrix. Since the exponential of a skew-Hermitian matrix is a unitary matrix, the time evolution, $exp(-i\hbar\mathcal{H})$, of a quantum system represented by the Hamiltonian $\mathcal{H}$ describes a unitary transformation. In general, a computation or an algorithm should be represented by a unitary matrix in order to be implemented on a quantum computer. If an algorithm can be formalized in terms of unitary matrices, then a possible implementation may yield a computational efficiency over the classical algorithms: e.g. Shor’s integer factoring algorithm Shor (1994). Many quantum algorithms use the exponential $exp(-i\hbar\mathcal{H})$ to map the description of the original problem $\mathcal{H}$ to a unitary matrix. An exact mapping through matrix exponential generally requires eigendecomposition of $\mathcal{H}$: In spite of the existence of various numerical methods, the eigendecomposition in the cases where $\mathcal{H}$ is a large-dense matrix is still a big computational challenge for classical computers. Because of this computational difficulty, it is common to use some order of Trotter-Suzuki approximationSuzuki (1976) by writing $\mathcal{H}$ as a sum of terms whose exponentials are known or easy to compute. In general case, the order of the approximation affects the correctness and the high order approximations dramatically increases the computational complexity Poulin et al. (2015). Therefore, the complexity of this exponentiation should be also embodied in the consideration of the computational complexity of a quantum algorithm (see Ref.Aaronson (2015) for a relevant discussion of the topic).

It is known that using additional qubits eases the implementation difficulty of a quantum circuit. A quantum operator (or a dynamic of a system) can be implemented inside a larger system, where some reduced part of the system represents the action of the original system. This idea is used in different contexts such as designing circuits from matrix elements Daskin et al. (2012, 2014), and generating exponential of a matrix written as a sum of unitary matrices through Taylor series (Berry et al., 2015). Moreover, in Ref.Daskin and Kais (2017), to use a matrix on a quantum computer directly, we considered converting a non-unitary matrix, $\mathcal{H}$, to a unitary by using its square root. In a more general fashion, Ref.Low and Chuang (2016) described a qubitization approach to complete a non-unitary matrix to a unitary matrix whose multiple applications generate Chebyshev polynomials of $\mathcal{H}$.

In this paper, we follow a different approach: for a Hamiltonian $\mathcal{H}$ of order $2^{n}$ with real eigenvalues; we first put the matrix in the form $\widetilde{\mathcal{H}}=\left(\mathbf{I}-i\mathcal{H}/\kappa\right)$, where $\mathbf{I}$ is an identity matrix and $\kappa$ is a coefficient. For a sufficiently large $\kappa$; the value of the angle, the phase, seen in the polar form of the eigenvalue of $\widetilde{\mathcal{H}}$ becomes approximately equal to $\lambda_{j}/\kappa$, where $\lambda_{j}$ is an eigenvalue of $\mathcal{H}$. Using this fact and assuming $\mathcal{H}$ is a sum of some simple terms which can be easily mapped to quantum circuits; we present a phase estimation framework for finding $\lambda_{j}/\kappa$. In this framework, we use $\widetilde{\mathcal{H}}$ directly and show that necessary powers of $\widetilde{\mathcal{H}}$ can be estimated through successive applications after the application of the oblivious amplitude amplification. This provides a way to estimate the eigenvalues of $\mathcal{H}$ by using $(n+1+logL)$ qubits and $O(2^{a}nL)$ number of quantum operations. Here, $a$ is the number of iterations in the phase estimation and $L$ is the number of terms in the sum yielding $\mathcal{H}$. We also show that the powers of $\widetilde{\mathcal{H}}$ can be more accurately estimated by using additional qubits and a permutation operator in the phase estimation algorithm, which requires also the same number of quantum gates and $O(a+n+logL)$ number of qubits. As an example system, we use the Hamiltonian matrix of the hydrogen molecule and estimated its ground state energy within the presented framework.

In the following sections, after first describing the proposed framework, we will present the complexity analysis. Then, in Sec.IV, the simulation results for the hydrogen molecule are presented and discussed. And in the final section the paper is concluded with a summary.

II Proposed Eigenvalue Estimation

Phase estimation algorithm (PEA) estimates the phase $\phi$ for a given operator with the eigenvalue $e^{i\phi}$. For a Hamiltonian $\mathcal{H}\in C^{\otimes n}$, we propose to consider the following matrix in the phase estimation algorithm:

[TABLE]

where $\mathbf{I}$ is an identity matrix, and $\kappa$ is a coefficient. The eigenvalues of this matrix are in the form $\widetilde{\lambda}_{j}=(1-i\lambda_{j}/\kappa)$, where $\lambda_{j}$ is the $j$th eigenvalue of $\mathcal{H}$. When $\kappa$ is sufficiently large; in the polar form of $\widetilde{\lambda}_{j}$, the value of the angle becomes $\lambda_{j}/\kappa$ since $sin(\lambda_{j}/\kappa)\approx\lambda_{j}/\kappa$:

[TABLE]

Also note that when $\kappa$ is large, $\widetilde{\mathcal{H}}^{\kappa}$ gives an approximation to $e^{i\mathcal{H}}$. To be able to use the matrix in Eq.(1) directly in the phase estimation, we will initially assume that we have an efficient mechanism to produce the following unitary matrix:

[TABLE]

where the first part of the matrix is equal to $\widetilde{\mathcal{H}}$ with/without some normalization, and a $\bullet$ indicates another part of the matrix. As also mentioned in the introduction, these types of matrices are used in various contexts: either to design quantum circuits from matrix elements Daskin et al. (2012, 2014) or to be able to use it with the oblivious amplitude amplificationDaskin and Kais (2017), or to estimate unitary dynamic of a Hamiltonian through truncated Taylor series Berry et al. (2015) and its generalized form Low and Chuang (2016). In the next section we will discuss how the above matrix can be generated by using an ancilla register. Now let us consider the input state \left|\mathbf{0}\right\rangle$$\left|\varphi\right\rangle, where $\left|\varphi\right\rangle$ represents an abstract quantum state on the system register and $\left|\mathbf{0}\right\rangle$ is the first vector in the standard basis and the state on the ancilla register. The application of $\mathbf{U}^{(1)}$ to this input yields the following state Berry et al. (2015):

[TABLE]

The application of $\widetilde{\mathcal{H}}$ to the input $\left|\varphi\right\rangle$ on the system register can be obtained from the above output with the probability given by $p_{0}||\widetilde{\mathcal{H}}\left|\varphi\right\rangle||$, where $||.||$ represents a vector norm. For a unitary $\widetilde{\mathcal{H}}$, the probability simply becomes $p_{0}$. The desired output can be singled out through the application of a projector $\mathbf{P}=\left|\mathbf{0}\right\rangle\left\langle\mathbf{0}\right|$ to the ancilla register.

II.1 Estimation of $\widetilde{\mathcal{H}}^{2^{j}}$

For the phase estimation, we need to be able to efficiently generate the set of operators $\{\mathbf{U}^{(2^{0})}\dots\mathbf{U}^{(2^{a})}\}$ which consist of $\widetilde{\mathcal{H}}^{2^{0}}\dots\widetilde{\mathcal{H}}^{2^{a}}$. This can be done through the successive applications of $\mathbf{U}^{(1)}$ after the application of the oblivious amplitude amplification described in Sec.II.3: The amplitude amplification process to maximize the probability of $\left|\mathbf{0}\right\rangle$ in the ancilla turns $\mathbf{U^{(1)}}$ into the following form:

[TABLE]

$k$ number of applications of the matrix $\mathbf{\widetilde{U}}^{(1)}$ can be used to obtain $k$th power of $\widetilde{\mathcal{H}}$:

[TABLE]

A more accurate estimation for the $k$th power of $\widetilde{\mathcal{H}}$ in $\mathbf{U}^{(k)}$ can be obtained from $\mathbf{U}^{(k/2)}$ by using an additional qubit in the ancilla and a permutation matrix: For instance, if $\widetilde{\mathcal{H}}^{2}$ is desired, we first add one more qubit to the ancilla. Then the matrix representation of the circuit becomes

[TABLE]

Using a permutation matrix, $\mathbf{\Pi}$, similar to the one below, we can obtain the square from a successive application of this matrix as shown below:

[TABLE]

where the dimension of the sub-matrices are assumed to conform to the multiplications. In more general sense, $\mathbf{U^{(k)}}=(\mathbf{I}^{\otimes logk}\otimes\mathbf{U^{(k/2)}})\mathbf{\Pi}(\mathbf{I}^{\otimes logk}\otimes\mathbf{U^{(k/2)}})$. Here, $\mathbf{\Pi}$ can be implemented through some swap operations on the ancilla register.

Consequently, we can use the set of operators $\{\mathbf{U}^{(2^{0})}\dots\mathbf{U}^{(2^{a})}\}$ generated by successive applications in the phase estimation algorithm and estimate the value of $\frac{\lambda_{j}}{\kappa}$ and hence the eigenvalue of the Hamiltonian with an accuracy affected by $a$. Because of the additional qubits, it may be easier to implement this method through the iterative version of the phase estimation algorithm as depicted in Fig.1.

II.2 The steps inside PEA

The phase estimation algorithm at the $x=(a-k-1)$th iteration is composed of three registers: vis., the phase with one-qubit, the ancilla with $x+logL$ qubits, and the system register with $n$ qubits. The initial input to the algorithm is given by the following:

[TABLE]

where $\left|\mathbf{\varphi}_{j}\right\rangle$ is an estimation to the $j$th eigenvector of $\mathcal{H}$. In PEA, we first apply the Hadamard gate to the phase qubit which simply puts this qubit into the superposition and yields the state:

[TABLE]

Then, $\mathbf{U^{(2^{k})}}$ controlled by the phase qubit is applied to the remaining qubits:

[TABLE]

With the help of Eq.(4), $\left|\mathbf{\psi}_{2}\right\rangle$ can be rewritten in the following form:

[TABLE]

where $\hat{p}=\sqrt{p_{0}}|\widetilde{\lambda}_{j}|^{2^{k}}$. The application of the operator $\mathbf{P}$ to the ancilla register causes the state to collapse into the following unnormalized state:

[TABLE]

In the iterative phase estimation, in each iteration we get the estimate $k$th bit value of the phase. This is achieved with the help of a $R_{z}$ gate whose rotation angle is determined from the previous iterations (please refer to Fig.1 for the value of the angle.). This gate converts the term $e^{i2\pi\phi 2^{k}}$ into the form $e^{i2\pi(0.x_{k})}$, where $x_{k}$ represents a binary bit: i.e., simply:

[TABLE]

In the final step, the Hadamard gate is applied again to the phase qubit:

[TABLE]

Measurement on the phase qubit in this unnormalized state yields the bit value of $x_{k}$:

• •

When $x_{k}=0$, the probability of $\left|0\right\rangle$ becomes greater in the final state:

[TABLE]

• •

In the case $x_{k}=1$, the probability of $\left|1\right\rangle$ becomes greater:

[TABLE]

Here, the probability difference between measuring 1 and 0 is determined by $\frac{2\hat{p}}{\mu}$, where $\mu$ is a normalization constant. Before the application of the projector $\mathbf{P}$, $\mu$ can also be eliminated by applying the oblivious amplitude amplification to $\mathbf{U}^{(1)}$ as described below.

II.3 Application of oblivious amplitude amplification

The amplitude amplification Mosca et al. (1998); Brassard et al. (2002); Kaye et al. (2006) is based on the Grover’s search algorithm Grover (1998) where one applies a sequence of the operators to increase the magnitude of the amplitudes of some desired states: Consider the following output state:

[TABLE]

where $\mathbf{A}$ is a quantum algorithm, $X_{good}$ and $X_{bad}$ are the sets of good (desired) states and bad(undesired) states of the first register, $x$ represents a standard basis vector, and $\left|\mathbf{\Phi}\right\rangle$ represents the states of the qubits in the second register. The probability of the good states in this output can be increased by the application of the iteration operator:

[TABLE]

Here, the operator $\mathbf{U_{f}}$ marks (multiply by -1) the amplitudes of the desired states and $\mathbf{U_{0}}^{\perp}=2\left|\mathbf{0}\right\rangle\left\langle\mathbf{0}\right|-\mathbf{I}$.

A version of the amplitude amplification called oblivious amplitude amplification Paetznick and Svore (2014); Berry et al. (2014, 2015) can be used to increase the probability of the part $\left({\sqrt{p_{0}}}\left|\mathbf{0}\right\rangle\widetilde{\mathcal{H}}\left|\varphi\right\rangle\right)$ in Eq.(3). In the oblivious amplitude amplification algorithm, $\mathbf{U_{f}}$ and $\mathbf{U_{0}}^{\perp}$ only operates on the ancilla and the probability change does not affect the state in the second register. In our case, we use the following iteration operator:

[TABLE]

where $\mathbf{U_{0}}^{\perp}$ acts on the ancilla register. The main difference from Eq.(19) is that Eq.(20) does not depend on the input to the system register and the marking and the amplifying operators, $\mathbf{U_{0}}^{\perp},$ are applied only on the ancilla register.

The oblivious amplitude amplification algorithm best works when $\widetilde{\mathcal{H}}$ is a unitary process. However, recently it is shown that the algorithm also works when $\widetilde{\mathcal{H}}$ is close to a unitary matrix Berry et al. (2015). In this case, the error in the output of one step of the oblivious amplification is bounded by the distance which in our case measured as: $\|\widetilde{\mathcal{H}}\widetilde{\mathcal{H}}^{\dagger}-\mathbf{I}\|$. Since $\widetilde{\mathcal{H}}\widetilde{\mathcal{H}}^{\dagger}=\mathbf{I}+\mathcal{H}^{2}/\kappa^{2}$, this distance becomes:

[TABLE]

Thus, the error in the output of the oblivious amplitude amplification is bounded by $O(\|\widetilde{\mathcal{H}}^{2}/\kappa^{2}\|)$. For a large $\kappa$, $\|\mathcal{H}^{2}/\kappa^{2}\|$ and so the order of the error can be expected to be very small. In addition, since $\widetilde{\mathcal{H}}^{\kappa}$ gives an approximation to $e^{i\mathcal{H}}$, in the iterations of PEA, $\widetilde{\mathcal{H}}^{2^{j}}$ for some $j$ can be expected to be nearly a unitary matrix. Therefore, in the case of $\widetilde{\mathcal{H}}$ with a large $\kappa$, it is possible to use oblivious amplitude amplification. To observe this, the simulation results in Sec.IV presented with and without the oblivious amplitude amplification.

II.4 Simulating sums of unitary matrices

Any Hamiltonian can be decomposed into a linear combination of unitary matrices:

[TABLE]

where $\alpha_{l}$ is some complex coefficient and $\mathbf{H_{l}}$ represents a unitary matrix.

Berry et al.Berry et al. (2015) have showed that the time evolution of $\mathcal{H}$ can be simulated as follows: The time is considered to be divided into $r$-segments so that the time evolution at segment $r$, $\mathbf{U_{r}}=e^{-i\mathcal{H}t/r}$, is approximated through the following Taylor expansion:

[TABLE]

Since $\mathcal{H}^{k}=\sum_{l1,\dots,l_{k}=1}^{L}\alpha_{l_{1}}\dots\alpha_{l_{k}}\mathbf{H_{l_{1}}}\dots\mathbf{H_{l_{k}}}$, this expansion is rewritten as:

[TABLE]

where $\mathbf{V_{j}}$ is in the form of $(-i)^{k}\mathbf{H_{l_{1}}}\dots\mathbf{H_{l_{k}}}$, $\beta_{j}$ is a complex coefficient, and $M$ is the number of terms in the summation and equal to $L^{k}$. Assuming a mechanism to implement each $\mathbf{V_{j}}$ exist, then $\mathbf{V_{j}}$s are combined in a quantum circuit $select(\mathbf{V})$ by using $m=logM$ ancilla qubits in a way that for each state $\left|\mathbf{j}\right\rangle$ on the ancilla, a $\mathbf{V_{j}}$ is applied to the main qubits. Here, $\left|\mathbf{j}\right\rangle$ represents the $j$th vector in the standard basis. For an arbitrary state $\left|\psi\right\rangle$, this can be shown as:

[TABLE]

The state on the ancilla can be put into a superposition state where each basis biased by the coefficients $\beta_{j}$s:

[TABLE]

where $p_{0}=\sum_{j=0}^{M}\beta_{j}$ and $\beta_{j}$s are assumed to be positive. The time evolution $\mathbf{U_{r}}$ can be obtained by combining the operator $select(\mathbf{V})$ with $\mathbf{B}$ in the following way:

[TABLE]

As similarly explained in the previous sections, this matrix can be used simulate action of $\mathbf{U_{r}}$ on some arbitrary state $\left|\psi\right\rangle$ with the probability determined by $p_{0}$.

When $\mathcal{H}$ given in the form of Eq.(22), the matrix described in Eq.(1) can also be represented as a sum of unitary matrices:

[TABLE]

where $\mathbf{H_{0}}=i\kappa\mathbf{I}$. Using the same methodology given in Eq.(27), we can form the matrix $\mathbf{U}^{(1)}$ in Eq.(3) easily.

Also note that since $\widetilde{\mathcal{H}}$ is a sum of unitary matrices, its powers are also similar to those of $\mathcal{H}$. For instance, its $k$th power reads as follows:

[TABLE]

Therefore,

[TABLE]

However, in terms of complexity, the required number of qubits is $klogL$ which grows linearly with the power. In the case of using this matrix inside PEA, the required number of qubits grows exponentially with the precision: That means if we need $\mathbf{U}^{2^{j}}$ at an iteration of PEA, we need to use $2^{j}logL$ qubits to represent $L^{2^{j}}$ number of terms in the matrix $select(\mathbf{V})$. Therefore, when the available number of qubits is limited, the power should be taken through successive applications of Eq.(5) or Eq.(8). In the case of Eq.(8), only one additional qubit is necessary to find the square of $\widetilde{\mathcal{H}}$ from the previous iteration of the phase estimation algorithm.

II.5 Simulating sum of rank one matrices

Any matrix can also be written as a sum of rank one matrices. Consider the matrices in the following form:

[TABLE]

where $\mathbf{H_{j}}=\left|\mathbf{x_{j}}\right\rangle\left\langle\mathbf{x_{j}}\right|$ is a rank one Hermitian matrix constructed by the normalized vector $\left|\mathbf{x_{j}}\right\rangle$. Such sums are frequently seen in the form $\mathbf{XX^{T}}$ in machine learning problems and statistical analysis of data points-where a column of the matrix $\mathbf{X}$ represents a data point. We can rewrite Eq.(31) as follows:

[TABLE]

where $\alpha_{j}$s are coefficients and $\mathbf{R_{0}}=\mathbf{I}$. The above equation describes $\mathcal{H}$ in terms of a sum of $(L+1)$ number of idempotent-unitary matrices: i.e., $\mathbf{R_{j}}^{2}=\mathbf{I}$. As a result, $\mathcal{H}$ can be simulated in the phase estimation after being turned into the form of Eq.(1). In the next section, we will analyze the complexity of the method in the case the Hamiltonian is given as a sum of unitary matrices.

III Computational Complexity

Computational complexity of a quantum algorithm involves two aspects: vis., classical and quantum complexities. The classical complexity covers the preprocessing time of circuits implementing the quantum algorithm, i.e. $e^{i\mathcal{H}}$, and post-processing of the output for obtaining a desired solution. In this paper, we have mainly focused on avoiding the necessity of $e^{i\mathcal{H}}$ in estimation of the eigenvalue of $\mathcal{H}$ on quantum computers. When the matrix $\mathcal{H}$ can be mapped to quantum circuits easily (e.g., it is a sum of ”simple” unitary matrices); then the classical complexity is bounded by the number of matrix elements: i.e., $O(N^{2})$, assuming all the matrix elements are needed to be processed and stored.

On the other hand, quantum complexity of a quantum algorithm is determined by the number of qubits and the number of one-and two-qubit quantum gates involved in the circuit implementing the quantum algorithm (Here, we only consider algorithmic complexity.). In this perspective, the complexity analysis of the method is given below.

III.1 The number of qubits

When the powers of $\widetilde{\mathcal{H}}$ is obtained through Eq.(5), the required number of qubits is $(n+1+logL)$. However, when obtained through Eq.(8), it becomes $(a+n+logL)$, where $a$ is the number of iterations in the phase estimation.

III.2 Circuit Implementation of $\mathbf{B}$

For a normalized $L$-dimensional column-vector $\left|\mathbf{x_{j}}\right\rangle\in C^{\otimes l}$ and the identity matrix $\mathbf{I}$, $\mathbf{R_{j}}=\mathbf{I}-2\left|\mathbf{x_{j}}\right\rangle\left\langle\mathbf{x_{j}}\right|$ is an Householder transformation describing a reflection operator around the vector $\left|\mathbf{x_{j}}\right\rangle$. On quantum computers, a Householder transformation Golub and Van Loan (2012) can be implemented by using $O(2^{l})$ total number of two- and one-qubit quantum gates Ivanov et al. (2006); Bullock et al. (2005); Urías and Quiñones (2015) ( In Ref.Ivanov et al. (2006), an implementation with the same complexity is also presented for a general version of the Householder transformation: i.e., $\mathbf{I}-(e^{i\varphi}-1)\left|\mathbf{x_{j}}\right\rangle\left\langle\mathbf{x_{j}}\right|$).

Note that a Householder matrix, $\mathbf{R}$, can be also implemented through $(2L-3)$ number of plane (Givens) rotationsGolub and Van Loan (2012) as follows: It is known that a Givens rotation $\mathbf{G_{j-1}}$ can be used to zero out $j$th entry in the column of a matrix. The Givens rotation $\mathbf{G_{L-1}}$ applied to $\mathbf{R}$ not only nullifies the last element in the first column of $\mathbf{R}$ but all the entries on the $L$th row upto the first sub-diagonal entry as well. Furthermore, the matrix product $(\mathbf{G_{L-1}RG_{L-1}^{T}})$ yields a matrix where all entries but the diagonal entry on the last row and column are zeros and the diagonal element is one. Because of this property, $\mathbf{R}$ can be shown as a product of $(2L-3)$ number of plane rotations: That is, $\mathbf{R}=\mathbf{G_{L-1}^{T}}\dots\mathbf{G{3}^{T}G_{2}^{T}G_{1}^{T}G_{2}G_{3}}\dots\mathbf{G_{L-1}}$. As a result, the circuit for $\mathbf{R}$ can be formed via plane rotations in this decomposition. Since a rotation matrix can be implemented through a multi controlled quantum gate Cybenko (2001); when considered together, the rotations in the product form a uniformly controlled gate-network. These networks are well-studied in Ref.Vartiainen et al. (2004) where it is shown that a uniformly controlled network acting on $l$ qubits can be simplified into $O(2^{l})$ number of one- and two-qubit gates. Also note that the sparsity of the vector $\left|\mathbf{x_{j}}\right\rangle$ directly affects the number of plane rotations in the decomposition and hence the number of gates required for the implementation.

III.3 Circuit Implementation of $select(\mathbf{V})$

The circuit for $select(\mathbf{V})$ is determined by the number of terms which is given by $L$ and the number of gates required to implement each term. If each term in the sum can be implemented through simple $O(n)$ quantum gates on different qubits, then the combination of the all terms form $O(n)$ different gray-coded networks controlled by $l$ number of qubits in the ancilla. The decomposition of this networks will require $O(nL)$ quantum gates.

Therefore, the total complexity to implement $\mathbf{U^{(1)}}=\mathbf{B}^{\dagger}\otimes\mathbf{I})select(\mathbf{V})(\mathbf{B}\otimes\mathbf{I})$ is bounded by $O(nL+2L)=O(nL)$ with the assumption that each term in the sum can be implemented with $O(n)$ quantum gates. In an iteration of the phase estimation, finding the power $\mathbf{U^{(2^{k})}}$ requires $2^{k}$ number of successive applications of $\mathbf{U^{(1)}}$. Therefore, the total required quantum gates for an iteration of the algorithm becomes $O(2^{k}nL)$, where $1\leq k\leq a$ with $a$ being the total number of iterations.

IV An Example Application to the Hamiltonian of $H_{2}$

The concept of the second quantization used in quantum chemistry to simplify the formalism of fermionic many particle systems. It represents the interacting systems of electrons and nuclei through the creation and annihilation operators. The molecular electronic Hamiltonian in electronic structure problem is expressed in the second quantization form as followsLanyon et al. (2010); Whitfield et al. (2011); Seeley et al. (2012):

[TABLE]

where $h_{pq}$ is the one-electron integrals including the electronic kinetic energy and the electron nuclear attraction terms. $h_{pqrs}$ represents the set of two-electron integrals with the electron-electron interactions. $a_{j}$ and $a_{j}^{\dagger}$ are the lowering and raising operators. This type of Hamiltonians can be represented in terms Pauli matrices by using the following Jordan-Wigner transforms:

[TABLE]

where

[TABLE]

Whitfield et al.Lanyon et al. (2010); Whitfield et al. (2011) and Seeley et al.Seeley et al. (2012) thoroughly studied this mapping and used the Hamiltonian for the hydrogen molecule as an example system: Using a minimal number of basis, only four spin orbitals indexed from [math] to $3$ are involved in the above sum. They found the values of the one- and two-electron integrals by using a restricted Hartree-Fock calculation at an internuclear separation of $7.414\times 10^{-11}m$. Because of the overlap in the integral values, the Hamiltonian is reduced to the following form:

[TABLE]

Using the mappings in Eq.(34), the Hamiltonian is rewritten into the terms of Pauli matrices with the values of the coefficients given in Table 1 (taken from Eq.(80) of Ref.Seeley et al. (2012)):

[TABLE]

The simulation of this Hamiltonian within the phase estimation algorithm is done after generating a circuit equvailent of $e^{i\mathcal{H}t}$ through the Trotter-Suzuki approximation. As mentioned before, the complexity of the generated circuit is determined by the order of the approximation which affects the accuracy of the obtained ground state energy of the Hamiltonian from the phase estimation algorithm.

In our case, we first convert the Hamiltonian into the form of Eq.(1) by using $\kappa=10\times||\mathcal{H}||_{1}$, which guaranties that $|\lambda_{j}/{\kappa}|\leq 0.1\approx sin(0.1)$. This results in total 16 terms for $select(\mathbf{V})$, where $\kappa$ is taken as the coefficient for the additional identity and the negative signs and $i$ are shifted to the terms inside $select(\mathbf{V})$ so that all $\beta_{j}$s become positive real numbers (This is because $\sqrt{\beta_{j}}$s are involved in the construction of $\mathbf{B}$).

In the circuit, an ancilla of 4 qubits is required to control each term separately to construct $select(\mathbf{V})$. As also mention in the previous section, in Ref.Möttönen et al. (2004), a circuit decomposition technique is described for circuits (called gray-coded networks) where each quantum gate on a qubit is controlled by $2^{l}$ distinct states of some other $l$ qubits. It is shown that, the network can be decomposed into $2^{l}$ number of $CNOT$ and $2^{l}$ number of single gates by using some form of Hadamard transformation. The circuit for $select(\mathbf{V})$ consists of 4 such networks (one network for per qubit in the system). Each network is controlled by $logL$ qubits in the ancilla. Since each decomposed network will require $L$ number of $CNOT$s and $L$ number of single quantum gates, the circuit complexity of $select(\mathbf{V})$ is bounded by $O(nL)$, where $n$ is the number of qubits in the system. As a result, the for circuit $H_{2}$ will require 64 $CNOT$ gates. In the previous section, the circuit complexity of $\mathbf{B}$ is given by $O(L)$, which makes the whole circuit complexity $\approx{100}$ gates. In the phase estimation, since we have also a phase qubit added to the control, the required number of gates for $select(\mathbf{V})$ will be $\approx 128$ gates.

The simulation is done in three different ways:

• •

In the first case, the powers are estimated through the successive application of $\mathbf{\widetilde{U}^{(1)}}$ given in Eq.(5). Here the amplitude amplification operator $\mathbf{Q}$ given in Eq.(20) is applied six times to obtain a maximum probability. The phase estimation algorithm is run for 25 iterations and the measured unnormalized probabilities for the phase qubit are drawn in Fig.2. The obtained phase from the bit values is 0.014603. From $\arcsin(\text{Im}(e^{-i2\pi 0.014603}))\times\kappa$, the value of the ground state energy is found as -1.845843. Since the system size is small, the ground state energy of the Hamiltonian has been also classically computed as $-1.851046$ via eigen-decomposition. This shows an error value $0.0052$ in our estimation caused by the error in the estimation of the powers of the matrices through Eq.(5) and the error in the approximation $sin(\lambda_{j}/\kappa)\approx\lambda_{j}/\kappa$.

• •

In the other cases, the power of the matrices is computed through Eq.(8) by adding one more qubit to the ancilla in each iteration. Here, because of the computational difficulty-the system size for 9 iteration requires 17 qubits, we present the output of the phase qubit for only 9 iterations of the phase estimation without and with the amplitude amplification applied to $\mathbf{U^{(1)}}$ respectively in Fig.3 and Fig.4. In comparison to classically computed value, the simulation results produce the correct eigenvalue if the previous bit values are provided; otherwise, it generates an estimate of the eigenvalue, as expected, with an error bounded by $O(\kappa/2^{-9})$. As seen in Fig.3, without the amplitude amplification, depending on the value of $k$ in $\mathbf{U}^{2^{k}}$, the probability difference between 1 and 0 on the phase qubit diminishes and the total success probability approaches to 0.5 (the probability of the part of the ancilla system simulating $\widetilde{\mathcal{H}}$). As shown in Fig.3 with the same $\kappa$ value, the total probability can be maximized by the application of the amplitude amplification only to $\mathbf{U^{(1)}}$. The probability difference now survives for 9 iterations. This can be further improved by using oblivious amplitude amplification algorithm in the subsequent iterations when the probability difference starts to drop and changing the value of $\kappa$. Note that a larger $\kappa$ will require more iterations to obtain a higher precision. Also note that bit values and the probabilities are affected by the number of iterations because the iterative phase estimation algorithm uses a rotation gate whose angle determined by the bit values measured in the previous iterations.

V Conclusion

In this paper, we have described a framework to estimate the eigenvalues of Hamiltonians in the phase estimation algorithm without using the time evolution operator, the exponential, of the Hamiltonian. We have shown how to find the powers of the matrices necessary in the phase estimation algorithm by successive applications. We have analyzed the circuit implementation of the whole framework in terms of classical and quantum complexities and showed that the framework provides an efficient way to estimate the eigenvalues of Hamiltonians which are written in terms of sums of simple unitary matrices (here, “simple” means a unitary matrix which can be implemented through a few number of one- and two-qubit quantum gates.). In addition, we have used the Hamiltonian of hydrogen molecule as an example system and showed how to estimate its ground state energy through the described method. We believe this framework can be used efficiently for many eigenvalue related problems such as the finding ground state energy of quantum systems.