Constructing Local Bases for a Deep Variational Quantum Eigensolver for Molecular Systems

The deep VQE [13] is an algorithm of the divide and conquer family. Its purpose is to calculate the ground state of a quantum system. We review the algorithm of the original deep VQE, which can treat spin Hamiltonians consisting of 2-local interactions. In deep VQE, we first divide the system into $M$ subsystems. The problem Hamiltonian with 2-local interactions can consequently be written in the form of

where $H_{i}$ acts on subsystem $i$ and $V_{ij}$ on subsystems $i$ and $j$. The interaction term $V_{ij}$ decompose into operators $V_{\mu,i}$ acting only on single subsystems $i$:

where $\mu$ indicates the different interaction terms. The deep VQE first finds ground states $\ket{G_{i}}$ of each subsystem Hamiltonian $H_{i}$ using a VQE algorithm. Then we build a subsystem basis $\{\ket{b_{i,k}}\}_{k=1}^{K_{i}}$ for each subsystem $i$ with dimensions $K_{i}$ by acting with certain excitation operators $\mathcal{B}_{i}=\{B_{i,k}\}$ on $\ket{G_{i}}$. The selection of the excitation operators $B_{i,k}$ is a crucial step which determines the performance of the deep VQE, and we will discuss it in detail in Sec. III. To ensure the created basis set $\{\ket{b_{i,k}}\}_{k=1}^{K_{i}}$ is orthogonal a Gram-Schmidt orthogonalization is applied.

In the next step, we construct an effective Hamiltonian $H^{\mathrm{eff}}$ by projecting the original Hamiltonian $H$ to the subspace spanned by $\bigotimes_{i}\{\ket{b_{i,k}}\}_{k=1}^{K_{i}}$. This can be achieved by measuring all matrix elements $\bra{b_{i,k}}H_{i}\ket{b_{i,l}}$ and $\bra{b_{i,k}}\bra{b_{j,l}}V_{i,j}\ket{b_{i,p}}\ket{b_{j,q}}=\sum_{\mu}\lambda^{\mu}_{i,j}\bra{b_{i,k}}V_{\mu,i}\ket{b_{i,p}}\bra{b_{j,l}}V_{\mu,j}\ket{b_{j,q}}$. Notice that calculating the expectation values only involves one subsystem, which can be calculated separately on a quantum computer. Now, we represent each subsystem by $N_{i}=\lceil\log_{2}K_{i}\rceil$ qubits, and we again utilize the VQE using $N_{\mathrm{tot}}=\sum_{i}\lceil\log_{2}K_{i}\rceil$ qubits to search for the ground state of $H^{\mathrm{eff}}$. Even though $H^{\mathrm{eff}}$ whose dimension is $d=\prod_{i}K_{i}$ could, in principle, be mapped to $\lceil\log_{2}d\rceil\leqslant N_{\mathrm{tot}}$ qubits, it is beneficial to map each subsystem individually since it preserves the locality of the Hamiltonian. The deep VQE algorithm can be repeated iteratively to treat larger and larger systems on a qubit-limited quantum computer.

The central problem in quantum chemistry is finding the first-principles Hamiltonian’s ground state, which describes interacting electrons. In second quantization, the Hamiltonian can be written in the form

where $N$ is the number of orbitals of the molecule. Here we consider how to apply the deep VQE to this Hamiltonian.

A frequent approach to treat the Hamiltonian of Eq. (3) on a quantum computer is to map the fermionic operators $a_{o}^{\dagger}$ and $a_{o}$ to qubit operators through, e.g., Jordan-Wigner transformation. This leads to a Hamiltonian in the form of,

where $P_{i}$ is a Pauli string, $P_{i}\in\{I,X,Y,Z\}^{\otimes N}$, and $h_{i}$ is a real coefficient. The obstacle to apply the deep VQE is that $P_{i}$ is not 2-local but can be nonlocal when using JW transformation as we will see below.

After partitioning the system into $M$ subsystems and $M_{\mathrm{int}}$ interaction terms, we can write the Hamiltonian in Eq. (4) in the form of,

where $H_{i}$ and $V_{\mu,i}$ is an operator acting only on the $i$-th subsystem, and $\mu$ indexes different interaction terms. To perform the deep VQE, we construct the effective Hamiltonian $H^{\mathrm{eff}}$ using the subsystem basis $\{\ket{b_{i,k}}\}_{k=1}^{K_{i}}$. To do this, we measure the matrix elements $(H_{i}^{\mathrm{eff}})_{k,l}=\bra{b_{i,k}}H_{i}\ket{b_{i,l}}$ and $(V_{\mu,1}^{\mathrm{eff}})_{k,l}=\bra{b_{i,k}}V_{\mu,i}\ket{b_{i,l}}$ for all combinations of $i$, $k$, $l$, and $\mu$. Let $H_{i}^{\mathrm{eff}}$ and $V_{\mu,i}^{\mathrm{eff}}$ be $n_{i}=\lceil\log_{2}K_{i}\rceil$-qubit operators with the evaluated matrix elements. The total effective Hamiltonian can be written as,

If we wish to perform the VQE of $H^{\mathrm{eff}}$, we must be able to efficiently measure the expectation value $\langle{H^{\mathrm{eff}}}\rangle=\bra{\psi}H^{\mathrm{eff}}\ket{\psi}$ for a given state $\ket{\psi}$ prepared on a quantum computer. A usual approach to measure the expectation value of a Hamiltonian is to expand it into Pauli operators. However, if we do so in Eq. (6), the number of Pauli operators grows exponentially to $M$ due to its nonlocality, which blocks efficient evaluation of $\langle{H^{\mathrm{eff}}}\rangle$.

We must avoid this exponential growth to apply the deep VQE to fermionic systems. To this end, we first diagonalize $V_{\mu,i}^{\mathrm{eff}}$ classically and obtain unitary $U_{\mu,i}$ such that

where $\bm{\upsilon}_{\mu,i}$ are eigenvalues of $V_{\mu,i}^{\mathrm{eff}}$. Note that this process can be performed in time $2^{O(n_{i})}$ and that it is natural to assume $n_{i}$ is a small constant in the deep VQE. We can also construct a quantum circuit to realize a $2^{n_{i}}$-dimensional unitary $U_{\mu,i}$ in time $2^{O(n_{i})}$. Therefore, we can measure $\langle{V_{\mu,1}^{\mathrm{eff}}\otimes\cdots\otimes V_{\mu,M}^{\mathrm{eff}}}\rangle$ for a given $\sum_{i=1}^{M}n_{i}$-qubit state $\ket{\psi}$ by first applying $U_{\mu,1}\otimes\cdots\otimes U_{\mu,M}$ to $\ket{\psi}$ and then measuring it in the computational basis. Therefore, the required measurements to determine the expectation value of the interaction term scale linearly with the number of interactions in the system. Thereby, we can efficiently measure $\langle{V_{\mu,1}^{\mathrm{eff}}\otimes\cdots\otimes V_{\mu,M}^{\mathrm{eff}}}\rangle$ and hence $\langle H^{\mathrm{eff}}\rangle$ as well.

Note that we prefer JW transformation for applying deep VQE to fermionic systems because it directly maps the $i$-th orbital to the $i$-th qubit; if the $i$-th qubit is $\ket{1}$, it indicates that an electron occupies the $i$-th orbital. This property allows us to split systems naturally into subsystems using localized orbitals. Although other techniques such as Bravyi-Kitaev (BK) transformation relaxes the locality of Pauli strings down to $O(\log N)$, it makes the correspondence between fermion occupation and qubit configurations not straightforward. Note that the above technique may also be used for BK-transformation-based deep VQE to obtain the expectation values efficiently.

Choosing a basis set impacts the performance of the deep VQE. It ultimately determines the number of qubits needed to run the algorithm and the accuracy of the result. It is therefore vital to make an educated decision for the used method.

Below we examine different procedures to generate the basis sets, which are summarized in Tab. 1 as well as their scaling in Table 2. In the single Pauli method, we use the set of the Pauli operators $X,Y,Z$ acting on each qubit as $\mathcal{B}_{i}$. In the interactions’ approach, we apply all $V_{\mu,i}$ in Eq. (6) to the subsystem $i$ to create a basis, this leads to $\mathcal{B}_{i}=\{V_{\mu,i}\}_{\mu=1}^{M_{\mathrm{int}}}$. For the particle conserving approach, we apply the excitation $a_{s}^{\dagger}$ or the de-excitation operators $a_{s}$ to every qubit $s$ of the $n_{i}$ qubits of the subsystem. Additionally, to capture effects inside the subsystems, we applied swap gates between different qubits belonging to the same subsystem to create a basis. Therefore, in this approch, $\mathcal{B}_{i}=\{a_{s}\}_{s=1}^{n_{i}}\cup\{a_{s}^{\dagger}\}_{s=1}^{n_{i}}\cup\{\mathrm{SWAP}_{s,s^{\prime}}\}_{s,s^{\prime}=1}^{n_{i}}$.

