Solving Quantum Master Equations with Deep Quantum Neural Networks

Deep quantum neural networks may provide a promising way to achieve quantum learning advantage with noisy intermediate scale quantum devices. Here, we use deep quantum feedforward neural networks capable of universal quantum computation to represent the mixed states for open quantum many-body systems and introduce a variational method with quantum derivatives to solve the master equation for dynamics and stationary states. Owning to the special structure of the quantum networks, this approach enjoys a number of notable features, including the absence of barren plateaus, efficient quantum analogue of the backpropagation algorithm, resource-saving reuse of hidden qubits, general applicability independent of dimensionality and entanglement properties, as well as the convenient implementation of symmetries. As proof-of-principle demonstrations, we apply this approach to both one-dimensional transverse field Ising and two-dimensional $J_1-J_2$ models with dissipation, and show that it can efficiently capture their dynamics and stationary states with a desired accuracy.

Every quantum system is inevitably coupled to its surrounding environment.In most cases, the coupling is of Markovian type and the dynamics of open quantum systems are governed by the Lindblad quantum master equation [23].Thus, solving this master equation plays a key role in studying open quantum systems.However, for quantum many-body systems this turns out to be a formidable challenge due to the exponential growing of the Hilbert space dimension with the size of the system.To combat this challenge, a number of prominent numerical approaches based on classical computers have been developed, such as these using tensor network representations [24][25][26][27][28][29][30][31] or quantum Monte Carlo methods [32,33].More recently, machine-learning inspired approaches based on artificial neural networks [34], especially the restricted Boltzmann machine (RBM), have also been introduced to solve the quantum master equation [35][36][37][38].Each of these approaches bears its own advantages and disadvantages, and the choice of which one to use is problem-specific.For instance, tensor-network approaches are typically very effective for one-dimensional (1D) systems involving small entangle-FIG.1.Schematic illustration of the DQFNN (deep quantum feedforward neural network) method in solving quantum master equations.This DQFNN has d layers in total, with an input (top), output (bottom), and d − 2 hidden layers.A quantum perceptron is defined as an arbitrary unitary operator acting on qubits from neighboring layers [22], and the DQFNN is parametrized by the product of these unitary operators.We use the output state as the variational ansatz state of a quantum open system.The qubits in the output layer are measured on different basis and the measurement results are used by a classical-quantum hybrid optimization algorithm to update the network parameters.
ment, but may face pronounced difficulties in higher dimensions or in situations where massive entanglement comes into play [39,40].Quantum Monte Carlo methods rely on efficient sampling of physical configurations [41] and are effective for certain open quantum systems [32,33], yet could suffer from a severe sign problem that is common in simulating dissipative dynamics [42].The neural-network approaches [35][36][37][38] are generally applicable to high dimensional systems and entanglement is not a limiting factor [43].However, like in most machine-learning tasks [44], their performance depends crucially on the suitable tuning of hyperparameters.In fact, the capability and limitation of neural-network methods in solv-ing quantum many-body problems still remain largely unclear and related studies are at the research forefront.
Here, inspired by the classical neural-network approaches and the exciting experimental progress in developing NISQ devices, we propose a variational DQFNN method to solve the quantum master equations.We note that variational algorithms running on NISQ devices for open quantum systems have been discussed in the literature [45][46][47][48], but to the best of our knowledge most existing works focus on straightforward variational quantum circuits that may suffer from the notorious barren plateau (i.e., vanishing gradient) problem [49,50].We utilize a deep quantum neural network, which is arranged in layers and convenient for the quantum analogue of the backpropagation algorithm, to serve as the ansatz density state of the open quantum system.We adopt a hybrid quantum-classical stochastic reconfiguration (SR) algorithm [51] to variationally solve the quantum master equation in the Lindblad form to obtain the dynamics and steady state.As a result of the special structure of the quantum networks, our approach escape barren plateaus [22] and allows for convenient implementation of translational symmetry and resourcesaving reuse of hidden qubits.In addition, it works for generic open quantum many-body systems, independent of dimensionality, the amount of entanglement involved, and the specific forms of interaction and dissipation.We benchmark our approach with both dissipative 1D transverse field Ising and 2D J 1 − J 2 models.Our results pave a way to explore the rich physics of open quantum systems with quantum neural networks that would be implemented using NISQ devices in the near future.
The general problem and DQFNN structure.-Weconsider the problem of solving the following master equation: where ρ is the density matrix of the system, H is the Hamiltonian governing the unitary part of the dynamics, F i are the jump operators describing the dissipative processes induced by the environment with dissipative rate denoted by γ i , the curly bracket represents the anticommutator, and L is the Liouvillian superoperator [52].
We use the final output state of a quantum neural network, the recently proposed DQFNN in particular [22], to serve as the variational ansatz for the density state of the open quantum system.The network structure is illustrated in Fig. 1.Each node represents a qubit and a quantum perceptron is an arbitrary unitary operator acting on several qubits from neighboring layers.Information propagates from the top (input) layer to the bottom (output) layer.At the i-th layer, the output state ρ o is determined by where ρ is the output state of the (i − 1)-th layer (which serves as the input state for the i-th layer), |ψ (i) ψ (i) | denotes the initial state of the i-th layer, U = ⊗ M j=0 U i,j represents the total unitary operation between the two layers with U i,j denoting the j-th perceptron at the i-th layer, and Tr i−1 denotes a partial trace of the (i − 1)-th layer.The Eq. ( 2) may also be regarded as a mapping M from ρ o , and the whole network features a mapping-composition structure, which is essential for the backpropagation algorithm [44]: Here, |ψ (1) ψ (1) | denotes the initial state of the first (input) layer and we suppose that the network has d layer in total.We mention that this quantum network is universal in the sense that it can map any input density state to an arbitrary output state due to its vast flexibility in designing the qubit perceptrons.In fact, it has been shown in Ref. [22] that such a DQFNN can carry out universal quantum computation even for perceptrons as simple as unitaries acting on two-input and one-output qubits.Here, we focus on two-input one-output qubit perceptrons for simplicity and experimental feasibility.
The general recipe.-Asmentioned, we use ρ to serve as the variational ansatz state and adapt a hybrid quantumclassical SR algorithm [51] to variationally solve Eq. ( 1) to obtain the dynamics and steady state of the open quantum system.More specifically, we solve the following optimization problem variationally: where Θ collectively denotes all the parameters used to describe the qubit perceptrons (unitary gates) and the Fubini-Study (FS) norm is used to measure the distance.For convenience, we omit the labels and rewrite the density matrix ρ o as a vector ρ.Similar to the SR algorithm in solving the quantum master equation based on the variational RBM representation [35,37], the above minimization problem results in the following system of equations: where Θ ν denotes the ν-th parameter, and We initialize the parameters to be small random values and update them iteratively according to the the following rule: where λ is a hyperparameter (called the learning rate in the machine learning literature) introduced to modulate the convergence of the iteration process.
In real experiments, both S and f might be obtained from measurements of observables or their linear combinations (see the Supplemental Material [53] for details).Then, by using Eq. ( 8) we update the parameters of the quantum neural network iteratively and obtain the dynamics of the open system , and σ x 1 σ x 2 are calculated by both the DQFNN and exact diagonalization (ED) method.Here, we set J = 1 as the energy unit and the dissipation ratio is γ = 1.For the DQFNN approach, the learning rate λ is chosen to be λ = 0.01 × 0.999 Ns , where Ns is the number of steps.consequently.At long time limit, the neural network state converges to the steady state of the system.As discussed in Ref. [22], a striking feature of this DQFNN approach is that the parameter matrices can also be calculated layer-by-layer and at any given time we need only access to two layers.This will greatly reduce the memory requirements of the algorithm and enables a resource-saving reuse of the qubits.
Numerical simulations.-Weapply the introduced DQFNN approach to two concrete models to benchmark how it works.The first example involves the dissipative transverse field Ising model in one dimension: where H = J j σ z j σ z j+1 + h j σ x j , γ denotes the dissipation rate, and σ ± j = 1 2 (σ x j ± iσ y j ), with σ x,y,z j being the Pauli matrices for the site j.The first term in the Hamiltonian represents the z-component spin-spin interaction in the longitudinal direction with strength J and the second term denotes a local uniform magnetic field with strength h along the transverse direction.We mention that the steady state of this model has been studied in Ref. [36] with a classical neural-network (i.e., RBM) approach.Here, instead we apply the DQFNN method to obtain both the dynamics and steady state.The quantum neural-network structure is shown in Fig. 2 (a).The system size is N = 5 and the quantum network has five layers and contains 15 qubits in total.Each perceptron consists of three qubits coupled by a series of unitary gates as shown by the quantum circuit in the dotted box in Fig. 2 (a).We consider the periodic boundary condition and implement the translational symmetry for reducing the number of parameters.In addition, owing to the special structures of the quantum perceptrons, the derivatives of the density states with respect to the parameter Θ ν can be obtained by [54]: We mention that the derivatives in Eq. ( 10) is exact, in sharp contrast to the finite-difference method that may inevitably induce a small differential error.This might be important for the convergence of the optimization procedure [55].
With classical computers, we numerically simulate the quantum neural networks and the whole process of using the DQFNN approach to obtain the dynamics and steady state for the dissipative Ising model defined in Eq. ( 9).We start with a DQFNN state where all the constituted qubits are initially set to be |+ = (|0 + |1 )/ √ 2 and all the involved parameters are initially set to be a small value near zero.Then the parameters are updated according to the optimization procedures shown in Eq. (8).Partial of our results are shown in Fig. 2(b) and Fig. 3.In Fig. 2(b), we compute the steady state with varying strength of the transverse field h through the DQFNN approach.We plot the magnetization values σ x,z and correlations σ x 1 σ x 2 for the stationary state, and compare them with the ED results.It is clear that the DQFNN results match the ED results excellently.In Fig. 3, we consider the dynamics of our DQFNN during the optimization process with fixed h = 0.6.Here, we add a random Gaussian noise into the quan- tum derivatives to account for the imperfections in real experiments and to verify the robustness of our DQFNN approach.In Fig. 3(a,b,c), we plot the magnetization values σ x,z and correlations σ x 1 σ x 2 with different noise strength, and compare their corresponding results obtained by ED.From these figures, it is evident that the noiseless DQFNN results match the ED results (with relative error smaller than 10 −2 ) and converge to the corresponding values for stationary state at long time.The added noise does cause small deviation during the dynamic process, while the evolution process eventually converges to the steady state as well.In Fig. 3(d), we show δL, which is the expectation value of the Lindblad operator to measure the convergence to the stationary state, as a function of time.We note that δL approaches zero as we increase the iterations steps, implying that the density state of the system represented by the quantum neural network converges to the stationary state indeed, despite the addition of strong noises.
Our second example concerns the dissipative J 1 −J 2 model defined on a square-lattice in 2D with Hamiltonian: where the first (second) term represents the nearest-neighbor (next-nearest-neighbor) z-component spin-spin interaction with strength J 1 (J 2 ).The dissipation considered here is the same as that for the 1D Ising case [see Eq. ( 9)].For this Hamiltonian, it is easy to observe that there is a geometric frustration to the antiferromagnetic state due to the competi-tion between the first and second terms.For simplicity, we set J 1 = 1 as the energy unit, J 2 = 1 2 J 1 , h = J 1 , and the dissipation rate γ = J 1 .Our numerically simulated results from the DQFNN approach are shown in Fig. 4. In Fig. 4(a,b,c), we plot the magnetization values σ x,z and the correlation function σ x 1 σ x 2 as a function of iteration steps.As shown, it is clear that all these quantities converges nicely to their corresponding exact values for the steady state (after 1000 iterations, their relative errors are all smaller than 10 −2 ).This indicates that after around 1000 iterations, the output state ρ of the DQFNN indeed converges to the steady state, which is also clearly manifested in Fig. 4 as δL converges to zero (after 1000 iterations, |δL| becomes less than 10 −2 ).
We remark that the accuracy of our DQFNN results can be improved in a number of ways, including increasing the number of layers of the quantum neural network or the sampling size at each iteration step, tuning hyperparameters, and designing more appropriate quantum perceptrons, etc. Due to limited classical computation resources, in this paper we only carry out simulations for small systems.However, this will not be a problem if quantum devices are used in practice.As discussed in the Supplementary Material [53], the time complexity of the DQFNN approach is about O(N 3 ), which indicates its scalability to larger systerms.We also stress the difference between the RBM [35][36][37][38] and DQFNN approaches: the RBM method is entirely classical and the gradients that are crucial in updating parameters rely on efficient Monte Carlo sampling; in contrast, for the DQFNN method the gradients can be obtained directly from measurements of observables in experiment [53].This might lead to an advantage in computational cost for the DQFNN approach for certain problems where no efficient sampling scheme is available.In the future, it is of both fundamental and practical interest to carry out a complete study on the advantages and limitations for the DQFNN approach.
Conclusion.-We have introduced a deep quantum neuralnetwork method to solve Lindblad master equations for open quantum many-body systems.Through concrete examples in both 1D and 2D, our results demonstrated that both the dynamics and stationary states of such systems can be efficiently obtained via this method with a desirable accuracy.Due to the special structures of DQFNNs, our approach is generally applicable to high dimensional systems and is independent of the amount of entanglement involved.In addition, it is robust to experimental imperfections and free from the barren plateau problem, and allows a convenient implementation of translational symmetry and a resource-saving reuse of qubits.This classical-quantum hybrid approach is particular suitable for running on NISQ devices in the near future.
Supplementary Material for: Solving Quantum Master Equations with Deep Quantum Neural Networks

