Quantum Davidson algorithm for excited states

In linear algebra, an order-r Krylov subspace generated by a matrix H and a reference vector b is the linear subspace spanned by the images of b under the first r powers of H [62]. This subspace is denoted as $\mathcal{K}_r(H,b) \equiv \{b, Hb, H^2b, \cdots, H^{r-1}b\}$. The Krylov subspace has been extensively utilized in numerical algorithms for finding solutions to high-dimensional matrices, such as the generalized minimal residual method, Davidson, and quasi-minimal residual algorithms [61].

In quantum chemistry, the iterative Krylov subspace $\mathcal{K}_r(\hat{H},\vert \psi\rangle)$ is especially valuable for identifying low-lying states in electronic structure theory since the number of such states is typically much smaller than the size of the solution space. Various versions of the Davidson algorithm have been developed to enhance convergence by designing efficient pre-conditioners [63–65]. Furthermore, the Davidson algorithm has been generalized to use non-orthogonal Krylov spaces [66]. The flowchart of a standard classical Davidson algorithm is presented in algorithm S1 of the supplementary information (SI). Despite its computational efficiency, the generation of the subspace $\hat{H}^r\vert \psi\rangle$ remains a significant bottleneck in the Davidson algorithm and becomes computationally intensive on classical computers as the dimension of $\hat{H}$ grows.

In this study, we introduce a quantum counterpart of the Davidson algorithm, termed QDavidson, that harnesses QCs to mitigate the scaling challenges of generating the Krylov subspace. Compared to the traditional Krylov method, the QKS scheme encodes arbitrarily complex states on quantum circuits, where the matrix's projection into the subspace is measured. The benefits of QKS methods over VQE-based techniques for computing low-lying excited states include (1) an independent ansatz for each reference state with a streamlined quantum circuit and (2) the elimination of the intricate optimization process. In [56], quantum imaginary time evolution (QITE) is conducted using the Trotter decomposition of the evolution operator $\mathrm{e}^{-\hat{H} t}$ and mapping each Trotter step into a unitary. Subsequently, the QLanczos algorithm is introduced in the Krylov space derived from QITE snapshots. However, the QKS arising from sequential time evolution lacks compactness. Here, we develop the QDavidson algorithm to create a compact QKS set, yielding a more efficient and numerically stable realization of the QKS approach.

Instead of pre-generating the subspace through real or imaginary-time evolution, the QDavidson framework adaptively expands the QKS, keeping the subspace growth closely with the true eigenspace. Initially, orthogonal states serve as the reference. HF $\vert \phi_0\rangle$, single-excitation configurations $\vert \phi^a_i\rangle = a^\dagger_a a_i\vert \phi_0\rangle$, and CI single (CIS) states $\{\sum_{ia}C_{ia}\vert \phi^a_i\rangle\}$, which are efficiently computed on classical computers, are natural choices for initial reference states. The set of these initial reference states is denoted as $\{\vert \psi_k\rangle\}$. Within the QKS framework, each basis (subsequently called a Krylov vector) within the Krylov subspace of the jth Davidson iteration assumes the following ansatz,

$$\begin{equation} \vert \psi_\mathrm{K}\rangle = \hat{U} ^j_\mathrm{K} \sum_\mathrm{I} C_\mathrm{KI} \vert \phi_\mathrm{I}\rangle \equiv \hat{U} ^j_\mathrm{K} \vert \Phi_\mathrm{K}\rangle. \end{equation} \tag{ 1 }$$

Here, multiple reference states (i.e. a linear combination of HF/CIS states) are utilized where $\vert \phi_I\rangle$ denote a single determinant. The entanglers $\hat{U}^j_\mathrm{K}$ arise from the Krylov space expansion and introduce correlations that extend beyond the initial states. In the case of QDavidson, the entanglers take the form $\hat{U}^{j}_\mathrm{K} = \mathrm{e}^{-\mathrm{i}\hat{A}^{j}_\mathrm{K}\delta\tau}$, where the operator $\hat{A}^{j}_\mathrm{K}$ is determined by mapping the non-unitary operator $\mathrm{e}^{-(\hat{H}-E_\mathrm{K}^j)\delta\tau}$ into a unitary form. The further details of $\hat{U}^j_\mathrm{K}$ will be elaborated later, with $\hat{U}^1_\mathrm{K} = 1$ for the initial iteration.

