Learning the tangent space of dynamical instabilities from data

The theory of dynamical systems has a long history of trying to elucidate one of the most important concepts in science and engineering: instability. Dynamical instabilities are important because they can give rise to a variety of phenomena with unexpected, and even disastrous, consequences. They occur in fluid mechanics,^1,2 climate dynamics,³ optics,⁴ and thermoacoustics⁵ and come in many shapes and forms. Perhaps the simplest of all are instabilities arising from a linear mechanism, whose investigation traditionally involves linearizing the governing equations around a fixed point and looking for unstable eigenvalues of the linearized problem.⁶ Generalization of this concept to periodic orbits led to the well-known Floquet theory, in which stability of periodic trajectories is ascertained by computing the spectrum of the monodromy matrix.⁶ The realization that episodes of transient growth cannot be predicted by linear or Floquet analyses came much later and in turn gave rise to the theory of non-normal and transient instabilities.⁷ This culminated with the introduction of the Lyapunov exponents and Lyapunov vectors, which provide a rigorous description of instability mechanisms in chaotic systems.^8,9 Since then, considerable effort has been devoted to the development of algorithms for computation of Lyapunov exponents and Lyapunov vectors from data. Great strides have been made by the likes of Eckmann et al.,¹⁰ Sano and Sawada,¹¹ Rosenstein et al.,¹² Kantz,¹³ and Wolfe and Samelson,¹⁴ who proposed new methods to compute Lyapunov exponents and Lyapunov vectors from experimental measurements or improved on existing ones.

The key feature shared by these algorithms is that they monitor the fate of perturbations around a reference orbit for sufficiently long times, sometimes resorting to orthonormalization in order to prevent blow-up of the perturbation vectors. However, tracking the evolution of perturbations along a trajectory is nothing more than a time-marching approach in disguise. More importantly, it does not utilize the fact that for a large class of dynamical systems, the $i$th Lyapunov vector $ui$ depends pointwise on the state of the system [with the mapping $x⟼ui(x)$ being uniquely determined by the phase-space point $x$], and consequently, the $i$th Lyapunov exponent—a quadratic form over $ui(x)$—can be recovered by measure-averaging over many measurement points, rather than time-averaging over a long trajectory. This was pointed out by Ershov and Potapov,¹⁵ who recognized the benefits of using measure-averaging in lieu of time-averaging, as the former allows for the measurement points to be arranged in arbitrary order and even for them to belong to different trajectories. These results show that if one were able to compute the pointwise function $ui(x)$ from data, then one would immediately have access to a complete picture of the directions of instabilities at any point in the phase space, as well as accurate estimates for the Lyapunov exponents.

As far as we know, the only attempt to compute the map $ui(x)$ was made in Ref. 16, where we used Fenichel's theory of slow invariant manifolds¹⁷ to derive analytical expressions for $ui(x)$ in situations exhibiting slow-fast dynamics. However, application of this method is limited because (a) the dimensionality of the system cannot be too large for the manifold analysis to be tractable and (b) the system must have a well-defined separation of time scales, short of which Fenichel’s theory becomes moot. Of course, the method being analytical, there is no data component to it. This is precisely what we set out to rectify in the present paper, with the introduction of a machine-learning algorithm that infers the map $ui(x)$ from a large collection of state measurements.

Machine-learning (ML) algorithms have pervaded virtually every area of science and engineering because of two reasons. First, the amount of data available from experiments, field measurements, or numerical simulations has reached unprecedented levels. Second, ML algorithms have proven to be extremely versatile and powerful from the standpoint of extracting patterns and information from tremendously complex systems that would otherwise remain inaccessible. Applications of ML to dynamical systems include dimensionality reduction and flow-feature extraction,^18–21 discovery of governing equations,^22–24 design of optimal control strategies,²⁵ turbulence closure modeling,^26,27 integration of partial differential equations,^28–30 and forecasting of dynamical instabilities in chaotic systems.^31,32 The recent emergence of deep learning has led to a push in the latter direction, with recurrent neural networks^33,34 and reservoir computing^35,36 leading the way.

The learning algorithm that we propose in this paper produces a cartography of directions associated with strongest instabilities in the phase space, from which the leading Lyapunov exponents can be extracted. These directions coincide with the backward Lyapunov vectors of Legras and Vautard⁹ and the optimally time-dependent (OTD) modes of Babaee and Sapsis.³⁷ (The equivalence between the two was established in Ref. 16.) The potential of the learning algorithm in problems related to prediction and control of transient instabilities and extreme events is considerable because the proposed method is fully data-driven, has no restricting assumptions other than invertibility, autonomy, ergodicity, and measure-invariance of the underlying dynamical system, and only requires state measurements as inputs.

The paper is organized as follows. We formulate the problem in Sec. II, introduce the learning algorithm in Sec. III, followed by results in Sec. IV, a discussion in Sec. V, and a conclusion in Sec. VI.

We consider the autonomous dynamical system

where $x$ belongs to a compact Riemannian manifold $X$ endowed with the Borel $σ$-algebra and a measure $μ$⁠, the map $F:X⟶X$ is sufficiently smooth to ensure existence and uniqueness of solutions, and overdot denotes differentiation with respect to the time variable $t$⁠. We assume that the transformation

sometimes referred to as the “flow map,” is invertible, measure-preserving, and ergodic. Measure-invariance is important from the standpoint of defining time-averages of scalar functions. (This is the well-known Birkhoff ergodic theorem.³⁸) Measure-invariance and ergodicity are important to guarantee that time-averages and measure-averages coincide,

for all $f∈Lμ2(X)$⁠. In words, these assumptions ensure that trajectories asymptotically settle into an attractor $A⊂X$ (which may be steady, time-periodic, quasiperiodic, or chaotic) and remain on that attractor (i.e., there is no “switching” between attractors).

Infinitesimal perturbations around a given trajectory obey the variational equation

