Bilinear dynamic mode decomposition for quantum control

Data-driven methods for establishing quantum optimal control (QOC) using time-dependent control pulses tailored to specific quantum dynamical systems and desired control objectives are critical for many emerging quantum technologies. We develop a data-driven regression procedure, bilinear dynamic mode decomposition (biDMD), that leverages time-series measurements to establish quantum system identification for QOC. The biDMD optimization framework is a physics-informed regression that makes use of the known underlying Hamiltonian structure. Further, the biDMD can be modified to model both fast and slow sampling of control signals, the latter by way of stroboscopic sampling strategies. The biDMD method provides a flexible, interpretable, and adaptive regression framework for real-time, online implementation in quantum systems. Further, the method has strong theoretical connections to Koopman theory, which approximates non-linear dynamics with linear operators. In comparison with many machine learning paradigms, it requires minimal data and the biDMD model is easily updated as new data is collected. We demonstrate the efficacy and performance of the approach on a number of representative quantum systems, showing that it also matches experimental results.


Introduction
Quantum optimal control (QOC) is a comprehensive mathematical framework for quantum control in which time-dependent control pulses are tailored to specific quantum dynamical systems and desired control objectives [1]. QOC algorithms are critical for emerging new quantum technologies in scientific and engineering disciplines, including computing, communications, simulation and sensing [2]. Despite the unique challenges of quantum theory [3,4], standard control optimization procedures have been developed and refined [1,5,6]. Many of the seminal innovations in QOC were developed for applications in Nuclear Magnetic Resonance (NMR) [7]; today there are many numerical methods and tools useful in a wide range of QOC application areas [8][9][10]. In QOC, practical control design ultimately depends on an experimentally-accurate model of the governing quantum dynamics and the action of controls. In this manuscript, we introduce the data-driven regression framework known as dynamic mode decomposition (DMD) to construct control models directly from time-series measurement data. Our DMD method is tailored to the bilinear structure of quantum control dynamics. Moreover, it is a completely data-driven approach that can accommodate both fast and slow (stroboscopic) sampling of control signals.
Many data-driven methods for characterizing quantum devices for purposes of quantum control have been developed under a variety of modelling assumptions [11]. For example, the

Quantum dynamical systems
The first complete formulation of quantum dynamical systems occurred in the 1930s with the operator-theoretic perspective of the Dirac-von Neumann axioms [40]. This description replaced the notion of a possibly non-linear equation of motion defined over a classical phase space with an infinite dimensional linear operator algebra acting on a Hilbert space. This formulation was not restricted to quantum mechanics. Indeed, an operator-theoretic linearization of classical dynamical systems was contemporaneously made using the Koopman-von Neumann theory [29,30]. We review the dynamics of quantum states here before transitioning to abstract states of measurement data in Section 3.

Autonomous systems
In quantum mechanics, the state of an N -dimensional physical system is represented by a unit vector, or ket, |ψ in a complex vector space C N . The dynamics of a quantum state are given by the Schrödinger equation, where the Hamiltonian operator, H 0 , has been assumed without loss of generality to be traceless , and Plank's constant has been set to unity. More generally, an ensemble of pure quantum states can be completely characterized, in the sense of its measurement statistics, by a density matrix ρ(t); that is, a non-negative self-adjoint operator in C N ×N with trace one. The Liouville-von Neumann, or quantum Liouville equation, describes the time evolution of a density matrix [41]: The quantum Liouville operator L 0 is sometimes known as a super-operator because of its linear action on the space of operators. There are a number of ways to vectorize the quantum Liouville equation. We will use the language of the familiar Bloch vector (also known as the vector of coherence) to describe passing to a vector of differential equations. Writing ρ(t) = 1/N + N 2 −1 j=1 Tr(ρσ j )σ j with {σ j } N 2 −1 j=1 as a complete and orthonormal basis for traceless Hermitian operators, we can then define x j := Tr(ρσ j ) so that (2) becomes [4,42] ∂ ∂t For a more complete description of this process, see Appendix A. Obtaining an expectation value x j (t) := Tr(ρ(t)σ j ) requires a collection of identically-prepared quantum states. In this paper, the "measurement" of a quantum state as represented by x(t) is understood to mean the measurement of such an ensemble.
In terms of the eigenvectors v j and eigenvalues λ j of L 0 , the solution to (3) is where the coefficients c j are the coordinates of x 0 in the eigenvector basis. For continuous dynamics as in (4), there exists an equivalent discrete time description of the system x(t) sampled at intervals of ∆t.

