Rodeo Algorithm for Quantum Computing

Quantum computing is a powerful paradigm with the potential to describe large complex systems and eventually perform computations beyond the reach of classical computing. Recently, there have been several exciting algorithmic advances in describing the time evolution of Hamiltonians on quantum computers using a variety of different tools Low and Chuang (2019, 2017); Childs et al. (2019); Campbell (2019); Childs et al. (2019); Meister et al. (2020). They can be broadly categorized as either Lie-Trotter-Suzuki product formulas Trotter (1959); Suzuki (1976) or linear combinations of unitaries Childs and Wiebe (2012). Unfortunately, the application of these techniques for quantum state preparation is limited by existing hardware capabilities. Quantum adiabatic evolution is one approach to quantum state preparation that starts with an eigenstate of a simple Hamiltonian that slowly evolves with an interpolating time-dependent Hamiltonian until reaching the desired target Hamiltonian Farhi et al. (2000, 2001). The problem is that calculations based on quantum adiabatic evolution require an extended time evolution that makes the computational cost prohibitive for large systems. To address this problem, we introduce a new framework for quantum state preparation and spectrum determination called the rodeo algorithm.

The rodeo algorithm employs a strategy that is opposite to quantum adiabatic evolution. As the name suggests, the rodeo algorithm operates by shaking off all other states until only the target eigenvector remains. In this regard, the rodeo algorithm is similar in character to the projected cooling algorithm Lee et al. (2020); Gustafson (2020). However, the rodeo algorithm has the advantage that it can be applied to any quantum Hamiltonian and is a recursive algorithm that achieves exponential convergence in the number of cycles. It can be used to compute the full energy spectrum as well as prepare any energy eigenstate. While the rodeo algorithm might appear similar to Kitaev’s iterative version of quantum phase estimation Kitaev (1995) and fixed-time energy band filtering methods Ge et al. (2017); Lu et al. (2020), none of these methods can be used efficiently to prepare individual eigenstates of a general quantum Hamiltonian.

We will refer to the Hamiltonian of interest as the object Hamiltonian, $H_{\rm obj}$, and the linear space which it acts upon the object system. By assumption, the object system starts in some initial state $\ket{\psi_{I}}$, which in general will update after each measurement. We will use auxiliary or ancilla qubits coupled to the object system. In the following we use the standard terminology, ancilla qubits. But we also mention that this collection of ancilla qubits is also informally called the “rodeo arena”.

If the quantum device allows for mid-circuit measurements, then only one ancilla qubit is needed. However, here we focus on the implementation using different ancilla qubits for each cycle of the rodeo algorithm. Each of the ancilla qubits is initialized in the state $\ket{1}$ and operated on by a Hadamard gate H. We use each ancilla qubit $n=1,\cdots,N$ to control the time evolution of $H_{\rm obj}$ for time $t_{n}$. In order to achieve the desired energy filtering, we operate on each ancilla qubit $n$ with the phase rotation gate P$(Et_{n})$, follow that with another Hadamard gate H, and then measure the qubit. We use the convention that P$(Et_{n})$ multiplies the phase $e^{iEt_{n}}$ to the $\ket{1}$ state and leaves the $\ket{0}$ state untouched. The circuit diagram for the rodeo algorithm is shown in Fig. 1.

In order to illustrate the effect of these gate operations, let us explicitly write out the operation for one cycle of the rodeo algorithm with one ancilla qubit. Starting from the initial state $\ket{1}\otimes\ket{\psi_{I}}$ and performing one rodeo cycle, we obtain

	$\displaystyle\begin{bmatrix}\left[\tfrac{I}{2}-\tfrac{1}{2}e^{-i(H_{\rm obj}-E)t_{1}}\right]\ket{\psi_{I}}\\ \left[\tfrac{I}{2}+\tfrac{1}{2}e^{-i(H_{\rm obj}-E)t_{1}}\right]\ket{\psi_{I}}\end{bmatrix}=$
	$\displaystyle\begin{bmatrix}\tfrac{I}{\sqrt{2}}&\tfrac{I}{\sqrt{2}}\\ \tfrac{I}{\sqrt{2}}&\tfrac{-I}{\sqrt{2}}\end{bmatrix}\begin{bmatrix}I&0\\ 0&Ie^{iEt_{1}}\end{bmatrix}\begin{bmatrix}I&0\\ 0&e^{-iH_{\rm obj}t_{1}}\end{bmatrix}\begin{bmatrix}\tfrac{I}{\sqrt{2}}&\tfrac{I}{\sqrt{2}}\\ \tfrac{I}{\sqrt{2}}&\tfrac{-I}{\sqrt{2}}\end{bmatrix}\begin{bmatrix}0\\ \ket{\psi_{I}}\end{bmatrix},$		(1)