The general ground and excited states, denoted $\vert \Psi_\mathrm{I}\rangle$, can be expressed as linear combinations of the Krylov vectors,

$$\begin{equation} \vert \Psi_\mathrm{I}\rangle = \sum_\mathrm{K} V_\mathrm{KI} \vert \psi_\mathrm{K}\rangle. \end{equation} \tag{ 2 }$$

Therefore, the challenge of finding the ground state and low-lying excited states, represented by the equation $\hat{H}\vert \Psi\rangle = E\vert \Psi\rangle$, can be recast as the generalized eigenvalue problem HV = ESV within the Krylov subspace. Here, H denotes the Hamiltonian matrix within the Krylov subspace, as described by equation (1),

$$\begin{equation} H_\mathrm{KL} = \langle \psi_\mathrm{K}\vert \hat{H} \vert \psi_L\rangle = \langle \Phi_\mathrm{K}\vert \hat{U}^\dagger_\mathrm{K} \hat{H} \hat{U} _\mathrm{L} \vert \Phi_\mathrm{L}\rangle. \end{equation} \tag{ 3 }$$

Additionally, S represents the overlap matrix among Krylov vectors, with individual elements given by,

$$\begin{align} S_\mathrm{KL} = \langle \psi_\mathrm{K} \vert \psi_\mathrm{L}\rangle = \langle \Phi_\mathrm{K}\vert \hat{U}^\dagger_\mathrm{K} \hat{U}_\mathrm{L} \vert \Phi_\mathrm{L}\rangle. \end{align} \tag{ 4 }$$

Considering each Krylov vector may contain distinct entanglers and reference states, the matrix elements of both H and S are measured on QCs using an ancillary qubit and the Hadamard test [51, 55], as illustrated in figure 1. Once the elements of H and S are measured, the generalized eigenvalue problem HV = ESV can be trivially solved on classical computers due to the small dimensionality of the Krylov subspace.

Zoom In Zoom Out Reset image size

After solving for the approximate excited states in the current subspace, the residues of these approximate states can be computed as

$$\begin{equation} \vert R_\mathrm{I}\rangle = \hat{H}\vert \Psi_\mathrm{I}\rangle - E_\mathrm{I} \vert \Psi_\mathrm{I}\rangle = \left(\hat{H} - E_\mathrm{I}\right) \sum_\mathrm{K} V_\mathrm{KI} \vert \psi_\mathrm{K}\rangle. \end{equation} \tag{ 5 }$$

The norm of the residue can be measured on QCs similarly to the $H_\mathrm{KL}$ measurement,

$$\begin{align} |\vert R_\mathrm{I}\rangle| = &\langle \Psi_\mathrm{I}\vert \left(\hat{H} - E_\mathrm{I}\right)^2 \vert \Psi_\mathrm{I}\rangle \nonumber\\ = &\sum_\mathrm{KL} V\,^*_\mathrm{KI} V_\mathrm{LI} \langle \psi_\mathrm{K}\vert \left(\hat{H} - E_\mathrm{I}\right)^2 \vert \psi_\mathrm{L}\rangle. \end{align} \tag{ 6 }$$

Hence, the measurement of $|\vert R_\mathrm{I}\rangle|$ is akin to the $H_\mathrm{KL}$ measurement but with the Hamiltonian $\hat{H}$ in equation (3) replaced by $(\hat{H}-E^2_\mathrm{I})^2$. The norm indicates the convergence of each eigenstate. If $|\vert R_\mathrm{I}\rangle|$ exceeds ε (where ε denotes the convergence criterion), a new Krylov vector, based on the normalized counterpart of $\vert R_\mathrm{I}\rangle$, should be incorporated into the Davidson algorithm, provided it is linearly independent within the existing Krylov subspace.

The classical Davidson employs the Krylov subspace $\mathcal{K}_r((\hat{H}-E),\vert \Phi_\mathrm{I}\rangle)$ to solve the eigenvalue problem. However, it is non-trivial to create the $\hat{H}\vert \Psi_\mathrm{I}\rangle$ state on quantum circuits with the non-unitary $\hat{H}$. As an alternative, within the QDavidson algorithm, a correction vector defined by

