A heuristic method for determining the minimum number of modes to be used for the reduced order model. The heuristic is based on the Shapiro-Wilkes hypothesis test comparing how well the empirical distribution is represented by a Gaussian distribution.
Koopman Reduced Order Modeling
Consider a stochastic dynamical system
\(\mathbf {x}(t+1) = T(\mathbf {x}(t), \xi (t)), \quad t\in  \mathbb {N}\) 
where \(\mathbf {x}(t) \in \Omega \subset  \mathbb {R}^d\) , \(\Omega \)  is a compact Riemannian manifold, and \(\xi (t)\)  is an IID noise process with \(\xi \sim \Xi \) . Let \(f_i \in L^2(M, d\mathbf {x})\)  for \(i=1,\dots , n\)  and \(\mathbf {f} = (f_1,\dots , f_n)\) . The Koopman operator induced by the stochastic dynamical system (and parameterized by the noise sequence) is denoted as \(U_\xi f(\mathbf {x}) = f(T(\mathbf {x}, \xi ))\) . With an abuse of notation, we will denote by \(U_\xi ^t f_i(\mathbf {x})\)  the composition of \(f_i\)  with the map (REF ) \(t\) -times with the that a new value for \(\xi (t)\)  is drawn from the distribution at each iteration; i.e.,
\(U_{\xi }^2 f_i(\mathbf {x}) = f_i ( T( \mathbf {x}(1), \xi (1)) ) = f_i ( T( T(\mathbf {x}, \xi (0)), \xi (1)) ).\) 
Furthermore, we define \(U_\xi \mathbf {f}(\mathbf {x}) = (U_\xi f_1(\mathbf {x}), \dots , U_\xi f_n(\mathbf {x}))\) . We assume there is a spectral decomposition for \(U_\xi \)  written as
\(U_{\xi }^t f_i(\mathbf {x}) = \sum _j c_j \lambda _j^t \phi _j(\mathbf {x}) + \int _z̏^t dE_c(z)f_i(\mathbf {x}),\) 
where \(c_j \in , \) j is an eigenvalue, \(\phi _j\)  is an eigenfunction, and \(dE_c\)  is a projection-valued measure corresponding to the continuous part of the spectrum. For the vector-valued observable case, we have
\(U_{\xi }^t \mathbf {f}(\mathbf {x}) = \sum _j \mathbf {m}_j \lambda _j^t \phi _j(\mathbf {x}) + \int _z̏^t dE_c(z) \mathbf {f}(\mathbf {x}),\) 
where \(\mathbf {m}_j\)  is the \(j\) -th Koopman mode.
Model reduction
Obviously, we cannot compute the infinite expansion, so we attempt to approximate the evolution of the observable due to the operator using a finite mode expansion. That is
\(U_\xi ^t \mathbf {f}(\mathbf {x}) = \sum _{j=1}^{J} \mathbf {m}_j \lambda _j^t \phi _j(\mathbf {x}) + \mathbf {r}(t),\) 
where here \(\mathbf {r}(t) \sim \mathfrak {R}\)  is a residual sequence. Note that if \(T\)  or any \(f_i\)  are nonlinear, then \(\mathbf {r}(t)\)  will not have the same distribution as \(\xi (t)\) . To make the notation more compact, let us write the deterministic part as
\(D_\xi ^t \mathbf {f}(\mathbf {x}) = \sum _{j=1}^{J} \mathbf {m}_j \lambda _j^t \phi _j(\mathbf {x}).\) 
So far we have split the evolution of the observable into a deterministic part and a stochastic part. Also note, that due to the truncation of the deterministic part from an infinite sum to a finite sum, the “coefficients”, \(\phi _j(x)\) , for the Koopman modes \(\mathbf {m}_j\)  may not give the optimal deterministic reconstruction of the evolution. We will come back to this when discussing algorithms to compute the decomposition.
Now we look into splitting the modal noise into a modal component and an innovation component. Let \(M = \operatorname{span}{\mathbf {m}_1, \dots , \mathbf {m}_J}\) , \(\mathbf {M}\)  the matrix having the Koopman modes as its columns, and \(P_M\)  the orthogonal projection onto \(M\) . The,
\(U_\xi ^t \mathbf {f}(\mathbf {x}) = \underbrace{\sum _{j=1}^{J} \mathbf {m}_j \lambda _j^t \phi _j(\mathbf {x}) + P_M \mathbf {r}(t)}_{\in M} + \underbrace{(I - P_M)\mathbf {r}(t)}_{\text{orthogonal to $M$}},\) 
Using the fact that \(P_M\)  can be computed as \(P_M = \mathbf {M} \mathbf {M}^+\) , where \(\mathbf {M}^+\)  is the Moore-Penrose pseudo-inverse of \(\mathbf {M}\) . Consider the terms lying in \(M\) :
\(\sum _{j=1}^{J} \mathbf {m}_j \lambda _j^t \phi _j(\mathbf {x}) + P_M \mathbf {r}(t)= \mathbf {M} \begin{bmatrix} \lambda _1^t \phi _1(\mathbf {x}) \\ \vdots \\ \lambda _J^t \phi _J(\mathbf {x})\end{bmatrix} + \mathbf {M} \mathbf {M}^+ \mathbf {r}(t)= \mathbf {M} \left(\begin{bmatrix} \lambda _1^t \phi _1(\mathbf {x}) \\ \vdots \\ \lambda _J^t \phi _J(\mathbf {x})\end{bmatrix} + \mathbf {M}^+ \mathbf {r}(t) \right).\) 