where $I$ is the identity operator on the object system. We note that $H_{\rm obj}$ commutes with all of our gates, and so we can describe the action of the rodeo algorithm for each individual eigenvector of $H_{\rm obj}$ with energy $E_{\rm obj}$. In that case, the probability of measuring the ancilla qubit $n$ in the $\ket{1}$ state is

$$\cos^{2}\left[(E_{\rm obj}-E)\tfrac{t_{n}}{2}\right]=\left|\tfrac{1}{2}+\tfrac{1}{2}e^{-i(E_{\rm obj}-E)t_{n}}\right|^{2}.$$

(2)

The success probability of measuring all $N$ ancilla qubits in the $\ket{1}$ state is given by product,

$$P_{N}=\prod_{n=1}^{N}\cos^{2}\left[(E_{\rm obj}-E)\tfrac{t_{n}}{2}\right].$$

(3)

If we now take random values of $t_{n}$, we have an energy filter for $E_{\rm obj}=E$. The geometric mean of $\cos^{2}\theta$ when sampled uniformly over all $\theta$ is equal to $\tfrac{1}{4}$. Therefore the spectral weight for any $E_{\rm obj}\neq E$ is suppressed by a factor of $\tfrac{1}{4^{N}}$ for large $N$.

In Fig. 2 we plot the probability $P_{N}$ of measuring the $\ket{1}$ state for all ancilla qubits versus $E_{\rm obj}-E$ for $1$ (dotted black line), $4$ (dashed blue line), and $8$ (solid red line) cycles with Gaussian random values of $t_{n}$ with root-mean-square value $t_{\rm RMS}=10$. If we use a Gaussian approximation for $P_{N}$ near its maximum value of $1$ at $E_{\rm obj}=E$, we find that the width of the peak scales as $O[1/(\sqrt{N}t_{\rm RMS})]$.

Let $\epsilon$ be the desired energy resolution of our rodeo algorithm such that all energy eigenvectors outside of the interval $[E-\epsilon,E+\epsilon]$ are exponentially suppressed. If we choose $t_{\rm RMS}$ to scale proportionally with $1/\epsilon$, then we achieve the desired energy filtering with energy resolution $\epsilon$. The actual peak width will be a factor of $1/\sqrt{N}$ narrower than $\epsilon$, but that is needed to get exponential suppression as a function of $N$ for all energies $E_{\rm obj}$ further than $\epsilon$ from $E$.

In Fig. 3 we plot $\ln P_{N}$ versus $N$ for Gaussian random values of $t_{n}$ using several values for $\theta_{\rm RMS}\equiv(E_{\rm obj}-E)\tfrac{t_{\rm RMS}}{2}$. We show $\theta_{\rm RMS}=0.5$ (open circles), $\theta_{\rm RMS}=1.0$ (open triangles), $\theta_{\rm RMS}=2.0$ (open diamonds), and $\theta_{\rm RMS}=3.0$ (open squares). We present the predicted asymptotic scaling, $\ln P_{N}=-N\ln 4$, with filled circles. We see that for $\theta_{\rm RMS}$ greater than $1$, the expected asymptotic scaling is achieved. Therefore, if $t_{\rm RMS}$ is larger than twice the inverse spacing between energy levels, then $P_{N}$ scales as $\tfrac{1}{4^{N}}$ for large $N$.