$$\begin{equation} \vert \delta_\mathrm{I}\rangle = \mathrm{e}^{-\Delta\tau \left(\hat{H}-E_\mathrm{I}\right)}\vert \Psi_\mathrm{I}\rangle \end{equation} \tag{ 7 }$$

is employed to expand the Krylov subspace. In other words, the QDavidson algorithm uses the subspace of $\mathcal{K}_r(\mathrm{e}^{-\Delta\tau(\hat{H}-E)},\vert \Phi_\mathrm{I}\rangle)$ to derive the eigenstates. Because $\vert \delta_\mathrm{I}\rangle = [1-\Delta\tau(\hat{H}-E)]\vert \Phi_\mathrm{I}\rangle+\mathcal{O}(\Delta\tau)$, it can be verified that $\mathcal{K}_r(\mathrm{e}^{-\delta\tau(\hat{H}-E)},\vert \Phi_\mathrm{I}\rangle)$ is equivalent to $\mathcal{K}_r((\hat{H}-E),\vert \Phi_\mathrm{I}\rangle)$ when $\Delta\tau\rightarrow 0$. The evolution detailed in equation (7) is subsequently mapped to unitary operators $\mathrm{e}^{-\mathrm{i}\hat{A}}$ as proposed in [56],

$$\begin{equation} \vert \delta_\mathrm{I}\rangle \simeq n_\mathrm{I} \mathrm{e}^{-\mathrm{i}\hat{A}}\vert \Phi_\mathrm{I}\rangle \end{equation} \tag{ 8 }$$

where $n_\mathrm{I}$ is the normalization factor and $\hat{A} = \sum_\alpha a_\alpha\hat{P}_\alpha$. The coefficients a_α are obtained by solving a system of linear algebra equations associated with the mapping (more details can be found in appendix A). An alternative method involves directly mapping the preconditioned residue $\frac{1}{\bar{H}-E_i}(\hat{H}-E)$ into $\mathrm{e}^{-\mathrm{i}\hat{A}}$, represented by $\vert \delta_\mathrm{I}\rangle = \frac{1}{\bar{H}-E_i}(\hat{H}-E) \vert \Phi_\mathrm{I}\rangle \approx C\mathrm{e}^{-\mathrm{i}\hat{A}} \vert \Phi_\mathrm{I}\rangle$ where $\bar{H}$ is the diagonal part of the Hamiltonian matrix in the Krylov subspace. We note that the investigation of the effects of different preconditioners on QDavidson's performance is beyond the scope of the current paper.

After mapping the residue operator to unitaries, the QDavidson algorithm determines if the new Krylov vectors (or the correction vector) are linearly dependent within the current subspace by introducing $\vert \delta^\prime_{\mathrm{K}^{\prime}}\rangle$

$$\begin{align} \vert \delta^{\prime}_\mathrm{I}\rangle&\equiv \vert \delta_\mathrm{I}\rangle-\sum_\mathrm{KJ} \vert \Psi_\mathrm{K}\rangle\left(S^{-1}\right)_\mathrm{KJ}\langle \Psi_\mathrm{J}\vert\delta_\mathrm{I}\rangle\nonumber\\ &\equiv \vert \delta_\mathrm{I}\rangle - \sum_\mathrm{KJ} \left(S^{-1}\right)_\mathrm{KJ} d_\mathrm{IJ} \vert \delta_\mathrm{I}\rangle. \end{align} \tag{ 9 }$$

Here, $d_\mathrm{KJ} = \langle \Psi_\mathrm{J}\vert\delta_\mathrm{I}\rangle$ and the norm $|\vert \delta^{\prime}_\mathrm{I}\rangle|$ is measured is measured on QCs. If the ratio $\frac{|\vert \delta^{\prime}_\mathrm{I}\rangle|}{|\vert \delta_\mathrm{I}\rangle|}\gt\epsilon$, it implies that the new Krylov vector $\vert \delta_\mathrm{I}\rangle$ is linearly independent within the current subspace, and hence, it should be incorporated into the subspace. As the inclusion of a new Krylov vector into the Krylov space introduces additional correlations, the next iteration will draw the results nearer to the exact solutions. The flowchart of the QDavidson algorithm is summarized in figure 2.

