Simulation-assisted learning of open quantum systems

Models for open quantum systems, which play important roles in electron transport problems and quantum computing, must take into account the interaction of the quantum system with the surrounding environment. Although such models can be derived in some special cases, in most practical situations, the exact models are unknown and have to be calibrated. This paper presents a learning method to infer parameters in Markovian open quantum systems from measurement data. One important ingredient in the method is a direct simulation technique of the quantum master equation, which is designed to preserve the completely-positive property with guaranteed accuracy. The method is particularly helpful in the situation where the time intervals between measurements are large. The approach is validated with error estimates and numerical experiments.


Introduction
Quantum learning has recently emerged as a field at the intersection of quantum physics and computer science [1,2,3].Its primary focus is on the estimation and characterization of energy levels and interaction coefficients within a quantum system.Such effort has become increasingly crucial for characterizing, optimizing, and controlling quantum systems.One particularly explored area is Hamiltonian learning, where the emphasis is to infer the Hamiltonians [4,5,6,7,8,9,10,11,12,13,14,15,16,17] , which ultimately determine the system's energy spectrum and dynamic behavior.Compared to direct simulation approaches, the learning process attempts to map observation data back to the model parameters.An important extension of Hamiltonian learning is the parameter identification of open quantum systems [18,19,20] , ones that continuously interact with their environment [21].The ultimate goal is to reconstruct the interactions in the Lindblad-Gorini-Kossakowski-Sudarshan quantum master equation (QME) [22,23]: where H is the system Hamiltonian and V j 's are jump operators that arise from the interactions with the environment [21].
In practice, the learning tasks are usually formulated as an optimization problem, which is then solved by an iterative method, and as such, one needs to frequently access Ke Wang: wangke.math@psu.eduXiantao Li: Xiantao.Li@psu.eduAccepted in Quantum 2024-07-03, click title to verify.Published under CC-BY 4.0.
the objective function, e.g., through the expectations of some observables.The ability to acquire such quantities might be limited by the measurement capabilities.Furthermore, the availability of the gradient of the objective function, a necessary ingredient for a fast optimization algorithm, might not be attainable through experiments either.In this paper, we propose a simulation-assisted method to address these issues.In contrast to the equation-based methods [18,19] to identify Lindbladians in (1), we formulate a trajectorybased framework, where the observations A(t) = tr Ae tL ρ(0) are fit by simulating the QME (1).The main departure from equation-based methods [19,24,20] is that we do not assume access to the quantum system.We designed a classical learning algorithm with a given time series measurements as input.Namely, we do not assume the flexibility of choosing the time steps or the utilization of classical shadow tomography to make further measurements.For example, the efficiency in the sample complexity in the approach [19] requires choosing the initial and end time carefully.Part of the challenges comes from the presence of noise in the data, which complicates the numerical differentiation procedure.In our approach, the accuracy is controlled by chopping measurement intervals ∆t into smaller steps δt, and in between we simulate the Lindblad dynamics with the proposed parameters.
Another notable difference from the equation-based approaches [18,19] is that our optimization problem does not require the expectations associated with the right-hand side of Eq. ( 1), i.e., tr OV j ρV † j , tr OV † j V j ρ , tr OρV † j V j , and thus fewer features are needed from the data set.
To enable such a simulation-assisted approach, we propose a simulation algorithm, denoted here by M δt (t), that has global error δt 2 , in that e Lt − M δt (t) = O δt 2 , we choose a semi-implicit algorithm, which is stable even when the coherent term in the QME (1) has large coefficients.This makes it more robust in practice than the second-order method in Breuer and Petruccione [21].Furthermore, we also show that the method has a completely positive property and it can be easily written in a Kraus form.This makes it possible to implement such a simulation method on a quantum computer as well.This is done by unraveling the Lindblad dynamics to a stochastic Schrödinger equation [21], followed by a stochastic expansion [25].
Another practical aspect of our approach is efficient optimization.With the Kraus representation of our numerical approximation, the calculation of the gradient is streamlined.This allows us to apply the Levenberg-Marquardt algorithm, which has very rapid convergence [26,27].With mild assumptions on the fitting error, we will show how the approximation error from the numerical simulation and the statistical error from the measurements affect the performance of the parameter identification.
The rest of the paper is organized as follows: We present a general learning framework in the form of a nonlinear least squares problem in Section 3.1, and explain the difference and connections to existing works in Section 3.2.To emphasize the algorithmic details, we present the methods without reference to specific open quantum systems, also with the hope that the methods can be applied to a broader class of problems.Then in Section 3.3 and Section 3.4, we present the specific simulation method and how it is integrated with the optimization procedure.In Section 3.5 and Section 4, we present some error analysis and results from some numerical experiments, where we detail the specific applications of open quantum systems.