We define the modal noise \(\rho (t) = (\rho _1(t), \dots , \rho _J(t))^T\)  as
\(\rho (t) = \mathbf {M} \mathbf {M}^+ \mathbf {r}(t).\) 
The innovation noise is then defined as
\(\eta (t) = (I - P_M)\mathbf {r}(t) = \mathbf {r}(t) - P_M \mathbf {r}(t) = \mathbf {r}(t) - \rho (t).\) 
The innovation noise \(\eta (t)\)  is orthogonal to \(M\) . With these definitions (REF ) becomes
\(U_\xi ^t \mathbf {f}(\mathbf {x}) = \left( \mathbf {M}\begin{bmatrix} \lambda _1^t \phi _1(\mathbf {x}) \\ \vdots \\ \lambda _J^t \phi _J(\mathbf {x})\end{bmatrix} + \rho (t) \right) + \eta (t).\) 

Prediction and Forecasting
Given the stochastic reduced order model (REF ) one could compute predictions and statistics for those predictions if we knew the distributions of \(\rho (t)\)  and \(\eta (t)\) . Note that the distribution of these can be computed since these are just linear transformations of \(\mathbf {r}(t)\) . However, \(\mathfrak {R}\)  is usually unknown. Even if we had the analytic form of the distribution \(\Xi \)  that the noise process \(\xi \)  follows (see eq. (REF )), an analytic form for \(\mathfrak {R}\)  can be very difficult to find since \(\Xi \)  passes through a couple of nonlinear transformations, namely \(T\)  and \(\mathbf {f}\) .
Let us assume that \(\mathfrak {R}\)  is not known analytically and we cannot numerically generate samples from it directly. Thus we have to estimate the distribution from the data. Algorithm REF  summarizes the basic computations for computing a reduced order model. It returns the deterministic part of the model, the residual noise sequence, the modal noise sequence, and the innovation noise sequences. We note that while we are working with i.i.d. stochastic dynamical systems, and therefore should be working with the stochastic Koopman operator [1]}, (REF ) computes the reconstruction coefficients along a single realization of evolution of the dynamical system. Wanner and Mezić [2]} show the convergence of these coefficients to the true reconstruction coefficients of the stochastic Koopman operator using an ergodicity argument in the large data limit \(T\rightarrow \infty \) .
[H]
\({U_\xi ^t \mathbf {f}(\mathbf {x})}_{t=0}^{T}\)  stochastic evolution of observable
Deterministic model: \(D_\xi = {(c_j, \lambda _j, \mathbf {m}_j) : j=1,\dots J}\)
Noise sequences: \({(\mathbf {r}(t), \rho (t), \eta (t)) : t=0,\dots , T-1 }\)
Compute KMD of data \(\rightarrow ({\mathbf {m}_j}_j^{M}, {\lambda _j}_j^{M} )\)  *Assume Koopman modes are normalized.
Choose truncation number \(J > 0\)  according to some criterion
Compute modal reconstruction coefficients:
\(\mathbf {c} = \operatornamewithlimits{arg\,min}_{\alpha \in n} \sum _{t=0}^{T-1} { U_\xi ^t \mathbf {f}(\mathbf {x}) - \sum _{j=1}^J \alpha _j \lambda _j^t \mathbf {m}_j}_2^2.\) 
Construct deterministic reconstruction model
\(D_\xi ^t \mathbf {f}(\mathbf {x}) = \sum _{j=1}^J c_j \lambda _j^t \mathbf {m}_j.\) 
Compute residual sequence
\(\mathbf {r}(t) = U_\xi ^t \mathbf {f}(\mathbf {x}) - D_\xi ^t \mathbf {f}(\mathbf {x}).\) 
Compute modal noise sequence
\(\rho (t) = \mathbf {M} \mathbf {M}^+ \mathbf {r}(t).\) 
Compute innovation sequence
\(\eta (t) = \mathbf {r}(t) - \rho (t).\) 
Koopman Reduced Order Model for Stochastic Systems (KROM)

Forecasting statistics
Once the noise sequences \(r(t), \rho (t), \eta (t)\)  are computed, probability density models need to be constructed for which drawing samples is easy, say with multivariate kernel density estimation techniques. Assume the densities for \(\rho \)  and \(\eta \)  are computed.
Prediction in the subspace M.
Recall that \(M\)  is the subspace spanned by the Koopman modes, \(M = \operatorname{span}{\mathbf {m}_1,\dots , \mathbf {m}_J}\) , and \(P_M\)  is the orthogonal projection onto \(M\) . By (REF ) and since the innovation sequence \(\eta \)  is orthogonal to \(M\)  we have
\(P_M U_\xi ^t \mathbf {f}(\mathbf {x}) = D_{\xi }^t \mathbf {f}(\mathbf {x}) + \rho (t).\) 
Then
\(\mathbb {E} P_M U_\xi ^t \mathbf {f}(\mathbf {x})= \mathbb {E} D_{\xi }^t \mathbf {f}(\mathbf {x}) + \mathbb {E} \rho (t)= D_{\xi }^t \mathbf {f}(\mathbf {x}) + \mathbb {E} \rho (t).\) 
Thus the expected prediction in \(M\)  is the deterministic prediction but with a constant offset depending on the expected value of \(\rho (t)\) . Similarly, the variance of the predictions in the plane are
\(\operatorname{Var}[P_M U_\xi ^t \mathbf {f}(\mathbf {x})] = \mathbb {E}( P_M U_\xi ^t \mathbf {f}(\mathbf {x}) - \mathbb {E} P_M U_\xi ^t \mathbf {f}(\mathbf {x}))^2 = \mathbb {E}( \rho (t) - \mathbb {E} \rho (t))^2 = \operatorname{Var}(\rho (t)).\) 

Both of these results make sense since we have structured our model to have additive noise.

