11institutetext: 1Institute of Aerodynamics and Flow Technology, German
Aerospace Center, Germany, 11email<EMAIL_ADDRESS>
2 TUM School of Engineering and Design, Technical University of Munich,
Germany
3 Institute for Hydromechanics, Karlsruhe Institute of Technology, Germany

# Direct numerical simulation of turbulent open channel flow: Streamwise
turbulence intensity scaling and its relation to large-scale coherent motions

Christian Bauer1 Yoshiyuki Sakai2 and Markus Uhlmann3

###### Abstract

We conducted direct numerical simulations of turbulent open channel flow (OCF)
and closed channel flow (CCF) of friction Reynolds numbers up to
$\mathrm{Re}_{\tau}\approx 900$ in large computational domains up to
$L_{x}\times L_{z}=12\pi h\times 4\pi h$ to analyse the Reynolds number
scaling of turbulence intensities. Unlike CCF, our data suggests that the
streamwise turbulence intensity in OCF scales with the bulk velocity for
$\mathrm{Re}_{\tau}\gtrsim 400$. The additional streamwise kinetic energy in
OCF with respect to CCF is provided by larger and more intense very-large-
scale motions in the former type of flow. Therefore, compared to CCF, larger
computational domains of $L_{x}\times L_{z}=12\pi h\times 4\pi h$ are required
to faithfully capture very-large-scale motions in OCF—and observe the reported
scaling. OCF and CCF turbulence statistics data sets are available at
https://doi.org/10.4121/88678f02-2a34-4452-8534-6361fc34d06b

##  1 Introduction

Plane Poiseuille flow, also known as closed channel flow (CCF), is one of the
most studied canonical flows by means of direct numerical simulations (DNSs,
Kim1987 ; Moser1999 ; Hoyas2006 ; Lee2015 ; Oberlack2022 ). The numerical
domain is defined by doubly-periodic boundary condition in the stream- and
spanwise directions, and impermeable no-slip walls at the bottom and the top.
Conversely, less attention has been paid to open channel flow (OCF), where one
of the no-slip walls is replaced by a free-slip plane, despite its even more
direct relevance in environmental flows. Most of the early OCF DNSs were
restricted to small friction Reynolds numbers $Re_{\tau}=u_{\tau}h/\nu$, based
on friction velocity $u_{\tau}$, channel height $h$ and kinematic viscosity
$\nu$, and/or small computational domains. Recently, Yao et al. Yao2022
compared CCF and OCF by means of DNS data in computational domains of
$L_{x}\times L_{z}=8\pi h\times 4\pi h$ for Reynolds numbers up to
$Re_{\tau}=2000$. In agreement with experiments Duan2020 ; Duan2021  and DNSs
Bauer2015 , they reported that so-called very-large-scale motions
(VLSMs)—which are a common feature of wall-bounded turbulent flows—are more
energetic and larger in OCF than in CCF. Generally, VLSMs become more
energetic with increasing Reynolds number Bauer2015 ; Duan2021 ; Yao2022 ,
leading to the failure of wall scaling of the streamwise intensity in CCF
Hoyas2006 . While for CCF the streamwise turbulence intensity scales neither
with the friction velocity nor with the bulk velocity, Afzal et al. Afzal2009
suggested that the streamwise turbulence intensity in OCF scales with the bulk
velocity based on an experimental data set. Other experimental studies, such
as Duan et al. Duan2020 ; Duan2021 , on the contrary, show the—relatively
scattered—streamwise turbulence intensity scaled in wall units together with
universal scaling functions proposed by Nezu and Rodi Nezu1986 . Thus,
additional highly accurate OCF data sets are required to answer the question
whether the streamwise turbulence intensity in OCF scales in bulk units.

Since the streamwise length of VLSMs is an order of magnitude larger than the
channel half width, exceptionally large domains are required to faithfully
capture VLSMs in DNS Lozano2014 ; Feldmann2018 . Moreover, Pinelli et al.
Pinelli2022  found that the small scales—responsible for the mass transfer
across the free surface in OCF—are spatially organised by VLSMs, which
requires an additional grid refinement towards the free surface. In the
present study, we investigate the scaling of the streamwise turbulence
intensity in OCF by means of DNSs. For this purpose, we generate new OCF and
CCF DNS datasets at comparable Reynolds numbers up to $\mathrm{Re}_{\tau}=900$
in identical computational domains of $L_{x}\times L_{z}=12\pi h\times 4\pi
h$, with grid refinements towards no-slip and free-slip walls. By removing
Reynolds number or domain size discrepancies between OCF and CCF, we are able
to quantify the impact of the free-slip boundary condition on turbulence one-
point statistics in the cleanest possible manner.

##  2 Methodology

We solve the incompressible Navier-Stokes equations in their wall-normal
velocity/vorticity formulation with a pseudo-spectral representation of the
flow variables Kim1987 . We apply periodic boundary conditions in the stream-
and the spanwise directions, the smooth no-slip wall at the bottom and the
free-slip boundary condition at the top of the channel. Although the numerical
method has shown its validity in numerous CCF simulations Kim1987 ; Moser1999
; DelAlamo2003 , the introduction of a free-slip surface has to be regarded
separately. A Chebyshev-Gauss-Lobatto grid, which is adopted in the surface-
normal direction provides a grid refinement both towards the wall and the free
surface. We set up six OCF and four CCF cases, where the Reynolds number and
the computational domain length are varied (cf. table 1).

Table 1: Turbulent channel flow simulation cases. $N_{x}$, $N_{y}$ and $N_{z}$ are the number of grid points in the streamwise, wall-normal and spanwise direction; $\Delta x^{+}$, the streamwise grid spacing; $\Delta z^{+}$, the spanwise grid spacing; $\Delta y^{+}_{max}$, the maximum grid spacing in the wall-normal direction. $\Delta T$ is the temporal averaging interval for turbulent statistics. Prefix “O” indicates OCF cases, “C” CCF cases. Case | O200 | O400 | O600 | O900L4 | O900L8 | O900L12 | C200 | C400 | C600 | C900
---|---|---|---|---|---|---|---|---|---|---
$Re_{\tau}$ | $200$ | $399$ | $596$ | $899$ | $895$ | $894$ | $200$ | $397$ | $593$ | $890$
$Re_{b}$ | $3170$ | $6969$ | $11\,047$ | $17\,512$ | $17\,512$ | $17\,512$ | $3170$ | $6969$ | $11\,047$ | $17\,512$
$L_{x}/h$ | $12\pi$ | $12\pi$ | $12\pi$ | $4\pi$ | $8\pi$ | $12\pi$ | $12\pi$ | $12\pi$ | $12\pi$ | $12\pi$
$L_{z}/h$ | $4\pi$ | $4\pi$ | $4\pi$ | $2\pi$ | $4\pi$ | $4\pi$ | $4\pi$ | $4\pi$ | $4\pi$ | $4\pi$
$N_{x}$ | $768$ | $1536$ | $1536$ | $1536$ | $2048$ | $3072$ | $768$ | $1536$ | $2048$ | $3072$
$N_{y}$ | $129$ | $193$ | $257$ | $385$ | $385$ | $385$ | $129$ | $193$ | $257$ | $385$
$N_{z}$ | $512$ | $1024$ | $1536$ | $1024$ | $2048$ | $2048$ | $512$ | $1024$ | $1536$ | $2048$
$\Delta x^{+}$ | $9.8$ | $9.8$ | $14.6$ | $7.4$ | $11.0$ | $11.0$ | $9.8$ | $9.8$ | $14.6$ | $11.0$
$\Delta z^{+}$ | $4.9$ | $4.9$ | $4.8$ | $5.5$ | $5.5$ | $5.5$ | $4.9$ | $4.9$ | $4.8$ | $5.5$
$\Delta y^{+}_{max}$ | $2.5$ | $3.3$ | $3.7$ | $3.7$ | $3.7$ | $3.7$ | $4.9$ | $6.5$ | $7.3$ | $7.3$
$\Delta Tu_{b}/h$ | $8660$ | $1925$ | $1460$ | $417$ | $570$ | $1054$ | $8600$ | $3260$ | $1757$ | $1013$
$\Delta Tu_{\tau}^{2}/\nu$ | $109\,700$ | $43\,900$ | $46\,980$ | $19\,080$ | $26\,100$ | $44\,700$ | $108\,200$ | $73\,730$ | $55\,960$ | $45\,820$

Hereinafter, the velocity fluctuation in $i$-direction is defined as
$u^{\prime}_{i}=u_{i}-\langle u_{i}\rangle$, where angular brackets denote
averaging over both homogeneous directions and time. Normalisation in viscous
or “wall units” is indicated by the $+$ superscript. Corresponding scales are
$u_{\tau}$ and the viscous length scale $\delta_{\nu}=\nu/u_{\tau}$.
Characteristic scales in “bulk units”, on the other hand, are the bulk
velocity $u_{b}$ and the channel height $h$, which forms the bulk Reynolds
number $Re_{b}=u_{b}h/\nu$.

##  3 Computational domain size

In the following, the effect of the constraints upon large-scale coherent
motions in the vicinity of the free surface due to the finite size of the
computational domain is investigated. Figure 1 displays instantaneous
realisations of the streamwise fluctuating velocity for the three open channel
simulation cases with different computational domain size O900L4, O900L8, and
O900L12 in the vicinity of the free-slip boundary ($y/h\approx 0.95$),
together with closed channel data for comparison.