where $v$ belongs to the tangent space of the manifold $X$ at point $x$⁠, denoted by $TxX$⁠, and $L(x;v)≜dF(x;v)$ is the Gâteaux derivative of $F$ evaluated at $x$ along the direction $v$⁠. In principle, the variational equation could be used to track directions of instabilities around trajectories. In practice, however, this is impossible, because any collection of perturbations ${vi}i=1r$ evolved with (4) for a sufficiently long time would see the magnitude of its members grow or decay exponentially fast, and the angle between them rapidly vanish.

To compute a set of meaningful directions (or “modes”) from the variational equation, Babaee and Sapsis³⁷ proposed to enforce orthonormality of the $vi$’s at all times. One way to achieve this is to continuously apply the Gram–Schmidt algorithm to the collection ${vi}i=1r$⁠, starting with $v1$ and moving down. Blanchard and Sapsis¹⁶ showed that the resulting vectors obey

for $i∈{1,…,r}$⁠, where the angle brackets denote a suitable inner product on $TxX$⁠. In the above, we recognize the variational equation (the left-hand side and the first term on the right-hand side), appended with terms enforcing the orthonormality constraint (the last two terms on the right-hand side). We also note that the summation index goes to $i−1$ rather than $r$ so that (5) assumes a lower-triangular form. (This reflects the very structure of the Gram–Schmidt algorithm.) The $ui$’s have been referred to as the “OTD modes” and the collection ${ui}i=1r$ as the “OTD subspace.”³⁷ The terms “subspace” and “modes” are appropriate because the collection ${ui}i=1r$ forms a real vector space, for which the $ui$’s are an orthonormal basis.

A key property of the OTD modes is that they and the $vi$’s span the same subspace. The first OTD mode aligns with the most unstable direction, just like $v1$ does. The second OTD mode is not free to align with the second-most unstable direction because it must remain orthogonal to $u1$⁠. However, the subspace spanned by $u1$ and $u2$ coincides with the two-dimensional subspace that is most rapidly growing (this is also the subspace spanned by $v1$ and $v2$⁠). For hyperbolic fixed points, the OTD subspace aligns with the most unstable eigenspace of the associated linearized operator.³⁷ For time-dependent trajectories, the OTD subspace aligns with the left eigenspace of the Cauchy–Green tensor, which characterizes transient instabilities.³⁹ As a result, the $i$th Lyapunov exponent $λi$ can be recovered from the $i$th OTD mode,

The fact that the OTD modes track directions of instabilities along any trajectory has been utilized on multiple occasions, including in the context of prediction of extreme events in chaotic systems⁴⁰ and design of low-dimensional controllers for stabilization of unsteady flows.^41,42

The OTD system (5) is a set of time-dependent differential equations, which must be solved concurrently with the dynamical system (1). The standard approach for infinite-dimensional systems is to discretize (1) and (5) in space and advance each using a time-stepping scheme. The dimension $d$ of the phase space after discretization may be quite large, with the discretized state potentially having thousands of millions of degrees of freedom. In such cases, computation of the first $r$ OTD modes involves solving $r$ $d$-dimensional differential equations (the OTD equations), plus a $d$-dimensional differential equation for the state itself. For very large $d$⁠, this procedure is computationally onerous, and alternative approaches might be desirable.

For dynamical systems satisfying the assumptions made earlier (autonomy, invertibility, ergodicity, and measure-invariance), the OTD modes asymptotically converge to a set of vectors defined at every point on the attractor.^15,16,43 In other words, in the post-transient regime, $ui$ only depends on the state $x$⁠, but not on the history of the trajectory or its own initial condition $ui(t0)$⁠. Hence, we may cease to view $ui$ as being parametrized by $t$ and instead view it as a graph from the phase space to tangent space,

In this context, the collection ${ui(x)}i=1r$ has been referred to as the “stationary Lyapunov basis” (SLB) at point $x$ (Ref. 15).

The existence of the SLB at every point $x$ of the attractor was established separately by Ershov and Potapov¹⁵ and Goldhirsch et al.⁴³ as a consequence of the Oseledec theorem.⁴⁴ The question of uniqueness and continuity with respect to $x$ was also addressed by Ershov and Potapov.¹⁵ They showed that for a given $x$⁠, more than one SLB may exist, but only one is stable with respect to perturbations of the underlying state. So, the OTD modes $ui(x)$ are uniquely determined by the point $x$ in the phase space. They also showed that if the Lyapunov spectrum is not quasidegenerate (i.e., all Lyapunov exponents are distinct, and there is no “crossing” of Lyapunov exponents under small perturbations), then the graph (7) is continuous in $x$⁠. Uniqueness and continuity are important because they allow for the possibility of representing the graph (7) as a superposition of smooth basis functions, or as a realization of a Gaussian process. Once the graph (7)—or an approximation of it—is known, the Lyapunov exponents can be computed by replacing time-averaging with measure-averaging in (6),

A promising approach is to learn the mapping (6) from data. This requires several ingredients. First, we assume that we have a large collection of “snapshots” ${xn}n=1N$ available for the state. Each $xn$ must belong to the attractor, but not necessarily to the same trajectory, a consequence of the use of measure-averaging. Second, we assume that we have a mechanism to compute the vector field $F(xn)$ and the action of the linearized operator $L$ at $xn$ in any direction $v$⁠. Third, we need to eliminate the dependence of the OTD system (5) on time. This is done by applying the chain rule to the left-hand side of (5), resulting in

where $dui(x;F(x))$ is the Gâteaux derivative of $ui$ evaluated at $x$ along the direction $F(x)$⁠. There is no explicit temporal dependence in (9) so that the variable $x$ should no longer be thought of as a point on a particular time-dependent trajectory, but rather as any point on the attractor. System (9) may also be thought of as a partial differential equation on $X$⁠.

For computational purposes, it is useful to consider the discretized counterpart of (9),