Non-autonomous systems
Many important classical non-autonomous control systems can be formulated as control-affine dynamical systems. Transforming to the operator theoretic perspective of the Koopman-von Neumann theory, these are bilinear dynamical systems [43]. In this section, we review nonautonomous dynamical systems and establish the bilinear formulation of the quantum control problem.

Direct actuation
Direct actuation is a common situation in classical control theory in which an autonomous dynamical system like (3) undergoes a linear response to control inputs u(t) ∈ R Nc [37]: Recall that if the control is input under a zero-order hold across ∆t, the discretization x n = x(t 0 + (n − 1)∆t) and u n = u(t 0 + (n − 1)∆t) transforms (5) into the discrete-time dynamical system Significantly, the dynamics remain control-affine for any span ∆t.

Bilinear quantum dynamics
The control of a quantum system can be modelled using N c real-valued control functions, u j (t), coupled to corresponding time-independent interaction Hamiltonians, H j , such that the dynamics are described by the bilinear Schrödinger equation, Following Appendix A as in (3), the bilinear Schrödinger equation induces the vectorized bilinear quantum Liouville equation [4,42], by identifying H n → L n for n = 0, 1, . . . , N c and x ∈ R N 2 −1 . This linear differential equation can be integrated in the usual way to obtain If a discretization is taken such that x n = x(t 0 + (n − 1)∆t) and u j;n = u j (t 0 + (n − 1)∆t), then to first order such that the discrete time dynamical system remains control-affine [39]. Note that in this first-order approximation the computation of derivatives via finite differences and the computation of the discrete time dynamical system are equivalent.

Dynamic Mode Decomposition
In recent years, a variety of practical computational tools have encouraged more widespread adoption of the Koopman-von Neumann perspective in the dynamical systems community [31,44]. The Dynamic Mode Decomposition (DMD) broadly refers to a suite of numerical methods originating in the fluid dynamics community for the purpose of studying coherent spatiotemporal structures in complex fluid flows [32]. DMD combines standard dimensionality reduction via the singular value decomposition with Fourier transforms in time. In application, all variations of the DMD algorithm involve collecting a time series of experimental measurements in order to compute a reduced-order dynamical model based on DMD modes and eigenvalues. The DMD modes identify coherent structures in the measured space and the DMD eigenvalues define the growth, decay, and oscillation frequency of each mode. The following sections will review the DMD framework.

DMD
The standard DMD is a data-driven algorithm defined via the observed trajectories from either experimental data or a numerical simulation of a dynamical system ∂ ∂t where x(t) ∈ R N 2 −1 to remain consistent with Section 2. (DMD ultimately deals with data and is agnostic regarding the underling system that produced it.) The first step of the algorithm is to assemble a sequential measurement record {x(t 1 ) = x 1 , x 2 , . . . , x M }, t 1 < t 2 < · · · < t m , into snapshot matrices: so X, X ∈ R N 2 −1×M −1 . In this paper we assume uniform sampling with t m := t 0 + (m − 1)∆t, but this assumption can be relaxed [33]. In its simplest form, DMD is a regression algorithm that estimates the propagator matrix X ≈ AX (13) by solving the optimization problem where M F = J j=1 K k=1 m jk is the Frobenius matrix norm defined for any given matrix M. The least-squares solution to (14) is A = X X + where + denotes the Moore-Penrose pseudoinverse of a matrix. From the eigendecomposition A = WΩW −1 , define the DMD modes w j as the columns of W. Future states can then be predicted as in a discrete analogue to (4) using where b are the coefficients of the initial condition in the eigenvector basis, x 0 = Wb.
The eigendecomposition of A ∈ R N 2 −1×N 2 −1 can easily become numerically intractable for large N . However, the big data limit is nothing unusual for the problems appearing in data science and fluid dynamics that DMD has already addressed with success. The DMD algorithm circumvents this problem by projecting the high-dimensional snapshot matrix X onto a low-rank subspace defined by R principle orthogonal decomposition (POD) modes computed from the singular value decomposition of X (note that because it is computed from the snapshot matrices, the rank of A is at most M − 1 when N 2 > M ). In this way, DMD determines a low-rank approximatioñ A ∈ R R×R to describe the dynamics of the measured trajectories. This paper introduces DMD for quantum control by way of simple low-rank examples. However, the POD-based DMD is also amenable to the study of many-body quantum states [45]. DMD accommodates tools like tensor networks [46,47] and compressed sensing [35,48].

