Quantum random power method for ground state computation

A main component in our algorithm is the estimation of the matrix element of the Hamiltonian matrix polynomial $p(H)$, for example, for any indices $i,j\in[N]$,

where $\ket{i}$ denotes the computational basis element corresponding to the integer $i$, and the same for $\ket{j}$. To estimate this quantity, one can consider two approaches: one based on the Hamiltonian evolution [23] and the other based on the quantum singular value transformation (QSVT) [12].

Hamiltonian evolution is built on the ability to implement $e^{iHt}$. We estimate a matrix polynomial of $H$ using a Fourier-series approximation. By the property of the Fourier expansion, we obtain the following approximation for any smooth function $p(x)$ on $[-\pi,\pi]$ and some sufficiently large $M$,

The coefficient $c_{m}$ depends on the choice of a kernel [43]. In our applciations, we used the Fej́er kernel as discussed in Section VIII. Using this approximation and Hamiltonian simulation, we obtain access to a matrix polynomial of $H$ approximately,

Especially, estimating matrix elements of $p(H)$ and the corresponding elements of the exponentials of $H$ are equivalent as follows,

To compute this matrix element, we use the Hadamard test in Fig. 2.

Alternatively, if we are given a block encoding of $H$, we construct a block encoding of $p(H)$ by combining quantum signal processing techniques [25, 12, 5]. Then, we are able to do the same task by replacing $e^{iHt}$ in Fig. 2 with $U_{p(H)}$, the block encoding of $p(H)$, and additional ancilla qubits,

For our algorithm, we let $p(x)$ be a Chebyshev polynomial, which was motivated by the Chebyshev filtering method [45].

This method requires a priori bounds to be properly implemented. For example, we need $\lambda_{lb}$ and $\lambda_{ub}$ such that

with which the given Hamiltonian is linearly transformed as

Applying this transformation and a Chebyshev polynomial, $p(x)$, to $H$ amplifies the contributions of small eigenvalues but dampens those of large eigenvalues of $H$ as in Fig. 3, say,

Notice that the high-lying eigenvalues in $[\lambda_{lb},\lambda_{ub}]$ are mapped into $[-1,1]$ where the Chebyshev polynomials value in $[-1,1]$, while the low-lying eigenvalues lying in $[\lambda_{1},\lambda_{lb})$ are located outside $[-1,1]$. Due to the fact that the Chebyshev polynomial is amplified outside $[-1,1]$, the contributions of low-lying eigenvalues of $p(H)$ increases rapidly.

A search of a priori bounds is an important step in the context of filtering algorithms [45, 22, 23]. In classical algorithms, the task can be completed with a combination of Rayleigh quotient computations and Lanczos algorithm as demonstrated in [45]. But this may yield an exponential arithmetic operation cost with system size. Alternatively, one may use the binary amplitude estimation with a block encoding of $H$ provided [22] or with an access to the Hamiltonian evolution $U=e^{iHt}$ [23].

Once we have access to the matrix polynomial $p(H)$ as discussed above, we build an iterative algorithm in which a gradient estimator is constructed from the elements of the matrix polynomial. To do so, we propose an estimator for the matrix-vector product $H\bm{x}$ by borrowing the ideas of randomized coordinate approach [30, 18] and the stochastic gradient descent [3]. A more precise definition is found in (16). For simple discussion, we look at the property that the estimator in (16) is proportional to $p(H)\bm{x}$ in the expectation. Essentially, we define an iterative algorithm as follows,

Here $\bm{g}_{k}$ satisfies that

For example, we sample two indices $i_{k},j_{k}\in[N]$ first, and then estimate the corresponding matrix element to define the random vector as

where $(\bm{x}_{k})_{i_{k}}$ is the $i_{k}$-th component of $x_{k}$, $\bm{e}_{j_{k}}$ is the classical computational basis vector corresponding to index $j_{k}$. For simplicity, we assume that there is no systematic error and let $\epsilon$ denote the sampling noise. Noting that the two independent random samplings occur for indices $i_{k},j_{k}$, we then observe that averaging $\bm{g}_{k}$ over the indices yields the desired matrix-vector product up to a scaling factor,