Notations and Terminologies
Throughout this paper, the Euclidean norm is denoted by ∥ . The density operator ρ ∈ C d×d is represented as a semi-positive definite matrix with trace 1, properties that will be expressed as ρ ⪰ 0 and tr(ρ) = 1.An emphasis will be placed on quantum evolutions that are trace-preserving and completely positive.These properties are formalized in terms of dynamic maps C d×d → C d×d , e.g.see [28].In particular, if tr(A(ρ)) = tr(ρ) for all ρ, then A is said to be trace-preserving.Similarly, A is said to be positive if ρ ⪰ 0 ⇒ Aρ ⪰ 0. Further, A is said to be completely positive if A ⊗ I m×m is a positive map for every m ≥ 1.A useful representation of trace preserving and completely positive (CPTP) maps is the Kraus form, which expresses A as follows, In this paper, τ m and t n respectively denote the simulation and measurement times.The corresponding time interval is represented by δt and ∆t, respectively, with the ratio denoted by L = ∆t δt , see Fig. 1.In particular, the simulation times are designed to be shorter or equal to the measurement times.The ratio L can be controlled to obtain the desired accuracy.3 The Learning Framework for Parameter Identification of Lindbladians

A Least Squares Formulation of the Learning Problem
In practice, the unknown parameters in an open quantum system can appear in both the Hamiltonian and the dissipative terms in the QME (1).To indicate the role of the parameters, we first rewrite Eq. ( 1) as follows, ( We can combine the parameters by writing θ : will be called the Lindbladian.We refer to the two sets of parameters as Hamiltonian and dissipative parameters, respectively.When L D = 0, the problem is reduced to a Hamiltonian learning problem.
We assume that the measurement data is a time series We have set t 0 = 0, and for simplicity we assume a uniform time-lapse ∆t between experiments, although the extension to general distributions of measurement times is straightforward.Also indicated in Eq. ( 4) is the dependence of the solution of (3) on the parameters, i.e., ρ(t n ; θ) = e tnL(θ) ρ(0).We now formulate the learning problem based on the consistency in Eq. ( 4), as a least squares problem, where the objective function is given by, Here, y * k,n is interpreted as the measurement data with respect to the true, but unknown parameter θ * ; θ * is the solution of the least squares problem, i.e., θ * = argminϕ(θ).Due to the typical nonlinear dependence of ρ(t n ; θ) on θ, the optimization problem is a nonlinear least squares problem [29].
An alternative viewpoint is a maximum likelihood estimate (MLE) from parameter estimation methods for dynamical systems [30].Namely, the measurement outcome y k,n comes with a Gaussian noise ϵ k,n , Then Eq. ( 6) is the log-likelihood function.

Related Works
The QME (1) is an important alternative to study electron transport properties, which are important in material science [31,32].Attempts have been made to integrate the model with existing parallel code [33] to simulate systems with many electrons.These simulation algorithms can be regarded as a bottom-up approach, where the bath operators are assumed in advance.Another increasingly important application is quantum computing, where the quantum circuits are subject to environmental noise, and many quantum error mitigation schemes require the knowledge of an error model, i.e., the channel interactions and strength [34,35,36].
Bairey et al. [18] proposed to learn the coefficients in the QME (1) by acting the steady-state equation on a collection of observables.These equations can thus be rewritten in terms of the expectations by applying the adjoint of the Lindbladians to the observables, leading to equations of Ehrenfest form.For example, by letting ⟨A⟩(t) := tr(Aρ(t)), one can derive from Eq. (1), The operators on the right-hand side can be expressed on a known basis with unknown parameters, which reformulate the equation in terms of measurement data and parameter values.At a steady state, the time derivative drops out, and the approach by Bairey et al.
[18] yields a system of equations, usually over-determined, for the unknown parameters.Franca et al. [19] extended this framework to dynamics problems, with the time derivatives estimated from polynomial approximations.As alluded to in the previous section, these two approaches can be summarized as equation-based: in that the objective function is set up so that the equation ( 7) holds the best it can.An important practical issue is that the time interval for the measurements, denoted in this article by ∆t, might be large, in which case, the accuracy of the inference is limited by the approximation of the time derivatives.If there is no means to decrease the sample time interval, we face inherent inaccuracies in the learning method.More importantly, the measurement data in practical applications always contain noise, leading to a subtle challenge: For the numerical differentiation procedure to be accurate, ∆t must be sufficiently small.However, a small ∆t can amplify the effect of measurement noise.This difficulty has been identified and analyzed in other applications [37,38].The approach in [19] mitigated this issue by choosing the measurement times based on the Chebyshev measure, which on the other hand, requires an upper bound on the time interval [0, t max ].Another challenge is that the right-hand side of Eq. (7) may involve many expectations that have to be obtained from measurements.It is also worthwhile to mention that Boulant et al. [39] also proposed numerical algorithms for reconstructing Lindblad operators and highlighted the importance of the CP property.But their method to ensure this property is quite involved.The more recent work [40] aims at reconstructing Lindbladians based on measurements at multiple times with an optimization problem.Their method is based on experimental measurements and the least squares problem is formulated in terms of the discrete probability induced by the observables.Numerical simulations are only used to determine a stopping criterion.Compared to this experiment-based approach, the simulation-assisted approach only works with the original set of measurement data.More importantly, the numerical simulation provides the gradient of the objective function, which can substantially speed up the optimization process.Another related approach to Lindbladian learning is the learning algorithms from Gibbs states [12,14], since the Gibbs state could be reached by certain Lindblad dynamics.
While the equation-based methods [18,19] fundamentally differ from the method presented in this paper, the two approaches can be combined to accelerate the parameter estimation procedure.Indeed, the current problem can be viewed, in the broader context, as parameter estimation of ODEs [30], where the polynomial approximation of derivatives is known as polynomial regression.It has been used as a preparation of a two-step estimation procedure [41,42].The parameters obtained from the linear regression, e.g., in the approach of Franca et al. [19] can be used as an initial guess for a trajectory-based approach.

Computing the Objective Function using Direct Simulations of the QME
Solving the optimization problem in (5) requires access to the expectation of A (k) for an approximate parameter θ, and such data are unavailable from experiments.This is treated by an efficient numerical algorithm for solving the QME (3).Here we outline a simple derivation of numerical methods, which can be implemented on either classical or quantum computers.
In a nutshell, to approximate the solution operator e Lt of the QME (1), a simulationassisted algorithm uses an approximation M δt (t) and minimizes the difference between e L † (tn) A and M † δt (t n )A, by taking the trace with ρ(0).Another advantage of this approach is that we no longer have to measure the terms on the right-hand side, as was done in [18,19].Thus the number of measurements can be significantly reduced.
The QME encodes a trace-preserving and completely positive (TPCP) map [22,23], which in the Kraus form, can be generally written as [28] (Choi-Kraus'theorem), where j F † j F j = I [21].Therefore, if such a Kraus form associated with the solution e tnL(θ) of the QME (3) can be found, then the objective function in Eq. ( 6) and its gradient can be directly evaluated to carry out the optimization task.Another interesting aspect is that this procedure also places the problem in the general framework of learning a quantum system [3].Boulant et al. [39] have demonstrated that the CP property increases the overall robustness of the parameter estimation procedure.
Our approach starts by first unraveling the Lindblad equation (3) into the stochastic Schrödinger equation [21], Here |ψ(t)⟩ represents a stochastic realization of a quantum state ρ(t), in the sense that, In the stochastic equation, W j;t 's are independent Brownian motions that incorporate the noise from the bath.
Although in practice the trace-preserving property can be ensured by simple scaling, the completely positive property is often destroyed by classical ODE methods [43].A simple example is the Euler's method, e.g., applied to a simple Lindblad equation: The right hand side is a Kraus form, but not in a diagonal form: It can be written as, But the corresponding matrix (a i,j ) is clearly not positive definite.Consequently, standard ODE methods may not induce a CP map.
The equivalence between the QME (1) and Eq. ( 10) is the key ingredient for constructing an approximation of the QME that preserves the CP property.As a quick demonstration, we first consider a time discretization of the stochastic Schrödinger equation ( 9) by the semi-implicit Euler method [25].Toward this end, we define numerical time steps,

and an approximate wave function |ψ
It is important to point out that the step size δt is a numerical parameter, and it can be much smaller than the measurement time ∆t.In order for the simulation to produce expectations at the same time steps as the experiments, choose δt such that, The semi-implicit Euler method follows a time-marching step and implements an iteration formula for |ψ m ⟩ as follows, The terminology "semi-implicit" comes from the treatment of the noise: It is only sampled from the current time step and if V j = 0, this method is the standard implicit Euler method for an ordinary differential equation.In Eq. ( 13) δW j,m 's are independent Gaussian random variables with mean zero and variance δt.In addition, the matrix G is given by, It also appears in the structure-preserving scheme [43].
We can express the method in a compact form At this point, an approximation method for the density operator ρ can easily be obtained.Specifically, let By taking expectations of (15), we arrive at, Here we have used E δW j,m ] = 0 and E δW j,m δW k,m = δtδ j,k .
Clearly, this can be written in the Kraus form (8). The corresponding Kraus operators are given by, This Kraus form from the semi-implicit Euler method is summarized in Lemma 1.
Lemma 1.The semi-implicit Euler method induces a Kraus form, i.e., The approximation property can be verified by direct expansions of the Kraus operators in Eq. (17).
In practice, such first-order methods have limited accuracy.To ensure better performance, we consider the second-order implicit approximation method for stochastic differential equations [25,Chapter 15].When applied to Eq. ( 9), the method can be written as, In this expression, δ Ŵ is an approximation of an increment of the Brownian motion.In addition, U j 1 ,j 2 are independent two-point distribution, defined as, We now state the approximation property of this method.
Lemma 2. The implicit second-order approximation induces a Kraus form, i.e.
In addition, the one-step error is given by, Proof.Similar to Lemma 1, we define Following the second-order implicit scheme (20), each iteration of the density operator is as follows, where the Kraus operators are given by, Notice that the one-step error is of the third order, so we can simplify the iteration formula without changing the error order by letting This will simplify the calculation of the gradient.The rest of the proof is given in Appendix A.2.
We chose an implicit method due to its numerical stability.For example, one can show that the spectral radius of each Kraus operator has an upper bound that is independent of H.This can be seen from the observation that the real part of the matrix G is positive definite.Therefore, the spectral radius of F 0 is less or equal to 1.For 0 < j ≤ N V , F j has a spectral radius less or equal to the spectral radius of V j √ δt.The matrix inverse of I − 1 2 Gδt will certainly complicate the numerical implementation.One robust approach is to solve the associated linear system of equations by bi-conjugate gradient method (Bi-CGSTAB) [44], which involves matrix-vector multiplication and orthogonalization.On the other hand, if H does not involve large eigenvalues, an explicit method can be used to replace Eqs. ( 13) and (20) to simplify the implement, e.g., using the second-order Itô-Taylor expansion [25].

Gradient Evaluations for Gradient-based Optimization
The problem of finding the optima of ϕ(θ) can be seen as the nonlinear least squares problem where we used the textbook notations [29]: is the residual error and R(θ) = (r k,n ) k,n will be referred to as the residual vector.We are interested in the scenario when the dimension of the residual is larger than the number of parameters, i.e., the overdetermined regime.
One practical iteration method for solving the least squares problem is the Levenberg-Marquardt (LM) algorithm, which, starting with an initial guess, θ (0) , involves updating the parameters iteratively, It can be seen as a combination of the Gauss-Newton method and the gradient descent algorithm.In Eq. ( 27), R ′ denotes the gradient of the residual with respect to the parameters.The parameter ν k serves as a regularization, and in practice, it can be chosen to be proportional to the norm of the residual error.
The convergence of the LM method has been established under very mild conditions.Here we follow [26,27] and pose the following local condition: There exists a constant C, such that, for all θ in a neighborhood of θ * .
Theorem 1. [26, Theorem 2.1] Assume that Eq. (28) holds and R ′ is Lipschitz continuous.Let the parameter ν k be If the initial guess θ (0) is sufficiently close to θ * , then there exists a constant, such that the iterations from (27) satisfy quadratic convergence, The convergence speed is faster than the super-linear convergence proved in [29].
The quadratic convergence property makes the LM method (27) an extremely useful algorithm.Global convergence has also been analyzed in [26,27].This is often accomplished by using a line search algorithm.(27).One alternative, which is particularly useful for large scale problems, is the gradient descent algorithm,

Remark 1. The LM algorithm's fast convergence comes with the cost of computing the inverse of a matrix in
where α k is the learning rate.Such an algorithm can often achieve linear convergence in a neighborhood of the minimizer [45].
In parameter estimation problems, the condition in Eq. ( 28) is viewed as a local identifiability condition [30].Meanwhile, the Lipschitz condition can be verified by assuming that the first order and second derivatives of L(θ) are bounded, which holds trivially if the Lindbladian has a linear dependence on the parameters.To elaborate on this, a partial derivative of the residual error is given by, This leads us to consider Γ α (t, θ) := ∂ θα ρ(t; θ), which from Eq. ( 3), follows the differential equation, Using the contraction property of e tL [46], we find that, To examine the Lipschitz continuity of the partial derivatives, we define, which follows the equation, and a simple bound follows, As a result, if the first and second order derivatives of L(θ) are bounded, R ′ fulfills the Lipschitz condition.
Meanwhile, in our learning task, the residual function is subject to approximation error.To incorporate these errors, we can follow the proof in [26], where the search direction at each step of the LM algorithm involves a linear least-square (LS) problem.Based on the classical sensitivity analysis for LS [47], an O(ϵ) perturbation of R(θ) and R ′ (θ) will introduce an O(ϵ) error in the iteration formula.Therefore we have Proposition 1.Let ϵ > 0. Suppose that the residual vector R and R ′ is subject to an ϵ perturbation, R → R + ∆R, with ∥∆R∥ ≤ ϵ and ∥∆R ′ ∥ ≤ ϵ, for all θ in a neighborhood of θ * .Then the LM iterations, under the same conditions on R as in the previous theorem, exhibit an approximate quadratic convergence, where C is independent of ϵ.
The LM algorithm in Eq. ( 27) requires access to the gradient of the residual function.To facilitate the calculation of the gradient, we first take the derivative of the Kraus form, Here θ refers to one parameter in θ.The entire gradient can be computed by visiting all the components in θ.For clarity, we express the results as an inner product in the space of Hermitian matrices.Namely, for any Hermitian matrices A and B, we define, Now we show how the gradient of the objective function can be computed from a back propagation procedure.
Lemma 3. Assume that the iteration ρ m → ρ m+1 follows a Kraus form, Then, for any Hermitian operator A, the following identify holds, where, K * stands for the adjoint of K. Namely, Proof.By the cyclic property of the trace operator, we have In light of this Lemma, we have, Theorem 2. When the density operator is simulated by the form of ρ m+1 = K[ρ m ] := j F j ρ m F † j for each simulation time step δt, the explicit form of the derivative of the objective function ϕ(θ) can be written as where ) ] can be viewed as a back-propagated operator.
Proof.The detailed proof can be found in Appendix A.1.
By using Theorem 2, the explicit expression of ∂ θα ϕ(θ) can be obtained once we know the derivative of F j .For instance, for the first-order semi-implicit Euler method, Similarly, for the second-order implicit approximation (23), we have, We now discuss the time complexity of the overall learning algorithm.Notice that simulating Lindblad dynamics for the time duration [0, T ], due to the second order accuracy, incurs complexity that is proportional to the number of time steps, i.e., LN T = O ϵ −1/2 T 3/2 .Meanwhile, the quadratic convergence property of the Levenberg-Marquardt algorithm implies that the number of iteration steps is of the order O(log log ϵ −1 ) such that the optimization error is within ϵ.In addition, with direct calculation, we find that the time complexity of calculating the objective function value and its gradient is , respectively, where N M M denotes the complexity associated with a matrix multiplication, e.g., F j ρF † j .Thus the overall complexity is where Õ ignores logarithmic factors.

