The majority of model reduction approaches use an efficient representation of the state and then derive equations to temporally evolve the coefficients that encode the state in the representation. In this paper, we instead employ an efficient representation of the entire trajectory of the state over some time interval and solve for the coefficients that define the trajectory on the interval. We use spectral proper orthogonal decomposition (SPOD) modes, in particular, which possess properties that make them suitable for model reduction and are known to provide an accurate representation of trajectories. In fact, with the same number of total coefficients, the SPOD representation is substantially more accurate than any representation formed by specifying the coefficients in a spatial (e.g., POD) basis for the many time steps that make up the interval. We develop a method to solve for the SPOD coefficients that encode the trajectories in forced linear dynamical systems given the forcing and initial condition, thereby obtaining the accurate representation of the trajectory. We apply the method to two examples, a linearized Ginzburg-Landau problem and an advection-diffusion problem. In both, the error of the proposed method is orders of magnitude lower than both POD-Galerkin and balanced truncation applied to the same problem, as well as the most accurate solution within the span of the POD modes. The method is also fast, with CPU time comparable to or lower than both benchmarks in the examples we present.

§ INTRODUCTION

The expense of many modern computational models can prohibit their use in applications where speed is required. In a design optimization problem, for example, many simulations at different boundary conditions or parameters must be performed. In control applications, simulations may need to be conducted in real time to inform actuation. Model reduction techniques strive to deliver the orders-of-magnitude speedup necessary to enable adequately fast simulation for these and other problems with only a mild accuracy sacrifice.

The great majority of model reduction methods employ the following two-step strategy: they (i) find an accurate compression of the state of the system at a particular time and (ii) find equations that evolve the coefficients that represent the state in this representation. The POD-Galerkin method [1, 25, 32], perhaps the most widely used starting point for model reduction, is representative of this approach. The proper orthogonal decomposition (POD) modes are an efficient means of representing the state in that with relatively few POD coefficients, the state can often be represented to high accuracy. In a POD-Galerkin reduced-order model (ROM), these coefficients are then evolved by projecting the governing equations into the space of POD modes, yielding a much smaller dynamical system to evolve. Many alternative choices have been explored for both steps. Examples of the compression step include using balanced truncation modes [24] and autoencoders [17, 12], and examples of deriving the equations in the reduced space include using Petrov-Galerkin projections [3, 4, 27] and learning the equations from data [30, 28]. All of these approaches to model reduction, however, fall within the two-step strategy outlined above.

We have investigated a different approach in this work: instead of representing the state (at a particular time) in a reduced manner, we instead employ a reduced representation for the entire trajectory, i.e., the state's evolution for some time interval. Whereas POD modes are the most efficient (linear) representation of the state, they are far from the most efficient representation of trajectories. This is true intuitively – to represent a trajectory with POD modes, one has to specify the POD coefficients for each time step along the trajectory, but from one time step to the next, the POD coefficients are highly correlated. The analog of POD for entire trajectories is space-time POD [19, 33, 11]. Space-time POD modes are themselves time- and space-dependent, so to represent a trajectory, they are weighted by static coefficients. These modes are the most efficient linear representation of trajectories in the sense that to represent a trajectory to some desired accuracy, fewer degrees of freedom are needed if the trajectory is represented with space-time POD modes than any other linear encoding scheme. Unfortunately, space-time POD modes have a number of characteristics that make them undesirable for model reduction; computing them requires much training data, storing them is memory intensive, and computing space-time inner products, which would be necessary in a space-time ROM method, is expensive.

Fortunately, an efficient space-time basis that does not share the undesirable properties of space-time POD modes exists. Spectral POD (SPOD) modes are most naturally formulated as a POD in the frequency domain. More precisely, at every temporal frequency, there exists a set of spatial modes that optimally represent the spatial structure at that frequency. These modes are the SPOD modes, and they may be thought of as space-time modes where each spatial mode $\boldsymbol{\psi}_{k,j}$ at frequency $\omega_k$ has the time dependence $e^{i\omega_k t}$. Each mode is associated with an energy, and these energies may be compared across frequencies; for example, the second mode at one frequency may have more energy than the first mode at another frequency. The fact that motivates this work is that the most energetic SPOD modes are also an excellent basis for representing trajectories. In fact, SPOD modes converge to space-time POD modes as the time interval becomes long, so for long intervals, the representation of a trajectory with SPOD modes is nearly as accurate (on average) as the space-time POD representation, which is optimal among all linear representations [11].

With this motivation, the goal of this work is to develop an algorithm to solve quickly for the SPOD coefficients that represent a trajectory in forced linear dynamical systems given the initial condition and forcing. If these coefficients can be obtained accurately, then the resulting error will be substantially lower than that of POD-Galerkin with the same number of modes. The method works as follows. The SPOD coefficients at a given frequency are related to the (temporal) Fourier transform of the state at the same frequency, which in turn is related to the forcing and initial condition. We derive these relations analytically and obtain an equation for the SPOD coefficients as a linear operation on the forcing and initial condition. We precompute the linear operators involved, leaving only small matrix-vector multiplications to be done online.

We demonstrate the method on two problems: a linearized Ginzburg-Landau problem with a spatial dimension of $N_x = 220$, and an advection-diffusion problem with $N_x = 9604$. We show that, indeed, we can solve for the SPOD coefficients accurately, resulting in two-orders-of-magnitude lower error than even the projection of the solution onto the same number of POD modes, which is itself a lower bound for the error in any time-domain Petrov-Galerkin method, such as balanced truncation (BT) [24]. We show that this accuracy improvement does not come with an increase in CPU time; the method is competitive with POD-Galerkin and balanced truncation in CPU time, as is predicted by scalings we derive, with a slight advantage in the two examples we present.

Though the space-time approach is uncommon, we are not the first to attempt it [18, 8, 29, 7, 16, 41, 37]. Previous methods have formed both spatial and temporal bases with simulation data. These methods have been used to solve large linear problems as well as small non-linear ones. For example, Ref. [7] solves a large linear Boltzmann transport problem, and Ref. [16] solves an advection-diffusion problem. In the nonlinear case, much of the effort has focused on solving simple nonlinear PDEs over relatively short time intervals [8, 29]. We believe that the representational advantage of SPOD modes relative to previous choices of space-time basis, as well as their analytic time dependence make them a more compelling choice for model reduction.

Our approach may also be related to harmonic balance [15, 14]. In this technique, the governing equations for a temporally periodic system are Fourier-transformed in time, resulting in a set of nonlinear equations to be solved for the transformed fields. This has been applied to the periodic flows arising in turbomachinery [9], resulting in significant computational speedup due to the relatively small number of relevant harmonics. Harmonic balance does not employ a spatial reduction, and our method may be viewed as a spatially reduced harmonic balance method for linear problems. Another important difference is that our method is applicable to non-periodic systems as well as periodic ones.

The remainder of this paper is organized as follows. In Section <ref>, we discuss the properties of POD, space-time POD, and SPOD that are relevant to the method. We present the main approach of the method in Section <ref>. However, there is a subtle issue with transforming to the frequency domain to be accounted for. We describe this issue in Section <ref> and then account for it in the method in Section <ref>. We then demonstrate the method on a linearized Ginzburg-Landau problem and a scalar transport problem in Section <ref>, and conclude the paper in Section <ref>.

§ SPACE-ONLY, SPACE-TIME, AND SPECTRAL POD

We review the space-only, space-time, and spectral forms of POD here. The most significant point for the purposes of this paper is the fact that spectral POD modes approach space-time POD modes as the time interval becomes long, and thus are very efficient in representing trajectories.

§.§ Space-only POD

Space-only POD aims to reconstruct snapshots of the state by adding together prominent modes weighted by expansion coefficients. The first mode is defined to maximize $\lambda[\vc{\phi}({\bf{x}})]$, the expected value of the energy captured by it,

\begin{equation} \label{Sop:maximization}
\lambda[\vc{\phi}({\bf{x}})] =  \frac{\mathbb{E}\big[ |\langle \vc{q}({\vd{x}}), \vc{\phi} ({\vd{x}})\rangle_{\vd{x}} |^2 \big]}{\| \vc{\phi}({\vd{x}})\|^2 }  \text{,}
\end{equation}

\begin{equation} \label{Sop:maximization_with_functional}
\vc{\phi}_1({\vd{x}})  = \argmax \lambda[\boldsymbol{\phi}({\bf{x}})] \text{.}
\end{equation}

The subsequent modes are defined to maximize the energy captured under the constraint that they are orthogonal to all previous ones,

\begin{equation} \label{Sop:maximization_with_functional_latter}
\vc{\phi}_{j}(\vd{x}) = \argmax_{\langle \vc{\phi}({\vd{x}}), \vc{\phi}_{k<j} ({\vd{x}})\rangle_{\vd{x}} = 0} \lambda[\boldsymbol{\phi}({\bf{x}})] \text{.}
\end{equation}

The space-only inner product $\langle \cdot, \cdot \rangle_{\vd{x}}$, which defines the energy captured, is defined as an integral over the spatial domain $\Omega$,

\begin{equation}
\langle \vc{q}({\vd{x}}), \vc{\phi} ({\vd{x}})\rangle_{\vd{x}} = \int_{\Omega} \vc{\phi}^* ({\vd{x}}) \tc{W}(\vd{x}) \vc{q}(\vd{x}) d\vd{x} \text{,}
\end{equation}

where $\tc{W}(\vd{x})$ is a weight matrix used to account for inter-variable importance or possibly to preference certain regions of the domain. One can show [38, 19] that the solution to the optimization problem (<ref>,<ref>) is modes which are eigenfunctions of the space-only correlation tensor,

\begin{equation}
\int_{\Omega}\tc{C}(\vd{x}_1,\vd{x}_2) \tc{W}(\vd{x}_2) \vc{\phi}_j(\vd{x}_2) d\vd{x}_2 = \lambda_{j}\vc{\phi}_j(\vd{x}_1) \text{,}
\end{equation}

where the eigenvalue is equal to the energy of the mode, i.e., $\lambda_j = \lambda[\vc{\phi}_j]$, and the correlation tensor is

\begin{equation}
\tc{C}(\vd{x}_1,\vd{x}_2) = \mathbb{E}[\vc{q}(\vd{x}_1)\vc{q}^*(\vd{x}_2)] \text{.}
\end{equation}

§.§ Space-time POD

Whereas space-only POD modes optimally represent snapshots, space-time POD modes optimally represent trajectories over the time window $[0,T]$. The formulation is much the same as in space-only POD; the modes optimize the expected energy (<ref>,<ref>), but in space-time POD, the inner product is over time as well as space,

\begin{equation}
\langle \vc{q}({\vd{x}},t), \vc{\phi} ({\vd{x}},t)\rangle_{\vd{x},t} = \int_{0}^{T} \int_{\Omega} \vc{\phi}^* ({\vd{x}},t) \tc{W}(\vd{x}) \vc{q}(\vd{x},t) d\vd{x} \ dt \text{.}
\end{equation}

The space-time POD modes that solve this optimization are eigenfunctions of the space-time correlation $\tc{C}(\vd{x}_1,t_1,\vd{x}_2,t_2) =  \mathbb{E}[\vc{q}(\vd{x}_1,t_1)\vc{q}^*(\vd{x}_2,t_2)]$, and the eigenvalues represent the energy (importance) of each mode. The most important property of space-time POD modes for the purpose of this paper is that the space-time POD reconstruction of a trajectory achieves lower error, on average, than the reconstruction with the same number of modes in any other space-time basis. More concretely, using the first $r$ space-time POD modes to reconstruct the trajectory,

\begin{equation}
\vc{q}^r(\vd{x},t) = \sum_{j = 1}^{r} \vc{\phi}_j ({\vd{x}},t) \langle \vc{q}(\vd{x},t) , \vc{\phi}_j ({\vd{x}},t) \rangle_{\vd{x},t} \text{,}
\end{equation}

yields lower expected error

\begin{equation}
\mathbb{E} [\|\vc{q}^r(\vd{x},t) - \vc{q}(\vd{x},t)\|_{\vd{x},t}^2 ]
\end{equation}

than would any other space-time basis. Above, expected error is measured over space and time using the norm $\| \cdot \|_{\vd{x},t}$ induced by the inner product.

§.§ Spectral POD

Spectral POD is most easily understood as the frequency domain variant of space-only POD for statistically stationary systems. In other words, SPOD modes at a particular frequency optimally reconstruct (in the same sense as above) the state at that frequency, on average. The property that makes them attractive for model reduction, however, is that SPOD modes are also the long-time limit of space-time POD modes for statistically stationary systems. These ideas are made precise below, but for a more complete discussion, see [38].

Spectral POD modes at frequency $k$ maximize

\begin{equation} \label{eq:POD:spectral_op}
\lambda_k[\vc{\psi}({\bf{x}})] =  \frac{\mathbb{E}\big[ |\langle \hat{\vc{q}}_k({\vd{x}}), \vc{\psi} ({\vd{x}})\rangle_{\vd{x}} |^2 \big]}{\| \vc{\psi}({\vd{x}})\|^2 }  \text{,}
\end{equation}

again, subject to the constraint that each mode $\vc{\psi}_{k,j}(\vd{x})$ is orthogonal to the previous ones at that frequency $\vc{\psi}_{k,i<j}(\vd{x})$. The Fourier-transformed state is defined as

\begin{equation}
\hat{\vc q}_k = \int_{-\infty}^{\infty} e^{-i\omega_k t}\vc{q}(\vd{x},t) dt \text{.}
\end{equation}

