Data-driven resolvent analysis

Resolvent analysis identifies the most responsive forcings and most receptive states of a dynamical system, in an input--output sense, based on its governing equations. Interest in the method has continued to grow during the past decade due to its potential to reveal structures in turbulent flows, to guide sensor/actuator placement, and for flow control applications. However, resolvent analysis requires access to high-fidelity numerical solvers to produce the linearized dynamics operator. In this work, we develop a purely data-driven algorithm to perform resolvent analysis to obtain the leading forcing and response modes, without recourse to the governing equations, but instead based on snapshots of the transient evolution of linearly stable flows. The formulation of our method follows from two established facts: $1)$ dynamic mode decomposition can approximate eigenvalues and eigenvectors of the underlying operator governing the evolution of a system from measurement data, and $2)$ a projection of the resolvent operator onto an invariant subspace can be built from this learned eigendecomposition. We demonstrate the method on numerical data of the linearized complex Ginzburg--Landau equation and of three-dimensional transitional channel flow, and discuss data requirements. The ability to perform resolvent analysis in a completely equation-free and adjoint-free manner will play a significant role in lowering the barrier of entry to resolvent research and applications.


Introduction
The resolvent is a linear operator that governs how harmonic forcing inputs are amplified by the linear dynamics of a system and mapped onto harmonic response outputs. Resolvent analysis refers to the inspection of this operator to find the most responsive inputs, their gains, and the most receptive outputs. The resulting low-rank approximation of the forcing-response dynamics of the full system is extremely valuable for modeling, controlling, and understanding the physics of fluid flows. Interest in the approach has continued to grow since McKeon and Sharma (2010) showed that, by interpreting the nonlinear term in the Fourier-transformed Navier-Stokes equations as an exogenous harmonic forcing, resolvent analysis can uncover elements of the structure in wall turbulence.
Two decades before coined as such, resolvent analysis was first used in the seminal work of Trefethen et al. (1993) to study the response of linearly stable flows to deterministic external disturbances, such as those coming from wall roughness, acoustic perturbations, body forces, or free-stream turbulence. A key result was to identify the non-normality of the linearized operator as the cause of transient energy amplification of disturbances, even for cases deemed as stable by eigenvalue analysis. During the 1990s, nonmodal stability theory emerged to provide a more complete picture of the linear perturbation dynamics for fluid flows using an initial-value problem formulation, as a complement to the eigenproblem from classic hydrodynamic stability theory (Schmid, 2007;Schmid and Brandt, 2014;Schmid and Henningson, 2001). This formulation allowed the study of the response of fluid flows to initial conditions (Butler and Farrell, 1992;Farrell and Ioannou, 1993a;Gustavsson, 1991;Herrmann et al., 2018b;Hwang and Cossu, 2010a;Reddy and Henningson, 1993), stochastic inputs (delÁlamo and Jiménez, 2006;Farrell and Ioannou, 1993b;Hwang and Cossu, 2010a,b), and harmonic forcing (Herrmann et al., 2018b;Hwang and Cossu, 2010a,b;Jovanović and Bamieh, 2005). Another landmark is the framework adopted by Jovanović and Bamieh (2005) that focuses on the response of certain outputs of interest to forcing of specific input components. This input-output viewpoint allows the examination of localized disturbances and provides mechanistic insight into multi-physics systems (Herrmann et al., 2018a;Jeun et al., 2016;Jovanović, 2020), making it particularly relevant for control applications. Reduced-order models based on resolvent modes have been studied for turbulent channel flows (Abreu et al., 2020;McKeon, 2017;Moarref et al., 2013Moarref et al., , 2014, laminar and turbulent cavity flows (Gómez et al., 2016;Sun et al., 2020), and turbulent jets (Lesshafft et al., 2019;. Other exciting recent developments include the design of an airfoil separation control strategy based on resolvent analysis by , the harmonic resolvent formalism to capture cross-frequency interactions in periodic flows by Padovan et al. (2020), and the application of a nonlinear input-output analysis to boundary layer transition by Rigas et al. (2020).
Even though there is a growing interest in the community to study flows with two and three inhomogeneous directions, global resolvent analysis is still far from commonplace. The main reasons for this are the requirement of a high-fidelity solver for the linearized governing equations, and the computational cost and memory allocation associated with handling a very large operator. The latter challenge has been addressed using matrix-free iterative techniques (Bagheri et al., 2009a;Loiseau et al., 2019;Martini et al., 2020;Monokrousos et al., 2010). However, this approach requires having access to a high fidelity solver for the adjoint equations, which adds to the first challenge. A promising alternative, first used by Moarref et al. (2013) and further investigated by Ribeiro et al. (2020), is the use of randomized numerical linear algebra techniques (Halko et al., 2011;Liberty et al., 2007). To simultaneously addresss both issues, we propose a purely data-driven approach to obtain the resolvent operator that doesn't rely on access to the governing equations and can be combined with randomized methods to alleviate the computational expense if needed.
The unprecedented availability of high-fidelity numerical simulations and experimental measurements has lead to incredible growth of research in data-driven modeling of dynamical systems during the past decade (Blanchard and Sapsis, 2019;Li et al., 2020;Loiseau and Brunton, 2018;Raissi et al., 2019Raissi et al., , 2020Rudy et al., 2017). In fluid dynamics, this has led to the development and application of machine learning (ML) algorithms to extract dominant coherent structures from flow data (Brenner et al., 2019;Brunton et al., 2020;Taira et al., 2017. Dynamic mode decomposition (DMD) is a particularly relevant technique introduced by Schmid (2010) to learn spatiotemporal patterns from time-resolved data that are each associated with a single frequency and growth/decay rate . During the last decade, considerable effort has been devoted to interpret its application to nonlinear dynamical systems, based on a deep connection to Koopman theory (Mezić, 2013;Rowley et al., 2009), and to develop numerous extensions to allow the use of non-sequential measurements (Tu et al., 2014), promote sparsity in the solution , work with streaming datasets (Hemati et al., 2014), improve its accuracy with nonlinear observables (Williams et al., 2015), incorporate the effect of control inputs , add robustness using nonlinear optimization (Askham and Kutz, 2018), and enable its application to massive datasets with randomized 2. Assemble data matrices   Figure 1: Schematic of the data-driven resolvent analysis algorithm demonstrated on the transitional channel flow example detailed in §4.2. Data is collected from time recordings of the system of interest, where one or more initial conditions are used to generate transient dynamics. Measurements are stacked into data matrices that are used to approximate an eigendecomosition of the underlying system via dynamic mode decomposition. A projection of the resolvent operator onto the span of the learned eigenvectors is analyzed, and, finally, the produced modes are lifted to physical coordinates. linear algebra techniques (Erichson et al., 2019). Recently, other data-driven modal decompositions have been derived, such as spectral proper orthogonal decomposition (SPOD) (Lumley, 1970;, and the bispectral mode decomposition (BMD) (Schmidt, 2020).  showed that, for turbulent and statistically stationary flows, SPOD modes are equivalent to the output resolvent modes if the nonlinear forcing exhibits no preferential direction. Nonetheless, it is the input resolvent modes that provide insight into the most amplified flow structures, the most sensitive actuator locations, and the most responsive control inputs.
In this work, we present an algorithm to obtain the leading input and output resolvent modes, and the associated gains, of linearly stable flows directly from data by following the procedure shown in Figure 1. The method relies on DMD to approximate the eigenvalues and eigenfunctions of the system based on snapshots from one or more transient trajectories of the flow. We show that, using an appropriate inner-product, the resolvent of the matrix of DMD eigenvalues is the resolvent of the system projected onto the span of the DMD eigenvectors, which are subsequently used to synthesize the resolvent modes in physical coordinates. Our method is able to find the optimal forcing, response, and gain of a dynamical system in an equation-free manner, and without access to data from adjoint simulations, therefore opening up the possibility of resolvent analysis purely based on experimental measurements.
The remainder of the paper is organized as follows. A brief and practical formulation of resolvent analysis is presented in §2, followed by a general description of our proposed method in §3. We demonstrate the approach on two examples and discuss data requirements and convergence in §4. Our conclusions and thoughts on possible extensions are offered in §5.

Resolvent analysis
In this section, we present a brief and practical formulation of resolvent analysis, followed by a description of a dimensionality reduction approach first used in the context of nonmodal stability theory by Reddy and Henningson (1993).
Let us consider a forced linear dynamical systeṁ where the dot denotes time-differentiation, x ∈ C n is the state vector, A ∈ C n×n is the linear dynamics matrix, and f ∈ C n is an external driving force. Such a system may arise from a semidiscretized partial differential equation, and in the case of fluid flows, the incompressible Navier-Stokes equations can be written in this form by projecting the velocity field onto a divergence-free basis to eliminate the pressure variable. The state x may either represent the deviation from a steady state of a laminar flow, or fluctuations about the temporal mean of a statistically stationary unsteady flow. In both cases the matrix A is the linearization of the underlying nonlinear system about the corresponding base flow, either the equilibrium or the mean. However, if we are dealing with an unsteady flow, it is important to note that a Reynolds decomposition yields a nonlinear term that typically cannot be neglected. In this scenario, the nonlinearity is lumped into f and considered as an external forcing, with the caveat that the system may exhibit preferred input directions due to internal feedback between linear amplification and nonlinear interactions (Beneddine et al., 2016;McKeon and Sharma, 2010).
In this work, we will focus on the case where x is the deviation from a stable steady state and f is an exogenous input with no preferential direction, representing disturbances from the environment, model discrepancy, or an open-loop control actuation. For a harmonic forcing f (t) =f e −iωt + c.c., where c.c. means complex conjugate, ω ∈ R is the angular driving frequency and t ∈ R the time variable, the long-term response is also harmonic, x(t) =xe −iωt + c.c., and is governed by the particular solution to (1) given bŷ where I is the n × n identity matrix. Let H(ω) = (−iωI − A) −1 be the matrix approximation of the resolvent operator, where the negative sign accompanying the ω-term follows the convention used for travelling waves. We seek the largest input-output gain, σ 1 (ω), optimized over all possible forcing vectorsf , or more formally where x 2 Q =x * Qx, with ( ) * denoting the Hermitian transpose, measures the size of the state based on a physically meaningful metric given by the positive-definite weighting matrix Q. This weighting accounts for integration quadratures, non-uniform spatial discretizations, and appropriate scaling of heterogeneous variables in multi-physics systems (Herrmann et al., 2018a;Jeun et al., 2016). The Cholesky decomposition is used to factorize Q = F * F and relate the physicallyrelevant norm to the standard Euclidean 2-norm via x 2 Q = Fx 2 2 . With this scaling of the states, and using the definition of the vector-induced matrix norm, the optimal gain optimization problem (3) is equivalent to The solution to (4), along with a hierarchy of optimal and suboptimal forcing and response vectors, is given by the singular value decomposition (SVD) of the weighted resolvent where Σ ∈ R n is a diagonal matrix containing the gains σ 1 ≥ σ 2 ≥ . . . ≥ σ n ≥ 0, also known as singular values, and are unitary matrices whose columns, when left-multiplied by F −1 , yield the input and output resolvent modes, φ j and ψ j , respectively.

Data-driven resolvent analysis
Our data-driven approach to perform resolvent analysis relies on dynamic mode decomposition (DMD) to approximate the eigenvalues and eigenvectors of the underlying dynamical system. Among the many choices of DMD variants, we use the exact DMD approach derived by Tu et al. (2014). In the absence of forcing, the evolution of measurements of the dynamical system of interest (1) is governed by where x k is the measurement at time t k = k∆t, and ∆t is the sampling time. As in the previous section, let the the weight matrix Q = F * F define a physically meaningful norm to quantify the size of the state vector. The following transformation allows us to work in the 2-norm framework, where Θ = F exp(A∆t)F −1 evolves the weighted measurements Fx k one time step into the future. Under this transformation, the adjoint of Θ based on the Q-norm is equivalent to the Hermitian adjoint. Thus, using the weighted measurements, we can proceed using readily available DMD codes. To begin, we collect snapshots of the state denoted by x (j) k , where the subscript k ∈ {1, . . . , m + 1} denotes the sample number, and the superscript j ∈ {1, . . . , p} denotes different trajectories started from p ≥ 1 initial conditions. The next step is to assemble the weighted data matrices where X and Y are of size n × pm. Based on these data matrices, the DMD framework with a rank-r truncated SVD yields the matrices D r = diag(ρ 1 , ρ 2 , . . . , ρ r ) ∈ C r×r containing the approximated eigenvalues, and V r,F ∈ C n×r and W r,F ∈ C n×r , whose columns are the approximated direct and adjoint eigenvectors of the underlying operator Θ. The reader is referred to the work of Tu et al. (2014) for the derivation of the adjoint DMD modes. These eigenvalues and vectors are related to those corresponding to A via λ j = log(ρ j )/∆t, V r = F −1 V r,F , and W r = F −1 W r,F . Hence, we obtain Λ r = diag(λ 1 , λ 2 , . . . , λ r ) ∈ C r×r , V r = [v 1 v 2 . . . v r ] ∈ C n×r , and W r = [w 1 w 2 . . . w r ] ∈ C n×r , that satisfy for an unknown underlying operator A, where A + = Q −1 A * Q is its Q-norm adjoint. In other words, we use DMD as a data-driven eigendecomposition, which is not surprising considering that its connection to Arnoldi methods and Krylov subspaces has been clear since the origins of the algorithm (Schmid, 2010). Next, we seek an approximation of the resolvent operator built on Λ r , V r , and W r . Our approach leverages an operator-based dimensionality reduction technique first used in the context of nonmodal stability analysis by Reddy and Henningson (1993). We now return our attention to the forced system (1), and consider an eigenvector expansion of x and f , as follows where a = [a 1 a 2 . . . a r ] T ∈ C r , and b = [b 1 b 2 . . . b r ] T ∈ C r are the vector of expansion coefficients in eigenvector coordinates. In the work of Reddy and Henningson (1993) V r are the eigenvectors associated with the first r eigenvalues with largest real part of a known operator A, whereas here they are the DMD modes. Substitution of (10) in (1), and taking the inner product with W r at both sides yieldsȧ = Λ r a + b, where we have used the bi-orthogonality property between the sets of direct and adjoint eigenvectors, and assumed that they have been normalized such that W * r QV r = I ∈ R r×r . Because we are now working in different coordinates, if we want to retain the physical meaning of the norm, we need to adjust our inner-product accordingly. The new weighting matrix is derived as where we have defined the new matrixF ∈ C r×r from the Cholesky factorization of V * r QV r =F * F . We are now ready to proceed with the resolvent analysis for the system (10). As presented in the previous section, the weighted resolvent modes and gains are obtained from the SVD of The final step is to synthesize the resolvent modes in physical coordinates, as follows A schematic that summarizes the entire procedure is shown in Figure 1. It is worth pointing out that an analogous procedure can be carried out to obtain data-driven transient growth modes, simply by replacing the resolvent operator with the matrix exponential propagator evaluated at a finite time-horizon. In addition, notice that the reduced-order resolvent matrix in (12) is of size r × r, meaning that its full SVD requires O(r 3 ) operations instead of O(n 3 ), and therefore is considerably cheaper to compute. This projection onto the span of the eigenvectors has been successfully used in operator-based frameworks to achieve computational speedups of a few orders of magnitude for non-modal stability analysis (Herrmann et al., 2018a;Reddy and Henningson, 1993), and several orders of magnitude for pseudospectra computations (Trefethen and Embree, 2005). Although the potential time savings are promising, the computational bottleneck of the presented data-driven method is the truncated SVD in the DMD step, which has an operation count of O(nl 2 m 2 ). However, our approach benefits from all past and future innovations to improve the accuracy, robustness, flexibility, and speed of DMD .

Examples and discussion
In this section, we demonstrate the application of data-driven resolvent analysis on two example problems.

Complex Ginzburg-Landau equation
Our first example is the linearized complex Ginzburg-Landau equation, which is a typical model for instabilities in spatially-evolving flows. The system is governed by the linear operator where x is the spatial coordinate, and D x and D 2 x are the 1 st and 2 nd -order spatial differentiation matrices with homogeneous boundary conditions at x → ±∞. We choose a quadratic spatial dependence for the parameter µ(x) = (µ 0 − c 2 µ ) + µ 2 2 x 2 , that has been used previously by several authors (Bagheri et al., 2009b;Chen and Rowley, 2011;. The other parameters are set to µ 0 = 0.23, µ 2 = −0.01, ν = 2 + 0.4i, and γ = 1 − i, giving rise to linearly stable dynamics. As in Bagheri et al. (2009b), we use spectral collocation based on Gauss-weighted Hermite polynomials to build the differentiation matrices D x and D 2 x and the integration quadrature Q (Weideman and Reddy, 2000). The spatial coordinate is discretized into n = 220 collocation points, and the domain is truncated to x ∈ [−85, 85], which is sufficient to enforce the far-field boundary conditions. Data is generated from 30 simulations that are each started from different initial conditions which we choose to be the first 30 Gauss-weighted Hermite polynomials. We record m = 100 snapshots that are sampled every ∆t = 0.5 time units. Data-driven resolvent analysis is performed using snapshot matrices assembled considering measurements from the first trajectory only. Subsequently, the method is applied on snapshot matrices where measurements from the other simulations are sequentially concatenated one by one, to investigate the convergence behavior in regards to the amount of data required. The Q-norm error between the operator-based and the data-driven leading resolvent modes as a function of the number of trajectories considered is shown in Figure 2. As more data is included, the abrupt drop in the error is expected, since DMD is able to accurately approximate a larger number of eigenvectors, therefore enriching the basis we use to represent the resolvent modes. The first four resolvent modes, as well as their gains as a function the input frequency are shown in Figure 2. For all results presented in this section, we used r = 24 for the rank truncation in both DMD and the data-driven eigenbasis V r . It is worth pointing out that, caution is needed when retaining more vectors in V r to not include spurious eigenvectors, which instead of enriching the basis can be detrimental to the performance of the method.

Transitional channel flow
Our second example is the three-dimensional flow in a plane channel of finite length and depth, and with periodic streamwise and spanwise boundary conditions. The system is governed by the incompressible Navier-Stokes equations, and we consider a Reynolds number Re = 2000 based on the channel half-height and the centerline velocity, and a domain size of 2π × 2 × 2π dimensionless length units along the x, y, and z coordinates that indicate the streamwise, wallnormal, and spanwise directions, respectively. The state vector in this case is composed of the three-dimensional flow field of disturbances about the base parabolic velocity profile.
In the case of single-wavenumber perturbations, the dynamics are described by the traditional Orr-Sommerfeld and Squire equations (Schmid and Henningson, 2001). The corresponding linear operator is built using Chebyshev spectral collocation (Trefethen, 2000) to discretize the wallnormal direction. The Orr-Sommerfeld/Squire operator is then used to compute the operatorbased spectrum and the leading resolvent gain of the three-dimensional flow, as shown in Fig. 1. This is achieved looping over wavenumber combinations that are compatible with the finite channel dimensions, i.e., integers for the current setup. We consider N y = 101 collocation points, and wavenumbers in the range |α| ≤ 7, for the streamwise component, and |β| ≤ 7 for the spanwise component. The operator-based leading resolvent modes, shown in Fig. 3, are computed using ω = 0, α = 0, β = 2, for which the maximum gain is observed to occur. The optimal forcing and response correspond to the familiar streamwise vortices that excite streamwise streaks (Trefethen et al., 1993).
In order to demonstrate our data-driven resolvent analysis, we need snapshots from the transient evolution of the full three-dimensional flow field. We use the spectral code Channelflow (Gibson, 2014;Gibson et al., 2008) to perform direct numerical simulations of the incompressible Navier-Stokes equations. The code uses Chebyshev and Fourier expansions of the flow field in the wall-normal and horizontal directions, and a 3 rd -order Adams-Bashforth backward differentiation scheme for the time integration. We find that a grid with N y = 65 and N x = N z = 32 points is sufficient to discretize the domain for the cases studied, and a time step of 0.01 time units is selected, keeping the CFL number at 0.32. All cases described below are simulated for 500 time units, and snapshots are saved every ∆t = 0.5 time units. The perturbation kinetic energy of all initial conditions simulated is set to 10 −5 to ensure that the effect of nonlinearity is negligible.
In addition to how much data is required for the data-driven resolvent, which was discussed in the previous section, here we investigate the more interesting question of which data. To do so, we learn the leading resolvent modes from three datasets obtained from qualitatively different initial perturbations. First, we consider an initial perturbation of the wall-normal velocity component with the following distribution v(x, y, z, 0) = α,β c α c β (cos(πy) + 1) e i(αx+βz) + c.c., α, β ∈ {−3, ..., 3}, where the y-dependence is selected to satisfy exactly the boundary conditions at the walls, c α and c β are randomly sampled real numbers between −1 and 1, and the horizontal velocity com-  The method is demonstrated using three datasets obtained from DNS initialized with: 1) small-wavenumber random disturbances, where three trajectories are considered, 2) the optimal forcing, and 3) localized actuation.

Operator-based
Operator-based results are also shown for comparison, including the resolvent modes obtained when the input and output are constrained to lie in the span of the snapshots from dataset 3.
ponents are computed from the continuity equation. The first dataset is then composed of the snapshots from three simulations, each initialized with these random small-wavenumber disturbances. Data-driven resolvent analysis is applied retaining r = 200 DMD eigenvectors, and the learned gains and modes are shown in Figure 3. The second dataset considers a single simulation started using the optimal forcing mode as the initial condition. In this scenario, the optimal gain and modes can be learned with only r = 20 DMD eigenvectors, as shown in Figure 3. It is interesting to look at the DMD eigenvalues, which for this case form a small subset of those learned from the first dataset, as shown in Figure 3. In the previous case we learned the eigenvalues with largest real part, now we discover ones from a very specific subset of the complex plane. This highlights that, although the resolvent modes can be accurately represented by very few eigenvectors, the amount of data required to learn those eigenvectors that form an efficient basis is highly-dependent on the dynamic trajectories sampled. Moreover, this opens up exciting research directions, emphasizing the importance of finding principled disturbance designs to effectively probe dynamical systems.
Lastly, the third dataset, considers a localized disturbance as initial condition, which was previously studied by Ilak and Rowley (2008). The exact form of the wall-normal velocity component is v(x, y, z, 0) = 1 − r 2 c 2 r (cos(πy) + 1) e (−r 2 /c 2 where r 2 = (x − π) 2 + (z − π) 2 , and the parameters are set to c r = 0.7 and c y = 0.6. This type of disturbance is close to what could be generated in experiments using a spanwise and streamwise periodic array of axisymmetric jets injecting fluid perpendicular to the wall. Data-driven resolvent analysis is performed using r = 200 DMD eigenvectors, and the resulting optimal forcing and response modes do not resemble the operator-based ones, as shown in Figure 3. In this scenario, the true resolvent modes of the system do not lie on the span of the learned DMD eigenvectors. More importantly, this occurs because some of the eigenvectors required to represent the true resolvent modes are not in the span of the data snapshots, and thus cannot be learned by DMD.
In fact, we show in Figure 3 that the learned resolvent modes for this dataset coincide with the operator-based ones projected onto the span of the data snapshots, which is as good as we can expect to do with a data-driven approach.

Summary and conclusions
In this work, we have developed an algorithm to perform resolvent analysis based on timeresolved data of dynamical systems. Unlike other modal decompositions, resolvent modes provide insight into the most amplified states, the most sensitive actuator locations, and the most responsive control inputs. Our method relies on dynamic mode decomposition to learn approximate eigenvalues and eigenvectors of the underlying linear operator from snapshots of transient dynamics of the system. Subsequently, we are able to compute the resolvent operator projected onto the learned eigenbasis. We perform data-driven resolvent analysis on numerical data of the linearized complex Ginzburg-Landau equation and of disturbances in a three-dimensional plane channel flow, demonstrating agreement between the leading data-driven and operator-based resolvent modes and gain distribution. A critical requirement is the design of initial disturbances to generate transients that are dynamically rich. We show that, using disturbances from a localized actuator, our method recovers the optimal forcing and response of the underlying system projected onto the span of the measured snapshots. This stresses the need for strategies to effectively explore the state space of a dynamical system. The proposed algorithm performs resolvent analysis in an equation-free and adjoint-free manner, therefore opening the possibility of only using experimental measurements. Data-driven resolvent analysis will play a significant role in lowering the barrier of entry to resolvent research and applications. Our results are encouraging for linearly stable flows; however, more work is required before applying this technique to turbulent flows, where the linear and nonlinear contributions to the transient dynamics of the system must be disambiguated.