Figure 1: Instantaneous streamwise velocity fluctuations $u^{\prime+}$ in a
wall-parallel plane directly below the free surface in OCF (a,b,c) and the
vicinity of the channel centreline in CCF (d, $y/h\approx 0.95$). (a) O900L4;
(b) O900L8; (c) O900L12; (d) C900.

A coherent structure with a streamwise extent comparable to the domain length
$l_{x}\approx 12\pi h$ and a spanwise extent of $l_{z}\approx 2h$ is clearly
visible in figure 1(c) for O900L12 at $z/h\approx 9$, but not for the smaller
boxes (a,b), where the flow domain interferes with the natural extensions of
such a structure. Moreover, the structures in the vicinity of the free-slip
boundary in OCF appear to be more intense and larger than the structures near
the centreline of CCF (figure 1d). Thus, larger computational domains than in
CCF are required to fully capture these scales. Figure 2 displays one-
dimensional pre-multiplied energy spectra of the streamwise velocity
fluctuations for $\mathrm{Re}_{\tau}\approx 900$ and different computational
domain length as a function of the wavelength and the wall distance.

Figure 2: Contour lines of the one-dimensional pre-multiplied energy spectra
of the streamwise velocity fluctuations as a function of the wavelength and
the wall distance. (a,c,e) Streamwise spectra
$\kappa_{x}\phi^{+}_{uu}(\lambda_{x})$; (b,d,f) spanwise spectra
$\kappa_{z}\phi^{+}_{uu}(\lambda_{z})$. (a,b) O900L4; (c,d) O900L8; (e,f)
O900L12. Spectra are normalised by their peak values. Values range from 0
(white) to 0.5 (orange) with an increment of 0.05.

Although at first glance both the streamwise and the spanwise pre-multiplied
energy spectra decay at large wavelengths for O900L4 (figure 2a,b), in fact
larger box sizes are required to physically capture the largest scales that
only appear in the larger domains. This can be seen in the streamwise spectra
for O900L8 (figure 2c), where a spectral pile-up at large wavelength occurs.
For O900L12 (figure 2e,f), on the other hand, both the streamwise and the
spanwise spectra decay at large wavelengths again. Thus, a domain size of
$L_{x}\times L_{z}=12\pi h\times 4\pi h$ appears to be sufficiently large to
capture a relevant amount of VLSMs. Moreover, a comparison of the spanwise
spectra (figure 2b,d,f) indicates an artificially large accumulation of energy
at spanwise wavelengths around $\lambda_{z}/h\approx 2$ in the smallest box
($L_{x}=4\pi h$, $L_{z}=2\pi h$). In summary, as in CCF, VLSMs appear in the
bulk flow, penetrating into the buffer layer and therefore through these
motions the whole flow domain is influenced by the free-surface layer. Due to
their large spatial extent, larger computational domains than in CCF are
needed to fully capture these scales. A domain size of $L_{x}=12\pi h$ was
found to be sufficient at $\mathrm{Re}_{\tau}=900$.

##  4 Turbulence intensity scaling

Let us now focus on the variation with Reynolds number of the turbulence
intensities normalised with the bulk velocity $u_{b}$. Both open and closed
channel turbulence intensities obtained from DNSs with the largest
computational domain ($L_{x}\times L_{z}=12\pi h\times 4\pi h$) are presented
in figure 3.

Figure 3: Turbulence intensities normalised by $u_{b}$ as function of the
distance from the wall $y/h$ for OCF (left) and CCF (right):  ,
$\mathrm{Re}_{\tau}\approx 200$;  , $\mathrm{Re}_{\tau}\approx 400$;  ,
$\mathrm{Re}_{\tau}\approx 600$;  , $\mathrm{Re}_{\tau}\approx 900$. The
insets show a zoom for the upper channel half. Left: The symbols ( ) indicate
a profile from OCF measurements at $\mathrm{Re}_{\tau}=2407$ Duan2020 . Right:
The gray lines indicate CCF DNS data at $\mathrm{Re}_{\tau}=2003$ ( ,
Hoyas2006 ), $\mathrm{Re}_{\tau}=5186$ ( , Lee2015 ),
$\mathrm{Re}_{\tau}=10049$ ( , Oberlack2022 ).

While the spanwise and wall-normal contributions to the turbulence intensity
scale well with the friction velocity $u_{\tau}$ in the bulk region of the
flow for $Re_{\tau}\geq 400$ for both flow configurations (figure omitted),
the streamwise contribution differs. As figure 3(a) indicates, the streamwise
turbulence intensity in OCF appears to scale with the bulk velocity $u_{b}$
for $Re_{\tau}\geq 400$ and $y/h\gtrsim 0.3$, whereas in the closed channel
case it neither scales with $u_{\tau}$ (figure omitted) nor with $u_{b}$
(figure 3b). This scaling property of OCF is in qualitative agreement with the
experimental study of OCF by Afzal2009 , who observed the scaling of the outer
flow streamwise turbulence intensity with $u_{b}$ for $1260\leq Re_{\tau}\leq
5280$. Besides, the streamwise turbulence intensity profile obtained from
recent OCF measurements by Duan2020 ; Duan2021  at $\mathrm{Re}_{\tau}=2407$
normalised with $u_{b}$ collapses well with our data (figure 3a inset). The
above observation differs significantly from CCF, where the spectral signature
of VLSMs first appears at $\mathrm{Re}_{\tau}\approx 600$ (figure omitted) and
the contribution of VLSMs to the streamwise turbulence intensity is much lower
than in OCF Duan2020 . Even at $\mathrm{Re}_{\tau}=10049$, which is currently
the highest Reynolds number ever been achieved by a CCF DNS Oberlack2022 , the
bulk scaling of the streamwise turbulence intensity has not been observed.
Hence, it is remarkable that OCF streamwise turbulence intensity scales with
$u_{b}$ at a Reynolds number as low as $\mathrm{Re}_{\tau}=400$, if a
sufficiently large numerical domain is used.

##  5 Conclusion

We performed DNSs of OCF and CCF in computational domains up to $L_{x}\times
L_{z}=12\pi h\times 4\pi h$ and Reynolds numbers up to
$\mathrm{Re}_{\tau}=900$. The enhancement of VLSMs in OCF leads to more
stringent requirements on the domain size than in the closed channel
counterpart. Since VLSMs appear at lower Re in OCF, a domain size $L_{x}\geq
12\pi h$ is already required at the moderate value $Re_{\tau}=400$. Unlike
CCF, the streamwise turbulence intensity in OCF scales with the bulk velocity
$u_{b}$ for $Re_{\tau}\geq 400$ and $y/h\gtrsim 0.3$, which is attributed to
the larger and more intense VLSMs in OCF. Although it was suggested by the
experiment of Afzal et al. Afzal2009 , this is the first time the bulk scaling
of the streamwise turbulence intensity has been observed in DNS data. In
addition, only with $L_{x}=12\pi h$ the bulk-scaling profile appears in the
$Re_{\tau}=400$ case (figure omitted), which indicates that VLSMs in OCF
already play a significant role at this Reynolds number. Thus, it is of
critical importance to use sufficiently large numerical domains. Finally, both
OCF and CCF turbulence statistics data sets are available online BauerData2023
.

###### Acknowledgements.

Many fruitful discussions with Genta Kawahara throughout this work are
gratefully acknowledged. The simulations were carried out at UC2 of SCC
Karlsruhe and at DLR’s CARA cluster. The computer resources and assistance
provided by these centres are gratefully acknowledged. YS and MU acknowledge
funding by DFG through grant UH 242/3-1.

## References

  * (1) Bauer C., Sakai Y., and Uhlmann M. (2023), _”Data underlying the publication: Direct numerical simulation of turbulent open channel flow”_ , 4TU.ResearchData, doi: 10.4121/88678F02-2A34-4452-8534-6361FC34D06B.V1
  * (2) Kim J., Moin P. and Moser R.D. (1987), _J. Fluid Mech., vol. 177, 133–166_
  * (3) Moser R.D., Kim J., and Mansour N.N. (1999), _Phys. Fluids, vol. 11, 943–945_
  * (4) Hoyas S. and Jiménez J. (2006), _Phys. Fluids, vol. 18, 11702_
  * (5) Lee M. and Moser R. D. (2015), _J. Fluid Mech., vol. 774, 395–415_
  * (6) Oberlack M., Hoyas S., Kraheberger S. V., Alcántara-Ávila F. and Laux J. (2022), _Phys. Rev. Lett., vol. 128, 024502_
  * (7) Yao J., Chen X. and Hussain F. (2022), _J. Fluid Mech., vol. 953, A19_
  * (8) Duan Y., Chen Q., Li D. and Zhong Q. (2020), _J. Fluid Mech., vol. 892, A3_
  * (9) Duan Y., Zhong Q., Wang G., Zhang P. and Li D. (2021), _J. Fluid Mech., vol. 918, A40_
  * (10) Bauer C. (2015), _Master’s thesis, KIT, Karlsruhe_
  * (11) Afzal B., Faruque M. A. and Balachandar, R. (2009), _J. Hydraul. Res., vol. 47(1), 66–81_ ,
  * (12) Nezu I. and Rodi W. (1986), _J. Hydraul. Eng., vol. 112,335–355_ ,
  * (13) Lozano-Durán A. and Jiménez J. (2014), _Phys. Fluids, vol. 26, 011702_ ,
  * (14) Feldmann D., Bauer C. and Wagner C. (2018), _J. Turbul., vol. 19(3), 274–295_
  * (15) Pinelli M., Herlina H., Wissink J.G. and Uhlmann, M. (2022), _J. Fluid Mech., vol. 933, A49_
  * (16) del Álamo J. C. and Jiménez J. (2003) _Phys. Fluids, vol. 15(6), L41–L44_


