Subspace Variational Quantum Simulator

Kentaro Heya, ∗ Ken M. Nakanishi, † Kosuke Mitarai, 4, ‡ and Keisuke Fujii 5, § Research Center for Advanced Science and Technology (RCAST), The University of Tokyo, Meguro-ku, Tokyo 153-8904, Japan Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan. Graduate School of Engineering Science, Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka 560-8531, Japan. QunaSys Inc., High-tech Hongo Building 1F, 5-25-18 Hongo, Bunkyo, Tokyo 113-0033, Japan. JST, PRESTO, 4-1-8 Honcho, Kawaguchi, Saitama 332-0012, Japan (Dated: April 19, 2019)


I. INTRODUCTION
We are now in the era of noisy intermediate-scale quantum (NISQ) devices [3]. A NISQ device is a controllable quantum system, which still has a non-negligible amount of errors in its operation. However, it is believed that such a quantum device with a sufficient number of qubits and a low-enough error rate is not simulatable on classical computers [4][5][6]. In this sense, NISQ devices might have higher computational power than a classical computer. Especially, we believe that they are likely to be advantageous in preparing a quantum state of interest or simulating a time evolution of a quantum system.
For such an application of NISQ devices, the variational quantum eigensolver (VQE) [7][8][9][10] has recently attracted much attention. The VQE was proposed originally as an algorithm to find the ground state of a given Hamiltonian H. There, we use a parametrized quantum circuit U (θ) and optimize the parameter θ to minimize an expectation value of the Hamiltonian H = ψ(θ)|H|ψ(θ) , where |ψ(θ) is an ansatz state obtained by applying the circuit U (θ) to an input state |ϕ in . Successfully optimizing the appropriate parame-terized quantum circuit, we can get the ground state of the Hamiltonian as the output of the circuit.
The dynamics of quantum systems can exhibit interesting phenomena such as dynamical quantum phase transition [11] or discrete time crystals [12,13], which are yet to be fully understood. Furthermore, dynamical properties of quantum systems are practically important, and hence researchers have used methods like the time-dependent density functional theory [14]. As Feynman stated [15], a controllable quantum system must be advantageous in simulating such dynamics over a classical computer.
As for the simulation of the time evolution of a quantum system, Ref. [1] has proposed a variational quantum simulator (VQS), which simulates the time evolution of quantum states under a given Hamiltonian H. In the VQS, time-dependent parameter θ(t + δt) is utilized as the parameter of the quantum circuit. The algorithm determines the parameter θ(t + δt) by minimizing the difference between the actual time-evolved state e −iHδt |ψ(θ(t)) and |ψ(θ(t + δt)) . For this minimization, one measures the gradient of energy expectation value V = . With these, θ(t + δt) is calculated according to the differential equation M dθ(t) dt = V . This method is a pioneering work for performing a quantum simulation without Trotterization on NISQ devices.
In this paper, we propose another approach to simulate arXiv:1904.08566v1 [quant-ph] 18 Apr 2019 a time evolution of a quantum system, which we call the subspace variational quantum simulation (SVQS), and demonstrate proof-of-principle experiments on Rigetti Quantum Cloud Service. The SVQS simulates dynamics that evolve in a space spanned by the low-lying excited states and the ground state. This is based on our recent work, subspace-search variational quantum eigensolver (SSVQE) [2], which is an extension of the VQE to a subspace search problem. Our algorithm works as follows. First, we run SSVQE to find a unitary U (θ * ) that maps k orthogonal input states {|ϕ i } k i=1 , which we choose from computational basis, to the excited states up to the k-th, {|E i } k i=1 of H, each having the corresponding eigenenergy E i . Note that the inverse of the unitary, U † (θ * ), maps each |E i to computational basis |ϕ i , where we can easily apply phase factors e iEit . Therefore, we can perform the simulation of e −iHt by applying U † (θ * ) to an input state |ψ in ∈ span{|E i } k i=1 , then applying the phase, and finally applying U (θ * ) to return the state to its initial basis.
Our protocol does not require the gradient or the metric tensor for simulation. The only thing we have to do is to find the unitary U (θ * ) that transforms the low energy subspace to a subspace spanned by computational basis. The simplicity of our protocol enables us to do the simulation with lower overhead than the VQS. Moreover, we do not need to discretize the simulation by a small timestep δt, that is, we can perform the simulation of an arbitrary time duration without increasing the depth of the circuit. Although the simulation is restricted into space spanned by the low-lying states in the proposed method, they are often the states of the interest. We successfully demonstrate the proof-of-principle experiments of our method on Rigetti's superconducting quantum computer. In the experiments, we employ error-mitigation techniques [1,16,17] to observe quantum dynamics accurately even in the presence of noise.
This work presents an efficient algorithm to simulate the quantum dynamics on a NISQ device. While the dimension of the simulatable subspace is limited, this algorithm shares the idea with the long-term quantum algorithms based on quantum phase estimation such as Harrow-Hassidim-Lloyd algorithm [18] for solving linear problem and the quantum recommendation systems [19]. More precisely, both approaches transform the input state to another basis on which we can manipulate their eigenvalues easily and then undo the basis change. From this respect, the SSVQE and the SVQS can be regarded as a NISQ version of quantum phase estimation and its application, respectively. This correspondence might provide a deeper insight into the relationship between the near-term quantum algorithms and long-term quantum algorithms.

