Subspace-search variational quantum eigensolver for excited states

The variational quantum eigensolver (VQE), a variational algorithm to obtain an approximated ground state of a given Hamiltonian, is an appealing application of near-term quantum computers. The original work [A. Peruzzo et al.; \textit{Nat. Commun.}; \textbf{5}, 4213 (2014)] focused only on finding a ground state, whereas the excited states can also induce interesting phenomena in molecules and materials. Calculating excited states is, in general, a more difficult task than finding ground states for classical computers. To extend the framework to excited states, we here propose an algorithm, the subspace-search variational quantum eigensolver (SSVQE). This algorithm searches a low energy subspace by supplying orthogonal input states to the variational ansatz and relies on the unitarity of transformations to ensure the orthogonality of output states. The $k$-th excited state is obtained as the highest energy state in the low energy subspace. The proposed algorithm consists only of two parameter optimization procedures and does not employ any ancilla qubits. The disuse of the ancilla qubits is a great improvement from the existing proposals for excited states, which have utilized the swap test, making our proposal a truly near-term quantum algorithm. We further generalize the SSVQE to obtain all excited states up to the $k$-th by only a single optimization procedure. From numerical simulations, we verify the proposed algorithms. This work greatly extends the applicable domain of the VQE to excited states and their related properties like a transition amplitude without sacrificing any feasibility of it.


I. INTRODUCTION
Supported by the world-wide active research for the development of quantum devices, quantum computers equipped with almost a hundred qubits are now within reach. Those near-term quantum computers are often called noisy intermediate-scale quantum (NISQ) devices [1], reflecting the fact that those quantum computers are not fault-tolerant, that is, they do not have the guaranteed accuracy of the computational result. However, such a NISQ device is believed not to be simulatable on classical computers if the gate fidelity is sufficiently high [2][3][4]. This fact encourages us to look for practical applications of them.
The variational quantum eigensolver (VQE) [5][6][7] is an attracting application of near-term quantum computers in the hope that the controllable quantum devices can simulate another quantum system more efficiently than classical devices. The VQE is an algorithm for finding an approximate ground state of a given Hamiltonian H. For this purpose, the VQE utilizes a parameterized quantum circuit U (θ), which is also called an ansatz circuit, to generate an ansatz state |ψ(θ) . The expectation value of the target Hamiltonian H(θ) = ψ(θ)| H |ψ(θ) is minimized by iterative optimization of the parameters θ. The circuit with the resultant optimal parameters θ * which minimizes H outputs the approximate ground state.
Not only the ground state, which the original VQE aims to find, but also excited states of molecules are responsible for many chemical reactions and physical processes. For example, the transition between a ground state and excited states is the origin of luminescence [8]. Intermediate states of a chemical reaction are, in general, not a ground state of a system, and therefore properties of such excited states are important for the analysis of them [9].
In spite of the importance of the excited states, classical computation suffers from the increasing computational cost and gives relatively poor results for them [9][10][11]. This motivates us to utilize quantum computers for the task of finding excited states and analyzing their property. A long-term quantum algorithm for a chemical reaction has been investigated in Ref. [12]. However, algorithms which we can run on NISQ devices are yet to appear.
In order to find the excited states using NISQ devices, we propose a method, which utilizes the conservation of orthogonality under the unitary transformation. We name the method as the subspace-search VQE (SSVQE). The SSVQE takes two or more orthogonal states as inputs to a parametrized quantum circuit, and minimizes the expectation value of the energy in the space spanned by those states. This method automatically imposes the orthogonality condition on the output states, and therefore allows us to remove the swap test [13], which has been employed in the previous works [14,15] to ensure the orthogonality. In principle, the proposed algorithm can find the k-th excited state by running optimization of the circuit parameters only twice. We also propose a generalized version of the SSVQE, which finds all excited states up to the k-th only one optimization procedure. As a possible application of the SSVQE, a method to measure a transition amplitude between two eigenstates is described. It can evaluate material properties such as permittivity and rate of spontaneous emission. We perform numerical simulations and show validity of proposed algorithms for fully connected random transverse Ising models and Helium hydride. This work greatly extends the practicability of the VQE by enabling it to find the excited states efficiently, and thereby, pushes the VQE as a candidate for a possible application of NISQ devices further.
The rest of the paper is organized as follows. In Sec. II we first propose the algorithm of the SSVQE and the extended version of it. Then, in Sec. III we briefly review the existing works addressing the same objective of finding the excited states in the framework of the VQE. An algorithm to obtain the transition amplitude is described in Sec. IV. Finally, we present the simple, proofof-principle numerical simulations in Sec. V.

