Variational quantum simulation of general processes
Suguru Endo, Jinzhao Sun, Ying Li, Simon Benjamin, Xiao Yuan
VVariational quantum simulation of general processes
Suguru Endo, ∗ Jinzhao Sun, Ying Li, Simon C. Benjamin, and Xiao Yuan † Department of Materials, University of Oxford, Parks Road, Oxford OX1 3PH, United Kingdom Clarendon Laboratory, University of Oxford, Parks Road, Oxford OX1 3PU, United Kingdom Graduate School of China Academy of Engineering Physics, Beijing 100193, China
Variational quantum algorithms have been proposed to solve static and dynamic problems ofclosed many-body quantum systems. Here we investigate variational quantum simulation of threegeneral types of tasks—generalised time evolution with a non-Hermitian Hamiltonian, linear alge-bra problems, and open quantum system dynamics. The algorithm for generalised time evolutionprovides a unified framework for variational quantum simulation. In particular, we show its appli-cation in solving linear systems of equations and matrix-vector multiplications by converting thesealgebraic problems into generalised time evolution. Meanwhile, assuming a tensor product struc-ture of the matrices, we also propose another variational approach for these two tasks by combiningvariational real and imaginary time evolution. Finally, we introduce variational quantum simulationfor open system dynamics. We variationally implement the stochastic Schr¨odinger equation, whichconsists of dissipative evolution and stochastic jump processes. We numerically test the algorithmwith a six-qubit 2D transverse field Ising model under dissipation.
Introduction.—
The variational method is a powerfulclassical tool for simulating many-body quantum sys-tems [1–5]. The core idea is based on the intuition thatphysical states with low energy belong to a small mani-fold of the whole Hilbert space. As quantum circuits canefficiently prepare states that may not be efficiently rep-resented classically, the variational method has been re-cently generalised to the quantum regime with trial statesefficiently prepared by a quantum circuit and informationextracted from a coherent measurement of the state [6–26]. The trial state in variational quantum algorithmscan be prepared with shallow quantum circuits [27–30],which is robust to a certain amount of device noise andis compatible with near-term Noisy Intermediate ScaleQuantum (NISQ) hardware [31]. Variational quantumalgorithms can be utilised for efficiently finding energyspectra [7–12, 19, 21–25, 32] and simulating real timeSchr¨odinger evolution [13, 33] of closed systems. Al-though quantum circuits are unitary operations, the vari-ational algorithm is not limited to energy minimisationand unitary processes and it can be used to simulate dissi-pative imaginary time evolution that cannot be straight-forwardly mapped to unitary gates [20, 34].In this work, we study the capability of variationalquantum algorithms and show that they are not lim-ited to these applications. First, we introduce a varia-tional quantum algorithm for simulating the generalisedtime evolution defined in Eq. (1) below. Our algorithmcan be regarded as a unified framework, which incorpo-rates the special cases of real and imaginary time evo-lutions [13, 20, 34], non-Hermitian quantum mechan-ics [35–37] that describes nonequilibrium processes [38],parity-time symmetric Hamiltonians [39–41], open quan-tum systems [42], general first order differential equa-tions, etc.Next we apply the variational method for solving lin-ear algebra problems, such as linear systems of equa- tions and matrix-vector multiplications, important tasksin machine learning and optimisation [43, 44]. Many al-gorithms have been developed for linear systems of equa-tions with universal quantum computers [45–52], whichhave profound applications in quantum machine learn-ing [53–57]. However they generally require deep circuitsthat rely on fault tolerant quantum computers. In thiswork, we introduce two types of variational quantum al-gorithms for the two linear algebra problems. For thefirst type, we consider general sparse matrices and showhow solutions of the problems can be converted into gen-eralised time evolutions, which can be variationally sim-ulated. For the second type, we consider special matricesthat are products of small matrices acting on a constantnumber of qubits, and use variational real and imaginarytime evolution to find solutions.Finally, we combine the developed variational al-gorithms to simulate the evolution of open quantumsystems [58–60]. Simulating the evolution of generalopen quantum systems is of great importance forunderstanding any quantum system that interacts withan environment. Existing quantum algorithms [61–66]for simulating open quantum systems generally requiredeep quantum circuits. In this work, we consider thedescription of open system dynamics via the stochasticSchr¨odinger equation, whose evolution can be regardedas an average of wave functions that undergo a continu-ous measurement induced from the environment [60, 67].The evolution of each wave function is composed of twoprocesses that can be both simulated with variationalalgorithms: It may continuously evolve under thegeneralised time evolution with the system Hamiltonianand the damping effect due to continuous measurement;alternatively the state discontinuously jumps accordingto the measurement results. The continuous process canbe described by the generalised time evolution, and thejump process is a matrix-vector multiplication process. a r X i v : . [ qu a n t - ph ] J un Therefore our algorithm is compatible with shallowcircuits and NISQ hardware.
Generalised time evolution.—
We first consider varia-tional quantum simulation of generalised time evolution, B ( t ) ddt | v ( t ) (cid:105) = | dv ( t ) (cid:105) . (1)Here | dv ( t ) (cid:105) = (cid:80) j A j ( t ) | v (cid:48) j ( t ) (cid:105) , A j ( t ) and B ( t ) are gen-eral time dependent sparse (non-Hermitian) operators, | v ( t ) (cid:105) is the system state, and each of | v (cid:48) j ( t ) (cid:105) can beeither | v ( t ) (cid:105) or any known state that can be efficientlyprepared with a quantum circuit. The states | v ( t ) (cid:105) and | v (cid:48) j ( t ) (cid:105) can be (un)normalised states | v ( t ) (cid:105) = α ( t ) | ψ ( t ) (cid:105) , | v (cid:48) j ( t ) (cid:105) = α (cid:48) j ( t ) | ψ (cid:48) j ( t ) (cid:105) with normalisation factors α ( t )and α (cid:48) j ( t ), respectively. In practice, we assume that A j ( t ) ( B ( t )) can be decomposed as a linear combina-tion of Pauli operators A j ( t ) = (cid:80) i λ ji ( t ) σ i with complexcoefficients λ i and a polynomial (with respect to the sys-tem size) number of tensor products of Pauli matrices σ i = (cid:78) i k σ i k with i k denoting the i k th qubit.In variational quantum simulation, instead of directlysimulating the dynamics, we assume that the state can berepresented by parameterised quantum states | v ( (cid:126)θ ( t )) (cid:105) = α ( (cid:126)θ ( t )) | ϕ ( (cid:126)θ ( t )) (cid:105) with (cid:126)θ := ( (cid:126)θ , (cid:126)θ ). Then we projectthe original evolution to the evolution of the parametersvia McLachlan’s principle [68],min (cid:13)(cid:13)(cid:13)(cid:13) B ( t ) ddt | v ( (cid:126)θ ( t )) (cid:105) − (cid:88) j A j ( t ) | v (cid:48) j ( t ) (cid:105) (cid:13)(cid:13)(cid:13)(cid:13) , (2)where (cid:107) | ψ (cid:105) (cid:107) = (cid:112) (cid:104) ψ | ψ (cid:105) and the minimisation is over thederivative of the parameters. By minimising the distancebetween the true evolution and the evolution of the pa-rameterised state, we find the equation of parameters as (cid:88) j ˜ M k,j ˙ θ j = ˜ V k , (3)where ˙ θ j = dθ j /dt and the coefficients are linear sumsof state overlaps that can be efficiently measured withquantum circuits [69]. We specify the detailed deriva-tion, expression of the coefficients, the quantum circuits,and a detailed resource estimation of the algorithm inSupplementary Materials [70]. Here we consider twoexamples with B ( t ) = 1 and | dv ( t ) (cid:105) = − iH | v ( t ) (cid:105) or | dv ( t ) (cid:105) = − ( H − (cid:104) v ( t ) | H | v ( t ) (cid:105) ) | v ( t ) (cid:105) , correspondingto real and imaginary time evolution [13, 20, 34],respectively. Compared to previous works studying realand imaginary time dynamics [13, 20, 34], our algorithmconsiders a much more general setting of first-orderdifferential equations. Therefore, it not only unifiesprevious results in the general setting, but provides thebasis for solving general problems as we discuss below. Variational algorithms for linear algebra.—
Now con-sider linear algebra problems of solving linear systems ofequations and matrix-vector multiplications. For a sparsesquare matrix M and a state vector | v (cid:105) , we aim to find | v M − (cid:105) = M − | v (cid:105) or | v M (cid:105) = M | v (cid:105) . (4)We introduce two types of algorithms where the first typeis more general and the second type is more efficient withassumptions of the matrix. Here we take matrix multipli-cation as an example and the derivation works similarlyfor linear equations, which can also be found in Supple-mentary Materials [70].For the first type, the algorithm is based on convert-ing the static algebraic problem into a dynamical pro-cess, evolving the initial vector | v (cid:105) to the target state | v M (cid:105) . A natural evolution path is via a linear ex-trapolation between | v (cid:105) and | v M (cid:105) as | v ( t ) (cid:105) = E ( t ) | v (cid:105) with E ( t ) = t/T · M + (cid:0) − t/T (cid:1) I , | v (0) (cid:105) = | v (cid:105) , and | v ( T ) (cid:105) = | v M (cid:105) . Different evolution paths can be alsoconsidered. For example, in the conventional Hamilto-nian simulation scenario, we have M = e − iHT and itcorresponds to an exponential extrapolation. We alsoconsider linear extrapolation between normalised statesin Supplementary Materials [70] and we leave the discus-sion of other possible evolution paths to future works.Given the evolution via linear extrapolation, the timederivative equation of | v ( t ) (cid:105) is ddt | v ( t ) (cid:105) = G | v (0) (cid:105) , with G = (cid:0) M − I (cid:1) /T . This corresponds to the case with A ( t ) = G , B ( t ) = 1, and | v (cid:48) ( t ) (cid:105) = | v (0) (cid:105) in Eq. (1),which can be variationally simulated.For the second type, we assume that M is given asa tensor product of matrices, M = M ⊗ · · · ⊗ M L ,with each M i acting on a small constant number ofqubits. We can thus sequentially act M i and focuson the realisation of each term. For each M ≡ M i ,we first consider a singular value decomposition as M = U DV , with unitary matrices U , V and diag-onal matrix D with non-negative entries. Now weshow how to realise the matrix multiplication of thethree matrices. Given a spectral decomposition of U = (cid:80) j e iλ j | λ j (cid:105) (cid:104) λ j | with λ j ∈ R , we can represent itas U = exp( − iH U T U ) with H U = − (cid:80) j λ j /T U | λ j (cid:105) (cid:104) λ j | and T U >
0, which can be realised by evolving the statewith Hamiltonian H U for time T U via variational realtime simulation [13]. The realisation is similar for V .To realise the diagonal matrix D = (cid:80) j a j | j (cid:105) (cid:104) j | ,we approximate it as D ≈ exp( − H D T D ) with − H D T = (cid:80) a j (cid:54) =0 log( a j ) | j (cid:105) (cid:104) j | − α (cid:80) a j =0 | j (cid:105) (cid:104) j | and aconstant α = O (log (1 /ε D ) + log Poly( n )) that ensuresan accuracy of ε D of the approximation with n qubits.We refer to Supplementary Materials for a detaileddiscussion [70]. Therefore, we can define an unnor-malised imaginary time evolution d | v ( τ ) (cid:105) dτ = − H D | v ( τ ) (cid:105) ,so that the initial vector | v (cid:105) at τ = 0 is evolved to D | v (cid:105) at τ = T . Note that even though the secondmethod assumes a tensor structure of M , the matrixmultiplication process can still be classically hard whenthe input state is a general multipartite entangled state.Such a case is practically relevant as we shortly discussfor its application in simulating jump processes of thestochastic Schr¨odinger equation. Open system simulation.—
Here we show quantum sim-ulation of open system dynamics. Conventional algo-rithms for this task generally require a deep circuit [61–66]. The recently proposed variational approach [34] alsoneeds two copies of the purified quantum state, thus re-quiring a number of qubits that is four times that ofthe system size. Our method, based on variational al-gorithms for generalised time evolution and matrix mul-tiplication and the stochastic Schr¨odinger equation, in-stead only requires to apply shallow circuits on a singlecopy of the state without purification. Our algorithmthus extensively alleviate the requirement of quantumsimulation of open systems.Consider the Lindblad master equation, ddt ρ = − i [ H, ρ ] + L ρ. (5)where H describes the system Hamiltonian and L ρ = (cid:80) k (2 L k ρL † k − L † k L k ρ − ρL † k L k ) describes the interac-tion with the environment with Lindblad operators L k .Instead of simulating the evolution of the density matrix,we consider its equivalent description by the stochasticSchr¨odinger equation, which averages trajectories of purestate evolution under continuous measurements [60, 67].Each single trajectory | ψ c ( t ) (cid:105) is described by d | ψ c ( t ) (cid:105) = (cid:32) − iH − (cid:88) k ( L † k L k − (cid:104) L † k L k (cid:105) ) (cid:33) | ψ c ( t ) (cid:105) dt + (cid:88) k (cid:20)(cid:18) L k | ψ c ( t ) (cid:105)|| L k | ψ c ( t ) (cid:105) || − | ψ c ( t ) (cid:105) (cid:19) dN k (cid:21) , (6)where d | ψ c ( t ) (cid:105) = | ψ c ( t + dt ) (cid:105) − | ψ c ( t ) (cid:105) , (cid:104) L † k L k (cid:105) = (cid:104) ψ c ( t ) | L † k L k | ψ c ( t ) (cid:105) , and dN k randomly takes either 0or 1, satisfing dN k dN k (cid:48) = δ kk (cid:48) dN k and E [ dN k ] = (cid:104) ψ c ( t ) | L † k L k | ψ c ( t ) (cid:105) dt . At each time t , we can as-sume a positive obervable valued measurement { O = I − (cid:80) k L † k L k dt, O k = L † k L k dt } happens. For measure-ment outcome O k , the state discontinuously jumps to L k | ψ c ( t ) (cid:105) / (cid:107) L k | ψ c ( t ) (cid:105) (cid:107) with probability E [ dN k ]. Andthe total jump probability is γ ( t ) = (cid:80) k E [ dN k ]. For out-come O with probability 1 − γ ( t ), we have dN k = 0 ∀ k ,and the state evolves under the generalised time evolu-tion Eq. (1) with operator A = − iH − (cid:88) k ( L † k L k − (cid:104) L † k L k (cid:105) ) , (7)and B ( t ) = 1. Here, − iH corresponds to the conven-tional real time Schr¨odinger evolution with Hamiltonian H and the other terms can be understood as a normaliseddamping process. Therefore, the whole process is com-posed of two parts: the continuous generalised time evo-lution and the quantum jump process.The stochastic Schr¨odinger equation can be simulatedwith the Monto Carlo method. Suppose the state jumpsat time t , then at time t + τ , the probability p ( t + τ )that the state does not jump is p ( t + τ ) = e − Γ( t,τ ) , withΓ( t, τ ) = (cid:82) t + τt γ ( t (cid:48) ) dt (cid:48) . When a jump happens at time t ,a uniform random number q ∈ [0 ,
1] is generated. Thenthe time of the next jump is determined by accumulat-ing time τ until we have p ( t + τ ) = q . When a jumphappens, a random number q (cid:48) ∈ [0 ,
1] is generated to de-termine which jump operator to apply at each timestep.The state is updated to L k | ψ c ( t ) (cid:105) / (cid:107) L k | ψ c ( t ) (cid:105) (cid:107) if q (cid:48) ∈ [˜ γ k − ( t ) , ˜ γ k ( t )], where˜ γ k ( t ) = (cid:80) kl =1 (cid:104) ψ c ( t ) | L † l L l | ψ c ( t ) (cid:105) (cid:80) N L l =1 (cid:104) ψ c ( t ) | L † l L l | ψ c ( t ) (cid:105) , (8)and N L is the number of the Lindblad operators. Con-sidering discretised time with initial state | ψ c (0) (cid:105) , thestochastic Schr¨odinger equation from time 0 to T can besimulated as follows. Algorithm 1
Stochastic evolution equation Set Γ = 0 and generate a random number q ∈ [0 , for t = 0 : dt : T do if e − Γ ≥ q then Evolve the state | ψ c ( t ) (cid:105) under A in Eq. (7). Calculate γ ( t ) = (cid:80) k (cid:104) ψ c ( t ) | L † k L k | ψ c ( t ) (cid:105) dt . Update Γ = Γ + γ ( t ). else Calculate ˜ γ k ( t ) in Eq. (8) Generate a random number q (cid:48) ∈ [0 , if q (cid:48) ∈ [˜ γ k − ( t ) , ˜ γ k ( t )] then Update | ψ c ( t ) (cid:105) to L k | ψ c ( t ) (cid:105) / (cid:107) L k | ψ c ( t ) (cid:105) (cid:107) . Reset Γ = 0 and randomly generate q ∈ [0 , Now we show how to variationally simulate thestochastic Schr¨odinger equation, Algorithm 1. Sup-pose the state | ψ c ( t ) (cid:105) at time t can be representedby the parametrised state | φ c ( (cid:126)θ ( t )) (cid:105) prepared by aquantum computer. We can simulate step 4, i.e.,the evolution under operator A defined in Eq. (7),with the algorithm for generalised time evolution.Specifically, we can evolve the parameters accord-ing to Eq. (3) with ˜ M k,j = Re (cid:16) ∂ (cid:104) ϕ ( (cid:126)θ ( t )) | ∂θ k ∂ | ϕ ( (cid:126)θ ( t )) (cid:105) ∂θ j (cid:17) ,˜ V k = Re (cid:16) (cid:104) ϕ ( (cid:126)θ ( t )) | ( − iH − ( L − (cid:104) L (cid:105) )) ∂ | ϕ ( (cid:126)θ ( t )) (cid:105) ∂θ k (cid:17) , and L = (cid:80) k L † k L k . The values of γ ( t ) and ˜ γ k ( t ) at step 5and 8 are measured as expectation values of quantumstates. The jump at step 11 is realised by the variationalalgorithms for matrix-vector multiplication. Especially,when considering L k as a product of operators of eachqubits, it can be efficiently realised with the singular FIG. 1. Geometry of the 2D Ising model. We considernearest-neighbor interactions for two connected qubits. value decomposition method. In practice, we considersparse Hamiltonian and Lindblad operators, thereforeall the measurements can be efficiently evaluated. Werefer to Supplementary Materials for the discussion ofthe resource estimation [70].
Numerical simulation.—
We numerically test the varia-tional algorithm for simulating a six-qubit 2D Ising modelin a transverse field coupled to a Markovian bath [71–74].The Hamiltonian is H I = J/ (cid:80) (cid:104) ij (cid:105) Z i Z j + h X (cid:80) i =1 X i ,where Pauli operators X i , Y i , and Z i act on the i th spinand (cid:104) ij (cid:105) represents nearest neighbor pairs in Fig. 1. TheLindblad term is L i = √ γσ + i with σ + i = | (cid:105) (cid:104) | i beingthe raising operator acting on the i th spin. In our sim-ulation, we set J = 1, h X = 1, and γ = 1. The ini-tial state is prepared in a product state | ϕ (0) (cid:105) = | (cid:105) ⊗ ,and then the Hamiltonian H I is suddenly turned on,which drives the qubits out of equilibrium. We simu-late both the ideal and dissipative evolution from t = 0to t = 6, and measure the normalised nearest-neighborcorrelations C = (cid:80) (cid:104) ij (cid:105) Z i Z j /
7. We use the Hamiltonianansatz (HA) [75] sandwiched with single qubit rotationsas the trial state, shown Fig. 2. The HA preserves thesymmetry of the Hamiltonian and the single qubit rota-tions are introduced to break the symmetry, ensuring itscapability for simulating the jump process.We apply the algorithm in Ref. [13] and our Algo-rithm 1 to simulate the ideal and dissipative evolution, | (cid:105) HA HA HA R X HA HA HA R X | (cid:105) R X R X | (cid:105) R X R X | (cid:105) R X R X | (cid:105) R X R X | (cid:105) R X R X FIG. 2. The ansatz consists of the Hamiltonianansatz (HA) and single qubit rotations. Each HA is e θ ( Z Z ) e θ ( Z Z + Z Z ) e θ ( Z Z ) e θ ( Z Z + Z Z ) e θ ( Z Z ) e θ ( X + X ) e θ ( X + X + X + X ) and each single qubit gate R X is e − iθX . With different parameters for different HA andsingle qubit gates, the ansatz has in total 54 parameters. Time C o rr e l a ti on Ideal & exactIdeal & variationalDissipative & exactDissipative & variational
Time -0.0100.01 E rr o r FIG. 3. Numerical simulation of ideal and dissipative evolu-tion of the 2D Ising model. We consider evolution from t = 0to 6 with δt = 0 .
005 and measure the nearest-neighbor corre-lations C = (cid:80) (cid:104) ij (cid:105) Z i Z j /
7. We repeat the stochastic processof Algorithm 1 with N trial = 2 × times. The results fromvariational simulation agree well with the exact ones, with amaximal error around 0.01 (see inset). respectively. We use the singular value decompositionmethod to simulate quantum jumps in Algorithm 1. Wedecompose the jump operator σ + as | (cid:105) (cid:104) | = | (cid:105) (cid:104) | X with | (cid:105) (cid:104) | = U DV , U = I , D = | (cid:105) (cid:104) | , and V = X . Torealise the V operator, we set H V = X and T V = π/ X = exp( − iH V T V ). Then we evolve thestate under Hamiltonian H V for time T V with time step δt V = 0 .
01 to have X | ϕ ( (cid:126)θ ) (cid:105) . To realise D = | (cid:105) (cid:104) | , weset H D = | (cid:105) (cid:104) | and T D = 6 so that D ≈ exp( − H D T D ).Then we realise D | ϕ ( (cid:126)θ ) (cid:105) / (cid:107) D | ϕ ( (cid:126)θ ) (cid:105) (cid:107) by using the nor-malised variational imaginary time evolution with totaltime T D and time step δt D = 0 . Discussion.—
To summarise, we extend the variationalquantum simulation method to general processes, includ-ing generalised time evolution, and its application forsolving linear algebra tasks and simulating the evolutionof open quantum systems. Our algorithm for simulat-ing the generalised time evolution can be applied to simu-late non-Hermitian quantum mechanics [35–37] includingnonequilibrium processes [38] and parity-time symmetricHamiltonians [39–41]. Especially, it is shown in Ref. [41]that a quantum state can evolve to the target state fasterwith non-Hermitian parity-time symmetric Hamiltoniansthan the case with Hermitian Hamiltonians. Therefore,our variational algorithm for simulating the generalisedtime evolution may also be useful for designing fasterquantum computing algorithms. Recently, several otheralgorithms for linear algebra tasks have been proposedto be compatible with near-term quantum hardware [77–80]. A thorough comparison of these algorithms and ourscould also be an interesting future work. Our algorithmfor simulating open systems only needs shallow circuitson a number of qubits matching the system size, thusit enables the possibility of investigating open systemphysics with near-term quantum computers. Finally theproposed algorithms are compatible with NISQ hardwareand can be further combined with quantum error miti-gation techniques [13, 17, 81–87].
Acknowledgements.—
This work is supported bythe EPSRC National Quantum Technology Hubin Networked Quantum Information Technology(EP/M013243/1). SCB acknowledges support from theEU AQTION project. We acknowledge the use of theUniversity of Oxford Advanced Research Computing(ARC) facility in carrying out this work. SE is supportedby Japan Student Services Organization (JASSO) Stu-dent Exchange Support Program (Graduate Scholarshipfor Degree Seeking Students). YL is supported byNational Natural Science Foundation of China (GrantNo. 11875050) and NSAF (Grant No. U1930403). JSacknowledges support from Chinese Scholarship Council. ∗ [email protected] † [email protected][1] R. Balian and M. V´en´eroni, Annals of Physics , 29(1988).[2] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari,Rev. Mod. Phys. , 463 (1999).[3] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Piˇzorn, H. Ver-schelde, and F. Verstraete, Phys. Rev. Lett. , 070601(2011).[4] T. Shi, E. Demler, and J. I. Cirac, Annals of Physics , 245 (2018).[5] L. Vanderstraeten, J. Haegeman, and F. Verstraete, Sci-Post Phys. Lect. Notes (2019).[6] E. Farhi, J. Goldstone, and S. Gutmann, arXiv preprintarXiv:1411.4028 (2014).[7] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung,Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien,Nature communications (2014).[8] Y. Wang, F. Dolde, J. Biamonte, R. Babbush,V. Bergholm, S. Yang, I. Jakobi, P. Neumann,A. Aspuru-Guzik, J. D. Whitfield, et al. , ACS nano ,7769 (2015).[9] P. J. J. O’Malley, R. Babbush, I. D. Kivlichan,J. Romero, J. R. McClean, R. Barends, J. Kelly,P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jef-frey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley,C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wen-ner, T. C. White, P. V. Coveney, P. J. Love, H. Neven,A. Aspuru-Guzik, and J. M. Martinis, Phys. Rev. X ,031007 (2016).[10] Y. Shen, X. Zhang, S. Zhang, J.-N. Zhang, M.-H. Yung,and K. Kim, Phys. Rev. A , 020501 (2017).[11] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New Journal of Physics , 023023 (2016).[12] S. Paesani, A. A. Gentile, R. Santagati, J. Wang,N. Wiebe, D. P. Tew, J. L. O’Brien, and M. G. Thomp-son, Phys. Rev. Lett. , 100503 (2017).[13] Y. Li and S. C. Benjamin, Phys. Rev. X , 021050 (2017).[14] J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok,M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A.de Jong, and I. Siddiqi, Phys. Rev. X , 011021 (2018).[15] R. Santagati, J. Wang, A. A. Gentile, S. Paesani,N. Wiebe, J. R. McClean, S. Morley-Short, P. J. Shad-bolt, D. Bonneau, J. W. Silverstone, D. P. Tew, X. Zhou,J. L. O’Brien, and M. G. Thompson, Science Advances (2018), 10.1126/sciadv.aap9646.[16] A. Kandala, A. Mezzacapo, K. Temme, M. Takita,M. Brink, J. M. Chow, and J. M. Gambetta, Nature , 242 (2017).[17] A. Kandala, K. Temme, A. D. C´orcoles, A. Mezzacapo,J. M. Chow, and J. M. Gambetta, Nature , 491(2019).[18] C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz,H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush,A. Aspuru-Guzik, R. Blatt, and C. F. Roos, Phys. Rev.X , 031022 (2018).[19] J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J.Love, and A. Aspuru-Guzik, Quantum Science and Tech-nology , 014008 (2018).[20] S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin,and X. Yuan, npj Quantum Information , 1 (2019).[21] T. Jones, S. Endo, S. McArdle, X. Yuan, and S. C.Benjamin, Physical Review A , 062304 (2019).[22] O. Higgott, D. Wang, and S. Brierley, Quantum , 156(2019).[23] R. Santagati, J. Wang, A. A. Gentile, S. Paesani,N. Wiebe, J. R. McClean, S. Morley-Short, P. J. Shad-bolt, D. Bonneau, J. W. Silverstone, et al. , Science ad-vances , eaap9646 (2018).[24] J. R. McClean, M. E. Kimchi-Schwartz, J. Carter, andW. A. de Jong, Physical Review A , 042308 (2017).[25] J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok,M. Kimchi-Schwartz, J. McClean, J. Carter, W. De Jong,and I. Siddiqi, Physical Review X , 011021 (2018).[26] C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K.Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F.Roos, et al. , Nature , 355 (2019).[27] I. Kassal, J. D. Whitfield, A. Perdomo-Ortiz, M.-H.Yung, and A. Aspuru-Guzik, Annual review of physi-cal chemistry , 185 (2011).[28] D. Lu, B. Xu, N. Xu, Z. Li, H. Chen, X. Peng, R. Xu,and J. Du, Phys. Chem. Chem. Phys. , 9411 (2012).[29] K. B. Whaley, A. R. Dinner, and S. A. Rice, Quantuminformation and computation for chemistry (John Wiley& Sons, 2014).[30] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin,and X. Yuan, Reviews of Modern Physics , 015003(2020). [31] J. Preskill, Quantum , 79 (2018).[32] K. M. Nakanishi, K. Mitarai, and K. Fujii, Physical Re-view Research , 033062 (2019).[33] K. Heya, K. M. Nakanishi, K. Mitarai, and K. Fujii,arXiv preprint arXiv:1904.08566 (2019).[34] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin,Quantum , 191 (2019).[35] N. Hatano and D. R. Nelson, Phys. Rev. Lett. , 570(1996).[36] N. Moiseyev, Non-Hermitian Quantum Mechanics (Cam-bridge University Press, 2011).[37] F. Bagarello, J.-P. Gazeau, F. H. Szafraniec, and M. Zno-jil,
Non-selfadjoint operators in quantum physics: Math-ematical aspects (John Wiley & Sons, 2015).[38] H. C. Fogedby, A. B. Eriksson, and L. V. Mikheev, Phys.Rev. Lett. , 1883 (1995).[39] C. M. Bender and S. Boettcher, Phys. Rev. Lett. , 5243(1998).[40] C. M. Bender, Reports on Progress in Physics , 947(2007).[41] C. M. Bender, D. C. Brody, H. F. Jones, and B. K.Meister, Physical Review Letters , 040403 (2007).[42] I. Rotter, Journal of Physics A: Mathematical and The-oretical , 153001 (2009).[43] C. E. Rasmussen, in Advanced lectures on machine learn-ing (Springer, 2004) pp. 63–71.[44] N. M. Nasrabadi, Journal of electronic imaging ,049901 (2007).[45] A. W. Harrow, A. Hassidim, and S. Lloyd, Phys. Rev.Lett. , 150502 (2009).[46] A. Childs, R. Kothari, and R. Somma,SIAM Journal on Computing , 1920 (2017),https://doi.org/10.1137/16M1087072.[47] A. Ambainis, in STACS’12 (29th Symposium on Theoret-ical Aspects of Computer Science) , Vol. 14 (LIPIcs, 2012)pp. 636–647.[48] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Phys.Rev. Lett. , 250504 (2013).[49] L. Wossnig, Z. Zhao, and A. Prakash, Phys. Rev. Lett. , 050502 (2018).[50] S. Chakraborty, A. Gily´en, and S. Jeffery, in , Leibniz International Pro-ceedings in Informatics (LIPIcs), Vol. 132, editedby C. Baier, I. Chatzigiannakis, P. Flocchini, andS. Leonardi (Schloss Dagstuhl–Leibniz-Zentrum fuer In-formatik, Dagstuhl, Germany, 2019) pp. 33:1–33:14.[51] A. Gily´en, Y. Su, G. H. Low, and N. Wiebe, in
Proceed-ings of the 51st Annual ACM SIGACT Symposium onTheory of Computing (2019) pp. 193–204.[52] Y. b. u. Suba¸s ı, R. D. Somma, and D. Orsucci, Phys.Rev. Lett. , 060504 (2019).[53] P. Rebentrost, M. Mohseni, and S. Lloyd, Phys. Rev.Lett. , 130503 (2014).[54] S. Lloyd, M. Mohseni, and P. Rebentrost, Nature Physics , 631 (2014).[55] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost,N. Wiebe, and S. Lloyd, Nature , 195 (2017).[56] C. Shao, arXiv e-prints , arXiv:1803.01601 (2018),arXiv:1803.01601 [quant-ph].[57] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Phys-ical Review A , 032309 (2018).[58] G. Lindblad, Communications in Mathematical Physics , 119 (1976). [59] H.-P. Breuer, F. Petruccione, et al. , The theory of openquantum systems (Oxford University Press on Demand,2002).[60] C. Gardiner and P. Zoller, “Quantum noise: A handbookof markovian and non-markovian stochastic process withapplications to quantum optics,” (2000).[61] H. Wang, S. Ashhab, and F. Nori, Physical Review A , 062317 (2011).[62] D. Bacon, A. M. Childs, I. L. Chuang, J. Kempe, D. W.Leung, and X. Zhou, Physical Review A , 062302(2001).[63] R. Sweke, I. Sinayskiy, D. Bernard, and F. Petruccione,Physical Review A , 062308 (2015).[64] R. Sweke, M. Sanz, I. Sinayskiy, F. Petruccione, andE. Solano, Physical Review A , 022317 (2016).[65] R. Cleve and C. Wang, arXiv preprint arXiv:1612.09512(2016).[66] A. M. Childs and T. Li, Quantum Information & Com-putation , 901 (2017).[67] H. Carmichael, Physical review letters , 2273 (1993).[68] A. McLachlan, Molecular Physics , 39 (1964).[69] K. Mitarai and K. Fujii, Physical Review Research ,013006 (2019).[70] See Supplemental Materials for resource estimations anddetailed derivations.[71] T. E. Lee, H. Haeffner, and M. Cross, Physical reviewletters , 023602 (2012).[72] C. Ates, B. Olmos, J. P. Garrahan, and I. Lesanovsky,Physical Review A , 043620 (2012).[73] I. Lesanovsky, M. van Horssen, M. Gut¸˘a, and J. P. Gar-rahan, Physical review letters , 150401 (2013).[74] D. C. Rose, K. Macieszczak, I. Lesanovsky, and J. P.Garrahan, Physical Review E , 052132 (2016).[75] D. Wecker, M. B. Hastings, and M. Troyer, Phys. Rev.A , 042303 (2015).[76] J. Raftery, D. Sadri, S. Schmidt, H. E. T¨ureci, and A. A.Houck, Physical Review X , 031043 (2014).[77] X. Xu, J. Sun, S. Endo, Y. Li, S. C. Benjamin, andX. Yuan, arXiv preprint arXiv:1909.03898 (2019).[78] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Sub-asi, L. Cincio, and P. J. Coles, arXiv preprintarXiv:1909.05820 (2019).[79] D. An and L. Lin, arXiv preprint arXiv:1909.05500(2019).[80] H.-Y. Huang, K. Bharti, and P. Rebentrost, arXivpreprint arXiv:1909.07344 (2019).[81] S. McArdle, X. Yuan, and S. Benjamin, Physical ReviewLetters , 180501 (2019).[82] X. Bonet-Monroig, R. Sagastizabal, M. Singh, andT. O’Brien, Physical Review A , 062339 (2018).[83] K. Temme, S. Bravyi, and J. M. Gambetta, Phys. Rev.Lett. , 180509 (2017).[84] S. Endo, S. C. Benjamin, and Y. Li, Physical Review X , 031027 (2018).[85] M. Huo and Y. Li, arXiv preprint arXiv:1811.02734(2018).[86] M. Otten and S. K. Gray, npj Quantum Information ,11 (2019).[87] J. Sun, X. Yuan, T. Tsunoda, V. Vedral, S. C. Bejamin,and S. Endo, arXiv preprint arXiv:2001.04891 (2020).[88] S. Lloyd, Science , 1073 (1996).[89] M. Suzuki, Journal of Mathematical Physics , 400(1991).[90] A. M. Childs, A. Ostrander, and Y. Su, Quantum , 182 (2019).[91] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Physical review letters , 090502 (2015).[92] G. H. Low and I. L. Chuang, Quantum , 163 (2019). SUPPLEMENTARY MATERIALS: VARIATIONAL QUANTUM SIMULATION OF GENERALPROCESSESREVIEW: VARIATIONAL QUANTUM SIMULATION OF REAL AND IMAGINARY TIME EVOLUTION
We first review the variational quantum algorithms for simulating real and imaginary time evolution introduced inRef. [13] and [20], respectively, as we generalise these algorithms to propose our new algorithm for general processes.We refer the reader to Ref. [34] for a comprehensive derivation of quantum variational algorithms from variousvariational principles — the Dirac and Frenkel variational principle, the McLachlan’s variational principle, and thetime-dependent variational principle.The real time evolution is described by the Schr¨odinger equation, d | ψ ( t ) (cid:105) dt = − iH | ψ ( t ) (cid:105) , (9)with Hermitian Hamiltonian H . Instead of directly simulating the real time dynamics with the Hamiltonian simulationalgorithms [88–92], the variational method assumes that the quantum state | ψ ( t ) (cid:105) is prepared by a parametrisedquantum circuit, | ϕ ( (cid:126)θ ( t )) (cid:105) = R N ( θ N ) . . . R k ( θ k ) . . . R ( θ ) | ¯0 (cid:105) with each gate R k ( θ k ) controlled by the real parameter θ k and the reference state | ¯0 (cid:105) . Here, we denote (cid:126)θ = ( θ , θ , . . . , θ N ). According to McLachlan’s variational principle [68],the real time dynamics of | ψ ( t ) (cid:105) can be mapped to the evolution of the parameters (cid:126)θ ( t ) by minimising the distancebetween the ideal evolution and the evolution induced of the parametrised trial state, δ (cid:107) ( ∂/∂t + iH ) | ϕ ( (cid:126)θ ( t )) (cid:105) (cid:107) = 0 , (10)where (cid:107) | ϕ (cid:105) (cid:107) = (cid:112) (cid:104) ϕ | ϕ (cid:105) . The solution is (cid:88) j M k,j ˙ θ j = V k , (11)with coefficients M k,j = Re (cid:18) ∂ (cid:104) ϕ ( (cid:126)θ ( t )) | ∂θ k ∂ | ϕ ( (cid:126)θ ( t )) (cid:105) ∂θ j (cid:19) ,V k = Im (cid:18) (cid:104) ϕ ( (cid:126)θ ( t )) | H ∂ | ϕ ( (cid:126)θ ( t )) (cid:105) ∂θ k (cid:19) . (12)For imaginary time evolution, the normalised Wick-rotated Shr¨odinger equation is obtained by replacing t = iτ inEq. (9), d | ψ ( τ ) (cid:105) dτ = − ( H − (cid:104) ψ ( τ ) | H | ψ ( τ ) (cid:105) ) | ψ ( τ ) (cid:105) . (13)Applying a similar procedure for real time evolution, the imaginary time evolution is mapped to the evolution of theparameters via McLachlan’s principle, δ (cid:107) ( ∂/∂τ + H − (cid:104) H (cid:105) ) | ϕ ( (cid:126)θ ( t )) (cid:105) (cid:107) = 0 . (14)The evolution of the parameters is (cid:88) j M k,j ˙ θ j = C k , (15)with M given in Eq. (12) and C defined by C k = − Re (cid:32) (cid:104) ϕ ( (cid:126)θ ( t )) | H ∂ | ϕ ( (cid:126)θ ( t )) (cid:105) ∂θ k (cid:33) . (16)The M , V , and C terms can be efficiently measured with quantum circuits. Considering gate based circuits, thederivative of the each parameterised gate can be expressed as ∂R k ∂θ k = (cid:88) i g k,i R k σ k,i , (17)where σ k,i are unitary operators and g k,i are complex coefficients. The derivative of the trial state can be written as ∂ | ϕ ( (cid:126)θ ( t )) (cid:105) ∂θ k = (cid:88) i g k,i R k,i | ¯0 (cid:105) , (18)where R k,i = R N R N − · · · R k +1 R k σ k,i · · · R R . (19)The M k,j terms can be expressed as M k,j = (cid:88) i,j (cid:60) (cid:16) g ∗ k,p g j,q (cid:104) ¯0 | R † k,p R j,q | ¯0 (cid:105) (cid:17) . (20)Similarly, considering sparse Hamiltonian with decomposition H = (cid:80) j λ j σ j , λ j ∈ R , we have C k and V k as V k = (cid:88) i,j Re (cid:16) ig ∗ k,i λ j (cid:104) ¯0 | R † k,i σ j R | ¯0 (cid:105) (cid:17) ,C k = − (cid:88) i,j Re (cid:16) g ∗ k,i λ j (cid:104) ¯0 | R † k,i σ j R | ¯0 (cid:105) (cid:17) , (21)All the M , C , and V terms can be written in the form a Re (cid:0) e iθ (cid:104) ¯0 | U | ¯0 (cid:105) (cid:1) , where a, θ ∈ R depend on the coefficients, and U is a unitary operator of either R † k,p R j,q or R † k,i σ j R . We can calculate M , C , and V by using the quantum circuit shown in Fig. 4. ( | (cid:105) + e iθ | (cid:105) ) / √ X • X • H... ... | ¯0 (cid:105) R R k − σ k,p R k R j − σ j,q (a)( | (cid:105) + e iθ | (cid:105) ) / √ X • X • H... ... | ¯0 (cid:105) R R k − σ k,i R k R N σ j (b)FIG. 4. Quantum circuits that evaluate (a) Re( e iθ (cid:104) ¯0 | R † k,p R j,q | ¯0 (cid:105) ) and (b) Re( e iθ (cid:104) ¯0 | R † k,i σ j R | ¯0 (cid:105) ). DERIVATION FOR VARIATIONAL SIMULATION OF GENERALISED TIME EVOLUTION EQUATION
Now, we consider variational simulation of the generalised time evolution equation, B ( t ) ddt | v ( t ) (cid:105) = (cid:88) j A j ( t ) | v (cid:48) j ( t ) (cid:105) . (22)By parametrising | v ( t ) (cid:105) and | v (cid:48) j ( t ) (cid:105) as | v ( (cid:126)θ ( t )) (cid:105) and | v (cid:48) ( (cid:126)θ (cid:48) j ( t )) (cid:105) , with McLachlan’s principle, we have δ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) B ( t ) (cid:88) i ∂ | v ( (cid:126)θ ( t )) (cid:105) ∂θ i ˙ θ i − (cid:88) j A j ( t ) | v (cid:48) ( (cid:126)θ (cid:48) j ( t )) (cid:105) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) = 0 . (23)This is equivalent to ∂∂ ˙ θ k (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) B ( t ) (cid:88) i ∂ | v ( (cid:126)θ ( t )) (cid:105) ∂θ i ˙ θ i − (cid:88) j A j ( t ) | v (cid:48) ( (cid:126)θ (cid:48) j ( t )) (cid:105) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) = ∂∂ ˙ θ k (cid:88) i ∂ (cid:104) v ( (cid:126)θ ( t )) | ∂θ i ˙ θ i B † ( t ) − (cid:88) j (cid:104) v (cid:48) ( (cid:126)θ (cid:48) j ( t )) | A † j ( t ) (cid:88) l B ( t ) ∂ | v ( (cid:126)θ ( t )) (cid:105) ∂θ l ˙ θ l − (cid:88) j A j ( t ) | v (cid:48) ( (cid:126)θ (cid:48) j ( t )) (cid:105) = 0 ∀ k. (24)Hence, we have (cid:88) j (cid:32) ∂ (cid:104) v ( (cid:126)θ ( t )) | ∂θ k B † ( t ) B ( t ) ∂ | v ( (cid:126)θ ( t )) (cid:105) ∂θ j + ∂ (cid:104) v ( (cid:126)θ ( t )) | ∂θ j B † ( t ) B ( t ) ∂ | v ( (cid:126)θ ( t )) (cid:105) ∂θ k (cid:33) ˙ θ j = (cid:88) j ∂ (cid:104) v ( (cid:126)θ ( t )) | ∂θ k B † ( t ) A j ( t ) | v (cid:48) j ( (cid:126)θ (cid:48) j ( t )) (cid:105) + h.c, (25)which leads to Eq. (3) in the main text. (cid:88) j ˜ M k,j ˙ θ j = ˜ V k . By substituting | v ( (cid:126)θ ( t )) (cid:105) = α ( (cid:126)θ ( t )) | ϕ ( (cid:126)θ ( t )) (cid:105) and | v (cid:48) ( (cid:126)θ (cid:48) j ( t )) (cid:105) = α (cid:48) ( (cid:126)θ (cid:48) j ( t )) | ϕ ( (cid:126)θ (cid:48) j ( t )) (cid:105) , we have˜ M k,j = Re (cid:32) | α ( (cid:126)θ ( t )) | ∂ (cid:104) ϕ ( (cid:126)θ ( t )) | ∂θ k B † ( t ) B ( t ) ∂ | ϕ ( (cid:126)θ ( t )) (cid:105) ∂θ j (cid:33) + Re (cid:32) ∂α ∗ ( (cid:126)θ ( t )) ∂θ k α ( (cid:126)θ ( t )) (cid:104) ϕ ( (cid:126)θ ( t )) | B ( t ) ∂ | ϕ ( (cid:126)θ ( t )) (cid:105) ∂θ j (cid:33) + Re (cid:32) ∂α ∗ ( (cid:126)θ ( t )) ∂θ j α ( (cid:126)θ ( t )) (cid:104) ϕ ( (cid:126)θ ( t )) | B ( t ) ∂ | ϕ ( (cid:126)θ ( t )) (cid:105) ∂θ k (cid:33) + Re (cid:32) ∂α ( (cid:126)θ ( t )) ∂θ k ∂α ∗ ( (cid:126)θ ( t )) ∂θ j (cid:104) ϕ ( (cid:126)θ ( t )) | B † ( t ) B ( t ) | ϕ ( (cid:126)θ ( t )) (cid:105) (cid:33) , ˜ V k = (cid:88) j Re (cid:32) ∂α ∗ ( (cid:126)θ ( t )) ∂θ k α (cid:48) ( (cid:126)θ (cid:48) j ( t )) (cid:104) ϕ ( (cid:126)θ ( t )) | B † ( t ) A j ( t ) | ϕ (cid:48) j ( (cid:126)θ (cid:48) j ( t )) (cid:105) (cid:33) . + (cid:88) j Re (cid:32) α ∗ ( (cid:126)θ ( t )) α (cid:48) ( (cid:126)θ (cid:48) j ( t )) ∂ (cid:104) ϕ ( (cid:126)θ ( t )) | ∂θ k B † ( t ) A j ( t ) | ϕ (cid:48) j ( (cid:126)θ (cid:48) j ( t )) (cid:105) (cid:33) . (26)0The first term of ˜ M k,j can written as (cid:88) i,q,l Re (cid:18) | α ( (cid:126)θ ( t )) | g ∗ k,i g j,q β l (cid:104) ¯0 | R † k,i σ l R j,q | ¯0 (cid:105) (cid:19) , (27)where we set B † ( t ) B ( t ) = (cid:80) l β l σ l , and σ l is a Pauli operator. Each term of ˜ M k,j can be written in the form of a Re( e iθ (cid:104) ¯0 | R † k,i σ l R j,q | ¯0 (cid:105) ), where a, θ ∈ R . The quantum circuit to compute this value is shown in Fig. 5. It is worthmentioning that this quantum circuit only necessitates three controlled operations from an ancilla qubit. The secondand third terms can be written as (cid:80) ij Re (cid:0) ∂α ∗ ( (cid:126)θ ( t )) ∂θ k α ( (cid:126)θ ( t )) γ j g ∗ k,i (cid:104) ¯0 | R † k,i σ j R | ¯0 (cid:105) (cid:1) , where we set B ( t ) = (cid:80) j γ j σ j .We can express each term in the form of a Re( e iθ (cid:104) ¯0 | R † k,i σ l R | ¯0 (cid:105) ), which can be evaluated by using the quantumcircuit shown in Fig. 4 (b). The fourth term can be simply computed by measuring the expectation value of B † ( t ) B ( t ) = (cid:80) l β l σ l for | ϕ ( (cid:126)θ ) (cid:105) .The term ˜ V can be computed as follows. We denote ˜ U t and ˜ U (cid:48) t to be the unitary circuit to prepare | ϕ ( (cid:126)θ ( t )) (cid:105) = ˜ U t | ¯0 (cid:105) and | ϕ (cid:48) j ( (cid:126)θ (cid:48) j ( t )) (cid:105) = ˜ U (cid:48) t | ¯0 (cid:105) . Replacing B † ( t ) A j ( t ) = (cid:80) l Λ jl ( t ) σ l , the first term of each˜ V k can be written as (cid:80) jl Re (cid:0) α ∗ ( (cid:126)θ ( t )) α (cid:48) ( (cid:126)θ (cid:48) j ( t ))Λ jl (cid:104) ¯0 | ˜ U † t σ l ˜ U (cid:48) t | ¯0 (cid:105) (cid:1) . Each term can be expressed in the form of a Re (cid:0) e iθ (cid:104) ¯0 | ˜ U † t σ l ˜ U (cid:48) t | ¯0 (cid:105) (cid:1) and computed using the quantum circuit shown in Fig. 6. Meanwhile, the second term canbe expressed as (cid:80) ijl Re( α ∗ ( (cid:126)θ ( t )) α (cid:48) ( (cid:126)θ j ( t ))Λ jl ( t ) g ∗ k,i (cid:104) ¯0 | R † k,i σ l ˜ U (cid:48) t | ¯0 (cid:105) ), and each term can be written in the form of a Re (cid:0) e iθ (cid:104) ¯0 | R † k,i σ l ˜ U (cid:48) t | ¯0 (cid:105) (cid:1) , which can be computed by using the quantum circuit shown in Fig. 7.Note that, although the quantum circuits to evaluate ˜ V generally need controlled U operations, in the case | v (cid:48) ( (cid:126)θ (cid:48) j ( t )) (cid:105) = | v ( (cid:126)θ ( t )) (cid:105) , e.g., open quantum simulation and solving linear equations, it can be computed by the quantumcircuit shown in Fig. 4 (b), which only needs two controlled operations. ( | (cid:105) + e iθ | (cid:105) ) / √ • X • • X H... ... ... | ¯0 (cid:105) R R k − σ k,i R k R j − σ j,q R j +1 R N σ l FIG. 5. Quantum circuits that evaluate Re (cid:18) g ∗ k,i g j,q (cid:104) ¯0 | R † k,i σ l R j,q | ¯0 (cid:105) (cid:19) . | (cid:105) + e iφ l | (cid:105) • • X • X H | ¯0 (cid:105) ˜ U t σ l ˜ U (cid:48) t FIG. 6. The quantum circuit for evaluating ˜ V k . | (cid:105) + e iφ l | (cid:105) • • X • X H | ¯0 (cid:105) R k,i σ l ˜ U t FIG. 7. The quantum circuit for evaluating ˜ V k . VARIATIONAL SIMULATION FOR LINEAR ALGEBRA TASKSMatrix evolution with normalised states
Here, we also consider the case where we are only interested in the normalised final state | ψ ( t ) (cid:105) = M | v (cid:105) / (cid:107)M | v (cid:105) (cid:107) .By extrapolating from | v (cid:105) / (cid:107) | v (cid:105) (cid:107) to | ψ ( t ) (cid:105) , we can similarly have an evolution of the state | ψ ( t ) (cid:105) as | ψ ( t ) (cid:105) = E (cid:48) ( t ) | ψ (cid:105) , (28)with E (cid:48) ( t ) = N ( t ) (cid:18) tT M + (1 − tT ) I (cid:19) , (29)and a normalisation factor N ( t ) = 1 (cid:115)(cid:13)(cid:13)(cid:13)(cid:13)(cid:18) tT M + (1 − tT ) I (cid:19) | ψ (cid:105) (cid:13)(cid:13)(cid:13)(cid:13) . (30)The normalisation factor N ( t ) can be measured from the expectation values of M † + M and M † M for | ψ (cid:105) . Giventhe definition of the state | ψ ( t ) (cid:105) at time t , the corresponding derivative equation is ddt | ψ ( t ) (cid:105) = ˙ N ( t ) N ( t ) | ψ ( t ) (cid:105) + N ( t ) G | ψ (0) (cid:105) . (31)Such an equation is also described by the generalised time evolution equation with | v ( t ) (cid:105) = | ψ ( t ) (cid:105) , A ( t ) = ˙ N ( t ) N ( t ) I , A ( t ) = N ( t ) G , | v (cid:48) ( t ) (cid:105) = | ψ ( t ) (cid:105) , | v (cid:48) ( t ) (cid:105) = | ψ (0) (cid:105) , and B ( t ) = 1. Solving linear equations
We now discuss how to solve linear equations, defined by
M | v M − (cid:105) = | v (cid:105) , with vector | v (cid:105) , invertible matrix M , and solution | v M − (cid:105) = M − | v (cid:105) . We also introduce two algorithms for thisproblem where the first does not assume the structure of the matrix M and the second one assumes tensor productstructure of M .For the first type of algorithm, we convert the static problem into a dynamical problem. We consider a path fromthe initial state | v (cid:105) to | v M − (cid:105) as E ( t ) | v ( t ) (cid:105) = | v (cid:105) with E ( t ) = t/T · M + (cid:0) − t/T (cid:1) I, | v (0) (cid:105) = | v (cid:105) , and | v ( T ) (cid:105) = | v M − (cid:105) . That is, we can evolve the state | v (cid:105) at time t = 0 to the state | v M − (cid:105) at time t = T . The derivation equation of | v ( t ) (cid:105) is given by E ( t ) ddt | v ( t ) (cid:105) = − G ( t ) | v ( t ) (cid:105) with G ( t ) = (cid:0) M − I (cid:1) /T . It is a special case of a generalised time evolution defined in Eq. (1), with B ( t ) = E ( t ), A j ( t ) = − G ( t ), and | v (cid:48) j ( (cid:126)θ (cid:48) j ( t )) (cid:105) = | v ( (cid:126)θ ( t )) (cid:105) . Therefore, we have the evolution equation as Eq. (3) with coefficients˜ M k,j = Re (cid:32) ∂ (cid:104) v ( (cid:126)θ ( t )) | ∂θ k E † ( t ) E ( t ) ∂ | v ( (cid:126)θ ( t )) (cid:105) ∂θ j (cid:33) , ˜ V k = − Re (cid:32) ∂ (cid:104) v ( (cid:126)θ ( t )) | ∂θ k E † ( t ) G ( t ) | v ( (cid:126)θ ( t )) (cid:105) (cid:33) . M − | v (cid:105) by measuring the coefficients and evolving the parameters accordingly. It is important tonote that ˜ V k can be efficiently computed with the quantum circuit shown in Fig. 4. Therefore, linear equations canbe solved with shallow quantum circuits which necessitate two or three controlled operations from an ancilla qubit.For the second type of algorithm, we assume that M = U DV and the solution is given by | v M − (cid:105) = V † D − U † | v (cid:105) . With similar methods for realising matrix multiplication, we can also realise the matrix inversion operation with realand imaginary time evolutions.
RESOURCE ESTIMATION
In this section, we discuss the complexity of the variational algorithms. We give the resource estimation for achievinga desired simulation accuracy.
Resource estimation for simulating general processes
We first discuss the resource required for implementing the variational algorithm for the generalised time evolution.Suppose we aim to simulate the evolution B ( t ) ddt | v ( t ) (cid:105) = (cid:88) j A j ( t ) | v (cid:48) j ( t ) (cid:105) , (32)from time t = 0 to t = T . We approximate | v ( t ) (cid:105) by | v ( (cid:126)θ ( t )) (cid:105) and evolve the parameters according to (cid:88) j ˜ M k,j ˙ θ j = ˜ V k , (33)with a time step δt . Here the definition of the coefficients can be found in the main text. The error of the algorithmcan be described by ε = D (cid:16) | v ( T ) (cid:105) , | v ( (cid:126)θ ( T )) (cid:105) (cid:17) , (34)where D ( ρ, σ ) = Tr [ | ρ − σ | ] is the trace distance between two states. Denote E ( t − δt, t ) as the evolution from time t − δt to t , then ε = D (cid:16) E ( T − δt, T )( | v ( T − δt ) (cid:105) ) , | v ( (cid:126)θ ( T )) (cid:105) (cid:17) , ≤ D (cid:16) E ( T − δt, T )( | v ( T − δt ) (cid:105) ) , E ( T − δt, T )( | v ( (cid:126)θ ( T − δt )) (cid:105) ) (cid:17) + D (cid:16) E ( T − δt, T )( | v ( (cid:126)θ ( T − δt )) (cid:105) ) , | v ( (cid:126)θ ( T )) (cid:105) (cid:17) , ≤ D (cid:16) | v ( T − δt ) (cid:105) , | v ( (cid:126)θ ( T − δt )) (cid:105) (cid:17) + D (cid:16) E ( T − δt, T )( | v ( (cid:126)θ ( T − δt )) (cid:105) ) , | v ( (cid:126)θ ( T )) (cid:105) (cid:17) . (35)Here the second line is due to the triangle inequality of the distance and the third line assumes that the evolution isdissipative in the sense that it can only decrease the distance. This is true for both real and imaginary time evolutionand the open system simulation. Whether this is true for linear algebra problems depends on the problem. As thiswork mainly focuses on the simulation of open quantum systems, we leave a detailed discussion of the resourceestimation for linear algebra problems in a future work.Following a recursive procedure, we thus arrive at the final upper bound for the error of the simulation algorithm ε ≤ (cid:88) t = δt : δt : T D (cid:16) E ( t − δt, t )( | v ( (cid:126)θ ( t − δt )) (cid:105) ) , | v ( (cid:126)θ ( t )) (cid:105) (cid:17) , (36)where we use D (cid:16) | v (0) (cid:105) , | v ( (cid:126)θ (0)) (cid:105) (cid:17) = 0 by assuming that the ansatz can perfectly represent the initial state at time t = 0. For each term, we can also further divide it into two parts D (cid:16) E ( t − δt, t )( | v ( (cid:126)θ ( t − δt )) (cid:105) ) , | v ( (cid:126)θ ( t )) (cid:105) (cid:17) = δε I + δε A , (37)3with δε I = D (cid:16) | v ( (cid:126)θ ( t )) (cid:105) , | v ( (cid:126)θ ( t )) (cid:105) (cid:17) ,δε A = D (cid:16) E ( t − δt, t )( | v ( (cid:126)θ ( t − δt )) (cid:105) ) , | v ( (cid:126)θ ( t )) (cid:105) (cid:17) . (38)Here | v ( (cid:126)θ ( t )) (cid:105) is the ideal state obtained with exact values of ˜ M and ˜ V . Therefore δε I characterises the implemen-tation error induced from imperfect ˜ M and ˜ V ; δε A characterises the algorithmic error due to finite time step andinsufficient ansatz. Overall, the total error is upper bounded by ε ≤ ε I + ε A = (cid:88) t = δt : δt : T ( δε I + δε A ) , (39)with ε I = (cid:80) t = δt : δt : T δε I and ε A = (cid:80) t = δt : δt : T δε A . In the following, we analyse these two terms separately. Implementation error
The implementation error mainly comes from the imprecise estimation of ˜ M and ˜ V owing to either physical erroror shot noise. Denote the exact value of ˜ M and ˜ V by ˜ M and ˜ V , respectively. Then we have the exact derivative˙ (cid:126)θ = ˜ M − ˜ V and the practically measure derivative ˙ (cid:126)θ = ˜ M − ˜ V . Denote ˜ M = ˜ M + δ ˜ M , ˜ V = ˜ V + δ ˜ V , and ˙ (cid:126)θ = ˙ (cid:126)θ + δ ˙ (cid:126)θ with δ ˜ M , δ ˜ V , and δ ˙ (cid:126)θ representing the small noise perturbation. In Ref. [13], it is shown that the implementationerror can be estimated by δε I = D (cid:16) | v ( (cid:126)θ ( t )) (cid:105) , | v ( (cid:126)θ ( t )) (cid:105) (cid:17) = (cid:113) δ ˙ (cid:126)θ T B δ ˙ (cid:126)θδt + O ( δt ) , (40)where B k,k (cid:48) = ∂ (cid:104) v ( (cid:126)θ ( t )) | ∂θ k ∂ | v ( (cid:126)θ ( t )) (cid:105) ∂θ k (cid:48) − ∂ (cid:104) v ( (cid:126)θ ( t )) | ∂θ k | v ( (cid:126)θ ( t )) (cid:105) (cid:104) v ( (cid:126)θ ( t )) | ∂ | v ( (cid:126)θ ( t )) (cid:105) ∂θ k (cid:48) . (41)Here we assumed that the evolution of Eq. (32) is normalised, i.e., the state vector | v ( (cid:126)θ ( t )) (cid:105) is normalised.The error of the parameters can be also represented as a function of δ ˜ M and δ ˜ V , δ ˙ (cid:126)θ = ˜ M − ˜ V − ˜ M − ˜ V , = (cid:0) ˜ M ( I + ˜ M − δ ˜ M ) (cid:1) − ( ˜ V + δ ˜ V ) − ˜ M − ˜ V = ( I + ˜ M − δ ˜ M ) − ˜ M − ( ˜ V + δ ˜ V ) − ˜ M − ˜ V ≈ ( I − ˜ M − δ ˜ M )( ˜ M − ˜ V + ˜ M − δ ˜ V ) − ˜ M − ˜ V ≈ ˜ M − δ ˜ V − ˜ M − δ ˜ M ˙ (cid:126)θ . (42)Here in the fourth line we consider a Taylor expansion with (cid:107) ˜ M − δ ˜ M (cid:107) ≤ (cid:107) ˜ M − (cid:107)(cid:107) δ ˜ M (cid:107) (cid:28) δε I = (cid:113) δ ˙ (cid:126)θ T B δ ˙ (cid:126)θδt + O ( δt ) (cid:46) (cid:112) (cid:107)B(cid:107)(cid:107) δ(cid:126) ˙ θ (cid:107) δt, (cid:46) (cid:112) (cid:107)B(cid:107) ( (cid:107) ˜ M − (cid:107)(cid:107) δ ˜ V (cid:107) + (cid:107) ˜ M − (cid:107) (cid:107) ˜ V (cid:107)(cid:107) δ ˜ M (cid:107) ) δt. (43)Suppose the dominant error is the shot noise, we can have the relationship between the implementation error andthe number of samples. Given the number of measurements N s for each term of ˜ M and ˜ V , we have (cid:107) δ ˜ M (cid:107) ≈ (cid:113)(cid:80) j,k ˜ m k,j √ N S , (cid:107) δ ˜ V (cid:107) ≈ (cid:112)(cid:80) k ˜ v k √ N S (44)4where ˜ m k,j = | α ( (cid:126)θ ( t )) | (cid:88) i,q | g ∗ k,i g j,q |(cid:107) B ( t ) (cid:107) + (cid:12)(cid:12)(cid:12)(cid:12) ∂α ∗ ( (cid:126)θ ( t )) ∂θ k (cid:12)(cid:12)(cid:12)(cid:12) | α ( (cid:126)θ ( t )) | (cid:88) i | g j,i |(cid:107) B ( t ) (cid:107) + (cid:12)(cid:12)(cid:12)(cid:12) ∂α ∗ ( (cid:126)θ ( t )) ∂θ j (cid:12)(cid:12)(cid:12)(cid:12) | α ( (cid:126)θ ( t )) | (cid:88) i | g k,i |(cid:107) B ( t ) (cid:107) + (cid:12)(cid:12)(cid:12)(cid:12) ∂α ∗ ( (cid:126)θ ( t )) ∂θ k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∂α ∗ ( (cid:126)θ ( t )) ∂θ j (cid:12)(cid:12)(cid:12)(cid:12) (cid:107) B ( t ) (cid:107) ˜ v k = (cid:88) j (cid:20)(cid:12)(cid:12)(cid:12)(cid:12) ∂α ∗ ( (cid:126)θ ( t )) ∂θ k (cid:12)(cid:12)(cid:12)(cid:12) | α ∗ ( (cid:126)θ (cid:48) j ( t )) | (cid:107) A j ( t ) (cid:107) + | α ∗ ( (cid:126)θ ( t )) | | α ∗ ( (cid:126)θ (cid:48) j ( t )) | (cid:88) i | g k,i | (cid:107) A j ( t ) (cid:107) (cid:21) (cid:107) B ( t ) (cid:107) , (45)and g k,i comes from ∂ | ϕ ( (cid:126)θ ( t )) (cid:105) ∂θ k = (cid:88) i g k,i R k,i | ¯0 (cid:105) . (46)Therefore, (cid:107) δ ˙ (cid:126)θ (cid:107) can be upper bounded by (cid:107) δ ˙ (cid:126)θ (cid:107) ≤ ∆ √ N S (47)with ∆ = (cid:107) ˜ M − (cid:107) − (cid:107) ˜ V (cid:107) (cid:115)(cid:88) j,k ˜ m k,j + (cid:107) ˜ M − (cid:107) (cid:115)(cid:88) k ˜ v k . (48)The implementation error at each step is thus upper bounded by δε I ≤ (cid:112) (cid:107)B(cid:107) ∆ √ N S δt, (49)and the total implementation error is ε I ≤ (cid:112) (cid:107)B(cid:107) max ∆ max √ N S T, (50)where (cid:107)B(cid:107) max and ∆ max denote the maximal possible value at different times. Therefore, in order to have animplementation error ε I , we find that we can write N S = (cid:107)B(cid:107) max ∆ T /ε I . Algorithmic error
The algorithmic error mainly comes from finite time step and imperfect ansatz. Following the definition of the timeevolution, we have δε A = D (cid:16) E ( t − δt, t )( | v ( (cid:126)θ ( t − δt )) (cid:105) ) , | v ( (cid:126)θ ( t )) (cid:105) (cid:17) = (cid:112) ∆ δt + ∆ δt + O ( δt ) (cid:46) (cid:113) ∆ (max)2 δt + (cid:113) ∆ (max)3 δtδt, (51)where the first term comes from imperfection of ansatz and the second term comes from finite time step [13]. Here∆ = (cid:104) δv ( (cid:126)θ ) | δv ( (cid:126)θ ) (cid:105) , ∆ = (cid:104) δv ( (cid:126)θ ) | δv ( (cid:126)θ ) (cid:105) + (cid:104) δv ( (cid:126)θ ) | δv ( (cid:126)θ ) (cid:105) , (52)with | δv ( (cid:126)θ ( t )) (cid:105) = | dv ( t ) (cid:105) − (cid:88) j ˙ θ j ∂∂θ j | v ( (cid:126)θ ( t )) (cid:105)| δv ( (cid:126)θ ( t )) (cid:105) = 12 (cid:0) ddt | dv ( t ) (cid:105) − (cid:88) jj (cid:48) ˙ θ j ˙ θ j (cid:48) ∂ ∂θ j ∂θ j (cid:48) | v ( (cid:126)θ ( t )) (cid:105) (cid:1) . (53)5And ∆ (max)2 and ∆ (max)3 are the maximal possible values at time t , respectively.Assuming that the ansatz can always represent the target state with ∆ (max)2 = 0, the total algorithmic error isupper bounded by ε A (cid:46) (cid:113) ∆ (max)3 δtT. (54)To suppress the effect of the second term to ε A , we need to have δt ≈ ε A / (∆ (max)3 T ), hence the number of stepsrequired is N A = T /δt ≈ ∆ (max)3 T /ε A . Resource estimation
Now, we combine the results and show the overall resource cost. At each step, the number of required individualquantum circuits N I is N I = N B † B N P N D + 2 N B N D N P + N B † B + N B N A j N (cid:48) A j ( N P N D + 1) , (55)where the first term and the second terms correspond to the number of individual quantum circuits to populate eachelement of ˜ M and ˜ V , respectively. Here, N P is the number of parameters, N D is the number of terms in Eq. (46), N A j , N B and N B † B are the number of Pauli terms to decompose the A j ( t ), B ( t ) and B † ( t ) B ( t ) matrices in Eq. (32),respectively, and N (cid:48) A j is the number of A j ( t ) matrices (the number of the terms in the right hand side of Eq. (32)).Thus the total number of measurements N tot required to suppress the implementation error by shot noise to ε I andthe effect of algorithmic error to ε A is N tot = N I × N A × N S ≈ (cid:107)B(cid:107) max ∆ ∆ (max)3 T ε A ε I ( N B † B N P N D + 2 N B N D N P + N B † B + N B N A j N (cid:48) A j ( N P N D + 1)) . (56)By taking ε I = ε A = ε/ ε , we have N tot ≈ N I × N A × N S = 16 (cid:107)B(cid:107) max ∆ ∆ (max)3 T ε ( N B † B N P N D + 2 N B N D N P + N B † B + N B N A j N (cid:48) A j ( N P N D + 1)) . (57)Note that this resource estimation is only a pessimistic asymptotic estimation, while the realistic resource requiredcan be much less for a given problem further with the optimisation of measurement schemes Resource estimation for matrix multiplication via singular value decomposition method
In this section, we show the resource required for realising matrix multiplication via the singular value decompositionmethod, which is used in the variational algorithm for simulating open quantum systems. We leave the detailedresource analysis for the other variational algorithms of linear algebra problems in a future work.Suppose the matrix M can be decomposed as M = U DV , where U and V are unitary matrices and D is adiagonal matrix with non-negative entries. Suppose that U , V and D can be decomposed as U = exp( − iH U T U )and V = exp( − iH V T V ) and D ≈ exp( − H D T D ), and these operations can be realised by variational real andimaginary time simulation algorithms. As we only approximate D = (cid:80) j a j | j (cid:105) (cid:104) j | by D ≈ exp( − H D T D ) with − H D T = (cid:80) a j (cid:54) =0 log( a j ) | j (cid:105) (cid:104) j | − α (cid:80) a j =0 | j (cid:105) (cid:104) j | , we first analyse the error introduced from this approximation.6 Accuracy of approximating D Denote D α = D + ∆ α with D = (cid:80) a j (cid:54) = a j | j (cid:105) (cid:104) j | and ∆ α = (cid:80) a j =0 e − α | j (cid:105) (cid:104) j | . We use D α to approximate D , for agiven vector | v (cid:105) , the error of the matrix can propagate to the state as ε D = (cid:13)(cid:13)(cid:13)(cid:13) D α | v (cid:105) (cid:104) v | D α (cid:104) v | D α | v (cid:105) − D | v (cid:105) (cid:104) v | D (cid:104) v | D | v (cid:105) (cid:13)(cid:13)(cid:13)(cid:13) = 1 (cid:104) v | D | v (cid:105) (cid:13)(cid:13)(cid:13)(cid:13) (cid:104) v | D | v (cid:105)(cid:104) v | D α | v (cid:105) D α | v (cid:105) (cid:104) v | D α − D | v (cid:105) (cid:104) v | D (cid:13)(cid:13)(cid:13)(cid:13) , = 1 (cid:104) v | D | v (cid:105) (cid:13)(cid:13)(cid:13)(cid:13) (cid:104) v | D | v (cid:105)(cid:104) v | D α | v (cid:105) ( D + ∆ α ) | v (cid:105) (cid:104) v | ( D + ∆ α ) − D | v (cid:105) (cid:104) v | D (cid:13)(cid:13)(cid:13)(cid:13) , = 1 (cid:104) v | D | v (cid:105) (cid:13)(cid:13)(cid:13)(cid:13)(cid:18) (cid:104) v | D | v (cid:105)(cid:104) v | D α | v (cid:105) − (cid:19) D | v (cid:105) (cid:104) v | D + (cid:104) v | D | v (cid:105)(cid:104) v | D α | v (cid:105) (( D + ∆ α ) | v (cid:105) (cid:104) v | ( D + ∆ α ) − D | v (cid:105) (cid:104) v | D ) (cid:13)(cid:13)(cid:13)(cid:13) , ≤ (cid:12)(cid:12)(cid:12)(cid:12) (cid:104) v | D | v (cid:105)(cid:104) v | D α | v (cid:105) − (cid:12)(cid:12)(cid:12)(cid:12) + 1 (cid:104) v | D | v (cid:105) (cid:107) ( D + ∆ α ) | v (cid:105) (cid:104) v | ( D + ∆ α ) − D | v (cid:105) (cid:104) v | D (cid:107) . (58)Here (cid:107) ρ (cid:107) = Tr | ρ | is the trace norm of ρ . The first term can be bounded as (cid:12)(cid:12)(cid:12)(cid:12) (cid:104) v | D | v (cid:105)(cid:104) v | D α | v (cid:105) − (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) (cid:104) v | D α | v (cid:105) − (cid:104) v | D | v (cid:105)(cid:104) v | D α | v (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) (cid:104) v | ∆ α | v (cid:105)(cid:104) v | D α | v (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) ≤ e − α C , (59)where C = (cid:104) v | D | v (cid:105) determines the norm after multiplying the matrix to the vector, which can be assumed to belower bounded by a constant. As the value of C can be experimentally measured, the cases with an exponential small C can be experimentally identified. In the case with an exponentially small C , it just indicate that the output vectoris a zero vector.For the second term, we can bound it as1 (cid:104) v | D | v (cid:105) (cid:107) ( D + ∆ α ) | v (cid:105) (cid:104) v | ( D + ∆ α ) − D | v (cid:105) (cid:104) v | D (cid:107) ≤ C (cid:107) ∆ α | v (cid:105) (cid:104) v | D + D | v (cid:105) (cid:104) v | ∆ α + ∆ α | v (cid:105) (cid:104) v | ∆ α (cid:107) ≤ C (cid:107) ∆ α | v (cid:105) (cid:104) v | D (cid:107) + (cid:107) D | v (cid:105) (cid:104) v | ∆ α (cid:107) + (cid:107) ∆ α | v (cid:105) (cid:104) v | ∆ α (cid:107) = e − α C , (60)where we used (cid:107) ∆ α | v (cid:105) (cid:104) v | D (cid:107) = (cid:107) D | v (cid:105) (cid:104) v | ∆ α (cid:107) = 0. Therefore, the total error induced from the approximation of D α is ε D = 2 e − α C , (61)and we can choose α = ln(2 /Cε D ) /
2. Therefore, when C ≥ / Poly( n ) with n denoting the number of qubits, we have α = O (cid:18) log (cid:18) ε D (cid:19) + log Poly( n ) (cid:19) . (62)We note that even when C is exponentially small with respect to n as C = 1 /c n with constant c , we can still choose α to be a constant. Nevertheless, we should always measure the value of C at the beginning of the algorithm anddirectly output a zero vector when C is too small. Resource analysis
According to our resource analysis for the generalised time evolution together with the error for approximating D ,the implementation error is ε I (cid:46) (cid:112) (cid:107)B(cid:107) max ∆ max T SV D √ N S , (63)7and the algorithmic error is ε A (cid:46) (cid:113) ∆ (max)3 δtT SV D + ε D , (64)where T SV D = T U + T V + T D . Therefore, the total number of measurements N SV D tot required to suppress the effectof shot noise to ε I and the effect of algorithmic error to ε A is N SV D tot ≈ (cid:107)B(cid:107) max ∆ ∆ (max)3 T SV D ( ε A − ε D ) ε I ( N P N D + N P N H N D ) , (65)where N H is the largest number of terms when H U , H V , and H D is decomposed as a linear combination of Paulioperators. Resource estimation for stochastic Sch¨odinger equation
The resource estimation for the variational algorithm of open quantum systems is more involved. The simulationerror consists of the following parts1. The algorithmic error from approximating the Lindblad master equation with the stochastic Schr¨odinger equa-tion with finite number of trajectory samples.2. For each trajectory, we approximate the continuous evolution with the generalised time evolution and the jumpprocess with the matrix multiplication algorithm.3. For each jump, we need to estimate the jump probability and determining the jump time, whose error can alsocause implementation errors.For the first type of error, the error can be upper bounded to a small value with a sufficiently large number of sampleproportional to a polynomial function of the inverse of the accuracy. For the second type of error, we can bound theerror with the analysis for the generalised time evolution. For the jump process, we can use the analysis for matrixmultiplication for each jump and we show in the following that the number of jumps is proportional to the evolutiontime. For the third type error, we can also bound it with a similar analysis for the generalised time evolution. In thefollowing, we show how to estimate the number of jumps and we leave the detailed analysis of resource estimation ina future work.
Number of jumps