II. METHODS
In this section, we first introduce the subspace-search variational quantum eigensolver (SSVQE) [2], which is the key ingredient for our proposed method. We then describe our proposal, the subspace variational quantum simulator (SVQS). In the following subsections, H denotes a given n-qubit target Hamiltonian, which can be written as follows, where E j and |E j are j-th eigenenergy and eigenstate of H, respectively. In application for quantum chemistry, Jordan-Wigner or Bravyi-Kitaev transformation can be utilized to map a fermionic Hamiltonian to an n-qubit Hamiltonian.

A. Subspace-search variational quantum eigensolver
The SSVQE is an algorithm for finding the k-th or up to the k-th excited states of Hamiltonian H. To find excited states up to the k-th, the SSVQE takes k orthogonal states as inputs of a parametrized quantum circuit and minimizes the weighted sum of the expected values of the energies of the output states. The output states are automatically orthogonal by the conservation of orthogonality under the unitary transformation, and therefore we are able to find all excited states up to the k-th by only one optimization procedure. A more detailed explanation is found in Ref. [2]. The procedure of the SSVQE is summarized as follows: 1. Construct a parameterized quantum circuit U (θ) and prepare k initial states {|ϕ j } k j=0 (k ≤ n), which are orthogonal with each other ( ϕ i |ϕ j = δ i,j ).

Minimize
Successfully optimizing θ of the appropriate parameterized quantum circuit U (θ) by the above procedure, each output state |ψ j (θ) ≡ U (θ) |ϕ j (j = 0, 1, · · · , k) converges to the following state where δ j is an unknown global phase factor, and θ * denotes the parameters theta that minimizes L ω (θ). Therefore, the obtained circuit U (θ * ) corresponds to a map between the computational subspace S com formed by the orthogonal initial states {|ϕ j } k j=1 and the eigensubspace S eig formed by the excited states {|E j } k j=1 . While the above unknown global phase factors are though to cause a problem to prepare a desired superposition of the energy eigenstate, this is not the case as we will see later.

B. Subspace variational quantum simulator
In this subsection, we describe our proposal, which we name the subspace variational quantum simulator (SVQS). The key idea of the SVQS is mapping low-energy eigensubspace S eig of the target Hamiltonian H to subspace S com spanned by computational basis which can be controlled easily. The procedure of the SVQS is summarized as follows: 1. Construct a parameterized quantum circuit U (θ) and prepare l input states {|ϕ j | |ϕ j = X j |0 ⊗n } l j=1 (l ≤ n) such that they are orthogonal with each other ( ϕ j |ϕ j = δ j,j ).

Minimize
where the weight vector ω is chosen such that ω i > ω j when i > j. Hereafter, θ * = arg min L ω (θ).
3. Prepare an initial state |ψ in in the eigensubspace S eig .
4. Encoding the input state |ψ in on the eigensubspace S eig into the computational subspace S com by applying the Hermitian conjugate of the obtained circuit U † (θ * ).
5. Apply single-qubit Z-rotation on each qubit, namely, are the eigenenergies of the eigenstates {|E j } l j=1 obtained by SSVQE.
6. Decoding the state T (t)U † (θ * ) |ψ in on the computational subspace S com into the eigensubspace S eig by applying U (θ * ). See also Fig. 1 which is the quantum circuit corresponding to the step 4-6 of the procedure. Let us explain how the above procedure works to simulate the time evolution of the quantum system. The obtained circuit U (θ * ) corresponds to a map between S eig and S com . Therefore, the circuit U (θ * ) can be denoted as follows, where {δ j } l j=1 are unknown phase offsets and U ⊥ is an unknown map between subspaces complementary to S eig and S com . The tensor product of single-qubit Z-rotations can be denoted as follows, where U ⊥ (t) is again a map between the complementary subspaces. Under the conjugation by Eq. (3), Eq. (4) transforms as follows: Equation (5) corresponds to the time evolution operator on the k-dimentional eigensubspace with unknown timedependent unitary operation on the complementary subspace. Note that the number of excited states that we can address here is restricted at most n because of the construction of the initial orthogonal states by one-hot states |ϕ j = X j |0 ⊗n . Note that we here assume that the experimentalist can prepare the initial state |ψ in of his/her interest appropriately without any knowledge about the global phase factors. This assumption would be unlikely, and hence we will relax it in the next subsection.

