Variational quantum simulation of general processes

Hybrid variational quantum algorithms have been proposed for simulating many-body quantum systems with shallow quantum circuits, and are therefore relevant to Noisy Intermediate Scale Quantum devices. These algorithms are often discussed as a means to solve static energy spectra and simulate the dynamics of real and imaginary time evolutions. Here we consider broader uses of the variational method to simulate general processes. We first show a variational algorithm for simulating the generalised time evolution with a non-Hermitian Hamiltonian, such as parity-time symmetric Hamiltonians. Then we consider linear algebra tasks including solving linear systems of equations and matrix-vector multiplications, vital components in many fields including machine learning and optimisation. We convert the linear algebra tasks into dynamics problems and show how they can be simulated with the algorithm for generalised time evolution. Furthermore, when considering matrices that are products of small matrices, we also propose an alternative variational algorithm that can realise matrix-vector multiplication with a simpler circuit. Finally, we focus on open quantum systems and apply the developed methods to the variational simulation of the stochastic Schr\"odinger equation. We numerically test our theory with a three-qubit 1D Ising model under relaxation.

Hybrid variational quantum algorithms have been proposed for simulating many-body quantum systems with shallow quantum circuits, and are therefore relevant to Noisy Intermediate Scale Quantum devices. These algorithms are often discussed as a means to solve static energy spectra and simulate the dynamics of real and imaginary time evolutions. Here we consider broader uses of the variational method to simulate general processes. We first show a variational algorithm for simulating the generalised time evolution with a non-Hermitian Hamiltonian, such as parity-time symmetric Hamiltonians. Then we consider linear algebra tasks including solving linear systems of equations and matrix-vector multiplications, vital components in many fields including machine learning and optimisation. We convert the linear algebra tasks into dynamics problems and show how they can be simulated with the algorithm for generalised time evolution. Furthermore, when considering matrices that are products of small matrices, we also propose an alternative variational algorithm that can realise matrix-vector multiplication with a simpler circuit. Finally, we focus on open quantum systems and apply the developed methods to the variational simulation of the stochastic Schrödinger equation. We numerically test our theory with a three-qubit 1D Ising model under relaxation.
Introduction.-Variational algorithms have been developed as a powerful classical tool for simulating manybody quantum systems [1][2][3][4][5]. The core idea is based on the intuition that physical states with low energy generally belong to an exponentially small manifold of the whole Hilbert space. As quantum circuits can efficiently prepare states that may not be efficiently represented classically, the variational method has been recently generalised to the quantum regime with trial states efficiently prepared by quantum hardware and information extracted from a coherent measurement of the state . The trial state in variational quantum algorithms can be prepared with shallow quantum circuits [27][28][29][30], which can be robust to a certain amount of device noise and is compatible with near-term Noisy Intermediate Scale Quantum (NISQ) hardware [31]. Variational quantum algorithms can be utilised for efficiently finding energy spectra [7-12, 19, 21-25, 32] and simulating real time Schrödinger evolution [13,33] of classically intractable many-body systems. Although quantum circuits are unitary operations, the variational algorithm is not limited to energy minimisation and unitary processes and it can be used to simulate dissipative imaginary time evolution that cannot be straightforwardly mapped to unitary gates [20,34].
In this work, we study the capability of variational quantum algorithms and show that they are not limited to these problems. First, we propose a variational quantum algorithm for simulating the generalised time evolution defined in Eq. (1) below. Our algorithm can be regarded as a unified framework for simulating general time evolutions, which incorporates the special cases of real and imaginary time evolutions [13,20,34]. This framework also works for non-Hermitian quantum mechanics [35][36][37] that describes nonequilibrium processes [38], parity-time symmetric Hamiltonians [39][40][41], open quan-tum systems [42], etc.
Next, we focus on problems of linear systems of equations and matrix-vector multiplications, which are important in many fields including machine learning and optimisation [43,44]. Various algorithms have been developed for linear systems of equations with universal quantum computers [45][46][47][48][49][50][51][52], which have profound applications in quantum machine learning [53][54][55][56][57]. However they generally require a long depth circuit that relies on a fault tolerant quantum computer. In this work, we introduce variational quantum algorithms for these two tasks. For general sparse matrices, we show how solutions of these two tasks can be converted into a generalised time evolution, which can be variationally simulated. For special matrices that are products of small matrices that only involve a constant number of qubits, the solutions can be even more easily obtained only with variational real and imaginary time evolution.
Finally, we combine the developed variational algorithms to simulate the evolution of open quantum systems [58][59][60]. Under the description of the stochastic Schrödinger equation, the evolution of open quantum systems can be regarded as an average of wave functions that undergo a continuous measurement induced from the environment [60,61]. The evolution of each wave function is composed of two processes that can be both simulated with variational algorithms: It may continuously evolve under the generalised time evolution with the system Hamiltonian and the damping effect due to continuous measurement; alternatively the state discontinuously jumps according to the measurement results. The continuous process can be described by the generalised time evolution, and the jump process is a matrix-vector multiplication process. Simulating the evolution of general open quantum systems is of great importance for understanding any quantum system that interacts with an environment. Existing quantum algorithms [62][63][64][65][66][67] for simulating open quantum systems generally require deep quantum circuits. As our algorithm is compatible with shallow circuits, it can be realised with NISQ hardwares.
Generalised time evolution.-We first introduce variational quantum simulation of generalised time evolution, with |dv(t) = j A j (t) |v ′ j (t) . Here, A j (t) is a time dependent general sparse (non-Hermitian) operators, |v(t) is the system state, and each of |v ′ j (t) can be either |v(t) or any known state that can be efficiently prepared with a quantum circuit. The states |v(t) and |v ′ j (t) can be (un)normalised states with normalisation factors α(t) and α ′ j (t), respectively. In practice, we assume that A j (t) can be decomposed as a linear combination of Pauli operators A j (t) = i λ i (t)σ i with complex coefficients λ i and a polynomial (to the system size) number of tensor products of Pauli matrices σ i .
In variational quantum simulation, instead of directly simulating the dynamics, we assume that the state can be represented by parameterised quantum states |v( θ(t)) = α( θ 0 (t)) |ϕ( θ 1 (t)) with θ := ( θ 0 , θ 1 ). Then we project the original evolution to the evolution of the parameters via McLachlan's principle [68], where |ψ = ψ|ψ . By minimising the distance between the true evolution and the evolution of the parameterised states, we find the equation of parameters as jM k,jθj =Ṽ k .
Here the coefficients are linear sums of state overlaps, where each term can be efficiently measured with quantum circuits. We specify the detailed derivation, expression of the coefficients and the quantum circuits in Supplementary Materials. We consider two examples: real and imaginary time evolution [13,20,34] [13,20,34].
Suppose |v(t) = |ϕ( θ(t)) , is the same for both real and imaginary time evolution; andṼ is Im ∂ ϕ( θ(t))| ∂θi H |ϕ( θ(t)) and for real and imaginary time evolution, respectively. To measureM orṼ , we generally need an extra ancillary qubit with controlled operations applied between the ancilla and the system state. Interestingly, it was recently shown that the ancilla and the constant number of controlled operations can be reduced [69].
Variational algorithms for linear algebra.-Next we focus on solving linear systems of equations and matrixvector multiplications. For a sparse square matrix M and a state vector |v 0 , we aim to find for these two approximations respectively. We first focus on matrix-vector multiplication by introducing a dynamical process that evolves the initial vector |v 0 to the target state |v M . The natural evolution path is via a linear extrapolation between |v 0 and |v M as |v(t) = E(t) |v 0 with E(t) = t/T · M + 1 − t/T I, |v(0) = |v 0 , and |v(T ) = |v M . Different evolution paths can be also considered. For example, in the conventional Hamiltonian simulation scenario, we have M = e −iHT and it corresponds to an exponential extrapolation. We also consider linear extrapolation between normalised states in Supplementary Materials and we leave the discussion of other possible evolution paths to future works. Given the evolution via linear extrapolation, the time (1), which can be variationally simulated.
We now introduce a second method for realising matrix-vector multiplication with only real and imaginary time evolution. This method assumes an efficient singular value decomposition of M as M = U DV , with unitary matrices U , V and diagonal matrix D with non-negative entries. Suppose the unitary matrices U and V can be represented by U = exp(−iH U T U ) and V = exp(−iH V T V ) with time T U and T V , respectively, then the multiplication of U and V can be realised by evolving the state with Hamiltonian H U and H V via variational real time simulation. Given a spectral decom- properly chosen large constant α satisfying α ≫ log(a j ). Therefore, we can define an unnormalised imaginary time evolution d|v(τ ) dτ = −H D |v(τ ) , so that the initial vector |v 0 is evolved to D |v 0 from τ = 0 to τ = T . This second method is easier to implement; it requires the singular value decomposition of M and evolution operator for the unitary and diagonal matrices. Although these assumptions are restrictive, they are natural when considering matrices M that are products of matrices that only involve a few qubits, i.e., M = M 1 ⊗ · · · ⊗ M L with M i acting on a small constant number of qubits. We will shortly show that this variational algorithm is particularly useful for simulating the jump process of the stochastic Schrödinger equation.
We now discuss the solution of the linear equation M |v M −1 = |v 0 with invertible matrix M. We consider the extrapolation between |v 0 and |v M −1 as Although this is slightly different from the generalised time evolution of Eq. (1), we can still simulate it with the variational method by assuming |v(t) = |ϕ( θ(t)) . With McLachlan's variational principle [68], we have the evolution equation as Eq. (3) with . Therefore we can calculate M −1 |v 0 by measuring the coefficients and evolving the parameters accordingly. We refer the reader to Supplementary Materials for the detailed derivation and the quantum circuits to implement this algorithm.
Open system simulation.-Now, we apply the developed variational algorithms to simulate the stochastic mater equation of open quantum systems. Note that this approach differs from that in Ref. [34] where the evolution of the Lindblad master equation is directly simulated using a method that requires two copies of the states. Our scheme here only require one copy of the state. Suppose the dynamics of open quantum systems is described by the Lindblad master equation, Here, the system Hamiltonian is H and the interaction with the environment Lρ is described by Lρ = , with Lindblad operators L i . Although the Lindblad master equation directly evolves the density matrix, it is equivalent to the stochastic Schrödinger equation by averaging the trajectory of each pure state evolved under continuous measurements [60,61]. Because the measurement process is stochastic, the wave function has a stochastic evolution.
Given this Lindblad master equation, the stochastic Schrödinger equation for each single trajectory |ψ c (t) is where d |ψ c (t) = |ψ c (t + dt) − |ψ c (t) , and dN k randomly takes either 0 or 1 which satisfies dN k dN k ′ = δ kk ′ dN k and E[dN k ] = ψ c (t)| L † k L k |ψ c dt. At each time t, we can assume a positive obervable valued measure- For measurement outcome O k , the state discontinuously jumps to L k |ψ c (t) / L k |ψ c (t) with probability E[dN k ]. And the total jump probability is γ(t) = k E[dN k ]. For outcome O 0 with probability 1 − γ(t), we have dN k = 0, ∀k and the state evolves under the generalised time evolution with operator Here, −iH corresponds to the conventional real time Schrödinger evolution with Hamiltonian H and the other terms can be understood as a normalised damping process. Therefore, the whole process is composed of two parts: the continuous process governed by the first term and the quantum jump process described by the second term. Now, we show how to simulate stochastic Schrödinger equation using the Monto Carlo method. Suppose the state jumps at time t, then at time t + τ , the probability p(t + τ ) that the state does not jump is with Γ(t + τ ) = t+τ t γ(t ′ )dt ′ . When jump happens at time t, a uniform random number q ∈ [0, 1] is generated. Then the time of the next jump is determined by accumulating time τ until we have p(t + τ ) = q. When jump happens, a random number q ′ ∈ [0, 1] is generated to determine which jump operator to apply at each timestep. The state is updated to L k |ψ and N L is the number of the Lindblad operators. Considering discretised time with initial state |ψ c (0) , the stochastic Schrödinger equation from time 0 to T can be simulated as follows. Evolve the state |ψc(t) under A in Eq. (7).

