Measurement-efficient quantum Krylov subspace diagonalisation

Finding the ground-state energy of a quantum system is of vital importance in many fields of physics [1, 2, 3]. The Lanczos algorithm [4, 5] is one of the most widely used algorithms to solve this problem. It belongs to the Krylov subspace methods [6], in which the solution usually converges to the true answer with an increasing subspace dimension. However, such methods are unscalable for many-body systems because of the exponentially-growing Hilbert space dimension [7]. In quantum computing, there is a family of hybrid quantum-classical algorithms that can be regarded as a quantum generalisation of the classical Lanczos algorithm [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. Following Refs. [11, 15, 16], we call them quantum Krylov subspace diagonalisation (KSD). These algorithms are scalable with the system size by carrying out the classically-intractable vector and matrix arithmetic on the quantum computer. They possess potential advantages as the conventional quantum algorithm to solve the ground-state problem, quantum phase estimation, requires considerable quantum resources [21, 22], while the variational quantum eigensolver is limited by the ansatz and classical optimisation bottlenecks [23, 24].

However, Krylov subspace methods are often confronted with the obstacle that small errors can cause large deviations in the ground-state energy. This issue is rooted in the Krylov subspace that is spanned by an almost linearly dependent basis [25]. Contrary to classical computing, in which one can exponentially suppress rounding errors by increasing the number of bits, quantum computing is inherently subject to statistical error. Since statistical error decreases slowly with the measurement number $M$ as $\propto 1/\sqrt{M}$, an extremely large $M$ can be required to reach an acceptable error. In this case, although quantum KSD algorithms perform well in principle, the measurement cost has to be assessed and optimised for realistic implementations [16, 20].

In this work, we present a general and rigorous analysis of the measurement cost in quantum KSD algorithms. Specifically, we obtain an upper bound formula of the measurement number that is applicable to all quantum KSD algorithms. Then, we propose an algorithm to construct the Hamiltonian-power Krylov subspace [6]. In our algorithm, we express the product of Hamiltonian power and a Gaussian function of Hamiltonian as an integral of real-time evolution. In this way, the statistical error decreases exponentially with the power, which makes our algorithm measurement-efficient. We benchmark quantum KSD algorithms by estimating their measurement costs in solving the anti-ferromagnetic Heisenberg model and the Hubbard model. Various lattices of each model are taken in the benchmarking. It is shown that the measurement cost in our algorithm is typically orders of magnitude fewer than other algorithms to reach an error comparable to the classical Lanczos algorithm. We also demonstrate the measurement efficiency of our algorithm by composing a ground-state projector.

First, we introduce the KSD algorithm and some relevant notations. The algorithm starts with a reference state $|{\varphi}\rangle$ (or a set of reference states [10, 11]; we focus on the single-reference case in this work). Then we generate a set of basis states

where $H$ is the Hamiltonian, and $f_{1},f_{2},\ldots,f_{d}$ are linearly-independent $(d-1)$-degree polynomials (to generate a $d$-dimensional Krylov subspace). For example, it is conventional to take the power (P) function $f_{k}(H)=H^{k-1}$ in the Lanczos algorithm. These states span a subspace called the Krylov subspace. We compute the ground-state energy by solving the generalised eigenvalue problem

where $\mathbf{H}_{k,q}=\langle{\phi_{k}}|H|{\phi_{q}}\rangle$ and $\mathbf{S}_{k,q}=\langle\phi_{k}|\phi_{q}\rangle$. Let $E_{min}$ be the minimum eigenvalue of the generalised eigenvalue problem. The error in the ground-state energy is

where $E_{g}$ is the true ground-state energy. We call $\epsilon_{K}$ the subspace error. One can also construct generalised Krylov subspaces in which $f_{k}$ are functions of $H$ other than polynomials; see Table 1. Appendix A provides a more detailed introduction to the KSD algorithm.

In the literature, subspaces generated by variational quantum circuits [26, 27, 28] and stochastic time evolution [29] are also proposed. Quantum KSD techniques can also be used to mitigate errors caused by imperfect gates [30, 31, 32]. In this work, we focus on Hamiltonian functions because of their similarities to conventional Krylov subspace methods.

In addition to $\epsilon_{K}$, the other error source is the statistical error. In quantum KSD algorithms, matrices $\mathbf{H}$ and $\mathbf{S}$ are obtained by measuring qubits at the end of certain quantum circuits. Measurements yield estimators $\hat{\mathbf{H}}$ and $\hat{\mathbf{S}}$ of the two matrices, respectively. Errors in $\hat{\mathbf{H}}$ and $\hat{\mathbf{S}}$ depend on the measurement number. Suppose the budget for each complex matrix entry is $2M$ measurements (the budget for each part is $M$). Variances of matrix entries have upper bounds in the form

and

where $C_{\mathbf{H}}$ and $C_{\mathbf{S}}$ are some factors depending on the measurement protocol. For example, suppose each entry of $\mathbf{H}$ and $\mathbf{S}$ can be expressed in the form $\sum_{s}q_{s}\langle{\varphi}|U_{s}|{\varphi}\rangle$, where $U_{s}$ are unitary operators, and we measure each term using the Hadamard test [33], then the variance of the entry has the upper bound $2C^{2}/M$, where $C=\sum_{s}|q_{s}|$.

We quantify the error in a matrix with the spectral norm. The distributions of $\|\hat{\mathbf{H}}-\mathbf{H}\|_{2}$ and $\|\hat{\mathbf{S}}-\mathbf{S}\|_{2}$ depend on not only the measurement number but also the correlations between matrix entries. We can measure matrix entries independently, then their distributions are uncorrelated. However, certain matrix entries take the same value, and we can measure them collectively. For example, for the P basis, all entries $\mathbf{H}_{k-q,q}$ take the same value as $\mathbf{H}_{k-1,1}$. We only need to measure one of them on the quantum computer, and then these entries become correlated. In Appendix B, we analyse the distributions of spectral norms for each basis in Table 1.

Errors in matrices cause an error in the ground-state energy in addition to $\epsilon_{K}$. It is common that the subspace basis is nearly linearly dependent, which makes the overlap matrix $\mathbf{S}$ ill-conditioned. Then, even small errors in matrices may cause a serious error in the ground-state energy. Thresholding is a method to overcome this problem [8, 16]. In this work, we consider the regularisation method; see Algorithm 1. An advantage of the regularisation method is that by taking a proper regularisation parameter $\eta$, the resulting energy $\hat{E}_{min}$ is variational, i.e. $\hat{E}_{min}\geqslant E_{g}$ up to a controllable failure probability.

We propose to choose the regularisation parameter $\eta$ according to the distributions of $\|\hat{\mathbf{H}}-\mathbf{H}\|_{2}$ and $\|\hat{\mathbf{S}}-\mathbf{S}\|_{2}$. Let $\kappa$ be the permissible failure probability. We take $\eta$ such that

In Table 2, we list $\eta$ satisfying the above condition for each basis. The value of $\eta$ depends on the failure probability $\kappa$ and the measurement number $M$. Taking $\eta$ in this way, an upper bound of the energy error is given in the following lemma.

See Appendix C for the proof. Notice that we can always subtract a sufficiently large positive constant from the Hamiltonian to satisfy the conditions. The upper bound of the energy error in the regularisation method is $\min_{\mathbf{a}\neq\bf{0}}E^{\prime}(\eta,\mathbf{a})-E_{g}$, up to the failure probability.

In this work, we rigorously analyse the number of measurements. Let $\epsilon$ be the target error in the ground-state energy. We will give a measurement number sufficient (i.e. an upper bound of the measurement number necessary) for achieving the target. For the regularisation method, we already have an upper bound of the error $\min_{\mathbf{a}\neq\bf{0}}E^{\prime}(\eta,\mathbf{a})-E_{g}$. Therefore, we can derive such a measurement number by solving the equation

in which $\eta$ is the unknown. With the solution, we can work out the corresponding measurement number $M$ according to Table 2. If matrix entries are measured independently, the total measurement number is $M_{tot}=4d^{2}M$; we have the following theorem.

See Appendix C for the proof. Similar to Lemma 1, we can always satisfy the condition $E_{g}+\epsilon<0$ by subtracting a sufficiently large positive constant from the Hamiltonian. The above theorem also holds (up to some approximations) for collective measurements if we replace functions $\alpha(\kappa)$ and $\beta(d)$ according to Table 2.

As we have shown, the measurement overhead $\gamma$ depends on the overlap matrix $\mathbf{S}$, i.e. how we choose the basis of Krylov subspace. The main advantage of our algorithm is that we utilise a basis resulting in a small $\gamma$.

To generate a standard Krylov subspace, we choose operators

where $E_{0}$ is a constant up to choice. We call it the GP basis. The corresponding subspace is a conventional Hamiltonian-power Krylov subspace, but the reference state $|{\varphi}\rangle$ has been effectively replaced by $e^{-\frac{1}{2}(H-E_{0})^{2}\tau^{2}}|{\varphi}\rangle$.

The spectral norm of $f_{k}$ operators for the GP basis [defined in Eq. (13)] has the upper bound

where $e$ is Euler’s constant. This upper bound is proven in Lemma 9 in Appendix F. Under the condition $e\tau^{2}>d-1$, the norm upper bound $(\frac{k-1}{e\tau^{2}})^{\frac{k-1}{2}}$ decreasing exponentially with $k$. This property is essential for the measurement efficiency of our algorithm.

We propose to realise the Gaussian-power basis by expressing the operator of interest as a linear combination of unitaries (LCU). The same approach has been taken to realise other bases like power [19, 12], inverse power [18], imaginary-time evolution [34], filter [15] and Gaussian function of the Hamiltonian [35].

Suppose we want to measure the quantity $\langle{\varphi}|A|{\varphi}\rangle$. Here $A=f_{k}^{\dagger}Hf_{q}$ and $A=f_{k}^{\dagger}f_{q}$ for $\mathbf{H}_{k,q}$ and $\mathbf{S}_{k,q}$, respectively. First, we work out an expression in the from $A=\sum_{s}q_{s}U_{s}$, where $U_{s}$ are unitary operators, and $q_{s}$ are complex coefficients. Then, there are two ways to measure $\langle{\varphi}|A|{\varphi}\rangle$. When the expression has finite terms, we can apply $A$ on the state $|{\varphi}\rangle$ by using a set of ancilla qubits; In the literature, LCU usually refers to this approach [36]. Alternatively, we can evaluate the summation $\langle{\varphi}|A|{\varphi}\rangle$ using the Monte Carlo method [37] by averaging Hadamard-test [33] estimates of randomly selected terms $\langle{\varphi}|U_{s}|{\varphi}\rangle$. The second way works for an expression with even infinite terms and requires only one or zero ancilla qubits [38, 39]. We focus on the Monte Carlo method in this work.

Now we give our LCU expression. Suppose we already express $H$ as a linear combination of Pauli operators, it is straightforward to work out the LCU expression of $A$ given the expression of $f_{k}$. Therefore, we only give the expression of $f_{k}$. Utilising the Fourier transformation of a modified Hermite function (which is proved in Appendix F), we can express $f_{k}$ in Eq. (13) as

where $H_{n}(u)$ denotes Hermite polynomials [40], $g_{\tau}(t)=\frac{1}{\tau\sqrt{2\pi}}e^{-\frac{t^{2}}{2\tau^{2}}}$ is the normalised Gaussian function, $x=H-E_{0}$, and $e^{-ixt}$ is the real-time evolution (RTE) operator. Notice that in Eq. (15) $f_{k}$ is expressed as a linear combination of RTE operators $e^{-iHt}$, and the summation in conventional LCU is replaced with an integral over $t$. See Appendix G for how to evaluate this integral using the Monte Carlo method.

With Theorem 1, we can benchmark the measurement cost in quantum KSD algorithms listed in Table 1. We focus on the overhead factor $\gamma$, while ignoring the slight differences in the other two factors in different algorithms (See Table 2). Given a target error $\epsilon$, first we work out the solution $\eta$ to Eq. (9), then $\gamma$ is calculated according to Eq. (12).

Two models of strongly correlated systems, namely, the anti-ferromagnetic Heisenberg model and Hubbard model [7, 48, 49, 50], are used in benchmarking. For each model, the lattice is taken from two categories: regular lattices, including chain and ladder, and randomly generated graphs. We fix the lattice size such that the Hamiltonian can be encoded into ten qubits. For each Hamiltonian specified by the model and lattice, we evaluate all quantum KSD algorithms with the same subspace dimension $d$ taken from $2,3,\ldots,30$. In total, $233$ instances of $(model,lattice,d)$ are generated to benchmark quantum KSD algorithms. The details of the numerical calculations including parameters taken in different algorithms are given in Appendix H.

Fig. 1 illustrates the result of the instance $(\mathrm{Heisenberg},\mathrm{chain},d=5)$, for example. We can find that the error decreases with the cost overhead. When the target error is relatively high, algorithms have similar performance, e.g. when $\epsilon\approx 0.01$, the cost $\gamma\approx 10$ is sufficient in most algorithms. However, in comparison with others, the error in the GP algorithm decreases much faster. When $d$ is sufficiently large, $M$ required to achieve the error $\epsilon$ always follows a scaling of $1/\epsilon^{2}$. However, the size of $d$ at which different algorithms converge to this behavior varies; see Fig. 9. The GP algorithm converges with a smaller $d$ compared to other algorithms. In this instance, the GP algorithm has already converged when $d=5$, whereas the other algorithms have not yet converged, resulting in poorer performance.

To confirm this advantage of the GP algorithm, we implement the simulation for a number of instances. We choose the power algorithm as the standard for comparison. We remark that the power algorithm and the classical Lanczos algorithm yield the same result when the implementation is error-free (i.e. without statistical error and rounding error), and the eventual error in this ideal case is the subspace error $\epsilon_{K}$. Given the model, lattice and subspace dimension, we compute $\epsilon_{K}$ of the power algorithm. Then, we take the permissible error $\epsilon=2\epsilon_{K}$ and compute the overhead factor $\gamma$ for each algorithm. Then, the empirical distribution of $\gamma$ is shown in Fig. 2.

From the distribution, we can find that our Gaussian-power algorithm has the smallest measurement overhead, and $\gamma$ is smaller than $10^{2}$ for almost all instances. The filter (F) algorithm is the most well-performed one among others. Part of the reason is that we have taken the optimal parameter found by grid search, of which the details are discussed in Appendix H. The median value of $\gamma$ is about $4$ for the GP algorithm, $13$ for the F algorithm and as large as $10^{12}$ for some other algorithms. Comparing the GP and F algorithms, although the median values are similar, the cost overhead is larger than $10^{4}$ in more than $20\%$ instances in the F algorithm, which never happens in the GP algorithm.

The performance of the GP algorithm depends on the choice of parameters: $\tau$, $E_{0}$ and $\eta$. The regularisation parameter $\eta$ is taken according to Table 2. To choose the value of $E_{0}$, we need prior knowledge about the true ground-state energy $E_{g}$, which can be obtained from certain classical algorithms computing the ground state or the quantum KSD algorithm taking a smaller subspace dimension. Although we prefer an $E_{0}$ close to $E_{g}$, the numerical result suggests that a relatively large difference $|E_{0}-E_{g}|$ ($10\%$ of the entire spectrum) is tolerable. The Gaussian factor in the GP basis acts as a filter on the spectrum. Therefore, if $|E_{0}-E_{g}|$ is finite and we take an over large $\tau$, we may have a vanishing overlap between the effective reference state $e^{-\frac{1}{2}(H-E_{0})^{2}\tau^{2}}|{\varphi}\rangle$ and the ground state. To avoid this problem, we need to assume an upper bound $|E_{0}-E_{g}|\leqslant\epsilon_{0}$, then we suggest taking $\tau$ in the range $\sqrt{\frac{d-1}{e}}<\tau\leqslant\frac{\sqrt{d-1}}{\epsilon_{0}}$, see Appendix K. Although we take the same $\tau$ for all $f_{k}$ in each instance in the benchmarking, one can flexibly choose different $\tau$ in practice.

We show that $\gamma$ is related to the spectrum of $\mathbf{S}$. Let’s consider a small $\eta$ and the Taylor expansion of $E^{\prime}$. The minimum value of the zeroth-order term $\frac{\mathbf{a}^{\dagger}\mathbf{H}\mathbf{a}}{\mathbf{a}^{\dagger}\mathbf{S}\mathbf{a}}$ is $E_{min}$ according to the Rayleigh quotient theorem [51]. Take this minimum value, we have $E^{\prime}\simeq E_{min}+s\eta$, where $s\propto\frac{\mathbf{a}^{\dagger}\mathbf{a}}{\mathbf{a}^{\dagger}\mathbf{S}\mathbf{a}}$ is the first derivative. The solution to $E^{\prime}=E_{g}+\epsilon$ is $\eta\simeq(E_{g}+\epsilon-E_{min})/s$, then

where $u\leqslant\frac{(C_{\mathbf{H}}+\|H\|_{2}C_{\mathbf{S}})\epsilon}{2\|H\|_{2}(\epsilon-\epsilon_{K})}\approx C_{\mathbf{S}}$ under assumptions $C_{\mathbf{H}}\approx\|H\|_{2}C_{\mathbf{S}}$ and $\epsilon-\epsilon_{K}\approx\epsilon$. Eq. (21) is derived in Appendix I. Therefore, when $\mathbf{S}$ is close to singular, the overhead can be large.

We can explain the measurement efficiency of our algorithm by composing a projector. In our algorithm with the rescaled basis, the solution is a state in the form $|{\Psi(\mathbf{a})}\rangle=\sum_{k=1}^{d}a_{k}c_{k}^{-1}f_{k}|{\varphi}\rangle$. Ideally, the linear combination of $f_{k}$ realises a projector onto the ground state, i.e. $\sum_{k=1}^{d}a_{k}c_{k}^{-1}f_{k}=|\psi_{g}\rangle\langle\psi_{g}|$. Then, $|{\Psi(\mathbf{a})}\rangle=\sqrt{p_{g}}|{\psi_{g}}\rangle$ up to a phase and $\mathbf{a}^{\dagger}\mathbf{S}\mathbf{a}=\langle\Psi(\mathbf{a})|\Psi(\mathbf{a})\rangle=p_{g}$. Then $\gamma\approx\left(\frac{p_{g}\mathbf{a}^{\dagger}\mathbf{a}}{\mathbf{a}^{\dagger}\mathbf{S}\mathbf{a}}\right)^{2}=\left(\mathbf{a}^{\dagger}\mathbf{a}\right)^{2}$ (notice that $C_{\mathbf{S}}=1$), we can find that $\mathbf{a}^{\dagger}\mathbf{a}$ determines $\gamma$. When $c_{k}$ is smaller, $|a_{k}|$ required for composing the projector is smaller. In our algorithm, $c_{k}$ decreases exponentially with $k$, and the speed of decreasing is controllable via the parameter $\tau$. In this way, our algorithm results in a small $\gamma$.

To be concrete, let’s consider a classic way of composing a projector from Hamiltonian powers using Chebyshev polynomials (CPs) $T_{n}$, which has been used for proving the convergence of the Lanczos algorithm [52], RTE and CP quantum KSD algorithms [16, 20]. Without loss of generality, we suppose $\|H\|_{2}=1$. An approximate projector in the form of Chebyshev polynomials is

where $Z=1-(H-E_{g}-\Delta)$, $\Delta$ is the energy gap between the ground state and the first excited state, $z_{1}=1+\Delta$, and $\Omega$ is an operator with the upper bound $\|\Omega\|_{2}\leqslant 2/(z_{1}^{n}+z_{1}^{-n})$. Notice that the error $\Omega$ depends on $n$, and its upper bound decreases exponentially with $n$ if $\Delta$ is finite. For simplicity, we focus on the case that $E_{0}=E_{g}$ in the Hamiltonian power. Then, in the expansion $T_{n}(Z)=\sum_{l=0}^{n}b_{l}(H-E_{g})^{l}$, coefficients have upper bounds $|b_{l}|\leqslant n^{l}T_{n}(z_{1})$ increasing exponentially with $l$. See Appendix J. In the LCU and Monte Carlo approach, large coefficients lead to large variance. Therefore, it is difficult to realise this projector because of the exponentially increasing coefficients.

In our algorithm, the exponentially decreasing $c_{k}$ can cancel out the large $b_{l}$. We can compose the CP projector with an additional Gaussian factor. Taking $d=n+1$ and $a_{k}=c_{k}b_{k-1}/T_{n}(z_{1})$, we have $\sum_{k=1}^{d}a_{k}c_{k}^{-1}f_{k}=\frac{T_{n}(Z)}{T_{n}(z_{1})}e^{-\frac{1}{2}(H-E_{g})^{2}\tau^{2}}$. Because the Gaussian factor is also a projector onto the ground state, the overall operator is a projector with a smaller error than the CP projector. The corresponding overhead factor is

When $\tau$ is sufficiently large, $\gamma\lesssim 4$. This example shows that the Gaussian-power basis can compose a projector with a small measurement cost.

By composing a projector with an explicit expression, we can obtain the worst-case performance of the algorithm. Similar analyses are reported in Refs. [16] and [20] for RTE and CP algorithms, respectively. Because the projector is approximate, it results in a finite error in the ground-state energy $\epsilon_{P}$. By solving the generalised eigenvalue problem, we can usually work out a better solution, i.e. the corresponding energy error is lower than $\epsilon_{P}$. Our result about the measurement cost is applicable to a target error that is not limited by $\epsilon_{P}$. Theorem 1 holds even when $\epsilon<\epsilon_{P}$, which is beyond the explicit projector analysis.

The method of constructing polynomial projection operators mentioned above can be applied not only to the GP algorithm but also to the P algorithm. Since the variance of the P algorithm does not decrease exponentially with $k$, it can be seen that our algorithm has an advantage over the P algorithm in this regard. We remark that for this particular projector, the decomposition into Hamiltonian powers is costly, and the existence of a better decomposition is possible. Compared to other algorithms, we have provided a numerical benchmark test in the previous section. Additionally, for the CP and RTE algorithms, it is possible to construct approximate projection operators with constant bounded expansion coefficients [16, 20]. In this case, further methods are needed to theoretically explain the improvement in measurement cost of our algorithm compared to these algorithms.

We can apply Theorem 1 to estimate the measurement cost of various quantum KSD algorithms for arbitrary target error $\epsilon$ and subspace dimension $d$. In this paper, we focus on a target error comparable to the Lanczos algorithm at the same $d$, i.e. $\epsilon\sim\epsilon_{K}$. The Lanczos algorithm is a well-known classical algorithm. It is natural to demand quantum KSD algorithms, as quantum counterparts of Lanczos algorithm, to reach a comparable target error. In this case, the measurement cost of the GP algorithm is observed much smaller than the others in the benchmarking.

In certain cases of practical applications, there is definite target error epsilon, e.g. the chemical accuracy. In these cases, these target errors are no longer depending on $\epsilon_{K}$. According to the Kaniel–Paige convergence theory [52], by increasing the subspace dimension $d$ and measurement number $M$, we can always achieve arbitrary finite target error. Our results focus on relatively small d, in which case the GP algorithm requires a much smaller $M$ than the others, as illustrated in Fig. 2. We show in Fig. 9 that $M$ converges to $O(1/\epsilon^{2})$ when $d$ is sufficiently large, as demonstrated in Refs. [20, 53]. Increasing $d$ can weaken the advantage of the GP algorithm in the measurement cost. However, increasing $d$ not only brings more measurement entries (i.e. a larger factor $\beta$), but also implies deeper quantum circuits and more gate noises.

In conclusion, we have proposed a regularised estimator of the ground-state energy (see Algorithm 1). Even in the presence of statistical errors, the regularised estimator is variational (i.e. equal to or above the true ground-state energy) and has a rigorous upper bound. Because of these properties, it provides a universal way of analysing the measurement cost in quantum KSD algorithms, with the result summarised in Theorem 1. This approach can be generalised to analyse the effect of other errors, e.g. imperfect quantum gates. The robustness to statistical errors in our algorithm implies the robustness to other errors.

To minimise the measurement cost, we have proposed a protocol for constructing Hamiltonian power with an additional Gaussian factor. This Gaussian factor is important because it leads to the statistical error decreasing exponentially with the power. Consequently, we can construct a Chebyshev-type projector onto the ground state with a small measurement cost. We benchmark quantum KSD algorithms with two models of quantum many-body systems. We find that our algorithm requires the smallest measurement number, and a measurement overhead $\gamma$ of a hundred is almost always sufficient for approaching the theoretical-limit accuracy $\epsilon_{K}$ of the Lanczos algorithm. In addition to the advantage over the classical Lanczos algorithm in scalability, our result suggests that the quantum algorithm is also competitive in accuracy at a reasonable measurement cost.

This paper primarily focuses on the ground-state problem. Notably, our quantum KSD algorithm holds promise for application to low-lying excited states as well. We leave it for future research.