We remark that the iterative scheme (11) is typically used for uncontrained optimization problems in machine learning, but still applicable for the problem (2) which is a constrained optimization problem with normalization as demonstrated in [26, 38]. To understand why it converges, we naively consider the dynamics of the averaged iteration. Let us denote $\bm{y}_{k}:=\mathbb{E}[\bm{x}_{k}]$. By taking expectation on the equation, we can obtain

which resembles the power method without normalization. With sufficiently small $a>0$, the iteration will converge to the eigenvector corresponding to the maximum eigenvalue of the matrix $I-ap(H)$ or the ground state of $H$. From this viewpoint, we can expect that the iteration (11) converges to the ground state of $H$ on average. In fact, a stronger convergence property holds as will be shown in Theorem 4 that contrasts to the previous works [26, 38].

We extend the estimator in (13) to the one defind with more than one matrix element sampling as follows,

We note that the quantity, $\xi(H,i_{r},i_{c})$, is an estimate of the element of a desired matrix. For our applications, the estimate $\xi$ will approximate the element of $p(H)$, such as (7). For the sake of simplicity, we first consider the case where $\xi(H,i_{r},i_{c})=\bra{\bm{i}_{r}}H\ket{\bm{i}_{c}}$. By the definition (16) and the definition of matrix 2-norm, we have

Note that the example (13) is a special case of this estimator when $m_{r}=m_{c}=1$ without any quantum noise. The following proposition shows that the estimator (16) is proportional to $H\bm{x}$ in the expectation, and has the variance scaling with $\norm{\bm{x}}^{2}$.

The proof of Proposition 2 is shown in Section VIII.1.

We remark that when $\xi$ in (16) estimates the element of $p(H)$, Proposition 2 can be restated with $p(H)$ instead of $H$.

As discussed in Section II, we incorporate the estimator (16) to the randomized power method (11). This outlines the algorithm as shown in Algorithm 1 that produces an approximation of a ground state and the corresponding energy estimate. As a stopping criterion, the Rayleigh quotient,

, is computed every iteration. To do so efficiently, we exploit the sparsity of $H$ and $\bm{g}$ in (16), otherwise the cost of a brute-force calculation of (21) scales with the dimension $N$. The sparisty of $H$ may inherit from that $H$ is a linear combination of $\mathcal{O}(poly(n))$ many Pauli strings. As a result, the estimation of (21) relies on the following recursive relations

Notice that the first relation requires $\mathcal{O}(m_{r}m_{c})$ cost, since $\bm{g}_{t-1}$ is $m_{c}$-sparse and the $m_{r}$-sparse columns of $H$ corresponding to $\bm{g}_{t-1}$ are added to the vector $H\bm{x}_{t-1}$. Similarly, the second relation needs $\mathcal{O}(m_{c})$ cost. The last relation takes $\mathcal{O}(m_{c}\min\{m_{r},m_{c}\})$ cost as the last term on the right side dominates overall cost. Therefore, the per-iteration cost for the estimation of (21) scales only as $\mathcal{O}(m_{r}m_{c})$.

We note that in Algorithm 1, the gradient estimate (16) may involve bias due to the use of the quantum filtering technique (10) and the sampling without repetition. This means that the matrix-element estimation is performed not for $p(H)$, but for a perturbed one $\widetilde{p(H)}$. As a consequence, Algorithm 1 converges to an approximation of ground state of $H$ as in the following theorem. In the next section, we provide a lower bound of the fidelity between perturbed and true ground states to ensure that the more accurate the estimator (16) is computed, the higher the achievable fidelity of Algorithm 1 is, as being arbitrarily close to one.

The following theorem guarantees convergence with probability one to a given precision if a sufficiently small step size is chosen. As usually observed in variants of the power method [26, 38], the iteration complexity has a logarithmic dependence on the overlap of the ground state. By combining Theorem 6, we present the overall query complexity for Algorithm 1.

The proof of Theorem 4 is shown in Section VIII.3.

In Algorithm 1, the elements of a matrix polynomial are estimated without repetition. The resulting matrix, $\widetilde{p(H)}$, will be one perturbed from the desired matrix $p(H)$ by at least the sampling noise, denoted as $E$, namely,

Here, we denote an error-free matrix and noisy one by $A$ and $\widetilde{A}$. Suppose that these matrices are Hermitian as in the situation (25). It may be challenging to derive a perturbative relation between the ground states of $A$ and $\widetilde{A}$ when the error $E$ is considered to be arbitrary. But when $E$ is smaller than some value, we can show that a lower bound for the fidelity between a perturbed ground state and the true ones depends only on the magnitude of noise, $\norm{E}$.

Let $\lambda_{1}(A)$ and $\lambda_{2}(A)$ be the first and second smallest eigenvalues of $A$, and similar for $\widetilde{A}$. Let $E=\widetilde{A}-A$ denote the error. The following lemma shows that when $\norm{E}$ is smaller than a certain value, the overlap of any ground state of $\widetilde{A}$ with the subspace of ground states of $A$ increases to $1$ linearly with $\norm{E}$.

The proof of this lemma is shown in Section VIII.

For the Hamiltonian simulation approach (i.e. (27) and (28)), we simulate Algorithm 1 using Qiskit [36]. To estimate the gradient, we incorporate the Chebyshev filtering method using the first-order Trotter formula, the Fourier-series approximation with the Fej́er kernel, and Hamiltonian simulation. Thus three types of errors are involved: numerical errors due to the Fourier expansion and the Trotter formula as well as the sampling noise from the Hadamard test. For the block-encoding approach (i.e. (27) and (30)), we simulate the quantum singular value transformation (QSVT) on MATLAB by computing the eigen-decomposition of a given Hamiltonian and constructing a corresponding matrix polynomial to mimic a block encoding of a matrix polynomial of $H$. In this approach, only sampling noise is considered. In all numerical experiments, we save the estimated matrix elements once and use the estimated matrix to carry out 100 independent simulations for Algorithm 1 with the same initial guess.

We consider the models (27) and (28) for different $D/J$ with $J=1$. For every $D/J$, we set the same parameter values, $n=2$, and the shot number to be $4\times 10^{5}$, the degree of the Fej́er kernel to be $30$, the degree of the Chebyshev filtering to be $3$, and $m_{r}=m_{c}=1$ in Algorithm 1. We use the first-order Trotter formula and set $\max\{k^{2}+20,50\}$ to be the number of Trotter steps when estimating matrix elements of $e^{iHk}$ where $k$ is a positive integer. In Fig. 4 that compares the performance of Algorithm 1 with and without polynomial filters, we measure the growing rate of the average fidelity for each $D/J$, namely,

where the constants $A$ and $B$ are

Note that this quantity measures how fast the average fidelity increases, according to the technical detail (101). In the following numerical results, $m_{r}=m_{c}=1$. Here we set $\epsilon=0.9$. Under this setting, the quantity (31) tells us the worst-case growing rate of the average fidelity unitl it reaches the precision $\epsilon$. For fair comparison, we compare the cases without scheduling decreasing step sizes in the bottom panels of Fig. 4, as the step size $a$ is a constant. For the polynomial filtering, we set the upper and lower bounds in (8) as $\lambda_{ub}=\lambda_{\max}+0.2(\lambda_{\max}-\lambda_{\min})$ and $\lambda_{lb}=\lambda_{\min}+0.2(\lambda_{\max}-\lambda_{\min})$ for the TFIM and the XXZ model.

The top panels in Fig. 4 show the convergence of our algorithm with and without polynomial filters, and with and without scheduling decreasing step sizes. The graphs show the average of fidelities obtained from the 100th iterates of 100 independent simulations with a fixed initial guess for each $D/J$. For the two models, we can see that in all cases, the optimization with the filtering method converges much faster than the ones without it. One reason may be that the convergence rate of the fidelity is related to the growing rate of the average fidelity (31). We observe that the overall growing rate is relatively higher when the filtering is applied as shown in the bottom panels of Fig. 4. In addition, the convergence result in Fig. 4 shows that Algorithm 1 with the filtering method approximates the true ground state well, despite the additional errors involved (e.g. the Trotter error, the Fourier approximation error, and the sampling error). Together with the convergence guarantee in Theorem 4, this can be verified by checking the high fidelities in the last column of Table 1, which are the fidelity of the ground state from the errorneous matrix contructed from Algorithm 1 with the true ground state.

Regarding Lemma 5, we observe that the perturbation error (LHS), $\norm{E}$, is smaller than half of the enlarged spectral gap (RHS), $\frac{1}{2}\absolutevalue{\lambda_{1}(p(H))-\lambda_{2}(p(H))}$, in all cases. This means that the condition in Lemma 5 is satisfied in the cases, and so is the lower bound in the lemma. Accordingly, the ordering of the first two smallest eigenvalues is preserved even in the presence of sampling noise, and moreover the interval $[\lambda_{1}-e,\lambda_{1}+e]$ filters out the smallest eigenvalue of the estimated matrix from the second smallest one.

In this section, we consider the TFIM (27) with $J=1,D=1.5$ and the Hubbard model (30) in weak ($t=1,U=1$) and strong coupling regimes ($t=1,U=6$). we set $n=10$, the parameters $m_{r}=20,m_{c}=20$ in Algorithm 1, and the degree of the Chebyshev polynomial to be $7$ for the TFIM and $8$ for the Hubbard model. We note that there is a non-degenerate ground state for the cases of (27) and (30) in the weak-coupling regime, while there are two degenerate ground states in the strong-coupling regime for (30). For this case, we measure the fidelity as follows,

where the $\psi_{GS,i}$’s for $i=1,2$ denote two orthogonal ground states. This is a special case of more general result in A lower bound of the fidelity between a perturbed ground state and the subspace of true ones, which considers the number of ground states to be $2$. For the polynomial filtering, we set the upper and lower bounds in (8) as $\lambda_{ub}=\lambda_{\max}+0.1(\lambda_{\max}-\lambda_{\min})$ and $\lambda_{lb}=\lambda_{\min}+0.03(\lambda_{\max}-\lambda_{\min})$ for the TFIM and the Hubbard model in the weak-coupling regime, and the same upperbound and $\lambda_{lb}=\lambda_{\min}+0.02(\lambda_{\max}-\lambda_{\min})$ for the HM in the strong-coupling regime. These lower bounds, $\lambda_{lb}$’s, seperate $\lambda_{1}$ and $\lambda_{2}$ for the TFIM model, $\lambda_{3}$ and $\lambda_{4}$ for the HM in the weak-coupling regime, and $\lambda_{15}$ and $\lambda_{16}$ for the HM in the strong-coupling regime, respectively.

Fig. 5 shows the effect of the polynomial filter and the shot number on the convergence of Algorithm 1. One can observe that the number of shots influences the performance of Algorithm 1. For example, the result with the largest shot number is much better than the others. Additionally, from the bottom panels of Fig. 5, it is evident that the use of the polynomial filter improves overall performance of Algorithm 1, exhibiting better convergence than the result without it. Importantly, the result in Fig. 5 shows high fidelities, which means that Algorithm 1 with the filtering method provides a fairly good approximation even with the sampling noise. By the convergence theorem Theorem 4, this can be understood by looking at the fidelity of approximate ground states with the true ones as in the last columns in Table 2.

As discussed in the previous section, we compare the perturbation error (LHS) to the enlarged spectral gap (RHS). For each of the three models, the condition in Lemma 5 is not satisfied in some cases where the LHS is larger than the RHS, but the interval $[\lambda_{1}-e,\lambda_{1}+e]$ filters the smallest perturbed eigenvalue $\widetilde{\lambda_{1}}$ from the second one, which means that the ordering of the first two smallest eigenvalues is preserved even in the presence of sampling noise. However, this observation implies that the condition in Lemma 5 can be restrictive in practice. It would be interesting to find a more relaxed condition, and we leave this to future work.

For the Hubbard model in a strong coupling regime, among 100 simulation results, Fig. 6 in Section VIII.4 plots 12 that display the respective fidelities with the two degenerate ground states and the mean-square sum with normalization as defined in (1). The graphs mean that our method converges to a ground state as shown in Theorem 4.

We use an outcome of the Weyl’s inequality, which states that for Hermitians $A,B\in\mathbb{C}^{N\times N}$ and any $k\in\{1,2,..,N\}$,

where the $\lambda_{k}(A)$ denotes the $k$-th smallest eigenvalue of $A$ and similarly for $A+B$ and $B$.

To prove Lemma 5, we claim that

where $\mathbb{I}$ is the indicator function that values one on the given interval and zero otherwise. To show this, we check that $\lambda_{1}(\widetilde{A})\in[\lambda_{1}(A)-e,\lambda_{1}(A)+e]$, and $\lambda_{2}(\widetilde{A})\not\in[\lambda_{1}(A)-e,\lambda_{1}(A)+e]$. By (34), $\lambda_{1}(\widetilde{A})\in[\lambda_{1}(A)-e,\lambda_{1}(A)+e]$. Again, by (34), $\lambda_{2}(\widetilde{A})\geq\lambda_{2}(A)-e$. By the assumption, $\lambda_{2}(A)-e>\lambda_{1}(A)+e$. Thus, the claim is proven.

Notice that the interval $[\lambda_{1}(A)-e,\lambda_{N}(A)+e]$ includes the spectra of $A$ and $\widetilde{A}$ by the Weyl’s inequality, where $\lambda_{N}(A)$ is the largest eigenvalue of $A$. Restricting to this interval, by the Stone–Weierstrass theorem, we can take a sequence of polynomials $\{p_{n}(x)\}_{n=1}^{\infty}$ that uniformly converges to the indicator function in (35). Then for any $\epsilon$, there exists a $n>0$ such that $p_{n}$ is $\frac{\epsilon}{2}$-close to the indicator function on the interval in the uniform norm. Thus, it immediately follows that

where $\mathbb{I}:=\mathbb{I}_{[\lambda_{1}(A)-e,\lambda_{1}(A)+e]}$. Together with the above claim, we observe

where $R_{n}(A,E):=p_{n}(\widetilde{A})-p_{n}(A)$. We note that $\mathbb{I}(A)=P_{\lambda_{1}(A)}$ by the selection of $e$. In addition, $\norm{R_{n}(A,E)}\leq C\norm{E}$, since all terms in $R_{n}(A,E)$ have at least one $E$ as a factor, and $C$ is determined by $\norm{A}$, $e$, and $\epsilon$, based on the triangle inequality argument. Let $\ket{\widetilde{\bm{q}}}=\sum_{j=1}^{g}\beta_{j}\ket{\widetilde{\bm{q}}_{j}}$ be a ground state of $\widetilde{A}$. By applying the braket of the state $\ket{\widetilde{\bm{q}}}$ to (37), we therefore obtain

which completes the proof.

To implement the Chebyshev filtering [45] with the Fourier expansion and Hamiltonian simulation, we first transform the Hamiltonian as (9), and normalize the transformed matrix such that its spectrum belongs to $[-\pi,\pi]$, which is required to simulate Hamiltonian. Denote the transformed matrix by $\overline{H}$. Let $\alpha$ be set as $\norm{\overline{H}/\alpha}\leq\pi$.

Our goal is to estimate the matrix elements of $p_{\ell}(H)$ for the Chebyshev polynomial of degree $\ell$, $p_{\ell}(x)$, using the Fourier expansion with the Fej́er kernel. The Fourier expansion with the Fej́er kernel of order $K$ for $p_{\ell}(\alpha x)$ is represented as a convolution form

where $c_{k}=\frac{K-k+1}{2\pi K}\int_{-\pi}^{\pi}p_{\ell}(\alpha x)e^{-ikx}dx$, and $S_{K}(x)=\frac{\sin^{2}\left(Kx/2\right)}{K\sin^{2}(x/2)}$. By the Fej́er’s theorem, the RHS converges uniformly to the LHS as $K\rightarrow\infty$.

The $c_{k}$ can be computed via a recursive relation for the following auxiliary variable,

By the integral of parts, we have