We also use additional low-lying energy eigenstates $\ket{G^{e}_{i}}$ to generate the basis sets instead of using only $\ket{G_{i}}$ as has been done in the original paper [13]. The superscript $e$ marks the $e$-th excited state and we define $\ket{G_{i}^{0}}:=\ket{G_{i}}$. Note that if the subsystem Hamiltonian has degenerate ground states, $\ket{G^{e}_{i}}$ does not necessarily have different energy than $\ket{G_{i}^{0}}$. Such excited states can be constructed with an algorithm like the subspace-search VQE [17]. After applying $\mathcal{B}_{i}$ to $\{\ket{G^{e}_{i}}\}_{e=1}^{l}$, we use the Gram-Schmidt algorithm to ensure that the created basis is orthogonal and minimal with $K_{i}$ elements.

The number of basis vectors determines how many qubits we need to run the deep VQE algorithm. We, therefore, also examine different approaches to reduce the dimensionality of the basis. For the interaction-based methods, we only consider interactions between subsystems with $\lambda^{\mu}$ exceeding a certain threshold $\varepsilon>0$. In the numerical experiments that follow, we examine the results for $\varepsilon=10^{-2}$ if not stated otherwise. This selection was necessary since otherwise, the interaction-based methods would lead to basis sets that resemble no reduction from the entire Hilbert space. This was caused due to the numerous interaction terms between the subsystems for the molecular Hamiltonian.

The particle conserving and single Pauli methods as well can produce basis sets with a dimensionality that is challenging to represent on a quantum machine. We, therefore, also examine methods to truncate their basis sets. Our approach aimed to reduce the number of participating qubits in the subsystems for the basis-creating step leading to a subset of the entire basis. To do so, we select the interaction $\bigotimes_{i=1}^{M}\lambda^{\mu}V_{\mu,i}$ which has the largest absolute $\lambda^{\mu}$. Let us denote the set of qubits involved in this interaction by $A_{i}$ for each subsystem $i$. We then apply the corresponding excitation operators of the particle conserving or single Pauli methods to the qubits in $A_{i}$ to create the basis. We call these methods particle conserving edges and single Pauli edges, respectively. Mathematically, the former method uses $\mathcal{B}_{i}=\{a_{s}\}_{s\in A_{i}}\cup\{a_{s}^{\dagger}\}_{s\in A_{i}}\cup\{\mathrm{SWAP}_{s,s^{\prime}}\}_{s,s^{\prime}\in A_{i}}$ and the latter uses $\mathcal{B}_{i}=\{X_{s},Y_{s},Z_{s}\}_{s\in A_{i}}$. The dimensionalities $K_{i}$ for the different strategies scale differently with increasing system sizes, which are summarized in Table 2. These method will be compared numerically in Sec. IV.3.

The above strategies can reduce the number of qubits but in a hardware-agnostic way. In practice, a quantum computer has a fixed number of qubits, and we wish to use as many qubits as possible within the hardware limitation. With this in mind, we also propose an adaptive interaction-based strategy to create a basis set. We first order the interactions between the subsystems based on the interaction strength. By gradually increasing the threshold $\varepsilon_{adapt}$, we find a basis set of the desired dimensionality. This selection of thresholds can be performed for each subsystem independently. We will present the numerical demonstration of this method in Sec. IV.4.

A subsystem Hamiltonian having degenerate ground states can lead to unexpected results. A VQE then cannot find a unique ground state. The found ground state would be semirandom and depends on multiple factors, such as the starting condition of the VQE and the noise in the quantum device. In principle, all linear combinations of the degenerate ground states are possible solutions to the VQE without any treatment. This results in an instability of the deep VQE since the random starting vectors, in general, lead to different basis sets and consequently to different $H^{\mathrm{eff}}$.

