# Simulation-assisted learning of open quantum systems

Ke Wang<EMAIL_ADDRESS>Department of Mathematics, Pennsylvania State
University  Xiantao Li<EMAIL_ADDRESS>Department of Mathematics,
Pennsylvania State University

###### Abstract

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.

##  I 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], 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 [12, 13], ones that
continuously interact with their environment [14]. The ultimate goal is to
reconstruct the interactions in the Lindblad-Gorini-Kossakowski-Sudarshan
quantum master equation (QME) [15, 16]:

$\frac{d}{dt}\rho=-i[H,\rho]+\sum_{j=1}^{N_{V}}(V_{j}\rho
V^{{\dagger}}_{j}-\frac{1}{2}V_{j}^{{\dagger}}V_{j}\rho-\frac{1}{2}\rho
V_{j}^{{\dagger}}V_{j}),\\\ $ (1)

where $H$ is the system Hamiltonian and $V_{j}$’s are jump operators that
arise from the interactions with the environment [14].

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 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 [12, 13] to identify Lindbladians in (1), we
formulate a _trajectory-based_ framework, where the observations
$A(t)=\tr(Ae^{t\mathcal{L}}\rho(0))$ are fit by simulating the QME (1). An
important consequence is that the dependence of accuracy on the measurement
time-lapse $\Delta t$ is replaced by the dependence on the time step $\delta
t$ of the simulation. This offers the flexibility to achieve an accuracy not
limited by experimental setups. For example, when the measurement time $\Delta
t$ is large, the numerical simulations fill the gap by generating the solution
of the QME (1) in between. Another advantage is that, unlike the equation-
based approaches [12, 13], our optimization problem does not require the
expectations associated with the right hand side of Eq. 1.

To enable such a simulation-assisted approach, we propose a simulation
algorithm, denoted here by $\mathcal{M}_{\delta t}(t)$, that has global error
$\delta t^{2}$, in that $e^{\mathcal{L}t}-\mathcal{M}_{\delta
t}(t)=\mathcal{O}\left(\delta t^{2}\right).$ 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 [14]. 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 [14], followed by a stochastic expansion
[17].

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 [18, 19]. 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 II.1,
and explain the difference and connections to existing works in Section II.2.
Section II.3 and Section II.4 present the specific simulation method and how
it is integrated with the optimization procedure. In Section II.5 and Section
III, we present some error analysis and results from some numerical
experiments.

##  II The Learning Framework for Parameter Identification of Lindbladians

###  II.1 A Least Squares Formulation of the Learning Problem

In practice, the unknown parameter 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,

$\displaystyle\frac{d}{dt}\rho$
$\displaystyle=\mathcal{L}_{H}(\bm{\theta}_{H})\rho+\mathcal{L}_{D}(\bm{\theta}_{D})\rho,$
(2) $\displaystyle\mathcal{L}_{H}(\bm{\theta}_{H})\rho=-i[H,\rho],$
$\displaystyle\quad\mathcal{L}_{D}(\bm{\theta}_{D})\rho=\sum_{j=1}^{N_{V}}\left(V_{j}\rho
V_{j}^{{\dagger}}-\frac{1}{2}V_{j}^{{\dagger}}V_{j}\rho-\frac{1}{2}\rho
V_{j}^{{\dagger}}V_{j}\right).$

We can combine the parameters by writing
$\bm{\theta}:=(\bm{\theta}_{H},\bm{\theta}_{D})$;
$\mathcal{L}(\bm{\theta}):=\mathcal{L}_{H}(\bm{\theta}_{H})+\mathcal{L}_{D}(\bm{\theta}_{D})$
will be called the _Lindbladian_. We refer to the two sets of parameters as
Hamiltonian and dissipative parameters, respectively. When
$\mathcal{L}_{D}=0$, the problem is reduced to a Hamiltonian learning problem.

We assume that the measurement data is a time series
$y_{k,n},n=1,2,\cdots,N_{T},$ which correspond to observables
$A^{(k)},k=1,\cdots,N_{O}$, at time instances $t_{n}=n\Delta t,\quad 0\leq
n\leq N_{T},$

$y_{k,n}:=\tr\big(A^{(k)}\rho(t_{n};\bm{\theta})\big{missing}).$ (3)

We have set $t_{0}=0,$ and for simplicity we assume a uniform time-lapse
$\Delta t$ between experiments, although the extension to general
distributions of measurement times is straightforward. Also indicated in Eq. 3
is the dependence of the solution of (2) on the parameters, i.e.,
$\rho(t_{n};\bm{\theta})=e^{t_{n}\mathcal{L}(\bm{\theta})}\rho(0).$

We now formulate the learning problem based on the consistency in Eq. 3, as a
least squares problem,

$\min_{\bm{\theta}}\phi(\bm{\theta}),$ (4)

where the objective function is given by,

$\phi(\bm{\theta})=\frac{1}{2N_{O}N_{T}}\sum_{\overset{k=1,\cdots,N_{O}}{n=1,\cdots,N_{T}}}\absolutevalue{y_{k,n}-y_{k,n}^{*}}^{2},\quad
y_{k,n}^{*}:=\tr\big(A^{(k)}\rho(t_{n};\bm{\theta}^{*})\big{missing}),$ (5)

Here, $y_{k,n}^{*}$ is interpreted as the measurement data with respect to the
true, but unknown parameter $\bm{\theta}^{*}$; $\bm{\theta}^{*}$ is the
solution of the least squares problem, i.e.,
$\bm{\theta}^{*}=\mathrm{argmin}\phi(\bm{\theta})$. Due to the typical
nonlinear dependence of $\rho(t_{n};\bm{\theta})$ on $\bm{\theta}$, the
optimization problem is a nonlinear least squares problem [20].

An alternative viewpoint is a maximum likelihood estimate (MLE) from parameter
estimation methods for dynamical systems [21]. Namely, the measurement outcome
$y_{k,n}$ comes with a Gaussian noise $\epsilon_{k,n}$,

$y_{k,n}=\text{tr}\big{(}A^{(k)}\rho(t_{n};\bm{\theta})\big{)}+\epsilon_{k,n}.$

Then Eq. 5 is the log-likelihood function.

###  II.2 Related Works

The QME (1) is an important alternative to study electron transport
properties, which are important in material science [22, 23]. 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 [24, 25, 26].

Bairey et al. [12] 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 $\langle A\rangle(t):=\text{tr}(A\rho(t))$, one can derive
from Eq. 1,

$\frac{d}{dt}\langle
A\rangle(t)=-i\langle[A,H]\rangle(t)+\sum_{j=1}^{N_{V}}\left(\langle
V_{j}^{{\dagger}}AV_{j}\rangle(t)-\frac{1}{2}\langle
AV_{j}^{{\dagger}}V_{j}\rangle(t)-\frac{1}{2}\langle
V_{j}^{{\dagger}}V_{j}A\rangle(t)\right).$ (6)

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. [12] yields a system of equations, usually
over-determined, for the unknown parameters. Franca et al. [13] 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 (6) holds the best it can. An
important practical issue is that the time interval for the measurements,
denoted in this article by $\Delta t$, might be large, in which case, the
accuracy of the inference is limited by the approximation of the time
derivatives. Another challenge is that the right hand side of (6) may involve
many expectations that have to be obtained from measurements. It is also
worthwhile to mention that Boulant et al. [27] 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 [28] 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.

While the equation-based methods [12, 13] 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 [21], where the
polynomial approximation of derivatives is known as polynomial regression. It
has been used as a preparation of a two-step estimation procedure [29, 30].
The parameters obtained from the linear regression, e.g., in the approach of
Franca et al. [13] can be used as an initial guess for a trajectory-based
approach.

###  II.3 Computing the Objective Function using Direct Simulations of the QME

