A Query-based Quantum Eigensolver

First, we briefly review how to implement the fixed-point quantum search. The advantage of fixed-point search over the non-fixed-point search is the asymptotic convergence of the system state to the target state after a sequence of fixed-point Grover iterations Grover (2005). Just like Grover’s original search design, the fixed-point search can also achieve the quadratic speedup in query complexity Yoder et al. (2014). Specifically, in an oracle-based search problem, there is a classical oracle function $f$ on the total set $A$ satisfying $f(x)=1$ if and only if $x\in Q\subset A$, with $|A|=N$ and $|Q|=M\leq N$. The goal is to find any of the $M$ solutions in $Q$. In standard Grover’s search algorithm, one can construct an $N$-dimensional quantum system to encode the classical data $x\in A$ into its basis vectors $\{|x\rangle\}$. The initial state is chosen as $|\varphi\rangle=\frac{1}{\sqrt{N}}\sum_{x\in A}|x\rangle$, and the target state as $|t\rangle=\frac{1}{\sqrt{M}}\sum_{x\in Q}|x\rangle$. By appending an ancilla qubit $|b\rangle$, the quantum oracle $O_{f}$ can be expressed as:

Based on $O_{f}$, we can construct the parametrized inversion operators $R_{t}(\beta_{k})$ and $R_{\varphi}(\alpha_{k})$ with respect to $|\varphi\rangle$ and $|t\rangle$:

where $\alpha_{k}$ and $\beta_{k}$ are angle parameters. It can be shown that by choosing appropriate values of $\alpha_{k}$ and $\beta_{k}$, $k=1,\cdots,l$, the following Grover iteration sequence

will make the final state converge to target state $|t\rangle$ asymptotically as $l$ increases, satisfying $\left|\langle t|U^{(l)}|\varphi\rangle\right|^{2}\geq 1-\delta$, where $G(\alpha_{j},\beta_{j})=-R_{t}(\beta_{j})R_{\varphi}(\alpha_{j})$ and $l$ is the number of iteration. When $\alpha=\beta=\pi$, $G$ is reduced to the original non-fixed-point Grover iteration; when $\alpha=\beta=\frac{\pi}{3}$ (or $-\frac{\pi}{3}$ at some positions of $R_{\varphi}$ and $R_{t}$ in the oracle sequence), $G$ is reduced to Grover’s $\pi/3$ fixed-point algorithm; when the value of $\{\alpha_{k},\beta_{k}\}$ is chosen according to the results by Yoder et al Yoder et al. (2014), $G$ is reduced to Yoder-Luo-Chuang(YLC)’s fixed-point algorithm, with a quadratic speedup. Since each $R_{t}(\beta_{k})$ contains two quantum oracles $O_{f}$,

Define $\mu=|\langle t|\varphi\rangle|^{2}=\frac{M}{N}$. Assuming that each iteration requires two queries, we find the following optimal value of $q$ which can achieve the error threshold $\delta$:

from which we can see that YLC’s method can achieve quadratic speedup: $q=O\big{(}\frac{1}{\sqrt{\mu}}\big{)}=O(\sqrt{N})$.

Next, we show how to use fixed-point Grover’s search to solve a Type II eigenvalue problem on a quantum processor. Given an $N\times N$ Hermitian matrix $A$, we assume it has the spectral decomposition: $A=\sum_{m}\lambda_{m}|u_{m}\rangle\langle u_{m}|$, with $A|u_{m}\rangle=\lambda_{m}|u_{m}\rangle$, which is unknown to us before the calculation. To solve it on a quantum processor, we consider $A$ as a Hermitian operator on an $n$-qubit system, with $n\equiv\left\lceil{\log N}\right\rceil$. Without loss of generality, we assume $n=\log N$ to get rid of the cumbersome ceiling notations.

The target Type II problem is to find the eigenvalues near a given point $\lambda_{0}\in\mathbb{R}$, and the corresponding eigenvectors. Mathematically, the problem can be formulated as searching for all $\lambda_{m}$ within the $\epsilon$-neighborhood of $\lambda_{0}$: $B_{\epsilon}(\lambda_{0})\equiv[\lambda_{0}-\epsilon,\lambda_{0}+\epsilon]$. If more than one solution is located in $B_{\epsilon}(\lambda_{0})$, our method will randomly output one of the solutions at the final measurement; if there is no solution in $B_{\epsilon}(\lambda_{0})$, we can always tell this from the final measurement outcome. When the latter happens, we can enlarge the value of $\epsilon$ and redo the whole process for multiple times until we find a sufficiently large $\epsilon$ such that at least one solution is located in $B_{\epsilon}(\lambda_{0})$. Hence, without loss of generality, in the following we assume there is one and only one eigenvalue located in $B_{\epsilon}(\lambda_{0})$. We further assume it is $\lambda_{1}$. Then for quantum search, we can define the target state