# Hamiltonian Learning using Machine Learning Models Trained with Continuous
Measurements

Kris Tucker<EMAIL_ADDRESS>Department of Applied Math,
University of Colorado Boulder, USA  Amit Kiran Rege<EMAIL_ADDRESS>Department of Computer Science, University of Colorado Boulder, USA  Conor
Smith  Center for Quantum Information and Control, University of New Mexico,
University of New Mexico, Albuquerque, NM, USA  Department of Electrical and
Computer Engineering, University of New Mexico, Albuquerque, NM, USA  Claire
Monteleoni  Department of Computer Science, University of Colorado Boulder,
USA  INRIA Centre de Recherche de Paris, France  Tameem Albash  Center for
Computing Research, Sandia National Laboratories, Albuquerque, New Mexico
87185, USA

###### Abstract

We build upon recent work on using Machine Learning models to estimate
Hamiltonian parameters using continuous weak measurement of qubits as input.
We consider two settings for the training of our model: (1) supervised
learning where the weak measurement training record can be labeled with known
Hamiltonian parameters, and (2) unsupervised learning where no labels are
available. The first has the advantage of not requiring an explicit
representation of the quantum state, thus potentially scaling very favorably
to larger number of qubits. The second requires the implementation of a
physical model to map the Hamiltonian parameters to a measurement record,
which we implement using an integrator of the physical model with a recurrent
neural network to provide a model-free correction at every time step to
account for small effects not captured by the physical model. We test our
construction on a system of two qubits and demonstrate accurate prediction of
multiple physical parameters in both the supervised and unsupervised context.
We demonstrate that the model benefits from larger training sets establishing
that it is in fact “learning,” and we show robustness to errors in the assumed
physical model by achieving accurate parameter estimation in the presence of
unanticipated single particle relaxation.

††preprint: SAND2024-04128O

##  I Introduction

As the scale, complexity, and availability of quantum devices continues to
grow over the coming decades [1], the ability to accurately characterize
device parameters, especially unwanted effects, will become ever-more
important. For example, non-local errors such as crosstalk in large systems
[2, 3] can be problematic for quantum error correction methods [4, 5, 6, 7,
8]. Detection of these unwanted effects is a critical first step in mitigating
them, with efforts to do so complicated by the exponential scaling of the
dimension of the Hilbert space with system size and the challenges of modeling
increasingly complex systems.

Recent advances in the application of machine learning (ML) tools to quantum
systems have shown promise in overcoming some of these difficulties. Scalable
ML models have been used to represent quantum states for state tomography [9,
10, 11, 12, 13, 14] and evolution [14, 15], for learning unknown dynamics [16,
17, 18], as well as for learning device characterization from measurements
[19, 20, 21, 22, 23]. Increasingly, a priori knowledge of a physical system
has been combined with ML models to improve accuracy and interpretability
[24]. The model-free nature of many ML solutions has made them scalable and
robust to many of the common pitfalls of physical models, such as non-
Markovian dynamics [19, 20, 21, 15]. There are trade-offs, however, as more
abstract models are generally less interpretable, which limits the physical
insights that can be gleaned from them. Finding a balance between these two
competing properties, representability and interpretability, thus becomes an
important challenge.

In this work, we build on recent advances in device characterization applying
ML models to the continuous measurement of qubits [20, 21, 19]. In contrast
with the parameter estimation approach of Ref. [19], where a stochastic master
equation is used as a trainable model, our ML model learns a direct map from
continuous measurement inputs to the system parameters of interest. In the
case where the measurement records are associated with known device parameters
(the case of supervised learning), this approach has two advantages. First,
the model is completely independent of any state representation. The
measurements need only provide enough information to estimate the parameters
of interest and not the complete quantum state, potentially freeing it from
the curse of dimensionality that comes from the exponential scaling of the
dimension of the Hilbert space and thus making it more scalable.

Second, once the ML model is trained, parameter estimation can be performed
quickly even for systems with parameters not seen during training. It can also
benefit from having been provided with many examples of noisy measurement
records associated with known true values in the training set. This makes it
possible for the model to learn to distinguish between those features of
measurement records that vary with parameters and those that are uncorrelated
noise, thus potentially requiring fewer measurements for estimation.

In the event that device parameters are not provided along with measurement
records, an unsupervised approach comparing the input to a measurement record
reconstructed from parameter estimates can be used. This requires a map from
parameter estimates to measurement output, which is accomplished by adding a
layer combining an integrator of the physical model with a recurrent neural
network (RNN) to provide a model-free correction at every time step to account
for unanticipated effects, as long as they are consistent and small relative
to the dynamics driven by the parameters of interest. This employs the
capabilities of neural ordinary differential equations found in other works
[25, 26, 15] but enhanced to provide a completely model-free correction not
bound by assumptions of linearity or Markovianity, with a projection step to
return the corrected state to the manifold of physical density operators [27].
The output of the model is then the estimated solution to the unconditioned
master equation with learned corrections for effects beyond the master
equation, which can be compared to the measurement records to update model
parameters. While this case is not completely model-free, it also does not
make rigid assumptions about the dynamics and remains capable of accurate
parameter estimation in the presence of unanticipated effects that would
otherwise severely impact accuracy. The ML model is not burdened with learning
all of the physics; it just has to correct for small effects the model may
have missed, an approach known as discrepancy modeling [28, 29].

The fully unsupervised model can be viewed as a denoising autoencoder (DNA)
[30] where the physical parameters being estimated take on the role of the
latent space, with interpretability enforced by the presence of the physical
model in the decoder. The encoder is the map from measurement records to
parameters, and is the desired product of the training to be used for fast and
accurate parameter estimation.

While the unsupervised model clearly has applications for parameter estimation
in cases where nothing is known about the parameters of interest, the
supervised approach using just the encoder could still find applications for
systems where parameters are known at the time training data is generated, but
a prediction routine is still required in other circumstances, as would be the
case for detecting drift away from device calibration over time. We apply both
approaches to the specific problem of learning one- and two-body Hamiltonian
parameters in a two qubit system.

The manuscript is organized as follows. In Sec. II, we describe the physical
system of interest to us, which will be given by two qubits that are weakly
measured. In Sec. III, we describe our ML models and their supervised and
unsupervised training. In Sec. IV, we show results for the performance of the
models. In Sec. V, we summarize our results and discuss possible future
extensions of this work.

##  II Physical System

Our physical system of interest is that of two qubits with fixed position in a
microwave cavity as illustrated in Fig. 1. A single common mode of the cavity
is coupled to the computational degree of freedom of the qubits, and they are
coherently driven on resonance with a Rabi drive of frequency $\Omega$. A two-
qubit interaction term is present with magnitude $\epsilon$. A weak
measurement tone [31, 32] is applied to the cavity to probe the qubit state in
one of the $\\{X,Y,Z\\}$ directions for each qubit, with a measurement back-
action dephasing rate of $\kappa$. Upon adiabatic elimination of the cavity
mode [33], which we assume has dynamics evolving at a rate much faster than
time scales relevant to the qubits, we can describe the system with the
stochastic master equation (SME) [34, 35]

$\displaystyle d\rho$ $\displaystyle=-i[H,\rho]dt$
$\displaystyle+\sum_{i=1}^{2}\mathcal{D}[L_{i}](\rho)dt+\sum_{i=1}^{2}\sqrt{\frac{\eta}{2}}\mathcal{H}[L_{i}](\rho)dW_{t}^{(i)},$
(1a) $\displaystyle H$
$\displaystyle=\sum_{i=1}^{2}\frac{\Omega}{2}X_{i}+\epsilon Z_{1}Z_{2},$ (1b)

where $\\{X_{i},Y_{i},Z_{i}\\}$ is the set of Pauli operators for qubit
$i\in\left\\{1,2\right\\}$, $\mathcal{D}[L](\rho)=L\rho
L^{\dagger}-\frac{1}{2}\left\\{L^{\dagger}L\rho+\rho L^{\dagger}L\right\\}$ is
the Lindblad super-operator, $\mathcal{H}[L](\rho)=L\rho+\rho
L^{\dagger}-\rho\mbox{Tr}\left[\rho\left(L+L^{\dagger}\right)\right]$ is the
measurement super-operator, $\eta=0.1469$ is the efficiency of the
measurement, and $L_{i}=\sqrt{\kappa}C_{i},$ where
$C_{i}\in\\{X_{i},Y_{i},Z_{i}\\}$ is the weak measurement operator with
$\kappa=3.326$ radians/$\mu$s. Here we have suppressed the dependence on $t$
for $\rho$, and we have adopted the normalization convention $\hbar=1$. The
parameter values are selected to be comparable to those used in Ref. [19].

The stochastic differential equations for the measurement records are given by

$dr_{i}=\sqrt{\frac{\eta}{2}}\mbox{Tr}\left[\rho\left(L_{i}+L_{i}^{\dagger}\right)\right]dt+dW_{t}^{(i)}$
(2)

where $r_{i}$ is the weak measurement record for qubit
$i\in\left\\{1,2\right\\}$ and the independent Wiener increments
$dW_{t}^{(i)}$ are the same as those appearing in Eq. (1a). This is comparable
to the system studied in Ref. [19], except generalized to two qubits and with
the addition of the two-qubit interaction term in Eq. (1b). We also consider
cases where the weak measurement operator $L_{i}$ is different for each qubit.

The unconditioned master equation is given by averaging over all possible
trajectories of the Wiener processes and is given by:

$\displaystyle d\overline{\rho}$
$\displaystyle=-i[H,\overline{\rho}]dt+\sum_{i=1}^{2}\mathcal{D}[L_{i}](\overline{\rho})dt$
(3a) $\displaystyle d\overline{r_{i}}$
$\displaystyle=\sqrt{\frac{\eta}{2}}\mbox{Tr}\left[\overline{\rho}\left(L_{i}+L_{i}^{\dagger}\right)\right]dt$
(3b)

Figure 1: The physical system of two qubits in a microwave cavity subject to a
Rabi drive $\Omega$, single particle relaxation with rate $\gamma_{s}$, and
dephasing due to weak measurement back-action with rate $\kappa$

##  III Machine Learning Model

###  III.1 Preliminaries

We briefly review the kinds of architectures used in this work. For a more
detailed overview, we refer the reader to Appendix A.

Recurrent Neural Networks (RNNs) [36] are a class of artificial neural
networks designed to recognize patterns in sequences of data. Unlike
traditional neural networks, which assume all inputs (and outputs) are
independent of each other, RNNs are characterized by their ability to retain
information from previous inputs in the sequence through internal memory,
making them ideal for tasks involving sequential data.

Traditional RNNs struggle to learn and retain information over long sequences
due to problems like the contribution of information decaying over time [37].
This makes it hard for the network to maintain long-term dependencies. Long
Short-Term Memory (LSTM) networks [38] are a type of RNN architecture designed
to address the issue of long-term dependencies in sequence data. In this work,
we utilize LSTM-based models in order to model trajectories of weak
measurements.

We now describe a common design recipe used in tasks performing unsupervised
translation of data from one domain to another, often referred to,
collectively, as the family of encoder-decoder models. These are also called
sequence-to-sequence models when working with sequential data. Such models
produce an arbitrary length, context-dependent sequence, given an input
sequence. The length of the input sequence need not match the output sequence,
and individual elements of the input sequence need not have one-to-one
correspondence with individual elements of the output sequence.

Figure 2: A high-level description of the encoder-decoder design recipe.

The key idea is to use the encoder network to obtain a compact representation
of the input sequence (often called the context vector). This representation
is then passed to the decoder network to produce an output sequence by
augmenting it to the hidden state in the decoder. This is an extremely general
framework and potentially any neural network can be used in such an
arrangement (see Fig 2).

In our unsupervised setup, the encoder is an LSTM that returns the master
equation parameters as the context. The decoder is a numerical ODE integrator
implemented as an RNN.

###  III.2 Model Description

Our objective is to train a ML model to estimate the physical master equation
parameters in Eq. (1) by observing the weak measurement records described by
Eq. (2). Input data is provided for training of the ML model in the form of a
collection of _averaged_ weak measurement records, where the average is over
some number of records $d$ to reduce the effects of diffusive noise. The
output of the model is either (a) the physical master equation parameters of
interest in the supervised case where labeled training data is available, or
(b) the predicted measurement record satisfying the unconditioned master
equation (Eq. (3)) in the unsupervised case where true parameters are not
available. The former configuration can be interpreted as an encoder mapping
noisy trajectories to estimated parameter values, while the latter
configuration adds a decoder that takes the estimated parameters as input and
produces a predicted solution to the unconditioned master equation. The full
architecture is shown in Fig. 3 and can be viewed as a denoising autoencoder
taking noisy input trajectories and producing output with the noise removed.

The encoder begins with an average pooling operation in the time dimension to
further smooth the trend component of the voltages, followed by a long short-
term memory (LSTM) layer [39] to process the sequential input data and
consolidate information from multiple qubits into a single sequence. This is
followed by a feed-forward neural network of dense layers ending with the
parameter layer. The LSTM and feed-forward layers comprise the neural network
shown in Fig. 3(ii).

In the unsupervised case, interpretability of the latent variables as
parameters is enforced by using them directly in an enhanced numerical ODE
integrator in the decoder for the unconditioned master equation. The
integrator is implemented as a RNN with a custom cell that combines a single
step of Euler’s method with a correction produced by a standard LSTM cell.
This is shown in Fig. 4. The correction term is designed to compensate for
unanticipated dynamics in the physical model used by the numerical integrator.
The design leverages the concept of neural ODE [25] or ODE-RNN [26, 40], as it
uses a neural network to approximate corrections to the unconditioned master
equation as:

$d\overline{\rho}=f_{\text{model}}\left(t,\overline{\rho};\theta\right)dt+f_{\text{LSTM}}\left(t,\overline{\rho};\theta\right)dt,$
(4)

where $\theta$ is a set of ODE parameters, $f_{\text{model}}$ is the drift
term of the form of Eq. (3), and $f_{\text{LSTM}}$ is an LSTM cell and as such
is a non-linear, model free function of the time, state, and parameters.

The carry-state used in the standard LSTM architecture is passed between
evaluations of $f_{\text{LSTM}}$ at each time point, which means information
from the full history of states is potentially available, allowing for the
modeling of non-Markovian effects We ensure that the output operator is a
quantum state at each time step by restricting the free parameters in the
correction to have a Hermitian form to enforce Hermiticity, by implementing a
normalization step to ensure that the trace is preserved, and by performing an
orthogonal projection of the operator back onto the state space according to
the algorithm presented in Ref. [27]. In this respect it is analogous to the
approach used by the time-dependent variational principle (TDVP) algorithm
[41, 42] for modeling dynamics of matrix product states (MPS) [43, 44].

We emphasize that the enhanced numerical integrator is primarily relying on a
known physical model given by $f_{\text{model}}$ and using the ML components
in $f_{\text{LSTM}}$ only to compensate for discrepancies between this model
and the experimental system, which are assumed to be small relative to the
known dynamics. As such, it could be categorized along with other approaches
in ML for discrepancy modeling [28, 29] that operate under the principle that
small corrections are easier to learn than the full dynamics. A similar
approach is taken to modeling non-Markovian dynamics via RNN in Ref. [15],
though our approach differs in that it combines discrepancy modeling with a
completely general, non-linear correction to the state evolution based on a
neural ODE. In addition, as suggested in the outlook section of Ref. [15], our
investigation focuses on the situation where only a subset of the measurement
information necessary for full state reconstruction is available, so training
data volume requirements are substantially reduced.

It should be noted that while the integrator used in the decoder for this work
is an Euler integrator for the unconditioned master equation, this piece is
completely modular in design and can be replaced with any integrator for
simulating quantum dynamics. For example, a TDVP integrator using MPS could be
used for larger systems that are approximated well by a MPS ansatz, with the
trainable decoder parameters accounting for errors introduced by this
approach.

###  III.3 Model Interpretation

The output of the decoder is the estimated solution to the unconditioned
master equation (Eq. (3b)), and the loss is calculated as the mean squared
error (MSE) between this estimate and provided approximations of clean
measurements used as label values in the case where true parameter estimates
are not available. Thus, it is desirable that these labels contain as little
noise as possible, although estimation can still take place with noisy labels,
as will be seen in Sec. IV. Clean label values are realized by averaging over
many input groups sharing the same parameter values and using this same
average for all of the contributing input trajectory groups during training.
See Sec. III.4 for more details. In this way the model can be viewed as a
denoising autoencoder [30], with de-noised output provided for every noisy
input example, and with a latent space corresponding to the physical
parameters being estimated.

In our model, the parameters learned during training are the various weights
and biases of the neural networks in the encoder and the decoder. Notably, we
are not directly estimating master equation parameters during training, but
rather learning a function that maps noisy measurement records to physical
parameters. This approach differs from that of Ref. [19], where the physical
parameters are the result of an optimization and no such function is learned.
One advantage of our approach is that it enables extremely fast parameter
estimation once the model is trained, even when observing systems with
parameter sets not included in the training data. Furthermore, in the case
where the training data includes true physical parameter values, it allows us
to learn the map with the encoder only, bypassing the need for a physical
model entirely. In this case, scaling is limited by the amount and type of
training data required, rather than system size. The amount of training data
required will vary by application and requires further investigation. Finally,
learning a map for physical parameters as a function of measurement records
allows for the possibility of fewer measurements being required for prediction
versus training as the model learns to account for noise in the input data, as
illustrated by results in IV showing how estimation accuracy improves with
training set size even when the number of input records is held constant.

A second difference between our approach and that of Ref. [19] is that when
labeled training data is unavailable and a physical decoder is necessary, the
flexible correction scheme in the decoder allows the model to compensate for
non-Markovian or nonlinear dynamics since it is not bound by a Lindblad form,
while still being informed by the master equation for a first order model of
the dynamics.

Figure 3: Diagram of the full ML model. Voltage records are measured from the
cavity system producing noisy averaged voltage records (i), these are sent
through a neural network encoder (ii) producing system parameters as output
(iii). These parameters are used by the flex integrator decoder (iv) to
produce noise-free voltage estimates (v) Figure 4: Details of the decoder RNN
cell combining a standard ODE integrator with an LSTM update.

###  III.4 Training and Data

In order to generate the data for training, validation, and testing, we use an
Euler-Maruyama integrator to solve Eq. (2) for both qubits for each of $40$
values of $\epsilon$ evenly spaced on $[0,2)$ with fixed $\Omega=1.395$
radians/$\mu$s, and the same number of $\Omega$ values evenly spaced on
$[1,5)$ with fixed $\epsilon=1.0$ radians/$\mu$s for a total of $K=80$
($\Omega$, $\epsilon$) pairs. For each pair of values ($\Omega_{k}$,
$\epsilon_{k}$), $N=32{,}000$ measurement trajectories are simulated for
$T=4\mu s$. Half of the collection of pair of values ($\Omega_{k}$,
$\epsilon_{k}$), corresponding to every other parameter pair, is used for
training. The other half, corresponding to parameters midway between training
values, is split evenly to be used for validation and testing. This ensures
the training set contains a disjoint set of parameter pairs from the
validation/test set. The two endpoints of the validation/test set are excluded
from each end to ensure that the training data extends slightly beyond the
domain of validation/test values. Only weak measurement records are being
used, and no strong measurements are being performed.

During training, for each true parameter set
$\theta_{k}=(\Omega_{k},\epsilon_{k})$, trajectory groups of a preset size $d$
are randomly selected from the full training set and their averages are
provided as input values to the model, such that each mini-batch is comprised
of $M=N/d$ averaged trajectories $\\{\tilde{x}_{j,k}(t)\\}_{j=1}^{M}$ where
$\tilde{x}_{j,k}(t)=\frac{1}{d}\sum_{i\in I}r_{i,k}(t)$, $r_{i,k}(t)$ is
measurement record $i$ for parameter set $k$, and $I$ is a set of trajectory
indices of size $d$ randomly selected from $1,2,...,N$ without replacement
until all $N$ trajectories are used, which defines one epoch of training. In
this way, each mini-batch consists of a different set of noisy trajectories as
input to maximize the diversity of training examples. In the unsupervised
case, for all $j=1,...,M$, the value $x_{k}(t)$ to be used in the loss
function is the average over the full training set of trajectories associated
with the true parameter set $\theta_{k}$, such that
$x_{k}(t)=\frac{1}{N}\sum_{i=1}^{N}r_{i,k}(t)$.

To make the interpretation of our model as a denoising autoencoder more
concrete, we note that DNAs are characterized by a corruption process
$C(\tilde{y}|y)$ whereby noisy inputs $\tilde{y}$ are generated for each
uncorrupted training example $y$ [30]. In our case, the full trajectory means
$x_{k}(t)$ take on the role of $y$, approximating the solution to the
unconditioned SME, while the random selection of much smaller trajectory
groups is the process by which the corrupted data elements $\tilde{y}$ are
generated.

In the supervised case, the loss function is the MSE between predicted
parameter sets $\tilde{\theta}_{k}$ and the provided true values $\theta_{k}$

$\mathcal{L}=\frac{1}{pK}\sum_{k}\left|\left|\tilde{\theta}_{k}-\theta_{k}\right|\right|_{2}^{2},$
(5)

where $p$ is the dimension of $\theta_{k}$, the number of parameters being
estimated. In the unsupervised case, the loss function is the MSE between each
$x_{k}$ and the estimated clean measurement record

$\mathcal{L}=\frac{1}{MKN_{t}}\sum_{j,k,t}\left|x_{k}(t)-\mathcal{M}(\tilde{x}_{j,k})(t)\right|^{2}$
(6)

where $\mathcal{M}$ denotes the model, $t$ is the time index, and $N_{t}\equiv
T/\Delta t$ is the number of time points excluding the initial condition.

During evaluation, validation and test error is evaluated in a similar manner
to how training mini-batches are selected, with groups of size $d$ randomly
chosen from the full validation and test sets, followed by the calculation of
the MSE. This process is repeated 100 times and the average is taken as the
calculated MSE for each set.

Training is performed for multiple runs of 100 epochs, with a learning rate of
$3\times 10^{-3}$ and a decay rate of 0.99 per epoch, resetting the learning
rate after each run, continuing until the validation loss drops by less than
five percent run-to-run, and for the last 20 epochs of the final run. Once
converged, parameter MSE is computed for the validation and test sets. This
entire process is performed 100 times, each with a different random
initialization of the model parameters. The best model is considered to be the
one with the smallest validation loss after the final epoch, and this model is
used to evaluate the MSE for the test set. Hyperparameter tuning for the model
layer sizes was performed using a grid search for a physical system with
parameters distinct from those in Sec. IV.

##  IV Evaluation

To fully assess the performance of the model, we evaluate two distinct cases:
the encoder alone to evaluate the case of labeled training data, and the full
model for the case where labels are not available. We consider a range of
trajectory group sizes $d$ as well as noise-free data derived from an Euler
integration of the unconditioned master equation. The noise-free case is
denoted by $\infty$ as the number of trajectories in all tables.

To minimize the impact of the weak measurement back-action via the parameter
$\kappa$, in Secs. IV.1 and IV.2 the initial state is chosen to be spin-up in
the directions of measurement, and the directions of $X$ and $Y$ were used for
the first and second qubit, respectively. A different configuration is used in
Sec. IV.3, where more diverse measurements are needed to correct for
unanticipated single particle relaxation not explicitly present in the
decoder’s physical model.

###  IV.1 Supervised Learning

First, we study the impact of noise in the training set and how it helps or
hinders the model’s performance when predicting parameters from noisy data. We
do this by performing estimation with two models, one trained using noise-free
measurement values, and another set trained using eight groups of $d=4{,}000$
averaged measurement records. Both models are then evaluated on a test set of
noisy, averaged measurement records with the same group size $d$ as the noisy
training set. In both cases, the model with the best validation loss on its
respective validation set out of 100 randomly initialized models was used.

Figure 5 shows an example of the error of the estimated parameter pair
$\tilde{\theta}_{k}$ compared to truth for each of the 32 $k$ values in the
trimmed test set. We see that for many specific pairs and for the overall MSE,
the model trained on the noise-free measurements is producing less accurate
estimates than the model trained on noisy data. This is consistent in both the
$\Omega$ and $\epsilon$ errors. This suggests that the model learns to account
for noise, as expected, and training with noise on the level of measurements
used for prediction is beneficial.

Figure 5: Estimations of the parameter pair
$\theta_{k}=(\Omega_{k},\epsilon_{k})$ (top) and the squared error (bottom)
averaged over four noisy trajectories for a model trained on clean data and
noisy data, then evaluated on noisy data. Shuffle evaluated test set MSEs are
(8.67e-3, 5.83e-4) and (4.94e-3, 4.87e-4) respectively

Next we examine the impact of the training group size $d$. Table 1 lists test
set MSEs for a number of training sets where evaluation is taking place with
the same group size $d$ used for training. The best and median were taken from
the list of results sorted according to validation loss, as this is what would
be known in practice during training. Here we see a steady improvement as the
number of trajectories in a group increases, and noise decreases, saturating
in the case where noise-free data is used for training and evaluation. We note
that this test differs from the results in Fig. 5 where the model trained
without noise was evaluated on noisy data, as here we consider only the case
where training and test group sizes are equal.

$d$ | Best MSE | Median MSE | Mean MSE
---|---|---|---
2,000 | 9.87e-3, 8.76e-4 | 9.72e-3, 8.71e-4 | 9.75e-3, 9.79e-4
4,000 | 4.94e-3, 4.87e-4 | 5.23e-3, 5.44e-4 | 5.62e-3, 5.75e-4
8,000 | 4.17e-3, 2.74e-4 | 2.94e-3, 2.74e-4 | 3.38e-3, 4.13e-4
16,000 | 2.19e-3, 2.06e-4 | 1.40e-3, 1.92e-4 | 1.80e-3, 2.94e-4
$\infty$ | 1.01e-5, 2.04e-5 | 4.51e-6, 4.82e-5 | 1.32e-5, 5.68e-5
Table 1: Each cell shows the MSE pair for ($\Omega$, $\epsilon$) estimates on
supervised training set using the encoder only.

###  IV.2 Unsupervised Learning

We now consider the unsupervised case where the labels of the training data
are not known, but it is assumed that the physical model in the decoder is
correct, i.e. we are still not learning parameters for the drift function
correction $f_{\text{LSTM}}$, which we will consider in Sec. IV.3. We are
therefore using the MSE in Eq. (6) as the loss function.

First we examine how parameter estimatation accuracy for a fixed group size
$d=4{,}000$ varies with the total training set size $N$. Table 2 shows
$(\Omega,\epsilon)$ test set MSEs for the case where $\Omega$ is unknown but
fixed at a true value of 1.395 radians/$\mu$s and $\epsilon$ is allowed to
vary, for various training set sizes $N$. The same test set containing
$16{,}000$ trajectories is used for each row. From the table, we see that
accuracy improves significantly as the amount of training data increases, even
though the groups being presented to the model for parameter estimation remain
the same. This indicates that a greater diversity of noisy measurement records
when training results in models that can produce more accurate parameter
estimates when presented with the same number of measurements when performing
prediction. This motivates using $N=32{,}000$ going forward.

$N$ | Best MSE | Median MSE | Mean MSE
---|---|---|---
4,000 | 1.51e-3, 1.12e-3 | 6.58e-4, 9.86e-4 | 8.33e-4, 1.08e-3
8,000 | 3.27e-4, 7.71e-4 | 5.20e-4, 7.94e-4 | 5.38e-4, 8.81e-4
16,000 | 1.83e-4, 7.49e-4 | 3.20e-4, 8.41e-4 | 3.07e-4, 8.38e-4
32,000 | 8.14e-5, 7.16e-4 | 9.82e-5, 8.15e-4 | 1.45e-4, 8.12e-4
Table 2: MSE of $(\Omega,\epsilon)$ when using $d=4{,}000$ trajectories to
estimate parameters for varying training set sizes $N$ for an unsupervised
training set with fixed $\Omega$ and varying $\epsilon$. The values of
$\kappa$ and $\eta$ are assumed to be known exactly.

Next, we consider the impact of group size $d$ and measurement record time
spacing $\Delta t$ on accuracy. Table 3 again shows $(\Omega,\epsilon)$ test
set MSEs for the case where $\Omega$ is unknown but fixed at a true value of
1.395 radians/$\mu$s and $\epsilon$ is allowed to vary. It illustrates how
extremely low test set MSEs are achievable for $\Omega$ in this case, which is
expected given the high volume of training data available for a single value.
The accuracy of $\epsilon$ estimates depends on both the number of
trajectories used to create each input sequence, as well as the time spacing
at which measurements are recorded. To evaluate measurement records with
varying time spacing, trajectories simulated with $\Delta t=2^{-8}\mu$s are
sub-sampled to avoid introducing numerical integration error associated with
simulations with a larger time step.

Table 4 shows test set MSEs for the training set where both $\Omega$ and
$\epsilon$ are allowed to vary. This is a harder task as the model has less
training data for each unique $\Omega$ value, hence the loss in accuracy for
that parameter. The error in $\epsilon$ remains roughly the same or better as
in the fixed $\Omega$ case, with a best case root mean-square error of around
one percent of the median test value of the $\epsilon$ parameter. Here we see
a roughly linear trend in the MSE versus $d$ at smaller time steps, as
doubling the number of trajectories in a group roughly halves the MSE. This
trend breaks down, however, for the $\Delta t=2^{-4}\mu$s case, suggesting a
permissive time step threshold around $\Delta t=2^{-6}\mu$s at or below which
the expected trend in accuracy versus input data size is realized.

$d$ | $\Delta t=2^{-8}\,\mu s$ | $2^{-6}$ | $2^{-4}$
---|---|---|---
2,000 | 5.78e-5, 1.40e-3 | 1.20e-4, 1.40e-3 | 1.22e-4, 2.83e-3
4,000 | 8.14e-5, 7.16e-4 | 1.66e-4, 7.90e-4 | 1.19e-4, 2.09e-3
8,000 | 8.51e-5, 3.42e-4 | 7.16e-5, 3.63e-4 | 8.02e-5, 1.82e-3
16,000 | 8.31e-5, 1.57e-4 | 1.66e-4, 1.57e-4 | 1.19e-4, 1.89e-3
$\infty$ | 1.66e-6, 6.08e-6 | 7.21e-6, 1.04e-4 | 1.52e-5, 1.71e-3
Table 3: MSE of ($\Omega$, $\epsilon$) estimates on unsupervised training set with fixed $\Omega$ and varying $\epsilon$ and $N=32,000$. The values of $\kappa$ and $\eta$ are assumed to be known exactly. $d$ | $\Delta t=2^{-8}\,\mu s$ | $2^{-6}$ | $2^{-4}$
---|---|---|---
2,000 | 1.12e-2, 8.01e-4 | 1.10e-2, 8.42e-4 | 1.23e-2, 2.12e-3
4,000 | 5.87e-3, 4.56e-4 | 6.46e-3, 5.07e-4 | 6.39e-3, 1.68e-3
8,000 | 3.08e-3, 2.32e-4 | 2.98e-3, 2.87e-4 | 3.34e-3, 1.46e-3
16,000 | 1.68e-3, 1.19e-4 | 1.11e-3, 1.63e-4 | 1.58e-3, 1.39e-3
$\infty$ | 1.25e-5, 7.75e-6 | 1.76e-5, 7.07e-5 | 1.14e-4, 1.23e-3
Table 4: MSE of ($\Omega$, $\epsilon$) estimates on unsupervised training set
with varying $\Omega$ and $\epsilon$ and $N=32,000$. The values of $\kappa$
and $\eta$ are assumed to be known exactly.

###  IV.3 Model Correction

In this section we demonstrate the ability of the decoder to correct for
dynamics not explicitly considered in the physical model (Eq. (1a)). This is
done by enabling training for the decoder LSTM parameters. For the datasets,
we simulate 30,000 trajectories with fixed $\Omega=1.395$ radians/$\mu$s but
with varying $\epsilon$, which is the parameter to be estimated, while adding
a dissipative term for each qubit corresponding to the Lindblad operator
$\mathcal{D}[\sqrt{\gamma_{s}}\sigma_{i}^{-}](\rho)$ in Eq. (1a), where
$\gamma_{s}=0.1$ and $\sigma_{i}^{-}=\frac{1}{2}(X_{i}-iY_{i})$ is the single
particle relaxation operator mapping the excited state to the ground state.

In this case, we take a more diverse set of measurements, simulating 10,000
trajectories measuring in each of the $X$, $Y$, and $Z$ directions for both
qubits, and relaxing $\kappa$ to one-fourth the value used in the last section
to reduce measurement back-action. Spin-up in the $Z$ direction is the initial
state for each qubit. This more cautious approach to measurement is warranted
if completely unknown effects are expected to be present. More information
about the type of phenomenon, but not necessarily the magnitude, could allow
for a more targeted measurement scheme, but here we keep it general. The
trajectory group size used for each input was $d=5{,}000$, and the best of 20
randomly initialized models was selected for the results in this section.

The results of the parameter estimation with and without the correction are
shown in Fig. 6 and Table 5. Here we see that single particle relaxation has
introduced a significant bias to the estimated $\epsilon$ parameters when
unaccounted for in the model, but the decoder LSTM has successfully corrected
for the effect, returning the MSE to a value much closer to where it would
have been had the physical model explicitly accounted for it.

Figure 6: Estimate of $\epsilon$ and $P(Y_{i}=0)$ when $\epsilon=1.7$ and measuring $X$ with spin-up in $Z$ as the initial state for both qubits for training data simulated with single particle relaxation rate $\gamma_{s}=0.1$ but a model that does not explicitly account for $\gamma_{s}$. Results are shown for the physical model alone and with a correction learned by the decoder LSTM $d$ | $\gamma_{s}=0.1$ | $\gamma_{s}=0.0$ |  $\gamma_{s}=0.0$ \+ correction
---|---|---|---
5,000 | 2.84e-3 | 0.142 | 3.82e-3
$\infty$ | 2.12e-6 | 0.145 | 2.71e-4
Table 5: MSE of $\epsilon$ estimates using different $\gamma_{s}$ values in
the model. The final column shows results when learning is enabled for the
free parameters in the decoder to account for the unanticipated term
$\gamma_{s}$. Measurement time spacing is $\Delta t=2^{-8}\mu s$

##  V Conclusions

We have proposed a machine learning model based on a denoising autoencoder
capable of direct estimation of physical parameters in a system modeled by a
stochastic master equation from weak measurement records. The model is capable
of learning in a supervised and an unsupervised context, and it can accurately
predict parameters for systems not seen in the training data. While leveraging
the use of a master equation integrator to enable unsupervised learning, the
autoencoder is robust to unanticipated dynamics not included in its physical
model. We have demonstrated this in the case of unanticipated Lindblad
dissipative terms. Previous work has established that LSTM models are capable
of learning dynamics beyond unknown Lindbladian dissipation [21, 15], and it
is an interesting subject for future work to examine the robustness of our
model to the presence of unanticipated non-Markovian and nonlinear dynamics in
the training data.

Another potential subject for future investigation is the ability of the model
to estimate parameters for much larger systems in the supervised context,
where it may not necessarily be subject to the exponential scaling of the
Hilbert space with system size that often limits other approaches to parameter
estimation.

###### Acknowledgements.

This material is based upon work supported by the National Science Foundation
the Quantum Leap Big Idea under Grant No. OMA-1936388. AKR and CM acknowledge
support from NSF QLCI OMA-2016244. This article has been co-authored by an
employee of National Technology & Engineering Solutions of Sandia, LLC under
Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The
employee owns all right, title and interest in and to the article and is
solely responsible for its contents. The United States Government retains and
the publisher, by accepting the article for publication, acknowledges that the
United States Government retains a non-exclusive, paid-up, irrevocable, world-
wide license to publish or reproduce the published form of this article or
allow others to do so, for United States Government purposes. The DOE will
provide public access to these results of federally sponsored research in
accordance with the DOE Public Access Plan
https://www.energy.gov/downloads/doe-public-access-plan.

##  Appendix A Recurrent Neural Networks

In this section, we describe RNNs in additional detail. Our presentation is
mostly standard and closely follows the treatment of [45].

RNNs are a class of models which employ a cycle within their network structure
such that output of a unit depends in some manner on its own output at some
previous step, so they have an underlying recursive structure [36]. This
structure allows RNNs to model sequences of data where information at the
current index depends on information processed in the past. While the term
RNNs refers to a general class of models which share this recursive pattern,
in the ML literature, it is used to refer to a specific kind architecture with
a simple cyclic pattern (sometimes also called Elman Networks) [46].

Figure 7: An illustration describing the basic architecture of a RNN. The
model on the left can be thought of as ‘unrolling’ itself in time.

The model operates as follows (see Fig. 7). Given an input vector $x_{t}$,
where $t$ is an index over the sequence to be modelled, we first multiply
$x_{t}$ with a weight matrix $W$ and pass the result through a non-linearity
$g$ to produce values in a hidden layer $h_{t}$. However, unlike feed-forward
networks we supplement this value with values from previous hidden layers
$h_{t-1}$, multiplied by a weight matrix $U$ as well. This can be summarized
as:

$h_{t}=g(Uh_{t-1}+Wx_{t}).$ (7)

Additionally, each unit of the RNN can produce an output $y_{t}$, which can be
parameterized as follows:

$y_{t}=f(Vh_{t}),$ (8)

where $f$ is usually taken to be a softmax function and $V$ is a weight
matrix. This network is then trained through an analog of backpropagation
which unrolls the sequence along its index referred to as Backpropagation
through Time (BPTT) [47, 48, 49].

In practice, it has been observed that such networks can be hard to train for
longer sequence lengths. This is mainly due to the “bottleneck” caused by the
previous hidden state vector $h_{t-1}$ having to summarize _all_ of the
relevant information at the current index inside a fixed size vector. Another
difficulty encountered is due to repeated multiplications of small values in
BPTT causing poor gradient flow to steps further back in time (referred to as
the ‘vanishing gradient problem’ [37].

LSTMs are a commonly used extension of RNNs to tackle the shortcomings pointed
out above [38]. The central idea is to gain more precise control over
information being allowed to flow through the network. In addition to the
hidden state, an additional so-called ‘memory cell’ (also called the cell
state) is computed, which, like the hidden state, allows the flow of
information along the length of the sequence. The memory cell is designed such
that information can be added to it or subtracted from it. We refer readers to
[30] for details of how this is done.

##  Appendix B Baseline: Random Forests

To justify our choice of employing LSTMs as the backbone of our proposed
architecture, we provide a comparison using Random Forests [50] for the
encoder-only experiments described in Section IV.1 with. We choose to compare
against Random Forests because they have been shown to have robust off-the-
shelf performance across a wide variety of tasks [51]. In this section we give
a brief overview of Random Forests and a comparison with results in Table 1,
which justify the use of Neural Network based architectures in our subsequent
experiments.

###  B.1 Model Description

Here we provide a brief overview of the Random Forest model, closely following
the treatment in Ref. [52]. Decision trees are simple tree-based models that
make decisions by recursively partitioning the input feature space into
regions and then returning the label or output associated with that region.
While decision trees are often robust to noisy features and exhibit a high
degree of interpretability, they can grow to large depths, thereby overfitting
the training set [52]. This effect is often addressed by learning a Random
Forest, i.e., a collection (ensemble) of decision trees, each of which is
constructed independently by training on a random subset of the data, in
addition to using some random choices in the algorithm. This introduces
diversity among the trees, preventing the model from becoming too dependent on
any single feature. For regression tasks (i.e. the output is continuous-
valued), as in this work, the final prediction is typically the average of the
predictions made by each tree.

The combination of random sampling and majority voting (or averaging) is a
technique known as Bagging (or Bootstrap Aggregating)[53]. Random Forests use
bagging to create an ensemble that is more robust (reducing uncertainty on the
output) and less prone to overfitting compared to a single decision tree.
Apart from their versatility and robustness, Random Forests provide insights
into feature importance, helping to identify the most influential variables in
the model.

###  B.2 Training Setup

Our aim is to compare the prediction accuracy for the parameter $\epsilon$ for
a fixed $\Omega$ in the case of supervised learning between the Random Forest
model and the model described in the main text. To do so, we use the same
procedure as described in Sec. IV.1 of the main text. The training set
contains averaged trajectories sampled randomly in groups of varying size $d$,
and the trajectories in the dataset correspond to $40$ values of $\epsilon$
evenly spaced on $[0,2)$ for a fixed $\Omega=1.395$ radians/$\mu$s. The
crucial hyperparameters for our method are the number of trees in the ensemble
(set to 100) and the maximum depth of each tree (set to 25). These were found
via 5-fold cross validation over the dataset.

A clear point of difference compared to the LSTM-based model is the form of
the input given to the Random Forest model. Unlike the LSTM model, input to
the Random Forest model _must_ be 1-dimensional. We use the concept of
‘dimensionality’ to refer to the number of weak measurement values at each
timestep. Our input is 2-dimensional since we have a weak measurement value
for each qubit. This dichotomy requires us to ‘flatten’ the input values into
a single dimension. This can be done in two ways: either by keeping
measurement values corresponding to a particular time step close together
(column-major order) or by concatenating the values of the second qubit after
_all_ the values of the first qubit (row-major order). We find that the choice
of ordering makes little difference in our results as we demonstrate next.

###  B.3 Results

We present the impact of the choice of input order in our Random Forest model
in Table 6. We find the difference in MSE values to be minimal.

$d$ | Best MSE | Median MSE | Mean MSE
---|---|---|---
2,000 | 2.35e-3, 2.63e-3 | 4.13e-3, 4.08e-3 | 4.41e-3, 4.52e-3
4,000 | 1.23e-3, 1.31e-3 | 3.44e-3, 3.52e-3 | 3.69e-3, 3.70e-3
8,000 | 7.49e-4, 7.34e-4 | 1.61e-3, 1.68e-3 | 1.86e-3, 1.86e-3
16,000 | 3.36e-4, 3.26e-4 | 1.47e-3, 1.50e-3 | 1.65e-3, 1.71e-3
$\infty$ | 2.34e-4, 3.48e-5 | 1.54e-3, 1.37e-3 | 1.75e-3, 1.79e-3
Table 6: Each cell shows MSEs for Random Forest models with row-major input
flattening and column-major input flattening respectively.

We now compare the results of the Random Forests model in column-major order
to the LSTM model in Table 7. We find about an order of magnitude difference
in the MSEs achieved by the Random Forest model compared to the LSTM model.
The trend of improved accuracy with increasing sizes of $d$, also observed in
the LSTM results, continues to hold in this case. To get a better sense of the
uncertainty, we also present the mean MSEs and their one standard deviation
error bars for both our models in Fig. 8. The fact that our choice of
‘flattening’ does not affect overall performance of the Random Forest model
suggests that the difference in performance between the Random Forest model
and the LSTM model could be attributed to the LSTM-based model taking the
sequential nature of the data into account versus the Random Forest model
treating each trajectory as just a high-dimensional vector. We leave a
thorough investigation of these issue to future work.

Figure 8: Mean MSEs and their 1 standard deviation error bars for the Random Forests model with different input orderings and the LSTM model. $d$ | Best MSE | Median MSE | Mean MSE
---|---|---|---
2,000 | 2.63e-3, 1.47e-3 | 4.08e-3, 1.51e-3 | 4.52e-3, 1.63e-3
4,000 | 1.31e-3, 7.37e-4 | 3.52e-3, 7.16e-4 | 3.70e-3, 8.56e-4
8,000 | 7.34e-4, 4.39e-4 | 1.68e-3, 6.00e-4 | 1.86e-3, 4.67e-4
16,000 | 3.26e-4, 1.83e-4 | 1.50e-3, 2.13e-4 | 1.71e-3, 2.68e-4
$\infty$ | 3.48e-5, 8.11e-6 | 1.37e-3, 5.59e-5 | 1.79e-3, 7.41e-5
Table 7: Each cell shows the MSE pair for $\epsilon$ estimates on supervised
training set using the encoder only. The first value is the MSE of the Random
Forest model using column-major ordering, and second value corresponds to the
MSE of the LSTM-based model.

We also provide a comparison of the number of parameters used in both our
models in Table 8. Since the number of parameters in a Random Forest model
depends on the depth of the trees, which in turn depend on the complexity of
the dataset, we report the maximum parameters across all models obtained after
training.

Method | Number of parameters
---|---
LSTM | 1,619,316
Random Forest | 129,380
Table 8: Number of parameters used in the two models studied.

## References

  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Sarovar  _et al._ [2020] M. Sarovar, T. Proctor, K. Rudinger, K. Young, E. Nielsen, and R. Blume-Kohout, Detecting crosstalk errors in quantum information processors, Quantum 4, 321 (2020).
  * Tripathi  _et al._ [2022] V. Tripathi, H. Chen, M. Khezri, K.-W. Yip, E. Levenson-Falk, and D. A. Lidar, Suppression of crosstalk in superconducting qubits using dynamical decoupling, Phys. Rev. Appl. 18, 024068 (2022).
  * Kitaev [1997] A. Y. Kitaev, Quantum computations: algorithms and error correction, Russian Mathematical Surveys 52, 1191 (1997).
  * Aharonov and Ben-Or [1997] D. Aharonov and M. Ben-Or, Fault-tolerant quantum computation with constant error, in  _Proceedings of the Twenty-Ninth Annual ACM Symposium on Theory of Computing_, STOC ’97 (Association for Computing Machinery, New York, NY, USA, 1997) p. 176–188.
  * Preskill [1998] J. Preskill, Reliable quantum computers, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454, 385 (1998).
  * Knill  _et al._ [1998] E. Knill, R. Laflamme, and W. H. Zurek, Resilient quantum computation, Science 279, 342 (1998).
  * Gottesman [1998] D. Gottesman, Theory of fault-tolerant quantum computation, Phys. Rev. A 57, 127 (1998).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Carleo  _et al._ [2018] G. Carleo, Y. Nomura, and M. Imada, Constructing exact representations of quantum many-body systems with deep neural networks, Nature Communications 9, 5322 (2018).
  * Carrasquilla  _et al._ [2019] J. Carrasquilla, G. Torlai, R. G. Melko, and L. Aolita, Reconstructing quantum states with generative models, Nature Machine Intelligence 1, 155 (2019).
  * Rocchetto  _et al._ [2018] A. Rocchetto, E. Grant, S. Strelchuk, G. Carleo, and S. Severini, Learning hard quantum distributions with variational autoencoders, npj Quantum Information 4, 28 (2018).
  * Carrasquilla and Torlai [2021] J. Carrasquilla and G. Torlai, How to use neural networks to investigate quantum many-body physics, PRX Quantum 2, 040201 (2021).
  * Torlai  _et al._ [2023] G. Torlai, C. J. Wood, A. Acharya, G. Carleo, J. Carrasquilla, and L. Aolita, Quantum process tomography with unsupervised learning and tensor networks, Nature Communications 14, 2858 (2023).
  * Banchi  _et al._ [2018] L. Banchi, E. Grant, A. Rocchetto, and S. Severini, Modelling non-markovian quantum processes with recurrent neural networks, New Journal of Physics 20, 123030 (2018).
  * Brunton  _et al._ [2016] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113, 3932 (2016).
  * Champion  _et al._ [2019] K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton, Data-driven discovery of coordinates and governing equations, Proceedings of the National Academy of Sciences 116, 22445 (2019).
  * Lusch  _et al._ [2018] B. Lusch, J. N. Kutz, and S. L. Brunton, Deep learning for universal linear embeddings of nonlinear dynamics, Nature Communications 9, 4950 (2018).
  * Genois  _et al._ [2021] E. Genois, J. A. Gross, A. Di Paolo, N. J. Stevenson, G. Koolstra, A. Hashim, I. Siddiqi, and A. Blais, Quantum-tailored machine-learning characterization of a superconducting qubit, PRX Quantum 2, 040355 (2021).
  * Flurin  _et al._ [2020] E. Flurin, L. S. Martin, S. Hacohen-Gourgy, and I. Siddiqi, Using a recurrent neural network to reconstruct quantum dynamics of a superconducting qubit from physical observations, Phys. Rev. X 10, 011006 (2020).
  * Koolstra  _et al._ [2022] G. Koolstra, N. Stevenson, S. Barzili, L. Burns, K. Siva, S. Greenfield, W. Livingston, A. Hashim, R. K. Naik, J. M. Kreikebaum, K. P. O’Brien, D. I. Santiago, J. Dressel, and I. Siddiqi, Monitoring fast superconducting qubit dynamics using a neural network, Phys. Rev. X 12, 031017 (2022).
  * An  _et al._ [2023] Z. An, J. Wu, M. Yang, D. L. Zhou, and B. Zeng, Unified quantum state tomography and hamiltonian learning using transformer models: A language-translation-like approach for quantum systems (2023), arXiv:2304.12010 [quant-ph] .
  * Castelano  _et al._ [2023] L. K. Castelano, I. Cunha, F. S. Luiz, M. V. de Souza Prado, and F. F. Fanchini, Physics informed neural networks learning a two-qubit hamiltonian (2023), arXiv:2310.15148 [quant-ph] .
  * Raissi  _et al._ [2017] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics informed deep learning (part ii): Data-driven discovery of nonlinear partial differential equations (2017), arXiv:1711.10566 [cs.AI] .
  * Chen  _et al._ [2018] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, Neural ordinary differential equations, in  _Advances in Neural Information Processing Systems_, Vol. 31, edited by S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (Curran Associates, Inc., 2018).
  * Habiba and Pearlmutter [2020] M. Habiba and B. A. Pearlmutter, Neural ordinary differential equation based recurrent neural network model, in  _2020 31st Irish Signals and Systems Conference (ISSC)_ (2020) pp. 1–6.
  * Smolin  _et al._ [2012] J. A. Smolin, J. M. Gambetta, and G. Smith, Efficient method for computing the maximum-likelihood quantum state from measurements with additive gaussian noise, Phys. Rev. Lett. 108, 070502 (2012).
  * Kaheman  _et al._ [2019] K. Kaheman, E. Kaiser, B. Strom, J. N. Kutz, and S. L. Brunton, Learning discrepancy models from experimental data (2019), arXiv:1909.08574 [cs.LG] .
  * de Silva  _et al._ [2020] B. M. de Silva, D. M. Higdon, S. L. Brunton, and J. N. Kutz, Discovery of physics from data: Universal laws and discrepancies, Frontiers in Artificial Intelligence 3, 10.3389/frai.2020.00025 (2020).
  * Goodfellow  _et al._ [2016] I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (MIT Press, 2016) http://www.deeplearningbook.org.
  * Jacobs and Steck [2006] K. Jacobs and D. A. Steck, A straightforward introduction to continuous quantum measurement, Contemporary Physics 47, 279 (2006).
  * Brun [2002] T. A. Brun, A simple model of quantum trajectories, American Journal of Physics 70, 719 (2002).
  * Bonifacio  _et al._ [1971] R. Bonifacio, P. Schwendimann, and F. Haake, Quantum statistical theory of superradiance. i, Phys. Rev. A 4, 302 (1971).
  * Blais  _et al._ [2021] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit quantum electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).
  * Gambetta  _et al._ [2008] J. Gambetta, A. Blais, M. Boissonneault, A. A. Houck, D. I. Schuster, and S. M. Girvin, Quantum trajectory approach to circuit qed: Quantum jumps and the zeno effect, Phys. Rev. A 77, 012112 (2008).
  * Rumelhart  _et al._ [1986] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, Learning internal representations by error propagation, in  _Parallel Distributed Processing: Explorations in the Microstructure of Cognition, Volume 1: Foundations_ , edited by D. E. Rumelhart and J. L. Mcclelland (MIT Press, Cambridge, MA, 1986) pp. 318–362.
  * Hochreiter  _et al._ [2001] S. Hochreiter, Y. Bengio, P. Frasconi, and J. Schmidhuber, Gradient flow in recurrent nets: the difficulty of learning long-term dependencies, in  _A Field Guide to Dynamical Recurrent Neural Networks_ , edited by S. C. Kremer and J. F. Kolen (IEEE Press, 2001).
  * Hochreiter and Schmidhuber [1997a] S. Hochreiter and J. Schmidhuber, Long short-term memory, Neural Computation 9, 1735 (1997a).
  * Hochreiter and Schmidhuber [1997b] S. Hochreiter and J. Schmidhuber, Long short-term memory, Neural Comput. 9, 1735–1780 (1997b).
  * Rubanova  _et al._ [2019] Y. Rubanova, R. T. Q. Chen, and D. K. Duvenaud, Latent ordinary differential equations for irregularly-sampled time series, in  _Advances in Neural Information Processing Systems_, Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).
  * Haegeman  _et al._ [2011] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Time-dependent variational principle for quantum lattices, Phys. Rev. Lett. 107, 070601 (2011).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94, 165116 (2016).
  * Perez-Garcia  _et al._ [2007] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Matrix product state representations, Quantum Info. Comput. 7, 401 (2007).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011).
  * Jurafsky and Martin [2024] D. Jurafsky and J. H. Martin,  _Speech and Language Processing An Introduction to Natural Language Processing, Computational Linguistics, and Speech Recognition_ , 3rd ed. (2024).
  * Elman [1990] J. L. Elman, Finding structure in time, Cognitive Science 14, 179 (1990).
  * Robinson  _et al._ [1987] A. Robinson, F. Fallside, and U. of Cambridge. Engineering Department,  _The Utility Driven Dynamic Error Propagation Network_ (University of Cambridge Department of Engineering, 1987).
  * Werbos [1988] P. J. Werbos, Generalization of backpropagation with application to a recurrent gas market model, Neural Networks 1, 339 (1988).
  * Mozer [1995] M. C. Mozer, A focused backpropagation algorithm for temporal pattern recognition, in  _Backpropagation: Theory, Architectures, and Applications_ (L. Erlbaum Associates Inc., USA, 1995) p. 137–169.
  * Breiman [2001] L. Breiman, Random forests, Machine Learning 45, 5 (2001).
  * Fernández-Delgado  _et al._ [2014] M. Fernández-Delgado, E. Cernadas, S. Barro, and D. Amorim, Do we need hundreds of classifiers to solve real world classification problems?, Journal of Machine Learning Research 15, 3133 (2014).
  * Hastie  _et al._ [2001] T. Hastie, R. Tibshirani, and J. Friedman,  _The Elements of Statistical Learning_ , Springer Series in Statistics (Springer New York Inc., New York, NY, USA, 2001).
  * Breiman [1996] L. Breiman, Bagging predictors, Machine Learning 24, 123 (1996).