For a specific n-qubit open quantum system described by Lindblad master equation, where the Hamiltonian H is k-local and the dissipation part L only contains single qubit terms, the overall complexity of the learning algorithm based on one and two-qubit Pauli observables is
Implicit in the final bound in the above theorem is a k-dependent prefactor, which depends on the structure of the locality terms.The dependence of the bound on n comes from the fact that for this specific quantum system, we have

Quantifying the Error in the Parameter Identification
Similar to modern machine learning problems, there are mainly three sources of error in a quantum learning problem.Specifically, as can be seen from the objective function Eq. ( 6), 1.The estimated values, tr A (k) ρ(t n ; θ) are replaced with tr A (k) ρ nL (θ) (Recall that ∆t = Lδt, and ρ(t n , θ) = ρ(τ Ln , θ) ≈ ρ Ln (θ)).This can be regarded as a function approximation error.To include this error, we define, Using the second-order method (23 ) perturbation of the residual error in the original objective function (6).
2. The data y * k,n , as indicated in (4), come from measuring a set of observables at different time instances, which are subject to statistical error.To clarify this perspective, we express the measurement values as random variables: N S indicates the number of times the measurements are repeated.The effect of this sampling error can be understood by considering the following objective function, Based on Chebyshev's inequality, R = ( r k,n ) k,n is a ϵ perturbation of residual error in the original objective function (6) with high probability if 3. The optimization problem (5) may lead to a nonlinear system of equations and we may not be able to find the exact solution after a finite number of iterations.This is the optimization error.
The perturbation caused by the numerical approximation (23) is a deterministic perturbation, which from (28), causes a perturbation of the identified parameter.Theorem 4.Under the assumption (28), there exists a δ 0 > 0, such that for all δt < δ 0 , there is a minimizer θ of the least squares problem defined by the residual vector in Eq. (40), for some constant C.
We now turn to the measurement error.Since the observables A (k) 's are bounded, the measurement noise is bounded, and they can be regarded as a sub-Gaussian random variable.The Hoeffding inequality [48] implies that, Lemma 4.There exists σ k,n , such that the following inequality holds for every ϵ > 0 Now by using a union bound, we can estimate the probability, Here we have set Theorem 5.Under the same assumption as in the previous theorem, and further assume that the observation data y k,n are sampled, times, then there is a minimizer θ of the nonlinear least squares with residual in Eq. (42), with high probability.
3.6 Quantum/classical Hybrid Algorithms for Lindblad Simulation We have thus far discussed a classical approach to simulating the Lindblad equation, which is used to evaluate the objective function ϕ(θ) and its gradients as shown in Eq. ( 37).This method is suitable for quantum systems of moderate size, such as those that can be efficiently treated using large-scale parallel algorithms.Alternatively, these dynamics can be simulated directly on quantum computers.Several such algorithms have already been developed [49,50,51,52].Due to the many technical elements involved in these algorithms, we will sketch the overall procedure and leave the detailed presentation of these methods in separate works.As illustrated in Fig. 2, the overall procedure can be regarded as a quantum/classical hybrid algorithm, where some approximation of the parameter θ is fed into the Lindbladian, from which the objective function ϕ(θ) is estimated as an expectation.In addition, a gradient estimation algorithm, such as the improved Jordan's algorithm [53], can be used to estimate the gradient of ϕ(θ).On the other hand, we run the LM algorithm on a classical device to provide an update of the parameter, which re-enters the quantum Lindblad simulation until it reaches convergence.

Numerical Tests
In this section, we present several numerical results to test the effectiveness of our learning algorithm.For the test problem, we consider a quantum system of qubits with dynamics described by Lindblad Eq. ( 1).In particular, we assume that H is linear combination of k−local operators in that each term is acting on at most k qubits.In addition, the jump operators V j are assumed to be 1−qubit Paulis.Specifically, we choose the number of qubits n = 6.The initial states are fixed as the product state with all spins up.Meanwhile, the observables are chosen to be 1 and 2-local Pauli Strings.The performance of our learning algorithm will be quantified with relative and absolute error with respect to 2-norm.The source code is available at [54].

Phase Damping Model with Linear Parameter Dependence
We investigate an n -spin system with dephasing and amplitude damping noise on every qubit [34].The dynamics are described by the QME in Eq. ( 1), with the parametric form of the Hamiltonian part given by, σ α j is the Pauli matrix σ α (α = x, y, z) applied to the jth qubit, and the condition c n,α,β = 0 models open boundaries.
Meanwhile, the dissipation term is parameterized as follows, where the operator T [V ] is defined as The operators σ ± j := 2 −1 (σ x j ± iσ y j ) are introduced.The parameters are expressed as a vector θ, consisting of Hamiltonian coefficients θ H = {e j,α , c j,α,β } and dissipative parameters θ D = (λ 1 , λ 2 ), with 65 unknown parameters in total.To initialize the learning algorithms, these parameters will be generated from a Gaussian distribution, and then held fixed as the exact values of the parameters.
We first test the accuracy of the simulation methods in Section 3.3.As a reference, we conducted direct simulations of the test model in Eqs. ( 47) and ( 48) by using very small step size: δt = 10 −4 .The measurement with 1-qubit local observable A := σ y 2 , will serve as the benchmark solution and is denoted by y exact (t).
Fig. 3 displays the evolution of this observable in the interval t ∈ [0, 10].It contrasts the "exact" solution with results from the semi-implicit Euler method (16) and the second-order implicit method (23) using step size δt = 10 −2 , referred to as y 1 (t), and y 2 (t) respectively.One can observe that the semi-implicit Euler's method yields very good accuracy in the initial period [0, 0.5], but only a qualitatively correct solution in the transient state [0.5, 4].In addition, it tends to smear out the solution in the long run, e.g., in the final period [6,10].But the second-order implicit method offers much better accuracy in the entire time duration.Due to the high accuracy, the measurement data in the following numerical experiments will be generated using the second-order implicit method with δt = 10 −2 .We now employ the Lindbladian simulation scheme to implement the optimization algorithms detailed in Section 3.4 .
Observations were made at time points t n = 0.1, 0.   Extending our analysis, we retained 1-local observables but increased measurement frequency, setting measurement times to τ = t = 0.01, 0.02, • • • , 1 (∆t = δt = 0.01).Fig. 5 indicates that higher measurement frequency slightly enhances parameter identification.The efficiency of the optimization can again be attributed to the rapid convergence of the LM algorithm.
Inspired by the study in [18], where the Lindbladians are learned from steady states , our last experiment uses observation times at a later stage t = 4.1, 4.2, • • • , 5. A separate numerical test verified that this is when the system is beginning to saturate (see Fig. 3 as well).As shown in Fig. 6, the optimization identified all the parameters.Despite the non-monotonic convergence of θ H , the last few iterations exhibited rapid convergence.

A Nonlinear Parametric Model
The model in the previous section involves dissipation terms that are linear with respect to the parameters.In this section, we assume that each jump operator in the dissipation term has a linear parametric form, effectively leading to a nonlinear parametric model.This coincides with the example in Bairey et al. [18].Specifically, the Hamiltonian is of the same linear form as in the previous section while the jump operators are expanded as a linear combination of local Pauli matrices with complex coefficients , Similarly, the parameters are randomly generated by the Gaussian distribution , ).
These parameters will then be fixed and considered to be exact.Subsequently, they are used to generate the observation data y * k,n in (4).We have noticed that the parametric form in Eq. ( 51) is not unique due to the simultaneous appearance of V j and V † j in the Lindbladian.To eliminate the redundancy from the global phase of d j for each j, we set the imaginary part of d j,1 to zeros .The parameter θ is composed by Hamiltonian part θ H = {e j,α , c j,α,β } as in Eq. ( 47) and the dissipation part θ D = {d The simulation time interval is much smaller: δt = 0.01.These experiments (results shown in Fig. 7) indicate that effective parameter identification could be achieved even with the nonlinear parametric model.Moreover, faster convergence was observed when employing all 1 and 2-local observables, which could be possibly caused by a larger constant in the condition (28).To further illustrate the effect of the choice of the observables, we choose 1-local observables, by only keeping σ x j and σ y j 's for each spin, while other configurations of the test remain the same.This modification led to a noticeable increase in the number of iterations to reach a plateau, as demonstrated in Fig. 8.The large optimization error suggests that the limited set of 1-local observables is insufficient for the Lindbladian learning problem.
So far we only considered a small initial guess to show the second-order convergence property of Levenberg-Marquardt algorithm.To illustrate the robustness of our algorithm, we implemented three additional numerical experiments with initial guesses chosen such that the difference between the initial parameter θ (0) and the minimizer θ * is O(1) for both linear and nonlinear cases.We repeat the numerical experiment in Fig. 5, but using a random initial guess with a much larger error ∥θ (0) − θ * ∥ = 2.3650.The convergence is within 10 iterations as shown in Fig. 9. Similarly, we repeat the experiment in Fig. 7 and the initial guess has error ∥θ (0) − θ * ∥ = 2.0359.Our algorithm still converged quite rapidly as depicted in Fig. 10.We now choose θ (0) with a much larger initial error ∥θ (0) − θ * ∥ = 7.8569 in the experiment in Fig. 5.As we can observe from Fig. 11, the algorithm still exhibits convergence but the parameter θ has settled to another minimum.Nevertheless, the objective function is almost zero, indicating an excellent fit to the observation data.In summary, our algorithm demonstrates robust performance even with a considerably larger initial guess.However, the algorithm might end up with another global minimum θ.

Summary and discussions
In this paper, we presented an algorithm to infer parameters in an open quantum system, specifically, in a Lindblad equation.Rather than working with the observation times, we introduce smaller time scales where the solution is obtained by direct numerical simulations.There are two advantages.On one hand, since the objective function is formulated based on trajectories, we no longer need the observables that correspond to the Lindbladian terms that appeared in the approach in [18,19].As a result, the number of observables can be much less than the number of jump operators.On the other hand, the algorithm enables flexible control of the accuracy by using smaller step sizes in the numerical simulation.
In this expression, δ Ŵ is an approximation of an increment of Brownian motion.It was suggested in [25] to sample with the three-point distribution, , P (δ Ŵ = 0) = 2 3 (58) Meanwhile, U j 1 ,j 2 's are independent two-point distributed random variables with In Eq. ( 57), the diffusion operator L j 1 's are defined as follows.
where the Kraus operator

Figure 1 :
Figure 1: Time steps associated with the simulations and measurements.The measurement time interval ∆t is chopped into smaller intervals with δt representing numerical step size.

Figure 2 :
Figure2: A flowchart describing the parameter estimation algorithm using a quantum algorithm for Lindblad simulation.The simulation scheme is implemented by three steps: 1) Block encoding the Hamiltonian term and jump operators in the Lindbladian with coefficients depending on the input parameter θ; 2) Evolve the density operator ρ(t) to time t; 3) Measure observables A(k) with the density operator ρ(t) and time t n , and calculate the objective function value ϕ(θ) and its derivatives ∂ϕ(θ).4) The optimization part is based on Levenberg-Marquardt algorithm on a classical device.

Figure 3 :
Figure 3: Comparison of the accuracy for the simulation methods in Section 3.3, applied to a 6 -spin chain with dephasing and amplitude damping noise.The density operator ρ is measured with 1-qubit local observables A := σ y 2 and the measurement at each time point denoted by y n = tr(Aρ(t n )).The blue solid line represents the "exact" solution generated with a very small step size.
2, • • • , 1 with measurement interval ∆t = 0.1 and in the learning algorithm, the underlying Lindblad dynamics were simulated with a smaller step δt = 0.01 (L = 10) .The performance of the algorithms was tested on the basis of all 1-local observables or all 1 and 2-local observables, generated by Pauli matrices.Numerical results in Fig. 4 include the objective function value and the relative error ∥θ − θ * ∥/∥θ * ∥.Notably, the Levenberg-Marquardt algorithm demonstrated rapid convergence.Furthermore, with all 1 and 2-local observable, the optimization yields slightly faster parameter identification.It demonstrates that an increased number of observables enhances identification.

( 61 )
By taking expectations, this time-marching scheme induces an iteration formula for the density operator ρ m = E[|ψ m ⟩ ⟨ψ m |], as follows.