where $x$ and $ui$ belong to $Rd$⁠, $L(x)≜∇xF(x)$ is the Jacobian of the vector field $F:Rd⟶Rd$⁠, and $∇xui$ is the Jacobian of $ui$ with respect to $x$⁠. Although not explicitly shown in (10), the vector $ui$ should be understood as $ui(x)$⁠. We are now in a position to state the learning problem:

In what follows, we solve the learning problem by a deep-learning approach based on neural networks.

We will find it convenient to operate in the “big-data” regime, so we assume that the dataset used in the learning algorithm contains a very large number of snapshots. Before we proceed, we reiterate the fundamental assumptions that are made about the data. The underlying dynamical system from which data are collected should be autonomous, invertible, measure-preserving, and ergodic, and should have a nonquasidegenerate Lyapunov spectrum. As discussed in Sec. II B, these assumptions are key to ensure existence, unicity, and continuity of the SLB in the phase space.

To learn the graphs ${ui}i=1r$ from the collection of snapshots, we employ a neural-network approach. This is appropriate because each $ui$ is a continuous function of $x$⁠. This allows us to leverage the universal approximation theorem,⁴⁵ which states that any function may be approximated by a sufficiently large and deep neural network.

Neural networks do suffer from several shortcomings. They come with heavy computational requirements, are sensitive to the choice of hyperparameters, and do not always deliver on noisy data. However, the pros largely outweigh the cons. First, neural networks are known to be more powerful and more expressive than simpler (neighborhood-based or autoregressive) methods. Second, as discussed in Sec. I, neural networks have been used successfully to solve partial differential equations, which is essentially what the learning problem proposed in Sec. II B amounts to. Third, neural networks are generally trained by some variant of stochastic gradient descent and thus do not suffer from the requirement that the entire dataset be loaded in memory at training time. With an eye on large-scale deep-learning applications, neural networks are, therefore, the de facto correct choice for our learning problem.

In what follows, we refer to the learned OTD modes as the “deep OTD (dOTD) modes.”

We find it natural to assign to each OTD mode its own neural network $ui(x;θi)$⁠, where $θi$ denotes the parameters (weights and biases) of the $i$th network. We use the same fully-connected feed-forward architecture with $L$ hidden layers for all OTD modes [Fig. 1(a)]. (The input and output layers are numbered $0$ and $L+1$⁠, respectively.) The loss function for the $i$th network is specified as

which is nothing more than the residual of the $i$th OTD equation in system (10). We note that any SLB is a global minimizer of $ℓipde$⁠, with $ℓipde$ being trivially zero. So, this choice of loss function gives rise to no estimation error, and the model is only limited by the hypothesis class (i.e., the set of functions within reach of the neural network for a given number of layers and neurons) and the tolerance specified for the optimization algorithm. Those are commonly referred to as “approximation error” and “optimization error,” respectively.⁴⁶ 

Equation (11) shows that the loss function for the $i$th dOTD mode depends on the first $i−1$ dOTD modes. This raises the question of the order in which the dOTD modes ought to be learned. In what follows, we opt for a “sequential” approach [Fig. 1(b)], whereby the parameters $θi$ are optimized sequentially (starting with $θ1$ and moving down), and the outputs of the first $i−1$ networks are fed into the $i$th network as dummy inputs. (By “dummy inputs,” we mean quantities that are fed into the neural network for the sole purpose of computing the loss function.) We find this architecture to be intuitive because it mimics the triangular structure of the Gram–Schmidt algorithm. Our numerical experiments suggest that this approach is more stable (compared to other approaches described below), in the sense that it facilitates convergence of the optimization algorithm. The only issue has to do with error accumulation, arising as a result of using the outcomes of the first $i−1$ optimizations to compute the $i$th dOTD mode. However, this is easily fixed by tightening the tolerance of the optimization algorithm or doing multiple passes of training with decreasing tolerance.

Of course, the “sequential” approach is not the only option. Instead, one could solve the $r$ optimization problems concurrently, whereby the optimization algorithm performs one iteration for each neural network before updating the parameters. Each iteration would need to be done sequentially (i.e., starting with $θ1$⁠, then $θ2$⁠, etc.), and “coupling” between the dOTD modes would be done by passing the first $i−1$ parameter vectors at the current iteration to the $i$th neural network. Alternatively, one could combine the $r$ loss functions (11) and solve for all of the dOTD modes in a single optimization pass using the combined loss function. However, these two approaches appear to cause difficulty for the optimization algorithm, both in terms of execution speed and its ability to converge.

Use of the loss function (11) alone, which is solely based on the residual of the OTD system, might be problematic for two reasons. First, as discussed in Sec. II B, any SLB is a global minimizer of (11), but only one of them is stable (this is the SLB to which the OTD modes converge when computed with a time-marching approach). So, we need a mechanism to ensure that the optimization algorithm converges to the stable SLB and not to any of the unstable ones. Second, there may be other (possibly local) minimizers of (11) besides SLBs, so we also need a mechanism that prevents the optimization routine from getting trapped in an irrelevant minimum.

We begin with the question of how to ensure that the optimization algorithm converges to an SLB, with no consideration yet for whether that SLB is stable or not. We first note that minimization of the loss function (11) does not guarantee orthonormality of the resulting dOTD modes. [For example, the trivial solution $ui=0$ is a global minimizer of (11).] The reason is that the terms responsible for orthonormalizing the OTD modes in the time-dependent problem (5) are no longer sufficient to enforce this constraint. So, orthonormality must be enforced explicitly in the neural network. This is important because our numerical experiments suggest that the SLBs are the only orthonormal minimizers of (11). In other words, explicitly enforcing orthonormality of the dOTD modes helps the optimization algorithm systematically converge to an SLB, rather than some other irrelevant minimum.

Enforcing orthonormality of the dOTD modes can be realized in two ways. One approach is to embed Gram–Schmidt orthonormalization immediately after the last layer of the network so that the dOTD modes are orthonormal by construction. Another approach is to append to the loss function (11) a regularization term,