Zoom In Zoom Out Reset image size

Since the size of the Krylov subspace $N_\mathrm{K}$ is small, computing the generalized eigenvalue problems on classical computers is cheap. The primary complexity of the QDavidson algorithm stems from mapping the residue vectors into unitaries. Hence, its main bottleneck is the formation of S, b, and the resolution of the linear system to obtain unitary $\hat{A}$. Earlier research has shown that mapping non-unitary exponential operators into unitaries scales exponentially with correlation domain D [56]. But, a local approximation can be applied to eliminate the exponential dependence on D [56], leading to polynomial complexity, which is also utilized in the QDavidson method.

Both long-range (LR) and short-range (SR) 1D Heisenberg models, analogous to the models described in [56], were tested in this work. The SR models account only for nearest-neighbor interactions between spins, characterized by $C_{ij}(\hat{X}_i\hat{X}_j+\hat{Y}_i\hat{Y}_j+\hat{Z}_i\hat{Z}_j)$ terms. In contrast, the LR models consider pairwise interactions among all spins, encompassing a larger number of terms in a qubit Hamiltonian. Explicitly, the LR and SR Hamiltonians are given by

$$\begin{equation} \hat{H}_\mathrm{SR} = -\sum_{i = 1}^{N}{\hat{S}_{i} \hat{S}_{i+1}}, \end{equation} \tag{ 10 }$$

and

$$\begin{equation} \hat{H}_\mathrm{LR} = -\sum_{j>i}^{N}{\frac{1}{D_{ij}}\hat{S}_{i} \hat{S}_{j}}, \end{equation} \tag{ 11 }$$

where $D_{ij} = \min(1+|j-i|, 1+N-|j-i|)$. N denotes the number of spins in the system. For the SR Hamiltonians, the index i is cyclic; thus, when it reaches N + 1, it reverts to 1. The QDavidson algorithm is state-dependent, necessitating the definition of initial reference states. We initiated the algorithm with the anti-ferromagnetic product state, which corresponds to the alternating $\vert 01..\rangle$ state in a computational basis. In order to benchmark the QDavidson algorithm against other QKS methods, a QLanczos algorithm with identical initial parameters was also implemented.

The results of the QDavidson and QLanczos algorithms for the 1D-Heisenberg models are shown in figure 3. For the initial four low-lying solutions, the QDavidson algorithm achieved convergence within 4 and 7 iterations for four-spin and six-spin systems, respectively (figures 3(E)–(H)). This performance surpasses the recently proposed QLanczos [56] algorithm (figures 3(A)–(D)), which not only converges slower but also displays numerical instability as the resultant states quickly become linearly dependent. By implementing root convergence criteria and a linear dependency check, the QDavidson algorithm exhibits enhanced numerical stability. The full results for the QLanczos algorithm are provided in the SI (figure S3). This algorithm was executed until an iteration produced a linearly dependent state. Notably, while the QLanczos managed to compute exact energy values for all four-spin models within six iterations (figures S3(A) and (B)), it did not converge for six-spin systems, and the final energy estimations showed a considerable deviation (figures S3(C) and (D)).

Zoom In Zoom Out Reset image size

As each iteration appends a new entangler after mapping the imaginary time evolution operator into unitaries, the circuit depth increases monotonically with iterations within the QITE, QLanczos, and QDavidson frameworks. Hence, faster convergence translates to more compact circuits. Compared to QITE algorithms, QLanczos converges faster in finding the ground state [56]. Our QDavidson algorithm substantially improves convergence by employing the residue operator to narrow the subspace search near the exact state, resulting in a significantly reduced circuit depth. Although this improvement may not be evident for smaller systems, it becomes remarkably advantageous for larger systems. To elucidate this, we examined how the lowest-state energy's accuracy varies with the resulting circuits' maximum depth for both QLanczos and QDavidson. The results are shown in figure 4. QDavidson reproduces the exact solution for six-spin systems when the subspace comprises circuits with maximum gate depths of approximately 400 and 900 for the SR and LR models, respectively. In contrast, QLanczos results are less accurate, and deeper circuits are required to perform the algorithm. For instance, the accuracy attained by QLanczos after eight iterations (gate depth of 480) for the six-spin SR model can be achieved by QDavidson using circuits with a peak gate depth of just 120. Given the iterative expansion of the subspace, each QDavidson algorithm iteration contains states characterized by circuits of various depths. Consequently, tracking solely the longest circuit within the subspace is sufficient. The circuit depth was analyzed based on the Pauli words present in the operator $\hat{A}$ in the $\mathrm{e}^{-\Delta\tau (\hat{H}-E)}\vert \Psi_I\rangle \rightarrow \mathrm{e}^{-\mathrm{i}\hat{A}}\vert \Psi_I\rangle$ expansion. Recognizing the elements of $\hat{A} = \sum_i \theta_i \hat\sigma_i$, a standard unitary evolution circuit could be constructed for $\mathrm{e}^{-\mathrm{i}\hat{A}}$ (assuming one Trotter step).