This dependency on the starting vector set is tricky as, without supplementary knowledge about the system, there is no justification for favoring one subvector set over the other without running it through the deep VQE first. Therefore, the deep VQE must consider all degenerate vectors as starting vectors to ensure that the result is independent of the randomness of the first VQE. Alternatively, we should design the first VQE to return a unique solution for the reproducibility of the experiment by, for example, unfolding the degeneracy with constrained VQE [18].

First, we discuss the dependency of the deep VQE on the starting vector. For the sake of readability, we only show the results using the particle conserving method to produce the basis in the deep VQE. We show the results for the 10-atom tree molecule in Fig. 4. The ground states of all the subsystems of the 10-atom tree molecule are doubly degenerate. The ground state found by the VQE or, in our case, the direct diagonalization is therefore not unique. We use the spin state to distinguish the degenerate states. We mark the spin-up ground state with $\uparrow$ and the spin-down ground state with $\downarrow$. We only show a selection of all possible spin configurations due to redundancy resulting from the symmetry in the molecule. We observe two distinct resulting energy levels, indicating that the stating vector $\ket{G_{i}}$ can considerably influence which Hilbert space can be searched for the overall ground state and, therefore, the performance of the deep VQE.

The choice of the starting vector is not trivial, and without further knowledge of the system, all possible starting vectors have to be considered. Additional starting vectors, however, come at the cost of requiring additional qubits to represent $H^{\mathrm{eff}}$.

Next, we compare the performance of different basis methods with different starting vectors introduced in Sec. III using the treelike molecules. We compare two choices of starting vectors. The first choice uses single ground states of each subsystem, whereas the second uses additional eigenstates of each subsystem. In the former case, we chose the spin configuration $\downarrow\downarrow\downarrow\downarrow$ for the starting vectors for 10-atom tree molecule since we had degenerate eigenstates. In the case of the 13-atom tree molecule, only the central subsystem had a twofold degeneracy, and we also chose the spin $\downarrow$ ground state as the starting vector. For the interaction method, an $\epsilon=10^{-2}$ is chosen for the 10-atom tree molecule and an $\epsilon=10^{-3}$ for the 13-atom tree molecule. The number of considered starting vectors is indicated in the methods name in the bracket. The first digit indicates the number of starting vectors for the central subsystem, whereas the three following digits count the starting vectors for the branchlike subsystems.

Figures LABEL:fig:all_methods_error_10_atom_distance and LABEL:fig:all_methods_mean_error_10_atom_distance show the corresponding results for the 10-atom tree system. We also show the number of qubits needed for each method in Table 3. We note that we are unable to perform the calculation for the particle conserving method with additional starting vectors due to its sizeable computational requirement. The particle conserving methods with only one starting vector per subsystem performed exceptionally well in terms of accuracy. For all tested basis creation methods, including an additional eigenstate significantly increased the accuracy of the deep VQE. This is a consequence of the expanded Hilbert space that can be explored.

Unlike the other methods, we observe that the interaction methods with a fixed truncation have a decreasing accuracy with increasing stretching factors. More concretely, the energy obtained by interaction methods jumps to the combined subsystem’s energy at the stretching factor of 1.4. The interaction methods also display different behavior regarding the saved qubits, showing a significant dependence on the stretching factor as shown in Table 3. This dependence is due to our selection process of the interactions to generate a basis set. We set a fixed cut-off of the interaction strength $\lambda^{\mu}$ for the interactions involved in the basis creation step. Consequently, as the interactions between subsystems become weaker by stretching the molecule, more interactions are disregarded. Lowering the dimensionality of the basis sets allows us to save more qubits but comes at the cost of lowering the method’s accuracy.

Fig. LABEL:fig:all_methods_error_13_atom_distance and Fig. LABEL:fig:all_methods_mean_error_13_atom_distance show the results for the 13-atom tree molecule. Here, only the central subsystem is doubly degenerate, and the first excited state of the branch subsystems for the 13-atom tree molecule is a triplet. To see the effect of adding additional excited states as stating vectors on the performance of the deep VQE, we consider all three states as additional starting vectors. Due to the high computational cost for the other methods, only the single Pauli edges and the particle conserving edges methods could be prepared with additional starting vectors. For both methods, the additional starting vectors are able to increase the accuracy of the deep VQE, but this comes at the expense of requiring extra qubits (see Table 4 for the exact number of qubits required for each method). The particle conserving edges methods performed slightly better than others. This is consistent with the case of the 10-atom tree molecule, indicating an advantage of using the particle conserving approach in deep VQE for chemistry problems.