DMD with Control (DMDc)
Dynamic Mode Decomposition with control (DMDc) [37] is an algorithm developed for modelling the special case of direct actuation: with u(t) ∈ R Nc . Like DMD in (12), the size-M measurement record associated with the trajectory of x is assembled into snapshot matrices X and X . This time, the control inputs, u n := u(t n ), are also collected: DMDc is a regression-based approach to system identification that disambiguates the intrinsic dynamics, A, and the effects of control, B: The DMDc algorithm achieves this disambiguation by interpreting the best-fit solution according to DMDc brings all the benefits of the DMD algorithm to a model-free framework for experimental control optimization.

Bilinear DMD (biDMD)
Bilinear Dynamic Mode Decomposition (biDMD) is a data-driven system identification framework for bilinear non-autonomous control systems: Construct the snapshot matrices X, X , and U as in (12) and (17) using the size-M measurement record of the trajectory x(t) and the corresponding input record u(t). Additionally, construct the Figure 1: The trajectory of a qubit driven by a linearly-polarized semi-classical drive u(t) (Hamiltonian: is shown on the Bloch sphere in Figure 1(a). The corresponding Pauli-spin measurements are shown in Figure 1(b). Measurements x j , j = 1, 2, . . . , are taken at discrete time steps and assembled into offset snapshot matrices X and X in Figure 1(c). The bilinear Dynamic Mode Decomposition (biDMD, Figure 1(d)) is a regression-based algorithm that uses the assembled data matrices and control input from sufficiently-resolved data to learn the intrinsic dynamics, A and the control, B, for the bilinear dynamics.
bilinear snapshot matrix where T is the Kronecker product and U X is the column-wise Kronecker product of two matrices (also known as the Khatri-Rao product). The discrete time dynamics of a bilinear dynamical system are well-approximated (for small enough ∆t) [39] by . In its simplest form, the biDMD algorithm is a regression algorithm in the spirit of DMD (13) and DMDc (18) such that  Figure 2(b). In Figure 2(c) a new control is provided using the resonance frequency, ω ≈ 1.001, estimated from the DMD eigenvalues in Figure 2(b). The biDMD model extrapolates using only an initial state and a desired control. Figure 2(c) compares the biDMD extrapolation to the known simulation.
Like DMDc, the biDMD algorithm disambiguates the effect of the intrinsic dynamics, A, from the effect of the control inputs, B, using a factorization of the best-fit solution:

Example 1
Consider a two-level quantum system, H(t) = πσ 3 + u(t)σ 1 , where {σ j } 3 j=1 are the standard Pauli matrices. This can be realized, to take an example from quantum optics, by applying a linearly polarized semi-classical drive to a dipole-allowed atomic transition ∆E = 2π [41]. Suppose the dynamics have been inaccurately characterized so that the system is driven in an experiment by a slightly off-resonant pulse, u(t) = cos(2πω D t) with ω D = 1.1. In Figure 2, the biDMD algorithm uses samples from such an experiment in order to disambiguate the effect of control. The phase of the DMD eigenvalues (15) provide an estimate of the resonance frequency, ω ≈ 1.001. The biDMD model can then be used to predict the behavior of the system under the influence of this estimated-resonance drive without a priori obtaining an accurate characterization.
Recall the trajectory samples are measured via the expected value of an ensemble of identically prepared quantum states. Advantageously, the DMD algorithm is based on a least-squares regression that is optimal with respect to Gaussian measurement noise. The algorithm does not require much data. Moreover, if the algorithm is applied to a system with dissipation, then the fit DMD eigenvalues are simply forced into the unit circle according to the modified drift operator. An otherwise identical result to Figure 2 is obtained without any modification to the framework.

DMD for Stroboscopic Sampling
Recall that the DMD algorithms in Section 3 are applicable only to systems in which the measurement resolution ∆t is small enough to accommodate the linear approximation to the desired accuracy. Higher-order approximations to (10) show that improvements to the biDMD model can be made by adding terms non-linear in the control [39]. This research direction may be appropriate for cases where the zero-order hold on the control is extended for the entirety of the interval ∆t. Instead, in this section we study the case of stroboscopic measurements of the trajectory x(t) separated by time steps T ∆t during which the control can change significantly. The observation of low-frequency dynamics ∼ T in a system under the influence of relatively high-frequency actuation ∼ ∆t is a familiar feature of bilinear systems. For example, the classroom experiment of an inverted pendulum stabilized by the high-frequency vertical drive of its pivot can be modelled as a bilinear control system (under an appropriate operator-theoretic transformation of Mathieu's equation). The two-level quantum system from Section 3.3.1 provides another example. The textbook analysis of the two-level system proceeds via a rotating frame transformation together with the rotating wave approximation (RWA). Consider H(t) = πσ 3 + u(t)σ 1 driven by a pure tone u(t) = u 0 cos(2πνt). Changing to a frame rotating with the drive frequency about the σ 3 axis of the Bloch representation, with U (t) = exp(2πνi(1 − σ 3 )t/2) where the RWA disregards the contribution of the emergent high-frequency terms to realize a constant Hamiltonian. The coefficient of σ x in the RWA is u 0 which can be small relative to 2πν. In this case, the characteristic Rabi cycle of the two-level system is a low-frequency oscillation divorced from the high-frequency drive. Motivated by these examples, Sections 4.1-4.2 discuss DMD strategies for the case of stroboscopic measurements.

DMD and Floquet Theory
In this section, we recall how periodically-driven dynamical systems are connected to the case of stroboscopic measurements (and rotating frames) using the framework of Floquet theory. We introduce Floquet theory and show how corresponding re-interpretations of DMD enable an efficient method for studying stroboscopically-measured, periodically-driven quantum systems.
In Floquet theory (known as Bloch theory in condensed matter physics), the discrete-time evolution of a periodic quantum dynamical system, H(t + T ) = H(t), is given by the timeindependent stroboscopic Floquet Hamiltonian [49], The dependence on the initial time, t 0 , is a gauge transformation known as the Floquet gauge. The Floquet Hamiltonian is a consequence of the observation U (t 0 + nT, t 0 ) = U (t 0 + T, t 0 ) n . A classical version of (26) also exists; for vectorized quantum systems, exchange −iH Take the eigendecomposition of the Floquet operator such that L F [t 0 ]ξ j (t 0 ) = ε j ξ j (t 0 ) for j = 1, 2, . . . , N 2 − 1. Floquet's theorem is the observation that the T -periodic fast-time dynamics can be separated from the slow-time dynamics governed by L F [t 0 ]. The theorem asserts that attaching the periodic fast-time propagation to the eigenvectors, ξ j (t + T ) = ξ j (t), provides a complete set of solutions for the dynamics. In this Floquet basis, the general time-evolution of an initial state x(t 0 ) is given by where the ξ(t) are known as Floquet modes, the ε j are known as quasi-energies, and the coefficients c j are the coordinates of x(t 0 ). Comparing (4) with (27), it is evident the DMD algorithm (13) can provide a framework for the numerical approximation of Floquet modes and quasi-energies using DMD modes and eigenvalues [50]. We refer to this re-interpretation of the DMD algorithm as Floquet DMD. Suppose a size-M measurement record of the trajectory x(t) is obtained from a non-autonomous system with a known period T . Further suppose that the first s measurements are contained within [0, T ), the second within [T, 2T ), and so on under the constraint t r ≡ t ns+r mod T for r = 1, 2, . . . , s. Consider the reshaped snapshot matrices for this stroboscopically-measured data, where X, X ∈ R s(N 2 −1)× (M −1)/s . In this way, the columns of X exist in the span of a discretization of the T -periodic Floquet modes, {ξ j (t)} N 2 −1 j=1 . These Floquet modes evolve with a single phase given by the quasi-energy. From the snapshot matrices, the DMD algorithm constructs an optimal propagator in terms of DMD modes and eigenvalues that provide numerical approximations to the Floquet modes and quasi-energies, respectively.

Example 2
Consider the two-level quantum system from Section 3.3.1 with H(t) = πσ 3 + u(t)σ 1 driven by a control u(t) = cos(2πω D t) with ω D = 1.1. Suppose the system is stroboscopically measured at 4 times per period over a total 4T (such that M = 16 and s = 4). Construct the snapshot matrices (28) from this measurement record. In Figure 3(a,b), the DMD modes and eigenvalues from applying the Floquet DMD algorithm are compared with the exact Floquet modes, exact quasi-energies, the rotating wave approximation (RWA) modes, and the RWA eigenvalues (see Equation 25). Because ω D = 1.1 is close to resonance, the RWA is valid and able to capture the quasi-energies. However, RWA does so by sacrificing the fast-time dynamics. In contrast, the Floquet DMD resolves both the fast and slow time-scales using the stroboscopic measurement data. Moreover, Floquet DMD is not limited to a regime around resonance. The Floquet DMD provides a linear reduced order model that approximates the exact linear dynamics of any periodically-driven quantum dynamical system. This implies the method is well-suited to extrapolation, as shown in Figure 3(c).
The Floquet DMD inherits all the model-reduction advantages of the DMD algorithm. However, one drawback of Floquet DMD is that control is internal to the constructed model. As such, the solution can assist control strategies involving switched systems, but cannot generalize to unseen control inputs.

biDMD and Average Hamiltonian Theory
Average Hamiltonian Theory (AHT) [51,52] is a method with origins in the NMR community that asserts that the dynamics of a quantum dynamical system driven by a time-dependent control Figure 3: A qubit with intrinsic period 1 is driven by a slightly off-resonance frequency with period T = 1/1.1. Equi-spaced stroboscopic measurements are sampled at a frequency of 4/T during a total sampling window of 4T to construct a model using Floquet DMD. Figures 3(a-b) compare Floquet modes and quasi-energies to the eigenmodes and eigenvalues from DMD and the rotating wave approximation (RWA). The columns in Figure 3(a) are indexed by the coordinate of the state x(t) = Tr(ρσ j ). The horizontal axes are labelled by τ ≡ t mod T . Measurement data from a fifth period is used to provide initial conditions for the model-based extrapolation in Figure 3(c).
can be described by the average effect of the field over a period T of its oscillation. In AHT, a known system Hamiltonian is expanded analytically according to the Magnus series [53,54]. In this section, we show how AHT can be combined with the ideas of Floquet DMD from Section 4.1 to provide a framework for biDMD (23) in the case of stroboscopic measurement data. This allows us to extend model-free optimal control to the limit of stroboscopic measurements. We refer to this interpretation of biDMD as AHT-biDMD.
Consider a periodic bilinear dynamical system such that L(t + T ) = L(t).
Define Ω := 2π/T . Recall the existence of the constant Floquet operator L F [t 0 ]x(t 0 ) = x(t 0 + T ). The main idea of this section is to approximate L F [t 0 ] using a highfrequency perturbative expansion which will motivate a way to include control as a parameter in the spirit of biDMD. There are a number of perturbative methods to find the Floquet operator in the  Figure 4(e)) is a regression-based algorithm that uses the assembled data matrices and control input from stroboscopic measurements to learn the intrinsic dynamics, A and the control, B, for the bilinear dynamics.
limit of a high-frequency drive; we will use the familiar Magnus expansion as used in AHT which allows us to represent the constant Floquet operator L F [t 0 ] as a series of constant operators [49], where L (n) means that the operator appears in the series expansion with a pre-factor Ω −n .
By inserting the Fourier series u(t) = K k=1 a k cos(kΩt) + b k sin(kΩt) into each L (j) (30), j = 0, 1, . . . , J − 1, we obtain a bilinear model up to terms of order T −J in which controls enter the Floquet operator as polynomials of the constant Fourier series coefficients. An analytic example is shown in Appendix B. The expansion presumes these coefficients remain constant across a single period [(n − 1)T, nT ] but allows for changes between periods. The bilinearization of the propagator through the basis decomposition of the control is motivated by analytic Floquet-Fourier methods [55][56][57].
To apply AHT-biDMD, construct the snapshot matrices X and X as in (12) using the size-M measurement record of the trajectory x(t) where x n are measured stroboscopically, t n = t 0 + (n − 1)T . Additional intra-stroboscopic measurements can be included using the methods outlined in Section 4.1. Next, suppose controls are restricted to u n (t) = K k=1 a nk cos(kΩt) + b nk sin(kΩt) where n indexes the relevant step [(n − 1)T, nT ] of the applied control. Define a control snapshot in terms of the coefficients,û n = [a n1 , a n2 , . . . , a nK , b n1 , b n2 , . . . , b nK ] T , such that the snapshot matrix is noŵ Furthermore, define the polynomial library matrix to be where e.g.Û P 2 denotes all quadratic non-linearities of the control coefficients, Finally, use the library to construct the bilinear operator as in (21), Θ(Û) X. Now, AHT-biDMD is the approximation of the discrete time dynamics by the model The algorithm proceeds as in biDMD in Section 3.3.

Example 3
Consider again the two-level quantum system with H(t) = πσ 3 + u(t)σ 1 . Restrict the control to the span a k cos(kΩt) + b k sin(kΩt), Ω = 1 2 Figure 5: Figure 5(a) shows the true trajectory compared to the biDMD predictions for a two-level system driven at resonance. Figure 5(b) plots the control over a single period and shows an exact match with the truncated Fourier series whose coefficients make up the control input of the biDMD model. Figure 5(c) shows the error of (a) and additional amplitudes of the same control. Figure 6: Figure 6(a) shows the true trajectory compared to the biDMD predictions for a two-level system driven by a random control. Figure 6(b) plots the control over a single period and shows an exact match with the truncated Fourier series whose coefficients make up the control input of the biDMD model. Figure 6(c) shows the error of (a) and additional amplitudes of the same control.
Leveraging the linearity of the biDMD model, assemble the measurement trajectory data as a horizontal stack of numerical experiments driven by pure tones. Suppose the trajectories are Figure 7: Figure 7(a) shows the true trajectory compared to the biDMD predictions for a two-level system driven by a sawtooth control. Figure 7(b) plots the control over a single period and compares to the approximation via the truncated Fourier series whose coefficients make up the control input of the biDMD model. Figure 7(c) shows the error of (a) and additional amplitudes of the same control.
measured over 5 periods of applied control in which each of a k and b k are taken separately to be 0, 0.1, . . . , 1. Include additive Gaussian noise with σ = 10 −2 . From this control data, construct the library snapshot matrix (32) up to second-order non-linearities. In the trajectories, we also include additional intra-stroboscopic measurements using the methods of Section 4.1. From all of these horizontally-stacked measurement trajectories and the polynomially-extended control coefficients, assemble the snapshot matrices according to (34), and apply the biDMD algorithm to obtain a model for the system. In Figures 5, 6, and 7 we look at the predictive capabilities of the biDMD model over 10 control periods (twice the training time) for certain representative control pulses. Part (b) of each figure shows the input control plotted alongside the Fourier-series signal (35) that best approximates the applied control. Part (c) quantifies the percent error of the measured versus predicted trajectories over a range of amplitudes assuming that the input control pulse remains in a fixed shape.
In particular, Figure 5 shows the case of a resonance drive. The model is successful up to a large amplitude limit. This can be understood by recalling that the validity of the Magnus expansion requires a control period that is "small enough"-a condition that exists in an inverse relationship with the amplitude of the applied drive. Figure 6 shows a successful prediction for the case of an arbitrary multi-frequency control in the span of the truncated Fourier series (for simplicity, we hold the same coefficients across all periods). Finally, Figure 7 shows a sawtooth drive that is only approximately captured by the truncated Fourier series input in the biDMD model; as a result, the model predicts with a reduced accuracy.

Conclusion
The success of quantum optimal control (QOC) is critically dependent on accurate quantum system identification. In cases where theoretical models are inaccurate or unknown, data-driven system identification is essential for practical QOC. Established ideas from machine learning and regression algorithms like the dynamic mode decomposition have advantages that can be useful to quantum technologies. Like the innovation of the bilinear dynamic mode decomposition (biDMD) in this paper, the ideas must first be adapted to accommodate the underlying mathematics and structure of quantum dynamical systems. Much like its forerunner, the DMDc algorithm [37] for the case of direct actuation in classical control systems, biDMD provides a data-driven and equation-free framework for use with QOC strategies that leverages the underlying structure of QOC. We have demonstrated the success of the biDMD formulation on a number of quantum systems, showing that it can successfully leverage time-series data alone to enact control.
The use of the DMD framework for quantum systems brings with it a variety of well-established extensions and improvements beyond the naive algorithms introduced in this paper. Studying the ways these improvements fit into the quantum world provides a myriad of future research directions. First, DMD is by construction a method for reduced-order modelling of high-dimensional systems; however, DMD can also accommodate tensor-network representations of high-dimensional data [39]. DMD can also utilize compressive measurements [35]. As such, biDMD should be studied in systems involving multiple qubits under control. Also, time-delays and their connection to random measurements [58,59] indicate ways DMD can increase system dimensionality in the case of sparse sampling and provide a path for DMD to capture non-Markovian dynamics [60]. Optimization-based DMD algorithms [33,61] can improve characterization of the DMD system by relaxing data-collection strategies and by incorporating additional physical constrains or modelling assumptions. Such strategies can be used to increase the fidelity of DMD-based models in experimental settings. Ultimately, the goal is to demonstrate experimentally that biDMD provides one path to data-driven and equation-free model-based feedback controllers for engineering future quantum technologies. Moreover, it provides a clear connection to Koopman theory for classical system identification and control, illustrating how non-linear classical dynamics are linearized via an infinite-dimensional Hilbert space of observables [29][30][31].

A Vectorization of open quantum dynamics
In quantum dynamics, a time-independent closed quantum system can be described by a Hermitian operator H 0 (the Hamiltonian) acting on a unit vector |ψ(t) (the wavefunction) in a complex Hilbert space H using the Schrödinger equation, In this paper we set = 1, choose H 0 to have zero trace, and assume that H = C N has finite dimension. The control of the system can be realized using control functions, u k (t), coupled to time-independent interaction Hamiltonians, H k , leading to a bilinear Schrödinger equation, More generally, open quantum dynamics involves mixtures of pure quantum states. The statistics of measurements of this kind of state can be completely described by a density operator, ρ, that is any element in the convex set of all non-negative (ρ ≥ 0) self-adjoint (ρ † = ρ) linear maps on H with Tr ρ = 1. The pure state density operator is ρ(t) = |ψ(t) ψ(t)|. The case of Markovian dynamics of ρ can be described by the action of completely-positive trace preserving maps generated by the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation [41], where H(t) is a trace-zero Hermitian operator (the system Hamiltonian), {D j } N 2 −1 j=1 is an orthonormal set of complex matrices with trace zero, and C := (c jk ) is positive semi-definite. For a closed quantum system C = 0, and the equation becomes the quantum Liouville equation.
The Lindbladian, L, is sometimes called a super-operator because of its linear action on the operators ρ (in the Liouville limit, L is instead called the Liouville operator). Several ways to vectorize the GKSL equation exist. We will use the physics language of the Bloch vector (also know as the vector of coherence) to describe passing to a vector ODE. Using the trace-one constraint of a density matrix, write Tr(ρσ j )σ j (39) with {σ j } N 2 −1 j=1 as a complete and orthonormal basis for traceless Hermitian operators. Take x j := Tr(ρσ j ) so that (38) (c) j = i N m,n c mn f mnj (43) where [σ j , σ k ] = i f jk σ and {σ j , σ k } = 2 N δ jk 1 + g jk σ are the structure constants of the basis [4,42]. Identifying H(t) = H 0 + k u k (t)H k with L H = L 0 + k u k (t)L k results in a bilinear GKSL equation for the vectorized density matrix, ∂ ∂t x(t) = (L 0 + k u k (t)L k + L D )x(t) + c .
Similarly, the bilinear quantum Liouville equation is

B Average Hamiltonian Theory example
Recall the two-level quantum dynamical system ( The control appears with the same multiplicity as its attached operator in each non-vanishing commutator. Note that: As such, the desired example expansion is computed to be: