Machine-aided guessing and gluing of unstable periodic orbits
Pierre Beck Jeremy P. Parker Tobias M. Schneider





pierre.beck@epfl.ch

Emergent Complexity in Physical Systems Laboratory (ECPS), École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland



jparker002@dundee.ac.uk

Division of Mathematics, University of Dundee, Dundee DD1 4HN, United Kingdom



tobias.schneider@epfl.ch

Emergent Complexity in Physical Systems Laboratory (ECPS), École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland





§ ABSTRACT 
Unstable periodic orbits (UPOs) are believed to be the underlying dynamical structures of spatio-temporal chaos and turbulence. Finding these UPOs is however notoriously difficult. Matrix-free loop convergence algorithms deform entire space-time fields (loops) until they satisfy the evolution equations. Initial guesses for these robust variational convergence algorithms are thus periodic space-time fields in a high-dimensional state space, rendering their generation highly challenging. Usually guesses are generated with recurrency methods, which are most suited to shorter and more stable periodic orbits. Here we propose an alternative, data-driven method for generating initial guesses: while the dimension of the space used to discretize fluid flows is prohibitively large to construct suitable initial guesses, the dissipative dynamics will collapse onto a chaotic attractor of far lower dimension. We use an autoencoder to obtain a low-dimensional representation of the discretized physical space for the one-dimensional Kuramoto-Sivashinksy equation, in chaotic and hyperchaotic regimes. In this low-dimensional latent space, we construct loops based on the latent POD modes with random periodic coefficients, which are then decoded to physical space and used as initial guesses. These loops are found to be realistic initial guesses and, together with variational convergence algorithms, these guesses help us to quickly converge to UPOs. We further attempt to `glue' known UPOs in the latent space to create guesses for longer ones. This gluing procedure is successful and points towards a hierarchy of UPOs where longer UPOs shadow sequences of shorter ones.


Machine-aided guessing and gluing of unstable periodic orbits
    Tobias M. Schneider
    September 9, 2024
=============================================================





§ INTRODUCTION


It is widely accepted that unstable periodic orbits (UPOs) play an important role in supporting chaotic dynamics in many driven dissipative nonlinear systems. 
Periodic orbits are believed to be dense in the chaotic attractor of such systems and are organized in a hierarchical fashion, where longer orbits shadow a sequence of shorter ones. 
A well-known example displaying this hierarchical organization structure is the chaotic Lorenz ODE (ordinary differential equation) system with its famous chaotic attractor, which resembles the shape of a butterfly. Periodic orbits and chaotic trajectories shadowing them can be encoded by their sequential passage from one wing to the other. This leads to a description in terms of symbolic dynamics and reveals that periodic orbits are related to each other, with longer ones shadowing shorter ones <cit.>. The symbolic encoding, moreover, allows one to enumerate all UPOs. 

For chaotic PDEs, including the Navier-Stokes equations in the turbulent regime, we expect the same hierarchical organization of UPOs as in low-dimensional ODEs; and formal periodic orbit theory aimed at describing ergodic averages in terms of expansions over UPOs or, alternatively, over prime cycle sequences characterizing those UPOs, assumes it <cit.>. However, at least in the context of fluid flows, we are not aware of any direct demonstration of the hierarchical UPO organization. This fact highlights the algorithmic and numerical challenges inherent in identifying UPOs in high-dimensional chaotic systems. These computational challenges are especially severe for long periodic orbits that may be shadowing several shorter ones.

Physically relevant nonlinear dissipative chaotic PDEs, including model equations and the full Navier-Stokes equations, are usually formulated in 1-3 spatial dimensions, but their solution space is a formally infinite-dimensional function space. The spatial and temporal dimensions are typically discretized with many points, yielding a high-dimensional set of coupled ODEs. The identification of UPOs of the Navier-Stokes equations <cit.> suggests that a similar dynamical systems approach as in ODEs can be applied to PDEs, with spatiotemporal chaos being viewed as a chaotic walk between invariant solutions, such as equilibria, UPOs and invariant tori <cit.>. 
However, for 3D fluid flows, only very few UPOs have been identified due to the difficulties in computing them efficiently. Specifically, an envisioned hierarchical organization of UPOs has not been directly and conclusively demonstrated thus far.  

Much research has been conducted on the identification of UPOs, which is usually done in two steps: first by defining an adequate guess for a UPO and secondly by converging this guess to a solution of the system. The obvious and most common methods for finding UPOs are Newton shooting methods <cit.>. Here, the initial guess for the UPO is represented by a point in state space, namely the initial condition, and a period T. These two are then varied until the time-integrated trajectory closes in on itself. Typically, the optimization step is solved with a Newton algorithm. However, the exponential error amplification encountered when time-integrating a chaotic dynamical system leads to convergence issues, particularly when searching for long UPOs. 

More recently, loop convergence algorithms <cit.> and their matrix-free variations <cit.> have shown to be effective in finding UPOs <cit.>. The guess is now a space-time field that is already periodic (a loop) with time T but it does not satisfy the flow equations. The matrix-free variational methods from <cit.> deform the loop until its tangent vectors align everywhere with the flow vectors prescribed by the equations. This removes the time-integration aspect and consequently the challenges associated with an exponential error amplification characteristic of chaotic systems.  

In order to converge to UPOs (and in particular many distinct and long ones) we require a method to construct good guesses. Conventionally, guesses are extracted from recurrency methods <cit.>, where one looks for sub-trajectories in a long DNS of the system that almost close in on themselves. The downside of this method is that a trajectory is required to follow a UPO for an entire period, which is unlikely due to their unstable nature. As a result, this method is biased towards the same few frequently visited UPOs (usually short and less unstable ones). Although short periodic orbits are expected to have larger contributions in periodic orbit theory <cit.>, longer and more unstable periodic orbits are still necessary to obtain more accurate statistics <cit.>. They are also interesting from a control-theoretic point of view as they capture the dynamics and can be tracked for varying control parameter values <cit.>. Moreover, some dynamics appear to only be captured by long periodic orbits. <cit.> study the Kuramoto-Sivashinsky PDE, and for their control parameter of choice, the shortest UPO they find has period 12.08, while orbits that connect dynamically different parts of the chaotic attractor have periods of around 355.34 or more. Identifying guesses for such long UPOs that converge in the context of shooting methods is extremely difficult and requires unrealistically precise and extremely rarely observed recurrences. However, loop-based convergence algorithms that formulate a guess as a loop representing an entire space-time field have much more robust convergence properties than shooting approaches. This may allow to extract guesses within the convergence algorithm's convergence radius using alternative methods instead of a recurrency analysis. 

While the solution space is formally infinite-dimensional, trajectories in nonlinear, chaotic, driven, dissipative systems (such as the Navier-Stokes equations) collapse on a chaotic attractor once transients have died down. This attractor can be embedded in a curved manifold of far lower dimension - often termed the inertial manifold. Consequently, the high-dimensionality of the system's state space that renders the identification of UPOs so challenging, may be interpreted as an artifact of not knowing the most appropriate coordinates for describing the lower-dimensional intrinsic dynamics within the inertial manifold. If one had access to coordinates that approximately parametrized the lower dimensional (as compared to the discretization dimension) manifold that the attractor is embedded in, one could use these reduced coordinates to construct initial guesses for UPOs. In combination with the robust loop-convergence algorithms, even randomly drawn closed curves that lie in the inertial manifold and match the statistical properties of the attractor may be sufficient to define realistic guesses and identify UPOs. In analogy to the analysis of low-dimensional ODEs such as the Lorenz system, concatenating sequences of short UPOs to formulate guesses for long UPOs within such reduced coordinates may further allow to construct the hierachical sets of UPOs for PDEs that are expected to exist in theory, but so far have not been demonstrated directly. 

Linear model order reduction methods such as Dynamic Mode Decomposition (DMD) <cit.> and Principal Component Analysis (PCA) <cit.> (more commonly known as Proper Orthogonal Decomposition (POD) in the fluid dynamics community) are very popular for dimensionality reduction. While they capture a great deal of information, they are known to generalise less well to highly nonlinear systems and are outperformed by nonlinear deep learning methods such as autoencoders <cit.> <cit.>. The fact that autoencoders in particular seem efficient in giving a low-dimensional representation of spatiotemporal chaos is demonstrated by <cit.>. They train an autoencoder to identify low-dimensional embeddings of monochromatically forced Kolmogorov flow. They find that even for very low latent dimensions, such as N_h = 3, they obtain small losses, and come to the conclusion that much of the dynamics, such as low-dissipation events, live in a low-dimensional space. Even high-dissipation events are captured by only slightly larger latent dimensions, such as N_h = 32. Within the low-dimensional latent space, <cit.> re-define a recurrence function to obtain guesses for periodic orbits. <cit.> also explore autoencoders and combine them with a neural network in the latent space for time-series prediction. In particular, they study how the quality of the autoencoder improves as the latent dimension approaches the manifold dimension of the chaotic attractor.  

We propose to use an autoencoder to obtain an approximation of the low-dimensional manifold coordinates. Inside the latent space defined by the autoencoder, we randomly define loops that are in statistical agreement with the attractor. To this end, we define the loops as linear combinations of the latent POD modes with random periodic coefficients chosen to match the moments of the latent flow. We can arbitrarily adjust the complexity and length of these loop guesses and target longer UPOs by increasing the number of `twists' or `turns' in the loop. Moreover, the low-dimensionality of the latent space allows us to concatenate or `glue' orbits together. By gluing orbits, we define longer and more accurate guesses in a hierarchical fashion as has been observed for ODEs but not for PDEs, to the best of our knowledge. 

The structure of this paper is as follows: in section <ref> we give a brief reminder of loop-based convergence methods and introduce our setup for illustrating the methods in the Kuramoto-Sivashinsky PDE in chaotic and hyperchaotic regimes. In section <ref> we describe our methods used for defining initial guesses and describe the complete convergence setup for finding UPOs. We also define the notion of latent gluing and explain how we generate new, longer guesses by concatenating shorter UPOs. In section <ref> we first apply these methods to Kuramoto-Sivashinsky for parameter value L = 39 (low-dimensional chaos), and then in the hyperchaotic case at L = 100. We discuss and conclude in section <ref>. 





§ BACKGROUND





 §.§ Loop convergence methods


<cit.> treat the general PDE ∂_tu = F(u) for a real field u(x,t) on a n-dimensional spatial domain 𝒳⊂ℝ^n with initial condition u_0. The flow function f^t advances the dynamical system in time u(t) = f^t(u_0) = u_0 + ∫_0^tFdt' where t is the time. A fixed point u^* is a solution that satisfies F(u^*) = 0. A periodic orbit is characterised by an instantaneous field u (the initial condition) and a period T>0 that satisfy f^T(u) - u = 0 such that for any T^* < T this equation is not satisfied.

In shooting methods, a guess for a periodic orbit consists of an initial condition u_0 and a period T. The pair (u_0, T) are then varied until the time-integrated curve closes in on itself. This is done by solving the equation f^T(u_0) - u_0 = 0, typically via Krylov subspace methods <cit.> including Newton GMRES hook-step methods <cit.> and variations <cit.>. More recently, <cit.> define a cost function that is the norm of this equation and use gradient-based optimization to minimize it until a root is found.


In the loop convergence algorithm of <cit.>, the guess consists of an entire space-time field u(x, t), defined on 𝒳× [0,T)_periodic that is time-periodic with a guess period T, but does not necessarily satisfy the evolution equations of the system. A priori, the period T is unknown, and hence the field is re-scaled such that s = t / T and ũ(x, s) := u(x, sT). Hence ũ is a function defined on 𝒳× [0,1)_periodic mapping to ℝ^n. A solution ṽ(x, s) of the system then satisfies the re-scaled equation


    1/T∂ṽ/∂ s = F(ṽ)


Defining the residual vector r for a loop ũ(x, s) by


    r = F(ũ) - 1/T∂ũ/∂ s


we obtain the cost function J of a loop:


    J := ∫_0^1∫_𝒳r·r drds


A loop with J = 0 satisfies the flow equations and is a closed curve, and hence a periodic orbit. Conceptually, the loop is deformed until it satisfies the flow equations. Geometrically, the tangent vectors of the loop are aligned with the velocity vectors of the flow.

<cit.> minimize the cost function with a variational approach. While the method is very robust <cit.>, it involves a large Jacobian and therefore does not scale to high-dimensional systems, a challenge laid out explicitly by <cit.> and <cit.>. Inspired by a similar approach on equilibria by <cit.>, <cit.> formulated a matrix-free adjoint-based method (which we will use in this paper) for minimizing J over the space of loops which scales to high-dimensional systems, like the Navier-Stokes equations. They use the Kuramoto-Sivashinsky equation as a test-bed, while <cit.> apply a similar method to 2D Kolmogorov flow and explicitly address incompressibility in different ways. <cit.> deal with the challenge of computing pressure in the presence of solid walls in the 3D Navier-Stokes equations and introduce a method to accelerate the adjoint-based variational method with DMD. <cit.> also show that a similar approach can be used to compute connecting orbits through the construction of an analogous cost-function that undergoes an adjoint-based minimization process.





 §.§ The Kuramoto-Sivashinsky equation


The 1D Kuramoto-Sivashinsky equation (KSE) is a nonlinear PDE which arises in the modelling of the evolution of viscous liquid films down a vertical plane <cit.>, reaction-diffusion systems <cit.>, and flame-fronts <cit.> and exhibits chaotic behaviour for certain parameter values. In this paper, we use the following non-dimensional formulation of the KSE:


    u_t + uu_x + u_xx  +u_xxxx = 0

where u is the velocity field and we assume that the spatial domain is L-periodic, such that u(x,t) = u(x+L,t). The domain length L is the control parameter. The system is invariant under spatial and temporal translation, as well as under reflection x→-x, u→-u. We will work in the anti-symmetric subspace u(x) = -u(-x) (denoted by 𝕌^+ in <cit.>), as was done for example in <cit.> and <cit.>. This discretizes the spatial translation invariance, and reduces it to x → x + L/2. The system is also invariant under Galilean transformations, however this is filtered out by the imposed anti-symmetry condition. 



  §.§.§ Low-dimensional chaos : L = 39

Initially, we set L = 39, as in <cit.>, for which low-dimensional chaos is observed. Although this is not the simplest chaos observed for the KSE <cit.>, it is simple enough for us to show-case our methods before we move to a more complicated system. In this case, we discretize the spatial dimension with N_x = 64 points, turning the scalar function u into a 64-dimensional state vector u (not to be confused with the n-dimensional continuous field u in section <ref>). The data of a long trajectory (dt = 0.1, T_max = 155,000) is generated using the ETDRK4 scheme <cit.>. To ensure that the data is indeed just from the chaotic attractor, we cut off the first 50,000 time-steps. A space-time plot of an example trajectory generated in this setup is shown at the top of figure <ref>. 








  §.§.§ Hyperchaos : L = 100

In the second instance, we set L = 100, when the system is hyperchaotic with 5 positive Lyapunov exponents according to <cit.>. At L = 39, the dynamics are still relatively simple (albeit chaotic) due to the restricted symmetry and the narrow spatial domain. The more complex, hyperchaotic system resembles true spatiotemporal chaos, which is more akin to turbulence in fluids. In this case, we discretize the spatial dimension with N_x = 170 points, for which we observe a sufficient drop between the largest and smallest frequency of the time-averaged energy spectrum, while not over-resolving the system and making computations too slow. A typical trajectory of the system is presented at the bottom of figure <ref>, showing the increased complexity. Again, we generate one long trajectory where we cut off the first 50,000 time-steps, in order to make sure that all our data is part of one chaotic attractor. 




§ METHODS


In this section we first introduce the data-driven dimensionality reduction technique we use, namely autoencoders. We then explain how we generate loop guesses inside the autoencoder's low-dimensional latent space and describe the algorithmic procedure which we employ to converge to periodic orbits. Finally, we set out the gluing procedure to connect two existing periodic orbits, which serves as a guess for longer UPOs and helps us identify a hierarchy of UPOs.



 §.§ Data-driven dimensionality reduction




  §.§.§ Architecture

We apply an autoencoder to reduce the physical discretization dimension N_x to a latent dimension N_h. Autoencoders are neural networks that consist of two parts, namely the encoder ℰ: ℝ^N_in→ℝ^N_h and the decoder 𝒟: ℝ^N_h→ℝ^N_in. Given input data y∈ℝ^N_in, we would like to train the network so that approximately 𝒟∘ℰ (y) ≈y with N_h ≪ N_in.

The architecture of the autoencoder we employ in this paper varies slightly with the control parameter L, but the general architecture is similar: the encoder ℰ consists of three dense layers, with the third one having N_h nodes. The decoder has the same setup as the encoder, just in reverse. For the L = 100 case, we add four linear layers after the third dense layer of ℰ in order to minimize the rank of the latent space as described by <cit.> and explored by <cit.>. Illustrations of the networks are shown in figure <ref>. Note that we chose these architectures because they worked well for our purposes, but we do not claim that they are perfectly optimized. In both variations of the KSE explored later, we use ReLU activation functions on all layers except the linear layers: ReLU(x) = max{0,x}.







  §.§.§ Training

The numerical data obtained from direct numerical simulation (as described in <ref>) is re-scaled before it is used for training. We subtract the mean flow and then use a min-max re-scaling:

    u^* = (u - u_mean) - u_min/u_max-u_min + ϵ∈ [0,1)^N_x

where division is applied component-wise and ϵ = 10^-8 is a small constant to avoid division by 0. The minimum and maximum are vectors and also taken component-wise over u - u_mean. We drop the ^* for what follows for convenience. Due to the imposed anti-symmetry in the KSE and one component always being zero, we only use N_in = N_x / 2 - 1 components of u as inputs y. As loss function ℒ, we use the mean relative difference between the input y and the output 𝒟∘ℰ(y) rather than the standard mean-squared error, as dividing by the norm of y scales the loss in an interpretable manner. For N data points {y_n}_n=1^N, the loss is:


    ℒ = 1/N∑_n=1^N||𝒟∘ℰ(y_n) - y_n||^2/||y_n||^2 + ϵ


The best choice for N_h is a priori not known other than that we would like N_h ≪ N_in. We take topological quantities such as the Kaplan-Yorke dimension D_KY as guidance <cit.>. Note, however, that D_KY is a global average dimension of the attractor, and locally the topology might be more complicated. Since we are not looking for perfect guesses, D_KY is a good starting point. One way to identify an appropriate N_h is to train the autoencoder for multiple values of N_h and compare how the loss of the network evolves <cit.>. For each parameter choice of L and a selection of values for N_h, we train 20 autoencoders with different initializations for each value of N_h. Of the 20 networks, we pick the one with the best test loss. We then decide on the value of N_h by considering D_KY and by observing for which value of N_h we have considerable drops in the test loss.

In this section, we introduced the data-driven dimensionality reduction techniques that we will use to obtain an approximate low-dimensional representation of the high-dimensional discretized system. In the next section, we explain how the autoencoder's latent space gives us effective low-dimensional coordinates appropriate for defining loop guesses that are time-periodic space-time fields, lie on the chaotic attractor by matching the statistics of the flow and can be adjusted in length.



 §.§ Loops based on POD modes


We construct the guesses for UPOs based on linear combinations of the proper orthogonal decomposition (POD) modes <cit.> in the latent space with periodic coefficients. In general, consider a long time-series stacked in a matrix U∈ℝ^p× N, where the rows are p time-steps {u_i}_i = 1^p, and u_i∈ℝ^N. To calculate the POD modes, we compute the covariance matrix of the zero-mean time-series: let ũ_i = u_i - u, where u is the mean flow, and let Ũ∈ℝ^p× N be the corresponding zero-mean time-series. The unbiased estimator C for the covariance matrix is then given by 

    C = 1/p-1Ũ^TŨ∈ℝ^N× N


The POD modes ϕ_1, ..., ϕ_N are the eigenvectors of C

    Cϕ_k = λ_kϕ_k

with corresponding eigenvalues λ_1, ...,λ_N. Without loss of generality, the modes ϕ_1, ..., ϕ_N are ordered such that the eigenvalues are in decreasing order λ_1 ≥ ... ≥λ_N ≥ 0 (note that since the covariance matrix C is symmetric and positive semi-definite, its eigenvalues are real and positive).

The POD modes can be interpreted as fluctuations around the mean flow. Thus, we define a loop L(x, s) via a linear combination of the ϕ_k (x) with periodic coefficients a_k(s,{X_m,k}_m = 0^M)

    L(x, s) = u + ∑_k = 1^N a_k(s, {X_m,k}) ϕ_k(x)

where s is a periodic parameter and the {X_m,k}_m = 0^M are a sequence of independent and identically distributed (iid) random variables with a distribution X_m,k∼ X to be determined. 
We want the distribution of X to be so that loops on average match the first and second moments of the flow. Matching the first moment then means that the loop-average over the distribution of loops should agree with the mean flow. We denote this averaging by 𝔼_X,s[...], meaning 𝔼_X[⟨ ... ⟩_s], where ⟨ ... ⟩_s = 1/2π∫_0^2π...  ds. This gives the following two moment matching conditions:

    𝔼_X,s[L]    = u
    cov_X,s(L) C^(L)   = C

For our purposes, we assume that the X_m,k do not depend on s, and thus the order of integration between 𝔼_X and ⟨ ... ⟩_s does not matter. 
We match the first (equation <ref>) and second moment (equation <ref>) by setting 

    𝔼_X,s[a_k]     = 0 
    var_X,s(a_k)     = λ_k


for k = 1, ..., N. The detailed derivation of equations <ref> and <ref> is given in the appendix.

We now define the coefficients a_k. Since we want them to be time-periodic, we write the coefficients as a sum of sines and cosines

    a_k(s, A_:,k, B_:,k) = ∑_m = 0^M α_m [A_m, kcos(ms) - B_m, ksin(ms)]

where s∈[0,2π), M is the number of sine/cosine modes to be included in the sum, and the coefficients A_m, k, B_m, k∼ X are iid. The α_m are constants that give different weights of choice to higher frequency terms, for example α_m = (m + 1) / (M + 1). 

By substituting equation <ref> into equations <ref> and <ref>, we find that

    𝔼_X[A_m,k]    = 𝔼_X[B_m,k] = 0 
    var_X(A_:,k)    = var_X(B_:,k) = λ_k ( ∑_m = 0^M α_m^2)^-1


Thus, setting A_:,k, B_:,k∼𝒩(0, λ_k ( ∑_m = 0^M α_m^2)^-1) fulfills our requirements. Again, the detailed derivation is given in the appendix. 

The parameter M allows us to define longer guesses. Geometrically, adding higher modes to the sum in equation <ref> introduces extra `twists' or `turns' in our loop. This can be seen intuitively by thinking of Poincaré sections: the shortest UPOs are expected to have only p=1 intersections with an adequately chosen Poincaré section. UPOs with p = 2 intersections would require an extra twist in their geometric loop representation to intersect twice with the hyperplane. Thus, if we want to target UPOs with p Poincaré intersections, we set M = p. We will verify this intuition in section <ref>.