where $M$ is the number of eigenvalues within $B_{\epsilon}(\lambda_{0})$. Under the unique-solution assumption, $|t\rangle=|u_{1}\rangle$.

The basic idea of solving a Type II eigenvalue problem through the quantum search is as follows. First we choose the initial state of the $n$-qubit quantum processor as a superposition of all eigenvectors of $A$: $|\varphi\rangle=\sum_{k}b_{k}|u_{k}\rangle$, such that the overlap $p\equiv|\langle t|\varphi\rangle|^{2}=|b_{1}|^{2}$ is nonzero (later we will show we can further guarantee that such $|\varphi\rangle$ with $p\geq\frac{1}{N}$, can be found through a fixed number of random selections). Then quantum search can be applied to amplify the overlap $p$ through a sequence of fixed-point Grover iterations until $p$ is sufficiently close to $1$. Finally a single measurement is enough to output the target eigenvector $|t\rangle=|u_{1}\rangle$ and the corresponding eigenvalue $\lambda_{1}$.

Specifically, we need to figure out how to design the quantum oracle $R_{t}(\beta)$ for this particular problem. The classical oracle $f$ should satisfy $f(u_{m})=1$ for $\lambda_{m}\in B_{\epsilon}(\lambda_{0})$, and $f(u_{m})=0$ for $\lambda_{m}\notin B_{\epsilon}(\lambda_{0})$. Hence, $R_{t}(\beta)$ should satisfy

The action of $R_{t}(\beta)$ is to add a relative phase $e^{i\beta}$ to all target solutions $|u_{m}\rangle$. We need a basic component in $R_{t}(\beta)$ whose input is an eigenvector $|u_{m}\rangle$, and whose output is the corresponding eigenvector $\lambda_{m}$, which determines whether to add the relative phase $e^{i\beta}$. A candidate gate to implement such component is a quantum phase estimation gate $U_{PE}=e^{2\pi iA}$. Indeed, $R_{t}(\beta)$ can be implemented by the circuit shown in Fig. 2, where $Z(\beta)=\text{diag}(1,e^{i\beta})$.

Hence, the total quantum processor is composed of three registers, the $r$-qubit eigenvalue register initialized as $|0\rangle$, the $n$-qubit eigenvector register initialized as $|\varphi\rangle$, and the ancilla single-qubit register initialized as $|1\rangle$. Then the total initial state becomes $|\Phi\rangle=|0\rangle|\varphi\rangle|1\rangle$ as in Fig. 2. In addition, due to the finite precision of the eigenvalue register, each eigenvalue $\lambda_{m}$ is denoted by a finite-precision state $|\tilde{\lambda}_{m}\rangle$. Define the conditional phase gate $CZ(\beta)$ on the eigenvalue register and the ancilla qubit to be:

Then the action of $R_{t}(\beta)\equiv U^{\dagger}_{PE}\cdot CZ(\beta)\cdot U_{PE}$ on $|\Phi\rangle$ gives:

where the approximate sign $\approx$ is due to the finite precision of $\tilde{\lambda}_{m}$ to represent $\lambda_{m}$. We can see that such construction of $R_{t}(\beta)$ is exactly what we need to implement the YLC’s fixed-point search method Yoder et al. (2014). On the other hand, given the initial state $\varphi$, there exists a unitary $V$ such that $|\varphi\rangle=V|0\rangle$. Then we can construct $R_{\varphi}(\alpha)\equiv VCZ^{\prime}(\alpha_{k})V^{\dagger}$, where $CZ^{\prime}(\alpha)$ is the controlled-phase gate acting on the eigenvector qubit and the ancilla, and has the form:

Based on $R_{\varphi}(\alpha)$ and $R_{t}(\beta)$, the fixed-point Grover iteration becomes $G=-R_{\varphi}(\alpha)R_{t}(\beta)$. Then we choose $\alpha_{k}$ and $\beta_{k}$ for each Grover iteration $G(\alpha_{k},\beta_{k})$ in the search sequence, according to the following rule Yoder et al. (2014):

where $\eta^{-1}=T_{1/L}(1/\sqrt{\delta})$, and $l$ is the number of iterations with $L=2l+1$. Here $T_{L}(x)=\cos[L\cos^{-1}(x)]$ is the $L^{th}$ Chebyshev polynomial of the first kind, and $\delta$ is the error requirement of the search results. So for error requirement $\delta$, we can obtain the optimal query number $q$:

where $q=2l$ represents the number of phase-estimation gates used in our proposal Yoder et al. (2014). If we assume, after $q$ number of fixed-point queries, the eigenvalue and the eigenvector registers becomes:

where the overlap $|c_{1}|^{2}$ between $|u_{1}\rangle$ and $|\psi\rangle$ is close to $1$, then applying one more QPE gate to $|\Psi\rangle$ will generate the final state

where $\tilde{\lambda}_{i}$ is the approximation of the eigenvalue $\lambda_{i}$ corresponding to the eigenstate $|u_{i}\rangle$. Hence, when we measure the first register, we can get the desired $|\tilde{\lambda}_{1}\rangle$, $|u_{1}\rangle$ with high probability. Thus, the Type II eigenvalue problem is solved. The algorithm flow chart is summarized in the following.

From Eqn. (6), we can see that the greater the overlap between the initial state $|\varphi\rangle$ and target state $|t\rangle$, the fewer number of queries $q$ is required. Hence, the choice of the initial state is crucial in order to keep $q$ small and keep the query-based method efficient. From the discussion in Appendix A, if the initial state $|\varphi\rangle$ is randomly chosen from the uniform distribution on the unit sphere of $\mathbb{C}^{N}$, with $|\varphi\rangle=\sum_{j=1}^{N}b_{j}|u_{j}\rangle$ and assuming the target state $|t\rangle=|u_{1}\rangle$, then the probability of the overlap $p=|b_{1}|^{2}=|\langle u_{1}|\varphi\rangle|^{2}$ to be larger than $\frac{1}{N}$ is given by

This implies that we can find the minimum $K=11$ such that $1-(1-1/e)^{K}\geq 0.99$. In other words, if we prepare a set of $11$ randomly-chosen initial states $\{|\varphi_{k}\rangle\}$, then the probability of at least one of $|\varphi_{k}\rangle$’s satisfying $|\langle\varphi_{k}|t\rangle|\geq 1/\sqrt{N}$ is larger than $0.99$. Such minimum value $K=11$ is independent of the size $N$ of the eigenvalue problem. For each $|\varphi_{k}\rangle$, if we implement the fixed-point query sequence with $q=\sqrt{N}\log(2/\sqrt{\delta})-1$ for $11$ times, then at least one of sequence after the final measurement will generate the final state $|\tilde{\lambda}_{1}\rangle|u_{1}\rangle$ with probability larger than $0.99$. Hence, the total query complexity of our method is $O(11q)=O(\sqrt{N})$.

In addition, as discussed in the appendix, the set of $11$ randomly-chosen initial states can be substituted by a set of $11$ computational basis states, which can be efficiently generated in the experiment. For example, if we can efficiently generate the ground state $|00\cdots 0\rangle$ of an $n$-qubit system, then all the other computational basis state can be efficiently prepared from the ground state with no more than $\log N=n$ number of bit-flips.

Analyzing the complexity of the quantum circuit in Fig. 2, we can see that the efficiency of the proposed query-based eigensolver algorithm depends on whether the unitary gate $U_{PE}$(or $U=e^{2\pi iA}$ in particular) can be efficiently generated. We can identify two cases when this is valid: (1) if $A$ is a $d$-sparse matrix Berry et al. (2007); (2) if $A$ is an inherent Hamiltonian that can be generated directly and efficiently on the quantum processor. In these cases, the complexity of generating $U_{PE}$ becomes $O(\operatorname{poly}(\log N))$. As we can prepare the initial state $|\varphi\rangle$ such that $p=|\langle u_{1}|\varphi\rangle|^{2}\geq\frac{1}{N}$, the query complexity of our method is $q=O(\sqrt{N})$. Since the total gate complexity is equal to the total query complexity multiplied by the oracle gate complexity, i.e. the complexity of $U_{PE}$, we have the total complexity of our algorithm to be $O(\sqrt{N}\text{poly}(\log N))$. Compared with the complexity of using QPE eigensolver to solve a Type II problem, our query-based eigensolver obtains a quadratic speedup.

The Heisenberg model is a well-known and useful testbed to study many-body physics and quantum information. The Hamiltonian of an $n$-qubit Heisenberg model is

where $J_{x},J_{y},J_{z}$ are coupling constants and $h$ indicates the external magnetic field. $\sigma^{k}$, $k=x,y,z$ are Pauli operators, and $\sigma_{j}^{k}$ denotes $\sigma^{k}$ acting on the $j$-th qubit. With further assumption of the periodic boundary condition, we assume $\sigma_{n+1}^{k}=\sigma_{1}^{k}$. Given $H$ and $\lambda_{0}$, our task is to find the eigenvalue $\lambda_{1}$ of $H$ near $\lambda_{0}$ and the corresponding eigenstate $|u_{1}\rangle$. In order to facilitate the simulation, we reconstruct a new Hamiltonian $\tilde{H}=H-\lambda_{0}I$, and then the original task is equivalent to solving the Type II eigenvalue problem for $\tilde{H}$ near the point $\tilde{\lambda}_{0}=0$. As mentioned earlier in the paper, we can use the query-based quantum eigensolver to solve this problem, using the quantum circuit shown in Fig. 2.