The first term in the curly brackets enforces normality of $ui$⁠, and the second term enforces orthogonality of $ui$ and $uk$ (⁠$k<i$⁠). The regularization parameters $α$ and $βk$ should be chosen on a case-by-case basis, which offers less flexibility than the Gram–Schmidt approach. We note that the regularization approach and the Gram–Schmidt approach have no consequence for the estimation error because their respective loss functions are exactly zero for any SLB.

In practice, we have found that the Gram–Schmidt approach is more robust than the regularization approach, in the sense that the former requires fewer iterations for the optimization algorithm to converge to an SLB. We note that each iteration in the Gram–Schmidt approach requires computing the gradients of the Gram–Schmidt layer with respect to the network parameters by backpropagation, which is significantly more expensive than computing the gradients of the regularizing terms in (12). However, this is a burden worth carrying, given that we have found cases in which the regularization approach failed to converge (for a given iteration budget) while the Gram–Schmidt approach succeeded, and no cases suggesting otherwise.

Now that we have designed a mechanism ensuring that the optimization algorithm converges to an SLB, we must address the question of how to isolate the stable SLB from all of the unstable ones. We begin with a simple example that provides insight into this issue. Consider a case in which the trajectory of interest is a hyperbolic fixed point, denoted by $xe$⁠. Theorem 2.3 in Babaee and Sapsis³⁷ states that any $r$-dimensional eigenspace of the operator $L(xe)$ is an SLB of $xe$ and that the only stable SLB is that associated with the $r$ most unstable eigenvalues of $L(xe)$⁠. In this work, “stability” was determined by examining the time-dependent problem governing the evolution of a perturbed SLB. In the learning problem, however, we have eliminated any notion of temporality from the OTD system. Hence, for the case of a hyperbolic fixed point, the neural network, in its current manifestation, may converge to any of the $d$-choose-$r$ eigenspaces of $L(xe)$ and not necessarily to the most unstable one. This is problematic because the OTD modes draw their power from their ability to track directions of greatest instabilities. Naturally, we wish our learning algorithm to also have this feature.

To make sure that the learning algorithm returns the SLB associated with directions of strongest instabilities, we use a criterion based on Lyapunov exponents. As discussed in Sec. II B, the $i$th Lyapunov exponent $λi$ can be computed as a measure-average of the Lagrange multiplier $⟨ui(x),L(x)ui(x)⟩$⁠. With a finite-size dataset, however, we can do no better than approximating $λi$ by a finite sum over the data points. If the dataset is composed of multiple long trajectories initialized on the attractor according to some probability distribution (in general, we want the initial conditions to be independent and identically distributed), we have that

where $xn$ is the state of the $n$th trajectory after some long time $T$⁠. The above limiting statement holds only when $T→∞$⁠, but in practice, we merely require that $T$ be large enough so as to ensure convergence of the distribution of initial conditions to an invariant one. Equation (13) also holds when the dataset is composed of uniformly-sampled snapshots collected in the asymptotic regime of a single long trajectory, in which case (13) is equivalent to standard time-averaging. (Note that in either case, the snapshots may be arranged in any arbitrary order.) Equation (13) can be modified to account for the fact that $ui$ is modeled as a neural network. We define

as the “learned” Lyapunov exponent associated with the $i$th dOTD mode. This is the best approximation of $λi$ available, given the constraints related to finiteness of the dataset and representability of the OTD modes with neural networks.

With (14) in hand, we can append to the loss function (11) a regularization term,

that penalizes small Lyapunov exponents. Here, $σ$ is a monotonically increasing function on $R$ that exacerbates differences between the $λ^i$’s. Possible choices for $σ$ include $σ(a)=a$⁠, $a3$⁠, $sinh⁡(a)$⁠, and $−exp⁡(−a)$⁠. Lyapunov regularization guarantees that the SLB to which the optimization algorithm converges is such that the $λ^i$’s come in decreasing order; that is, the $i$th dOTD mode picks up the $i$th-largest Lyapunov exponent. Lyapunov regularization thus ensures that the dOTD modes learned by the neural network coincide with the unique stable solution of the time-dependent OTD system (5).

Adding Lyapunov regularization to the loss function (11) has the effect of introducing an estimation error because the value of $ℓilyap(θi)$ for the optimal $θi$ is generally nonzero, except in very specific cases [for example, if $σ(a)=a$ and $λ^i$ is zero]. No estimation error is a feature worth preserving because it allows us to focus our attention on the approximation error and the optimization error, thereby facilitating design and optimization of the neural network. To this effect, we specify an optimization schedule so that Lyapunov regularization is “switched off” after a certain number of iterations. This allows us to “steer” the optimization algorithm into a favorable direction initially, while ensuring no estimation error for iterations subsequent to relaxation of Lyapunov regularization.

The last ingredient needed to make the method fully data-driven is a mechanism to reconstruct the vector field $F(xn)$ and the action of the Jacobian matrix $L(xn)$ from the collection of snapshots ${xn}n=1N$⁠. We note that if the governing equations of the dynamical system are known, then there is no need for such a mechanism because $F(xn)$ and $L(xn)$ can be evaluated directly from the equations of motion. We also note that if we can only record some observable $f(x)$ of the state, but not the state itself, then we can use delay embedding to reconstruct the attractor and compute the dOTD modes in the embedded space. (This case will not be considered in this work.)

As discussed in Sec. I, discovery and reconstruction of governing equations from data is an active field of research. Any of the state-of-the-art methods could be applied to the present problem, each introducing its own level of complexity. The key issue is that reconstruction of $F(x)$ can be done offline, regardless of the dimensionality of the system. In other words, computation of $F(xn)$ for each $xn$ may be viewed as a preprocessing step and, therefore, does not add to the computational burden related to optimizing the parameters of the neural network. In what follows, we opt for perhaps the simplest of all approaches. We assume that the snapshots are sampled along a single long trajectory with a uniform and sufficiently small sampling time-step $Δts$ so that we may approximate $F(xn)$ as a standard Euler-forward finite difference,

Higher-order finite-difference formulas (e.g., Adams–Bashforth, Adams–Moulton, or backward differentiation formulas) may be used if higher accuracy is desired. Finite differences have the advantage of being extremely cheap to compute, even for high-dimensional systems. Implementation is straightforward, and the requirement that $Δts$ be small is far from drastic.

To compute the action of the Jacobian matrix $L(x)$ from data, we employ the classical algorithm proposed independently by Eckmann et al.¹⁰ and Sano and Sawada.¹¹ First, we scan the dataset to identify the $K$ nearest neighbors of each datapoint $xn$⁠. The nearest neighbors of $xn$ are defined as those points $xk$ of the orbit (past or future) that are contained in a ball of radius $γ$ centered at $xn$⁠,

If $γ$ is sufficiently small, then each vector $vkn=xn−xk$ may be viewed as a perturbation vector from the reference orbit. We, therefore, have

which allows us to compute the action of the Jacobian matrix $L(xn)$ on the vectors ${vkn}k=1K$⁠. Now, the critical step is to note that the vectors ${vkn}k=1K$ belong to the tangent space at point $xn$ and so do the OTD modes. (In fact, the OTD modes form a basis of this space when $r=d$⁠.) So, if we stack up the vectors ${vkn}k=1K$ into a matrix $Vn∈Rd×K$⁠, then the least-square fit

should be a reasonably good approximation for the dOTD modes. Here, $Vn†$ is the pseudoinverse of $Vn$⁠. (We note that the least-square approach allows for $K$ exceeding the dimension of the phase space $d$⁠.) From (19), we can compute the action of the linearized operator on the dOTD modes as

where $ΔVn$ is a $d$-by-$K$ matrix with columns $(vk+1n+1−vkn)/Δts$⁠. Equation (20) requires no information other than the snapshot data and can be used to evaluate the loss function (11).

For low-dimensional systems, it is possible to form the matrix $ΔVnVn†∈Rd×d$ explicitly and store it in memory. This computation can be done offline for every datapoint $xn$⁠, and consequently, the computational burden associated with reconstruction is nil for the neural network. For high-dimensional systems, however, forming and storing $ΔVnVn†$ is not possible, so (20) must be computed online [i.e., every time the loss function (11) is evaluated], which introduces an additional cost. Another aggravating issue is that the complexity of most nearest-neighbor-search algorithms grows exponentially with the dimension of the state, making the above approach intractable for $d$ greater than about 25 (Ref. 47). The curse of dimensionality thus requires that we pursue a different strategy.

For high-dimensional systems, the neural-network approach is intractable for two reasons. First, there is the issue of reconstructing the Jacobian matrix (or, equivalently, its action on the dOTD modes), which was discussed in Sec. III B. Second, to evaluate the term $∇xui(xn;θi)$ appearing in (11), we must compute the gradient of the dOTD modes with respect to the state vector, resulting in a $d$-by-$d$ matrix. For large $d$⁠, this computation is virtually hopeless. So, to make the neural-network approach applicable to high-dimensional systems, we proceed to an order-reduction of the phase space.

If $x$ arises from discretizing a partial differential equation (PDE) defined in a domain $Ω$⁠, then one approach is to randomly select $M$ points in that domain according to some probability distribution. Each snapshot $xn$ then has dimension $M$⁠, with $M$ presumably much smaller than $d$⁠. When $x$ contains nodal values of the state, this approach amounts to randomly excising $M$ entries from each $xn$⁠. (The excised entries need not be the same for all the snapshots.) Random sampling has been used successfully in a number of problems related to deep learning of PDEs,^29,48 largely because it has the advantage of being a “mesh-free” approach. However, applicability of this method to the present problem is limited because the algorithm for attractor reconstruction requires that the sampled points be the same for all snapshots. As a result, we may need a large number of sampled points (with $M$ possibly on the order of $d$⁠) to faithfully capture the dynamics. If $M$ is reasonably small, however, the dOTD modes can be learned at the sampled points and subsequently reconstructed over the entire domain using any standard interpolation algorithm.

To avoid these difficulties, we use an approach based on the Galerkin projection. We assume that any state on the attractor can be represented as a superposition of proper-orthogonal-decomposition (POD) modes,

where $x¯$ is the mean flow, $Φ∈Rd×ns$ contains the first $ns$ POD modes, and $ξ∈Rns$ contains the corresponding POD coefficients. Measure-invariance and ergodicity of the attractor allow us to view $ξ$ as a function of time or as a function of the state, so that we may use $ξ(t)$ and $ξ(x)$ interchangeably. Each POD mode can be computed directly from data by the method of snapshots,¹⁸ and the $i$th POD coefficient can be obtained by projecting $x−x¯$ onto the $i$th POD mode. The number of retained POD modes is determined by examining the cumulative energy of the POD eigenvalues. In general, $ns$ is selected so as to account for at least 95% of the total energy. We also note that implicit in (21) is a one-to-one correspondence between snapshots in the $x$-space (⁠$xn$⁠) and snapshots in the $ξ$-space (⁠$ξn$⁠). This allows us to view any function of $xn$ as a function of $ξn$ and vice versa.

In the POD subspace, the dynamics obeys

so, in principle, we could use the neural-network approach to learn the OTD modes in the $ξ$-subspace. We would simply apply the method proposed in Secs. III A and III B with $ξ$ in place of $x$⁠, $G$ in place of $F$⁠, and $∇ξG$ in place of $L$ and then project the dOTD modes learned in the $ξ$-subspace back to the full space, resulting in $Φui(ξ)$⁠. However, this approach is flawed because by construction, it assumes that the OTD modes live in the same subspace as the state itself (that is, the POD subspace spanned by the columns of $Φ$⁠). In reality, the OTD modes belong to the tangent space at point $x$⁠, whose principal directions have no reason to coincide with that of the state on the attractor. This inconsistency persists regardless of the number of POD modes used in (21) so that resorting to very large $ns$ does not solve the problem.

All is not lost, though, since we can use different projection subspaces for the state and the OTD modes,