Performance and Comparisons
Linear Modal System
We first test the Koopman Reduced Order Model (KROM) method for a system with known modes and eigenvalues. The general model of the system is given by
\(\mathbf {x}(t) = \left(\sum _{j=1}^{J} \mathbf {m}_j \lambda _j^t \right) + \xi (t)\) 
where \(\mathbf {m}_j \in  \mathbb {R}^n\) , \(\lambda _j \in , \) t \(\mathbb {N}\)\(, and \) (t)  \(\mathbb {R}\) n\( is a noise sequence with each coordinate independent with a normal distribution. We note that the evolution of this system is in \) n\( rather than \)  \(\mathbb {R}\) n\( due to the choice of complex \) j\(^{\prime }s; therefore, there will not be complex conjugate pairs.\) For the following numerical results, the system parameters were chosen as \(J = 10\) , \(n=20\) , and \(\xi _i(t) \sim N(0, 0.25^2)\) . The values of the coordinates of each mode \(\mathbf {m}_j\)  were chosen from a uniform distribution between -1 and 1 and afterwards the mode was normalized to have norm 1. The associated eigenvalues were chosen uniformly at random in the box \([-1, 1] \times [-1i , 1i]\)  in the complex plane, then scaled to be on the unit circle.
The trajectory was simulated for 500 time steps resulting in a data matrix \(\mathbf {D} \in {20\times 500}\) . Only the first 20 snapshots are used for computing the reduced order model. The rest of the snapshots are used to compute the modal noise and innovation noise.
Figure REF  compares the true eigenvalues of the model which generated the data and the eigenvalues computed when constructing the reduced order model. Overall, the computed ROM eigenvalues are close with the discrepancy due to both the ROM having half the modes as the true model and the additive noise. Figures REF  – REF  show the corresponding modal noise distributions (noise in plane to the ROM modes), the innovation noise distribution (orthogonal to the ROM modes), and a comparison of the reconstruction of the signal for just the first coordinate of the state vector. Note that the distributions for both the modal and innovation noises are approximately Gaussian. Additionally, the ROM signal matches the general trend of the true signal, with the true signal lying within the 95% confidence interval of the ROM model. Here, we have shown only a single coordinate of the ROM and true signal, but these results are typical for the other coordinates. Given the close matching of the eigenvalues and the true signal lying in the 95% confidence interval of the ROM signal, we are confident in the methodology and move onto different examples.
<FIGURE><FIGURE><FIGURE><FIGURE>
Switched Anharmonic Oscillators with additive noise
We apply the KROM methodology to a switched anharmonic oscillator example. Additionally, we incorporate an additive Gaussian noise term to the action variables. The anharmonic oscillator was chosen due to its lack of point spectra other than the eigenvalue at 1 corresponding to the constant eigenfunctions. It can be shown [3]} that the rest of the spectrum is purely continuous with “eigenmeasures” \(\phi _j(I, \theta ) = e^{ij\theta }\delta (I - c)\)  supported on level-sets \({(c, \theta ) : \theta \in S^1}\)  of the action variable having eigenvalues \(e^{ict}\) . In this example, we wire a collection of these anharmonic oscillators together, periodically changing their wiring topology, and investigate our prediction algorithm on such a system. Below, we give the construction of the model.
Let there be \(j=0, \dots , n-1\)  oscillators with action/angle dynamics given by
\(I_j(t) - I_j(k) &= \sum _{i} \mathbb {I}_{E(t)}((i,j))\, c_{j,i}( I_i(k) - I_j(k) ) + \eta _j(t), \qquad ( k < t \le k+1)  \\\dot{\theta }_j(t) &= f( I_j(t) ) \) 

where \(f(0) = 0\)  and \(f\)  is a nondecreasing function, \(0 < c_{j,i} < 1\)  and \(c_{i,j} = c_{j,i}\) , \(E(t)\)  is the set of undirected edges denoting the connections between oscillators at time \(t\) , and \(\mathbb {I}_{E(t)}\)  is the indicator function for the set \(E(t)\) , and \(\eta _j(t) \sim N(0, \sigma ^2)\)  is a normally distributed noise sequence . For example, \(f(I) = I\)  gives the anharmonic oscillator. Physically, the dynamics are as follows. We have a collection of nonlinear oscillators whose frequency increases as energy is injected into the oscillator. At any given time there are random connections between the oscillators, and energy flows from the oscillator of higher action (energy) to the one with lower action. Note that this is an instantaneous injection of energy at time \(t = k + 0^+\) . In everything that follows, all edges are undirected edges.
Edge set dynamics.
Here, the edge set dynamics are given by permutation dynamics, although (REF ) allows for more random network connections. Let \(\mathbf {p}(t), \mathbf {q}(t) \in \Theta (n)\)  where \(\Theta (n)\)  is the set of all permutations of the vector \((0,\dots , n-1)\) . Write \(\mathbf {p}(t) = (p_0(t),\dots , p_{n-1}(t))\)  (similarly for \(\mathbf {q}(t)\) ). Let \(\mathcal {P}(n) = \lbrace  \mathbf {P}:\Theta (n) \rightarrow \Theta (n)\rbrace \)  be the set of all permutation matrices on \(\Theta (n)\) . Fix \(\mathbf {P}, \mathbf {Q} \in \mathcal {P}(n)\) . Then
\(\mathbf {p}(t) &= \mathbf {P}^{{t}}\mathbf {p}  \\\mathbf {q}(t) &= \mathbf {Q}^{{t}}\mathbf {q}  \\E(t) &= { (p_0(t), q_0(t)) } \) 

Here, the edge set changes at integer multiples in time and there is one random connection between oscillators (with it being possible an oscillator is connected to itself). The edge set consists of a single, undirected edge which connects oscillators \(p_0(t)\)  and \(q_0(t)\)  at time \(t\) .

Experiments.
We show results for a range of the values so that the dynamics have a separation between scales between the speed of rotation of the oscillators and the time scale of switching the topology. Since the action variables \(I\)  directly control \(\dot{\theta }\) , then these variables can be changed to accomplish this separation of scale. The action variables \(I\)  will be distributed according to an exponential distribution. The \(\lambda \)  parameter for the exponential distributions (see eq. (REF )) will be varied.
Let the state space be
\(X = { (\vec{I}, \vec{\theta }) \times (\mathbf {p}, \mathbf {q}) : I \in ( \mathbb {R}^+)^n, \vec{\theta }\in ( \mathbb {R}/ \mathbb {Z})^n, \mathbf {p}, \mathbf {q} \in \Theta (n)}.\) 
and let the dynamics be given by Eqs. (REF )–().
We will specify a distribution for the initial conditions. The exponential distribution is defined as
\(h(y; \lambda ) = {\left\lbrace \begin{array}{ll}\lambda e^{-\lambda y}, & y \ge 0 \\0, & y < 0.\end{array}\right.}\) 
The exponential distribution has mean \(1/\lambda \) . For \(\vec{I} = (I_1,\dots , I_n) \in ( \mathbb {R}^+)^n\) , we will take the \(n\) -product of exponential distributions
\(H(\vec{I}; \lambda ) = \prod _{j=1}^{n} h(I_j; \lambda )\) 
The distribution for the initial conditions will be
\(d\mathbb {P} = H(\vec{I}; 1) \otimes d\vec{\theta }\otimes d\mathbf {p} \otimes d\mathbf {q}.\) 
where \(d\vec{\theta }\) , \(d\mathbf {p}\) , and \(d\mathbf {q}\)  are the uniform distributions on their respective spaces. Draw initial conditions independently from the distribution:
\(x^{(s)} = (\vec{I}^{(s)}, \vec{\theta }^{(s)}, \mathbf {p}^{(s)}, \mathbf {q}^{(s)}) \sim \mathbb {P}, \qquad (s=1,\dots , S). \) 
Fix the observable \(\vec{f} \in L^2(X; n)\)  as
\(\vec{f}(\vec{I}, \vec{\theta }, \mathbf {p}, \mathbf {q}) = ( I_1, \dots , I_n, \theta _1, \dots , \theta _n) \in ( \mathbb {R}^+)^n \times ( \mathbb {R}/ \mathbb {Z})^n \subset  \mathbb {R}^{2n}.\) 
\(\vec{f}\)  does not directly observe the permutation variables \(\mathbf {p}, \mathbf {q}\) . For this observable, compute a KMD for each sample of the initial conditions \(x^{(s)}\) . The following sections specify the lambda parameters used. When applying the DMD algorithm, we use a delay embedding of the observable \(\vec{f}\)  . For all simulations, we use a delay embedding of 300 time steps resulting in a data matrix of size \(600\times 101\)  for the DMD algorithm to act on. We therefore will be able to compute a maximum of 100 modes.
<TABLE><FIGURE><FIGURE><FIGURE><FIGURE><FIGURE><FIGURE><FIGURE>
Kuramoto Model
As our final example for applying the KROM, we turn to the Kuramoto model of coupled oscillators. In the classical model (see eq.REF ), the strength of the connections between the oscillators is fixed and uniform for all oscillators. As the strength of the connection increases, the oscillators can display behaviors such as synchronization.
\(\frac{d\theta _i}{dt} = \omega _i + \xi _i + \frac{K}{N} \sum _{j=1}^{N} \sin (\theta _j - \theta _i).\) 
Here, we introduce an enhancement where instead of the fixed coupling strength \(K/N\)  we introduce randomness. Our enhanced model looks like
\(\frac{d\theta _i}{dt} = \omega _i + \xi _i + \sum _{j=1}^{N} \frac{K + \zeta _{i,j}}{N} \sin (\theta _j - \theta _i).\) 
Here, each \(\zeta _{i,j}\)  is distributed according to \(N(0, \sigma ^2)\) , with \(\sigma = 1\)  and \(K=5\)  in our examples. Additionally, \(N=10\)  is the examples. The \(N\times N\)  coupling matrix \(Z = [\zeta _{i,j}]_{1\le i,j \le N}\)  is restricted to be symmetric with \(\zeta _{i,j} = \zeta _{j,i}\) . In our experiments, we have chosen to set \(\xi _i \sim N(0, 0.25)\)  for all \(i\) . The natural frequencies were drawn uniformly at random from the interval \([0.25, 0.75]\) , whereas the initial condition for each coordinate of \(\theta \)  was chosen uniformly at random from the interval \([0, 2\pi ]\) . When constructing the data matrix we use a complex representation of the angles (\(\theta \mapsto z = \exp (i \theta )\) ).
<TABLE>Figure REF  shows the reduced order models' eigenvalues for models with 20, 50, 70, and 100 modes used for reconstruction. The set of eigenvalues form a strictly nested sequence of sets parameterized by the number of modes used for reconstruction. Specifically, the set of eigenvalues of the 20 mode model is a strict subset of the set of eigenvalues for the 50 mode model which is a strict subset of 70 mode model's eigenvalues, etc. This is a general property of the Koopman ROM algorithm since the method uses the data to compute all of the modes, orders them according to their norm, and then selects a subset in order to build a reduced order model. Another property to note is the distribution of the eigenvalues in the complex plane. The Type 2 models are distributed in a wedge \(W = {re^{is} : \forall r \in [0,1], \forall s\in (-\phi , \phi ), \phi \in [0,\pi ]}\)  with \(\phi \)  increasing with increasing number of modes used for reconstruction.
<FIGURE>Figure REF  shows the reconstruction of the signal of the first oscillator. The ROMs predictions converge rapidly to the true signal as the number of modes increases. Additionally, the computed variance rapidly converges to 0. In fact, one cannot distinguish a \(\pm 2\)  standard deviation band when using 60 modes or greater. See the appendix  for a comparison of these reconstructions when using real angles \(\theta \in [0, 2\pi )\)  and when using the complexified angles \(z = \exp (i\theta )\) . That discussion shows that since the KMD algorithms minimize loss functions using Euclidean distances, a representation for a variable needs to be determined where the Euclidean distance in that new representation closely approximates the natural norm for the original representation of the variable.
<FIGURE>We quantify the KROM algorithm using two metrics. The first is the reconstruction error. This is simply the average “geodesic” distance between the trajectories. If the true trajectory is \(\Theta = {\theta (t) : t \in [0, T], \theta (t) \in [0, 2\pi )}\)  and the ROM trajectory is \(\Phi = {\phi (t) : t \in [0, 1, \dots , T], \phi (t) \in [0, 2\pi )}\) , then the “geodesic” distance between the trajectories is computed as
\(d(\Theta , \Phi ) = \frac{1}{T+1} \left(\sum _{t=0}^{T} {e^{i\theta (t)} - e^{i\phi (t)}}^2 \right)^{1/2}.\) 
Note this is not the true geodesic distance for \([0, 2\pi ) \mod {2}\pi \)  but is a good approximation.
Figure REF  shows the natural logarithm of the reconstruction errors for each oscillator, parameterized by the number of modes used for reconstruction. The error between the signal and the ROM is a smoothe function of the number of modes used for reconstruction. There is (generally) a significant decrease in the log error with increasing number of modes.
<FIGURE>While the log error between the signal and the reduced order model is a good metric, we also have access to the modal noise which we use to put error bounds on the predictions given by the ROM. We compute the percentage of time the true trajectory is within \(\pm 2\)  standard deviations of the ROM's prediction. This is computed as
\(P = \frac{1}{T+1} \sum _{t=0}^{T} \mathbb {I}_{[0, 2\sigma )}( d_g(\theta (t), \phi (t)) )\) 
where \(d_g(a, b)\)  is the trueAs opposed to the approximate geodesic distance (REF ) geodesic distance on \([0, 2\pi ) \mod {2}\pi \)  and \(\sigma \)  is the modal noise standard deviation.
Even when using a lower number modes to construct the model (50 modes or less), the ROM shows good residence times, (REF ), for most of the oscillators. With higher numbers of modes, the residence times remain above \(0.8\) . One thing to note is that the residence time is not monotonic with respect to the number of modes used for reconstruction. We see non-monotonic behavior for models utilizing 50 modes up to 100 modes (full reconstruction). We believe part of this is due to the way the standard deviation of the modal noise is computed. The modal noise is computed using the residual between the signal and the reconstruction. As a consequence, as the reconstruction becomes more accurate, the variance of the residual sequence reduces. As can be seen, especially in Fig. REF , the standard deviation in the modal noise (the green bands) rapidly tightens as the number of modes increases. This may account for the non-monotonic behavior. A true analysis requires computing the rate of decrease in the variance of the modal noise versus the number of modes used for reconstruction contrasted with the rate of reduction of the error between the signal and the reconstruction with respect to the number of modes used for reconstruction.
<FIGURE><FIGURE>
A heuristic for estimating the minimum number of modes to use for the ROM
Above, we have computed a sequence of reduced order models for 3 examples, from a minimum of 10 modes up to a full reconstruction (100 modes in all examples). The question remains how to choose the final number of modes to be used in the reduced order model (ROM); too few and the we get poor prediction accuracy, too many results in a larger and more expensive model.
In this section, we present a heuristic method to determine the minimum number of modes that must be taken. This method is based on testing the normality of the modal distribution that is computed for each ROM. The rationale is as follows. In our examples, we have structured dynamics plus added noise and a changing network topology. When we construct our model, we compute a deterministic model and then estimate the modal noise distribution. We assume that once we have adequately recovered the true structured dynamics with the compute deterministic model the modal noise will have a Gaussian distribution. If too few modes are used for the deterministic model, much of the structured part of the dynamics will be modeled as noise, destroying the normality of the modal distribution. By using too many modes in the deterministic model, the noisy part of the dynamics are subsumed into a deterministic model, possibly destroying the normality of the distribution.
Algorithm  specifies the heuristic. We note that this test can be computed automatically and does not require manually looking at histograms or quartile-quartile plots to determine this although we show such plots in Figures REF , REF , and REF  below.
[H]
The set modal noise distributions \(\rho _i\)  fro model \(i\) .
Minimum number of modes to use for ROM.
Let \(\rho _i = {\rho _{i,j}}_{j=1}^{J}\)  be the set of modal distributions for reduced order model \(i\) . \(\rho _{i, j}\)  is the modal distribution for the \(j\) -th coordinate of the system. The \(\rho _i\) 's should be ordered from smallest model to largest (e.g., \(i = 10, 20, \dots , 100\)  in our examples.)
For \(j=1,\dots , J\) , perform a Shapiro-Wilkes hypothesis test to compute \(J\)  p-values, \(p_i = {p_{i,1}, \dots , p_{i,J}}\)  (one for each modal distribution).
For each \(i\) , compute the mean and median of the set \(p_i\) .
Heuristic: The smallest ROM where both the mean and median are above 0.05 is the smallest number of modes one should take for the model.
Heuristic for determining the minimum number of modes used in the ROM.
<FIGURE><FIGURE><FIGURE>
Discussion and Conclusions
We have developed a method for reduced order modeling using Koopman operator theory that gives confidence bounds on the predictions of the reduced order model. This is accomplished in the following way. While the reduced order model is by necessity a finite process, the spectral expansion of the Koopman operator is infinite. The reduced order model represents the process with a finite number of Koopman modes. The rest of the dynamics are modeled as a noise process. The part of the noise process that is in the subspace modeled by the ROM's modes is called the modal process. This noise process is modeled as a Gaussian noise process and the estimated standard deviation of the noise noise is used to compute a confidence bound on the ROM's predictions.
We have applied this modeling to a sequence of examples. The first was a synthetic example where we could directly specify the Koopman eigenvalues and modes, so the computed Koopman spectrum could be compared to ground truth. The spectrums did not match precisely since Gaussian noise was introduced during the evolution of the model. The next two examples represented networked dynamical systems that had changing network topology. The first was a network of noisy, anharmonic oscillators whose connections were switch at each integer time. This system was chosen due to the fact that a single, non-noisy, anharmonic oscillator has a purely continuous spectrum, except the eigenvalue at 1, and thus does not have an expansion in terms of Koopman modes and eigenfunctions. Despite this, the methodology was able to model the system's evolution with a small number of computed modes plus the noise process. The actual trajectory of the system stayed within the confidence bounds of the ROM's prediction a large percentage of the time. The last example consisted of a noisy Kuramoto model with random coupling strength between the oscillators. As with the anharmonic oscillator example, the true signal stayed within the confidence bounds of the ROM's prediction the majority of the time.
In this paper, we assumed a Gaussian distribution of the modal noise when constructing the confidence bounds for the ROM predictions. The variance of a Gaussian function fitted to the modal noise distribution was then used to construct the confidence bounds. While a fairly good approximation for the examples chosen, other examples might not have a modal noise distribution that can be well-approximated with a Gaussian function. In future work, we could apply a more general kernel density estimation of the modal noise and use this more general form of the noise to get the confidence bounds.

Appendix: KMD and discontinuous representations of continuous variables.
In our work with the Kuramoto model, we noticed that the reduced order models had difficulty resolving discontinuities in the angle variables \(\theta \in [0, 2\pi ) \subset  \mathbb {R}\) , which is indicative of a more general pitfall when applying KMD to problems. These discontinuities were artificial due to a poor choice of representation. The angles truly live on the torus \( \mathbb {R}/2\pi  \mathbb {N}\equiv [0, 2\pi ) \mod {2}\pi \)  so that 0 and \(2\pi \)  are identified. Any reconstruction method needs to use a metric that is natural to this space. However, the internals of the Koopman ROM method (and the internal KMD algorithm) use the Euclidean metric. This unfortunately results in an artificial discontinuities in the angle variables; e.g. two angles \(\theta _1 = 0\)  and \(\theta _2 = 2\pi - 0.1\)  would be have a distance of \(2\pi - 0.1\) , whereas it should have the true distance of 0.1. These artificial discontinuities produced problems, both in the interval KMD algorithm and other areas in the Koopman ROM code.
It is more natural to use a complex representation of the angles (\(\theta \mapsto z = \exp (i \theta )\) ) due to the trajectory being continuous in the Euclidean metric of the complex plane. The Koopman ROM algorithm would then be approximating a continuous complex, trajectory and the Euclidean metric \({z} = {z}\)  approximates the geodesic distance on \([0,2\pi ) \mod {2}\pi \) .
this was tested for reconstruction models using \(10, 20, \dots , 100\)  modes for reconstruction. The performance of these models were compared to the performance of the analogous models built using the non-complexified \(\theta \)  variables. We will use the following terminology to differentiate between these different types of Kuramoto models

Type 1 models: Kuramoto models built using the non-complexified \(\theta \)  variables.

Type 2 models: Kuramoto models built using the complexified variables \(z = \exp (i\theta )\) .

Figures REF  and REF  compare the reconstructions of Type 1 and Type 2 models. As the number of modes used for the ROM increases, the Type 2 models show much faster convergence of the ROM predictions to the true signal and much tighter confidence bounds as compared to the Type 1 models.
<FIGURE><FIGURE>