The solution to the optimization (<ref>) is that the modes are eigenvectors of the cross-spectral density
$\tc{S}_{k}(\vd{x}_1,\vd{x}_2) = \mathbb{E}[\hat{\vc q}_k(\vd{x}_1) \hat{\vc q}_k^*(\vd{x}_2)]$

\begin{equation}
\int \limits_{\Omega}\tc{S}_{k}({\vd{x}}_1,{\vc{x}}_2)\tc{W}({\vd{x}}_2)\vc{\psi}_{k,j}({\vd{x}}_2)d{\vd{x}}_2 = \lambda_{k,j}\vc{\psi}_{k,j}({\vd{x}}_1) \text{.}
\end{equation}

The eigenvalue is again equal to the energy of the mode, i.e., $\lambda_{k,j} = \lambda_{k}[\vc{\psi}_{k,j}]$. Modes at frequency $k$ have an implicit time dependence of $e^{i\omega_k t}$.

SPOD modes and their energies become identical to space-time POD modes as the time interval on which the latter are defined becomes long [19, 38, 11]. Thus, the SPOD modes with the largest energies among all frequencies are the dominant space-time POD modes (for long times) and are most efficient for reconstructing long-time trajectories. We denote by $\tilde{\lambda}_j$ the $j$-th largest SPOD eigenvalue among all frequencies, which may be compared to the space-time POD eigenvalues. The convergence in the energy of the trajectory they capture is relatively fast, so for time intervals beyond a few correlation times, the SPOD modes capture nearly as much energy as the space-time modes [11].

If the simulation time of a reduced-order model is long enough for this convergence to be met, the ability of the SPOD modes to capture structures is not diminished relative to that of space-time POD modes. SPOD modes also have two properties that make them more suitable for model reduction than space-time modes: they have analytic time dependence, and they are separable in space and time. The former makes some analytic progress possible in writing the equations that govern the modes and enables Fourier theory to be applied. The latter means that storing the modes requires $N_t$ times less memory, where $N_t$ is the number of times steps in the simulation.

Figure <ref> shows the convergence in the representational ability of space-time POD and SPOD as the time interval becomes long. Specifically, to represent trajectories with some level of accuracy, $98 \%$, say, one needs the same number of SPOD coefficients as space-time POD coefficients if the interval is long compared to the correlation time in the system. This convergence in representation ability occurs because the SPOD modes themselves converge to space-time POD modes in the limit of a long time interval. With space-only POD, one must specify the coefficients for every time step, which leads to a far less efficient encoding of the data because the coefficients are highly correlated from one time step to the next. That SPOD modes are near-optimal in representing trajectories, and that they are substantially more efficient than space-only POD modes motivate this work. If one can efficiently solve for some number of the SPOD coefficients of a trajectory, then these coefficients will lead to substantially lower error than solving for the same number of space-only POD coefficients.

011[tc][tc]Time interval $T$
012[bc][bc]# DOFs for $98 \%$ accuracy

Number of degrees of freedom (DOFs) required to achieve $98\%$ representation accuracy of trajectories as a function of the length of the time interval of the trajectory $[0,T]$. For POD, one must specify all the mode coefficients at every time step, whereas spectral and space-time POD modes are themselves time-dependent. Thus, by leveraging spatiotemporal correlations, fewer DOFs are needed to represent a trajectory to a given accuracy by specifying the SPOD or space-time POD coefficients. As the time interval becomes long, the SPOD and space-time POD modes become equally efficient at representing trajectories.

§.§ Discretization of modes

Upon discretization, the continuous tensors become discrete, and the modes become $N_x$-component vectors in $\mathbb{C}^{N_x}$. Integration over space becomes matrix-vector multiplication. For example, the discrete SPOD modes at frequency $\omega_k$ are defined as the eigenvectors of the (weighted) cross-spectral density matrix,

\begin{equation}
\td{S}_k \td{W} \td{\Psi}_k = \td{\Psi}_k\td{\Lambda}_k \text{.}
\end{equation}

$\td{S}_k$ is the cross-spectral density, $\td{\Psi}_k$ is the matrix with the discrete SPOD modes as its columns, and $\td{\Lambda}_k$ is the diagonal matrix of eigenvalues which are the SPOD mode energies. It is important to note that in practice, the matrix $\td{S}_k$ is not formed; the SPOD modes are calculated by the method of snapshots [36, 38]. In particular, given $r_d$ trajectories from which to obtain SPOD modes, each $N_{\omega}$ time steps in length, the discrete Fourier transform (DFT) of each trajectory is taken. This yields $r_d$ realizations of the $k$-th frequency, for every frequency $k$. These can be formed into a data matrix

\begin{equation}
\td{Q}_k = [\hat{\ket{q}}_k^1, \hat{\ket{q}}_k^2, \dots,  \hat{\ket{q}}_k^{r_d}] \text{,}
\end{equation}

where $\hat{\ket{q}}_k^i \in \mathbb{C}^{N_x}$ is the $k$-th frequency of the DFT of the $i$-th trajectory. The SPOD modes at frequency $\omega_k$ and the associated energies may then be obtained by first taking the singular value decomposition $\td{U}\td{\Sigma}\td{V}^* = 1/\sqrt{r_d}\td{W}^{1/2}\td{Q}_k$. The SPOD modes are then given by $\td{W}^{-1/2}\td{U}$ and the energies by $\td{\Sigma}^2$ [38].

A finite number of evenly spaced frequencies are retained. The lowest one, $\omega_1$, induces a time $T = 2\pi/\omega_1$, which determines the interval $[0,T]$ on which the modes are periodic. The trajectories themselves are, of course, not periodic on this interval, so $T$ is the longest we may use SPOD modes for prediction, though, if a longer prediction is needed, the method may be repeated. Using the SPOD modes, the trajectory may be written

\begin{equation} \label{eq:disc:reconstruction}
\ket{q}(t) = \frac{1}{N_\omega}\sum_{j = 0}^{N_{\omega} - 1} \td{\Psi}_j \ket{a}_j e^{i\omega_j t} \text{,}
\end{equation}

where $\ket{a}_k = \td{\Psi}_k^* \td{W} \hat{\ket{q}}_k \in \mathbb{C}^{N_x}$ is the vector of expansion coefficients at the $k$-th frequency. The factor of $\frac{1}{N_\omega}$ makes (<ref>) an (interpolated) inverse discrete Fourier transform with $\td{\Psi}_k \ket{a}_k$ as the $k$-th coefficient in the DFT of $\ket{q}(t)$. The fully reduced representation of the trajectory is

\begin{equation}
\ket{q}^r(t) = \frac{1}{N_\omega}\sum_{j = 0}^{N_{\omega} - 1} \td{\Psi}^r_j\ket{a}_j^r e^{i\omega_j t} \text{.}
\end{equation}

In this paper, an $r$ superscript indicates a reduced quantity. The mean number of modes retained at each frequency is $r$, and $N_\omega r$ modes are retained in total. The same number is not retained at all frequencies because the most energetic space-time POD modes are the most energetic SPOD modes over all frequencies; the number of modes retained at the $k$-th frequency $r_k$ is the number of modes at this frequency that are among the $N_\omega r$ most energetic overall,

\begin{equation}
r_k = | \{ l : \lambda_{k,l} \geq \tilde{\lambda}_{N_\omega r} \} | \text{.}
\end{equation}

With these notations established, $\ket{a}_k^r \in \mathbb{C}^{r_k}$ and $\td{\Psi}_k^r \in \mathbb{C}^{N_x \times r_k}$. Finally, the Fourier transformed trajectory is defined using the DFT as

\begin{equation}
\hat{\ket{q}}_k = \sum_{j = 0}^{N_\omega - 1} \ket{q}(j \Delta t) e^{-i\omega_k j \Delta t} \text{.}
\end{equation}

§ SPOD PETROV-GALERKIN METHOD

Our goal is to derive a SPOD-based method to solve the linear ordinary differential equation

\begin{equation}
\dot{\ket q}(t) = \td{A} \ket{q}(t) + \td{B}\ket{f}(t) \text{}
\end{equation}

\begin{equation}
\ket{y}(t) = \td{C} \ket{q}(t)
\end{equation}

on the interval $t \in [0,T]$, where $\ket{q}(t) \in \mathbb{C}^{N_x}$, $\ket{f}(t) \in \mathbb{C}^{N_f}$ is some known forcing, $\td{B} \in \mathbb{C}^{N_x \times N_f}$ maps the forcing to the state evolution, $\ket{y}(t) \in \mathbb{C}^{N_y}$ is some linear observable of the state given by $\td{C} \in \mathbb{C}^{N_y \times N_x}$. Given the forcing and initial condition ${\bf q}_0$, our goal is to find the retained SPOD coefficients for the trajectory $\ket{q}(t)$, thereby obtaining the near-optimal rank-$N_\omega r$ space-time representation of the trajectory. With these coefficients, $\ket{y}(t)$ can easily be obtained taking the inverse DFT of $\hat{\ket{y}}_k = \td{C}\td{\Psi}_k^r \ket{a}_k^r$.

§.§ Method

The starting point is the following equation for $\ket{a}^r_k$ in terms of $\hat{ \ket{q}}_k$:

\begin{equation} \label{eq:Uncor:a from qhat}
\ket{a}^r_k = \td{\Psi}^{r*}_k \td{W} \hat{\ket q}_k \text{.}
\end{equation}

This equation is exact and follows from the fact that $\hat{\ket q}_k = \td{\Psi}_k \ket{a}_k$. The Fourier transformed state $\hat{\ket q}_k$ must be obtained from the known forcing and initial condition, and there is some subtlety here. For the sake of clarity in describing the model reduction approach, we ignore the subtlety in this section. The correct relation between the Fourier-transformed forcing and initial condition and the Fourier-transformed state is derived in Section <ref> and is given in (<ref>), and the model reduction method is adjusted for this correction in Section <ref>.

To write the equations in the frequency domain, we take the Fourier transform, replacing $\ket{q}(t)$ and $\ket{f}(t)$ with $\hat{\ket{q}}_k$ and $\hat{\ket{f}}_k$ and (naively, it turns out) replacing $\dot{\ket q}(t)$ with $i\omega_k \hat{\ket{q}}_k$, yielding

\begin{equation}
i\omega_k \hat{\ket{q}}_k = \td{A} \hat{\ket{q}}_k + \td{B}\hat{\ket{f}}_k \text{.}
\end{equation}

Then, assuming the stability of $\td{A}$, the state is related to the forcing by

\begin{equation} \label{eq:Uncorr:full_order_freq}
\hat{\ket{q}}_k = \td{R}_k \td{B}\hat{\ket{f}}_k \text{,}
\end{equation}

where the resolvent operator is defined as

\begin{equation}
\td{R}_k = (i\omega_k \td{I} - \td{A})^{-1} \text{.}
\end{equation}

Plugging (<ref>) into (<ref>), we have

\begin{equation} \label{eq:uncorr_PG_ROM}
\ket{a}^r_k = \td{\Psi}^{r*}_k \td{W}\td{R}_k \td{B}\hat{\ket{f}}_k \text{.}
\end{equation}

This method may be derived as a Petrov-Galerkin method in the following way. Starting from the equations in the frequency domain, $\td{L}_k \hat{\ket q}_k= \td{B} \hat{\ket{f}}_k$, where $\td{L}_k = i\omega_k \td{I} - \td{A}$, and plugging in the representation for $\hat{\ket{q}}_k$, we have

\begin{equation}
\td{L}_k  \td{\Psi}_k^r \ket{a}_k^r= \td{B} \hat{\ket{f}}_k \text{.}
\end{equation}

By left-multiplying the equations by $\td{\Psi}_k^{r*} \td{W}$ and inverting the matrices, the method would be Galerkin [18, 37]. Instead, (<ref>) may be recovered by multiplying the equations on the left by the test basis $\td{\Psi}_k^{r*} \td{R}_k \td{W} $, yielding

\begin{equation}
\td{\Psi}_k^{r*} \td{W} \td{R}_k \td{L}_k  \td{\Psi}_k^r \ket{a}_k^r = \td{\Psi}_k^{r*} \td{W} \td{R}_k \td{B} \hat{\ket{f}}_k \text{.}
\end{equation}

The operator $\td{\Psi}_k^{r*} \td{W} \td{R}_k \td{L}_k  \td{\Psi}_k^r$ is the identity, so (<ref>) is recovered and the method is Petrov-Galerkin.

We denote the matrix that maps the forcing at frequency $k$ to the SPOD mode coefficients at that frequency as $\td{M}_k = \td{\Psi}^{r*}_k \td{W}\td{R}_k\td{B} \in \mathbb{C}^{r_k \times N_f}$. For small problems (i.e., where $\mathcal{O}(N_x^3)$ is not prohibitive) the resolvent, and thus $\td{M}_k$, may be computed directly. More care is needed in larger problems, and, indeed, significant research has been devoted to approximating the resolvent operator in fluid mechanics [10, 21], where it is often used in the context of stability theory [39, 23, 22, 35]. These techniques may be used here to approximate $\td{M}_k$ using some number of the dominant resolvent modes. Implementing these methods, however, may require substantial effort, and an alternative is possible due to the availability of data (the same data used to calculate the SPOD modes).

§.§ Approximating $\td{M}_k$

An accurate and simple-to-implement approximation of $\td{M}_k$ can be obtained by leveraging the availability of data as follows. Defining

\begin{equation}
\ket{g}_k^i  = \td{L}_k\hat{\ket q}_k^i \text{,}
\end{equation}

where, again, $\hat{\ket q}_k^i$ is the $i$-th realization of $\hat{\ket q}_k$ in the training data, we have

\begin{equation}
\hat{\ket q}_k^i = \td{R}_k \ket{g}_k^i \text{.}
\end{equation}