Let us now consider what happens for an arbitrary initial state $\ket{\psi_{I}}$. Let $E_{j}$ be the energy eigenvalue nearest to $E$, and let $\ket{E_{j}}$ be the corresponding eigenvector. In the limit of large $N$, the probability that we measure the $\ket{1}$ state $N$ times in a row is $pP_{N}$, where $p$ is the overlap probability of the initial state with $\ket{E_{j}}$, and $P_{N}$ is the success probability for $E_{\rm obj}=E_{j}$. By tuning $E$ equal to $E_{j}$, this probability becomes $p$. If we require that the spectral weights of all other energy eigenvectors outside the interval $[E-\epsilon,E+\epsilon]$ are suppressed by a factor $\delta$, then the computational effort scaling for the rodeo algorithm is $Nt_{\rm RMS}/p=O[|\log\delta|/(p\epsilon)]$.

In order to successfully determine any given energy eigenvalue $E_{j}$ with error $\epsilon$, the computational cost scales as $O[(\log\epsilon)^{2}/(p\epsilon)]$. The search process involves $O(\log\epsilon)$ sequential scans of the energy, each scan sweeping over an energy range that is some constant factor $K$ smaller than the previous scan. Each scan is performed for several evenly spaced values of $E$, with a fixed number of rodeo cycles, and $t_{\rm RMS}$ a factor of $K$ larger than that used for the previous scan. The total time evolution required will scale as $O(1/\epsilon)$, and the factor of $(\log\epsilon)^{2}/p$ comes from the required statistics needed to perform the energy scans successfully with high probability. The resulting performance as a function of $\epsilon$ is close to the $O(1/\epsilon)$ bound set by the Heisenberg uncertainty principle. For comparison, the computational effort for phase estimation is $O[1/(p\epsilon)]$ plus an additional cost that is $O[(\log\epsilon)^{2}]$ associated with the quantum Fourier transform. Iterative phase estimation eliminates the need for the quantum Fourier transform, but is suitable only for finding the energy of a pure eigenstate. We note that the direct calculation of the expectation value of the Hamiltonian for a pure eigenstate requires $O(1/\epsilon^{2})$ measurements due to statistical errors.

In order to successfully prepare any given eigenstate $\ket{E_{j}}$ with a residual orthogonal component that has magnitude $\Delta$, the computational cost is $O(\log\Delta/p)$. For this case, we keep $t_{\rm RMS}$ fixed but large enough that we are filtering out only the desired eigenstate. We perform $N=O(\log\Delta)$ cycles of the rodeo algorithm, and this must be multiplied by $1/p$ for the number of measurements required. In preparing the eigenstate, it is important to keep $E$ centered on the peak maximum associated with $E_{j}$. Re-centering $E$ with each cycle requires only a constant overall factor in the computational cost that is independent of $p$ and $\Delta$. In contrast, the computational cost for the same task using phase estimation requires an effort that scales as $O[1/(p\Delta)]$. Adiabatic evolution requires an effort that is $O(1/\Delta)$ times a function of $p$ which depends on the adiabatic path connecting the initial and final Hamiltonians Wiebe and Babcock (2012). We see that for eigenstate preparation, the rodeo algorithm is exponentially faster than both phase estimation and adiabatic evolution in the limit $\Delta\rightarrow 0$.

From Eq. (3), we can also derive two useful estimates for $\Delta$, the magnitude of the residual orthogonal component, as a function of the number of rodeo cycles, $N$,

	$\displaystyle F_{A}\equiv\sqrt{2^{-N}(1-p)/[p+2^{-N}(1-p)]},$
	$\displaystyle F_{G}\equiv\sqrt{4^{-N}(1-p)/[p+4^{-N}(1-p)]}.$		(4)

$F_{A}$ is appropriate for the case when $N$ is much smaller than the number of eigenstates with nonzero overlap with our initial state. This corresponds with a spectral suppression factor of $1/2$ at each order for the undesired eigenstates, corresponding with the arithmetic mean of $\cos^{2}\theta$. $F_{G}$ is appropriate for the case when $N$ is much larger than the number of eigenstates with nonzero overlap with our initial state. This corresponds with a spectral suppression factor of $1/4$ at each order for the undesired eigenstates, corresponding with the geometric mean of $\cos^{2}\theta$.