Zoom In Zoom Out Reset image size

In addition to the 1D Heisenberg models, we also examined our algorithm on chemical systems. As expected, molecular Hamiltonians are notably more complex than the 1D Heisenberg models due to the intricate entanglement between numerous orbitals. To highlight the advantage of the QDavidson algorithm for molecular systems, we studied three molecular systems: ethylene, cyclopropene cation, and benzene. The Cartesian coordinates of the investigated molecular systems are provided in the SI (see table S1). An active space representing the π-conjugated system was selected for all the systems. The detailed active orbitals can be found in the SI (see figure S1). Therefore, the (2e, 2o) active space was considered for the $\mathrm{C_2H_4}$ molecule, (2e, 3o) for $\mathrm{C_3H_3^+}$, and (6e, 6o) for $\mathrm{C_6H_6}$. The Jordan–Wigner transformation [67] was utilized to map the second-quantized Hamiltonian into the qubit Hamiltonian. Consequently, the largest system examined in this study is the benzene molecule, comprising six electrons and 12 spin-orbitals in the active space (which corresponds to a 12 qubit Hamiltonian with 407 terms).

Taking into account that the electronic structure problems in molecular systems maintain particle-preserving and total spin-preserving symmetries, it is beneficial to establish a symmetric pool of Pauli terms for the correction vector $\vert \delta_\mathrm{I}\rangle$ mapping. Therefore, a set of all possible spin projection-preserving single and double excitation operators was selected and then transformed into qubit operators employing the Jordan–Wigner (JW) scheme [67]. All unique Pauli terms with an odd number of $\hat{Y}$ operators were selected into the pool. This results in 12, 40, and 828 unique Pauli terms for the ethylene, cyclopropene cation, and benzene, respectively.

The performance of the QDavidson algorithm depends on the defined initial state. Multiple simulations were conducted using varying initial configurations ($\vert \phi_0\rangle, \vert \phi^a_i\rangle = a^\dagger_a a_i\vert \phi_0\rangle$). In particular, the HF, the first excited singlet, and the first triplet product states were taken into consideration. The QDavidson algorithm converges to exact energy values for the smallest molecular systems within several iterations. For $\mathrm{C_2H_4}$, starting from the HF product state, the QDavidson algorithm obtained the exact ground state energy (−77.115 18 Hartree) and the exact second excited singlet state energy (−76.375 41 Hartree) after just one iteration. Similar results were acquired for the $\mathrm{C_3H_3^+}$ system, where the exact ground state (−113.649 29 Hartree) and the first three singlet excited states were determined. Table 1 presents results for other initial states. Intriguingly, the eigenstates identified for the $\mathrm{C_3H_3^+}$ system maintain particle preservation and do not coincide with a neutral state solution with three electrons instead of two, even though they are lower in energy.

Table 1. Results of QDavidson algorithm for $\mathrm{C_2H_4}$ and $\mathrm{C_3H_3^+}$ molecules. The initial state and number of iterations for the algorithm convergence are given. Different eigenstates are denoted in parenthesis as the following: S₀—ground singlet state, T₁—the lowest triplet state, S_i —ith excited singlet state, T_i —ith excited triplet state. The S_z denotes the spin projection of the obtained state.

Molecule	Initial	Number of	Energies
	state	iterations