Using the many realizations of the training data (the same ones used to generate the SPOD modes), we have

\begin{equation}
\td{Q}_k = \td{R}_k \td{G}_k \text{,}
\end{equation}

where $\td{Q}_k = [\hat{\ket{q}}_k^1, \hat{\ket{q}}_k^2, \dots] $ and $\td{G}_k = [\ket{g}_k^1, \ket{g}_k^2, \dots] $. Multiplying this equation by its Hermitian transpose gives

\begin{equation}
\td{Q}_k\td{Q}_k^* = \td{R}_k \td{G}_k \td{G}_k^* \td{R}_k^* \text{.}
\end{equation}

Both $\td{Q}_k\td{Q}_k^*$ and $\td{G}_k \td{G}_k^*$ may be represented by their eigendecompositions, i.e., their POD decompositions, giving

\begin{equation} \label{eq:uncorr:first_w_g}
\td{\Psi}_k^r \td{\Lambda}_k^r \td{\Psi}_k^{r*} = \td{R}_k \td{\Psi}_k^{gr} \td{\Lambda}_k^{gr}\td{\Psi}_k^{gr*} \td{R}_k^* \text{,}
\end{equation}

where $\td{\Psi}_k^{gr}$ and $\td{\Lambda}_k^{gr}$ come from the POD of $\td{G}_k$ (the $r$ superscript is a reminder that $\td{\Psi}_k^{gr}$ is not full-rank). By left-multiplying by $\td{\Psi}^{r*}_k \td{W}$ and right-multiplying by $\td{L}^*_k \td{\Psi}_k^{gr} \td{\Lambda}_k^{gr-1} \td{\Psi}_k^{gr*}\td{B}$, we have

\begin{equation} \label{eq:aprxM:almost_there}
\td{\Lambda}_k^r(\td{L}_k \td{\Psi}_k^r )^* \td{\Psi}^{gr}_k \td{\Lambda}_k^{gr-1} \td{\Psi}^{gr*}_k \td{B} = \td{\Psi}_k^{r*} \td{W} \td{R}_k \td{P}_k^{gr} \td{B} \text{,}
\end{equation}

where $\td{P}_k^g = \td{\Psi}_k^{gr} \td{\Psi}_k^{gr*}$ is the projection onto the modes of $\ket{g}_k$. The right-hand-side of (<ref>) is an approximation of $\td{M}_k$, which we label $\td{E}_k$

\begin{equation}
\td{E}_k = \td{\Psi}_k^{r*} \td{W} \td{R}_k \td{P}_k^g \td{B} \approx \td{M}_k \text{.}
\end{equation}

If $\td{P}_k^g$ were the identity, the approximation would be exact. The operator $\td{M}_k$ is used to operate on forcing vectors, so the approximation is accurate to the extent that the forcing is within the span of the realizations of $\ket{g}_k$. The forcing is likely to be near this span because it is closely related to $\ket{g}$, which will become clear in Section <ref>.

§ FULL-ORDER FREQUENCY DOMAIN EQUATIONS

§.§ Problem with the naive equations

The $k$-th Fourier coefficient of the time derivative $\hat{\dot{\ket{q}}}_k$ is not $i\omega_k \hat{\ket{q}}_k$. It is shown in appendix <ref> that in fact,

\begin{equation}
\hat{\dot{\ket{q}}}_k = i\omega_k \hat{\ket{q}}_k + \frac{\Delta \ket{q}} {T} \text{,}
\end{equation}

where $\Delta \ket{q} = \ket{q}(T) - \ket{q}(0)$ is the change of the state on the time interval [20, 26]. If the state were periodic, this term would be zero, and the usual equation would hold, but the solution to the linear ODE has no reason to be periodic. How, then, does one solve in the frequency domain? Modifying (<ref>) with $\frac{\Delta \ket{q}}{T}$ to give

\begin{equation}
\hat{\ket{q}}_k = \td{R}_k \left( \hat{\ket{f}}_k - \frac{\Delta \ket{q}}{T} \right)
\end{equation}

is not useful as $\Delta \ket{q}$ is now known a priori (and, in some sense, is what we are solving for)!

§.§ Solution: derivation of correction formula

To circumvent the problem described above, we use the analytic solution in time of (<ref>). We then take the discrete Fourier transform (DFT) of this solution analytically to obtain $\hat{\ket{q}}$ in terms of $\hat{\ket{f}}$. We use the DFT rather than the Fourier series because, from the exact DFT, one can easily obtain the exact solution at a set of points in the time domain by taking the inverse DFT. On the other hand, with the first $N_\omega$ Fourier series coefficients, point values cannot be recovered exactly.

The analytic solution of (<ref>) is [13]

\begin{equation} \label{eq:Corr:analytic_sol_time}
\ket{q}(t) = e^{\td{A}t}\td{q}(0) + \int_0^{t} e^{\td{A} (t - t')} \td{B} \ket{f}(t') \ dt' \text{.}
\end{equation}

We want the (analytic) DFT of (<ref>) where the $k$-th component of the DFT $\hat{\ket{q}}_k$ of $\ket{q}(t)$ is defined

\begin{equation} \label{DFT}
\hat{\ket{q}}_k = \sum_{j = 0}^{N_\omega-1} \ket{q}_j e^{-\frac{2\pi i}{N_\omega}jk} \text{,}
\end{equation}

where $\ket{q}_j = \ket{q}(j\Delta t)$. Plugging (<ref>) into (<ref>), we have

\begin{equation} \label{DFT lin sys 1}
\hat{\ket{q}}_k = \sum_{j = 0}^{N_\omega - 1} \left( e^{\td{A}j\Delta t}\ket{q}_0 + \int_0^{j \Delta t} e^{\td{A} (j \Delta t - t')}  \td{B} \ket{f}(t') \ dt' \right)  e^{-\frac{2\pi i}{N_\omega}jk} \text{.}
\end{equation}

For later convenience, we refer to the initial condition and forcing terms in (<ref>) as $\hat{\ket{q}}_{k,ic}$ and $\hat{\ket{q}}_{k,force}$, respectively.

First, we evaluate $\hat{\ket{q}}_{k,ic}$. It may be rewritten as

\begin{equation} \label{just IC term DFT}
\hat{\ket{q}}_{k,ic} = \sum_{j = 0}^{N_\omega - 1}  e^{(\td{A} - i \omega_k)j\Delta t} \ket{q}_0 \text{,}
\end{equation}

where we have defined $\Delta t = \frac{T}{N_\omega} $ and substituted $\omega_k = \frac{2\pi k}{T}$. The key to evaluating (<ref>) is to notice that it is a matrix geometric sum, i.e., each term is the previous one multiplied by $e^{(\td{A} - i \omega_k)\Delta t}$. The solution, entirely analogous to the scalar geometric sum, is

\begin{equation}
\hat{\ket{q}}_{k,ic} = (\td{I} - e^{(\td{A} - i \omega_k)\Delta t})^{-1} (\td{I} - e^{(\td{A} - i \omega_k)N_\omega\Delta t}) \ket{q}_0 \text{.}
\end{equation}

Because $e^{i\omega_k N_\omega \Delta t} = 1$, this simplifies to

\begin{equation}
\hat{\ket{q}}_{k,ic} = (\td{I} - e^{(\td{A} - i \omega_k)\Delta t})^{-1} (\td{I} - e^{\td{A}T})\ket{q}_0 \text{.}
\end{equation}

We view the first factor $ (\td{I} - e^{(\td{A} - i \omega_k)\Delta t})^{-1}$ as a time-discrete resolvent operator; expanding the matrix exponential to first order in $\Delta t$, one recovers a multiple of the resolvent. Note, however, that truncating the expansion to first order is always invalid for the higher frequencies because, for these frequencies, $\omega_k \Delta t$ is order unity.

To evaluate $\hat{\ket{q}}_{k,force}$, we assume that $\ket{f}(t)$ can be written as

\begin{equation}
\ket{f}(t) = \frac{1}{N_\omega} \sum_{l = 0}^{N_\omega - 1}\hat{\ket{f}}_l e^{i\omega_l t} \text{.}
\end{equation}

Though this is likely not exactly true in practice, it introduces minimal error (as we show in our numerical examples), and is necessary for evaluating the integral in (<ref>). After pulling out the constant matrix exponential term and including the Fourier interpolation of $\ket{f}(t)$, the integral in (<ref>) is

\begin{equation}
\hat{\ket{q}}_{k,force} = \frac{1}{N_\omega} \sum_{j = 0}^{N_\omega - 1} e^{-i\omega_k j \Delta t} e^{\td{A} j \Delta t} \int_0^{j \Delta t} \sum_{l = 0}^{N_\omega - 1}  e^{(i\omega_l - \td{A})t'} \td{B} \hat{\ket{f}}_l \ dt' \text{.}
\end{equation}

Integrating (and hence inverting the matrix exponential within the integral) brings out a resolvent, and the integral evaluates to

\begin{equation}
\hat{\ket{q}}_{k,force} = \frac{1}{N_\omega} \sum_{j = 0}^{N_\omega - 1} e^{-i\omega_k j \Delta t} \sum_{l = 0}^{N_\omega - 1} \td{R}_l \left( e^{i\omega_l j \Delta t} - e^{\td{A} j \Delta t} \right) \td{B} \hat{\ket{f}}_l \text{,}
\end{equation}

where, again, $\td{R}_l = (i\omega_l\td{I} - \td{A})^{-1}$ is the resolvent operator at the $l$-th frequency. Now, using this result, the forcing term in (<ref>) is

\begin{equation}
\hat{\ket{q}}_{k,force} = \frac{1}{N_\omega} \sum_{j = 0}^{N_\omega - 1}  \sum_{l = 0}^{N_\omega - 1} \td{R}_l \left( e^{i(\omega_l - \omega_k) j \Delta t} - e^{(\td{A} - i\omega_k \td{I}) j \Delta t} \right) \td{B} \hat{\ket{f}}_l
\end{equation}

The frequency difference term evaluates to zero for $\omega_l \neq \omega_k$, and the other term may be evaluated by the same geometric sum argument as before. The entire forcing term in (<ref>) becomes

\begin{equation}
\hat{\ket{q}}_{k,force} = \td{R}_k \td{B} \hat{\ket{f}}_k - \frac{1}{N_\omega} (\td{I} - e^{(\td{A} - i \omega_k)\Delta t})^{-1} (\td{I} - e^{\td{A}T}) \sum_{l = 0}^{N_\omega - 1} \td{R}_l \td{B} \hat{\ket{f}}_l \text{.}
\end{equation}

Adding $\hat{\ket{q}}_{k,ic}$ and $\hat{\ket{q}}_{k,force}$, we obtain the following equation for the $k$-th component of the DFT of the state

\begin{equation} \label{solution}
\hat{\ket{q}}_k = \td{R}_k \td{B}  \hat{\ket{f}}_k + \left(\td{I} -e^{(\td{A} - i \omega_k)\Delta t} \right)^{-1} \left(\td{I} - e^{\td{A}T} \right) \left(\ket{q}_0 - \frac{1}{N_\omega} \sum_{l = 0}^{N_\omega - 1} \td{R}_l \td{B} \hat{\ket{f}}_l \right) \text{.}
\end{equation}

The first term is the naive one. For a stable system, it may be shown that the correction term is inversely proportional to the length of the time interval. Therefore, while the error from ignoring it does decrease as the time interval increases in this case, it is not negligible unless the interval is extremely long. Intuitively, the correction term accounts for the transient in the solution, which decays exponentially at the rate prescribed by the slowest decaying eigenvector of $\td{A}$. The decay time is (of course) independent of the time interval $T$ on which we take the transform, so the fraction of the time interval occupied by this decay is inversely proportional to $T$. The term $\frac{1}{N_\omega}\sum_{l = 0}^{N_\omega - 1} \td{R}_l \hat{\ket{f}}_l$ is the initial condition that is `in sync' with the forcing. If $\ket{q}_0$ were equal to this, then there would be no transient, $\ket{q}(t)$ would be exactly periodic on the interval, and the naive equation (<ref>) would be correct. For long time intervals relative to the slowest decaying eigenvector of $\td{A}$, the term $(\td{I} - e^{\td{A}T})$ may be ignored. We also note that summing the operator multiplying the initial condition and forcing sum over all frequencies can be shown to give $N_\omega \td{I}$

\begin{equation}
\sum_{k = 0}^{N_\omega - 1} \left(\td{I} -e^{(\td{A} - i \omega_k)\Delta t} \right)^{-1} \left(\td{I} - e^{\td{A}T} \right) = N_\omega \td{I} \text{.}
\end{equation}

Error from uncorrected and corrected full-order frequency domain solutions for a stable system (a) and unstable system (b). In the stable case, the error is due to the transient, and the decay is determined by the largest eigenvalue of the system. In both cases, the corrected equations are substantially better. If the forcing is exactly periodic, there is only machine precision error for the solution to the corrected equations.

Figure <ref> shows the inadequacy of the uncorrected method for capturing trajectories in initial value problems. The results shown are for the Ginzburg-Landau system described in Section <ref>, with a stable (a) and unstable (b) value of $\mu_0$. While the corrected frequency domain equations (<ref>) are correct regardless of the stability of the system, we only recommend the method in the stable case; eigensystem methods are likely superior in the unstable case. The remaining error in the corrected method is due to the fact that the forcing is not periodic, and is not composed solely of retained frequencies, meaning that its Fourier series representation used above is approximate. However, this error is indeed small in both the stable and unstable cases.

§ CORRECTED METHOD

The Petrov-Galerkin method is the full-order frequency domain equations (<ref>) left-multiplied by $\td{\Psi}^{r*}_k \td{W}$,

\begin{equation} \label{eq:corr_pgp_toapprox}
\ket{a}_k = \td{\Psi}^{r*}_k \td{W} \td{R}_k \td{B} \hat{\ket{f}}_k +  \td{\Psi}^{r*}_k \td{W} \left(\td{I} -e^{(\td{A} - i \omega_k \td{I})\Delta t} \right)^{-1} \left(\td{I} - e^{\td{A}T} \right) \left(\ket{q}_0 - \frac{1}{N_\omega} \sum_{l} \td{R}_l \td{B} \hat{\ket{f}}_l \right) \text{.}
\end{equation}

As in the uncorrected case, for small systems, the inverses and matrix exponentials above may be computed directly, but this is not possible for larger systems where $\mathcal{O}(N_x^3)$ operations are infeasible. The first term is the same as in the uncorrected case, so it may be approximated with $\td{E}_k$, as before. We approximate the latter terms at each frequency by using the basis of SPOD modes at that frequency, as described below.

§.§ Operator approximations

We define the number of SPOD modes that are to be used in approximating the operators as $r_d$, and we assume this is the same for each frequency, though this is not necessary. At a maximum, this number is the number of SPOD modes available in the data, but it can be fewer if desired. To approximate the operators accurately, it will likely be larger than the number of modes retained to represent the state, $r_d > r_k$ for all $k$. To accomplish the approximations, we use the available $r_d$ SPOD modes as follows. To approximate the inverse $\left(\td{I} -e^{(\td{A} - i \omega_k \td{I})\Delta t} \right)^{-1}$ and matrix exponential $ \left(\td{I} - e^{\td{A}T} \right)$ terms, we first multiply by the identity in various places,

\begin{equation} \label{eq:corr_pgp:intermediate1}
\td{\Psi}^{*}_k \left(\td{I} -e^{(\td{A} - i \omega_k \td{I})\Delta t} \right)^{-1}  \left(\td{I} - e^{\td{A}T} \right) = \td{\Psi}^{*}_k \left(\td{I} -e^{(\td{A} - i \omega_k)\Delta t} \right)^{-1} \td{\Psi}_k \td{\Psi}^{*}_k \td{W} \left(\td{I} - e^{\td{A}T} \right) \td{\Psi}_k \td{\Psi}^{*}_k \td{W} \text{.}
\end{equation}

Note that, at this point, $\td{\Psi}_k$ is full rank, so $\td{\Psi}_k \td{\Psi}^{*}_k \td{W}$ is the identity. It is straightforward to show that (<ref>) may be rewritten with the SPOD bases inside the inverse and matrix exponentials as

\begin{equation} \label{eq:corr_pgp:intermediate2}
\td{\Psi}^{*}_k \left(\td{I} -e^{(\td{A} - i \omega_k \td{I})\Delta t} \right)^{-1}  \left(\td{I} - e^{\td{A}T} \right) = \left(\td{I} -e^{(\td{\Psi}^{*}_k \td{W} \td{A}\td{\Psi}_k - i \omega_k \td{I})\Delta t} \right)^{-1}  \left(\td{I} - e^{\td{\Psi}^{*}_k \td{W} \td{A}\td{\Psi}_kT} \right)  \td{\Psi}^{*}_k \td{W} \text{.}
\end{equation}

Finally, truncating the operators in (<ref>), i.e., $\td{\Psi}_k \to \td{\Psi}_k^{r_d}$, and denoting $\tilde{\td{A}} = \td{\Psi}^{r_d*}_k \td{W} \td{A}\td{\Psi}^{r_d}_k$, the approximated term is

\begin{equation}
\td{\Psi}^{r*}_k \td{W} \left(\td{I} -e^{(\td{A} - i \omega_k \td{I})\Delta t} \right)^{-1} \left(\td{I} - e^{\td{A}T} \right) \approx \td{P}_k \left(\td{I} -e^{(\tilde{\td{A}}_k - i \omega_k \td{I})\Delta t} \right)^{-1} \left(\td{I} - e^{\tilde{\td{A}}_kT} \right) \td{\Psi}^{r_d*}_k \td{W} \text{.}
\end{equation}

Here, $\td{P}_k = \begin{bmatrix}
\td{I}_{r_k} && \td{0}
\end{bmatrix} \in \mathbb{R}^{r_k \times r_d}$ selects the first $r_k$ rows of the matrix it multiplies and $\td{\Psi}_k^{r_d}$ is the first $r_d$ columns of $\td{\Psi}_k$. As desired, all matrix exponentials and inverses are of size $r_d \times r_d$, which makes them tractable.

The forcing sum in (<ref>), which is difficult to compute directly because it involves resolvents, must also be approximated. In the unreduced equations, this term is the same at each frequency, so the sum over frequencies only needs to be computed once. Any approximation of this term should be the same for each frequency to avoid quadratic scaling in $N_\omega$. The natural choice for approximating the forcing sum is to approximate each term by

\begin{equation}
\td{R}_l \td{B} \hat{\ket{f}}_l \approx \td{\Psi}_l^{r}\td{\Psi}_l^{r*} \td{W} \td{R}_l \td{B} \hat{\ket{f}}_l \text{.}
\end{equation}

This approximation is accurate because the SPOD modes at the $l$-th frequency are the best basis for $\hat{\ket{q}}_l$ and are thus a very good basis for $\td{R}_l\hat{\ket{f}}_l$, which is $\hat{\ket{q}}_l$ without the correction. The operator $\td{\Psi}_l^{r*} \td{W} \td{R}_l \td{B}$ may again be approximated with $\td{E}_l$.

The equations can, at this point, be written as

\begin{equation} \label{eq:simpler_pgp_noib}
\ket{a}_k = \td{E}_k \hat{\ket{f}}_k + \td{F}_k  \left(\ket{q}_0 - \frac{1}{N_\omega} \sum_{l} \td{\Psi}_l^r \td{E}_l \hat{\ket{f}}_l \right) \text{,}
\end{equation}

where $\td{F}_k = \left(\td{I} -e^{(\tilde{\td{A}}_k - i \omega_k \td{I})\Delta t} \right)^{-1} \left(\td{I} - e^{\tilde{\td{A}}_kT} \right) \td{\Psi}^{r*}_k \td{W} \in \mathbb{C}^{r_k \times N_x}$. The operators have been approximated, so far, to avoid  $\mathcal{O}(N_x^3)$ scaling in the offline phase of the algorithm. To avoid $\mathcal{O}(N_x)$ scaling in the online phase, one final approximation must be made. The term in parentheses in (<ref>) is multiplied on the left by $\td{F}_k$ for every frequency $\omega_k$, leading to $N_x N_\omega r$ scaling. This can be avoided by storing the term in parentheses in (<ref>) in a rank-$p$ reduced basis $\td{\Phi} \in \mathbb{C}^{N_x \times p}$ and precomputing the product of this basis with each $\td{F}_k$. This basis should represent the initial condition and forcing sum terms accurately, and in practice, we choose POD modes of the state. With this approximation, the corrected SPOD Petrov-Galerkin method becomes

\begin{equation} \label{eq:corr_pgp:final_eq}
\ket{a}_k = \td{E}_k \hat{\ket{f}}_k + \td{H}_k   \left(\td{\Phi}^*\td{W} \ket{q}_0 - \frac{1}{N_\omega} \sum_{l} \td{T}_l \td{E}_l \hat{\ket{f}}_l \right) \text{,}
\end{equation}

where $\td{H}_k = \td{F}_k \td{\Phi} \in \mathbb{C}^{r_k \times p}$ and $\td{T}_l = \td{\Phi}^* \td{W} \td{\Psi}_l^r \in \mathbb{C}^{p \times r_l}$.

§.§ Formal statement of the algorithm and scaling

The offline and online phases of the SPOD Petrov-Galerkin method are shown in Algorithms <ref> and <ref>, respectively.

SPOD Petrov-Galerkin method (offline)
[1]

Inputs: $\{ \td{\Psi}^{r_d}_i \}$, $\{ \td{\Lambda}^{r_d}_i \}$, $\{ \td{\Psi}^{gr}_i \}$, $\{ \td{\Lambda}^{gr}_i \}$, $\td{\Phi}$, $\{ r_i \}$, $\td{A}$, $\td{B}$, $\td{C}$, $\td{W}$
$k \in \{1,2,\dots, N_{\omega} \}$
$\td{\Psi}_k^{r} \gets [\ket{\psi}_{k,1}, \dots , \ket{\psi}_{k,r_k}]$ Retained SPOD modes
$\td{L}_k \gets i\omega_k \td{I} - \td{A}$ Inverse of resolvent

$\td{E}_k \gets (\td{\Lambda}_k^r (\td{L}_k \td{\Psi}_k^r)^* \td{\Psi}_k^{gr} ) (\td{\Lambda}_{k}^{gr-1} \td{\Psi}_k^{gr*}) \td{B}$ Precomputation of first operator

$\tilde{\td{A}}_k \gets \td{\Psi}^{r_d*}_k \td{W} \td{A} \td{\Psi}_k^{r_d}$ Reduced $\td{A}$ to compute matrix exponentials
$\td{P}_k \gets \begin{bmatrix}
\td{I}_{r_k} && \td{0}
\end{bmatrix}$ Row selector matrix
$\td{H}_k \gets \td{P}_k (\td{I} - e^{(\tilde{\td{A}}_k - i\omega_k \td{I})\Delta t} )^{-1} (\td{I} - e^{\tilde{\td{A}}_k T}) (\td{\Psi}_k^{r_d*} \td{W} \td{\Phi}$) Precomputation of second operator
$\td{T}_k \gets \td{\Phi}^* \td{W} \td{\Psi}^r_k$ Precomputation of third operator
$\td{C}^{\td{\Psi}}_k \gets \td{C}\td{\Psi}^r_k$ $\td{C}$ in SPOD basis

Inputs: $\{ \td{\Psi}^{r_d}_i \}$, the $r_d$ SPOD modes at each frequency; $\{ \td{\Lambda}^{r_d}_i \}$, the $r_d$ SPOD energies at each frequency; $\{ \td{\Psi}^{gr}_i \}$, the POD modes of $\ket{g}_i$ for each frequency (see (<ref>)); $\{ \td{\Lambda}^{gr}_i \}$, the energies for the POD modes of $\ket{g}_i$ for each frequency; $\td{\Phi}$, the basis for reducing the initial condition and forcing terms, $\{ r_i \}$, the number of modes to be kept at each frequency; $\td{A}$, $\td{B}$, $\td{C}$, the system matrices; $\td{W}$, the weight matrix.

Outputs: $\{ \td{E}_i \}$, $\{ \td{H}_i \}$, $\{ \td{T}_i \}$, the operators in (<ref>) for each frequency; $\{ \td{C}^{\td{\Psi}}_i \}$, the operators that map the SPOD mode coefficients to $\ket{y}$ for each frequency.

SPOD Petrov-Galerkin method (online)
[1]
Input parameters: $\{ \td{\Psi}^r_i \}$, $\td{\Phi}$, $\td{W}$, ${\ket f}$, $\ket{q}_0$, $\{ \td{E}_i \}$, $\{ \td{H}_i \}$, $\{ \td{T}_i \}$, $\{ \td{C}^{\td{\Psi}}_i \}$
$\hat{\ket f} \gets \mathtt{FFT}({\ket f})$ FFT of forcing
$\ket{a}^\td{\Phi}_0 = \td{\Phi}^* \td{W} \ket{q}_0$ Reduced initial condition
$\tilde{\ket{a}}^{\td{\Phi}}_0 \gets \td{0}_{p \times 1}$ Initializing forcing sum term
$k \in \{1,2,\dots, N_{\omega} \}$

$\ket{b}_k \gets \td{E}_k \hat{\ket f}_k$
$\tilde{\ket{a}}^{\td{\Phi}}_0  \gets \tilde{\ket{a}}^{\td{\Phi}}_0  + \frac{1}{N_\omega}\td{T}_k \ket{b}_k$ Forcing sum

$k \in \{1,2,\dots, N_{\omega} \}$

$\ket{a}^r_k \gets \ket{b}_k + \td{H}_k(\ket{a}^\td{\Phi}_0 - \tilde{\ket{a}}^{\td{\Phi}}_0 )$ Assigning SPOD coefficients
$\hat{\ket y}^r_k \gets \td{C}^{\td{\Psi}}_k \ket{a}^r_k$ Constructing observable in frequency domain

$\ket{y}^r \gets \mathtt{IFFT}(\hat{\ket y}^r)$ Observable in time domain

Inputs: $\{ \td{\Psi}^r_i \}$, the retained SPOD modes for each frequency; $\td{\Phi}$, the basis for reducing the initial condition and forcing terms; $\td{W}$, the weight matrix; ${\ket f}$, the forcing as a function of time, $\ket{q}_0$, the initial condition; $\{ \td{E}_i \}$, $\{ \td{H}_i \}$, $\{ \td{T}_i \}$, the operators in (<ref>) for each frequency; $\{ \td{C}^{\td{\Psi}}_i \}$, the operators that map the SPOD mode coefficients to $\ket{y}$ for each frequency.

Output: $\ket{y}^r$, the observable in the time domain.

So long as the offline time is feasible, which we show next, the online scaling is the salient cost. In calculating the online cost, we count the operations necessary to go from the time domain forcing and initial condition to the SPOD coefficients. In practice, $\ket{y}$ is likely small, thus multiplying the SPOD coefficients by $\td{C} \td{\Psi}_k$ and taking the inverse DFT contributes insignificantly to the scaling. The complexity of the online algorithm is

\begin{equation} \label{eq:alg:scaling}
\mathcal{O} ((r p + r N_f + N_f \log N_\omega ) N_\omega) \text{.}
\end{equation}

This time should be compared with the POD-Galerkin complexity, which is $\mathcal{O}( (rN_f + r^2 )N_t)$. The two are similar, and the differences are due to how $p$ compares with $r$, how $N_t$ compares with $N_\omega$, and the constants involved. In our numerical experiments, we find that the SPOD method is slightly faster for equal rank (but much more accurate). In Appendix <ref> we detail a further approximation that removes the $N_f$ scaling from (<ref>) by employing the discrete empirical interpolation method (DEIM) [5] to approximate the forcing, and other $N_x$-tall vectors, via sparse samplings of them. This can lead to significant speed-up in cases where $N_f$ is large.

In Algorithm <ref>, we have assumed that the SPOD modes for the $\hat{\ket{q}}$ and $\hat{\ket{g}}$ have already been obtained. These modes can be obtained using the techniques described in, e.g., Refs. [38, 34], and the cost for this step is $\mathcal{O}(r_d^2N_xN_\omega)$, where again $r_d$ is the number of SPOD modes obtained from the data, which is the same as the number of temporal blocks formed in Welch's algorithm. If the full-order model (FOM) is sparse, the scaling of the remaining steps is within the time required to obtain the modes, so, not including the time to generate the FOM data, the offline scaling is  $\mathcal{O}(r_d^2N_xN_\omega)$. If the FOM is dense, the scaling is $\mathcal{O}(r_d^2N_xN_\omega + r_dN_x^2N_\omega)$.

§ EXAMPLES

Here, we demonstrate the proposed method on two examples, a linearized Ginzburg-Landau problem and a scalar transport problem. The former is a dense system of dimension $N_x = 220$, and the latter is a sparse system of dimension $N_x = 9604$. For the Ginzburg-Landau system, we compare the accuracy and cost of the proposed method to those of POD-Galerkin and balanced truncation. The error is roughly two orders of magnitude lower for the proposed method than the other two (depending on the case), and the CPU time is similar for all methods. For the scalar transport case, the offline time for balanced truncation makes the method (in its unapproximated form) infeasible, so we compare only to POD-Galerkin. The proposed method is again orders of magnitude more accurate at similar CPU cost.

§.§ Linearized Ginzburg-Landau problem

In continuous space, the complex linearized Ginzburg-Landau equation is

\begin{equation} \label{eq:GL_def}
\dot{q}(x,t) =  \mathcal{A}q(x,t) + f(x,t) \text{,}
\end{equation}

\begin{equation}
\mathcal{A} = -\nu \frac{\partial}{\partial x} + \gamma \frac{\partial^2}{\partial x^2} + \mu(x) \text{,}
\end{equation}

and $f(x,t)$ is a forcing. Following Ref. [2], we set $\nu = 2 + 0.2i$ and $\gamma = 1 - i$. The parameter $\mu(x)$ takes the form

\begin{equation}
\mu(x) = (\mu_0 - c_\mu^2) + \frac{\mu_2}{2}x^2 \text{,}
\end{equation}

with $c_\mu = 0.2$, $\mu_2 = -0.01$ [2], and with $\mu_0 = 0.229$ [38]. The system can be interpreted as an advection-diffusion equation with a local exponential term. The equation supports traveling wave behavior in the positive $x$-direction and is stable in the sense that all the eigenvalues of the linear operator (discretized or continuous) are negative, so all solutions to the unforced equations decay asymptotically. Whether the exponential term promotes local growth or local decay depends on the sign of $\mu(x)$. With the parameters used, $\mu(x)$ is positive when $x \in [-6.15,6.15]$ and negative elsewhere, so as waves move through this region, they grow substantially before decaying once again after passing through it.

017[bc][bc]Gaussian forcing
019[bc][bc]White forcing
023[bc][bc]Gaussian-forced state
025[bc][bc]White-forced state

000[ct][ct]$\matlabfragNegXTick 50$
003[ct][ct]$\matlabfragNegXTick 50$

Ginzburg-Landau state and forcing trajectories. (a) the state resulting from the white forcing in (c). (b) the state resulting from the Gaussian forcing in (d). Both forcings have the same spatial correlation, but the short temporal correlation in the white forcing leads to more jagged structures in the corresponding state. Both trajectories consist of waves traveling in the positive $x$-direction that are amplified in a region near $x=0$.

When the equation is discretized in space, it takes the general form in (<ref>). Following [2, 6], we use a pseudo-spectral Hermite discretization with $N_x = 220$ collocation points [38], and solve the discretized equations using MATLAB's . This full-order model is run for $12000$ time steps with $\Delta t = 0.2$, generating training data from which to educe the modes. We use the procedure described in [38], segmenting the single long trajectory into $142$ (overlapping) trajectories, each of length $N_\omega = 1024$ time steps.

Figure <ref>(a,b) shows two space-time trajectories of the state $q$, each $1024$ time steps in length. The diagonally oriented structures demonstrate the traveling wave behavior of the system, and it is clear that the waves are amplified and then attenuated as they pass through $x=0$. These space-time trajectories are to space-time POD (and thus SPOD) as snapshots are to POD: the more structure there is in the trajectories, the fewer space-time modes are needed to accurately represent them. For example, in Figure <ref>(a) the state is forced with band-limited temporally white noise and a Gaussian spatial correlation while the state in Figure <ref>(b) is forced with a temporally, as well as spatially, Gaussian noise. The resulting state from the white forcing has more detailed structures – a good proxy for higher rank behavior – so trajectories with this forcing require more space-time modes to be accurately represented. The temporally white and temporally Gaussian forcings are shown in Figure <ref>(c,d). Note that the forcing occupies the entire domain in this example, i.e., $N_f = N_x$. We thus choose $p = N_x$, because this will not have a large impact on the CPU time in this case. The SPOD method, therefore, scales like $N_x$, as do POD-Galerkin and balanced truncation, and there is no scaling benefit gained by using an intermediary basis, so we set it to the identity in this example. Despite this spatially extensive forcing, all three ROMs maintain a significant advantage over the FOM due to the latter being dense owing to its pseudo-spectral discretization (the ROMs are also dense but of substantially lower dimension).

032[tc][tc]Modes (per time step)
034[tc][tc]Modes (per time step)
035[bc][bc]Excluded $\lambda$ fraction

037[bc][bc]Gaussian forcing
042[bc][bc]White forcing

009[ct][ct]$\matlabfragNegXTick 20$
010[ct][ct]$\matlabfragNegXTick 10$
020[ct][ct]$\matlabfragNegXTick 20$
021[ct][ct]$\matlabfragNegXTick 10$

Top panel: retained and unretained modes for $r = 10$ modes. The number of modes retained is frequency dependent, as shown in green on the right axis. The retained modes (red) are the overall highest-energy modes and the threshold (dashed) is determined as the energy of the $N_\omega r = 10240$-th most energetic mode. Bottom panel: the fraction of excluded energy as a function of the number of modes retained. This quantity vanishes faster in the case of Gaussian forcing, indicating that the trajectories are more accurately represented with a given number of modes in this case relative to the white-forcing case.

000[ct][ct]$\matlabfragNegXTick 20$
002[ct][ct]$\matlabfragNegXTick 20$
004[ct][ct]$\matlabfragNegXTick 20$
006[ct][ct]$\matlabfragNegXTick 20$
013[ct][ct]$\matlabfragNegXTick 20$
015[ct][ct]$\matlabfragNegXTick 20$
017[ct][ct]$\matlabfragNegXTick 50$

A FOM trajectory with the white forcing and the $2$-mode POD-Galerkin, balanced truncation, and SPOD Petrov-Galerkin approximations thereof, along with the errors in these ROM approximations. The error fields shown are the absolute value of the difference of the FOM and ROM trajectories. The peak error value (and the upper limit on the error color scale) is $87 \%$ of the peak absolute value of the state.

Figure <ref> illustrates various features of the mode energies. The top panel shows the SPOD mode energies $\lambda$ as a function of $\omega$. Each curve is a particular mode number as a function of frequency. The decision to retain $N_\omega r$ total modes in the ROM selects an energy threshold below which modes are not retained. After ordering the energies of all mode numbers at all frequencies, the threshold is given by the $N_\omega r$-th energy $\tilde{\lambda}_{N_\omega r}$. At frequencies where a given mode number is above (red) the energy threshold (dashed), it is retained. Where it is below (blue) the energy threshold, it is truncated. The green curve shows the number of modes that meet or exceed this threshold as a function of $\omega$. For example, at the dominant frequency in the white noise case, $22$ modes are retained, whereas at the highest and lowest frequencies, only $3$ are retained. The highlighted mode numbers are the lowest mode that is always retained, and the highest that is always excluded. The bottom panel shows the fraction of energy that is excluded depending on the number of modes retained. For the Gaussian forcing case, this quantity drops more steeply initially, indicating that the state is more accurately represented with a given number of SPOD modes in the Gaussian case relative to the white case.

The test data comprises $173$ trajectories, again with $\Delta t = 0.2$ and with each trajectory of length $1024$ time steps. For each ROM, we calculate the error at the $1024$ points for each trajectory as the square norm of the difference with the FOM solution. Throughout this section, we compare the performance of the proposed SPOD Petrov-Galerkin method to that of POD-Galerkin and an enhanced version of balanced truncation. In its usual form, balanced truncation does not make use of (i.e., does not require) data or knowledge of the statistics of the problem it is applied to. It may be improved with this information by `whitening' the forcing, i.e., transforming the system to one where the forcing is spatially white before performing the usual balanced truncation algorithm. In this application, where the forcing is far from spatially white, we observe that this variant of balanced truncation substantially outperforms the standard version, so we use it as a benchmark along with POD-Galerkin. We solve the reduced equations with  for both methods.

Figure <ref> shows the three ROM approximations of a trajectory, along with the errors in these approximations. All ROMs here use $r = 2$ modes, and the error field shown here is the absolute value of the difference between the FOM and ROM trajectories. Balanced truncation produces a better result than POD-Galerkin, but the error of the proposed SPOD Petrov-Galerkin method is much lower than both benchmarks.

Figure <ref> shows the error for the three ROMs as a function of time. The error reported is the square norm of the difference with the FOM averaged over all $173$ trajectories and normalized by the mean square norm of the state. For POD-Galerkin and balanced truncation, $10$ spatial modes are used for each of the $1024$ time steps. For the SPOD Petrov-Galerkin method, the $10 \times 1024 = 10240$ most energetic SPOD modes are used, and are distributed over the frequencies as shown in Figure <ref>. Most notably, the error is nearly two orders of magnitude smaller for the SPOD Petrov-Galerkin method than it is for the other two methods. Given the analysis in the previous sections, this is not surprising: the SPOD modes are (nearly) optimal in that the representation error with some number of SPOD modes is smaller than (nearly) every other space-time basis. Again, this representation is recovered by the Petrov-Galerkin method up to the errors introduced from approximating the operators and the non-periodicity of the forcing on the temporal interval. The error from all methods is larger in the white-noise forcing case. This is to be expected because the resulting behavior of the state is higher rank in this case relative to the Gaussian forcing. The error of the SPOD Petrov-Galerkin method decreases by more in the Gaussian forcing case because it takes explicit advantage of the additional spatiotemporal coherence relative to the white forcing case.

015[bc][bc]Gaussian forcing
018[bc][bc]White forcing

Error with $r = 10$ modes relative to the FOM, averaged over $173$ trajectories. The large difference in accuracy is due, in large part, to the ability of the SPOD modes to represent trajectories more accurately than space-only modes. This difference is larger in the Gaussian forcing case.

Next, we investigate the dependence of the error on the number of modes retained. Figure <ref> shows, as one might expect, that the error in all methods decreases with the number of modes retained. The gulf between the SPOD Petrov-Galerkin method and the two space-only ROMs is roughly maintained over the range of modes shown for both the white and Gaussian forcing cases. As SPOD-PG approaches $20$ modes, the error shallows because it approaches the limit of the full-order frequency domain error of roughly $10^{-6}$, which can be seen in Figure <ref>. The dashed lines are the projection of the exact solution onto each respective set of modes, which we refer to as the representation error. For example, the dashed green line is the SPOD mode representation error, i.e., the error of the FOM solution projected onto the span of the SPOD modes. We emphasize that the motivation for this work is the fact that the SPOD representation error is substantially below the POD representation error. The SPOD-PG solution error and representation error are nearly identical before the accuracy of the former is limited by the full-order frequency domain error at around $16$ modes and $10^{-6}$ error. Until this point, the SPOD-PG method indeed achieves the lowest error possible using SPOD modes, which is nearly the lowest error with any set of space-time modes. The SPOD-PG solution error is not only lower than the POD-Galerkin and balanced truncation solution errors, but also the respective representation errors of these bases. The POD representation error is a lower bound for the error for any time-domain Petrov-Galerkin method, such as balanced truncation or least-squares Petrov-Galerkin, because this is precisely the quantity that POD modes minimize. Indeed, the balanced truncation error is within this bound. We view the fact that the SPOD-PG solution error is significantly below the POD representation error as one of the major achievements of this work.

016[tc][tc]$r$ (# modes)
019[tc][tc]$r$ (# modes)

015[bc][bc]Gaussian forcing
018[bc][bc]White forcing

Average error as a function of number of space-time modes for the proposed SPOD Petrov-Galerkin method, POD-Galerkin and balanced truncation. The $x$-axis is the total number of space-time modes divided by the number of times steps. For the latter two methods, this is the number of spatial modes.

017[tc][tc]$r$ (# modes)
020[tc][tc]$r$ (# modes)
021[bc][bc]CPU Time (fraction of FOM)

016[bc][bc]Gaussian forcing
019[bc][bc]White forcing

Average CPU time as a function of the number of space-time modes used as a fraction of the FOM CPU time.

In Figure <ref>, we show the CPU time as a function of the number of modes retained. All values are reported as a fraction of the FOM time of $3543$ seconds in total for the $173$ runs. The SPOD Petrov-Galerkin CPU time scales linearly with the number of modes retained, but here, there are too few modes to see this scaling. Nonetheless, the SPOD Petrov-Galerkin method is substantially faster than the two benchmarks and runs in roughly two thousandths of the FOM time. This time includes the time to take the Fourier transform of the forcing and the inverse Fourier transform of the response. The times for the white forcing cases and Gaussian forcings are comparable for the SPOD Petrov-Galerkin method, but the FOM takes substantially longer with the white forcing.

013[bc][bc]$\nabla \times  {\bf u}$

The vorticity of the velocity field $\ket{u}(\ket{x})$ that transports the scalar, and a snapshot of the scalar. The region that is forced is shown in the red circle.

016[bc][bc]SPOD-PG $r=20$
019[bc][bc]POD-G $r=20$
023[bc][bc]SPOD-PG $r=10$
025[bc][bc]POD-G $r=10$
028[bc][bc]SPOD-PG $r=2$
030[bc][bc]POD-G $r=2$

The error fields for snapshots of the POD-Galerkin and SPOD Petrov-Galerkin ROMs for $2$,$10$, and $20$ modes. Note that there is less error in the $2$-mode SPOD Petrov-Galerkin snapshot than the $20$-mode POD-Galerkin snapshot. The extremes of the color scale are $31 \%$ of the maximum absolute value of the FOM solution, and the maximum in the $2$-mode POD-Galerkin error field is $62\%$ of the maximum absolute value of the FOM solution.

§.§ Scalar transport problem

Next, we demonstrate the proposed algorithm on an advection-diffusion system modeling the transport of a scalar quantity in a steady fluid flow. The flow profile is the mean of a lid-driven cavity flow simulation at $Re = 30,000$. This problem differs from the Ginzburg-Landau example in three important ways: it is substantially larger – $N_x = 9604$ as opposed to $220$ in the Ginzburg-Landau case – the matrix $\td{A}$ is sparse, and the forcing occupies only a subset of the domain. The former means that the model is too large to compute the matrix operations without the approximations we described earlier. With this large $N_x$, computing the Gramians in balanced truncation is too costly as well, so we do not compare to it here, though we do note that there are effective data-driven approximations of it as well [40, 31]. Balanced truncation must have greater error than the POD representation error, which we do show, and can be expected to share the CPU time of POD-Galerkin. That the forcing does not occupy the entire domain means that there is a CPU time savings in using the intermediary basis.

The continuous governing equations for the scalar transport case may be written in the same form as (<ref>), where $\mathcal{A}$ is now defined as

\begin{equation}
\mathcal{A} = -\ket{u}(\td{x}) \cdot \nabla  + \eta \nabla^2  \text{,}
\end{equation}

and where $\ket{u}(\td{x})$ is the mean flow in the lid-driven cavity. We take $\eta = 0.001$ and the velocity of the lid to be Gaussian in $x$ (as opposed to a constant) to avoid discontinuities at the upper corners. The problem is nondimensionalized such that the maximum speed, occurring in the center of the lid, is $1$. We prescribe a forcing that is stochastic with Gaussian spatial and temporal autocorrelation in a region of Gaussian support centered at $\overline{\ket{x}} = [0.75,0.25]^T$; its statistics are given by

\begin{equation} \label{eq:ST forcing stats}
C_{ff}(\ket{x}_1 \ket{x}_2,t_1,t_2) = \exp \Bigg[- \left( \frac{|\ket{x}_1 - \overline{\ket{x}}|^2 + |\ket{x}_2 - \overline{\ket{x}}|^2}{l^2} + \frac{|\ket{x}_2 - \ket{x}_1|^2}{\xi^2} + \frac{(t_2 - t_1)^2}{\tau^2} \right) \Bigg] \text{,}
\end{equation}

where $| \ \cdot \ |$ is Euclidean distance. Here, $l = 0.1$ is the spatial width of the support of the forcing, $\xi = 0.07$ is the spatial correlation length, and $\tau = 1$ the temporal correlation length. We use a second-order finite difference discretization with $98$ points in both directions and Dirichlet boundary conditions so as to mimic a heat bath. The FOM is solved using MATLAB's . The vorticity of the underlying velocity field and a snapshot of the transported scalar are shown in Figure <ref>. The red dashed circle in the latter indicates the region in which the equations are substantially forced – it is the radius at which the autocorrelation in the forcing drops by a factor of $e$ from its peak value (see (<ref>)).

For this example, we gather $50,000$ time steps of FOM data spaced $\Delta t = 0.5$ apart from which to calculate SPOD modes and approximate operators. Time is nondimensionalized such that it takes one time unit for the lid to cross the cavity. In this example, the time step used in the FOM integration was much smaller – more than an order of magnitude for most times – due to the stiffness of the system. The data was then segmented into $648$ overlapping blocks each $256$ time steps in length, so $N_\omega = 256$, and $T = 128$. The test data is $128$ trajectories of the system, and error is defined in the same way as before – as the square norm of the difference of the ROM and FOM solutions averaged over all test trajectories and normalized by the mean square norm of the solution itself. Finally, we use $p = 50$ for the intermediary basis.

Figure <ref> shows error field snapshots for the two ROMs for $2$, $10$, and $20$ modes. The error field for POD-Galerkin with $20$ modes is larger than that for the proposed method with $2$ modes.

Figure <ref> shows the accuracy over a range of mode numbers. Once again, the proposed method produces error two orders of magnitude lower than does POD-Galerkin. The dashed lines are again the error just due to projecting the solution into the space spanned by the respective modes. The error of the POD projection is a lower bound for the POD-Galerkin error, as well as any spatial Petrov-Galerkin method, such as balanced truncation. We see that in this problem, like in the Ginzburg-Landau problem, the SPOD-method far undershoots even this bound, and is indeed fairly close to its own error bound.

017[tc][tc]$r$ (# modes)

Accuracy of the proposed method compared to that of POD-Galerkin applied to the scalar transport problem. The dashed lines are the error of the full-order solution projected on the respective bases and are lower bounds for the error of the methods.

Figure <ref> shows the CPU time required to solve the scalar transport problem. The SPOD Petrov-Galerkin method is slightly faster in this case than POD-Galerkin. Too few modes are used to see the asymptotic linear scaling with the number of modes. Timing results from using DEIM to remove the $N_f$ scaling are also shown for both methods, and the DEIM-augmented version of the method is described in Appendix <ref>. In the appendix, we also show that there is no noticeable decrease in error, and discuss some drawbacks of the augmented version of the method.

015[tc][tc]$r$ (# modes)
016[bc][bc]CPU time (fraction of FOM)

Average CPU time to solve the scalar transport problem as a function of the number of space-time modes used as a fraction of the FOM CPU time. The timing results from using DEIM to remove the $N_f$ scaling from the method (and from POD-Galerkin) are also shown, and the details for this method are given in Appendix <ref>.

§ CONCLUSIONS

Space-time bases allow for a more accurate representation of a trajectory than do space-only bases with the same number of coefficients. In particular, the SPOD encoding of a trajectory with some number of coefficients may be orders of magnitude more accurate than the POD encoding of the same trajectory with the same number of coefficients. The obvious objectives are, therefore, to solve for these coefficients quickly and accurately, and we have pursued these for linear time-invariant systems.

We show, via two examples, that the method we derive can indeed accomplish both these objectives: the SPOD method takes comparable CPU time to benchmark time-domain methods like POD-Galerkin and balanced truncation, and the solution is roughly two orders of magnitude more accurate in the SPOD case than both benchmarks. In fact, the SPOD solution is nearly two orders of magnitude more accurate than the projection of the FOM solution onto the POD modes, which is the lower bound on error for any time-domain Petrov-Galerkin method.

A few negative aspects of the method are worth mentioning. The most limiting is that the method requires the entire forcing over the time interval of interest to be known before beginning the computation; the first step of the method is to take a FFT of the forcing, which cannot be done without the entire forcing in time. For some applications, this prevents the method from being applicable, while for others, it is not a problem. Second, the method works on a preprescribed interval $[0,T]$. If one wishes to obtain the solution longer than this interval, one can repeat the method with the value at the end of the interval as the initial condition. This is more cumbersome than extending the solution in a time-domain method. Finally, SPOD modes require more training data to obtain than POD modes, which limits the applicability of the proposed method in cases where training data is scarce. Where these disadvantages are not obstacles, however, we have shown the proposed method to be substantially more accurate at the same CPU time compared to standard methods for linear model reduction. We hope that this will aid in applications of linear model reduction and increase interest in space-time methods.

§ ACKNOWLEDGMENTS

P.F. and A.T. gratefully acknowledge funding from the National Science Foundation grant No. 2237537.

§ FOURIER REPRESENTATION OF THE TIME DERIVATIVE

Here, we illustrate why substituting $i\omega_k \hat{\ket{q}}_k$ for $\dot{\ket{q}}$ when going to the frequency domain is only correct if $\ket{q}$ starts and ends at the same value on the interval. It is easier to do this using the continuous definition of the Fourier coefficients, rather than the DFT, so we compute the integral

\begin{equation}
\hat{\dot{\ket{q}}}_k = \fint \dot{\ket q} e^{-i \omega_k t} \ dt\text{.}
\end{equation}

This integral may be solved by parts,

\begin{equation}
\hat{\dot{\ket{q}}}_k = \frac{1}{T} \ket{q} \big|_{0}^{T} - \fint -i \omega_k \ket{q}(t) e^{-i \omega_k t} dt  \text{.}
\end{equation}

The boundary term gives $\frac{\Delta \ket{q}}{T}$ [21], where $\Delta \ket{q} = \ket{q}(T) - \ket{q}(0)$, and the integral is the naive term $i\omega_k \hat{\ket q}_k$. The Fourier representation of the derivative is thus

\begin{equation}
\hat{\dot{\ket{q}}}_k = i\omega_k \hat{\ket q}_k + \frac{\Delta \ket{q}}{T} \text{.}
\end{equation}

§ DEIM-AUGMENTED ALGORITHM

Here, we present a means of decreasing the scaling in cases where the forcing is spatially structured, but $N_f$ is large. The idea is to use a sparse sampling of the forcing vectors, which are size $N_f$, using the discrete empirical interpolation method (DEIM) [5]. We also sparse sample the two vectors that are size $N_x$, the initial condition and the forcing sum terms in (<ref>), though these terms can also be handled using an intermediary basis, as before. The rank-$p$ DEIM approximation of a vector $\ket{v} \in \mathbb{C}^{N_x}$ is $\ket{v} \approx \td{U}_{\ket{v}} (\td{P}^T_{\ket{v}}\td{U}_{\ket{v}})^{-1}\td{P}^T_{\ket{v}} \ket{v}$, where the columns of $\td{U}_{\ket{v}} \in \mathbb{C}^{N_x \times p}$ are the POD modes for the ensemble from which  $\ket{v}$ is a sample, and $\td{P}^T_{\ket{v}} \in \{0,1\}^{p \times N_x}$ samples $p$ elements from $\ket{v}$ and is formed via the DEIM algorithm.

The DEIM algorithm is run for the forcing in the time domain, giving a set of sample points $\td{P}^T_{\ket{f}} \in \mathbb{C}^{p \times N_f}$ from which the forcing can be reconstructed accurately. The structures in the forcing at different frequencies will, in general, be different, so it is best to use a different spatial basis at each frequency to complete the DEIM approximation. We label the spatial basis for the forcing at the $k$-th frequency $\td{U}_{\hat{\ket{f}}_k}$. The approximation for the $k$-th forcing is $\hat{\ket{f}}_k \approx \td{U}_{\hat{\ket{f}}_k}(\td{P}^T_{\ket{f}}\td{U}_{\hat{\ket{f}}_k})^{-1} \td{P}^T_{\ket{f}} \hat{\ket{f}}_k$, so the first term in (<ref>) is now

\begin{equation}
\td{E}_k \hat{\ket{f}}_k \approx \td{K}_{\hat{\ket{f}}_k}  \td{P}^T_{\ket{f}} \hat{\ket{f}}_k \text{,}
\end{equation}

where the matrix $\td{K}_{\hat{\ket{f}}_k} = \td{E}_k \td{U}_{\hat{\ket{f}}_k}(\td{P}^T_{\ket{f}}\td{U}_{\hat{\ket{f}}_k})^{-1} \in \mathbb{C}^{r_k \times p}$ is precomputed for each frequency. Note that the same sampling is used for each frequency; if the sampling were different for each frequency, one would need to take the DFT of the entire forcing (or the union of all the samplings), which would negate the scaling benefit of sparse sampling. This approach for the forcing term can be used in conjunction with the intermediary basis approach for the initial condition and forcing sum terms in (<ref>). These latter terms can also be handled with DEIM, which we show now.

The initial condition and forcing sum terms can be approximated with DEIM matrices $\td{U}_{\ket{q}_0}$ and $\td{P}^T_{\ket{q}_0}$ for the initial condition, and  $\td{U}_{fs}$ and $\td{P}^T_{fs}$ for the forcing sum. For the former, these matrices are formed by gathering all initial conditions within the training data and running the DEIM algorithm to obtain $\td{U}_{\ket{q}_0}$ and $\td{P}^T_{\ket{q}_0}$. From these matrices, the initial condition multiplied by $\td{F}_k$ is approximated as

\begin{equation}
\td{F}_k \ket{q}_0 \approx \td{K}_{\ket{q}_0} \td{P}^T_{\ket{q}_0} \ket{q}_0 \text{,}
\end{equation}

where  $\td{K}_{\ket{q}_0,k} = \td{F}_k \td{U}_{\ket{q}_0} (\td{P}^T_{\ket{q}_0} \td{U}_{\ket{q}_0})^{-1} \in \mathbb{C}^{r_k \times p}$ is precomputed. Similarly, $\td{U}_{fs}$ and $\td{P}^T_{fs}$ are obtained by running DEIM on the set of forcing sums, which must be calculated from each trajectory in the training data. With these matrices, the forcing sum multiplied by $\td{F}_k$ is approximated as

\begin{equation}
\td{F}_k \frac{1}{N_\omega} \sum_{l} \td{\Psi}_l^r \td{E}_l \hat{\ket{f}}_l \approx \frac{1}{N_t} \td{K}_{fs,l} \sum_{l}  \td{T}_{l,fs} \td{E}_l \hat{\ket{f}}_l \text{,}
\end{equation}

where $\td{K}_{fs,l} = \td{F}_k \td{U}_{fs} (\td{P}^T_{fs} \td{U}_{fs})^{-1} \in \mathbb{C}^{r_k \times p}$ and $\td{T}_{l,fs} = \td{P}^T_{fs}\td{\Psi}_l \in \mathbb{C}^{p \times r_l}$ are both precomputed. These operators are kept separate (as opposed to multiplied as a precomputation step) to avoid $N_\omega^2$ scaling. With these approximations, the DEIM-augmented equation is

\begin{equation}
\ket{a}_k =  \td{K}_{\hat{\ket{f}}_k}  \td{P}^T_{\ket{f}} \hat{\ket{f}}_k + \td{K}_{\ket{q}_0} \td{P}^T_{\ket{q}_0} \ket{q}_0 - \frac{1}{N_t} \td{K}_{fs,l} \sum_{l}  \td{T}_{l,fs} \td{E}_l \hat{\ket{f}}_l \text{.}
\end{equation}

This method removes the $N_f$ scaling and replaces it with $p$, thus, the DEIM-augmented version of the algorithm scales like $\mathcal{O}( (rp + p\log N_\omega)N_\omega )$.

017[tc][tc]$r$ (# modes)

Accuracy of the DEIM-augmented version of the method compared to POD-Galerkin (and the DEIM version thereof) applied to the scalar transport problem. As long as enough sample points are used, DEIM does not introduce additional error. Again, the dashed lines are the error of the full-order solution projected on the respective bases.

To demonstrate the augmented version of the method, we apply it to the scalar transport problem described in the main text. The forcing in this problem meets the conditions for a large improvement: $N_f = 2050$ is large, but the forcing is spatially structured, coming from (<ref>), so it can be approximated via sparse sampling effectively. Figure <ref> shows the error of the DEIM-augmented version with $p = 200$ along with that of the non-augmented version. There is almost no error sacrifice relative to the non-approximated method with this number of sample points. The timing for the method applied to the scalar transport problem is shown in Figure <ref>, along with the timing for a DEIM-augmented version of POD-Galerkin. Indeed, DEIM offers substantial speedup, both for the proposed method and for the POD-Galerkin method.

Two drawbacks of the DEIM-augmented version lead us to favor the non-augmented method in most cases. First, the DEIM-augmented version of the method relies on the initial condition and forcing being accurately approximated by sparse samplings, and these sparse samplings will only be accurate if the initial condition and forcing are similar in character to those in the training data. Second, the DEIM-augmented version is cumbersome to implement, requiring more precomputation of modes and matrices, relative to the non-augmented method.

[1]
N. Aubry, P. Holmes, J. L. Lumley, and E. Stone.
The dynamics of coherent structures in the wall region of a turbulent boundary layer.
Journal of Fluid Mechanics, 192:115–173, July 1988.

[2]
S. Bagheri, D. S. Henningson, J. Hœpffner, and P. J. Schmid.
Input-output analysis and control design applied to a linear model of spatially developing flows.
Applied Mechanics Reviews, 62(2), Feb. 2009.

[3]
K. Carlberg, C. Bou‐Mosleh, and C. Farhat.
Efficient non‐linear model reduction via a least‐squares Petrov–Galerkin projection and compressive tensor approximations.
International Journal for Numerical Methods in Engineering, 86(2):155–181, Oct. 2010.

[4]
K. Carlberg, C. Farhat, J. Cortial, and D. Amsallem.
The GNAT method for nonlinear model reduction: Effective implementation and application to computational fluid dynamics and turbulent flows.
Journal of Computational Physics, 242:623–647, June 2013.

[5]
S. Chaturantabut and D. C. Sorensen.
Nonlinear model reduction via discrete empirical interpolation.
SIAM Journal on Scientific Computing, 32(5):2737–2764, Jan. 2010.

[6]
K. K. Chen and C. W. Rowley.
Optimal actuator and sensor placement in the linearised complex Ginzburg–Landau system.
Journal of Fluid Mechanics, 681:241–260, June 2011.

[7]
Y. Choi, P. Brown, W. Arrighi, R. Anderson, and K. Huynh.
Space–time reduced order model for large-scale linear dynamical systems with application to Boltzmann transport problems.
Journal of Computational Physics, 424:109845, 2021.

[8]
Y. Choi and K. Carlberg.
Space–time least-squares petrov–galerkin projection for nonlinear model reduction.
SIAM Journal on Scientific Computing, 41(1):A26–A58, 2019.

[9]
K. Ekici and K. C. Hall.
Nonlinear analysis of unsteady flows in multistage turbomachines using harmonic balance.
AIAA Journal, 45(5):1047–1057, May 2007.

[10]
A. Farghadan, E. Martini, and A. Towne.
Scalable resolvent analysis for three-dimensional flows, 2023.

[11]
P. Frame and A. Towne.
Space-time POD and the Hankel matrix.
PLOS ONE, 18(8):e0289637, Aug. 2023.

[12]
F. J. Gonzalez and M. Balajewicz.
Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems.
arXiv preprint arXiv:1808.01346, 2018.

[13]
M. D. Grigoriu.
Linear Dynamical Systems.
Springer International Publishing, 2021.

[14]
K. C. Hall, K. Ekici, J. P. Thomas, and E. H. Dowell.
Harmonic balance methods applied to computational fluid dynamics problems.
International Journal of Computational Fluid Dynamics, 27(2):52–67, Jan. 2013.

[15]
K. C. Hall, J. P. Thomas, and W. S. Clark.
Computation of unsteady nonlinear flows in cascades using a harmonic balance technique.
AIAA Journal, 40(5):879–886, May 2002.

[16]
Y. Kim, K. Wang, and Y. Choi.
Efficient space–time reduced order model for linear dynamical systems in python using less than 120 lines of code.
Mathematics, 9(14), 2021.

[17]
K. Lee and K. T. Carlberg.
Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.
Journal of Computational Physics, 404:108973, Mar. 2020.

[18]
C. Lin.
Model order reduction in the frequency domain via spectral proper orthogonal decomposition.
Master's thesis, University of Illinois at Urbana-Champaign, 2019.

[19]
J. L. Lumley.
The structure of inhomogeneous turbulent flows.
Atmospheric Turbulence and Radio Wave Propagation, 1967.

[20]
E. Martini, A. V. G. Cavalieri, P. Jordan, and L. Lesshafft.
Accurate frequency domain identification of ODEs with arbitrary signals, 2019.

[21]
E. Martini, D. Rodríguez, A. Towne, and A. V. Cavalieri.
Efficient computation of global resolvent modes.
Journal of Fluid Mechanics, 919, May 2021.

[22]
B. J. McKeon.
The engine behind (wall) turbulence: perspectives on scale interactions.
Journal of Fluid Mechanics, 817:P1, 2017.

[23]
B. J. McKeon and A. S. Sharma.
A critical-layer framework for turbulent pipe flow.
Journal of Fluid Mechanics, 658:336–382, July 2010.

[24]
B. Moore.
Principal component analysis in linear systems: Controllability, observability, and model reduction.
IEEE Transactions on Automatic Control, 26(1):17–32, Feb. 1981.

[25]
B. R. Noack, K. Afanasiev, M. Morzyński, G. Tadmor, and F. Thiele.
A hierarchy of low-dimensional models for the transient and post-transient cylinder wake.
Journal of Fluid Mechanics, 497:335–363, Dec. 2003.

[26]
P. A. S. Nogueira, P. Morra, E. Martini, A. V. G. Cavalieri, and D. S. Henningson.
Forcing statistics in resolvent analysis: application in minimal turbulent Couette flow.
Journal of Fluid Mechanics, 908, Dec. 2020.

[27]
S. E. Otto, A. Padovan, and C. W. Rowley.
Optimizing oblique projections for nonlinear systems using trajectories.
SIAM Journal on Scientific Computing, 44(3):A1681–A1702, June 2022.

[28]
A. Padovan, B. Vollmer, and D. J. Bodony.
Data-driven model reduction via non-intrusive optimization of projection operators and reduced-order dynamics, 2024.

[29]
E. J. Parish and K. T. Carlberg.
Windowed least-squares model reduction for dynamical systems.
Journal of Computational Physics, 426:109939, 2021.

[30]
B. Peherstorfer and K. Willcox.
Data-driven operator inference for nonintrusive projection-based model reduction.
Computer Methods in Applied Mechanics and Engineering, 306:196–215, July 2016.

[31]
C. W. Rowley.
Model reduction for fluids, using balanced proper orthogonal decomposition.
International Journal of Bifurcation and Chaos, 15(03):997–1013, Mar. 2005.

[32]
C. W. Rowley, T. Colonius, and R. M. Murray.
Model reduction for compressible flows using POD and Galerkin projection.
Physica D: Nonlinear Phenomena, 189(1–2):115–129, Feb. 2004.

[33]
O. T. Schmidt and P. J. Schmid.
A conditional space–time POD formalism for intermittent and rare events: example of acoustic bursts in turbulent jets.
Journal of Fluid Mechanics, 867, Apr. 2019.

[34]
O. T. Schmidt and A. Towne.
An efficient streaming algorithm for spectral proper orthogonal decomposition.
Computer Physics Communications, 237:98–109, Apr. 2019.

[35]
O. T. Schmidt, A. Towne, G. Rigas, T. Colonius, and G. A. Brès.
Spectral analysis of jet turbulence.
Journal of Fluid Mechanics, 855:953–982, Sept. 2018.

[36]
L. Sirovich.
Turbulence and the dynamics of coherent structures. ii. symmetries and transformations.
Quarterly of Applied Mathematics, 45(3):573–582, 1987.

[37]
A. Towne.
Space-time Galerkin projection via spectral proper orthogonal decomposition and resolvent modes.
In AIAA Scitech 2021 Forum. American Institute of Aeronautics and Astronautics, Jan. 2021.

[38]
A. Towne, O. T. Schmidt, and T. Colonius.
Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis.
Journal of Fluid Mechanics, 847:821–867, 2018.

[39]
L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll.
Hydrodynamic stability without eigenvalues.
Science, 261(5121):578–584, 1993.

[40]
K. Willcox and J. Peraire.
Balanced model reduction via the proper orthogonal decomposition.
AIAA Journal, 40(11):2323–2330, Nov. 2002.

[41]
M. Yano, A. T. Patera, and K. Urban.
A space-time hp-interpolation-based certified reduced basis method for Burgers' equation.
Mathematical Models and Methods in Applied Sciences, 24(09):1903–1935, May 2014.

[1]
N. Aubry, P. Holmes, J. L. Lumley, and E. Stone.
The dynamics of coherent structures in the wall region of a turbulent boundary layer.
Journal of Fluid Mechanics, 192:115–173, July 1988.

[2]
S. Bagheri, D. S. Henningson, J. Hœpffner, and P. J. Schmid.
Input-output analysis and control design applied to a linear model of spatially developing flows.
Applied Mechanics Reviews, 62(2), Feb. 2009.

[3]
K. Carlberg, C. Bou‐Mosleh, and C. Farhat.
Efficient non‐linear model reduction via a least‐squares Petrov–Galerkin projection and compressive tensor approximations.
International Journal for Numerical Methods in Engineering, 86(2):155–181, Oct. 2010.

[4]
K. Carlberg, C. Farhat, J. Cortial, and D. Amsallem.
The GNAT method for nonlinear model reduction: Effective implementation and application to computational fluid dynamics and turbulent flows.
Journal of Computational Physics, 242:623–647, June 2013.

[5]
S. Chaturantabut and D. C. Sorensen.
Nonlinear model reduction via discrete empirical interpolation.
SIAM Journal on Scientific Computing, 32(5):2737–2764, Jan. 2010.

[6]
K. K. Chen and C. W. Rowley.
Optimal actuator and sensor placement in the linearised complex Ginzburg–Landau system.
Journal of Fluid Mechanics, 681:241–260, June 2011.

[7]
Y. Choi, P. Brown, W. Arrighi, R. Anderson, and K. Huynh.
Space–time reduced order model for large-scale linear dynamical systems with application to Boltzmann transport problems.
Journal of Computational Physics, 424:109845, 2021.

[8]
Y. Choi and K. Carlberg.
Space–time least-squares petrov–galerkin projection for nonlinear model reduction.
SIAM Journal on Scientific Computing, 41(1):A26–A58, 2019.

[9]
K. Ekici and K. C. Hall.
Nonlinear analysis of unsteady flows in multistage turbomachines using harmonic balance.
AIAA Journal, 45(5):1047–1057, May 2007.

[10]
A. Farghadan, E. Martini, and A. Towne.
Scalable resolvent analysis for three-dimensional flows, 2023.

[11]
P. Frame and A. Towne.
Space-time POD and the Hankel matrix.
PLOS ONE, 18(8):e0289637, Aug. 2023.

[12]
F. J. Gonzalez and M. Balajewicz.
Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems.
arXiv preprint arXiv:1808.01346, 2018.

[13]
M. D. Grigoriu.
Linear Dynamical Systems.
Springer International Publishing, 2021.

[14]
K. C. Hall, K. Ekici, J. P. Thomas, and E. H. Dowell.
Harmonic balance methods applied to computational fluid dynamics problems.
International Journal of Computational Fluid Dynamics, 27(2):52–67, Jan. 2013.

[15]
K. C. Hall, J. P. Thomas, and W. S. Clark.
Computation of unsteady nonlinear flows in cascades using a harmonic balance technique.
AIAA Journal, 40(5):879–886, May 2002.

[16]
Y. Kim, K. Wang, and Y. Choi.
Efficient space–time reduced order model for linear dynamical systems in python using less than 120 lines of code.
Mathematics, 9(14), 2021.

[17]
K. Lee and K. T. Carlberg.
Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.
Journal of Computational Physics, 404:108973, Mar. 2020.

[18]
C. Lin.
Model order reduction in the frequency domain via spectral proper orthogonal decomposition.
Master's thesis, University of Illinois at Urbana-Champaign, 2019.

[19]
J. L. Lumley.
The structure of inhomogeneous turbulent flows.
Atmospheric Turbulence and Radio Wave Propagation, 1967.

[20]
E. Martini, A. V. G. Cavalieri, P. Jordan, and L. Lesshafft.
Accurate frequency domain identification of ODEs with arbitrary signals, 2019.

[21]
E. Martini, D. Rodríguez, A. Towne, and A. V. Cavalieri.
Efficient computation of global resolvent modes.
Journal of Fluid Mechanics, 919, May 2021.

[22]
B. J. McKeon.
The engine behind (wall) turbulence: perspectives on scale interactions.
Journal of Fluid Mechanics, 817:P1, 2017.

[23]
B. J. McKeon and A. S. Sharma.
A critical-layer framework for turbulent pipe flow.
Journal of Fluid Mechanics, 658:336–382, July 2010.

[24]
B. Moore.
Principal component analysis in linear systems: Controllability, observability, and model reduction.
IEEE Transactions on Automatic Control, 26(1):17–32, Feb. 1981.

[25]
B. R. Noack, K. Afanasiev, M. Morzyński, G. Tadmor, and F. Thiele.
A hierarchy of low-dimensional models for the transient and post-transient cylinder wake.
Journal of Fluid Mechanics, 497:335–363, Dec. 2003.

[26]
P. A. S. Nogueira, P. Morra, E. Martini, A. V. G. Cavalieri, and D. S. Henningson.
Forcing statistics in resolvent analysis: application in minimal turbulent Couette flow.
Journal of Fluid Mechanics, 908, Dec. 2020.

[27]
S. E. Otto, A. Padovan, and C. W. Rowley.
Optimizing oblique projections for nonlinear systems using trajectories.
SIAM Journal on Scientific Computing, 44(3):A1681–A1702, June 2022.

[28]
A. Padovan, B. Vollmer, and D. J. Bodony.
Data-driven model reduction via non-intrusive optimization of projection operators and reduced-order dynamics, 2024.

[29]
E. J. Parish and K. T. Carlberg.
Windowed least-squares model reduction for dynamical systems.
Journal of Computational Physics, 426:109939, 2021.

[30]
B. Peherstorfer and K. Willcox.
Data-driven operator inference for nonintrusive projection-based model reduction.
Computer Methods in Applied Mechanics and Engineering, 306:196–215, July 2016.

[31]
C. W. Rowley.
Model reduction for fluids, using balanced proper orthogonal decomposition.
International Journal of Bifurcation and Chaos, 15(03):997–1013, Mar. 2005.

[32]
C. W. Rowley, T. Colonius, and R. M. Murray.
Model reduction for compressible flows using POD and Galerkin projection.
Physica D: Nonlinear Phenomena, 189(1–2):115–129, Feb. 2004.

[33]
O. T. Schmidt and P. J. Schmid.
A conditional space–time POD formalism for intermittent and rare events: example of acoustic bursts in turbulent jets.
Journal of Fluid Mechanics, 867, Apr. 2019.

[34]
O. T. Schmidt and A. Towne.
An efficient streaming algorithm for spectral proper orthogonal decomposition.
Computer Physics Communications, 237:98–109, Apr. 2019.

[35]
O. T. Schmidt, A. Towne, G. Rigas, T. Colonius, and G. A. Brès.
Spectral analysis of jet turbulence.
Journal of Fluid Mechanics, 855:953–982, Sept. 2018.

[36]
L. Sirovich.
Turbulence and the dynamics of coherent structures. ii. symmetries and transformations.
Quarterly of Applied Mathematics, 45(3):573–582, 1987.

[37]
A. Towne.
Space-time Galerkin projection via spectral proper orthogonal decomposition and resolvent modes.
In AIAA Scitech 2021 Forum. American Institute of Aeronautics and Astronautics, Jan. 2021.

[38]
A. Towne, O. T. Schmidt, and T. Colonius.
Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis.
Journal of Fluid Mechanics, 847:821–867, 2018.

[39]
L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll.
Hydrodynamic stability without eigenvalues.
Science, 261(5121):578–584, 1993.

[40]
K. Willcox and J. Peraire.
Balanced model reduction via the proper orthogonal decomposition.
AIAA Journal, 40(11):2323–2330, Nov. 2002.

[41]
M. Yano, A. T. Patera, and K. Urban.
A space-time hp-interpolation-based certified reduced basis method for Burgers' equation.
Mathematical Models and Methods in Applied Sciences, 24(09):1903–1935, May 2014.

[1]
N. Aubry, P. Holmes, J. L. Lumley, and E. Stone.
The dynamics of coherent structures in the wall region of a turbulent boundary layer.
Journal of Fluid Mechanics, 192:115–173, July 1988.

[2]
S. Bagheri, D. S. Henningson, J. Hœpffner, and P. J. Schmid.
Input-output analysis and control design applied to a linear model of spatially developing flows.
Applied Mechanics Reviews, 62(2), Feb. 2009.

[3]
K. Carlberg, C. Bou‐Mosleh, and C. Farhat.
Efficient non‐linear model reduction via a least‐squares Petrov–Galerkin projection and compressive tensor approximations.
International Journal for Numerical Methods in Engineering, 86(2):155–181, Oct. 2010.

[4]
K. Carlberg, C. Farhat, J. Cortial, and D. Amsallem.
The GNAT method for nonlinear model reduction: Effective implementation and application to computational fluid dynamics and turbulent flows.
Journal of Computational Physics, 242:623–647, June 2013.

[5]
S. Chaturantabut and D. C. Sorensen.
Nonlinear model reduction via discrete empirical interpolation.
SIAM Journal on Scientific Computing, 32(5):2737–2764, Jan. 2010.

[6]
K. K. Chen and C. W. Rowley.
Optimal actuator and sensor placement in the linearised complex Ginzburg–Landau system.
Journal of Fluid Mechanics, 681:241–260, June 2011.

[7]
Y. Choi, P. Brown, W. Arrighi, R. Anderson, and K. Huynh.
Space–time reduced order model for large-scale linear dynamical systems with application to Boltzmann transport problems.
Journal of Computational Physics, 424:109845, 2021.

[8]
Y. Choi and K. Carlberg.
Space–time least-squares petrov–galerkin projection for nonlinear model reduction.
SIAM Journal on Scientific Computing, 41(1):A26–A58, 2019.

[9]
K. Ekici and K. C. Hall.
Nonlinear analysis of unsteady flows in multistage turbomachines using harmonic balance.
AIAA Journal, 45(5):1047–1057, May 2007.

[10]
A. Farghadan, E. Martini, and A. Towne.
Scalable resolvent analysis for three-dimensional flows, 2023.

[11]
P. Frame and A. Towne.
Space-time POD and the Hankel matrix.
PLOS ONE, 18(8):e0289637, Aug. 2023.

[12]
F. J. Gonzalez and M. Balajewicz.
Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems.
arXiv preprint arXiv:1808.01346, 2018.

[13]
M. D. Grigoriu.
Linear Dynamical Systems.
Springer International Publishing, 2021.

[14]
K. C. Hall, K. Ekici, J. P. Thomas, and E. H. Dowell.
Harmonic balance methods applied to computational fluid dynamics problems.
International Journal of Computational Fluid Dynamics, 27(2):52–67, Jan. 2013.

[15]
K. C. Hall, J. P. Thomas, and W. S. Clark.
Computation of unsteady nonlinear flows in cascades using a harmonic balance technique.
AIAA Journal, 40(5):879–886, May 2002.

[16]
Y. Kim, K. Wang, and Y. Choi.
Efficient space–time reduced order model for linear dynamical systems in python using less than 120 lines of code.
Mathematics, 9(14), 2021.

[17]
K. Lee and K. T. Carlberg.
Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.
Journal of Computational Physics, 404:108973, Mar. 2020.

[18]
C. Lin.
Model order reduction in the frequency domain via spectral proper orthogonal decomposition.
Master's thesis, University of Illinois at Urbana-Champaign, 2019.

[19]
J. L. Lumley.
The structure of inhomogeneous turbulent flows.
Atmospheric Turbulence and Radio Wave Propagation, 1967.

[20]
E. Martini, A. V. G. Cavalieri, P. Jordan, and L. Lesshafft.
Accurate frequency domain identification of ODEs with arbitrary signals, 2019.

[21]
E. Martini, D. Rodríguez, A. Towne, and A. V. Cavalieri.
Efficient computation of global resolvent modes.
Journal of Fluid Mechanics, 919, May 2021.

[22]
B. J. McKeon.
The engine behind (wall) turbulence: perspectives on scale interactions.
Journal of Fluid Mechanics, 817:P1, 2017.

[23]
B. J. McKeon and A. S. Sharma.
A critical-layer framework for turbulent pipe flow.
Journal of Fluid Mechanics, 658:336–382, July 2010.

[24]
B. Moore.
Principal component analysis in linear systems: Controllability, observability, and model reduction.
IEEE Transactions on Automatic Control, 26(1):17–32, Feb. 1981.

[25]
B. R. Noack, K. Afanasiev, M. Morzyński, G. Tadmor, and F. Thiele.
A hierarchy of low-dimensional models for the transient and post-transient cylinder wake.
Journal of Fluid Mechanics, 497:335–363, Dec. 2003.

[26]
P. A. S. Nogueira, P. Morra, E. Martini, A. V. G. Cavalieri, and D. S. Henningson.
Forcing statistics in resolvent analysis: application in minimal turbulent Couette flow.
Journal of Fluid Mechanics, 908, Dec. 2020.

[27]
S. E. Otto, A. Padovan, and C. W. Rowley.
Optimizing oblique projections for nonlinear systems using trajectories.
SIAM Journal on Scientific Computing, 44(3):A1681–A1702, June 2022.

[28]
A. Padovan, B. Vollmer, and D. J. Bodony.
Data-driven model reduction via non-intrusive optimization of projection operators and reduced-order dynamics, 2024.

[29]
E. J. Parish and K. T. Carlberg.
Windowed least-squares model reduction for dynamical systems.
Journal of Computational Physics, 426:109939, 2021.

[30]
B. Peherstorfer and K. Willcox.
Data-driven operator inference for nonintrusive projection-based model reduction.
Computer Methods in Applied Mechanics and Engineering, 306:196–215, July 2016.

[31]
C. W. Rowley.
Model reduction for fluids, using balanced proper orthogonal decomposition.
International Journal of Bifurcation and Chaos, 15(03):997–1013, Mar. 2005.

[32]
C. W. Rowley, T. Colonius, and R. M. Murray.
Model reduction for compressible flows using POD and Galerkin projection.
Physica D: Nonlinear Phenomena, 189(1–2):115–129, Feb. 2004.

[33]
O. T. Schmidt and P. J. Schmid.
A conditional space–time POD formalism for intermittent and rare events: example of acoustic bursts in turbulent jets.
Journal of Fluid Mechanics, 867, Apr. 2019.

[34]
O. T. Schmidt and A. Towne.
An efficient streaming algorithm for spectral proper orthogonal decomposition.
Computer Physics Communications, 237:98–109, Apr. 2019.

[35]
O. T. Schmidt, A. Towne, G. Rigas, T. Colonius, and G. A. Brès.
Spectral analysis of jet turbulence.
Journal of Fluid Mechanics, 855:953–982, Sept. 2018.

[36]
L. Sirovich.
Turbulence and the dynamics of coherent structures. ii. symmetries and transformations.
Quarterly of Applied Mathematics, 45(3):573–582, 1987.

[37]
A. Towne.
Space-time Galerkin projection via spectral proper orthogonal decomposition and resolvent modes.
In AIAA Scitech 2021 Forum. American Institute of Aeronautics and Astronautics, Jan. 2021.

[38]
A. Towne, O. T. Schmidt, and T. Colonius.
Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis.
Journal of Fluid Mechanics, 847:821–867, 2018.

[39]
L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll.
Hydrodynamic stability without eigenvalues.
Science, 261(5121):578–584, 1993.

[40]
K. Willcox and J. Peraire.
Balanced model reduction via the proper orthogonal decomposition.
AIAA Journal, 40(11):2323–2330, Nov. 2002.

[41]
M. Yano, A. T. Patera, and K. Urban.
A space-time hp-interpolation-based certified reduced basis method for Burgers' equation.
Mathematical Models and Methods in Applied Sciences, 24(09):1903–1935, May 2014.