where $Ψ(x)∈Rd×nt$ is a reduced orthonormal basis of the tangent space at $x$ and $μi(x)∈Rnt$ contains the basis coefficients. We note that the number of POD modes (⁠$ns$⁠) and columns of $Ψ$ (⁠$nt$⁠) need not be the same, which allows for the possibility of learning the OTD modes in a subspace bigger than the POD subspace. If $Ψ(xn)$ is a good approximation of the tangent space at $xn$⁠, then the $K$ nearest neighbors of $xn$ satisfy

where $ank$ is a matrix of coefficients. To compute $Ψ(xn)$ from data, we solve the minimization problem

which is tantamount to a principal component analysis (PCA) of the set of nearest neighbors ${xk}k=1K$⁠. In other words, the columns of $Ψ(xn)$ are the leading $nt$ POD modes of the $K$ nearest neighbors of $xn$⁠; we refer to them as the “tangent POD (tPOD) modes.” We note that this approach has been used as the basis for a number of manifold learning algorithms that involve reconstruction of the tangent space from data.^49,50

Now that we have two low-dimensional representations available for the state and the OTD modes, the learning problem reduces to finding the collection of graphs,

given a dataset of reduced states ${ξn}n=1N$ and tangent bases ${Ψ(ξn)}n=1N$⁠. To derive an equation for $μi$⁠, we substitute (23a) and (23b) into (10) and arrive at

Three remarks are in order. First, we note that (27) is a system of $nt$-dimensional differential equations, much less expensive to solve the original system of $d$-dimensional differential equations. Second, we note the presence of an additional term on the right-hand side of (27) arising from the dependence of the tPOD modes on the state $x$⁠. Third, we have defined the reduced operator

which is the projection of the high-dimensional operator $L(x)$ onto the reduced basis $Ψ(x)$⁠. Equation (28) provides additional insight as to why the “naive” approach described two paragraphs earlier was not a good one. This approach was equivalent to having $Ψ(x)=Φ$ for every point $x$ on the attractor, leading to $LΨ=ΦTL(x)Φ$⁠. However, this is a poor approximation for $L(x)$ because $Φ$ is a reduced basis of the phase space but not the tangent space.

To compute the reduced operator $LΨ(xn)$ from data, we use an approach similar to that proposed in Sec. III B. We first recall that (18) provides a mechanism to compute the action of $L(xn)$ on a collection of perturbation vectors ${vkn}k=1K$ at point $xn$⁠. However, the columns of $Ψ(xn)$ are linear combinations of these perturbation vectors, as per the POD construction. Therefore, we may write $Ψ(xn)=Vnκn$⁠, where $κn∈RK×nt$ is a matrix of coefficients. This leads to

The operator $LΨ(xn)$ has the advantage of being low-dimensional, so it can be computed offline and stored in memory, along with the POD-reduced vector field $G(ξn)$⁠.

The last term on the right-hand side of (27), however, is problematic because it involves the temporal derivative of the reduced tangent basis $Ψ$⁠. We could use the chain rule and write $Ψ˙=∇ξΨG(ξ)$⁠, but the gradient $∇ξΨ$ is expensive to compute. Another option is to use a finite-difference formula in the spirit of (16),

If $Δts$ is sufficiently small, then the tPOD modes smoothly deform from $xn$ to $xn+1$⁠, so (30) is a good approximation. [We note in passing that continuity of $ui(x)$ requires continuity of both $Ψ(x)$ and $μi(x)$⁠.] With this approach, the term $Ψ⊺Ψ˙$ can be computed offline and passed to the neural network as a dummy input.

Finally, we address the question of whether the “local” tangent bases $Ψ(x)$ could be combined into a larger “global” subspace that does not depend on $x$⁠. Mathematically, this is equivalent to seeking a reduced basis of the tangent bundle,

where $⨆$ is the disjoint union operator. Such a construct would have the benefit of eliminating the last term on the right-hand side of (27). Also, it would be esthetically more appealing because the OTD modes would be learned in a common subspace, regardless of the point at which they are computed. To construct a reduced basis of $TX$⁠, one can perform POD on the set of tangent bases ${Ψ(xn)}n=1N$⁠. This results in a set of “bundle modes,” denoted by $Π∈Rd×nb$⁠, whose span is the best $nb$-dimensional approximation of the tangent bundle $TX$⁠. The expectation is that the number of bundle modes $nb$⁠, although generally greater than the number of local tPOD modes $nt$⁠, will still be much smaller than $d$⁠. The dOTD modes may then be sought as vectors in the bundle space; that is, $ui(x)=Πρi(x)$⁠, where the coefficients $ρi∈Rnb$ satisfy

and $LΠ$ is defined as

For $nb$ not too large, the reduced operator $LΠ$ may be computed offline and stored in memory.

Applicability of the tangent-bundle method is limited to cases in which the governing equations of the dynamical system are available. This is because $LΠ$ cannot be computed solely from snapshot data, except when the linearized operator is known $L(x)$ explicitly. To see this, we first recognize that by construction, each bundle mode is a linear combination of the columns of $Υ=V1,…,Vn$⁠. So, evaluation of (33) from data is conditioned on our ability to compute $L(xn)Υ$ or, equivalently, each member of the set ${L(xn)Vm}m=1n$⁠. However, there is no mechanism to compute $L(xn)Vm$ from data when $m≠n$⁠. [Equation (18) holds only for $m=n$⁠.] Thus, evaluation of (33) requires explicit knowledge of the linearized operator and, in turn, the governing equations. Because of this restriction, this approach will not be pursued any further in this work.

We conclude this section with a few words on implementation. The neural network is built from scratch in Python. We use Autograd⁵¹ for automatic differentiation. The activation function between hidden layers is the hyperbolic tangent, although other choices (e.g., sigmoid function or swish function) are possible.⁵² The activation function for the output layer is a linear function. As discussed in Sec. III A 2, the last layer of the neural network is followed by a Gram–Schmidt layer. The weights are initialized according to Xavier initialization.⁵³ For Lyapunov regularization, we use $σ=sinh$⁠.