$\mathrm{C_2H_4}$	$\vert 0011\rangle$	1	−77.115 18 (S0), Sz = 0
			−76.375 41 (S2), Sz = 0
	$\vert 1001\rangle$	1	−76.922 18 (T1), Sz = 0
			−76.582 76 (S1), Sz = 0
	$\vert 0101\rangle$	0 a	−76.922 18 (T1), Sz = 1
$\mathrm{C_3H_3^+}$	$\vert 000011\rangle$	2	−113.649 29 (S0), Sz = 0
			−113.198 54 (S1), Sz = 0
			−112.759 64 (S2), Sz = 0
			−112.681 41 (S3), Sz = 0
	$\vert 001001\rangle$	2	−113.352 09 (T1), Sz = 0
			−113.198 54 (S1), Sz = 0
			−112.968 18 (T2), Sz = 0
			−112.759 64 (S2), Sz = 0
	$\vert 000101\rangle$	0 a	−113.352 09 (T1), Sz = 1

^a The initial product state is already an exact solution.

For the more extensive $\mathrm{C_6H_6}$ system, the QDavidson algorithm also accurately located the energies of several low-lying states (within chemical accuracy). However, as depicted in figures 5(A)–(C), many more iterations were needed. This is presumably due to the limited operator pool chosen, which only includes Pauli words from JW-mapped single and double excitation operators. Consequently, multi-electron systems (with electron count exceeding 2) might lack sufficient flexibility in the $\mathrm{e}^{-\Delta\tau (\hat{H}-E)}$ operator mapping. Regardless, the algorithm yielded meaningful results even with this restricted operator pool. For the benzene molecule, the algorithm found the first excited singlet and triplet states within 10⁻² Hartree after approximately 30 iterations (see figures 5(B) and (C)). For comparison, the QLanczos algorithm tends to get trapped in linearly dependent solutions after approximately 30 iterations, and the precision of the energies obtained is inferior to that in the QDavidson case. Consequently, when initiating the algorithm with the $\vert \mathrm{HF}\rangle$ state, the energy error for the ground state was found to be $5\times10^{-3}$ Hartree, while the first excited singlet and triplet states could not be accurately reproduced, as illustrated in figure S2. With initialization with the excited state determinant, the algorithm successfully reproduced the energy of the first excited triplet within 10⁻¹ Hartree. However, this accuracy is an order of magnitude lower than that achieved with QDavidson.

Zoom In Zoom Out Reset image size