The rodeo algorithm depends heavily on the initial-state overlap probability $p$. It is, therefore, very helpful to improve the quality of the initial state. One possible approach is to make a variational ansatz based on physical intuition about the nature of the eigenvector. Another strategy is to use domain decomposition to define a variational ansatz as a tensor product of wave functions on smaller subsystems. Yet another approach is some combination of variational methods and adiabatic evolution such as the quantum approximate optimization algorithm Farhi et al. (2014).

As a first application of the rodeo algorithm, we consider the spin-$\tfrac{1}{2}$ Heisenberg model in a uniform magnetic field with $10$ sites forming a closed one-dimensional chain Heisenberg (1928). The Hamiltonian has the form

$$H_{\rm obj}=J\sum_{\braket{j,k}}\vec{\sigma}_{j}\cdot\vec{\sigma}_{k}+h\sum_{j}\sigma^{z}_{j},$$

(5)

where $J$ is the exchange coupling, $\vec{\sigma}_{j}$ are the Pauli matrices on site $j$, $\braket{j,k}$ indicates nearest neighbors, and $h$ is the coupling to a uniform magnetic field in the $z$ direction. We consider the antiferromagnetic case with values $J=1$ and $h=3$. For our initial state we use an alternating tensor product state,

$$\ket{\psi_{I}}=\ket{0101010101}.$$

(6)

Since our initial state has a high degree of symmetry, we expect our initial state to have nonzero overlap with a relatively small number of energy eigenstates.

Let us label the energy eigenstates of $H_{\rm obj}$ as $\ket{E_{j}}$. We define the initial-state spectral function as $S(E)=|\braket{E_{j}}{\psi_{I}}|^{2}$ for $E=E_{j}$ and $S(E)=0$ otherwise. For the case of exact degeneracy, we sum the contribution from all degenerate energy states. In Fig. 4, we plot the initial-state spectral function using the rodeo algorithm for the Heisenberg spin chain with $N=3$ (thin blue line), $6$ (thick green line), and $9$ (medium red line) cycles. We have averaged over $20$ sets of Gaussian random values for $t_{n}$ with $t_{\rm RMS}=5$. This averaging over sets of random values for $t_{n}$ decreases the stochastic noise and results in a roughly constant background that can be distinguished from the spectral signal. For comparison, we show the exact initial-state spectral function with black open circles. We see that the agreement obtained using the rodeo algorithm is excellent.

The real challenge will be to perform these calculations on quantum computing devices with gate errors, measurement errors, and short decoherence times. But it is promising that we can obtain good results even though neither $t_{\rm RMS}$ nor $N$ are very large. The resulting short gate depth is crucial for implementation on noisy quantum devices. In the Supplemental Materials, we show the corresponding results with $t_{\rm RMS}=1$. Even in that case, we can clearly identify the spectrum of energy states with strong overlap with the initial state.

In addition to computing the initial-state spectral function, we can also prepare any energy eigenstate that has nonzero overlap with our initial state. In the Supplemental Materials, we show the overlap probability with energy eigenvector $\ket{E_{j}}$ after $N$ cycles of the rodeo algorithm. All of the energy eigenvectors with nonzero overlap with our initial state can be prepared with a relatively small number of rodeo cycles. After an energy eigenstate has been prepared using the rodeo algorithm, we can measure any properties of that eigenstate that we choose, such as expectation values of observables, transition matrix elements to other energy eigenstates, and linear response functions. Several examples of such calculations using prepared energy eigenstates will be discussed in a future publication. See also Ref. Roggero and Carlson (2019); Roggero et al. (2020); Roggero (2020) for some approaches to computing linear response functions and spectral densities.

In Fig. 5, we show results for the logarithm of the wave function error for the rodeo algorithm when preparing the energy eigenstate $\ket{E_{j}}$ corresponding to $E_{j}=-18.1$. We plot the logarithm of the magnitude of the orthogonal complement, $\log\Delta$, versus the total time evolution, $T$, for $t_{\rm RMS}=1$. For comparsion, we show the estimates $\log F_{A}$ and $\log F_{G}$ defined in Eq. (4). As expected, for small $T$, $\log F_{A}$ provides a good estimate, while for very large $T$, $\log F_{G}$ provides a better estimate. For comparison, we also show the analogous results obtained with the same initial state but instead using phase estimation and adiabatic evolution. For the adiabatic evolution calculation, we use the initial Hamiltonian $H_{I}=\sum_{j=1}^{10}(-1)^{j}\sigma^{z}_{j}$, with an interpolating function $H(t)=\cos^{2}[\pi t/(2T)]H_{I}+\sin^{2}[\pi t/(2T)]H_{\rm obj}$. We note that phase estimation and adiabatic evolution are similar in performance, while the rodeo algorithm is exponentially faster than both.