II. METHODS
The VQE is a quantum-classical hybrid algorithm to find a ground state of a given Hamiltonian H using NISQ devices. For this purpose, the VQE utilizes a parameterized quantum circuit U (θ), also called an ansatz circuit, to generate an ansatz state |ψ(θ) . The expectation value of the target Hamiltonian H(θ) = ψ(θ)| H |ψ(θ) is minimized by iterative optimization of the parameters θ.
Our objective here is to find excited states of the Hamiltonian H. Since the eigenstates of the Hamiltonian H are mutually orthogonal, a straightforward construction of the algorithm to find k-th excited state is to minimize H(θ) imposing an orthogonality condition between the ansatz state |ψ(θ) and all of the ground/excited states up to (k − 1)-th. By inductively repeating this, we can find the excited states of interest. The swap test [13], which can measure the inner product between the ground state and the ansatz state, has been employed to ensure the orthogonality in the previous works [14,15]. In contrast, the SSVQE and the weighted SSVQE we propose here utilize the conservation of orthogonality under the unitary transformation in an effort to satisfy the orthogonality condition. These methods automatically impose the orthogonality on the output states, and therefore remove the swap test.

A. Subspace-search variational quantum eigensolver
The key idea is to ensure the orthogonality at the input of the quantum circuit, not at the output. Below we describe the algorithm to find the k-th excited state that works on an n-qubit quantum computer. We define the ground state as the 0-th excited state. The algorithm, which we refer to as the subspace-search VQE (SSVQE), runs as follows. Algorithm: 1. Construct an ansatz circuit U (θ) and choose in- We denote the optimal θ by θ * .

Construct another parametrized quantum circuit
V (φ) that only acts on the space spanned by 4. Choose an arbitrary index s ∈ {0, · · · , k}, and maximize We note that, in practice, the input states {|ϕ j } k j=0 will be chosen from a set of states which are easily preparable, such as the computational basis.
Let the set of eigenstates of Then, the circuit optimized by the step 2 of the above algorithm is a unitary that best approximates the mapping from the space spanned by {|ϕ j } k j=0 to one spanned by {|E j } k j=0 . Therefore, in step 2, we can find the subspace which includes |E k as the highest energy state, using a carefully constructed ansatz U (θ). The unitary V (φ) is responsible for searching in that subspace. By maximizing L 2 (φ), we find the k-th excited state |E k .
In the case of k ≥ 2 n−1 , it is faster to choose 2 n − k of orthogonal input states |ϕ j and maximize L 1 (θ) instead of minimizing it in the step 2, then minimize L 2 (θ) instead of maximizing it in the final step.