C. Variational initial state prepraration for SVQS
In Sec. II B, we assume that the experimentalist can find and prepare the desired initial quantum states in the low-enegy eigensubspace S eig which we call a simulatable eigensubspace. However, since there are unknown global phase factors on the obtained energy eigenstates, and hence it is not straightforward for the experimentalist to prepare the desired initial quantum state. In this subsection, we relax this issue by variationally preparing the initial state of interest without any knowledge on the global phase factors.
The dimension of the simulatable eigensubspace S eig in SVQS is restricted by the number of qubits used. Therefore, in general cases, there is no guarantee that the desired initial state is included in the simulatable eigensubspace. Thus, we redefine the initial state preparation as a minimization problem of some cost function F (|ψ in ), which reflects the experimentalist's requests. For example, the cost function can be defined as follows, where O i is some observable, and O opt i is its expectation value that the experimentalist wants. In such cases, the procedure of searching the desired initial quantum states in the simulatable eigensubspace can be further regarded as VQE in the simulatable eigensubspace.
The initial states |ψ in for SVQS are generated by parameterized quantum circuits U ST (φ). If the parameterized quantum gates closed in the computational subspace S com , such as R Z (θ) and exp (−iθ (|01 10| + h.c.)), are feasible, the parameterized quantum circuit U ST (φ) can be efficiently composed of them.
After the successful optimization of φ to minimize F (|ψ in ), the parameterized initial quantum state converges to the |ψ in (φ * ) , where φ * denotes the optimal value of φ. Then, the time-developed state |ψ in (φ * , t) can be written as follows: D. Extension of simulatable eigensubspace using ancilla qubits In the methods described in Sec. II B and II C, the dimension of the simulatable eigensubspace S eig is limited by n because of the dimension of the computational subspace. In this section, we propose an additional method to extend the dimension from n to n+a by using a ancilla qubits. The Hamiltonian of the (n + a)-qubit system is defined in terms of the H problem as follows: where I and Z denote the identity and Pauli-Z operator on a single-qubit respectively, and B is a constant set much larger than the (n + a)-th excited eigenenergy. The first term of the target Hamiltonian corresponds to the problem Hamiltonian on the n-qubit system and the second term of the Hamiltonian corresponds to the static magnetic field B applied on Z-axis for the ancilla qubits. If the magnetic field B is larger than the eigenenergy of the (n + a)-th excited eigenstate of the problem Hamiltonian, the i-th excited eigenstate of the target Hamiltonian becomes where |ψ i problem denote the i-th excited eigenstate of the problem Hamiltonian. Notice that such a method does not disturb the low-lying eigenstates of the problem Hamiltonian. The a-ancilla qubits make it possible to use the (n + a)-orthogonal input states X j |0 This allows us to search up to (n + a)-th excited states of the n-qubit Hamiltonian.

III. EXPERIMENT
Here, we demonstrate our algorithm with 2-qubit experiments on Rigetti Quantum Cloud Service. The anzats circuit U D (θ), which we employed in following experiments, consists of X-rotation R x (π/2), parametrized Z-rotation R z (θ), and control-Z gate CZ as shown in Fig. 3. Here, we adopt the virtual-Z gate as the parameterized Z-rotation R z (θ), whose gate fidelity is known as much higher than other quantum gates, R x (π/2) and CZ [20].
The initial values of the parameters are randomly sampled from a uniform distribution [0, 2π). For the optimization of the parameters in the weighted-SSVQE, we used an optimization method recently reported in Ref. [21]. For the optimization of the parameters of the variational state preparation in the simulatable eigensubspace, we used the BFGS method [22,23] implemented in the SciPy library [24].