Defining loop guesses via equation <ref> allows us to generate random guesses by ad-hoc loops that are time-periodic and statistically lie on the attractor by matching moments up to second order. Moreover, the number of sine/cosine modes M allows us to adjust the length of the guess and therefore to directly target long UPOs. Nevertheless, the guesses themselves do not know anything about the dynamics, making it a crude approach for loop definition. While we will see that we can easily generate promising guesses this way and converge to UPOs, we will also present the method's current limitations. 




 §.§ Algorithm for guessing and converging loops to UPOs


In sections <ref> and <ref>, we laid out the methods for obtaining a low-dimensional representation of the high-dimensional discretized system and for generating guesses that are time-periodic space-time fields that match the flow statistics up to second order. Together with a loop convergence algorithm, we can now devise a general scheme for guessing loops and converging them to periodic orbits:

    
  * Obtain data {u_m}_m=1^M of the PDE of interest via direct numerical simulation.
    
  * Train an autoencoder for an adequate choice of N_h to get ℰ and 𝒟.
    
  * Obtain the latent POD modes ξ_1, ... ,ξ_K with eigenvalues γ_1, ..., γ_K based on a timeseries {h_i}_i = 1^p, where h_i∈ℝ^N_h, and K ≤ N_h is such that ∀ i = 1, ..., K we have γ_i > 0.
    
  * Define a loop L in the latent space following the approach from section <ref>:
    
    L(s) = h + ∑_k = 1^K a_k(s)ξ_k

    where 
    
    a_k(s) = ∑_m = 0^M α_m [A_m, kcos(ms) - B_m, ksin(ms)]

    
    and s∈ [0,2π). The coefficients A_m,k, B_m,k are randomly drawn from a normal distribution 𝒩(0, λ_k ( ∑_m = 0^M α_m^2)^-1).

    
  * Decode the loop to physical space 𝒟(L).
    
  * Use the adjoint solver from Azimi et al. <cit.> to converge the cost function J of the loop to J≈10^-4. Small weight adjustments are made to the gradient of the period T for stability.
    
  * Use a Newton solver to converge the loop to machine precision.

As a post-processing step, every UPO's time-resolution is increased to 256 time-steps and re-converged.



 §.§ Latent gluing of UPOs


Given two periodic orbits 𝒫_1 and 𝒫_2 with respective periods T_1, T_2, we will explore the possibility of concatenating or `gluing' them together in latent space and using this as an initial guess for a longer periodic orbit. The motivation behind this is symbolic dynamics, where a trajectory is described symbolically by its sequential passage through different parts of state space <cit.>. This appeals to a hierarchy of periodic orbits, where long UPOs shadow shorter ones <cit.>. Since this hierarchy appears to exist in ODE systems, we expect that this also applies in the KSE's physical and latent spaces. 

Given discretizations P_1∈ℝ^N_t_1× N_x and P_2∈ℝ^N_t_2× N_x of these orbits and latent representations L_i = ℰ(P_i)∈ℝ^N_t_i× N_h, we define the latent glued orbit G∈ℝ^ (N_t_1 + N_t_2) × N_x by first finding the indices I, J that minimize the distance between the two orbits in the latent space


    I,J = _i,j ||L_1^(i) - L_2^(j)||_2


where L_1^(i) and L_2^(j) are the i-th and j-th rows (or time-steps) of L_1, L_2 respectively. We define this minimal distance to be

    ℓ_2 = ||L_1^(I) - L_2^(J)||_2


The number of time-steps of each discretization are adequately chosen such that N_t_1/N_t_2≈ T_1/T_2. Define the naively glued orbit G_0 by vertically stacking L_1 and L_2 at the points of closest passage


    G_0 = 
        [                       L_1^(1:I);                 L_2^((J+1):end);                       L_2^(1:J);                 L_1^((I+1):end) ]


Since this introduces a jump discontinuity, we smooth G_0 in the latent space to get the new guess G: We set the high-frequency temporal modes in Fourier space to zero and keep only the lowest 1/6 positive (and lowest 1/6 negative) Fourier modes. The guess is then defined as 𝒟(G) with guess period T = T_1 + T_2, and then follows steps 6 and 7 of algorithm <ref>.




§ RESULTS


We apply the methods described in section <ref> to the KSE for two different parameter regimes: first L=39, for which low-dimensional chaos is observed, and then for the hyperchaotic case at L = 100.


 §.§ Low-dimensional chaos L = 39

Due to the imposed anti-symmetry in the system, half of the N_x = 64 discretization components are redundant. One of the remaining components is always zero. Thus the input vectors for the autoencoder have N_in = 31 dimensions. We train 20 autoencoders for each of N_h = 1,...,5 and find the final test losses shown in table <ref> and figure <ref>. As expected, the loss decreases when we increase N_h. The most significant drop is observed up to N_h = 3, which is in agreement with the Kaplan-Yorke dimension D_KY≈ 2.3 <cit.>. One might interpret this as measuring the dimension of the attractor <cit.>. 

We want to be able to visualize our methods in 3D, avoid overfitting and also test the effectiveness of the method. Although it might be too small to obtain an exact representation of the system, we decide to continue from here on with N_h = 3. Figure <ref> shows the performance of the autoencoder on test data. This trajectory was not part of the training data and is hence entirely new to the network. We observe that the network is able to identify the key structures of the trajectory. 

Figure <ref> visualises the attractor in the autoencoder's latent space together with a UPO in 3D (top) and in 2D projections (bottom). The points plotted in the figure are part of the data-sets used to train and test the autoencoder. 













  §.§.§ Periodic orbit searches

We generate multiple loops for different ranges of periods. We use all 3 latent POD modes to generate the loops as all 3 eigenvalues γ_1, γ_2, γ_3 are non-zero. As defined in section <ref>, we set the number of sine/cosine modes M to be equal to the targeted number of intersections p with an adequate Poincaré section. <cit.> uses û_1 = 0, (û_1)_t > 0 (where û_1 is the first component of the Fourier transform of u), in which case the dynamics appear to have a return time (the time between two consecutive intersections with the Poincaré section) of approximately 25 time units. Therefore, when we target short orbits with p = 1, we let M = 1 and choose guess period T = 25. For orbits with p = 2, we let M = 2 to introduce a `twist' in the loop and set the guess period to T = 50. For p = 3, we introduce two twists by choosing M = 3, and so on. For longer UPO searches (M≥ 3), we pick a range of periods between 25M and 25(M+1). Table <ref> in the appendix shows nicely how the number of Poincaré intersections of the final UPOs scales with M.

Figure <ref> compares decoded guess loops (of various lengths) to the UPOs and periods they converged to. We note that the guesses are realistic as they look like they could be trajectories of the KSE. This is emphasized when comparing the guesses to the final orbit they converge to: the guesses have similar sequences of patterns as the UPOs, and thus they look alike. In general, we find that many of our loops converge to periodic orbits, confirming that they are good initial guesses for loop convergence algorithms. The guesses look realistic and are already close to the UPOs that they eventually converge to. The detailed outputs of the runs for M = 1,2 and 3 are given in the tables <ref>-<ref>, where we generated 200, 500 and 700 loops respectively. We also indicate the number of times we converge to fixed points and when the algorithm does not converge. The latter happens either if we stopped the convergence algorithm too early or when the minimization of J gets caught in a local minimum, where J > 0 and ∇ J =0, instead of converging to a global minimum with J = 0. In cases where we converge to multiples of a short orbit (for example twice the 25.37 UPO) in higher-period runs, we cut the orbit into its shortest periodic component and re-converge. We conduct such searches for orbits with Poincaré intersections up to p = 4. Every UPO that we find is verified to still exist at temporal resolution N_t = 256. We note that for short orbits, we converge very often, with around 70% or 76% of guesses converging for targeted Poincaré intersections p = 1 and p = 2 respectively. As orbits get longer, the success rate naturally drops. While for p = 3, over half of the guesses still converge, when we target p = 4 this falls to around one quarter. These drops are to be expected as longer UPOs are generally harder to converge to, however we note that when we do latent gluing in section <ref>, the success rate for guesses with large periods shoots up significantly. A detailed overview of all runs is given in table <ref> in the appendix.










  §.§.§ Latent gluing


We create multiple new guesses using the methodology laid out in section <ref> by gluing the orbits 𝒫_i found above. We limit ourselves to gluing orbits with period T_i < 100 for computational efficiency reasons. We run through all possible combinations of gluing two orbits together. Every 𝒫_i has a symmetric counter-part 𝒫_i^s obtained by the shift x↦ x + L/2. Thus, we glue 𝒫_i with both 𝒫_j and 𝒫_j^s. Since we find 18 orbits with T<100, this gives a total of 306 combinations. Figure <ref> compares the distribution of random distances of a long time-series in latent space to the distribution of distances of closest passage ℓ_2 between two UPOs in latent space and confirms that the points of gluing of two UPOs are indeed close together.




For illustration purposes, figure <ref> shows a 2D projection of the gluing process in the latent space between two short orbits with periods T_1 ≈ 24.908 and T_2 ≈ 25.371. In this case, the gluing is easy to follow visually, and the glued guess with initial period T_1 + T_2 ≈ 50.279 converges quickly to a periodic orbit with T ≈ 50.368. One can easily see that the converged orbit shadows the initial two. Figure <ref> shows the same process for two longer initial periodic orbits with periods T_1 ≈ 83.804 and T_2 ≈ 85.559. The glued guess with initial period T_1 + T_2 ≈ 169.363 looks very similar to the converged orbit with period T ≈ 169.467. This is even more apparent in figure <ref>, which shows the decoded, physical plots of these orbits, comparing the glued loop to the final periodic orbit. Again, the converged long orbit appears to shadow the initial two short ones. Figure <ref> confirms that for the glued guesses that converged, the final period is approximately equal to the sum of the initial periods.









Out of the 306 total guesses, 227 converge to UPOs, 160 of which are distinct, and 79 do not converge (either the optimization of J gets stuck in local minima or we need to run the convergence for longer). The largest period found this way is T ≈ 171.096 and we find a general success rate of 74.2%. Impressive is also the success rate for long UPO guesses: loops with guessed periods T_1 + T_2 > 100 converge in 210 / 284 ≈ 73.9% of cases, compared to the approximately 25% from purely random guesses (with M = 4) observed in the previous section. We note that some UPOs struggle with being glued to other orbits, such as those with periods 53.135, 57.227, and 57.627 (see table <ref> in the appendix). 

Since we expect this hierarchy of periodic orbits, where long orbits shadow shorter ones, we indeed also expect the glued guesses to perform much better than the random ones. Doing the gluing process directly in physical space might also work just as well. The key take-away is that the expected hierarchy of UPOs appears to be present in this PDE system, and also in the latent attractor. This is a good confirmation that the autoencoder is able to capture a coherent, low-dimensional representation of the chaotic attractor of the KSE obtained through a nonlinear dimensionality reduction technique. Moreover, when the ℓ_2 distance between the two points of closest passage of two UPOs in latent space is small, then the convergence rate is larger. We obtain a convergence rate of 82.2% for guesses with ℓ_2 < 0.07. These cover over 50% of the guesses we attempted. Guesses with ℓ_2 ≥ 0.07 converge in only 65.8% of cases, showing that this ℓ_2 is a good indicator of whether orbits are gluable. Figure <ref> shows the cumulative convergence rate against ℓ_2.



Out of the total 153 glued symmetry pairs, in 14 cases neither glued guess converges. In 51 cases, we find that while one glued guess does not converge the other converges to a periodic orbit. We also observe that while 39 pairs converge to the same orbit, 49 pairs converge to two distinct ones. For the full output details, see table <ref> in the appendix.


In summary, in this section we applied the methods described in section <ref> to the KSE in the case of low-dimensional chaos. First, we generated guesses for UPOs by sampling random closed curves in the latent space that on average match the latent flow statistics up to second moments. By varying the number of sine/cosine modes in the linear combination of POD modes, namely the parameter M, we observed that the resulting UPOs usually had p = M intersections with the Poincaré section. This allows us to directly target longer UPOs by increasing M. Next, we glued UPOs together in the latent space at their closest points of passage. The resulting new guesses had very high convergence rates, indicating that the hierarchy of UPOs observed in ODEs is also present here. This also gives a method to search for longer UPOs. In the next section, we will apply these methods to the hyperchaotic case with L = 100.




 §.§ Hyperchaos L = 100

Taking the Kaplan-Yorke dimension D_KY≈ 9.2 <cit.> of the system as guidance, we train 20 autoencoders for each of the latent dimensions N_h = 8, ..., 14. Figure <ref> shows the best test losses for each value of N_h. We note that there is a strong drop until N_h = 11, followed by another strong drop at N_h = 13. It is likely that, when training more N_h = 12 models, we would find a test loss that continues the exponential trend observed for N_h = 9, 10, 11, 13. By considering the test losses, as well as the D_KY, we decided to continue with the more parsimonious N_h = 11. Figure <ref> shows the performance of the autoencoder applied to a physical trajectory part of the test set. We note that the yellow spots (where the relative error is ≥ 2) are more frequent and larger, indicating that the autoencoder struggles more with the hyperchaotic system. Nevertheless, the autoencoder is able to reconstruct the general shapes and structures of the system, and thus shows satisfactory results. Moreover, the latent POD decomposition only has 9 non-zero eigenvalues, resulting in an effective 9 latent dimensions that we use to define our guesses.








  §.§.§ Guesses for the period

In the L = 39 case, we related the period guesses to the system's approximate average return time of the Poincaré section. This becomes less straight-forward in higher-dimensional systems, such as the L = 100 case, as a Poincaré section of an n-dimensional system should be an (n-1)-dimensional subspace. Defining a guess period in this case is non-trivial, however we require one to be able to apply the convergence algorithm from <cit.>.
For a UPO of the system, it is possible to integrate around the UPO and obtain the time taken to traverse it. While a guess for a UPO is not a solution of the KSE, we do have access to ∂_t u at each point of the loop guess by evaluating the right-hand side of the KSE. Therefore, we decided to estimate the period guess of the loop by treating it as if it was a solution of the system and integrating around the loop the time taken to traverse it. More concretely, for a discretized loop time-series {u_t}_t = 1^N_t with N_t time-steps, the guess for the period is


    T_guess = ∑_i = 1^N_t || Δu_i||/1/N_t∑_i = 1^N_t || ∂_t u_i||


where we assume that the time-gap dt between time-steps is constant. The numerator can be interpreted as the length of the loop, and the denominator as the average velocity around the loop. Figure <ref> shows the distribution of T_guess for different values of M. Similar to the L = 39 system, we can aim for longer or shorter UPOs T_guess by introducing more "twists" in our loop i.e. by changing M.

Note that it would also have been a reasonable choice to use ∂_t u_i ·t̂ in the denominator of equation <ref>, where t̂ is the normalized loop tangent at u_i, instead of ||∂_t u_i||. The sign of this dot product would tell us whether the guess loop is going in the correct direction at u_i. However, since random vectors in high-dimensional spaces are likely to be orthogonal <cit.>, this results in T_guess blowing up and over-estimating the period.








  §.§.§ Periodic orbit searches

For the initial test searches that we conducted, we observed a sharp increase in T during the adjoint looping procedure related to the increase in parameter L. To counteract this and to better target the guess period, we reduce the weight of the period gradient and to gradually increase it again as the loss approaches 0. Moreover, we looked at different values for α_m = [(m+1)/(M+1)]^β by varying the exponent β that distributes the weighting between smaller and higher modes. We tried β = 1/2, 1, 2 and M, and found the best performance again with β = 1.

We conduct searches for UPOs based on guesses generated with M = 1,2 and 3. For each value of M, we generate 1,000 guesses. A summary of these three searches and their results is given in table <ref>. Examples of guesses generated from the autoencoder and the periodic orbits they converge to are shown in figure <ref>. As for the L = 39 case, we observe that the guesses and the periodic orbits that they converge to share similar shapes. A noticeable difference to L = 39 case is that during the convergence, the loops undergo a stretching and squeezing process with certain parts of the orbit being traversed faster or slower than for the guess, and the initial period often being an underestimate. This highlights the challenge of correctly parametrizing the initial guess inside the latent space. T_guess appears to underestimate the period, which may be due to a bias from our Riemann-sum approach in the definition of T_guess.





In general, we observe good convergence rates for short orbits, with 153 orbits converged for M = 1, of which 130 are distinct. We note that the converged periods create more of a continuous spectrum (rather than multiples of an average return time as in the L = 39 system), again reflecting the increased complexity of this system. The success rate naturally drops for M = 2 and 3, going to 49 and 11 orbits respectively (all of which are distinct). We attribute this drop to various factors, namely the much increased complexity of the system at L = 100, the ad-hoc definition for guesses based on sines and cosines without taking into account the dynamics apart from the moment matching, the difficulty in correctly parametrizing the speed of traversing the guess in the latent space, and how efficient the algorithm that we used is at converging long UPOs. Indeed, since long UPOs require more time-steps for accurate convergence, it is likely that if we let the convergence run for longer that more of our guesses would converge. Nevertheless, the ability to find many short orbits with such ease is a success in itself: the generation of these guesses only requires a one-time up-front cost in training the network. Once this is done, guesses can be generated cheaply and instantaneously.



From the three main searches we obtained 213 UPOs, 190 of which are distinct, with periods ranging from 10.03 to 110.11. As mentioned earlier, we also conducted other searches to examine different values of α_m, different discretization methods, and also initial searches without an adequately weighted period gradient. During these less successful searches, we also found other periodic orbits, giving us a total of 492 distinct orbits, with periods ranging from 9.96 to 110.11. All UPOs that we found were verified to still exist at temporal resolution N_t = 256 and spatial resolution N_x = 256. Figure <ref> shows the number of periodic orbits that we have found up to a given period. Based on <cit.> we expect the number of UPOs to increase exponentially for large periods. On figure <ref> we plot the exponential trend for orbits with periods T∈ [20,30]. For T > 30 we have not found enough UPOs to continue the trend. For such long UPOs, we expect the structure of the loop state space to be much more complex, implying the existence of many more local minima and thus less successful convergence rates (as also observed in <cit.>). 






  §.§.§ Latent gluing

 Since we have a set of 492 UPOs available to us, an exhaustive gluing approach would give us a total of over 241,000 glued guesses. We shall not attempt all of these for computational efficiency reasons. Using the results from section <ref>, we will restrict ourselves to the following two searches:
 
     
  * Search A: Guesses where the initial two orbits are close in the latent space, with ℓ_2 < 0.1, and where T_1 + T_2 < 125. This gives a total of 877 guesses.
     
  * Search B: A random selection of 1,000 glued guesses among those with T_1 + T_2 < 125.
  
 
We restrict ourselves to T_1 + T_2 < 125 for computational efficiency again, since the larger the expected period T, the higher the time-resolution should be and the longer it takes to compute these orbits. Physical plots of a typical example of a successful gluing between two UPOs are shown in figure <ref>. Figure <ref> shows the 2D projections of the initial orbits and the final one onto the first three latent POD modes. We note that the initial UPOs appear embedded within the glued UPO, in the sense that the long glued UPO shadows the shorter ones. 

Search A results in a success rate of approximately 9.8%, with the cumulative convergence rates shown in figure <ref>. In particular, looking at the 88 closest orbits (the closest 10%), the convergence rate is 18.2%. Search B, which consists of 1,000 random glued guesses, only has a success rate of 3.8%, clearly showing that a small ℓ_2 increases the likelihood of two orbits being able to be glued together. This indicates that also for the hyperchaotic system there seems to be a hierarchy of UPOs, where long UPOs shadow shorter ones, both in the physical space and in the autoencoder's latent space.

The success rate for L=100 is noticeably smaller than for L = 39. While this may be partly attributed to the higher complexity of the system, it also indicates that there might be potential for optimizing such gluing procedures. Of course, since glued orbits have a larger expected period, attempting to converge these orbits for longer might also increase the success rate.












§ CONCLUSION & DISCUSSION


In this paper we have introduced a new method for generating initial guesses for periodic orbits by randomly drawing loops in the low-dimensional latent space defined by an autoencoder. The autoencoder's latent dimension was chosen in a process where multiple networks were trained for various latent dimensions. We picked those networks for which we observed a strong drop with respect to the latent dimension (one can interpret this as autoencoders approximating the intrinsic coordinates of the manifolds that describe the chaotic attractor) and that performed the best among those. The loop guesses are constructed by linear combinations of the latent POD modes (that have non-zero eigenvalues) with periodic coefficients drawn from a random distribution that match the latent flow statistics. These loops are then decoded back to the physical space and together with a guess period serve as an initial guess for a loop convergence algorithm. We apply this method to the Kuramoto-Sivashinsky PDE in regimes of low-dimensional chaos and hyperchaos. The decoded loops lie close to the chaotic attractor, look realistic and prove to be good guesses for periodic orbits in loop convergence algorithms, with many of them converging to periodic orbits. This provides an alternative to recurrence methods, where guesses are based on near recurrences in a long DNS, and are thus rare and expensive to generate, while also being biased towards short and less unstable periodic orbits. We note that the derivation for loop guesses based on POD modes holds for both the physical and the latent space. For systems exhibiting low-dimensional chaos, defining guesses based on the physical POD modes is a valid approach and gives acceptable guesses. However, we observed that with guesses purely based on the physical POD modes we tend to converge to the same few UPOs. As the system becomes more complicated, it appears necessary to rely on a nonlinear order reduction method, like autoencoders, and define loops in the latent space via the latent POD modes to enforce a larger variation of loop guesses that approximately lie on the latent attractor.

Motivated by a hierarchy of UPOs that is present in ODE systems, such as the Lorenz equations, we explored the method of latent gluing, where we concatenate two periodic orbits based on where they are closest to each other in the latent space. After smoothing this new loop inside the latent space and decoding it back to the physical space, we use this as an initial guess for longer periodic orbits, as we expect the hierarchy of UPOs to carry over from ODEs to PDEs. Many of these guesses converge to UPOs, indicating that this hierarchy does indeed hold in the specific PDE system, with long UPOs shadowing shorter ones. This motivates further research into whether this hierarchy is also present in more complicated spatiotemporally chaotic PDE systems, such as the Navier-Stokes equations. Additionally, the gluing provides a method for generating new, longer UPOs from known shorter ones. Importantly, in both regimes of low-dimensional chaos and hyperchaos, the gluing is much more successful if the distance of closest passage between the two initial UPOs is smaller. This indicates that the hierarchy is also present in the autoencoder's latent space and that the autoencoder is able to capture a coherent low-dimensional representation of the system. For both the chaotic and hyperchaotic systems we conclude that a small distance between the points of closest passage of two periodic orbits in latent space significantly increases the likelihood of whether two UPOs are `gluable'. 

These results are a step forward for using loop convergence methods to find unstable periodic orbits, as we have shown that we can easily and cheaply generate good initial guesses. As expected, there are limitations when moving to long UPOs in hyperchaos. However, this also constitutes the outlook for the future: we believe that optimizing the autoencoder, as well as the guesses themselves so that they are built on knowledge of the dynamics will significantly improve this method. 
Having identified the hierachical organization of UPOs with long orbits shadowing a sequence of shorter ones moreover fuels the hope to generate symbolic encodings of those UPOs not only for simple ODEs but for formally infititely dimensional PDE systems and thereby - at least approximately - enumerate UPOs which is of key relevance when trying to transfer ergodic theory concepts from chaotic ODEs to dissipative nonlinear chaotic PDEs underlying physically relevant phenomena such as fluid turbulence.







§ ACKNOWLEDGEMENTS

The authors thank Omid Ashtari for his invaluable insights on loop convergence methods as well as his code. This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant no. 865677).






§ APPENDIX


tocsectionAppendix



 §.§ Training specifications

The networks were created with TensorFlow 2.11.0 and trained on a NVIDIA RTX A4000 GPU. The L = 39 networks had inputs of size 31 and the encoder part consisted of three dense layers with 32, 32 and N_h nodes respectively. The decoder has the inverse setup and consists of three dense layers of sizes 32, 32 and 31. Each layer has ReLU activation functions. We trained 20 networks for each N_h∈{1,2,3,4,5} for 500 epochs. 
The L = 100 networks had inputs of size 84 and the encoder part consisted of three dense layers with 256, 128 and N_h nodes respectively with ReLU activation functions, followed by 4 linear layers with N_h nodes for implicit rank minimization. The decoder consists of three dense layers of sizes 128, 256 and 84, each followed by a ReLU activation function. We trained 20 networks for each N_h∈{8,9,10,11,12,13,14} for 1,000 epochs.

For all networks we used the AdamW optimizer and an initial learning rate of 7*10^-4 with exponential decay. Each network also had l_2 regularization.



 §.§ Loops based on POD modes



  §.§.§ Matching moments


The a_k are independent as we assume the X_m,k to be iid. Since 𝔼_X,s is linear and the POD modes form an orthonormal basis, equation <ref> is satisfied if and only if 𝔼_X,s[a_k] = 0. Using this in equation <ref>, we find that


    C^(L)_ij = cov_X,s(L_i, L_j) 
       = cov_X,s(u_i + ∑_k = 1^N a_k(ϕ_k)_i, u_j + ∑_k = 1^N a_k(ϕ_k)_j) 
       = cov_X,s(∑_k = 1^N a_k(ϕ_k)_i, ∑_k = 1^N a_k(ϕ_k)_j) 
       = 𝔼_X,s[ ∑_p = 1^N a_p(ϕ_p)_i ∑_q = 1^N a_q(ϕ_q)_j] 
       = 𝔼_X,s[ ∑_p, q = 1^N a_p a_q(ϕ_p)_i (ϕ_q)_j] 
       = ∑_p, q = 1^N 𝔼_X,s[a_p a_q] (ϕ_p)_i (ϕ_q)_j 
       C^(L) =  ∑_p, q = 1^N cov_X,s(a_p, a_q) ϕ_p ϕ_q^T

Since the a_k are independent from each other, then cov_X,s(a_p, a_q) = 0 when p≠ q and cov_X,s(a_p, a_q) = var_X,s(a_p) if p = q, giving 

    C^(L)_ij = ∑_p = 1^N var_X,s(a_p) (ϕ_p)_i (ϕ_p)_j

Now note that C = VDV^T, where V has columns ϕ_1, ..., ϕ_N and D = diag(λ_1, ..., λ_N). Hence

    C_ij   = ∑_p,q=1^N V_ip D_pq V_qj^T 
       = ∑_p=1^N λ_p V_ip V_jp
     C_ij   = ∑_p=1^N λ_p(ϕ_p)_i (ϕ_p)_j

By comparing this expression to equation <ref>, we can match second moments by setting 

    var_X,s(a_k) = λ_k

for k = 1, ..., N.



  §.§.§ Deriving coefficients


We define the coefficients to be a sum of sines and cosines

    a_k(s, A_:,k, B_:,k) = ∑_m = 0^M α_m [A_m, kcos(ms) - B_m, ksin(ms)]

where M is the number of sine/cosine modes to be included in the sum, and the coefficients A_m, k, B_m, k∼ X are iid. The α_m are constants that give different weights of choice to higher frequency terms, for example α_m = (m + 1) / (M + 1). Then

    𝔼_X,s[a_k] = ∑_m = 0^M α_m{𝔼_X[A_m, k] ⟨cos(ms)⟩_s 
     - 𝔼_X[B_m, k] ⟨sin(ms)⟩_s}

Since the s-integrals are 0 for m>0, this only requires 𝔼_X[A_0,k] = 0. To simplify the variance calculation, we will set 𝔼_X[A_m,k] = 𝔼_X[B_m,k] = 0 for all m. The variance is given by

    var_X,s(a_k) = 𝔼_X,s[(∑_m = 0^M α_m(A_m, kcos(ms) 
                        - B_m, ksin(ms) ))^2] 
        = ∑_m,n = 0^M α_mα_n( 𝔼_X[A_m, kA_n, k] ⟨cos(ms)cos(ns)⟩_s 
            - 𝔼_X[A_m, kB_n, k] ⟨cos(ms)sin(ns)⟩_s 
            - 𝔼_X[B_m, kA_n, k] ⟨sin(ms)cos(ns)⟩_s 
            + 𝔼_X[B_m, kB_n, k] ⟨sin(ms)sin(ns)⟩_s)
        = α_0^2𝔼_X[A_0,k^2] + ∑_m=1^M1/2α_m^2[𝔼_X[A_m,k^2] + 𝔼_X[B_m,k^2]] 
        = var_X(A_:,k)∑_m = 0^Mα_m^2 
       var_X(A_:,k) = λ_k ( ∑_m = 0^M α_m^2)^-1


Thus, letting A_:,k, B_:,k∼𝒩(0, λ_k ( ∑_m = 0^M α_m^2)^-1), the loops are random, periodic and on average match the first and second moments of the flow.



 §.§ Supplementary tables