B. Weighted SSVQE for finding the k-th excited state
Here we extend the algorithm described in the previous section to find the k-th excited state of a given Hamiltonian which requires only a single optimization procedure. It runs as follows. Algorithm: 1. Construct an ansatz circuit U (θ) and choose input states {|ϕ j } k j=0 which are orthogonal with each other ( ϕ i |ϕ j = δ ij ).
When the cost L w reaches its global optimum, the circuit U (θ) becomes a unitary which maps |ϕ k to the k-th excited state |E k of the Hamiltonian and others to the subspace spanned by {|E j } j=k−1 j=0 . Therefore, by minimizing the cost L w , we can find the k-th excited state by a single optimization process. Note that the overall time required for the optimization might increase, due to the more complicated landscape of the cost function.
C. Weighted SSVQE for finding up to the k-th excited states We further generalize the above argument and propose an algorithm for finding all excited states of a given Hamiltonian up to the k-th with only one optimization procedure. Algorithm: 1. Construct an ansatz circuit U (θ) and choose input states The weight vector introduced here has the effect of choosing which |ϕ j is converted to which excited state. It is easy to see the circuit U (θ) when the cost L w reaches its global optimum becomes a unitary which maps |ϕ j to the j-th excited state |E j of the Hamiltonian for each j ∈ {0, 1, · · · k}. In this case, too, note that the overall time required for the optimization might increase due to the same reason as the previous section.

III. RELATED WORKS
In this section, we first overview previous works, and then point out the advantages of our methods over them.
Ref. [16] has proposed a method which hybridizes the quantum phase estimation algorithm and the VQE. Although it is experimentally demonstrated [16], the method is unlikely to be implemented on a NISQ device, due to the need for the controlled time evolution.
In Ref. [17], a method called quantum subspace expansion has been proposed. The algorithm first finds the ground state |E 0 by the usual VQE protocol, and then measures the matrix elements of the Hamiltonian with respect to the space spanned by {O α |E 0 }, where {O α } is a set of excitation operators. The diagonalization of the matrix, which is done classically, can determine the approximate eigenvalue spectra. They have used a set of one electron excitation operators {a † v a o }, where a † i and a j are the fermion creation and annihilation operators respectively and i, j running all possible indices, for {O α }.
The constrained VQE proposed in Ref. [18] can also be used for finding a certain set of excited states. They proposed a way to introduce constraints, such as the number of electrons or the overall spin of the system, on the VQE. The introduction of the constraints is done by adding the penalty term to the cost function. Their method finds the lowest energy state under the constraints. Since the difference in the constraints, such as the difference in the number of electrons or the overall spins, generally changes the energy of the system, it can be utilized to find a certain set of excited states.
Ref. [14] has recently proposed an inductive method which adds a penalty term to ensure the orthogonality of the ansatz state with respect to the low-lying state. To be more concrete, to find the k-th excited state, they use H where θ * i is the optimal parameters for the i-th excited state and β i is a hyperparameter that determines the strength of the penalty, as the target cost function to be minimized by tuning θ k . To estimate the overlap, their method uses the swap test, which requires us to double the number of qubits with additional gates. Their method works well when the hyperparameter β i is set properly as shown in Ref. [14]. Ref. [19] has enabled the optimization by the imaginary time evolution of the parameters in this approach.
The advantages of our methods, when compared to the methods above, are as follows.
1. The energy spectrum found by the SSVQE or the weighted SSVQE is exact when U (θ) and V (φ) have the ability to represent the exact unitary which maps the k input states to the k eigenstates of the Hamiltonian.
2. The swap test is not employed and thus easily implementable on the NISQ devices.
3. In the SSVQE, there are no hyperparameters.
4. In the weighted SSVQE, the results are unique regardless of the value of the hyperparameters if they meet the conditions. 5. Optimization runs only twice for the SSVQE and only once for the weighted SSVQE.

IV. CALCULATION OF TRANSITION MATRIX ELEMENTS
It is possible to measure a transition amplitude of an operator A, E i |A|E j , using the result of the SSVQE. Note that E i |A|E j = ϕ i |U † (θ * )AU (θ * )|ϕ j where U (θ * ) is the optimized unitary. We can measure this by expanding it as: where Recall that in practice the input states are chosen from simple states such as the computational basis, and therefore we assume that the superpositions like |+ x ij can easily be prepared. Each term of the above equation are measured separately on the NISQ device and then are summed up on a classical computer.