We show the results for the fixed qubit numbers approach in Figs. LABEL:fig:Results_flex_coef_10_atoms_compare_fix_qubits to LABEL:fig:Results_flex_coef_13_atoms_compare_fix_qubits_mean_error. We fix the number of qubits to 11, 14, and 17. For the 13-atom tree molecule, we could not produce an 11-qubit version as already including only the strongest interaction strength results in a basis set requiring more than 11 qubits. The qubit number is increased by three per step to take into account that there are three equivalent branchlike subsystems of the molecules.

In Figs. LABEL:fig:Results_flex_coef_10_atoms_compare_fix_qubits and LABEL:fig:Results_flex_coef_13_atoms_compare_fix_qubits, we see that if we define a fixed number of qubits, we can avoid the decreasing accuracy for increasing stretching factors. However, if the number of qubits is too restrictive, we cannot improve the result from the ”combined subsystems”. We see this behavior for the interaction(2222) method with 11 qubits for the 10-atom tree molecule in Fig. LABEL:fig:Results_flex_coef_10_atoms_compare_fix_qubits.

Additional to the toy model of a 10 and 13-atom hydrogen trees, we also apply the deep VQE algorithm to retinal to examine the performance on a natural molecule. We apply different basis forming strategies and compared them based on their accuracy and the number of qubits in Tab. 5. The calculation is performed with no stretching factors applied. The molecule got separated into two 10-qubit subsystems. With such a division, the eigenstates of the subsystems are double degenerate. In the case when only one starting vector for the basis forming step was considered, we choose the spin-down starting vector.

Deep VQE proves to be effective in treating such a molecule as retinal. Especially the addition of additional states as starting vectors seems to be a valid strategy. For both the 12 and the 14 qubit interaction treatment, the addition of a second starting vector resulted in a significantly better ground state energy approximation. This can be a consequence of the doubly degenerate ground state of the individual subsystems. The particle conserving edges method with each 10 starting vectors performed again exceptionally well and is able to approximate the ground state energy in the STO-3G basis within chemical accuracy while saving 8 qubits.

Overall we achieve similar accuracy for the 10 and the 13-atom tree molecules (see Figs. LABEL:fig:all_methods_error_10_atom_distance and LABEL:fig:all_methods_error_13_atom_distance). We expect this from the minor influence the newly added outer atoms have on the overall ground state. This behavior would also explain the remarkable success the edge methods provide, saving up to 15 qubits with comparable accuracy to the other methods Tab. 4. All methods are able to produce lower energies than the Hartree-Fock method. Multiple methods are able to approximate the ground state energies within an error of below 1% of the electron correlation energy of the molecule. The electron correlation energy is equivalent to the error of the Hartree-Fock method shown in Figs. LABEL:fig:all_methods_error_10_atom_distance and LABEL:fig:all_methods_error_13_atom_distance and represents the error due to the mean-field approximation of Hartree-Fock. This improvement helped us reach chemical accuracy for some stretching factors. The particle conserving edges method performed exceptionally well both for the treelike molecules as well for retinal. Especially the use of additional starting vectors proved to be beneficial for retinal. Using an edge method allows us to focus on changes affecting the orbitals involved in the most substantial interaction. We expect electrons occupying such orbitals to experience the most dramatic changes from the individual subsystem solutions when forming a bond with the other subsystems. To use the computational resources offered most effectively, it is advisable to use basis sets that focus on exploring the changes to the occupation of the most involved orbitals.

In contrast to the other methods, the interaction methods with a fixed truncation for the participating interactions significantly depended on the stretching factor. The number of qubits decreased with increasing distance. This reduction allows us to save more and more qubits but comes at the cost of decreasing accuracy. However, we consider this is not preferable if we have access to a quantum device with a fixed number of qubits. It seems unreasonable to not use all of them. The fixed-number qubit-interaction method we propose in this paper seems a more reasonable approach as it allows the use of all qubits. It also performs more stable for less interacting systems, not suffering from decreasing accuracy.