We use Adam⁵⁴ to solve the optimization problem. In its current manifestation, the code supports mini-batching, although caution must be exercised when selecting the batch size in order not to deteriorate the accuracy of (14). (The results presented in Sec. IV do not use mini-batching.) We have found no need for specifying learning rate schedules in the optimizer. Different stopping criteria may be used, depending on how the vector field and linearized operator are computed. If $F(xn)$ and $L(xn)$ are evaluated from the governing equations, and the hypothesis class is reasonably large, then optimization may be terminated when the loss function (11) is smaller than a user-specified tolerance. (As a result, the overall error is dominated by the approximation error.) If $F(xn)$ and $L(xn)$ are reconstructed from data, then it is preferable to terminate optimization after a user-specified number of iterations, in recognition of the fact that the reconstruction process is approximate and, therefore, introduces an estimation error.

To generate the dataset, we consider a single long trajectory of the dynamical system, rather than multiple shorter trajectories with distinct initial conditions. The dataset is constructed by collecting equally-spaced snapshots (with sampling time $Δts$⁠) along that long trajectory. (For the results presented in Sec. IV, $Δts$ coincides with the time-step size $Δt$ used to advance the dynamical system.) For each snapshot, we reconstruct the vector field and linearized operator using the KNeighborsClassifier implemented in scikit-learn. We note that the nearest-neighbor search is carried out over the entire collection of snapshots, including those collected in the transient regime, so as to improve accuracy of the reconstruction algorithm.¹⁰ For high-dimensional systems, the nearest-neighbor search is conducted in the POD subspace to alleviate computational cost. (By this, we mean that we merely extract time stamps for the nearest neighbors in the $ξ$-space and then use the corresponding snapshots in the $x$-space in the reconstruction algorithm.) In the examples presented in Sec. IV, the vector field and linearized operator are formed explicitly, stored in memory, and passed to the neural network as dummy inputs. This is possible either because the dynamical system is low-dimensional or because we first proceeded to a reduction of the dynamics using the Galerkin approach described in Sec. III C.

Upon reconstruction of $F(xn)$ and $L(xn)$⁠, we discard from the dataset those snapshots that were collected in the transient regime of the long trajectory. As discussed in Sec. III A 1, this step is critical to ensure validity of the measure-averaging operation. The resulting “truncated” dataset is then split into training, validation, and test sets. The training set is used to optimize the weights of the neural network and the test set to evaluate its predictive capability. The validation set is used to tune the hyperparameters of the neural network, as described in Ref. 55. (For a given dynamical system, we consider a range of random combinations of hyperparameters and select the architecture for which the validation error is the smallest.) The validation set is also used to tune the learning rate for the Adam optimizer, the schedule for Lyapunov regularization, and the number of nearest neighbors used in the reconstruction step.

The code is available on GitHub.⁵⁶ 

To evaluate the accuracy of the learning algorithm, we consider two performance metrics. The first metric is the PDE loss $ℓipde$ defined in (11), which quantifies the extent to which the neural network satisfies the OTD equations (10). We use this metric as our primary guide in our search for appropriate hyperparameters because a small PDE loss on training and validation data is a key prerequisite for generalizability on test data.

However, as discussed in Sec. III A 1, the PDE loss does not discriminate between stable and unstable SLBs. So, we supplement it with

where $uideep$ denotes the $i$th dOTD mode and $uinum$ denotes the $i$th OTD mode computed by direct numerical integration of (5). The distance $di$ takes values between 0 and 1, with the former indicating that $uideep$ and $uinum$ coincide and the latter indicating that $uideep$ and $uinum$ are orthogonal. We could also use the distance between the subspaces ${uideep}i=1r$ and ${uinum}i=1r$⁠, but that measure, unlike (34), assigns the same score to all the SLBs, regardless of stability (when $r=d$⁠).

If the OTD modes are learned in a reduced subspace, then it is useful to compute

where $Θ$ is a placeholder for $Ψ(x)$ or $Π$⁠, depending on whether the “local” or “bundle” approach is used. The above is equivalent to computing the distance function in the reduced subspace in which the modes are learned. This quantity is important because dimensionality reduction of the tangent space (or tangent bundle) introduces an additional error over which the neural network has no control. Thus, situations may arise in which $diΘ$ is small, but $di$ is large. Should this occur, a quick fix is to increase the dimension of the reduced subspace $Θ$⁠.

For each example considered below, we report training, validation, and testing error as measured by the two metrics $ℓipde$ and $di$⁠. (For the latter, we report its measure-average, $d¯i$⁠, over the dataset considered.) These numbers have been averaged over ten learning experiments. Each learning experiment corresponds to a different (random) initialization of the weights $θi$⁠. The hyperparameters, training set, validation set, and test set are kept the same from one learning experiment to the next. The results were found to be robust to small changes in hyperparameters and dataset selection.

We begin with a three-dimensional nonlinear system that was proposed by Noack et al.⁵⁷ as a testbed for investigating Hopf bifurcations in laminar bluff-body flows. We choose this system because it provides a good illustration for many of the comments made in Secs. II and III. The governing equations are given by

with $μ$ being a positive constant. The system admits a linearly unstable fixed point (⁠$x=y=z=0$⁠) and a stable periodic solution,

which defines a limit cycle of radius $μ$ in the $z=μ$ plane. As noted by Noack et al.,⁵⁷ the limit cycle is asymptotically and globally stable.

For $μ≤1/8$⁠, it is possible to derive analytical expressions for the OTD modes on the limit cycle,

where

These are the modes to which solutions of the OTD equations converge when computed with a time-stepping approach. [To the best of our knowledge, this is the first time that exact (nonasymptotic) expressions have been reported for unsteady OTD dynamics.] Each $ui$ may be expressed as a function of the state $x=xyzT$ as follows:

The ordered set ${u1(x),u2(x),u3(x)}$ is the unique stable SLB at point $x$ on the limit cycle, although other unstable SLBs exist. If we define

then any of the ordered triples ${u1,u2±,u3±}$⁠, ${u2±,u1,u3±}$⁠, and ${u2±,u3±,u1}$ are solutions to the OTD equations, but only ${u1,u2+,u3+}$ is stable.

We consider the case $μ=0.1$ and use the neural-network approach to learn the graphs ${x⟼ui(x)}i=13$ from data. We begin by generating a long trajectory initiated close to the limit cycle. This trajectory is computed by a third-order Adams–Bashforth method with time-step size $Δt=0.01$ for a total duration of $T=50$ time units. We record snapshots at every time step, resulting in an initial dataset comprising 5000 points. For each datapoint, we reconstruct the vector field and linearized operator using its seven nearest neighbors. We then discard the first 1000 datapoints corresponding to the transient regime. The resulting dataset is comprised of 4000 snapshots recorded in the interval $10≤t≤50$⁠, spanning about six periods in the asymptotic regime of the long trajectory. We then divide up this dataset into three. The training set is comprised of 10 equally-spaced points over the period $20≤t≤20+2π$⁠; the validation test of 629 equally-spaced points over the period $10≤t≤10+2π$⁠, and the training set of the remaining points. (We use the same dataset splitting to train the three neural networks.) The neural network is composed of one hidden layer with 40 neurons, and the learning rate for the Adam optimizer is set to 0.04. Lyapunov regularization is active for the first 1000 iterations and turned off for the remainder of the optimization. The dOTD modes are learned sequentially, with the maximum number of iterations specified as 3000. For these parameters, the training, validation, and test errors are shown in Table I.

Figure 2(a) shows the distance function $di$ computed on part of the test set. [Time series for the dOTD modes $uideep$ and the numerically integrated OTD modes $uinum$ are shown in Figs. S1(a)–S1(c) of the supplementary material.] The key points are that (a) the neural network is able to learn the graphs ${x⟼ui(x)}i=13$ from the limited number of training points supplied to it and (b) error accumulation is inconsequential because we were careful not to terminate the optimization prematurely. We note that $d3$ closely follows $d2$⁠, a consequence of the fact that $u3deep$ is completely determined by $u1deep$ and $u2deep$⁠. We also note that Lyapunov regularization is instrumental in the optimizer converging to the stable SLB (39). The Lyapunov exponents learned by the algorithm (⁠$λ^1=−0.005$⁠, $λ^2=−0.272$⁠, and $λ^3=−0.655$⁠) are reasonably close to the analytical values (⁠$λ1=0$⁠, $λ2=−0.276$⁠, and $λ3=−0.724$⁠).

As discussed in Sec. I, the novelty of our approach is that the neural network provides, locally for each $x$⁠, not only the directions of strongest instabilities, but also the degree of instability along each of these directions. This is made visually clear in Fig. 2(b), which shows all the points from the test set in the three-dimensional phase space. Each point is colored according to its local leading Lyapunov exponent $δ1(x)=⟨u1(x),L(x)u1(x)⟩$⁠, as learned by the neural network. The neural network correctly learns that $δ1(x)=0$ for all $x$ on the limit cycle. Figure 2(b) also shows, for two of the test points, the three dOTD modes (i.e., the directions of strongest instabilities) learned by the algorithm. Consistent with (39), the first mode is tangent to the trajectory. Figure 2(b) thus reinforces the claim made in Sec. I that our approach provides a “cartography” of instabilities in the phase space.

Next, we consider a modified version of the classical Charney–DeVore model, which describes atmospheric circulations at midlatitudes. We consider a six-dimensional truncation of the original system, which models a barotropic flow in a plane channel with topography.⁵⁸ The governing equations are given by

with parameters

where $m=1$ or 2. The parameters $αm$ and $δm$ account for zonal advection in the $z1$ and $z4$ directions, respectively; $βm$ for the so-called $β$ effects; $γm$ and $γm∗$ for topographic interactions; $C$ for Ekman damping; and $z1∗$ and $z4∗$ for zonal forcing in the $z1$ and $z4$ directions, respectively. We set $z1∗=0.95$⁠, $z4∗=−0.76095$⁠, $C=0.1$⁠, $β=1.25$⁠, $γ=0.2$⁠, and $b=0.5$⁠. These values of the parameters give rise to significant transitions between regimes of the “zonal” and “blocked” flow, resulting from the nonlinear interaction between barotropic and topographic instabilities.⁵⁸ The extreme episodes of blocked flow are of main interest to us because they are the manifestation of transient instabilities. Thus, it is in those intervals where we attempt to learn the OTD modes from data.

We use a third-order Adams–Bashforth scheme with time-step size $Δt=0.05$ to generate a long trajectory spanning 4000 time units. We use zero initial conditions, except for $z1(0)=1.14$ and $z4(0)=−0.91$⁠. Snapshots are recorded at every time step, resulting in an initial dataset with 80 000 points. For each datapoint, we reconstruct the vector field and linearized operator using its 60 nearest neighbors and then discard the first 10 000 data points corresponding to the transient regime (⁠$0≤t≤500$⁠). The resulting dataset is then divided up into training, validation, and test sets. For the training set, we consider the interval $1075≤t≤1165$ during which the trajectory passes through a regime of blocked flow [see Figs. S2(a)–S2(f) in the supplementary material] and select 50 uniformly-spaced training points in this interval. The remaining 1750 points in this interval form the validation set. The test set comprises all data points in the original dataset, except for those used for training and validation. The test data, therefore, contains multiple episodes of the blocked flow [Fig. 3(a)]. We only attempt to learn the first OTD mode $u1$⁠. The neural network has two hidden layers, each with 128 neurons. The learning rate for the Adam algorithm is 0.001. Lyapunov regularization is used for the first 2000 iterations and switched off thereafter. Optimization is terminated at 5000 iterations. For these parameters, the training, validation, and test errors are shown in Table II.