Solving the optimization problem in (4) requires access to the expectation of
$A^{(k)}$ for an approximate parameter $\bm{\theta}$, and such data are
unavailable from experiments. This is treated by an efficient numerical
algorithm for solving the QME (2). 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^{\mathcal{L}t}$ of the
QME (1), a simulation-assisted algorithm uses an approximation
$\mathcal{M}_{\delta t}(t)$ and minimizes the difference between
$e^{\mathcal{L}^{\dagger}(t_{n})}A$ and $\mathcal{M}_{\delta
t}^{\dagger}(t_{n})A$, by taking the trace with $\rho(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 [12, 13]. Thus the number of measurements can be
significantly reduced.

The QME encodes a trace-preserving and completely positive (TPCP) map [15,
16], which in the Kraus form, can be generally written as [31] (Choi-
Kraus’theorem),

$\rho_{n+1}=\mathcal{K}[\rho_{n}]:=\sum_{j}F_{j}\rho_{n}F_{j}^{\dagger},$ (7)

where $\sum_{j}F_{j}^{{\dagger}}F_{j}=I$[14]. Therefore, if such a Kraus form
associated with the solution $e^{t_{n}\mathcal{L}(\bm{\theta})}$ of the QME
(2) can be found, then the objective function in Eq. 5 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. [27] have demonstrated that
the CP property increases the overall robustness of the parameter estimation
procedure.

Our approach starts by first unraveling the Lindblad equation (2) into the
stochastic Schrödinger equation [14],

$\mathrm{d}\ket{\psi(t)}=\left(-iH\ket{\psi(t)}-\frac{1}{2}\sum_{j=1}^{N_{V}}V_{j}^{{\dagger}}V_{j}\ket{\psi(t)}\right)\mathrm{d}t+\sum_{j=1}^{N_{V}}V_{j}\ket{\psi(t)}\mathrm{d}W_{j;t}$
(8)

Here $\ket{\psi(t)}$ represents a stochastic realization of a quantum state
$\rho(t)$, in the sense that,

$\rho(t)=\mathbb{E}\big{[}\outerproduct{\psi(t)}{\psi(t)}\big{]}.$ (9)

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 [32]. The equivalence between the QME (1) and Eq. 9 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 (8) by the semi-implicit Euler method
[17]. Toward this end, we define numerical time steps, $\tau_{m}=m\delta
t,m=0,1,\cdots,M$, and an approximate wave function $\ket{\psi_{m}}$,

$\ket{\psi(\tau_{m})}\approx\ket{\psi_{m}},\;m\geq 0.$ (10)

It is important to point out that the step size $\delta t$ is a numerical
parameter, and it can be much smaller than the measurement time $\Delta t.$ In
order for the simulation to produce expectations at the same time steps as the
experiments, choose $\delta t$ such that,

$L=\frac{\Delta t}{\delta t}\in\mathbb{N},\quad M=LN_{T}.$ (11)

The semi-implicit Euler method follows a time-marching step and implements an
iteration formula for $\ket{\psi_{m}}$ as follows,

$\ket{\psi_{m+1}}=\ket{\psi_{m}}+G\ket{\psi_{m+1}}\delta
t+\sum_{j=1}^{N_{V}}V_{j}\ket{\psi_{m}}\delta W_{j,m}.$ (12)

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.
12 $\delta W_{j,m}$’s are independent Gaussian random variables with mean zero
and variance $\delta t.$ In addition, the matrix $G$ is given by,

$G=-iH-\frac{1}{2}\sum_{j=1}^{N_{V}}V_{j}^{{\dagger}}V_{j}.$ (13)

It also appears in the structure-preserving scheme [32].

We can express the method in a compact form

$\ket{\psi_{m+1}}=(I-G\delta t)^{-1}\left[I+\sum_{j=1}^{N_{V}}V_{j}\delta
W_{j,m}\right]\ket{\psi_{m}}.$ (14)

At this point, an approximation method for the density operator $\rho$ can
easily be obtained. Specifically, let
$\rho_{m}=\mathbb{E}\big{[}\ket{\psi_{m}}\bra{\psi_{m}}\big{]},m\geq 0$;
$\rho(\tau_{m})\approx\rho_{m}.$ By taking expectations of (14), we arrive at,

$\rho_{m+1}=(I-G\delta t)^{-1}\rho_{m}(I-G\delta t)^{-{\dagger}}+(I-G\delta
t)^{-1}\sum_{j=1}^{N_{V}}V_{j}\rho_{m}V_{j}^{{\dagger}}(I-G\delta
t)^{-{\dagger}}\delta t.$ (15)

Here we have used $\mathbb{E}\big{[}\delta W_{j,m}]=0$ and
$\mathbb{E}\big{[}\delta W_{j,m}\delta W_{k,m}\big{]}=\delta t\delta_{j,k}.$

Clearly, this can be written in the Kraus form (7). The corresponding Kraus
operators are given by,

$\displaystyle F_{0}=$ $\displaystyle(I-G\delta t)^{-1},$ (16) $\displaystyle
F_{j}=$ $\displaystyle(I-G\delta t)^{-1}V_{j}\sqrt{\delta
t},j=1,\cdots,N_{V}.$ (17)

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.,
$\rho_{m+1}=\mathcal{K}[\rho_{m}]=\sum_{j}F_{j}\rho_{m}F_{j}^{{\dagger}}$
where $\sum_{j}F_{j}^{{\dagger}}F_{j}=I+\mathcal{O}(\delta t^{2}).$ In
addition,

$e^{\delta t\mathcal{L}(\bm{\theta})}\rho-\mathcal{K}\rho=\mathcal{O}(\delta
t^{2}).$ (18)

The approximation property can be verified by direct expansions of the Kraus
operators in Eq. 16.

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 [17, Chapter 15]. When applied to Eq. 8, the
method can be written as,

$\displaystyle\ket{\psi_{m+1}}$ $\displaystyle=(I-\frac{1}{2}G\delta
t)^{-1}(I+\frac{1}{2}G\delta
t)\ket{\psi_{m}}+\sum_{j=1}^{N_{V}}(I-\frac{1}{2}G\delta
t)^{-1}(V_{j}+\frac{1}{2}V_{j}G\delta t)\ket{\psi_{m}}\delta\hat{W}_{j}$ (19)
$\displaystyle+\frac{1}{2}(I-\frac{1}{2}G\delta
t)^{-1}\sum_{j_{1},j_{2}=1}^{N_{V}}V_{j_{2}}V_{j_{1}}\ket{\psi_{m}}\left(\delta\hat{W}_{j_{1},m}\delta\hat{W}_{j_{2},m}+U_{j_{1},j_{2},m}\right).$

In this expression, $\delta\hat{W}$ is an approximation of an increment of the
Brownian motion. In addition, $U_{j_{1},j_{2}}$ are independent two-point
distribution, defined as,

$\displaystyle P(U_{j_{1},j_{2}}=\pm\delta t)=\frac{1}{2},\forall
j_{2}=1,\cdots,j_{1}-1,$ (20) $\displaystyle U_{j_{1},j_{1}}=-\delta t,$
$\displaystyle U_{j_{1},j_{2}}=-U_{j_{2},j_{1}},\text{ for
}j_{2}=j_{1}+1,\cdots,N_{V}.j_{1}=1,\cdots,N_{V},$

We now state the approximation property of this method.

######  Lemma 2.

The implicit second-order approximation induces a Kraus form, i.e.
$\rho_{m+1}=\sum_{j}F_{j}\rho_{m}F_{j}^{{\dagger}}$, where
$\sum_{j}F_{j}F_{j}^{{\dagger}}=I+\mathcal{O}(\delta t^{3})$. In addition, the
one-step error is given by,

$e^{\delta t\mathcal{L}}\rho=\sum_{j}F_{j}\rho
F_{j}^{{\dagger}}+\mathcal{O}(\delta t^{3}).$ (21)

###### Proof.

Similar to Lemma 1, we define
$\rho_{m}=\mathbb{E}[\ket{\psi_{m}}\bra{\psi_{m}}]$. Following the second-
order implicit scheme (19), each iteration of the density operator is as
follows,

$\displaystyle\rho_{m+1}$ $\displaystyle=(I-\frac{1}{2}G\delta
t)^{-1}(I+\frac{1}{2}G\delta t)\rho_{m}(I+\frac{1}{2}G\delta
t)(I-\frac{1}{2}G\delta t)^{-1}$ (22) $\displaystyle+(I-\frac{1}{2}G\delta
t)^{-1}\sum_{j=1}^{N_{V}}V_{j}(I+\frac{1}{2}G\delta
t)\rho_{m}(I+\frac{1}{2}G\delta t)V_{j}^{{\dagger}}(I-\frac{1}{2}G\delta
t)^{-1}\delta t$ $\displaystyle+\frac{1}{2}(I-\frac{1}{2}G\delta
t)^{-1}\sum_{j_{1},j_{2}=1}^{N_{V}}V_{j_{1}}V_{j_{2}}\rho_{m}V_{j_{2}}^{{\dagger}}V_{j_{1}}^{{\dagger}}(I-\frac{1}{2}G\delta
t)^{-1}(\delta t)^{2}$
$\displaystyle=\sum_{j=0}^{N_{V}^{2}+N_{V}}F_{j}\rho_{m}F_{j}^{{\dagger}}=:\mathcal{K}(\rho_{m})$

where the Kraus operators are given by,

$\displaystyle F_{0}$ $\displaystyle=(I-\frac{1}{2}G\delta
t)^{-1}(I+\frac{1}{2}G\delta t),$ (23) $\displaystyle F_{j}$
$\displaystyle=(I-\frac{1}{2}G\delta t)^{-1}V_{j}(I+\frac{1}{2}G\delta
t)\sqrt{\delta t},j=1,\cdots,N_{V},$ $\displaystyle F_{j_{1}+N_{V}j_{2}}$
$\displaystyle=\frac{1}{\sqrt{2}}(I-\frac{1}{2}G\delta
t)^{-1}V_{j_{1}}V_{j_{2}}\delta t,\quad j_{1},j_{2}=1,\cdots,N_{V}.$

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

$\displaystyle F_{j_{1}+N_{V}j_{2}}=\frac{1}{\sqrt{2}}V_{j_{1}}V_{j_{2}}\delta
t,\quad j_{1},j_{2}=1,\cdots,N_{V}.$ (24)

This will simplify the calculation of the gradient. The rest of the proof is
given in Section A.2.

∎

We choose 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\leq N_{V}$, $F_{j}$ has a
spectral radius less or equal to the spectral radius of $V_{j}\sqrt{\delta
t}$.

###  II.4 Gradient Evaluations for Gradient-based Optimization

The problem of finding the optima of $\phi(\bm{\theta})$ can be seen as the
nonlinear least squares problem

$\phi(\bm{\theta})=\frac{1}{2N_{O}N_{T}}\sum_{\overset{k=1,\cdots,N_{O}}{n=1,\cdots,N_{T}}}|r_{k,n}(\bm{\theta})|^{2}=\frac{1}{2N_{O}N_{T}}R(\bm{\theta})^{T}R(\bm{\theta}),$
(25)

where we used the textbook notations [20]: $r_{k,n}=y_{k,n}-y_{k,n}^{*}$ is
the residual error and $R(\bm{\theta})=(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 over-
determined regime.

One practical iteration method for solving the least squares problem is the
Levenberg-Marquardt (LM) algorithm, which, starting with an initial guess,
$\bm{\theta}^{(0)}$, involves updating the parameters iteratively,

$\bm{\theta}^{(k+1)}-\bm{\theta}^{(k)}=-\left(\nu_{k}I+R^{\prime}(\bm{\theta}^{(k)})^{T}R^{\prime}(\bm{\theta}^{(k)})\right)^{-1}R^{\prime}(\bm{\theta}^{(k)})^{T}R(\bm{\theta}^{(k)}).$
(26)

It can be seen as a combination of the Gauss-Newton method and the gradient
descent algorithm. In Eq. 26, $R^{\prime}$ denotes the gradient of the
residual with respect to the parameters. The parameter $\nu_{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 [18, 19] and pose the following local condition:
There exists a constant $C,$ such that,

$\norm{R(\bm{\theta})}\geq C\norm{\bm{\theta}-\bm{\theta}^{*}},$ (27)

for all $\bm{\theta}$ in a neighborhood of $\bm{\theta}^{*}$.

######  Theorem 1.

[18, Theorem 2.1] Assume that Eq. 27 holds and $R^{\prime}$ is Lipschitz
continuous. Let the parameter $\nu_{k}$ be

$\nu_{k}=\norm{R(\bm{\theta}^{(k)})}^{2}.$

If the initial guess $\bm{\theta}^{(0)}$ is sufficiently close to
$\bm{\theta}^{*}$, then there exists a constant, such that the iterations from
(26) satisfy,

$\norm{\bm{\theta}^{(k+1)}-\bm{\theta}^{*}}\leq
C\norm{\bm{\theta}^{(k)}-\bm{\theta}^{*}}^{2}.$ (28)

The convergence speed is faster than the super-linear convergence proved in
[20]. The quadratic convergence property makes the LM method (26) an extremely
useful algorithm. Global convergence has also been analyzed in [18, 19]. This
is often accomplished by using a line search algorithm.

In parameter estimation problems, the condition in Eq. 27 is viewed as a local
identifiability condition [21]. Meanwhile, the Lipschitz condition can be
verified by assuming that the first order and second derivatives of
$\mathcal{L}(\bm{\theta})$ 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,

$\frac{\partial}{\partial{\theta_{\alpha}}}r_{k,n}(\bm{\theta})=\text{tr}(A^{(k)}\partial_{\theta_{\alpha}}\rho(t_{n};\bm{\theta})).$
(29)

This leads us to consider
$\Gamma_{\alpha}(t,\bm{\theta}):=\partial_{\theta_{\alpha}}\rho(t;\bm{\theta}),$
which from Eq. 2, follows the differential equation,

$\frac{d}{dt}\Gamma_{\alpha}(t,\bm{\theta})=\mathcal{L}(\bm{\theta})\Gamma_{\alpha}(t,\bm{\theta})+\frac{\partial}{\partial{\theta_{\alpha}}}\mathcal{L}(\bm{\theta})\rho(t,\bm{\theta}).$

Using the contraction property of $e^{t\mathcal{L}}$ [33], we find that,

$\norm{\Gamma_{\alpha}(t,\bm{\theta})}\leq
t\norm{\partial_{\theta_{\alpha}}\mathcal{L}(\bm{\theta})}.$

To examine the Lipschitz continuity of the partial derivatives, we define,

$\chi_{\alpha,\beta}:=\partial_{\theta_{\beta}}\Gamma_{\alpha}(t,\bm{\theta}),$

which follows the equation,

$\frac{d}{dt}\chi_{\alpha,\beta}(t,\bm{\theta})=\mathcal{L}(\bm{\theta})\chi_{\alpha,\beta}(t,\bm{\theta})+\frac{\partial^{2}}{\partial\theta_{\beta}\partial\theta_{\alpha}}\mathcal{L}(\bm{\theta})\rho(t,\bm{\theta})+\frac{\partial}{\partial{\theta_{\alpha}}}\mathcal{L}(\bm{\theta})\Gamma_{\beta}(t,\bm{\theta})+\frac{\partial}{\partial{\theta_{\beta}}}\mathcal{L}(\bm{\theta})\Gamma_{\alpha}(t,\bm{\theta}),$

and a simple bound follows,

$\norm{\chi_{\alpha,\beta}(t,\bm{\theta})}\leq
t\norm{\frac{\partial^{2}}{\partial\theta_{\alpha}\partial\theta_{\beta}}\mathcal{L}(\bm{\theta})}+\frac{t^{2}}{2}\big{(}2\norm{\partial_{\theta_{\alpha}}\mathcal{L}(\bm{\theta})}+\norm{\partial_{\theta_{\beta}}\mathcal{L}(\bm{\theta})}\big{)}.$

As a result, if the first and second order derivatives of
$\mathcal{L}(\bm{\theta})$ are bounded, $R^{\prime}$ fulfills the Lipschitz
condition.

Meanwhile, in our learning task, the residual function is subject to
approximation and statistical error. To incorporate these errors, we first
state a generalization of the convergence result, which can be proved by
following the proof in [18],

######  Proposition 1.

Let $\epsilon>0.$ Suppose that the residual vector $R$ is subject to an
$\epsilon$ perturbation,

$R\to R+\Delta R,$

with $\norm{\Delta R}\leq\epsilon$, for all $\bm{\theta}$ in a neighborhood of
$\bm{\theta}^{*}$, and assume that $\Delta R^{\prime}$ is Lipschitz continuous
as well. Then the LM iterations, under the same conditions on $R$ as in the
previous theorem, exhibit an approximate quadratic convergence,

$\norm{\bm{\theta}^{(k+1)}-\bm{\theta}^{*}}\leq
C\norm{\bm{\theta}^{(k)}-\bm{\theta}^{*}}^{2}+\epsilon.$ (30)

The LM algorithm in Eq. 26 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,

$(\partial_{\theta}\mathcal{K})[\rho]=\sum_{j}\partial_{\theta}F_{j}\rho
F_{j}^{{\dagger}}+F_{j}\rho\partial_{\theta}F_{j}^{{\dagger}}$ (31)

Here $\theta$ refers to one parameter in $\bm{\theta}$. The entire gradient
can be computed by visiting all the components in $\bm{\theta}$. 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,

$\big{\langle}A,B\big{\rangle}:=\tr(AB).$ (32)

Now we show how the gradient of the objective function can be computed from a
back propagation procedure.

######  Lemma 3.

Assume that the iteration $\rho_{m}\to\rho_{m+1}$ follows a Kraus form,

$\rho_{m+1}=\mathcal{K}\rho_{m}:=\sum_{j}F_{j}\rho_{m}F_{j}^{\dagger}.$ (33)

Then, for any Hermitian operator $A$, the following identify holds,

$\big{\langle}A,\rho_{m+1}\big{\rangle}=\big{\langle}\mathcal{K}^{*}A,\rho_{m}\big{\rangle},$
(34)

where, $\mathcal{K}^{*}$ stands for the adjoint of $\mathcal{K}$. Namely,

$\mathcal{K}^{*}[A]:=\sum_{j}F_{j}^{\dagger}AF_{j}.$ (35)

###### Proof.

By the cyclic property of the trace operator, we have

$\displaystyle\tr(A\rho_{m+1})=\sum_{j}\tr(AF_{j}\rho_{m}F_{j}^{{\dagger}})=\sum_{j}\tr(F_{j}^{{\dagger}}AF_{j}\rho_{m})=\tr(\mathcal{K}^{*}[A]\rho_{m})$

∎

In light of this Lemma, we have,

######  Theorem 2.

When the density operator is simulated by the form of
$\rho_{m+1}=\mathcal{K}[\rho_{m}]:=\sum_{j}F_{j}\rho_{m}F_{j}^{{\dagger}}$ for
each simulation time step $\delta t$, the explicit form of the derivative of
the objective function $\phi(\bm{\theta})$ can be written as

$\partial_{\theta_{\alpha}}\phi(\bm{\theta})=\frac{1}{N_{O}N_{T}}\sum_{k=1}^{N_{O}}\sum_{n=1}^{N_{T}}\sum_{l=1}^{nL}\left(y_{k,n}-y_{k,n}^{*}\right)\Big{\langle}\partial_{\theta_{\alpha}}\mathcal{K}^{*}A^{(k)}_{nL-l},\;\rho_{nL-l}\Big{\rangle},$
(36)

where $A^{(k)}_{nL-l}:=\big{(}\mathcal{K}^{*}\big{)}^{l-1}[A^{(k)}]$ can be
viewed as a back-propagated operator.

###### Proof.

The detailed proof can be found in Section A.1. ∎

By using Theorem 2, the explicit expression of
$\partial_{\theta_{\alpha}}\phi(\bm{\theta})$ can be obtained once we know the
derivative of $F_{j}$. For instance, for the first-order semi-implicit Euler
method,

$\displaystyle\partial_{\theta}F_{0}$
$\displaystyle=F_{0}\partial_{\theta}GF_{0}\delta t$ (37)
$\displaystyle\partial_{\theta}F_{j}$
$\displaystyle=\partial_{\theta}F_{0}V_{j}\sqrt{\delta
t}+F_{0}\partial_{\theta}V_{j}\sqrt{\delta t},j=1,\cdots,N_{V}$

Similarly, for the second-order implicit approximation (22), we have,

$\displaystyle\partial_{\theta}F_{0}=\frac{1}{2}\delta t(I-\frac{1}{2}G\delta
t)^{-1}\partial_{\theta}G(F_{0}+I)$ (38)
$\displaystyle\partial_{\theta}F_{j}=\frac{1}{2}\delta t(I-\frac{1}{2}G\delta
t)^{-1}\partial_{\theta}GF_{j}+\sqrt{\delta t}(I-\frac{1}{2}G\delta
t)^{-1}\partial_{\theta}V_{j}(I+\frac{1}{2}G\delta t)$
$\displaystyle+\frac{1}{2}\delta t^{3/2}(I-\frac{1}{2}G\delta
t)^{-1}V_{j}\partial_{\theta}G,\quad j=1,\cdots,N_{V}$
$\displaystyle\partial_{\theta}F_{j_{1}+N_{V}j_{2}}=\frac{\delta
t}{\sqrt{2}}(\partial_{\theta}V_{j_{1}}V_{j_{2}}+V_{j_{1}}\partial_{\theta}V_{j_{2}}),\quad
j_{1},j_{2}=1,\cdots,N_{V}$

###  II.5 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. 5,

  1. 1.

The estimated values, $\text{tr}\big{(}A^{(k)}\rho(t_{n};\bm{\theta})\big{)}$
are replaced with $\text{tr}\big{(}A^{(k)}\rho_{nL}(\bm{\theta})\big{)}$
(Recall that $\Delta t=L\delta t$, and
$\rho(t_{n},\bm{\theta})=\rho(\tau_{Ln},\bm{\theta})\approx\rho_{Ln}(\bm{\theta})$).
This can be regarded as a function approximation error. To include this error,
we define,

$\widehat{r}_{k,n}(\bm{\theta})=\text{tr}\big{(}A^{(k)}\rho_{nL}(\bm{\theta})\big{)}-y_{k,n}^{*}.$
(39)

Using the second-order method (22), $\widehat{R}=(\widehat{r}_{k,n})_{k,n}$ is
a $\mathcal{O}(\delta t^{2})$ perturbation of the residual error in the
original objective function (5).

  2. 2.

The data $y_{k,n}^{*}$, as indicated in (3), 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:

$y_{k,n}^{*}\approx\hat{y}_{k,n}:=\frac{1}{N_{S}}\sum_{\ell=1}^{N_{S}}a_{k,n,\ell}.$
(40)

$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,

$\widetilde{r}_{k,n}(\bm{\theta})=\text{tr}\big{(}A^{(k)}\rho(t_{n},\bm{\theta})\big{)}-\hat{y}_{k,n}.$
(41)

Based on Chebyshev’s inequality, $\widetilde{R}=(\widetilde{r}_{k,n})_{k,n}$
is a $\epsilon$ perturbation of residual error in the original objective
function (5) with high probability if $N_{S}=\mathcal{O}(1/\epsilon^{2})$.

  3. 3.

The optimization problem (4) 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 (22) is a deterministic
perturbation, which from (27), causes a perturbation of the identified
parameter.

######  Theorem 3.

Under the assumption (27), there exists a $\delta_{0}>0,$ such that for all
$\delta t<\delta_{0}$, there is a minimizer $\hat{\bm{\theta}}$ of the least
squares problem defined by the residual vector in Eq. 39, such that

$\norm{\hat{\bm{\theta}}-\bm{\theta}^{*}}\leq C\delta t^{2}.$ (42)

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 [34] implies that,

######  Lemma 4.

There exists $\sigma_{k,n}$, such that the following inequality holds for
every $\epsilon>0$

$\mathbb{P}\Big{(}\absolutevalue{\hat{y}_{k,n}-{y}_{k,n}^{*}}>\epsilon\Big{)}\leq
2e^{-\frac{N_{S}\epsilon^{2}}{\sigma_{k,n}^{2}}}.$ (43)

Now by using a union bound, we can estimate the probability,

$\mathbb{P}\Big{(}\sqrt{\sum_{k,n}\absolutevalue{\hat{y}_{k,n}-{y}_{k,n}^{*}}^{2}}>\epsilon\Big{)}\leq
2N_{O}N_{T}e^{-\frac{N_{S}\epsilon^{2}}{N_{O}N_{T}\sigma^{2}}}.$

Here we have set $\sigma^{2}=\frac{1}{N_{O}N_{T}}\sum_{k,n}\sigma_{k,n}^{2}.$

######  Theorem 4.

Under the same assumption as in the previous theorem, and further assume that
the observation data $y_{k,n}$ are sampled,

$N_{S}=\Omega\left(\frac{\sigma^{2}N_{O}N_{T}}{\epsilon^{2}}\log\big(N_{O}N_{T}\big{missing})\right),$
(44)

times, then there is a minimizer $\widetilde{\bm{\theta}}$ of the nonlinear
least squares with residual in Eq. 41, such that

$\norm{\widetilde{\bm{\theta}}-\bm{\theta}^{*}}\leq\epsilon.$ (45)

with high probability.

##  III Numerical Tests

###  III.1 Phase Damping Model with Linear Parameter Dependence

We first consider an $N$-spin system with dephasing and amplitude damping
noise on every qubit [24]. The dynamics are described by the QME in Eq. 1,
with the parametric form of the Hamiltonian part given by,

$H=\sum_{j=1}^{N}H_{j},\quad
H_{j}=\sum_{\alpha=1}^{3}e_{j,\alpha}\sigma_{j}^{\alpha}+\sum_{\alpha,\beta=1}^{3}c_{j,\alpha,\beta}\sigma_{j}^{\alpha}\sigma_{j+1}^{\beta},\quad
j=1,\cdots,N.$ (46)

$\sigma_{j}^{\alpha}$ is the Pauli matrix $\sigma^{\alpha}$ ($\alpha=x,y,z$)
applied to the $j$th qubit, and we choose $c_{N,\alpha,\beta}=0$ to model an
open boundary condition.

Meanwhile, the dissipation term is parameterized as follows,

$\mathcal{L}_{D}\rho=\lambda_{1}\sum_{j=1}^{N}\mathcal{T}[\sigma_{j}^{-}]\rho+\lambda_{2}\sum_{j=1}^{N}\mathcal{T}[\sigma_{j}^{z}]\rho,$
(47)

where the operator $\mathcal{T}[V]$ is defined as

$\mathcal{T}[V]\rho=V\rho V^{{\dagger}}-\frac{1}{2}\\{V^{{\dagger}}V,\rho\\}.$
(48)

In Eq. 47 we defined operators $\sigma_{j}^{\pm}:=2^{-1}(\sigma_{j}^{x}\pm
i\sigma_{j}^{y})$. The density operator $\rho(0)$ is initialized as a product
state with all spins up.

In this example, the parameter $\bm{\theta}$ is composed of Hamiltonian
coefficients $\bm{\theta}_{H}=\\{e_{j,\alpha},c_{j,\alpha,\beta}\\}$ and
dissipative parameters $\bm{\theta}_{D}=(\lambda_{1},\lambda_{2})$.
Specifically, there are 65 unknown parameters when $N=6$. 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 II.3. As a
reference, we conducted direct simulations of the test model in Eqs. 46 and 47
by using very small step size: $\delta t=10^{-4}$, and we compute the 1-qubit
local observable $A:=\sigma_{2}^{y}$. This solution will be regarded as the
exact solution $y_{\text{exact}}(t)$ for comparison purposes. Fig. 1 shows the
history of this observable in the time interval $t\in[0,10].$ Also shown is
the solution computed from the semi-implicit Euler method (15) and the second-
order implicit method (22) using step size $\delta t=10^{-2}$, with solutions
denoted 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 $\delta t=10^{-2}$.

Figure 1: Comparison of the accuracy for the simulation methods in Section
II.3. We consider a chain of $6$ spins with dephasing and amplitude damping
noise in this section. The density operator $\rho$ is measured with 1-qubit
local observables $A:=\sigma_{2}^{y}$ and we denoted the measurement as
$y_{n}=\tr(A\rho(t_{n}))$. The blue solid line represents the ”exact” solution
generated with a very small step size.

Now we test the optimization algorithms in Section II.4 with the objective
function (5) computed from simulations of the Lindblad equation. We consider
observable time $t_{n}=0.1,0.2,\cdots,1$ with measurement interval $\Delta
t=0.1$ and simulation interval $\delta t=0.01$: $L=10$. In addition, we test
the algorithms by choosing $A^{(k)}$’s as the basis of all 1-local observables
or all 2-local observables generated by Pauli matrices. The optimization
results are displayed in LABEL:fig:exp30.3&4. To indicate the progress of the
optimization, we show the value of the objective function as well as the
relative error: $\norm{\bm{\theta}-\bm{\theta}^{*}}/\norm{\bm{\theta}^{*}}.$
We observe that the optimization using the Levenberg-Marquardt algorithm (26)
converges very quickly. In addition, the optimization problem (5) using all
2-local observables leads to a slightly faster parameter identification.

Figure 2: Reconstruction of the Lindbladians in Eqs. 46 and 47 by solving the
optimization problem (5). In terms of $y_{k,n}$ in (3), we consider two cases:
one with all 1-local observables ($N_{O}=19$), and the other with all 2-local
observables ($N_{O}=154$). The initial parameter $\bm{\theta}^{(0)}$ is
randomly selected with initial error
$\|\bm{\theta}^{(0)}-\bm{\theta}^{*}\|=0.3658$. The observable time is
$t_{n}=0.1,0.2,\cdots,1$ with measurement interval $\Delta t=0.1$ and
numerical simulation time difference $\delta t=0.01$.

We see from the previous experiment that having more observables can improve
the identification of the parameters. In the next experiment, we keep the
1-local observables, but increase the number of measurements in time instead.
Specifically, we set the observation time to $t=0.01,0.02,\cdots,1$, i.e.,
$\Delta t=\delta t=0.01$. From Fig. 3, we observe that these more frequent
measurements also improve the identification of the parameters slightly: The
optimization converged in fewer steps. The efficiency of the optimization can
again be attributed to the rapid convergence of the LM algorithm.

Figure 3: Reconstruction of the Lindbladians in Eqs. 46 and 47 by solving the
optimization problem (5). In terms of $y_{k,n}$ in (3), we consider all
1-local observables ($N_{O}=19$) and set the observation time to
$t=0.01,0.02,\cdots,1$, i.e. the measurement interval $\Delta t=0.01$
($N_{T}=100$). The simulation time interval $\delta t=0.01$. The initial
parameter $\bm{\theta}^{(0)}$ is the same as LABEL:fig:exp30.3&4.

Our last experiment is motivated by the study in [12], where the coefficients
in the Lindbladians are learned after the system reaches a steady state. To do
this, we choose the observation time at a later stage $t=4.1,4.2,\cdots,5$,
and a separate numerical test verified that this is when the system is
beginning to saturate (see LABEL:{fig:exp19Simu} as well). We see from Fig. 4
that the optimization is still able to identify all the parameters. In
particular, the parameters $\bm{\theta}_{H}$ do not converge monotonically.
But the last few iterations converged very fast.

Figure 4:  Reconstruction of the Lindbladians in Eqs. 46 and 47 by solving the
optimization problem (5). In terms of $y_{k,n}$ in (3), we consider all
1-local observables ($N_{O}=19$) and set the observation time to
$t=4.1,4.2,\cdots,5$, ($N_{T}=10$). And the simulation time interval $\delta
t=0.01$. The initial parameter $\bm{\theta}^{(0)}$ is the same as
LABEL:fig:exp30.3&4.

###  III.2 A Nonlinear parametric model

In this section, we consider a test problem similar to the ones considered in
Bairey et al. [12]. The model has the same parametric form for the Hamiltonian
terms (46). However, the jump operators are expanded as a linear combination
of local Pauli matrices,

$V_{j}=\sum_{\alpha=1}^{3}d_{j,\alpha}\sigma_{j}^{\alpha}=\sum_{\alpha=1}^{3}d_{j,\alpha}^{(1)}\sigma_{j}^{\alpha}+id_{j,\alpha}^{(2)}\sigma_{j}^{\alpha},\quad
j=1,\cdots,N_{V}.$ (49)

This parametric form is clearly more general. In particular, the Lindbladian
is a _nonlinear_ function of the parameters.

We first choose these dissipative parameters randomly according to a Gaussian
distribution, specifically,

$e_{j,\alpha},c_{j,\alpha,\beta}\sim\mathcal{N}(0,1),\quad\text{Re}(d_{j,\alpha}),\text{Im}(d_{j,\alpha})\sim\mathcal{N}(0,\frac{1}{2}).$
(50)

These parameters will then be fixed and considered to be exact. Subsequently,
they are used to generate the observation data $y_{k,n}^{*}$ in (3). We have
noticed that the parametric form in Eq. 50 is not unique due to the
simultaneous appearance of $V_{j}$ and $V_{j}^{\dagger}$ in the Lindbladian.
To eliminate the redundancy from a global phase of $d_{j}$ for each $j$, we
let

$d_{j,1}=d_{j,1}^{(1)},\quad j=1,\cdots,N_{V},$ (51)

which means that the imaginary part $d_{j,1}^{(2)}=0$. The parameter
$\bm{\theta}$ is composed by Hamiltonian part
$\bm{\theta}_{H}=\\{e_{j,\alpha},c_{j,\alpha,\beta}\\}$ as in Eq. 46 in the
previous test problem, as well as the dissipation part of the parameter is
$\bm{\theta}_{D}=\\{d_{j,\alpha}^{(1)},d_{j,\alpha}^{(2)}\\}$.

We consider the optimization problem (5) with values $y_{k,n}^{*}$ computed
from 1-local and 2-local observables at time $t=0.1,0.2,\cdots,1$. The
measurement time interval is thus $\Delta t=0.1$, while the simulation time
interval is much smaller: $\delta t=0.01$. In our numerical test, the randomly
chosen initial guess for the optimization has an error
$\|\bm{\theta}^{(0)}-\bm{\theta}^{*}\|=0.4529$. We then run the optimization
algorithm. The results are summarized in Fig. 5. Although the Lindblad
operator is a nonlinear function of the parameters, the optimization problem
still identifies the parameters rather quickly. Also observed is that with all
2-local observables, the optimization converges faster. One speculation is
that the constant in the condition (27) is larger for the choice with all
2-local observables.

To further test the effect of the choice of the observables, we choose 1-local
observables, by only keeping $\sigma_{j}^{x}$ and $\sigma_{j}^{y}$’s for each
spin, while other configurations of the test remain the same. In this case, we
observe that the optimization takes many more iterations to reach a plateau,
as shown in Fig. 6, and the optimization error is still quite large,
indicating that this collection of observables are inadequate for the
Lindbladian learning problem.

Figure 5: Reconstruction of the Lindbladians in Eqs. 46 and 49 by solving the
optimization problem (5). In terms of $y_{k,n}$ in (3), we consider two cases:
one with all 1-local observables ($N_{O}=19$), and the other with all 2-local
observables ($N_{O}=154$). The measurement time interval $\Delta t=0.1$ and
the step size for the simulation is $\delta t=0.01$. The initial parameter
$\bm{\theta}^{(0)}$ is randomly selected with initial error
$\|\bm{\theta}^{(0)}-\bm{\theta}^{*}\|=0.4529$.  Figure 6:  Reconstruction of
the Lindbladians in Eqs. 46 and 49 by solving the optimization problem (5). In
terms of $y_{k,n}$ in (3), we consider fewer 1-local observables
$\sigma^{x}_{j},\sigma^{y}_{j}$ ($N_{O}=12$) and measurement time
$t=0.1,0.2,\cdots,1$ ($N_{T}=10$) with simulation time interval $\delta
t=0.01$. The initial parameter $\bm{\theta}^{(0)}$ is the same as the previous
test.

##  IV 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 [12, 13]. 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.

Due to the close relation with quantum computing, one natural question is
whether we can use quantum algorithms to perform the simulation of the
Lindblad equation, especially when the dimension is too large for simulations
on classical computers. The idea of simulating Lindblad dynamics on quantum
computers has recently been explored by several groups [35, 36, 37, 38]. The
existing algorithms can be combined with gradient estimation algorithms to
update the parameters in a least squares algorithm, e.g., see [39, 40] for
examples of such hybrid schemes. Another alternative is to concatenate Eq. 19
and form a linear system of equations, which can be solved by quantum linear
solvers (QLS) [41, 42]. Due to the many technical elements involved in these
algorithms, we will present these methods in separate works.

## Acknowledgements

We thank Dr. Di Fang for providing some important references. This work is
supported by NSF Grant DMS-1819011 and DMS-2111221.

## References

  * da Silva  _et al._ [2011] M. P. da Silva, O. Landon-Cardinal, and D. Poulin, Practical characterization of quantum devices without tomography, Physical Review Letters 107, 210404 (2011).
  * Wang  _et al._ [2017a] J. Wang, S. Paesani, R. Santagati, S. Knauer, A. A. Gentile, N. Wiebe, M. Petruzzella, J. L. O’brien, J. G. Rarity, A. Laing, _et al._ , Experimental quantum Hamiltonian learning, Nature Physics 13, 551 (2017a).
  * Huang  _et al._ [2022] H.-Y. Huang, S. T. Flammia, and J. Preskill, Foundations for learning from noisy quantum experiments, arXiv preprint arXiv:2204.13691 (2022).
  * Granade  _et al._ [2012] C. E. Granade, C. Ferrie, N. Wiebe, and D. G. Cory, Robust online Hamiltonian learning, New Journal of Physics 14, 103013 (2012).
  * Rudinger and Joynt [2015] K. Rudinger and R. Joynt, Compressed sensing for Hamiltonian reconstruction, Physical Review A 92, 052322 (2015).
  * Wiebe  _et al._ [2015] N. Wiebe, C. Granade, and D. G. Cory, Quantum bootstrapping via compressed quantum Hamiltonian learning, New Journal of Physics 17, 022005 (2015).
  * Burgarth and Ajoy [2017] D. Burgarth and A. Ajoy, Evolution-free Hamiltonian parameter estimation through Zeeman markers, Physical Review Letters 119, 030402 (2017).
  * Sone and Cappellaro [2017] A. Sone and P. Cappellaro, Hamiltonian identifiability assisted by a single-probe measurement, Physical Review A 95, 022335 (2017).
  * Wang  _et al._ [2017b] Y. Wang, D. Dong, B. Qi, J. Zhang, I. R. Petersen, and H. Yonezawa, A quantum Hamiltonian identification algorithm: Computational complexity and error analysis, IEEE Transactions on Automatic Control 63, 1388 (2017b).
  * Zubida  _et al._ [2022] A. Zubida, E. Yitzhaki, N. Lindner, and E. Bairey, Optimal short-time measurements for Hamiltonian learning, Bulletin of the American Physical Society 67 (2022).
  * Huang  _et al._ [2023] H.-Y. Huang, Y. Tong, D. Fang, and Y. Su, Learning many-body Hamiltonians with Heisenberg-limited scaling, Physical Review Letters 130, 200403 (2023).
  * Bairey  _et al._ [2020] E. Bairey, C. Guo, D. Poletti, N. H. Lindner, and I. Arad, Learning the dynamics of open quantum systems from their steady states, New Journal of Physics 22, 032001 (2020).
  * Franca  _et al._ [2022] D. S. Franca, L. A. Markovich, V. Dobrovitski, A. H. Werner, and J. Borregaard, Efficient and robust estimation of many-qubit Hamiltonians, arXiv preprint arXiv:2205.09567 (2022).
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Communications in Mathematical Physics 48, 119 (1976).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of N-level systems, Journal of Mathematical Physics 17, 821 (1976).
  * Kloeden and Platen [2011] P. Kloeden and E. Platen,  _Numerical Solution of Stochastic Differential Equations_, Stochastic Modelling and Applied Probability (Springer Berlin Heidelberg, 2011).
  * Yamashita and Fukushima [2001] N. Yamashita and M. Fukushima, On the rate of convergence of the Levenberg-Marquardt method, in  _Topics in Numerical Analysis: With Special Emphasis on Nonlinear Problems_ (Springer, 2001) pp. 239–249.
  * Fan and Yuan [2005] J.-y. Fan and Y.-x. Yuan, On the quadratic convergence of the Levenberg-Marquardt method without nonsingularity assumption, Computing 74, 23 (2005).
  * Kelley [1999] C. T. Kelley,  _Iterative methods for optimization_ (SIAM, 1999).
  * Brunel [2008] N. J. Brunel, Parameter estimation of ODE’s via nonparametric estimators, Electronic Journal of Statistics 2, 1242 (2008).
  * Biele and D’Agosta [2012] R. Biele and R. D’Agosta, A stochastic approach to open quantum systems, Journal of Physics: Condensed Matter 24, 273201 (2012).
  * Di Ventra and D’Agosta [2007] M. Di Ventra and R. D’Agosta, Stochastic time-dependent current-density-functional theory, Physical Review Letters 98, 226403 (2007).
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Error mitigation for short-depth quantum circuits, Physical review letters 119, 180509 (2017).
  * Endo  _et al._ [2018] S. Endo, S. C. Benjamin, and Y. Li, Practical quantum error mitigation for near-future applications, Physical Review X 8, 031027 (2018).
  * Cai  _et al._ [2022] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean, and T. E. O’Brien, Quantum error mitigation, arXiv preprint arXiv:2210.00921 (2022).
  * Boulant  _et al._ [2003] N. Boulant, T. F. Havel, M. A. Pravia, and D. G. Cory, Robust method for estimating the lindblad operators of a dissipative quantum process from measurements of the density operator at multiple time points, Physical Review A 67, 042322 (2003).
  * Av  _et al._ [2020] E. B. Av, Y. Shapira, N. Akerman, and R. Ozeri, Direct reconstruction of the quantum-master-equation dynamics of a trapped-ion qubit, Physical Review A 101, 062305 (2020).
  * Strebel [2013] O. Strebel, A preprocessing method for parameter estimation in ordinary differential equations, Chaos, Solitons & Fractals 57, 93 (2013).
  * Bhaumik and Ghosal [2015] P. Bhaumik and S. Ghosal, Bayesian two-step estimation in differential equation models, Electronic Journal of Statistics 9, 3124 (2015).
  * Watrous [2018] J. Watrous,  _The theory of quantum information_ (Cambridge university press, 2018).
  * Cao and Lu [2021] Y. Cao and J. Lu, Structure-preserving numerical schemes for Lindblad equations, arXiv preprint arXiv:2103.01194 (2021).
  * Ruskai [1994] M. B. Ruskai, Beyond strong subadditivity? improved bounds on the contraction of generalized relative entropy, Reviews in Mathematical Physics 6, 1147 (1994).
  * Jin  _et al._ [2019] C. Jin, P. Netrapalli, R. Ge, S. M. Kakade, and M. I. Jordan, A short note on concentration inequalities for random vectors with subgaussian norm, arXiv preprint arXiv:1902.03736 (2019).
  * Kliesch  _et al._ [2011] M. Kliesch, T. Barthel, C. Gogolin, M. J. Kastoryano, and J. Eisert, Dissipative quantum Church-Turing theorem, Physical Review Letters 107, 10.1103/physrevlett.107.120501 (2011).
  * Childs and Li [2017] A. M. Childs and T. Li, Efficient simulation of sparse Markovian quantum dynamics, Quantum Information & Computation 17, 0901 (2017).
  * Cleve and Wang [2017] R. Cleve and C. Wang, Efficient quantum algorithms for simulating Lindblad evolution, in  _44th International Colloquium on Automata, Languages, and Programming, (ICALP 2017)_ (2017) pp. 17:1–17:14.
  * Xiantao Li [2022] C. W. Xiantao Li, Simulating Markovian open quantum systems using higher-order series expansion, arXiv preprint arXiv:2212.02051 (2022).
  * Gilyén  _et al._ [2019] A. Gilyén, S. Arunachalam, and N. Wiebe, Optimizing quantum optimization algorithms via faster quantum gradient computation, in  _Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms_ (SIAM, 2019) pp. 1425–1444.
  * Li and Wang [2023] X. Li and C. Wang, Efficient quantum algorithms for quantum optimal control, arXiv preprint arXiv:2304.02613 (2023).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for solving linear systems of equations, Physical Review Letters 103, 150502 (2009), arXiv: 0811.3171.
  * Childs  _et al._ [2017] A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM Journal on Computing 46, 1920 (2017).

##  Appendix A Appendix

###  A.1 Proof of Theorem 2

###### Proof.

The derivative of the objective function $\phi$ is of the form

$\partial_{\theta_{\alpha}}\phi(\bm{\theta})=\frac{1}{N_{O}N_{T}}\sum_{k=1}^{N_{O}}\sum_{n=1}^{N_{T}}(y_{k,n}-y_{k,n}^{*})\big{\langle}A^{(k)},\partial_{\theta_{\alpha}}\rho(\tau_{nL};\bm{\theta})\big{\rangle}$
(52)

Let $\chi_{nL}^{\alpha}=\partial_{\theta_{\alpha}}\rho(\tau_{nL};\bm{\theta})$
and $\rho_{nL}:=\rho(\tau_{nL};\theta)$. The iteration formula for the density
operator is of the Kraus form
$\rho_{m+1}=\mathcal{K}[\rho_{m}]:=\sum_{j}F_{j}\rho_{m}F_{j}^{{\dagger}},m=1,\cdots,nL$.
It reveals
$\chi_{m+1}^{\alpha}=\mathcal{K}[\chi_{m}^{\alpha}]+\partial_{\theta_{\alpha}}\mathcal{K}[\rho_{m}]$,
then

$\displaystyle\big{\langle}A^{(k)},\chi_{nL}^{\alpha}\big{\rangle}$
$\displaystyle=\big{\langle}A^{(k)},\mathcal{K}\chi_{nL-1}^{\alpha})+\big{\langle}A^{(k)},\partial_{\theta_{\alpha}}\mathcal{K}\rho_{nL-1}\big{\rangle}$
(53)
$\displaystyle\big{\langle}A^{(k)},\mathcal{K}\chi_{nL-1}^{\alpha}\big{\rangle}$
$\displaystyle=\big{\langle}A^{(k)},\mathcal{K}^{2}\chi_{nL-2}^{\alpha}\big{\rangle}+\big{\langle}A^{(k)},\mathcal{K}(\partial_{\theta_{\alpha}}\mathcal{K})\rho_{nL-2}\big{\rangle}$
$\displaystyle\vdots$
$\displaystyle\Rightarrow\big{\langle}A^{(k)},\chi_{nL}^{\alpha}\big{\rangle}$
$\displaystyle=\big{\langle}A^{(k)},\mathcal{K}^{nL}\chi_{0}^{\alpha}\big{\rangle}+\sum_{l=1}^{nL}\big{\langle}A^{(k)},\mathcal{K}^{l-1}(\partial_{\theta_{\alpha}}\mathcal{K})\rho_{nL-l}\big{\rangle}$
$\displaystyle=\sum_{l=1}^{nL}\big{\langle}A^{(k)},\mathcal{K}^{l-1}(\partial_{\theta_{\alpha}}\mathcal{K})\rho_{nL-l}\big{\rangle}$

as $\chi_{0}^{\alpha}=0$. By taking the adjoint operator,

$\displaystyle\partial_{\theta_{\alpha}}\phi(\bm{\theta})$
$\displaystyle=\frac{1}{N_{O}N_{T}}\sum_{k,n}(y_{k,n}-y_{k,n}^{*})\sum_{l=1}^{nL}\big{\langle}(\partial_{\theta_{\alpha}}\mathcal{K})^{*}(\mathcal{K}^{*})^{l-1}A^{(k)},\rho_{nL-l}\big{\rangle}$
(54)
$\displaystyle=\frac{1}{N_{O}N_{T}}\sum_{k,n}(y_{k,n}-y_{k,n}^{*})\sum_{l=1}^{nL}\big{\langle}(\partial_{\theta_{\alpha}}\mathcal{K})^{*}A^{(k)}_{nL-l},\rho_{nL-l}\big{\rangle}$

where $A^{(k)}_{nL-1}:=(\mathcal{K}^{*})^{l-1}[A^{(k)}]$ is a back-propagated
operator and the adjoint Kraus form

$(\partial_{\theta_{\alpha}}\mathcal{K})^{*}(\rho)=\sum_{j}(\partial_{\theta_{\alpha}}F_{j})^{{\dagger}}\rho
F_{j}+F_{j}^{{\dagger}}\rho(\partial_{\theta_{\alpha}}F_{j}),\quad\mathcal{K}^{*}(\rho)=\sum_{j}F_{j}^{{\dagger}}\rho
F_{j}$ (55)

for real parameters $\bm{\theta}$. ∎

###  A.2 Second Order Implicit Approximation Algorithm

Similar to Lemma 1, we start with the unraveled SDE,

$\mathrm{d}\ket{\psi}=(-iH_{S}\ket{\psi}-\frac{1}{2}\sum_{j=1}^{N_{V}}V_{j}^{{\dagger}}V_{j}\ket{\psi})\mathrm{d}t+\sum_{j=1}^{N_{V}}V_{j}\ket{\psi}\mathrm{d}W_{j;t}=G\ket{\psi}\mathrm{d}t+\sum_{j=1}^{N_{V}}V_{j}\ket{\psi}\mathrm{d}W_{j;t}$
(56)

According to the second order implicit weak scheme[17, Chapter 15], we can
write down a time-marching scheme for the wave function,

$\displaystyle\ket{\psi_{m+1}}$
$\displaystyle=\ket{\psi_{m}}+\frac{1}{2}(G\ket{\psi_{m+1}}+G\ket{\psi_{m}})\delta
t+\sum_{j=1}^{N_{V}}V_{j}\ket{\psi_{m}}\delta\hat{W}_{j}$ (57)
$\displaystyle+\frac{1}{2}\sum_{j=1}^{N_{V}}L^{0}V_{j}\ket{\psi_{m}}\delta\hat{W}_{j}\delta
t+\frac{1}{2}\sum_{j_{1},j_{2}=1}^{N_{V}}L^{j_{1}}V_{j_{2}}\ket{\psi_{m}}(\delta\hat{W}_{j_{1}}\delta\hat{W}_{j_{2}}+U_{j_{1},j_{2}}).$

In this expression, $\delta\hat{W}$ is an approximation of an increment of the
Brownian motion. It was suggested in [17] to sample with the three-point
distribution,

$P(\delta\hat{W}=\pm\sqrt{3\Delta})=\frac{1}{6},P(\delta\hat{W}=0)=\frac{2}{3}$
(58)

Meanwhile, $U_{j_{1},j_{2}}$’s are independent two-point distributed random
variables with

$\displaystyle P(U_{j_{1},j_{2}}=\pm\Delta)=\frac{1}{2},\forall
j_{2}=1,\cdots,j_{1}-1$ (59) $\displaystyle U_{j_{1},j_{1}}=-\Delta$
$\displaystyle U_{j_{1},j_{2}}=-U_{j_{2},j_{1}},\text{ for
}j_{2}=j_{1}+1,\cdots,N_{V},j_{1}=1,\cdots,N_{V}.$

In Eq. 57, the diffusion operator $L^{j_{1}}$’s are defined as follows.

$\displaystyle L^{0}V_{j}\ket{\psi}$
$\displaystyle=\sum_{k=1}^{d}(G\psi)^{k}\frac{\partial}{\partial\psi^{k}}V_{j}\ket{\psi}=\sum_{k=1}^{d}(G\psi)^{k}V_{j}(:,k)=V_{j}G\ket{\psi}$
(60) $\displaystyle L^{j_{1}}V_{j_{2}}\ket{\psi}$
$\displaystyle=\sum_{k=1}^{d}(V_{j_{1}}\psi)^{k}\frac{\partial}{\partial\psi^{k}}V_{j_{2}}\ket{\psi}=\sum_{k=1}^{d}(V_{j_{1}}\psi)^{k}V_{j_{2}}(:,k)=V_{j_{2}}V_{j_{1}}\ket{\psi}$

Here the notation $(:,k)$ indicates the $k$th column of the matrix. Thus, the
equation 19 or 22 can be rewritten as

$\displaystyle\ket{\psi_{m+1}}$ $\displaystyle=(I-\frac{1}{2}G\delta
t)^{-1}(I+\frac{1}{2}G\delta
t)\ket{\psi_{m}}+\sum_{j=1}^{N_{V}}(I-\frac{1}{2}G\delta
t)^{-1}(V_{j}+\frac{1}{2}V_{j}G\delta t)\ket{\psi_{m}}\delta\hat{W}_{j}$ (61)
$\displaystyle+\frac{1}{2}(I-\frac{1}{2}G\delta
t)^{-1}\sum_{j_{1},j_{2}=1}^{N_{V}}V_{j_{2}}V_{j_{1}}\ket{\psi_{m}}(\delta\hat{W}_{j_{1}}\delta\hat{W}_{j_{2}}+U_{j_{1},j_{2}}).$

By taking expectations, this time-marching scheme induces an iteration formula
for the density operator $\rho_{m}=\mathbb{E}[\ket{\psi_{m}}\bra{\psi_{m}}]$,
as follows.

$\displaystyle\rho_{m+1}$ $\displaystyle=(I-\frac{1}{2}G\delta
t)^{-1}(I+\frac{1}{2}G\delta t)\rho_{m}(I+\frac{1}{2}G\delta
t)(I-\frac{1}{2}G\delta t)^{-1}$ (62) $\displaystyle+(I-\frac{1}{2}G\delta
t)^{-1}\sum_{j=1}^{N_{V}}V_{j}(I+\frac{1}{2}G\delta
t)\rho_{m}(I+\frac{1}{2}G\delta t)V_{j}^{{\dagger}}(I-\frac{1}{2}G\delta
t)^{-1}\delta t$ $\displaystyle+\frac{1}{2}(I-\frac{1}{2}G\delta
t)^{-1}\sum_{j_{1},j_{2}=1}^{N_{V}}V_{j_{1}}V_{j_{2}}\rho_{m}V_{j_{2}{\dagger}}V_{j_{1}{\dagger}}(I-\frac{1}{2}G\delta
t)^{-1}(\delta t)^{2}$
$\displaystyle=\sum_{j=0}^{N_{V}}F_{j}\rho_{m}F_{j}^{{\dagger}}+\sum_{j_{1},j_{2}=1}^{N_{V}}A_{j_{1},j_{2}}\rho_{m}A_{j_{1},j_{2}}^{{\dagger}}$
$\displaystyle=\sum_{j=0}^{N_{V}^{2}+N_{V}}F_{j}\rho_{m}F_{j}^{{\dagger}}=\mathcal{K}[\rho_{m}]$

where the Kraus operator

$\displaystyle F_{0}$ $\displaystyle=(I-\frac{1}{2}G\delta
t)^{-1}(I+\frac{1}{2}G\delta t)$ (63) $\displaystyle F_{j}$
$\displaystyle=(I-\frac{1}{2}G\delta t)^{-1}V_{j}(I+\frac{1}{2}G\delta
t)\sqrt{\delta t},j=1,\cdots,N_{V}$ $\displaystyle A_{j_{1},j_{2}}$
$\displaystyle=\frac{1}{\sqrt{2}}(I-\frac{1}{2}G\delta
t)^{-1}V_{j_{1}}V_{j_{2}}\delta t,j_{1},j_{2}=1,\cdots,N_{V}$ $\displaystyle
F_{j_{1}+N_{V}j_{2}}$
$\displaystyle=A_{j_{1},j_{2}},j_{1},j_{2}=1,\cdots,N_{V}$

To show that the method has second order accuracy, we first expand the
solution of the Lindblad equation (2),

$\displaystyle\rho(\tau_{m+1})=\rho(\tau_{m})+\delta
t\mathcal{L}\rho(\tau_{m})+\frac{1}{2}\delta
t^{2}\mathcal{L}^{2}\rho(\tau_{m})+\mathcal{O}(\delta t^{3}).$ (64)

In particular, we have,

$\mathcal{L}\rho=-i[H,\rho]+\sum_{j=1}^{N_{V}}\left(V_{j}\rho
V_{j}^{{\dagger}}-\frac{1}{2}\\{V_{j}^{{\dagger}}V_{j},\rho\\}\right).$

We now show that the method (19) has a consistent expansion at (64) up to
$\mathcal{O}(\delta t^{2}).$ Toward this end, we expand the Kraus operators,

$\displaystyle F_{0}=I+G\delta t+\frac{1}{2}G^{2}\delta
t^{2}+\mathcal{O}(\delta t^{3}),$ $\displaystyle F_{j}=V_{j}\sqrt{\delta
t}+\frac{1}{2}\\{G,V_{j}\\}\delta t\sqrt{\delta
t}+\frac{1}{2}\\{G,GV_{j}\\}\delta t^{2}\sqrt{\delta t}+\mathcal{O}(\delta
t^{7/2}),\qquad j=1,\cdots,N_{V}$ $\displaystyle
F_{j_{1}+N_{V}j_{2}}=\frac{1}{\sqrt{2}}V_{j_{1}}V_{j_{2}}\delta
t+\mathcal{O}(\delta t^{3/2}),\qquad j_{1},j_{2}=1,\cdots,N_{V}$

By substituting these expansions of the jump operators, and only keeping terms
of order up to $\mathcal{O}(\delta t^{2})$, we find that,

$\displaystyle\rho_{m+1}=$ $\displaystyle\rho_{m}+\delta
t(G\rho_{m}+\rho_{m}G^{\dagger})+\frac{1}{2}\delta
t^{2}(G^{2}\rho_{m}+2G\rho_{m}G^{\dagger}+\rho_{m}{G^{\dagger}}^{2})$
$\displaystyle+\delta
t\sum_{j=1}^{N_{V}}V_{j}\rho_{m}V_{j}^{\dagger}+\frac{1}{2}\delta
t^{2}\sum_{j=1}^{N_{V}}\left(V_{j}\rho_{m}\\{G,V_{j}\\}^{\dagger}+\\{G,V_{j}\\}\rho_{m}V_{j}^{\dagger}\right)$
$\displaystyle+\frac{1}{2}\delta
t^{2}\sum_{j_{1},j_{2}=1}^{N_{V}}V_{j_{1}}V_{j_{2}}\rho_{m}V_{j_{2}}^{\dagger}V_{j_{1}}^{\dagger}+\mathcal{O}(\delta
t^{3}).$

In light of Eq. 13, the $\mathcal{O}(\delta t)$ terms are consistent with that
in (64). With lengthy calculations, we can also see that the
$\mathcal{O}(\delta t^{2})$ terms are also consistent. From the calculation of
the consistency of $\mathcal{O}(\delta t)$, we know that
$G\rho_{m}+\rho_{m}G^{{\dagger}}=\mathcal{L}\rho_{m}-\sum_{j}V_{j}\rho_{m}V_{j}^{{\dagger}}$.
The detailed calculation can be summarized as follows:

$\displaystyle
G^{2}\rho_{m}+2G\rho_{m}G^{\dagger}+\rho_{m}{G^{\dagger}}^{2}=G(G\rho_{m}+\rho_{m}G^{{\dagger}})+(G\rho_{m}+\rho_{m}G^{{\dagger}})G^{{\dagger}}$
$\displaystyle=$
$\displaystyle\mathcal{L}^{2}\rho_{m}-\sum_{j=1}^{N_{V}}\mathcal{L}(V_{j}\rho_{m}V_{j}^{{\dagger}})-\sum_{j=1}^{N_{V}}V_{j}\mathcal{L}\rho_{m}V_{j}^{{\dagger}}+\sum_{j_{1},j_{2}=1}^{N_{V}}V_{j_{2}}V_{j_{1}}\rho_{m}V_{j_{1}}V_{j_{2}}$
$\displaystyle\sum_{j=1}^{N_{V}}\left(V_{j}\rho_{m}\\{G,V_{j}\\}^{\dagger}+\\{G,V_{j}\\}\rho_{m}V_{j}^{\dagger}\right)=\sum_{j}(V_{j}\rho_{m}V_{j}^{{\dagger}}G^{{\dagger}}+GV_{j}\rho_{m}V_{j}^{{\dagger}})+V_{j}(\rho_{m}G^{{\dagger}}+G\rho_{m})V_{j}^{{\dagger}}$
$\displaystyle=$
$\displaystyle\sum_{j}\mathcal{L}(V_{j}\rho_{m}V_{j}^{{\dagger}})-2\sum_{j_{1},j_{2}}V_{j_{2}}V_{j_{1}}\rho_{m}V_{j_{1}}^{{\dagger}}V_{j_{2}}^{{\dagger}}+\sum_{j}V_{j}\mathcal{L}\rho_{m}V_{j}^{{\dagger}}$
$\displaystyle\Rightarrow$ $\displaystyle
G^{2}\rho_{m}+2G\rho_{m}G^{\dagger}+\rho_{m}{G^{\dagger}}^{2}+\sum_{j=1}^{N_{V}}\left(V_{j}\rho_{m}\\{G,V_{j}\\}^{\dagger}+\\{G,V_{j}\\}\rho_{m}V_{j}^{\dagger}\right)+\sum_{j_{1},j_{2}}V_{j_{2}}V_{j_{1}}\rho_{m}V_{j_{1}}^{{\dagger}}V_{j_{2}}^{{\dagger}}=\mathcal{L}^{2}\rho_{m}.$

This completes the proof.