(S5)
To sample the S and f efficiently, we generate a Markov chain of the right and left state configuration from step zero ( r 0 , l 0 ) to step s ( r s , l s ).In each step, a new trial right and left state configurations are randomly generated, and the new configuration is accepted with a probability min(1, | ρ lnew, rnew ρ l old , r old | 2 ).We improve our sampling efficiency by generating each new configuration with slightly modifying the strategies discussed in [35,37] In our numerical simulations, the probabilities of the first three types of moves are about 30 times larger than the last two ones.To approximate the expectation value of a general hermitian observable, we consider the accepted probability as min(1, and generate a list of state configuration l 0 → • • • → l s by following possible moves: 1.One lattice site is flipped.2. Two neighboring lattice sites are flipped.
3. All the lattice sites are flipped.4. A new random configuration is generated.
In our numerical simulations, the probabilities of the first two types of moves are about 30 times larger than the last two types of moves.

II. QUANTUM UPDATING SCHEME
With quantum devices, S and f can be obtained efficiently from linear combinations of measurements of proper observables in the experiment.A straightforward way is to do tomography of the output density matrix, which, however, is exponentially expensive as the system size increases.Thus, this method is limited to small systems.To overcome this difficulty, we use the following methods to obtain S and f efficiently.We rewrite S and f as follows: Thanks to the special structures of DQFNNs studied in this paper, the derivative ∂ρ ∂Θµ of ρ with respect to the parameter Θ µ can be written as Eq.( 10) in the main text.For the convenience, we write the l l| • • • | l as a trace operation and the Eqs.(S6-S7) then reduce to: As we can see in Eqs.(S8-S9), each term in f µ and S µ,ν only has two general form, namely, Tr[ρ(Θ)ρ(Θ )] and Tr[ρ(Θ)Lρ(Θ )], where Θ and Θ are two parameter sets that can be different or identical.
We adopt the quantum circuit in Fig. S1(a) to measure Tr[ρ(Θ)ρ(Θ )], which has been discussed in [57].We firstly generate the output density matrices ρ(Θ) and ρ(Θ ) from two DQFNNs, and introduce another ancillary qubit which is initialized to (|0 + |1 )/ √ 2 state.We decompose the spectrum of the whole system: r, l , r ρ l, r ρ l , r | l, l r, r |.After that, we apply a controlled-swap gate, which is defined as on the whole system Finally, we apply another Hardmard gate on the ancilla qubit and measure the spin probability in state |0 : l, r, l , r , m, n ρ l, r ρ l , r m, n| l, l r , r| m, n + 1 4 l, r, l , r , m, n ρ l, r ρ l , r m, n| l , l r, r | m, n l, r, l , r , m, n ρ l, r ρ l , r δ m, l δ n, l δ m, r δ n, r + 1 4 l, r, l , r , m, n ρ l, r ρ l , r δ m, l δ n, l δ m, r δ n, r where P (0) is the probability that the ancillary qubit is in state |0 , so that we can obtain Tr [ρ(Θ)ρ(Θ )] via measure the ancillary state population.
To measuring Tr[ρ(Θ)Lρ(Θ )], we expand the Lindblad equation into a explicit form: Taking into account the one dimensional case in main text Eq.( 9), the Hamiltonian H only involves the local Pauli operators, so that we only need to consider Tr [ρ(Θ)σ x i ρ(Θ )] and Tr ρ(Θ)σ z i σ z i+1 ρ(Θ ) .Here, we write the local Pauli operators in a uniform form Ô = σ x i , σ z i σ z i+1 .Utilizing the quantum circuits in Fig. S1 and obtain Tr ρ(Θ) To calculate the term Tr ρ(Θ)σ These two quantites can be measured by using the method above.We expand the last term σ and measure the real part and imaginary part separately.To calculate Tr [ρ(Θ)σ x,y i ρ(Θ )σ x,y i ], as shown in Fig. S1(c), we apply two controlled-Ô1,2 gates, where Ô1,2 = σ x,y i , and apply the controlled-swap gate and (S20) Then, similar to our previous discussion, we measure the ancillary qubit where w i is the number of perceptrons between i-th layer and ((i + 1))-th layer.II.4.Partial trace the i-th layer, namely, ρ i+1 o = Tr i (ρ i+1 ) and set i = i + 1. II.5.Reuse the qubits.If i < d, the partial traced qubits in i-th layer is re-initialized as (|0 + |1 )/ √ 2 and they can be used in the rest of the network.II.6.Repeat II.2-II.4 until i = d.III Update parameters: III.1.Compute S and f .As the discussion in section I and section II, we can obtain the S and f via a stochastic Markov-chain sampling or directly measure these two quantities.III.2.Update parameters.According to the system equation in the main text, the parameters can be updated as Θ new µ = Θ old µ + λ ν (S −1 ) µ,ν f ν .Here we choose λ to be small enough to monitor the convergence of the optimization process.IV Repeat step II and III until δL reaches minimal and the process get converged.

IV. TOTAL NUMBER OF PARAMETERS AND SYMMETRY IMPLEMENTATION SCHEME
In our DQFNN, the j-th perceptron in i-th layer contains three qubits and a series of coupling unitary operation belongs to SU (2 3 ) group with the most universal form is: where Θ i,j, α denote the parameters for the j-th perceptron in the i-th layer.The index α 1,2,3 takes 0, x, y, z to indicate the pauli matrix I, σ x , σ y , σ z respectly.Hence, to fully describe a quantum perceptron in this case needs 64 parameters at most.As mentioned in the main text, we use an experimentally more practical quantum circuit antasz, which contains 12 Euler rotations and 6 CNOT gates, to implement a quantum perceptron with only 36 parameters.For a fully connected DQFNN, the perceptrons are composed by any possible combination of two qubits from the i-th layer and one qubit from the (i + 1)-th layer, where the total number of perceptrons between two neighbor layers is C 2 ni n i+1 .Hence, the total number of parameters of the DQFNN is: where d is the total number of layers, n i denotes the number of qubits in the i-th layer, N p is the number of parameters for each perceptron, and N All denotes the total number of qubits in the DQFNN.Here we reduce the number of parameters by using the symmetry structure in the Hamiltonian we discussed in our main text.We restrict each perceptron to only involve local connection between i-th and (i + 1)-th layers.The j-th qubit in the (i + 1)-th layer is only grouped with the (j mod n i )-th and [(j + 1) mod n i ]-th qubits in the i-th layer into the same perceptron, where mod is the module operation.So that there are n i+1 perceptrons between i-th layer and (i + 1)-th layer.The total number of parameters is which means that the number of parameters grows linearly with the system size.This leads to the computational complexity O(N 3 All ) to the DQFNN approach.

FIG. 2 .
FIG. 2. (a) A DQFNN with layer structure (2, 2, 3, 3, 5) is used to solve the dissipative transverse field Ising model [see Eq. (9)].Here, each quantum perceptron consists of a unitary acting on twoinput and one-output qubits and this unitary is parametrized by a simple quantum circuit (black-dotted box), which contains twelve Euler rotations and six controlled-not gates.(b)Numerical results for the steady state of the dissipative Ising model with five spins.The expectation values of three different observables σx = 1 5 ( k σ x k ), σ z = 1 5 ( k σ z k ), and σ x 1 σ x 2 are calculated by both the DQFNN and exact diagonalization (ED) method.Here, we set J = 1 as the energy unit and the dissipation ratio is γ = 1.For the DQFNN approach, the learning rate λ is chosen to be λ = 0.01 × 0.999 Ns , where Ns is the number of steps.

FIG. 3 .
FIG. 3. Numerical results for the time dynamics of the dissipative Ising model with five spins under the periodic boundary condition.For the DQFNN approach, random Gaussian noises with varying strength R are added to the quantum derivatives to account for the imperfections in real experiments.The time step is chosen to be 5 × 10 −3 and the sampling size for every step is 5 × 10 4 .(a), (b), and (c) plot the time dependences of σ z , σ x , and σ x 1 σ x 2 , respectively.(d) plots the real and imaginary parts of δL, which measures the convergence to the steady state.

FIG. 4 .
FIG. 4. Numerical results for the dissipative 2D J1 − J2 model.Here, the simplest two-by-two square lattice is considered and we use a DQFNN with layer structure (4, 4, 4) to obtain the results.The parametrization of the DQFNN is similar to the case for the 1D dissipative Ising model.(a), (b), and (c) show respectively σ z , σ x , and σ x 0 σ x 1 as a function of iteration steps with the dissipation rate chosen to be γ = 1.0.The red-dotted lines represent their corresponding values for the steady state (SS) obtained from exact diagonalization.(d) plots both the real and imaginary parts of δL, which converges to zero after around 1000 iteration steps (|δL| < 10 −2 ), indicating that the output state of the DQFNN indeed converges to the steady state.

4 l 4 l 4 l 4 l 4 l 4 l
(b), we can measure the real and imaginary part of Tr ρ(Θ) Ôρ(Θ ) separately by setting the initial state of the ancillary qubit as (|0 + |1 )/ √ 2 for measuring real part and (|0 + i|1 )/ √ 2 for measuring the imaginary part.Then, we introduce a controlled-Ô gate, which is defined as: Ôc |x |ψ = |x Ôx |ψ ., r, l , r , m, n ρ l, r ρ l , r m, n| l, l r| n r | Ô † | m + i , r, l , r , m, n ρ l, r ρ l , r m| Ô| l n| l r, r | m, n + 1 , r, l , r , m, n ρ l, r ρ l , r m| Ô| l n| l r| n r | Ô † | m , r, l , r , m, n ρ l, r ρ l , r Ô † r , m δ n, l δ m, l δ n, r + i , r, l , r , m, n ρ l, r ρ l , r Ô m, l δ n, l δ m, r δ n, r + 1 , r, l , r , m, n : 1.One lattice site is flipped both in the left state and right state.2. One lattice site is flipped either in the left state or right state.3. Two neighboring lattice sites are flipped either in the left state or right state.4. All the lattice sites in the right state and left state are flipped.5.A new random configuration of the right state and left state are generated.