12:
Reset Γ = 0 and randomly generate q ∈ [0, 1]. Now we show how to simulate the stochastic Schrödinger equation, Algorithm 1, with the variational algorithms developed in this work. Suppose the state |ψ c (t) at time t can be represented by the parametrised state |φ c ( θ(t)) prepared by a quantum 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 according to Eq.
∂θ k , and The values of γ(t) andγ k (t) at step 5 and 8, respectively, can be efficiently measured. The jump at step 11 can be realised by the variational algorithms for matrix-vector multiplication. Especially, when considering L k as a product of operators of each qubits, it can be efficiently realised also with the singular value decomposition method. In practice, we consider sparse Hamiltonian and Lindblad operators L k , therefore all the measurements can be efficiently evaluated. We refer the reader to Supplementary Materials for a detailed resource estimation of our algorithm.
Numerical simulation.-Now we show an example of variational quantum simulation of a three-qubit dissipative 1D Ising model with a transverse field and open boundary, discussed in Ref. [70][71][72][73]. With the Lindblad master equation as in Eq. (5), the Hamiltonian is 1| j being the lowering operator acting on the i th spin. In our simulation, we set J = 1, h X = 1, h Z = 0, and γ = 1. The initial state is |ϕ(0) = |0 1 |0 2 |0 3 , and we simulate the evolution from t = 0 to t = 10. In our simulation, we make use of the Hamiltonian ansätz [74] as shown in Fig. 1 with ten parameters.
Hamiltonian ansätz in our numerical simulation. The single qubit gate is defined by RX (θi) = e −iθ i X and the two qubit gate is RZZ (θi) = e −iθ i Z⊗Z . The last gate U (θ6, . . . , θ10) is a repetition of the first five gates with five different parameters. In total, there are ten parameters.
To simulate quantum jumps induced by Lindblad operators, we use the singular value decomposition method. We decompose the jump operator σ − as |0 1| = |0 0| X. Suppose |0 1| = U DV , we can see that U = I, D = |0 0|, V = X. To realise the V operator, we set H V = X and T V = π/2 such that X = exp(−iH V T V ). In Fig. 2, we show the simulation results and compare it with the exact solution. We measure the average value of Z 1 and we can see that the variational simulation result agrees well with the exact solution with a deviation less than 10 −2 .
Discussion.-To summarise, we extend the variational quantum simulation method to general processes, including the generalised time evolution, matrix-vector multiplication, and the evolution of open quantum systems. Our algorithm for simulating the generalised time evolution can be applied to simulate non-Hermitian quantum mechanics [35][36][37] including nonequilibrium processes [38] and parity-time symmetric Hamiltonians [39][40][41]. Especially, it is shown in Ref. [41] that a quantum state can evolve to the target state faster with non-Hermitian parity-time symmetric Hamiltonians than the case with Hermitian Hamiltonians. Therefore, our variational algorithm for simulating the generalised time evolution may also be useful for designing faster quantum computing algorithms. Meanwhile, the proposed algorithms are compatible with NISQ hardwares and can be further combined with the recently proposed quantum error mitigation techniques [13,17,[75][76][77][78][79][80].

ACKNOWLEDGEMENTS.
This work is supported by the EPSRC National Quantum Technology Hub in Networked Quantum Information Technology (EP/M013243/1). SCB acknowledges support from the EU AQTION project. SE is supported by Japan Student Services Organization (JASSO) Student Exchange Support Program (Graduate Scholarship for Degree Seeking Students). YL is supported by NSAF (Grant No. U1730449).

SUPPLEMENTARY MATERIALS: VARIATIONAL QUANTUM SIMULATION OF GENERAL PROCESSES VARIATIONAL QUANTUM SIMULATION OF REAL AND IMAGINARY TIME EVOLUTION
We first review the variational quantum algorithms for simulating real and imaginary time evolution introduced in Ref. [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 various variational principles -the Dirac and Frenkel variational principle, the McLachlan's variational principle, and the time-dependent variational principle.
The real time evolution is described by the Schrödinger equation, with Hermitian Hamiltonian H. Instead of directly simulating the real time dynamics with the Hamiltonian simulation algorithms [81][82][83][84][85], the variational method assumes that the quantum state |ψ(t) is prepared by a parametrised quantum circuit, |ϕ( θ(t)) = R N (θ N ) . . . R k (θ k ) . . . R 1 (θ 1 ) |0 with each gate R k (θ k ) controlled by the real parameter θ k and the reference state |0 . Here, we denote θ = (θ 1 , θ 2 , . . . , θ N ). According to McLachlan's variational principle [68], the real time dynamics of |ψ(t) can be mapped to the evolution of the parameters θ(t) by minimising the distance between the ideal evolution and the evolution induced of the parametrised trial state, where |ϕ = ϕ|ϕ is the norm of |ϕ . The solution is with coefficients For imaginary time evolution, the normalised Wick-rotated Shrödinger equation is obtained by replacing t = iτ in Eq. (10), Applying a similar procedure for real time evolution, the imaginary time evolution is mapped to the evolution of the parameters via McLachlan's principle, The evolution of the parameters is with M given in Eq. (13) and C defined by The M , V , and C terms can be efficiently measured with quantum circuits. Considering gate based circuits, the derivative of the each parameterised gate can be expressed as state |ψ(t) as with and a normalisation factor The normalisation factor N (t) can be measured from the expectation values of M † + M and M † M for |ψ 0 . Given the definition of the state |ψ(t) at time t, the corresponding derivative equation is Such an equation is also described by the generalised time evolution equation with |v(t) = |ψ(t) , A 1 (t) =Ṅ (t) N (t) I, A 2 (t) = N (t)G, |v ′ 1 (t) = |ψ(t) and |v ′ 2 (t) = |ψ(0) .

DERIVATION OF THE TIME DERIVATIVE EQUATION OF PARAMETERS FOR SOLVING LINEAR EQUATIONS
Here, we derive the time derivative equation of parameters for solving linear equations. Firstly, we consider the extrapolation between |v(0) = |v 0 and |v M −1 = M −1 |v 0 as E(t) |v(t) = |v 0 , where We differentiate both sides with t, to obtain Now, we parametrise the state |v( θ(t)) . According to McLachlan's principle, we have = δ E(t) jθ j ∂ ∂θ j |v( θ(t)) + G(t) |v( θ(t)) = 0, which is equivalent to Then we obtain It is important to note thatṼ k in Eq. (41) can be efficiently computed with the quantum circuit shown in Fig. 3. M k,j can written asM where we set E † (t)E(t) = l β l σ l , and σ l is a Pauli operator. The quantum circuit to compute this value is shown in Fig. 5

RESOURCE ESTIMATION FOR SIMULATING STOCHASTIC SCHRÖDINGER EQUATION
One trajectory of the Stochastic Schrödinger equation is composed of the continuous evolution and jumps processes. The resource cost of the continuous evolution is similar to the one for real time evolution discussed in Ref. [13], which is shown to be polynomial to the evolution time and system size. For the jump processes, as each jump is simulated with real and imaginary time evolution, the resource cost of each jump is therefore also polynomial to the evolution time and system size.
Then we discuss how many jump process occurs on average in the simulation of the Stochastic Schrödinger equation. The averaged number of jump events during time from t to t + dt is Therefore, the average number of jump events from t = 0 to t = T where L ∞ = max eig(L) is the operator norm. For physical systems, we generally have L † k L k ∞ and the number of Lindblad terms equal O(Poly(n)), where n is the system size and Poly(n) is a polynomial function of n. Therefore, the averaged number of jumps is The number of jumps is much fewer when considering the case where each Lindblad operator only locally acts on a constant subsystem. That is, we assume that L † k L k has orthogonal support to each other and L † k L k = O(1). In this case, we have k L † k L k ∞ = O(1) and hence To simulate the Stochastic Schrödinger equation, we also need to sample different random trajectories. When we measure an observable O and hope to suppress the sampling error to ǫ = 1/ √ M , we need M random samples. Therefore, the overall cost should also be multiplied by M . However, it is worth noting that every trajectory is exactly parallel so the overhead M can be also constantly reduced.