A. Single-qubit Rabi oscillation
First, we demonstrate our idea with a single-qubit Rabi oscillation as the simplest quantum dynamic in a twolevel system. We consider the Hamiltonian given as follows: where Ω = 0.612 is a sample from a uniform distribution on [0, 1). In this case, the system that we want to simulate is two dimensions, but the dimension of simulatable subspace is one because the system consists of one qubit. Therefore, we employ one ancilla qubit to expand the simulatable subspace from one to two dimensions, which allows us to simulate the Rabi oscillation on a two-level system. An additional Hamiltonian for the ancilla qubit is given as follows: where the coefficient h = 5 is selected to be sufficiently larger than the maximum eigenenergy of the system Hamiltonian. The total Hamiltonian becomes Here, we adopt the ansatz circuit U D (θ) with D = 1 as shown in Fig. 3. The optimization process of SSVQE is shown in the appendix. After the convergence of the numerical optimization of the circuit with the weighted-SSVQE on a classical computer, we construct the quantum circuit for the time evolution    Here, we perform the simulation for several initial states of the data qubit given by R Y (i/8π) |0 with an integer i (−4 ≤ i ≤ 4) and the ancilla qubit as |0 , which are each shown by different colors in Fig. 5.
To obtain an ideal result even in the presence of noise, we use the following error mitigation method [25][26][27]. Using the result from the error-increased quantum circuits, we apply a linear extrapolation to restore the noiseless data. In order to construct such an error-increased circuit U (E) (θ), we replace quantum circuit elements R X (π/2) and CZ in U (θ) with redundant quantum circuit elements R (E) x (π/2) and CZ (E) for E ∈ [1, 3, 5, · · · ]. They can be denoted as follows, In the ansatz circuit, the fidelity is decreased mainly by R x (π/2) and CZ gates. Therefore, the error in U (E) (θ) is roughly E-times larger than that of U (1) (θ). From Figs. 5 and 6, one can see that the states of the data qubit rotate in X-axis and the states of the ancilla qubit stays in |0 state, which is the behavior expected from the Hamiltonian.

B. Hydrogen molecule
Next, we demonstrate SVQS with variational initial state preparation for hydrogen molecule using Rigetti QCS. According to Ref. [28], the Hamiltonian of the hydrogen molecule in STO-3G basis can be converted into 2-qubit Hamiltonian as follows, where the coefficients c i are calculated by openfermion [29] and Psi4 [30]. Here, the occupancy of the coupled and anti-coupled orbital can be calculated from the expectation values of the Pauli operators as follows, Figure 7 shows the eigenvalue dependence of the Hamiltonian of hydrogen molecules with respect to the interatomic distance. We simulate the dynamics of the Hamiltonians at the interatomic distance 0.2Å and 1.0Å. The mapping from the computational subspace to the eigensubspace of the hydrogen molecule is obtained from the numerical calculation of weighted-SSVQE as shown in the appendix. In the weighted-SSVQE, we adopt the ansatz circuit U D (θ) with D = 2 as shown in Fig. 3. Next, we prepare the initial quantum state for SVQS. Here, we adopt the cost function for the state preparation as follows, which means that the occupancy of the coupled and anticoupled orbital need to be close to 0.5 to reduce the cost value. As shown in Figs This is because, when the interatomic distance is around 0.61Å, the first excited state changes, and thus the quantum state included in the eigensubspace formed by the ground and the first excited state changes.
After the convergence of the numerical optimization of the circuit with weighted-SSVQE on a classical computer, we construct the quantum circuit for the simulation of the time evolution of the system as shown in Fig. 10, where θ * and φ * are converged parameters. Figures 11 and 12 shows the actual experimental results of the quantum simulation with 10 4 shots.  From the Figs. 11 and 12, one can see that when the hydrogen atoms are close to each other, the orbital occupancy rate does not change. This is because, when the interatomic distance is closer, the Hartree-Fock orbitals nicely approximates the eigenstate of the Hamiltonian. On the other hand, When the hydrogen atoms are separated from each other, the orbital occupancy rate changes, because of the poorness of the approximation performed by the Hartree-Fock orbitals.

IV. CONCLUSION
We proposed an efficient algorithm to simulate the quantum dynamics on a NISQ device by replacing the QPE with the SSVQE. The circuit depth of the proposed method is at most twice as large as that of VQE. The simplicity of the proposed method enables us to demonstrate it on the existing experimental quantum devices. The proposed method is the first attempt to apply the SSVQE as a near-term version of the QPE, and provides a deeper insight into the relationship between the near-term quantum algorithms and long-term quantum algorithms. For example, it might be able to do more complex processing, such as the HHL algorithm and the simulation of the relaxation process in the low-energy eigensubspace. Also, it would be interesting to do quantum signal processing with simulatable eigensubspace.

V. ACKNOWLEDGEMENTS
This work was supported by QunaSys and Rigetti Inc. The numerical simulations were performed with Quantum Toolbox in Python (QuTiP) [31], Open Source Scientific Tools for Python (SciPy) [24], and QuPy [32].
Appendix A: Numerical SSVQE in SVQS procedures 1. Numerical SSVQE for single-qubit Rabi model From the Fig. 13, one can see the formation of the energy subband caused by the magnetic field imposed on the ancilla qubit, and weighted-SSVQE can actually find desired excited states in the lowest energy sub-band all at once.