V. NUMERICAL SIMULATION
Here we numerically simulate our algorithms with 4qubit Hamiltonians. Figure 1 shows the variational ansatz used in the simulations. We chose the input states as {|ϕ j } = {|0000 , |0001 , |0010 , |0011 }. The depth D 1 is set to D 1 = 2 for all of them. D 2 is set to D 2 = 6 for the SSVQE and the weighted SSVQE for finding the k-th excited states, and D 2 = 8 for the weighted SSVQE for finding the all excited states up to the k-th. The initial values of the parameters were randomly sampled from a uniform distribution [0, 2π). For each simulation, the optimization was run for 10 times starting from different initial values. The results shown in the following sections are the ones which achieved the lowest value of the cost function among those 10 results. We used the BFGS method [20] implemented in the SciPy library [21] for the optimization of the parameters. A. Transverse Ising model First, we demonstrate our idea with a Hamiltonian of the fully connected transverse Ising model: with N = 4. The coefficients a i and J ij are sampled randomly from a uniform distribution on [0, 1). In this subsection, we use one Hamiltonian with the same coefficients as an example. All experiments were conducted on the case of k = 3.

SSVQE
The SSVQE can find the k-th excited state with only two optimization procedures. Figure 2 shows the first optimization process of θ to minimize L 1 (θ). In Fig. 2, the fidelity is defined by the overlap between the space spanned by {|E j } 3 j=0 and the output of the quantum circuit {U (θ) |ϕ j } 3 j=0 , namely, We see that, as the cost function gets close to its global minimum, the fidelity approaches unity as expected. Figure 3 shows the second process of optimizing φ to minimize L 2 (φ). Here the fidelity is defined by One can see the subspace-search approach works well from  2. Weighted SSVQE for finding the k-th excited state The method described in Sec. II B can find the k-th excited state by only one optimization sequence. Here, we chose w = 0.5 as the weight. Figure 4 shows the optimization process of θ to minimize L w (θ). Here the fidelity is defined by | E 3 |U (θ)|ϕ 3 | 2 . In this case, too, the algorithm succeeds in finding the third excited state of the Hamiltonian. However, the number of iterations to the convergence is larger than the number of overall iterations of the simple SSVQE. It might be attributed to the more complicated landscape existing in the cost function.

Weighted SSVQE
The weighted SSVQE described in Sec. II C can find 0, 1, · · · , k-th excited states all at once. Here, we chose w = (4, 3, 2, 1) as the weight vector. Figure 5 shows the optimization process of θ to minimize L w (θ). From Fig. 5, one can see that this approach can actually find the desired excited states all at once. The number of iterations to the convergence is almost equivalent to the one presented in the previous section.

B. Helium hydride
Next, we apply our idea for the molecular Hamiltonians of HeH with a fixed distance between two atoms. Our ansatz (Fig. 1) does not consider the conservation of the number of electrons, and therefore, the calculated excited states can have the different number of them. The molecular Hamiltonians are calculated with Open-Fermion and OpenFermion-Psi4 [22,23]. We used the STO-3G minimal basis set, and therefore, obtained the 4-qubit Hamiltonian. We calculated the Hamiltonians at the 24 different bond lengths and performed the VQE at each point. In the weighted SSVQE simulation, we used the same weights as in the previous section.
The result of SSVQE is shown in Fig. 6 and one from using weighted SSVQE for finding the k-th excited state is shown in Fig. 7. Both of the results agree nicely with the exact values of the third excited state at each bond length. Next, we used the weighted SSVQE described in Sec. II C to find all excited states up to the third. The result is shown in Fig. 8. One can see that the energy eigenvalues are well approximated by the optimized output of the weighted SSVQE. In this work, we proposed efficient algorithms for finding excited states of a given Hamiltonian, extending the framework of the VQE. The proposed method assures the orthogonality of the states at the input of the ansatz circuit. Minimizing a carefully designed cost function by optimizing the parameters of the quantum circuit, we can map each of the orthogonal states onto one of the energy eigenstates. Our algorithms require us to run, in principle, optimization only once or twice, and find one or more arbitrary excited states. We believe that this work greatly extends the practicability of the VQE for finding the excited states.