Even though the rodeo algorithm is more efficient than adiabatic evolution for eigenstate preparation, we can use adiabatic evolution as a preconditioner for the rodeo algorithm, in order to amplify the overlap of the initial state with the desired eigenvector. In Table S1, we show the overlap probability for energy eigenvector $\ket{E_{j}}$ with $E_{j}=-18.1$ after preconditioning with adiabatic evolution for time $t_{\rm AE}$ and applying $N$ cycles of the rodeo algorithm with $t_{\rm RMS}=5$. We see that by preconditioning with $t_{\rm AE}=5$, we achieve a more than sevenfold increase in the initial state overlap probability. We gain a significant computational advantage when using adiabiatic evolution as a preconditioner for the rodeo algorithm and expect this strategy to be very useful for larger system calculations.

As a second application of the rodeo algorithm, we consider the Anderson localization model in one dimension, which describes the transition between extended and localized electronic states in the presence of spatially-varying disorder Anderson (1958); Abrahams et al. (1979). Our object Hamiltonian $H_{\rm obj}$ describes a single particle on a periodic, one-dimensional lattice with $100$ sites. Let us denote the position basis states as $\ket{k}$ with $k=0,\cdots,99$. The matrix elements of the object Hamiltonian are

$$[H_{\rm obj}]_{k^{\prime},k}=-\delta_{k^{\prime},k+1}-\delta_{k^{\prime},k-1}+c_{k}\delta_{k^{\prime},k},$$

(7)

where the coefficients $c_{k}$ provide diagonal disorder that control the amount of localization of the electronic wave functions. In this example, we take the diagonal terms $c_{k}$ to be Gaussian random numbers with zero mean and root-mean-square value equal to $\tfrac{1}{2}$. With this amount of diagonal disorder, the resulting orbitals are rather localized in space. In the Supplemental Materials, we also present results for the case where diagonal disorder is much weaker and the energy eigenstates are delocalized in space.

As a variational ansatz for the ground state, we take our initial state to be the basis state $\ket{k_{\rm min}}$ such that $c_{k_{\rm min}}$ is the lowest diagonal element of the Hamiltonian. This choice reflects our physical intuition that the ground state is spatially localized. In Fig. 6, we plot the initial-state spectral function using the rodeo algorithm for the Anderson localization model with the root-mean-square diagonal disorder equal to $\tfrac{1}{2}$. We consider $N=3$ (thin blue line), $6$ (thick green line), and $9$ (medium red line) cycles. We have also averaged over $20$ sets of Gaussian random values for $t_{n}$ with $t_{\rm RMS}=10$. For comparison, we show the exact initial-state spectral function with black open circles. We see that the exact spectral function is well reproduced. The fact that our point-like initial state has significant overlap with only a small fraction of energy eigenstates is an indication of the localized character of the orbitals.

In this letter, we have presented a new method called the rodeo algorithm for preparing quantum eigenstates and determining spectral properties. It uses a tunable energy filter and stochastic methods to prepare any eigenstate of a given quantum Hamiltonian and can in fact be used to prepare the eigenstates of any quantum observable. It is an efficient algorithm with exponential convergence that can be performed using circuits with relatively short gate depth. In particular, the speed for eigenstate preparation is exponentially faster than that for phase estimation or adiabatic evolution. The rodeo algorithm can be combined with variational methods and/or quantum adiabatic evolution to provide a tool for solving the quantum many-body problem even when there is no a priori information about the target eigenvector. It has the potential for wide applicability across many different fields, including combinatorial optimization problems, strongly-correlated electrons, nuclear structure and dynamics, and lattice quantum chromodynamics. The results presented here were obtained from classical computation, but we are now working to implement the rodeo algorithm on quantum devices using several quantum Hamiltonians. The results will be presented in future publications.