In the following, we present the simulation results for two cases of the Heisenberg model with $n=4,5$. For each case, we randomly choose the parameters $J_{x},J_{y},J_{z},h$ from the uniform distribution on $[0,1]$. In order to verify the efficiency of our proposed initial-state preparation method, we apply two methods to generate the initial state $|\psi\rangle$: (1) randomly choosing $11$ computational basis states $|x_{k}\rangle$, $x_{k}\in\{1,\cdots,N\}$, (2) randomly generating $11$ states $|y_{k}\rangle$ from the uniform distribution on the unit sphere in $\mathbb{C}^{n}$. According to the previous discussion, there exists at least one out of the 11 initial states satisfying the overlap $p=b_{1}\geq\frac{1}{2^{n}}$ with a high probability($\geq 0.99$), assuming the target state is $|u_{1}\rangle$, and $|b_{1}|^{2}=|\langle u|_{1}\varphi\rangle|^{2}$. Substituting $p=\frac{1}{2^{n}}$ and error threshold $\delta=0.01$ into Eqn. (6), we obtain the minimum values of the query number $q_{\min}=11$ for $n=4$, and $q_{\min}=16$ for $n=5$. In the circuit design, for each case of $n=4$ and $n=5$, our query-based eigensolver circuit encloses $4$ or $5$ qubits as the eigenvector register, $7$ qubits as the phase register and one more qubit as the ancilla. Next we implement our the query-based eigensolver $11$ times for each initial state $|\psi\rangle=|x_{k}\rangle$, or $|\psi\rangle=|y_{k}\rangle$, $k=1,\cdots,11$.

Table 1 and 2 summarize the simulation results for $n=4$ and $n=5$. For $n=4$, after $q=11$ fixed-point Grover’s queries, there are four $|x_{k}\rangle$’s satisfying $p_{1}\geq 1/16$ with final fidelity $F\geq 0.9933$, and three $|y_{k}\rangle$’s satisfying $p_{2}\geq 1/16$ with $F\geq 0.9892$(Table 1). Analogously, for $n=5$, after $q=16$ fixed-point Grover’s queries, there are four $|x_{k}\rangle$ satisfying $p_{1}\geq 1/32$, with $F_{1}\geq 0.9937$, and six $|y_{k}\rangle$ satisfying $p_{2}\geq 1/32$, with $F_{2}\geq 0.9724$(Table 2). Thus, our query-based method is able to solve the Type II eigenvalue problem for the Heisenberg model, within the expected accuracy and computational complexity.

To illustrate our proposed algorithm, we simulate the hydrogen molecule to find an eigenvalue near a given point and its corresponding eigenstate. After Born-Oppenheimer approximation, the Hamiltonian of the hydrogen molecule is as follows:

where the coefficients $h_{ij}$ and $h_{ijkl}$ are one- and two-electron overlap integrals Aspuru-Guzik et al. (2005); Whitfield et al. (2011); McWeeny (1992). The operators $a^{\dagger}_{j}$ and $a_{j}$ are the creation operator and the annihilation operators, respectively. After the Jordan-Wigner transformation Whitfield et al. (2011), the Hamiltonian $H$ becomes

We aim to solve the eigenvalue problem for $H_{JW}$ near a given point $\lambda_{0}=-0.8837$. For convenience, we reconstruct a new Hamiltonian $\tilde{H}=H_{JW}-\lambda_{0}I$. In order to apply our method to $\tilde{H}$, we enclose $4$ qubits in the eigenvector register, $7$ qubits in the phase register and one more qubit as the ancilla. With the error threshold chosen as $\delta=0.01$, we find the minimum query number required is $q_{\min}=11$. Analogous to the Heisenberg model, we randomly choose $11$ $|x_{k}\rangle$’s and $11$ $|y_{k}\rangle$’s as the initial states, which guarantees that at least one of them satisfies $p\geq\frac{1}{2^{n}}=1/16$. Indeed, simulation results demonstrate that, after $q=11$ queries, there are two $|x_{k}\rangle$’s satisfying $p\geq 1/16$, with $F\geq 0.9917$, and there are four $|y_{k}\rangle$’s satisfying $p\geq 1/16$ with $F\geq 0.9870$ (as shown in Table 3), in line with our expectation.