For larger systems, it becomes evident that initiating the algorithm from a singular initial state is not the optimal strategy to capture the first few low-lying excited states, as illustrated in figures 5(A)–(C). Even though the algorithm can locate states of interest, the accuracy is far from the desired precision of 10⁻³ Hartree. Even initiating the algorithm with the first excited singlet (figure 5(F)) or triplet (figure 5(G)) determinants leads to quicker convergence compared to using the HF initial state. However, the energy convergence remains slow and often diverges from the exact values. Additionally, the initialization of the QDavidson with the entangled triplet CIS state was performed. Although the initial rapid convergence toward the exact excited energy was observed as illustrated in figure S3 (the accuracy of 10⁻² Hartree was achieved already at the third iteration), the convergence slows down and the desired accuracy of 10⁻³ Hartree was not achieved even after 30 iterations. In turn, the generalization of the algorithm with multiple initial states results in faster convergence and higher accuracy (see figure 5(H), highlighting the significance of multireference states. Notably, chemically accurate energies for S₀, T₁, and T₂ states could be obtained after just 11 iterations. Therefore, to efficiently find the low-lying excited states of a large molecular system, it is advised to begin with multiple initial states to expedite energy convergence. Nevertheless, it is anticipated that the convergence will be affected by different initial states, and we expect a more significant overlap between the initial states and true ground state will accelerate the convergence as is in classical Davidson algorithms.

Owing to the finite number of quantum measurements (shots), exact measurements of the Hamiltonian matrix and overlap matrix elements within the Krylov subspace are unattainable. A limited shot count may introduce numerical instabilities into the QDavidson procedure. In general, to estimate the expectation value of a Hamiltonian with respect to a single state $\vert \psi\rangle$ with precision p, it is required to perform $O(|h_{\max}|^2 Mp^{-2})$ measurements, where $h_{\max}$ is the largest coefficient in the Hamiltonian decomposition $\hat{H} = \sum_i{h_i\hat{P}_i}$, and M denotes the number of Hamiltonian terms [18]. Similar conclusions apply when evaluating Hamiltonian matrix elements within the Krylov subspace. Since $\langle \Psi_\mathrm{I}\vert \hat{H} \vert \Psi_\mathrm{J}\rangle = \sum_i{h_i\langle \Psi_\mathrm{I}\vert \hat{P}_i \vert \Psi_\mathrm{J}\rangle} = \sum_i\sum_\mathrm{KL}V^*_\mathrm{KI}V_\mathrm{LJ}$ $\langle \psi_K\vert \hat{P}_i\vert \psi_\mathrm{L}\rangle$ and each element $\langle \psi_\mathrm{K}\vert \hat{P}_i\vert \psi_\mathrm{L}\rangle$ can be evaluated according to figure 1, the cost of evaluating the Hamiltonian matrix element is $O(|h_{\max}|^2 Mp^{-2}N^2_\mathrm{K})$ where N_K is a size of the Krylov subspace. Although this measurement process can be complex, in practice, only a small subspace is typically spanned, and the growth of the subspace is linear with respect to the number of iterations, making the $N^2_\mathrm{K}$ factor relatively insignificant compared to $|h_{\max}|^2 Mp^{-2}$. Likewise, the overlap matrix's measurement incurs a cost of $O(p^{-2}N^2_K)$ shots. The procedure for mapping residual vectors to unitaries also requires quantum measurements. The required shot count here depends on the size of the operator pool (P) chosen for the mapping and is in order of $\max[O(|h_{\max}|^2 Mp^{-2}N^2_\mathrm{K} P),O(p^{-2}N^2_\mathrm{K} P^2)]$. Another procedure necessitating quantum measurements involves checking for root convergence and assessing linear dependency. The root convergence check, as introduced in equation (6), incurs a cost of $O(|h_{\max}|^4 M^2p^{-2}N^2_\mathrm{K})$ measurements. This is considered affordable since size of Krylov subspace $N_\mathrm{K}$ is small and M grows polynomially with the system size. For linear dependency, if it is realized through the ratio $\frac{|\vert \delta^{\prime}_\mathrm{I}\rangle|}{|\vert \delta_\mathrm{I}\rangle|}\gt\epsilon$ the number of measurements required would be $O(p^{-2}(N_\mathrm{K}+1)^2)$. Alternatively, linear dependency can be assessed by computing the determinant of the new overlap matrix (checking $\det(S) \lt \epsilon$) obtained by including $\vert \delta^{\prime}_\mathrm{I}\rangle$ in the Krylov subspace, with the same number of measurements. It is noteworthy that QLanczos formally requires fewer measurements per iteration as it does not involve checking for root convergence and linear dependency. However, its performance is significantly hampered by numerical instability, which might be deteriorated with the presence of shot noise. Besides, according to our numerical comparison, the QLancocs require more iterations to converge, eventually leading to more measurements. Additionally, it is noteworthy that the shots needed for QDavidson to converge are fewer than for the VQE procedure since the latter's ansatz optimization involves a substantial number of energy evaluations.

To illustrate the number of shots required to perform the algorithm on a real-world example, we incorporated the shot-noise into the calculations for the $\mathrm{C_2H_4}$ molecule. The desired precision p was determined after examining the algorithm's noise robustness by introducing the random error to each element of Hamiltonian matrix, (off-diagonal) overlap matrix element, S_αβ and b_α (appendix A). It was found that an evaluation precision of 10⁻⁴ is sufficient to guarantee the algorithm's numerical stability. With this level of precision, energy values of $-77.115\pm0.001 (S_0)$ and $-76.375\pm0.001 (S_2)$ were replicated. The algorithm, when initiated with a HF reference state, converged after one iteration. The estimated shot count per circuit evaluation was set to 10⁸, and the total number of shots was found to be ${\sim}10^{10}$. This experiment underscores the QDavidson algorithm's capability to reproduce the low-lying spectra of the evaluated Hamiltonian even in the presence of statistical shot noise.