Harmonic decompositions of multivariate time series are considered for which we adopt an integral operator approach with periodic semigroup kernels. Spectral decomposition theorems are derived that cover the important cases of two-time statistics drawn from a mixing invariant measure. The corresponding eigenvalues can be grouped per Fourier frequency and are actually given, at each frequency, as the singular values of a cross-spectral matrix depending on the data. These eigenvalues obey, furthermore, a variational principle that allows us to define naturally a multidimensional power spectrum. The eigenmodes, as far as they are concerned, exhibit a data-adaptive character manifested in their phase which allows us in turn to define a multidimensional phase spectrum. The resulting data-adaptive harmonic (DAH) modes allow for reducing the data-driven modeling effort to elemental models stacked per frequency, only coupled at different frequencies by the same noise realization. In particular, the DAH decomposition extracts time-dependent coefficients stacked by Fourier frequency which can be efficiently modeled—provided the decay of temporal correlations is sufficiently well-resolved—within a class of multilayer stochastic models (MSMs) tailored here on stochastic Stuart-Landau oscillators. Applications to the Lorenz 96 model and to a stochastic heat equation driven by a space-time white noise are considered. In both cases, the DAH decomposition allows for an extraction of spatio-temporal modes revealing key features of the dynamics in the embedded phase space. The multilayer Stuart-Landau models (MSLMs) are shown to successfully model the typical patterns of the corresponding time-evolving fields, as well as their statistics of occurrence.

This article introduces a new approach for the stochastic inverse modeling of multivariate time series issued from a mixing dynamical system. For such time series, a new framework, relying on semigroup and integral operator theory, is first presented, and data-adaptive spectral theorems are derived. This framework allows us in turn to extract—for a broad class of time-evolving datasets—reduction coordinates that can be efficiently modeled within a family of stochastic differential equations (SDEs) built from a fixed set of predictor functions. The data-driven modeling effort is thus reduced to elemental stochastic models identified as multilayer Stuart-Landau (SL) oscillators.

The success of machine learning algorithms generally depends on data representation, and it is often hypothesized that this is because different representations can entangle and hide more or less the different explanatory factors of variations behind the data.1 To be reasonably addressed, the problem of finding relevant data-adaptive basis functions for a successful inverse modeling must be considered within a specific domain of applications.2 This article is concerned with addressing this problem for time-evolving datasets produced by a dynamical system, either deterministic or stochastic.

Among the numerous techniques useful for the data representation or decomposition of time-evolving datasets, we may distinguish several seminal techniques like those based on variance's decomposition such as principal component analysis (PCA)3 and its probabilistic formulation4 or nonlinear extensions;5–7 techniques using eigenfunctions of Markov matrices reflecting the (local) geometry and density of the data;8–11 or other recent approaches exploiting the Koopman operator theory.12–15 

These methods often tied to dimensionality reduction, extract data-adaptive modes that come with reduction coordinates, which for datasets evolving in time, correspond to time series known as principal components (PCs) and the like. The effective derivation of reduced evolution equations to characterize the long-term dynamics of the underlying dynamical system16 based on these coordinates is then a central issue.

For some systems, the form of the master equations, prior experience, or some underlying physical intuition may help to determine good reduced systems. In some special cases, their exact form can even be rigorously found by adopting, for instance, the Mori-Zwanzig (MZ) projection approach17–22 or techniques rooted in the approximation theory of (stochastic) local invariant manifolds; see Refs. 23, 24, and references therein. In more complex cases where a rigorous derivation of the reduced dynamics is mathematically intractable, or when the master equations are even not known, a data-driven inverse modeling approach is often adopted.

Many methods have been proposed in the literature to address this task, such as those using nonlinear autoregression moving average with exogenous inputs (NARMAX), artificial neural networks (ANNs), stochastic linear inverse models (LIMs), empirical model reduction (EMR), and many others; see Refs. 25–33 and references therein. In a data-driven context, the MZ formalism also provides a theoretical guidance for the determination of nonlinear ingredients and other dependences on the past; see Ref. 34. Nevertheless, whether the approach retained uses an a priori set of predictor functions or allow for dictionaries of such functions,35 the modeler is often left, either with an explanatory deficit or with a selection problem of the appropriate (class of) predictors for a reliable emulation of the dynamics. In the latter case, a certain level of familiarity with the dataset is typically required to envision the right form of the predictors.

This article proposes to overcome such hurdles for inverse modeling, by introducing a data-adaptive spectral representation of the original data, aimed at providing—for a broad class of time-evolving datasets—reduction coordinates that can be efficiently modeled within a family of stochastic differential equations (SDEs) constituted from a fixed set of predictor functions; see Secs. VIII and IX. Roughly speaking, the proposed modeling approach operates well for finite multivariate time series for which decays of correlations are sufficiently well-resolved for the amount of available data. It thus covers a broad class of datasets issued from dynamical systems that are mixing.

The approach adopted in this article for the data representation of time-evolving datasets relies on the spectral analysis of integral operators L for which the kernels are built from periodic semigroups applied to the data's correlations. The resulting class of operators constitutes one of the two main contributions of this article, and a substantial portion of the latter is devoted to the analysis of their spectral properties; see Secs. II–V.

The other main contribution is concerned with the aforementioned class of SDEs, aimed at the modeling of the reduction coordinates obtained as projection of the original dataset onto the eigenfunctions of L . These SDEs fall into the class of networks of linearly coupled Stuart-Landau oscillators,36 that may be embedded within the class of multilayer stochastic models (MSMs),34 allowing thus for the inclusion of memory effects in their formulation; see Sec. VIII.

Note that readers who might be more motivated by the applications may skip Secs. II–V, at a first reading, and focus on Secs. VI–IX; the essential material from the theoretical results of Secs. II–V being recalled as needed in Secs. VI–IX.

Another originality of the framework lies in its harmonic flavor. The framework allows indeed for extracting modes naturally ranked per Fourier frequency, with reduction coordinates that are, for mixing dynamical systems, typically narrowband and modulated in amplitude; see Sec. VIII A. The corresponding reduced SDEs are also organized per Fourier frequency and the modeling of a specific frequency band is thus naturally made available and facilitated by the resulting framework.

The article is organized as follows. In Sec. II, we first recall basic elements of the spectral theory of periodic semigroups on a Banach space. In Sec. III, this spectral theory is applied to analyze the spectral properties of a class of integral operators L ϕ with periodic kernels built from periodic semigroups acting on a given (one-dimensional) data-function ϕ; see Theorem III.1. The important case for applications of kernels built from left circular translation groups is analyzed in Sec. IV. Here, one recovers the Fourier spectrum (power and phase) from the spectral elements of L ϕ for which the eigenfunctions are phase-shifted sinusoids and their phase relates to the phase spectrum of ϕ; see Theorem IV.1. This representation of the Fourier spectrum via the operator L ϕ allows us to propose in turn an innovative generalization for multidimensional signals in Sec. V. There we show, for multidimensional data, that the eigenvalues of an operator matrix generalization, L , can be grouped per Fourier frequency, and are actually given, at each frequency, as the singular values of a cross-spectral matrix depending on the data; see Theorem V.1. These eigenvalues obey furthermore a variational principle that allows us to define naturally a multidimensional power spectrum; see Remark V.1–(ii). The eigenmodes, as far as they are concerned, exhibit a data-adaptive character manifested in their phase which allows us in turn to define a multidimensional phase spectrum; see (80).

In Sec. VI, we then consider the case for which ϕ is—or a collection of such ϕ's are—obtained from lagged correlations estimated from time series issued from a mixing dynamical system. After discussing a natural periodization to apply the theory (see Sec. VI B), we introduce in Secs. VI C and VI D, the Hankel matrices and block-Hankel matrices resulting, respectively, from the discretization of L ϕ in the one-dimensional case, and of L in the multidimensional case. First applications of the resulting data-adaptive harmonic (DAH) methodology are then discussed for analytical examples of modulated and traveling waves; see Sec. VII.

In Sec. VIII, we introduce the class of MSMs tailored on stochastic Stuart-Landau oscillators, aimed at the modeling of the aforementioned reduction coordinates corresponding here to the projection of the dataset onto the DAH modes (DAHM); see Secs. VIII A and VIII B. These elemental models—the multilayer Stuart-Landau models (MSLMs)—are stacked per frequency and are only coupled at different frequencies by the same noise realization; see Sec. VIII C. In Sec. IX, we finally show the flexibility of the DAH-MSLM modeling approach, by its ability to provide skilled stochastic inverse models for data issued either from the nonlinear chaotic Lorenz 96 model or a stochastic heat equation driven by a space-time white noise. Section X discusses the directions for future research.

Let us first recall the definition of a strongly continuous semigroup also known as C0-semigroup; see Ref. 37 (p. 36). For an introduction to semigroup theory, we refer to Refs. 38 and 39. More advanced treatments and applications of semigroup theory can be found in, e.g., Refs. 37 and 40–42.

Definition II.1. A family ( T ( t ) ) t 0 of bounded linear operators on a Banach space W is called a strongly continuous (one-parameter) semigroup if

T ( t + s ) = T ( t ) T ( s ) , for all t , s 0 , T ( 0 ) = Id ,
(1)

and for every φ in W , the orbit

+ W t T ( t ) φ
(2)

is continuous.

The generator A of a strongly continuous semigroup ( T ( t ) ) t 0 is then defined as the operator

A : D ( A ) W W ,

such that

A φ = lim h 0 + 1 h ( T ( h ) φ φ )
(3)

for every φ in the domain

D ( A ) : = { φ W | lim h 0 + 1 h ( T ( h ) φ φ ) exists } .
(4)

Our approach relies on the following important characterization of spectral properties of periodic semigroups, e.g., Ref. 37 (Chap. IV, Sec. 2.25). Recall that a strongly continuous semigroup ( T ( t ) ) t 0 is periodic if there exists t 0 > 0 such that T ( t 0 ) = Id . The period τ of the semigroup is then

τ = inf { t 0 > 0 : T ( t 0 ) = Id } .

Periodic semigroups are always groups with inverses T ( t ) 1 = T ( n τ t ) , for 0 t n τ , and the spectrum of their generators lies always on the imaginary axis, more precisely we have

Theorem II.1. Ref. 37 (Chap. IV, Lemma 2.25) Let T(t) be a strongly continuous τ-periodic semigroup on a Banach space W , and A its generator, then

σ ( A ) 2 π i τ · ,

and the resolvent R(μ,A) ( μ ) of A is given by

R ( μ , A ) = ( 1 e μ τ ) 1 0 τ e μ s T ( s ) d s , μ 2 π i τ · .

The above representation of the resolvent

R ( μ , A ) : = ( A μ Id ) 1

shows that the resolvent of the generator of a τ-periodic semigroup is a meromorphic function having only poles of maximal order, with poles μn given by

μ n = 2 i π n τ .
(5)

The residue of R(μ,A) at the pole μn is given by the bounded linear operator of W [Ref. 37, (p. 267)]

P n : = lim μ μ n ( μ μ n ) R ( μ , A ) = 1 τ 0 τ e μ n s T ( s ) d s .
(6)

Hereafter, we denote by L ( W ) the space of bounded linear operators of the Banach space W .

From Ref. 37 (Chap. IV, Sec. 1.17), it follows furthermore that Pn coincides with Πn, the Riesz projector associated with the pole μn, namely, with

Π n : = 1 2 i π Γ c R ( μ , A ) d μ ,
(7)

where Γc denotes a rectifiable closed curve in the complex plane that separates μn from the rest of the poles. Furthermore, the rank of Πn, rg Π n , satisfies

rg Π n = ker ( A μ n Id ) ,

i.e., the algebraic multiplicity is equal to the geometric multiplicity. In particular, this implies that the spectrum σ ( A ) of A consists of eigenvalues only (i.e., no continuous spectrum) and that

A P n = μ n P n .
(8)

Recall that for any strongly continuous T(t) with generator A, for every λ in , t > 0, and φ in D(A), the following identity holds [Ref. 37 (Chap. II, Lem. 1.9)]:

e λ t T ( t ) φ φ = 0 t e λ s T ( s ) ( A λ Id ) φ d s ,
(9)

which here due to (8) leads to

T ( t ) P n = e μ n t P n , t 0.
(10)

Moreover,

P m P n φ = 1 τ 0 τ e μ m s T ( s ) d s · P n φ = 1 τ 0 τ e μ m s T ( s ) P n d s · φ = 1 τ 0 τ e ( μ n μ m ) s d s · P n φ = 0 , for n m .
(11)

With these properties at hand, it is then not difficult to infer the following useful decomposition result of any strongly continuous τ-periodic semigroup and of its generator; see Ref. 37 (Chap. IV, Thm. 2.27).

Theorem II.2. Let T(t) be a strongly continuous τ-periodic semigroup on a Banach space W , and A its generator, then

T ( t ) φ = m e μ m t P m φ , φ D ( A )
(12)

and

A φ = m μ m P m φ , φ D ( A 2 ) .
(13)

Let us introduce the one-dimensional torus T : = / ( τ ) which we define by identifying points in that differ by n τ for each n in . We denote by L 2 ( T ) the space of complex-valued measurable functions on that are τ-periodic, and square-integrable with respect to an arbitrary reference interval of length τ, I = [ α , α + τ ] for some α in . We endow L 2 ( T ) with the natural inner product

f , g = 1 τ T f ( r ) g ( r ) ¯ d r , f , g L 2 ( T ) .
(14)

Here, the integral over T is the integral with respect to r taken over any interval of length τ. In what follows, we indeed often identify T with its interval representation I, making precise the distinction when necessary.

Given a τ-periodic C0-semigroup T(t) acting on a function space W intersecting the set of τ-periodic functions, and a τ-periodic function ϕ in W , our goal is to analyze the spectral property of the following operator:

L ϕ [ Ψ ] ( r ) : = 1 τ T [ T ( r ) ϕ ] ( s ) Ψ ( s ) d s , Ψ L 2 ( T ) .
(15)

In practice, the space W may differ from L 2 ( T ) (by requiring, for instance, other regularities on ϕ), but to simplify the presentation, we restrict ourselves to the case

W = L 2 ( T ) .

The following lemma results directly from the theory of integral operators:

Lemma III.1. Assume ϕ lies in L 2 ( T ) and let T(t) be aτ-periodic C0-semigroup on L 2 ( T ) . If

  • (A1) ( r , s ) [ T ( r ) ϕ ] ( s ) belongs to L 2 ( T × T ) ,

    then the operator L ϕ defined by(15), maps L 2 ( T ) into L 2 ( T ) and iscompact.If furthermore

  • (A2) [ T ( r ) ϕ ] ( s ) ¯ = [ T ( s ) ϕ ] ( r ) ,

    then the operator L ϕ is self-adjoint.

Proof. The proof is standard and consists of noting that under the integrability condition (A1), the operator L ϕ is a Hilbert-Schmidt operator on L 2 ( T ) and thus compact; see, e.g., Ref. 43 (Thm. 6.12) or Ref. 44 (Chap. 16, Theorem 3). Finally, condition (A2) allows us to write

L ϕ [ Ψ ] , Ψ = T T [ T ( r ) ϕ ] ( s ) Ψ ( s ) Ψ ( r ) ¯ d s d r = T T [ T ( s ) ϕ ] ( r ) Ψ ( r ) Ψ ( s ) ¯ d s d r ¯ = L ϕ [ Ψ ] , Ψ ¯ ,
(16)

which shows that L ϕ [ Ψ ] , Ψ lies in for any Ψ in L 2 ( T ) , and thus L ϕ is self-adjoint. ◼

As a consequence of the spectral theorem for self-adjoint compact operators, we have the following theorem, in which either J = { 1 , , N } with N = dim ( rg ( L ϕ ) ) if L ϕ has finite rank, or J = otherwise.

Theorem III.1. Let L ϕ be defined in (15), with T(t) and ϕ satisfying the assumptions of Lemma III.1 including(A1)and(A2). Then, 0 is an eigenvalue of L ϕ and there exist:

  • countably many nonzero real numbers { λ n } n J either finitely many, or such that λ n 0 if infinitely many, and

  • an orthonormal basis { E n } n J of cl ( rg ( L ϕ ) ) , the topological closure of the range of L ϕ in L 2 ( T ) ,

for which

L ϕ [ Ψ ] = n J λ n Ψ , E n E n .
(17)

Each λn is an eigenvalue and each En is an eigenfunction.

Due to the convergence to zero, we usually order the eigenvalues λn in decreasing order according to the absolute value

| λ 1 | | λ 2 | .
(18)

The multiplicity of a given eigenvalue λ is the number of times it is repeated in the list of eigenvalues given above, or, equivalently, it is the dimension of the λ-eigenspace ker ( L ϕ λ Id ) . The multiplicity of any given eigenvalue is finite.

So far, we have applied classical results from the theory of integral operators to deal with the spectral properties of L ϕ . We can derive furthermore a key identity satisfied by the eigenfunctions of L ϕ [see (21)], by exploiting the spectral theory of periodic semigroups recalled in Sec. II and requiring more regularity on ϕ.

For that purpose, we focus here on the spectral properties of L ϕ in L 2 ( T ) , i.e., the subspace of L 2 ( T ) constituted by real-valued functions. We have then the following result concerning the spectral elements of L ϕ in L 2 ( T ) , for which we make use of the following version of the Fourier transform Ψ ̂ of a function Ψ in L 2 ( T ) :

Ψ ̂ ( f ) : = T Ψ ( s ) e 2 i π f s d s , f .
(19)

Theorem III.2. Assume ϕ lies in D(A), where A denotes the generator of a τ-periodic C0-semigroup T(t) on L 2 ( T ) . Assume that(A1)and(A2)hold. Assume

  • (A3)there exists an orthogonal basis { e n } of L 2 ( T ) for which Pn defined in (6) is rank-one, i.e.,
    P n f = f , e n e n .
    (20)

Let us denote by En an eigenfunction of L ϕ in L 2 ( T ) associated with an eigenvalue λn, then for any frequency f k = k / τ , we have

λ n E ̂ n ( f k ) = ϕ , e k e k , E n .
(21)

Proof. Let us consider an eigenfunction En associated with an eigenvalue λn of the operator L ϕ given in (15). Then, the identity (12) of Theorem II.2 gives

τ λ n e μ k r E n ( r ) = T e μ k r T ( r ) ϕ ( s ) E n ( s ) d s , = m e ( μ m μ k ) r T ( P m ϕ ) ( s ) E n ( s ) d s ,
(22)

which due to (5) leads, by integration in the r-variable, to

τ λ n T e μ k r E n ( r ) d r = T ( P k ϕ ) ( s ) E n ( s ) d s .
(23)

The latter identity gives (21) by applying (19) with f = k / τ , combined with Assumption (A3) and the definition of the inner product · , · in (14), recalling that En is real-valued. ◼

Remark III.1. Under the conditions of Theorem III.1, one has furthermore

λ 1 = sup Ψ W , Ψ 0 L ϕ Ψ , Ψ | Ψ | L 2 2 ,
(24)

with W = L 2 ( T ) , e.g., Ref. 43 (Prop. 6.9).

We actually have a variational characterization of any eigenvalue of the self-adjoint compact operator L ϕ by relying on the Courant-Fischer min-max principle; see Ref. 43 (problem 37) and Ref. 45. However, since the operator L ϕ is not positive, the Courant-Fischer min-max principle needs to be amended. The positive eigenvalues of L ϕ listed in decreasing order and denoted by 0 λ 2 + λ 1 + satisfy thus the relationships

λ k + 1 + = min V codim ( V ) = k max Ψ V Ψ 0 L ϕ Ψ , Ψ | Ψ | L 2 2 , k 0.
(25)

That is, the max is taken over a linear subspace V W of co-dimension k, and the min is taken over all such subspaces. Moreover, the minimum is attained on the subspace

V = span { E 1 + , , E k + } ,
(26)

where E j + denotes the eigenfunction associated with a positive eigenvalue λ j + .

Similarly, for 0 λ 2 λ 1 , we have

λ k + 1 = max V codim ( V ) = k min Ψ V Ψ 0 L ϕ Ψ , Ψ | Ψ | L 2 2 .
(27)

with the maximum attained for

V = span { E 1 , , E k } ,
(28)

where E j denotes the eigenfunction associated with a negative eigenvalue λ j . Such min-max characterizations of eigenvalues of L ϕ -like operators are made more explicit in Sec. Vfor the multidimensional case; see Remark V.1–(ii).

Let g be a τ-periodic continuous differentiable function such that g / g L , and ϕ in L 2 ( T ) . Let us consider

A φ = d φ d s + g g φ ,
(29)

whose domain is given by

D ( A ) : = { φ L 2 ( T ) : φ AC ( T ) , d φ d s + g g φ L 2 ( T ) } ,
(30)

where AC ( T ) denotes the collection of absolutely continuous functions on T .

Then, an adaptation of the arguments of Ref. 46 (Lem. 3.5) shows that A generates a weighted left translation semigroup S(t) on L 2 ( T ) . The corresponding operator L ϕ is then given by

L ϕ [ Ψ ] ( r ) = 1 τ T g ( s + r ) g ( s ) ϕ ( s + r ) Ψ ( s ) d s .

By assumptions on g and ϕ, assumption (A1) is satisfied and L ϕ is compact; see Lemma III.1. If g and ϕ are furthermore real-valued, then assumption (A2) is satisfied and thus L ϕ is self-adjoint as well; see again Lemma III.1.

An analysis of the resolvent of A reveals that the corresponding operators Pn defined in (6) are each rank-one, and that the En's are real. Theorem III.2 can be thus applied here when ϕ lies in D(A). Instead of providing a proof of this statement in this rather general setting, we focus on an important special case for applications, namely, the case in which g 1 , i.e., the case of the left circular translation group T(t) defined by

T ( t ) ϕ ( s ) = ϕ ( t + s ) , mod ( τ ) .
(31)

We have in this case

L ϕ [ Ψ ] ( r ) = 1 τ T ϕ ( s + r ) Ψ ( s ) d s ,
(32)

and the corresponding Fredholm integral equation of the second kind47 possesses explicit solutions such as those described by the following theorem. These solutions actually provide a spectral representation of the Fourier transform in terms of both its power and phase spectra. This is the object of the next theorem.

Theorem IV.1. Let A be the generator of the left circular group ( T ( t ) ) t on L 2 ( T ) defined by (31). Assume ϕ to be a real-valued function that lies in D(A) given by (30) with g 1 . Then, the eigenvalues of the operator L ϕ are real-valued, form a discrete set { λ k } , where k in characterizes the Fourier frequency 2 i π k / τ . More precisely, to each frequency f k = k τ corresponds an eigenvalue λk such that

| λ k | = | ϕ ̂ ( f k ) | .
(33)

Furthermore, the corresponding eigenfunctions in L 2 ( T ) are given by

E k ( t ) = 2 cos ( 2 π f k t + θ k )
(34)

with

{ θ k = 1 2 arg ϕ ̂ ( f k ) , mod ( 2 π ) if λ k 0 , θ k = 1 2 arg ϕ ̂ ( f k ) + π 2 , mod ( 2 π ) if λ k < 0 ,
(35)

where arg ( z ) denotes the principal value48 of the argument of the complex number z.

Remark IV.1. The eigenfunctions are thus phase-shifted sinusoids, for which the phase θk relates to the phase spectrum of the data ϕ. The eigenvalues provide the power spectrum. So far, for a one-dimensional signal ϕ, one can recover the Fourier spectrum (power and phase) from the spectral elements of L ϕ . The operator representation of the Fourier spectrum via the operator L ϕ allows us to propose an innovative generalization for multidimensional signals (see Sec. V) which turns out to be particularly useful for the inverse stochastic-dynamic modeling of spatio-temporal scalar fields; see Secs. VIII and IX.

Proof. As explained earlier, assumptions (A1) and (A2) are satisfied. Let us now check assumption (A3) of Theorem III.2. For that let us determine the resolvent of the generator A of the left circular group ( T ( t ) ) t defined by (31). Note that the operator A is the differentiation operator defined in (29) with g 1 .

The resolvent R ( λ , A ) of A exists for all λ { 0 , ± i ω 1 , ± i ω 2 , } , with ω k = 2 π k / τ , and is obtained by solving the differential equation

λ u u ̇ = ψ , ψ L 2 ( T ) .
(36)

The variation-of-constant formula gives

u ( s ) = e λ ( s t ) u ( t ) t s e λ ( s r ) ψ ( r ) d r .
(37)

Taking s = t + τ , and requiring u ( t + τ ) = u ( t ) , we get

( 1 e λ τ ) u ( t ) = t t + τ e λ ( t + τ r ) f ( r ) d r .
(38)

The factor on the left hand side is invertible if and only if λ does not coincide with 2 π i n / τ . In such a case

u ( t ) = ( 1 e λ τ ) 1 t t + τ e λ ( t r ) ψ ( r ) d r = ( 1 e λ τ ) 1 0 τ e λ τ ψ ( t + r ) d r .
(39)

The right-hand side of this formula maps L 2 ( T ) into the Sobolev space W 1 , 2 ( T ) .

Now let p be arbitrary in and let us take ψ ( r ) : = e i ω p r . Then

u ( t ) = e i ω p t ( 1 e λ τ ) 1 0 τ e ( i ω p λ ) r d r = ( e ( i ω p λ ) τ 1 ) 1 e λ τ e i ω p t i ω p λ = e i ω p t i ω p λ ,
(40)

since ω p τ 2 π , by definition. Thus, by expanding now any arbitrary Ψ in L 2 ( T ) into its Fourier expansion, we get

R ( λ , A ) ψ = p = ψ ̂ ( f p ) i ω p λ e i ω p θ
(41)

with f p = p / τ .

Therefore, the residue Pk at the pole μ k = i ω k [see (6)] is given, for any k in , by

P k ψ ( t ) = ψ ̂ ( f k ) e i ω k t ,
(42)

and the assumption (A3) of Theorem III.2 is satisfied with e k ( t ) = e i ω k t .

The identity (21) becomes in this case

λ k E ̂ k ( f k ) = ϕ ̂ ( f k ) E ̂ k ( f k ) ¯ .
(43)

Taking the modulus on both sides of (43) gives

| λ k | = | ϕ ̂ ( f k ) | .
(44)

Let us seek solutions to (43) under the form

E k ( t ) = cos ( 2 π f k t + θ k ) = cos ( 2 π f k ( t θ k 2 π f k ) ) .
(45)

Now by using the Fourier time-shift property, i.e., that the Fourier transform of x ( t θ k ) is given by x ̂ ( f ) e 2 π i f θ , we obtain from the Fourier transform of cos ( 2 π f k t ) that

E k ̂ ( f ) = 1 2 δ ( f f k ) exp ( i f θ k f k ) + 1 2 δ ( f + f k ) exp ( i f θ k f k )
(46)

and

E k ̂ ( f k ) = 1 2 exp ( i θ k ) .
(47)

With this expression substituted into (43), one finds

λ k = ϕ ̂ ( f k ) exp ( 2 i θ k )
(48)

and since λk is real, necessarily

{ θ k = 1 2 arg ϕ ̂ ( f k ) , mod ( 2 π ) if λ k 0 , θ k = 1 2 arg ϕ ̂ ( f k ) + π 2 , mod ( 2 π ) if λ k < 0.
(49)

Finally, the factor 2 comes as a normalization constant of Ek, for the norm subordinated to the inner product (14). ◼

Remark IV.2.

  • The negative part of the spectrum of L ϕ comes with eigenmodes of sine-type, namely, if λ k < 0 then
    E k ( t ) = cos ( ω k t 1 2 arg ϕ ̂ ( f k ) + π 2 ) = sin ( ω k t 1 2 arg ϕ ̂ ( f k ) ) .
    (50)
  • Sometimes instead of (35), we will make use of the following equivalent expression of the phase spectrum in terms of the eigenmodes of L ϕ , namely,
    arg ϕ ̂ ( f k ) = arg ( λ k E ̂ k ( f k ) ) arg ( E ̂ k ( f k ) ¯ ) , mod ( 2 π ) .
    (51)

The latter identity is a direct consequence of Eq. (43).

In this section, we consider a block operator matrix generalization of what precedes. We assume that we are given a collection { ϕ i , j ( t ) , 1 i , j d } of τ-periodic functions in L 2 ( T ) . Examples of such functions for the spectral analysis of, e.g., spatio-temporal fields, are discussed in Sec. VI D.

For the moment, let us consider the following operator matrix L defined for each u in ( L 2 ( T ) ) d by:

L u : = ( M 1 Ψ 1 , , M d Ψ d ) , u = ( Ψ 1 , , Ψ d ) ( L 2 ( T ) ) d
(52)

with

M p Ψ p : = q = 1 p 1 L ϕ q , p ( Ψ q ) + q = p d L ϕ p , q ( Ψ q ) , 1 p d ,
(53)

where each L ϕ i , j is given by (32). We have then the following Lemma.

Lemma V.1. The operator L given by (52) is self-adjoint and compact on X : = ( L 2 ( T ) ) d endowed with the inner product given for any u = ( Ψ 1 , , Ψ d ) and v = ( v 1 , , v d ) in X by

u , v X : = n = 1 d Ψ n , v n = τ 1 n = 1 d T Ψ n ( r ) v n ( r ) ¯ d r .
(54)

Proof. Since for each p, q in { 1 , , d } , the function ϕ p , q lies in L 2 ( T ) by assumption, each operator L ϕ p , q given by (32) is compact (see Sec. IV), and thus each operator M p is a (real-valued) compact operator on L 2 ( T ) , as the finite sum of compact operators. The operator L is then compact on X as the finite cartesian product of compact operators, for X endowed with the product topology induced by the inner product (54).

The self-adjoint property of L results from its definition. Indeed,

L u , u X = p = 1 d q = 1 p 1 L ϕ q , p ( Ψ q ) + q = p d L ϕ p , q ( Ψ q ) , Ψ p ,
(55)

which, after similar manipulations as in (16), leads to

L u , u X = p = 1 d ( q = 1 p 1 L ϕ q , p ( Ψ p ) , Ψ q ¯ + q = p d L ϕ p , q ( Ψ p ) , Ψ q ¯ ) = q = 1 d ( p = 1 q 1 L ϕ p , q ( Ψ p ) , Ψ q ¯ + p = q d L ϕ q , p ( Ψ p ) , Ψ q ¯ ) = L u , u X ¯ .
(56)

We are now in position to formulate the main result that extends Theorem IV.1 to the multidimensional case. For that purpose, we introduce for each frequency f k = k / τ , the d × dsymmetrized cross-spectral matrix S ( f k ) whose entries are given by

S p , q k : = { ϕ p , q ̂ ( f k ) if q p , ϕ q , p ̂ ( f k ) if q < p .
(57)

The next theorem shows that the spectrum of the operator L relates naturally to the cross-spectral matrix S ( f ) . At every frequency line, the singular values of S ( f ) gives actually d pairs of eigenvalues of L , and as f is varied one recovers the full spectrum of L .49 

Theorem V.1.Let A be the generator of the left circular group ( T ( t ) ) t on L 2 ( T ) defined by (31). Assume that each ϕ p , q is a real-valued function lying in D(A) given by (30) with g 1 . Then, the eigenvalues of the operator L defined in (52) are real-valued, form a discrete set { λ n } that possesses the following characterization.

For each frequency f k 0 , one can extract 2d eigenvalues (counting multiplicities) from the set { λ n } , d positive and d negative ones that satisfy the following property.

  • (S) For each pair of negative-positive eigenvalues denoted by ( λ p ( f k ) , λ + p ( f k ) ) , there exists a singular value σ p ( f k ) of S ( f k ) such that
    λ + p ( f k ) = λ p ( f k ) = σ p ( f k ) , 1 p d .
    (58)

Furthermore, the eigenfunctions Wk of L in ( L 2 ( T ) ) d possess the following representation:

W k = ( E 1 k , , E d k ) tr
(59)

with

E p k ( t ) = B p k cos ( 2 π f k t + θ p k ) , B p k , θ p k , with f k = k τ , k , p { 1 , , d } .
(60)

Finally, introducing the notations

ϕ a , b ̂ ( f k ) = X a b k e i α a b k , with X a b k = | ϕ a , b ̂ ( f k ) | ,
(61)

the amplitudes B p k and angles θ p k satisfy for each frequency fk and each p in { 1 , , d } , the relation

B p k λ + p ( f k ) = q = 1 p 1 B q k X q p k e i ( θ p k + θ q k + α q p k ) + q = p d B q k X p q k e i ( θ p k + θ q k + α p q k ) .
(62)

Proof. Step 1. We first start with the derivation of an analogue of Eq. (43) in the multidimensional case considered here. The resulting extension of (43) is key to the derivation of the main properties of the spectral elements of the block operator matrix L defined by (52) and (53).

For that purpose, let us consider Wk in ( L 2 ( T ) ) d under the form given in (59) with the E p k to be determined. Due to the operator matrix form of L given by (52) and (53), the eigenvalue problem L W n = λ W n can be equivalently rewritten as

q = 1 p 1 L ϕ q , p ( E q k ) + q = p d L ϕ p , q ( E q k ) = λ k E p k , p { 1 , , d } .
(63)

Let p be fixed in { 1 , , d } and k be in . By reproducing the arguments of the proof of Theorem III.2, we have

λ T e μ k r E p n ( r ) d r = q = 1 p 1 T ( P k ϕ q , p ) ( s ) E q n ( s ) d s + q = p d T ( P k ϕ p , q ) ( s ) E q n ( s ) d s .
(64)

Here again, Pk denotes the residue of the resolvent R ( ξ , A ) of A at the pole μ k = 2 π i k / τ and is thus given by

P k = lim ξ μ k ( ξ μ k ) R ( ξ , A ) .

Now let us recall the expression of the resolvent of R ( λ , A ) derived in the proof of Theorem IV.1. Namely, R ( ξ , A ) exists for all ξ { 0 , ± i ω 1 , ± i ω 2 , } , with ω k = 2 π k / τ , and is given for ψ in L 2 ( T ) by

R ( ξ , A ) ψ = p = ψ ̂ ( f p ) i ω p ξ e i ω p θ .
(65)

Therefore, Pk is given by

P k ψ ( t ) = ψ ̂ ( f k ) e i ω k t
(66)

and thus from (64) we obtain, since E q n ( t ) is a real-valued function [since Wn lies in ( L 2 ( T ) ) d ],

λ E p n ̂ ( f k ) = q = 1 p 1 ϕ q , p ̂ ( f k ) E q n ̂ ( f k ) ¯ + q = p d ϕ p , q ̂ ( f k ) E q n ̂ ( f k ) ¯ ,
(67)

i.e., the aforementioned extension of (43) to the multidimensional case.

Seeking solutions to Eq. (67) of the form (60), and recalling that

E p n ̂ ( f ) = 1 2 δ ( f f n ) exp ( i f θ p n f n ) + 1 2 δ ( f + f n ) exp ( i f θ p n f n ) .
(68)

[see also (46)], we conclude that

( E p n ̂ ( f k ) 0 and E q n ̂ ( f k ) 0 ) ( k = n ) .
(69)

We consider thus k = n, hereafter. Replacing λ + p by λ in Eq. (62) for a moment, it is easy to see that Eq. (62) is a direct consequence of Eq. (67). The former is indeed obtained by taking the real and imaginary parts of Eq. (67) [in which E p k is given by (60)] and forming the complex exponentials by summation. To complete the proof, we are thus left with the proof of property (S) which we deal with next.

Step 2. Note first that the symmetrized cross-spectral matrix S ( f k ) whose entries are defined in (57) is complex and symmetric. The Autonne-Takagi factorization [see Ref. 50 (Chap. VIII, Théorème 92) and Ref. 51] applied to S ( f k ) ensures thus the existence of a unitary matrix Uk such that

S ( f k ) = U k Σ k U k T ,
(70)

where Σ k = diag { σ 1 ( f k ) , , σ d ( f k ) } is a diagonal matrix with real and nonnegative entries. These entries are the singular values of S ( f ) at the frequency f = fk.

In particular, denoting by U the conjugate transpose of U, we have

S ( f k ) S ( f k ) ¯ = U k Σ U k T U ¯ Σ k U k = U k Σ 2 U k
(71)

and therefore the matrix S ( f k ) S ( f k ) ¯ is positive semi-definite, and its eigenvalues are given by the σ j ( f k ) 2 , for 1 j d .

On the other hand, Eq. (67) may be rewritten using (59) as

λ W k ̂ ( f k ) = S ( f k ) W k ̂ ( f k ) ¯ ,
(72)

which leads to

S ( f k ) S ( f k ) ¯ W k ̂ ( f k ) = S ( f k ) ( S ( f k ) W k ̂ ( f k ) ¯ ¯ ) = λ S ( f k ) W k ̂ ( f k ) ¯ = λ 2 W k ̂ ( f k ) ,
(73)

and thus λ 2 is an eigenvalue of S ( f k ) S ( f k ) ¯ .

To summarize, λ is (up to a sign) the square root of an eigenvalue of S ( f k ) S ( f k ) ¯ and is thus also a singular value σ j ( f k ) of S ( f k ) , up to a sign. Now since k = n in Eq. (67) (see step 1), there is for each frequency fk, d relations (67), one for each p. Considering the positive and negative eigenvalues that coexist at each frequency f k 0 , we infer thus (58). ◼

Remark V.1.

  • It should be noted that an a singular value decomposition (SVD)of S ( f k ) is not automatically its Autonne-Takagi factorization. The proof of the Autonne-Takagi factorization provides the details of how an SVD of a complex symmetric matrix can be turned into a symmetric SVD; see the proof of Ref. 52 (Theorem 3.1).

  • It is also worth mentioning that due to a general result dealing with the characterization of singular values of complex symmetric matrices,53the singular values of S ( f k ) possess a variational characterization in terms of the quadratic form
    Q ( x ) : = Re ( x T S ( f k ) x ) , x d .
    (74)
    In particular, by application of Ref. 53 (Corollary 5)and Theorem V.1, we obtain for 0 p d 1 the following variational characterization by maximizing over real subspaces, namely:
    max V dim ( V ) = p + 1 min x V x = 1 Re ( x T S ( f k ) x ) = λ + 1 + p ( f k ) ,
    (75)
    and
    max V codim ( V ) = p min x V x = 1 Re ( x T S ( f k ) x ) = λ 1 + p ( f k ) .
    (76)

Theorem V.1 and Remark V.1–(ii) ensure that the eigenvalues of the operator L (at the frequency f = fk) relate naturally to an energy via the quadratic form Q given in (74). Contrarily to the one-dimensional case [see (35)] and as (62) shows, the θ p k 's do not relate trivially to the multidimensional phase information contained in symmetrized cross-spectral matrix S ( f k ) .

Nevertheless, the relation (67) allows us to extract a useful phase information from S ( f k ) by noting that (provided that B p k 0 )

λ + p ( f k ) E p k ̂ ( f k ) ( E p k ̂ ( f k ) ¯ ) 1 = ϕ p , p ̂ ( f k ) + R p ( f k ) ,
(77)

with

R p ( f k ) = q = 1 p 1 ϕ q , p ̂ ( f k ) γ p q k + q = p + 1 d ϕ p , q ̂ ( f k ) γ p q k
(78)

and γ p q k : = E q k ̂ ( f k ) ( E p k ̂ ( f k ) ) 1 ¯ . Obviously, a similar relation holds for λ p ( f k ) (involving the appropriate E p k ).

By taking the argument on both sides of (77), we obtain (up to a multiple of 2 π )

arg ( λ + p ( f k ) E p k ̂ ( f k ) ) arg ( E p k ̂ ( f k ) ¯ ) = arg ( ϕ p , p ̂ ( f k ) + R p ( f k ) ) .
(79)

We denote by η + p ( f k ) the LHS of (79).

On the other hand, (68) gives E p k ̂ ( k ) = exp ( i θ p k ) / 2 , which leads to

2 θ p k = η + p ( f k ) = arg ( ϕ p , p ̂ ( f k ) + R p ( f k ) ) .
(80)

We understand thus that η + p ( f k ) gives (twice) the phase shift, θ p k , contained in the eigenfunctions [see (60)], and also provides a measure about the perturbation brought by the “cross-functions” ϕ p , q ( p q )—as encapsulated into R p ( f k ) —to the phase spectrum associated with ϕ p , p when considered as isolated [i.e., arg ( ϕ p , p ̂ ( f k ) ) ].

This perturbation is dependent on the eigenfunctions of the operator L as the expression of R p ( f k ) given in ( 78 ) shows. When R p ( f k ) = 0 , on recovers the phase spectrum of ϕ p , p as in the unidimensional case; see Remark IV.2–(iii). In what follows, we will look at the quantities, λ + p ( f k ) and η + p ( f k ) (resp. η + p ( f k ) ) as fk is varied; the former providing, in the multidimensional case, a notion of power spectrum (as associated with the form Q ) and the latter a notion of (perturbed) phase spectrum.

Section VI below focuses on the important case for applications in which the functions ϕ p , q are built from correlation functions. The concepts of DAH power and DAH phase spectra are further detailed in Sec. VI, addressing their numerical computation, while the determination and visualization of the underlying modes are illustrated by analytical examples in Sec. VII.

To simplify the presentation, we consider throughout this section a finite-dimensional Euclidean space H .

In the case of a deterministic flow ( S t ) t 0 acting on H and having an ergodic invariant measure μ supported by an attractor A ,54,55 let us recall that the correlation function ρ f , g ( t ) associated with two (sufficiently smoothed) observables f , g : H is given by

ρ f , g ( t ) = H f ( x ) g ( S t x ) d μ , t 0 , x H .
(81)

In the case of an ergodic stochastic system, the correlation function is given by

ρ f , g ( t ) = H f ( x ) P t g ( x ) d μ , t 0 , x H ,
(82)

where ( P t ) t 0 denotes the Markov semigroup possessing μ as (unique) ergodic invariant measure.56 In each case, the function ρ f , g ( t ) is called a correlation function subordinated to the observables f and g. Depending on these observables, various higher-order moment statistics can be considered.

In the deterministic setting, if μ is a Sinaï-Ruelle-Bowen measure,54,55,57 ρ f , g ( t ) corresponds to the more familiar cross-correlation coefficient at lag t (for the observables f and g) which takes then the equivalent form58 

ρ f , g ( t ) = lim T 1 T 0 T f ( S s x ) g ( S t + s x ) d s
(83)

for almost every x in the basin of attraction of A .55 

For an ergodic stochastic system, it takes the form

ρ f , g ( t ) = lim T 1 T 0 T f ( X t x ) g ( X t + s x ) d s .
(84)

P -almost surely and for every x in H , when, e.g., the Markov semigroup Pt associated with the Markov process X t x is strong Feller and irreducible; see Ref. 56.

We place ourselves in one of the working framework described in Sec. VI A above—i.e., either stochastic or deterministic—and apply the theory of Sec. IV to ϕ ( t ) = ρ f , g ( t ) for a choice of two observables, f and g. Note that if one chooses f = g, one recovers the familiar notion of lagged autocorrelation of, e.g., the univariate time series t f ( S t x ) . For this reason, we will sometimes, when f g , name ρ f , g ( t ) as a cross-correlation function, to emphasize the difference with the notion of autocorrelation.

Obviously, the main assumption on which the theory in Secs. III and IV relies on is the periodicity of ϕ ( t ) , which is not guaranteed, in practice, except for periodic trajectories of course. However, motivated by physical applications, we are interested with mixing systems, i.e., systems for which a decay of correlations takes place. For stochastic (resp. deterministic) mixing systems, see, e.g., Ref. 56 (resp. Refs. 55 and 58) and references therein.

For such systems, one should, in practice, choose τ large enough so that the decay of correlations has been sufficiently well-resolved. To apply the theory, one then needs to periodize ϕ ( t ) , and the resulting periodization φ plays a role in the definition of the operator L φ . To understand it, let us consider an interval I = [ α , α + τ ] with α in . We define

φ ( t ) : = { ϕ ( t ) , t [ α , α + τ ) , ϕ ( t n τ ) , t [ α + n τ , α + ( n + 1 ) τ ) , n + ϕ ( t + n τ ) , t ( α + ( n 1 ) τ , α + n τ ] , n .
(85)

The left circular shift T(t) defined in (31) applied to φ thus gives the following expression in terms of ϕ, for r lying in the interval I:

T ( r ) φ ( s ) : = { ϕ ( r + s ) , s [ α , α + τ r ) , ϕ ( r + s τ ) , s [ α + τ r , α + τ ) .
(86)

In other words, the kernel of the integral operator L φ is modified compared with its expression given in (32).

As a choice of the base interval I representing T , a reasonable one consists of taking α = τ / 2 , i.e.,

I = [ τ 2 , τ 2 ] .
(87)

This choice as a representation of T relies in part on the symmetry of ϕ ( t ) = ρ f , g ( t ) , i.e., ρ f , g ( t ) = ρ f , g ( t ) . Thus, ϕ ( τ / 2 ) = ϕ ( τ / 2 ) . Other choices, like α = 0, results—due to the assumed decay of correlations—with “jumps” that occur at the junction points in the periodization of ϕ given in (85).

These jumps do not constitute, however, an obstruction to apply the theory. Theorem IV.1 allows indeed for a discontinuity of this type for φ used in the definition of L φ . The reason is tied to the characterization of D(A) given in (30) (with g 1 ) that shows that a function φ with a discontinuity is allowed, as long as φ defined in (85) is absolutely continuous and φ lies in L 2 ( T ) .

In Sec. VI D, we consider the choice of the base interval I given by (87) in the case of a multidimensional DAH spectrum obtained from a collection of correlation functions. We turn next to the case I = [ 0 , τ ] , i.e., α = 0, for which the framework adopted in this article allows us to recover instructive results from the literature on Hankel matrices. In that respect, we address first the discretization of the operator L φ .

Given a uniform partition t 0 = 0 < < t N 1 < τ = t N , we consider the family of functions

χ j ( s ) : = { 1 , s [ t j , t j + 1 ) , 0 , else . 0 j N 1.
(88)

We simply denote by ρ, the correlation function ρ f , g . By the definition (86) of the left circular shift acting on the periodization φ in (85) (with α = 0), we have from (32)

L φ ( Ψ ) ( r ) = 1 τ ( 0 τ r ρ ( s + r ) Ψ ( s ) d s + τ r τ ρ ( r + s τ ) Ψ ( s ) d s ) .
(89)

Now let us approximate ρ by a piecewise constant function f, as follows:

f ( s ) = ρ ( t j ) , if s [ t j , t j + 1 ) .
(90)

Take t j = j δ τ with δ τ = τ / N , then

L φ [ χ j ] ( r ) δ τ τ ρ ( ( j + k ) δ τ ) mod ( τ ) , if r [ t k , t k + 1 ) .
(91)

Denoting ρ ( l δ ) by ρl, the operator L φ can be thus approximated by the following Hankel N × N matrix:

H = 1 N ( ρ 0 ρ 1 ρ 2 ρ N 1 ρ 1 ρ 2 ρ 3 ρ 0 ρ 2 ρ 1 ρ N 1 ρ 0 ρ 1 ρ N 2 ) .
(92)

Without any surprise, H is a left circulant matrix, which we denote sometimes as

H = l circ ( ρ 0 , ρ 1 , , ρ N 1 ) ,
(93)

or in other words, the rows of H are obtained by successive shifts to the left of one position starting from the row r = ( ρ 0 , ρ 1 , , ρ N 1 ) as a first row.

Let T be the right circulant matrix which has the same first row than the left circulant matrix H. Then, T is Toeplitz and H = P T for the permutation matrix P = 1 J n 1 , where J N 1 is the reversal ( N 1 ) × ( N 1 ) matrix obtained by flipping the ( N 1 ) × ( N 1 ) identity matrix I N 1 from left to right. Hence, H and T have identical singular values. However, circulant matrices and real anticirculant matrices (which are also real symmetric) are normal matrices, so they can be unitarily diagonalized and their singular values are the moduli of their eigenvalues. Therefore, the eigenvalues H and T have identical moduli. Finally, since real anticirculant matrices are real symmetric, H has real eigenvalues.

The spectral theory of circulant matrices is well established although somehow scattered in the literature; see, e.g., Refs. 59–61. The inherent periodicity of circulant matrices relates them to Fourier analysis and group theory. In contrast to the aforementioned standard circulant Toeplitz matrix T arising in Fourier analysis (see, e.g., Ref. 62 (Def. 5.12), the left-circulant matrix H is not diagonalized by the unitary Fourier matrix Ref. 61 (Theorem 3.5) and, as Proposition VI.1 summarizes, its eigenvectors actually depend on the data and more exactly relate to the phase spectrum as a counterpart of Theorem IV.1 and a consequence of formula (43).

Let us introduce the Fourier frequency ω k = 2 π k / N , and the periodogram of { ρ j } , namely,

I N ( ω k ) : = 1 N | j = 0 N 1 ρ j e i j ω k | 2 , 0 k N 1
(94)

with i 2 = 1 .

We have then the following:

Proposition VI.1. The eigenvalues { λ k } k { 0 , , N 1 } of the Hankel matrixHgiven by (92) can be arranged in terms of Fourier frequencies so that

λ 0 = N 1 n = 0 N 1 ρ n , λ N 2 = N 1 n = 0 N 1 ( 1 ) n ρ n , if N is even ,
(95)
and for
1 k [ N 1 2 ] ,
λ k = λ N k = N 1 / 2 I N ( ω k ) ,
(96)

where [ x ] denotes the largest integer less than or equal to x.

Furthermore, for each pair ( λ k , λ N k ) corresponding thus to a Fourier frequency ω k = 2 π k / N , the corresponding pair of eigenvectors ( v k , v N k ) satisfies

arg ( ρ ̂ ( ω k ) ) = arg ( λ j v j ̂ ( ω k ) ) arg ( v j ̂ ( ω k ) ¯ ) , with j { k , N k } , 1 k [ N 1 2 ] ,
(97)

and where u ̂ ( ω ) denotes the discrete Fourier transform—assessed at the frequency ω—of a vector u.

Proof. The identity (96) results from (33) of Theorem IV.1. For a direct proof of (95) and (96), we refer to Ref. 60 (Lemma 1). The identity (97) is a consequence of (51) in the discrete setting. ◼

Remark VI.1. Note that to the best of the authors' knowledge, the identity (97) has not been previously derived in the literature concerned with Hankel matrices.

In what follows, the eigenvalues of the matrix H when ranked from the lowest to the highest resolved frequency form what we call the DAH power spectrum. The DAH phase spectrum is formed by the phase of the corresponding eigenvectors, also ranked per frequency. Due to Proposition VI.1, these concepts coincide for one-dimensional signals with the standard notions of power and phase spectra [see also Remark IV.2–(ii)], whereas their multidimensional generalization following Sec. V and application to a collection of cross-correlation functions lead to more subtle but still useful interpretations. This generalization is discussed in further details in Sec. VI D.

Remark VI.2. Let us introduce the standard cyclic permutation matrix

K : = ( 0 1 0 0 0 0 0 0 0 1 0 1 1 0 0 0 ) .
(98)

Then,His naturally associated with the Krylov space K ( K , ρ ) generated by the matrixKgiven by (98) and the vector ρ = ( ρ 0 , , ρ N 1 ) tr , namely,

K ( K , ρ ) = span ( ρ , K ρ , K 2 ρ , , K N 1 ρ ) .
(99)

Indeed, the j-th column ofHis exactly given by K j 1 ρ .

Remark VI.3.

  • Proposition VI.1 remains valid when other base intervals I than [ 0 , τ ] are used. The modification consists essentially of replacing accordingly the summation bounds in the definition of the periodogram IN given by (94).

  • Let t X t x denote a stochastic process emanating from x in H : = q ( q > 1 ) and generated by a system of SDEs. Given two observables f , g : H , and recalling that ρ denotes the correlation function ρ f , g given by (82), we infer from the formula (96) that the λk's provide the(one-sided) cross-power spectrum, Γ f , g , of the real-valued stochastic processes t f ( X t x ) and t g ( X t x ) .

    The DAH eigenvectors provide furthermore the(one-sided) cross phase spectrumvia the formula (97), in the sense that
    Γ f , g ( ω k ) N | λ j | = exp ( arg ( λ j v j ̂ ( ω k ) ) arg ( v j ̂ ( ω k ) ¯ ) ) , with j { k , N k } , 1 k [ N 1 2 ] .
    (100)

Given a base interval I to be [ α , α + τ ] , one considers the periodization φ f , g = φ given in (85) for ϕ = ρ f , g , with ρ f , g denoting a correlation function given either by (83) or by (84), depending on the context.

The operator L ϕ defined in (32) then becomes

L φ f , g ( Ψ ) ( r ) = 1 τ ( α α + τ r ρ f , g ( s + r ) Ψ ( s ) d s + α + τ r α + τ ρ f , g ( r + s τ ) Ψ ( s ) d s ) ,
(101)

where Ψ is any arbitrary function of L 2 ( T ) .

By choosing I = [ τ / 2 , τ / 2 ] (i.e., α = τ / 2 ), a uniform partition

τ 2 = t M + 1 t 1 t 0 = 0 , 0 t 1 t M 1 τ 2 = t M ,
(102)

and elementary functions such as in (88), each operator L φ i , j is approximated by the following Hankel ( 2 M 1 ) × ( 2 M 1 ) matrix:

H ( i , j ) = ( ρ M + 1 ( i , j ) ρ M + 2 ( i , j ) ρ 0 ( i , j ) ρ 1 ( i , j ) ρ M 1 ( i , j ) ρ M + 2 ( i , j ) ρ M + 1 ( i , j ) ρ M + 2 ( i , j ) ρ 0 ( i , j ) ρ M + 1 ( i , j ) ρ 1 ( i , j ) ρ M + 2 ( i , j ) ρ 0 ( i , j ) ρ M 1 ( i , j ) ρ M + 1 ( i , j ) ρ M 1 ( i , j ) ρ M + 1 ( i , j ) ρ M + 2 ( i , j ) ρ 0 ( i , j ) ρ M 2 ( i , j ) ) .
(103)

The operator L defined in (52)(53)—in which each operator L ϕ p , q therein is replaced by L φ p , q given in (101) above—is then approximated by the following block-Hankel matrix C formed by d2 blocks, C ( i , j ) , each of size ( 2 M 1 ) × ( 2 M 1 ) and given as follows:

C ( i , j ) = H ( i , j ) , if i j , C ( i , j ) = H ( j , i ) , if j < i
(104)

for (i, j) in { 1 , , d } 2 .

Hereafter, we use M = 2 M 1 for concision, reindexing possibly from 1 to M the string { M + 1 , , M 1 } . By construction and in agreement with Lemma V.1, the matrix C is real symmetric, not circulant with each block that is Hankel.

Since the matrix C is the discretization of the operator L [defined in (52) and (53) and in which each operator L ϕ p , q therein is replaced by L φ p , q given in (101)], Theorem V.1 is enlightening concerning the spectrum of C . In particular, the eigenvalues of C relate according to (58) to the singular values of the symmetrized cross-spectral matrix S δ ( f k ) whose entries are still given by (57) except that ϕ p , q ̂ ( f k ) is obtained here as the discrete Fourier transform (evaluated at the frequency fk) of the correlation function ρ p , q ( t ) , after its periodisation following Sec. VI B.

As noted earlier for the operator L and applied here to the matrix C , Theorem V.1 and Remark V.1–(ii) ensure that the eigenvalues of C relate naturally to the quadratic form Q δ given by

Q δ ( x ) : = Re ( x T S δ ( f k ) x ) , x d
(105)

via the variational principle given in (75)(76), in which S δ ( f k ) replaces S ( f k ) . The eigenvalues of C give thus account on an energy distribution per frequency, as measured through the quadratic form Q δ .

Grouping per frequency the absolute value of the eigenvalues of the matrix C , and displaying those eigenvalues from the lowest to the highest resolved frequency, we obtain what we call the multidimensional DAH power spectrum, omitting the adjective “multidimensional” when this is clear from the context.

Due to Theorem V.1, the eigenvalues of the grand matrix C come in pairs organized per (non-zero) frequencies according to (58). However, the formula (58) is not used in practice and instead eigenvalues of the d M × d M matrix C are computed directly, and listed as the set of eigenvalues ( λ j ) 1 j d M . To group these eigenvalues per Fourier frequency and to determine the DAH power spectrum, we then proceed as follows:

  1. Compute the eigenvalues of C listed as ( λ j ) 1 j d M , and the associated eigenvectors Wj.

  2. Theorem V.1 ensures that for each j in { 1 , , d M } , the eigenvector Wj possesses the representation (59) where each E p j is a M -dimensional row vector that is explicitly associated with a Fourier frequency f according to (60) which translates in its discrete form as
    E p j ( s ) = B p j cos ( f s + θ p j ) , 1 s M ,
    (106)

    and where p is in { 1 , , d } .

    Note that f in the above formula takes value in the following set of Fourier frequency:
    f = 2 π ( 1 ) M 1 , = 1 , , M + 1 2 .
    (107)
  3. Compute the discrete Fourier transform E ̂ p j of E p j and group the j's for which the average power spectral density, Γ ¯ j : = d 1 p = 1 d | E ̂ p j | 2 , exhibits a dominant peak at a frequency f .

  4. From step 3 above, we have thus formed for each Fourier frequency f the following subset of indices in { 1 , , d M } :
    J ( f ) : = { j : s . t Γ ¯ j is peaked at f = f } .
    (108)
  5. Form for each ranging from 1 to ( M + 1 ) / 2 , the discrete set
    P : = { | λ j | , : j J ( f ) } .
    (109)

    The collection of the P for ranging from 1 to ( M + 1 ) / 2 denotes the DAH power spectrum.

Note that due to (58), J ( f ) is composed of 2d indices when 1 and of d indices if = 1 such that the J ( f ) 's form a partition of the total set of indices, { 1 , , d M } .

In applications, the eigenvectors of C are called the DAH modes (DAHMs), and the E p j are called DAHM snippets.

The multidimensional DAH phase spectrum is obtained by grouping per frequency the η + p ( f k ) [resp. η p ( f k ) ] as discussed at the end of Sec. V. Note that, in practice, the η + p ( f k ) [resp. η p ( f k ) ] are calculated by evaluating the LHS of (79).

More precisely, we form

Φ k j ( f ) : = arg ( λ j E ̂ k j ( f ) ) arg ( E ̂ k j ¯ ( f ) ) ,
(110)

where E ̂ k j and E ̂ k j ¯ denote the discrete Fourier transform of the DAHM snippet, E k j , and its conjugate, respectively. The DAH phase spectrum is then obtained as the collection of the following discrete set:

Φ : = { Φ k j ( f ) : j J ( f ) }
(111)

as varies in from 1 to ( M + 1 ) / 2 .

Finally, given a d-dimensional time series X ( t ) from which the matrix C is formed, we describe below another natural byproduct that can be computed from the spectral elements of C , namely, the DAH coefficients (DAHCs) obtained by projection of X ( t ) onto the orthogonal set formed by the DAHMs.

More precisely, a DAHC is a time series, ξ j ( t ) , obtained according to the formula

ξ j ( t ) : = s = 1 M p = 1 d X p ( t + s 1 ) E p j ( s ) , 1 s M , 1 t N , with N = N M + 1.
(112)

As shown in Theorem V.1, for each frequency f 0 corresponds d pairs of DAH eigenvalues; see (58). For each such pair, one can associate a pair of DAHCs that is, as explained in Sec. VIII A for the continuous time setting, constituted of time series that are narrowband about the frequency f. As we will see in Sec. VIII, this property shared by the DAHCs plays a key role in the inverse modeling approach proposed in this article. The fact that the DAHCs constitute a key ingredient for modeling relies also on their usage in natural reconstruction formulas allowing for a reconstruction of the full dataset or a part of it, possibly over a targeted frequency band; see  Appendix.

The modeling aspects exploiting DAHCs will be discussed later. For the moment, we turn next to basic examples aimed at providing first intuitions about DAH power and phase spectra, as well as about the related DAHMs and DAHCs.

We consider the following one-dimensional time-periodic scalar field, given for x in [ 0 , L ] and t > 0 by

q ( x , t ) = Re ( exp ( i ( k ( x ) x ω c t ) ) )
(113)

with a wave length k(x) depending on the spatial location. Given a uniform discretization [ 0 , L ] by Nx points for a mesh size δx, we choose hereafter k(x) to be given according to

k ( x ) = x δ x p α , α = N x m , 1 m N x .
(114)

Here, p denotes the integer for which the quotient x / ( α δ x ) is within the roundoff error of that integer. The parameter ωc is the Fourier frequency given by

ω c = 2 π ( l 1 ) K

with K = 2 K 1 and 1 l M . For the experiments reported below, we set L = 5 π , Nx = 42, K = 15, l = 3, and m = 19. For simplicity, a time-stepping, δ t = 1 , is chosen to discretize the time-variable t in (113). The resulting patterns are shown in Fig. 1(a). As a benchmark, a traveling wave with k(x) replaced by its mean value k ¯ is shown in Fig. 1(b). The effect of k(x) is thus to distort a spatial cosine wave, cos ( k ¯ x ω c t ) , by introducing some local oscillations occurring at scales smaller than 1 / k ¯ .

FIG. 1.

Spatio-temporal fields. Panel (a) shows a time evolution of q(x, t) with modulated wavelength, k(x), given by (114). Panel (b) shows the field obtained when replacing k(x) by its mean value k ¯ ; see text for more details. (a) Modulated wave and (b) travelling wave.

FIG. 1.

Spatio-temporal fields. Panel (a) shows a time evolution of q(x, t) with modulated wavelength, k(x), given by (114). Panel (b) shows the field obtained when replacing k(x) by its mean value k ¯ ; see text for more details. (a) Modulated wave and (b) travelling wave.

Close modal

Compared to the case of a pure traveling wave shown in Fig. 1(b), some channels in Fig. 1(a) still exhibit collectively some traveling wave patterns (between, e.g., x 14 and x = 5 π ), while others exhibit now standing-like waves patterns as it can be observed between x 4 and x 7 .

The corresponding DAH power spectrum and DAH phase spectrum are shown in Figs. 2(a) and 2(b), respectively. Their computations have been obtained by diagonalizing the matrix C given in (104) and following the procedure described in Sec. VI D, in which the cross correlations are estimated from q ( x j , t n ) obtained from (113), with x j = ( j 1 ) L / N x , 1 j N x , and tn = n, for 0 n 3 × 10 3 . The corresponding DAHMs and DAHCs are shown in Fig. 4.

FIG. 2.

Multidimensional DAH power and phase spectra. The DAH spectra shown here are those associated with the modulated wave shown in Fig. 1(a). On each vertical line in panel (a) [resp. (b)] above a frequency f / ( 2 π ) [with f given by (107)] is shown the discrete set P (resp. Φ ) given by (109) [resp. (111)]. For the underlying calculations, we chose M = 43 to estimate the cross correlations and form the corresponding matrix C given in (104). Each non-zero frequency is associated with 42 eigenpairs that correspond to the total of d = N x = 42 channels used to generate the field q(x, t) from (113). (a) DAH power spectrum and (b) DAH phase spectrum.

FIG. 2.

Multidimensional DAH power and phase spectra. The DAH spectra shown here are those associated with the modulated wave shown in Fig. 1(a). On each vertical line in panel (a) [resp. (b)] above a frequency f / ( 2 π ) [with f given by (107)] is shown the discrete set P (resp. Φ ) given by (109) [resp. (111)]. For the underlying calculations, we chose M = 43 to estimate the cross correlations and form the corresponding matrix C given in (104). Each non-zero frequency is associated with 42 eigenpairs that correspond to the total of d = N x = 42 channels used to generate the field q(x, t) from (113). (a) DAH power spectrum and (b) DAH phase spectrum.

Close modal

We first comment on the DAH power spectrum and DAHMs. The main noticeable feature of the DAH power spectrum shown in Fig. 2(a) [resp. Fig. 3(a)] is the presence of a peak located at the (temporal) frequency close to f c = ω c / ( 2 π ) (due to our discretization), the characteristic frequency of the spatio-temporal field q(x, t) given by (113) [resp. (113) with k(x) replaced by k ¯ ]. Recalling the variational characterization of the eigenvalues of C , one has thus here an illustration of the relevance of the energy Q δ given by (105) for indicating the presence of a dominant oscillation; the latter being expressed by the presence of a dominant pair of eigenvalues [i.e., with the highest max-min for Q δ ; see Remark V.1–(ii)] at that frequency.

FIG. 3.

Multidimensional DAH power and phase spectra. The DAH spectra shown here are those associated with the pure traveling wave shown in Fig. 1(b). On each vertical line in panel (a) [resp. (b)] above a frequency f / ( 2 π ) [with f given by (107)], is shown the discrete set P (resp. Φ ) given by (109) [resp. (111)]. For the underlying calculations, we chose M = 43 to estimate the cross correlations and form the corresponding matrix C given in (104).

FIG. 3.

Multidimensional DAH power and phase spectra. The DAH spectra shown here are those associated with the pure traveling wave shown in Fig. 1(b). On each vertical line in panel (a) [resp. (b)] above a frequency f / ( 2 π ) [with f given by (107)], is shown the discrete set P (resp. Φ ) given by (109) [resp. (111)]. For the underlying calculations, we chose M = 43 to estimate the cross correlations and form the corresponding matrix C given in (104).

Close modal

This dominant pair of eigenvalues come with an energetic oscillatory pair of DAHCs as shown in the panel located at the intersection of the bottom row and third column of Fig. 4 (resp. Fig. 5) for f = 0.071 ω c / ( 2 π ) ; compare with the DAHCs of the third column for other frequencies.

FIG. 4.

DAHMs and DAHCs from the modulated wave of Fig. 1(a). The DAHMs (1st and 2nd columns) are represented here according to their visualization adopted in this article; see text. Each row of the 1st and 2nd columns, shows a pair of DAHMs at a given frequency f such as indicated. There, the horizontal axes represent the temporal embedding dimension M = 2 M 1 (here M = 43), while the vertical ones indicate the spatial channel 1 k d with d = N x = 42 . The DAHCs shown in the 3rd column, are computed according to (112) where the E p j therein are selected for the indicated frequency f, i.e., j lies in J ( f ) given by (108). The x-axis represents time [t in (112)] and the y-axis represents the amplitude of the DAHCs. Note that the DAHMs and DAHCs shown here correspond, for each f, to the dominant pair of DAH eigenvalues located above f in Fig. 2(a).

FIG. 4.

DAHMs and DAHCs from the modulated wave of Fig. 1(a). The DAHMs (1st and 2nd columns) are represented here according to their visualization adopted in this article; see text. Each row of the 1st and 2nd columns, shows a pair of DAHMs at a given frequency f such as indicated. There, the horizontal axes represent the temporal embedding dimension M = 2 M 1 (here M = 43), while the vertical ones indicate the spatial channel 1 k d with d = N x = 42 . The DAHCs shown in the 3rd column, are computed according to (112) where the E p j therein are selected for the indicated frequency f, i.e., j lies in J ( f ) given by (108). The x-axis represents time [t in (112)] and the y-axis represents the amplitude of the DAHCs. Note that the DAHMs and DAHCs shown here correspond, for each f, to the dominant pair of DAH eigenvalues located above f in Fig. 2(a).

Close modal
FIG. 5.

DAHMs and DAHCs from the traveling wave of Fig. 1(b). Same as in Fig. 4 but for the traveling wave of Fig. 1(b). The DAHMs and DAHCs shown here correspond also, for each f, to the dominant pair of DAH eigenvalues located above f in Fig. 3(a).

FIG. 5.

DAHMs and DAHCs from the traveling wave of Fig. 1(b). Same as in Fig. 4 but for the traveling wave of Fig. 1(b). The DAHMs and DAHCs shown here correspond also, for each f, to the dominant pair of DAH eigenvalues located above f in Fig. 3(a).

Close modal

Visualization of DAHMs: The eigenvectors of the grand matrix C (i.e., the DAHMs) are visualized by exploiting the representation (59). More precisely, given a frequency f of interest and a DAHM, Wj, associated with this frequency, we extract d DAHM snippets, E p j , each of length M , and we form the array in which each row for 1 p d is given by E p j ( s ) when s varies from 1 to M ; the index p referring to channels and s to (embedding) time. The resulting “spatio-temporal” array with d rows and M columns gives thus a natural way to visualize a d M -dimensional DAHM. Such a visualization of DAHMs is adopted for the case of modulated and pure traveling waves considered here, as shown in the 1st and 2nd columns of Figs. 4 and 5. Such a representation of DAHMs is not limited to this example and is adopted throughout this article.

By adopting this representation, the resulting DAHM pattern provides information on, e.g., any out-of-phase or strong collinearity that would exist between different channels composing the original datasets, at a given frequency. Note that due to the explicit form of the DAHM snippet given in (106), this information is encapsulated into the θ p j 's and each DAHM visualized this way thus provides quantitative information about the phase shifts between channels as associated with a DAH eigenvalue.

For the modulated-wave case, the mixture of traveling/standing-like wave patterns as well as the local small-scale oscillations are outstandingly well captured, at each frequency f 0 , by the pair of DAHMs associated with the pair of dominant eigenvalues; see first and second columns of Fig. 4. We emphasize indeed that this mixture of patterns present in the field q(x, t) shown in Fig. 1(a) is not only found for the (few) mode(s) associated with the DAHCs of largest variance(s), but rather is distributed across a broad frequency band of Fourier frequencies f.

Noticeably, these patterns scale with the frequency: at each frequency f, the pair of dominant DAHMs exhibits the same patterns that repeat with a 1 / f -temporal periodicity; see 1st and 2nd columns of Fig. 4. The same elemental spatial pattern is thus repeated for a pair of (dominant) DAHMs across the frequencies.

The same feature is observed in the case of a pure traveling wave shown in Fig. 1(b). Here again, the (dominant) DAHMs capture the dominant pattern displayed by the field q(x, t) across a broad range of frequencies; see 1st and 2nd columns of Fig. 5.

Such a scaling property (in time)—of the principal pattern contained in the original dataset—as revealed by the dominant pairs of DAHM across the frequencies is in sharp contrast with other modes that would be obtained by techniques relying on a decomposition of the variance.

Such methods applied to q(x, t) would obviously identify another type of leading modes that would capture most of the variance, without however displaying a scaled-version of the dominant pattern across the frequencies, in particular, for the modes that capture a small fraction of the variance. The third column of Figs. 4 and 5 shows to the contrary that even for pairs of DAHMs associated with DAHCs of small amplitude, these modes exhibit the dominant pattern, repeated at the associated frequency f; see the cases f = 0.024 and f = 0.048 in the 1st and 2nd columns of Figs. 4 or 5.

We comment now the DAH phase spectra shown in Figs. 2(b) and 3(b). In both cases, the most striking feature about these phase spectra is the distinction appearing in their shape by comparing the frequency band 0 < f < f c with f > f c . For f > f c , one can indeed distinguish that all the lines converge towards a single point as one approaches the Nyquist frequency fN.

By relying on the discussion conducted at the end of Sec. V, this convergence corresponds to a convergence of the η + p ( f ) towards a same value, expressing thus a form of phase synchronization between the channels of a DAHM as f f N , by recalling (80) and (60).

One could think that such a phase synchronization revealed by the DAH phase spectrum is limited to the case of pure or modulated traveling waves, but actually, several conducted experiments on other spatio-temporal fields issued from various dynamical systems (not shown) seem to indicate that this feature occurs for a broad set of deterministic dynamics, including chaotic ones; see also Sec. IX A 2. Roughly speaking, the poor resolution that comes with the Nyquist frequency explains that differences between the phases across the channels—subtle by nature—must vanish as f f N . However, here lies maybe the most useful aspect of the DAH phase spectrum. As we will see in Sec. IX B 2, dataset issued from a stochastic system does not share this convergence property of the phases across the channels of a DAHM, as f f N . The DAH phase spectrum appears thus to be potentially a useful diagnostic tool to distinguish between spatio-temporal deterministic or stochastic dynamics.

Let X ( t ) = ( X 1 ( t ) , , X d ( t ) ) be a d-dimensional time series from which the operator L defined in (52)(53) is built from the cross-correlation functions ρ i , j ( t ) , after (possible) periodization such as described in Sec. VI B. We adopt the notations of Theorem V.1 and form the following time-continuous version of the DAHC associated with an eigenfunction Wk which—in the continuous time setting adopted here—gives

ξ k ( t ) : = p = 1 d T X p ( t + s ) E p k ( s ) d s = p = 1 d T X p ( t + s ) B p k cos ( ω k s + θ p k ) d s = p = 1 d T T ( t ) X p ( s ) B p k cos ( ω k s + θ p k ) d s ,
(115)

i.e., the continuous counterpart of (112).

If Xp lies in D(A) given in (30) (with g = 1), then

d T ( t ) X p d t = A T ( t ) X p = T ( t ) A X p ,
(116)

see e.g., Ref. 39 (Chap. II, Lemma 1.3).

In particular,

d ξ k ( t ) d t = p = 1 d T T ( t ) d X p d s B p k cos ( ω k s + θ p k ) d s .
(117)

As for (49), one can prove due to (58) that to a pair of DAH eigenvalues corresponds a pair of DAHMs, for which their corresponding phases differ by π / 2 , i.e., in each DAHM pair, the modes are shifted by one fourth of the period.

As a consequence, the DAHC, ξ ̃ k , associated with the mode in phase quadrature with Wk is given by

d ξ ̃ k ( t ) d t = p = 1 d T T ( t ) d X p d s B p k sin ( ω k s + θ p k ) d s .
(118)

Forming the complex time series z k ( t ) = ξ k ( t ) + i ξ ̃ k ( t ) where i 2 = 1 , we obtain

d z k ( t ) d t = p = 1 d T T ( t ) d X p d s B p k exp ( i ( ω k s + θ p k ) ) d s .
(119)

An integration by parts gives

d z k ( t ) d t = i ω k p = 1 d T T ( t ) X p ( s ) B p k exp ( i ( ω k s + θ p k ) ) d s ,
(120)

that is

d z k ( t ) d t = i ω k z k ,
(121)

supplemented by the initial condition

z k ( 0 ) = p = 1 d T X p ( s ) B p k exp ( i ( ω k s + θ p k ) ) d s .
(122)

If Xp is not in the domain D(A) (such as for X ( t ) issued from a mixing dynamical system) and the base interval I for T is taken to be [ τ / 2 , τ / 2 ] , then the boundary terms appear and we are left (formally) with

d z k ( t ) d t = i ω k z k + p = 1 d B p k e i θ p k ( X p ( t + τ 2 ) X p ( t τ 2 ) ) .
(123)

An integration of (123) gives

z k ( t ) = z k ( 0 ) e i ω k t + p = 1 d B p k 0 t e ( i ω k ( t r ) + i θ p k ) ( X p ( r + τ 2 ) X p ( r τ 2 ) ) d r ,
(124)

and thus

ξ k ( t ) = Re ( z k ( 0 ) e i ω k t ) + p = 1 d B p k 0 t cos ( ω k ( t r ) + θ p k ) ( X p ( r + τ 2 ) X p ( r τ 2 ) ) d r .
(125)

The boundary terms introduce exogenous forcing from the time series components, themselves. Since the convolution term in the RHS of (125) involves a cosine function oscillating at the frequency ωk, the frequencies g ω k contained in the time series X p ( t + τ / 2 ) X p ( t τ / 2 ) are filtered out and the resulting DAHC, ξ k ( t ) , is narrowband about the frequency ωk. This latter property is shared by all the DAHCs associated with the frequency ωk. The resulting DAHC pair ( ξ k ( t ) , ξ ̃ k ( t ) ) exhibits thus narrowband oscillations (about the frequency ωk) that are superimposed to the pure oscillation Re ( z k ( 0 ) e i ω k t ) .

Although the RHS of (125) provides an exact representation (in the case of an infinite time series) of a DAHC, its usage in practice is limited since it requires the knowledge of the original dataset for future time instants, which is incompatible with, e.g., a prediction purpose. Furthermore, we desire a model of the ξ k ( t ) 's which avoids an explicit dependence on the data, i.e., the Xp's appearing in (125).

Since the oscillatory and modulated behavior of time series that are nearly in phase quadrature and carried out about one dominant frequency may be reasonably emulated by Stuart-Landau (SL) oscillators,63 we propose closures of Eq. (123) built from SL oscillators. By closure, we mean that we want to mimic the effects of the term

p = 1 d B p k e i θ p k ( X p ( t + τ 2 ) X p ( t τ 2 ) )

in Eq. (123), by means of equations involving only the z k ( t ) -variables plus stochastic ingredients to be determined. In Sec. VIII B below we propose a solution to this problem exploiting both the properties of DAH spectra (such as summarized by Theorem V.1) and ideas borrowed from the theory of multilayer stochastic models (MSMs) used for inverse modeling of time-evolving datasets.34 

We turn now to the use of DAHMs and DAHCs for deriving inverse stochastic-dynamic models. The purpose is to show that these objects allows for a reduction of the data-driven modeling effort to the learning of elemental MSMs34 stacked per frequency. These elemental models fall into the class of networks of linearly coupled Stuart-Landau oscillators36 that may include memory terms as described below; see also Ref. 64. In a certain sense, the goal is to show that, given a sequence (of possibly partial) observations issued from a dynamical model, the DAHCs allow for recasting these observations so that they can be modeled within a universal parametric family of simple stochastic models, provided that, as mentioned earlier, the decay of correlations are sufficiently well resolved.65 We refer to Sec. IX for examples supporting this statement, and describe hereafter the aforementioned elemental models.

Let (xj, yj) be a pair DAHCs such as shown, for instance, in the third column of Figs. 9 or 12. Except at the zero-frequency and as motivated in Sec. VIII A, DAHCs are typically oscillatory time series that are narrowband about a dominant frequency f and whose oscillations are modulated in time. Such a pair of time series can be reasonably expected to be modeled via nonlinear oscillating systems near a Hopf bifurcation to capture the frequency f, and subject to noise, to capture the amplitude modulations; see Ref. 63. We assume that we have a collection of d such pairs, each associated with the same frequency f and exhibiting modulations occurring in different time locations.

If the signal is made of a concatenation of d different time series such as for the L96 flow, such a pair has a simple interpretation: it provides the manifestation, in the time domain, of the frequency f contained in the system's flow as captured by the corresponding pair of DAHMs.

For a narrowband pair of DAHCs, ( x j ( t ) , y j ( t ) ) , given at a frequency f 0 , we envision thus its modeling by some stochastic perturbation of a Stuart-Landau model, namely, by

z ̇ = ( μ + i γ ) z ( 1 + i β ) | z | 2 z + ϵ t , z ,
(126)

where z ( t ) = x j ( t ) + i y j ( t ) , and μ , γ , and β are real parameters to be estimated. The system is here driven by some “reddish noise,” ϵ t = ( ϵ j x , ϵ j y ) , to be estimated, also, from the time history of z(t). The issue is that the pairs, when modeled this way, would fail in reproducing the phase coherences across the channels, that would be contained in the signal such as in the case of a wave propagation, for instance. The estimation and simulation of ϵt requires thus a special care due to the presence of typical non-zero cross correlations among the channels.

The MSM framework34 allows us to address this issue, while keeping the advantages of the modeling framework provided by Eq. (126). It allows, as we will see in applications (see Sec. IX), for providing stochastic-dynamic inverse models able to respect the phase coherence of the collective behavior of the DAHC pairs ( x j ( t ) , y j ( t ) ) at a given frequency f 0 . In that respect, one proposes to model each such a DAHC pair by the following type of MSMs, which in the case of two layers used to model the main noise residual, ϵt, can be written as the following system of SDEs:

x j ̇ = β j ( f ) x j α j ( f ) y j + σ j ( f ) x j ( x j 2 + y j 2 ) + b i j x i j i J d ( f ) ( f ) x i + a i j x i j i J d ( f ) ( f ) y i + ϵ j x , y j ̇ = α j ( f ) x j + β j ( f ) y j + σ j ( f ) y j ( x j 2 + y j 2 ) + a i j y i j i J d ( f ) ( f ) x i + b i j y i j i J d ( f ) ( f ) y i + ϵ j y , ϵ j x ̇ = K 11 j ( f ) x j + K 12 j ( f ) y j + M 11 j ( f ) ϵ j x + M 12 j ( f ) ϵ j y + ζ j x , ϵ j y ̇ = K 21 j ( f ) x j + K 22 j ( f ) y j + M 21 j ( f ) ϵ j x + M 22 j ( f ) ϵ j y + ζ j y , ζ j x ̇ = L 11 j ( f ) x j + L 12 j ( f ) y j + P 11 j ( f ) ϵ j x + P 12 j ( f ) ϵ j y + N 11 j ( f ) ζ j x + N 12 j ( f ) ζ j y + Q 11 j ( f ) B ̇ 1 j + Q 12 j ( f ) B ̇ 2 j + i j i J d ( f ) k = 1 2 Q 1 k i ( f ) B ̇ k i , ζ j y ̇ = L 21 j ( f ) x j + L 22 j ( f ) y j + P 21 j ( f ) ϵ j x + P 22 j ( f ) ϵ j y + N 21 j ( f ) ζ j x + N 22 j ( f ) ζ j y + Q 21 j ( f ) B ̇ 1 j + Q 22 j ( f ) B ̇ 2 j + i j i J d ( f ) k = 1 2 Q 2 k i ( f ) B ̇ k i .
(MSLM)

Such a model is named multilayer stochastic Stuart-Landau models (MSLM). In (MSLM), the set J d ( f ) denotes the subset of indices i in J ( f ) [see (108)] that correspond to d distinct pairs of DAH eigenvalues associated with the frequency f 0 . Recall that the existence of such pairs when f 0 is ensured by Theorem V.1 and, in particular, (58).

The B k i 's with k in {1, 2} and i in J d ( f ) form 2d independent Brownian motions. At a given frequency f, the d pairs are linearly coupled as indicated by the terms in the summation terms appearing in the xj- and yj-equations. As depicted in the schematic of Fig. 6, the main and last layers (here the xj-, yj-, and ζ j x / y -equations) are the only layers in which these pairs are explicitly coupled. Without these coupling terms, we are left with a collection of d uncoupled normal forms of a Hopf bifurcation perturbed by a “reddish” noise, which are the analogue in real domain of the aforementioned Stuart-Landau model (126). The linear coupling at the main level is here aimed to enhance some coherence between the stochastic nonlinear oscillations.

FIG. 6.

Schematic of Eq. (MSLM). Each circle represents either a DAHC pair (xj, yj) or the corresponding pairs of (successive) noise residuals present in Eq. (MSLM). The vertical or horizontal dashed lines connecting two pairs indicate the existence of a mutual dependence between these pairs in the formulation of Eq. (MSLM). Note that these mutual dependences are only considered here for the pairs associated with a same frequency f, i.e., for j in J ( f ) given by (108).

FIG. 6.

Schematic of Eq. (MSLM). Each circle represents either a DAHC pair (xj, yj) or the corresponding pairs of (successive) noise residuals present in Eq. (MSLM). The vertical or horizontal dashed lines connecting two pairs indicate the existence of a mutual dependence between these pairs in the formulation of Eq. (MSLM). Note that these mutual dependences are only considered here for the pairs associated with a same frequency f, i.e., for j in J ( f ) given by (108).

Close modal

The noise residual ( ϵ j x , ϵ j y ) in the main layer of (MSLM) is modeled in the intermediate layer by means of linear dependences involving only ( ϵ j x , ϵ j y ) on one hand, and the j-th pair (xj, yj) on the other. The noise residual ( ζ j x , ζ j y ) of that intermediate layer is then modeled still via linear dependences that now include this noise itself, the noise residual of the previous layer, and the DAHC pair (xj, yj). The terms in the summations appearing in the ζ j x - and ζ j y -equations of the last layer involve the other pairs as mentioned earlier. The presence of such terms is aimed to take into account—at the noise level—of cross-correlations between the d DAHC pairs corresponding to a same frequency. In practice, the corresponding coefficient Q l k j ( l , k { 1 , 2 } , j { 1 , , d } ) are estimated from the vector of size 2d formed after the concatenation of the last level residuals for each pair.

Note that for the frequency f = 0, and as mentioned in Sec. VI D, there are exactly d DAHMs that are not paired, and they are modeled by a linear MSM, following Ref. 34; see also Ref. 66 [Eq. (9)].

An MSLM is not limited to the two-layer case shown earlier for the modeling of the main noise residual ( ϵ j x , ϵ j y ) . Following Ref. 34, extra layers may be added, with intermediate layers that only involve a linear dependence on the j-th DAHC pair (xj, yj) and on the corresponding noise residuals from the previous layers. In other words, the extra layers must respect the overall structure of (MSLM), with coupling between the other pairs (associated with the same frequency f) only present in the main and last layers. The addition of an extra layer would thus correspond to the addition of another intermediate plane in the schematic of Fig. 6, with only vertical connections to the rest of the structure.

The procedure of adding layers is controlled by a stopping criterion as described in Ref. 34 (Appendix A) which allows thus, in principle, for a different number of layers per frequency. In general, the layers with their surrogate variables, here ( ϵ j x , ϵ j y ) and ( ζ j x , ζ j y ) , allow for a learning of the temporal correlations contained in the main noise residual ϵt, as well as of its dependences on the history of ( x j ( t ) , y j ( t ) ) often required for an appropriate modeling of ϵt; see Ref. 34 for a general discussion on this topic. Typically, the more layers, the more complex are these correlations and memory dependences; see Ref. 34 (Proposition 3.3), Ref. 34 (Sec. 7), and Ref. 66 for examples. Sections IX A 3 and IX B 3 also illustrate this point.

As emphasized in (MSLM), the coefficients to be learned depend on the frequency f. An MSLM seeks for the modeling of the DAHC pairs ( x j , y j ) j J ( f ) associated with the frequency f, such as extracted from the original (or preprocessed67) dataset by application of the DAH decomposition. Once this extraction is computed, an MSLM can be learned in parallel for each frequency of interest, by application of, e.g., simple successive regressions. In practice, linear constraints must be incorporated in these regressions such as σ j ( f ) 0 to ensure stability, as well as antisymmetry, to estimate the coefficients β j ( f ) and α j ( f ) .

Once all the resulting (few) MSLM coefficients have been estimated, the parallelization advantages of the MSLMs across the frequencies rely also on its simulation: no need of extra coupling across the frequencies is indeed needed other than running the MSLMs by the same noise realization, for each frequency. The simulation procedure can be thus also parallelized at the convenience of the modeler to obtain the collective behavior of the simulated DAHCs. These simulated DAHCs can be then multiplied with the corresponding DAHMs such as described in the  Appendix, to obtain various possible frequency-band modeling of the original dataset.

Of course, one condition to the success of the resulting modeling approach relies on the ability of the DAH decomposition in extracting modulated DAHC pairs that come in nearly phase quadrature, and that are narrowband in the frequency domain, an MSLM being naturally dedicated to the modeling of such time series. The tradeoff between a good resolution of the decay of temporal correlations and the temporal length of the available dataset, seem to play a key role in that respect. Various conducted experiments (not shown) have indeed corroborated the idea that the DAHC pairs are expected to be narrowband and to come in nearly phase quadrature when the decay of temporal correlations has been sufficiently resolved.

In what follows, the combination of (i) the extraction of DAHCs by the DAH method of Sec. VI D and (ii) their modeling according the the MSLM approach described earlier is named a DAH-MSLM modeling. We present Sec. IX below applications to examples for which this criterion is met and a successful DAH-MSLM modeling is obtained in each case.

We apply hereafter the DAH-MSLM inverse modeling approach to two complex flows, obtained, respectively, from a chaotic ODE system and an stochastic partial differential equation (SPDE), with a linear drift.

1. The dynamical model

In this section, we apply the DAH approach to the so-called Lorenz 96 (L96) model.68 The evolution equations of the L96 model, in its simplified version,69 is written as

d X k d t = X k 1 ( X k + 1 X k 2 ) X k + F 1 , k = 1 , , K ,
(127)

supplemented with the boundary conditions

X k K = X k + K = X k .
(128)

In the absence of forcing and dissipation, the sum of the squares of the variables (the energy of the system) is conserved. This model has become popular as a toy model of the atmosphere to test various methods in weather and climate sciences as displaying a large set of chaotic behaviors with various mixing properties, as K and F1 are varied.

2. DAH spectra, DAHMs, and DAHCs

The model is set to a standard chaotic regime corresponding to K = 40 and F 1 = 6 , and is integrated by using a 4th-order Runge-Kutta scheme with a time step δ t = 5 × 10 3 . Excluding transient behavior and subsampling every 10th time step, we store N = 11 900 time instances of the K-dimensional vector X ( t ) = ( X 1 ( t ) , , X K ( t ) ) .

As shown in Fig. 7(a), the L96 flow is dominated by irregular spatio-temporal traveling wave-like patterns. We estimated the cross-correlations ρ i , j ( t ) between the variables X i ( t ) and X j ( t ) for i, j in { 1 , , 40 } , and we formed the ( 2 M 1 ) × ( 2 M 1 ) -matrix C given in (104), with M = 160. The corresponding absolute values of its eigenvalues (ranked per frequency)—forming thus the corresponding DAH power spectrum—are shown in Fig. 8(a). The corresponding DAH phase spectrum is shown in Fig. 8(b).

FIG. 7.

Spatio-temporal fields. Panel (a) shows a time evolution the Lorenz 96's flow solving Eqs. (127) and (128), and panel (b) shows the solution field of the stochastic heat equation (129). See text for the parameter values and the numerical schemes. (a) L96 model and (b) Stochastic heat equation.

FIG. 7.

Spatio-temporal fields. Panel (a) shows a time evolution the Lorenz 96's flow solving Eqs. (127) and (128), and panel (b) shows the solution field of the stochastic heat equation (129). See text for the parameter values and the numerical schemes. (a) L96 model and (b) Stochastic heat equation.

Close modal
FIG. 8.

Multidimensional DAH power and phase spectra for the L96 flow. We chose M = 160 to estimate the cross correlations and form the corresponding the matrix C given in (104). The DAH power spectrum and DAH phase spectrum are obtained as described in Sec. VI D. Each non-zero frequency is associated with 40 eigenpairs that correspond to the total of d = K = 40 variables in the L96 model. (a) DAH power spectrum and (b) DAH phase spectrum.

FIG. 8.

Multidimensional DAH power and phase spectra for the L96 flow. We chose M = 160 to estimate the cross correlations and form the corresponding the matrix C given in (104). The DAH power spectrum and DAH phase spectrum are obtained as described in Sec. VI D. Each non-zero frequency is associated with 40 eigenpairs that correspond to the total of d = K = 40 variables in the L96 model. (a) DAH power spectrum and (b) DAH phase spectrum.

Close modal

Interestingly, these spectra as well the underlying DAHMs share some common features with those displayed in the case of a pure traveling wave shown in Secs. VII A and VII B. First, the DAHMs shown in the 1st and 2nd columns of Fig. 9 share clearly common patterns with those shown in the respective columns of Fig. 5.

FIG. 9.

DAHMs and DAHCs from the L96 flow. Same as in Fig. 4 but for the L96 flow shown in Fig. 7(a). Here, M = 2 M 1 (here M = 160), and d = 40. The DAHMs and DAHCs shown here correspond also, for each frequency f, to the dominant pair of DAH eigenvalues located above f in Fig. 8(a).

FIG. 9.

DAHMs and DAHCs from the L96 flow. Same as in Fig. 4 but for the L96 flow shown in Fig. 7(a). Here, M = 2 M 1 (here M = 160), and d = 40. The DAHMs and DAHCs shown here correspond also, for each frequency f, to the dominant pair of DAH eigenvalues located above f in Fig. 8(a).

Close modal

Second, a phase synchronization between the channels of a DAHM as f approaches the Nyquist frequency is also observed in Fig. 8(b) as in Fig. 3(b).

Third, this synchronization is in contrast with the disordered distribution of phases shown in a frequency band 0 < f < f c with f c 0.1 ; the frequency band for which the most energetic pair of DAH eigenvalues (for the quadratic form Q δ ) are located; see Fig. 8(a).

Nevertheless, the DAH power spectrum exhibits much broader peaks in the case of the L96 flow than in the case of a pure traveling wave; compare Fig. 8(a) with Fig. 3(a).

Finally, the most noticeable difference is shown by the DAHCs. Whereas the latter were all periodic in the case of traveling wave as expected [see (121)], they are—as predicted by (125)70—clearly modulated in amplitude and narrowband in frequency in the case of the L96 flow; compare the 3rd column of Fig. 5 with that of Fig. 9. Furthermore, the pairs of DAHCs are typically constituted by time series (for f 0 ) that are nearly in phase quadrature; see the panels corresponding to f = 0.016, f = 0.031, and f = 0.047 in the 3rd column of Fig. 12. We are thus in a situation described in Sec. VIII B. Section IX A 3 below addresses the inverse-modeling problem of the L96 flow, following the approach described in Sec. VIII B.

3. DAH-MSLM modeling of the Lorenz 96 model

In this section, we report on the inverse modeling skills obtained by application of the DAH-MSLM approach described in Sec. VIII. For that purpose, the original field obtained by integrating Eq. (127), as mentioned earlier, has been sampled every 10 time steps in time. The frequency band that is modeled is from f = 0 up to f = 6.29 × 10 2 , which represents a fraction of the variance of 86 % , and 20 Fourier frequencies according to the resolution set by M = 2 × 160 1 = 319 .

The corresponding DAHCs are modeled for each of these frequencies according to an MSLM of the form described in Sec. VIII B, in which the number of layers per frequency is 11, including the main one. After multiplication by the corresponding DAHMs, the simulated DAHCs allow us, as described in the  Appendix, to form a simulated DAH-MSLM field. The resulting field is shown in Fig. 10 and exhibits striking common features with the original field, shown in the same figure for a better comparison.

FIG. 10.

The L96 flow and its simulated DAH-MSLM field. The L96 flow shown in the top panel corresponds to the sampling units used at the learning stage, i.e., every 10th time steps; see text for more details.

FIG. 10.

The L96 flow and its simulated DAH-MSLM field. The L96 flow shown in the top panel corresponds to the sampling units used at the learning stage, i.e., every 10th time steps; see text for more details.

Close modal

Interested by the ability of the simulated DAH-MSLM field to reproduce the time-variability of the L96 flow, we computed for each of these fields the autocorrelation function (ACF) (in time) of each channel, followed by an averaging of the resulting individual ACFs over the channels. The result is shown in Fig. 15(a) and as one can see the ability of the simulated DAH-MSLM field in reproducing the time-variability of the L96 flow is very satisfactory, according to this metric.

It should be emphasized that beyond the embedding window size, here M = 319 , these reproduction skills of the (mean) ACF, deteriorate due in particular, to the periodization involved in the determination of the matrix C , as in any other method exploiting time-embedding techniques.

So far we have been dealing in applications with deterministic systems. The Sec. IX B illustrates that striking inverse modeling skills can also be achieved by applying the DAH-MSLM approach for the inverse modeling of stochastic systems, still with the same elemental MSLMs as described in Sec. VIII B.

1. The dynamical model

The dataset analyzed hereafter is produced from a simple stochastic partial differential equation (SPDE), namely, the following stochastic heat equation posed on the interval [ 0 , L ] (L > 0) with periodic boundary conditions:

d u ( t , x ) = ( λ u ( t , x ) + D x x 2 u ( t , x ) ) d t + 2 λ d W ( t , x ) , x T L = / ( L ) .
(129)

Here, λ and D are positive parameters, and d W denotes a “space-time” white noise process on the circle T L , roughly speaking a Gaussian random field with spatio-temporal correlations for t , s 0 , given by

E ( d W ( t , x ) d W ( s , y ) ) = δ 0 ( t s ) δ 0 ( y x ) , x , y T L .

There are (at least) two approaches to frame mathematically the notion of space-time white noise and related concepts of solutions; see, e.g., Ref. 71 for a comparison. Both, the random-field approach based on Green functions72 and the functional analysis approach,73,74 give the same notion of weak/mild solutions, and stochastic heat equations such as Eq. (130) serve as a fundamental object to study various properties or formulas about such solutions; see e.g., Refs. 75–78.

The functional analysis approach teaches us that the stochastic process u ( t , · ) solving Eq. (129) is actually an infinite-dimensional Ornstein-Uhlenbeck process (OU) in the Hilbert space H = L 2 ( T L ) ; see Ref. 75. Our goal is to analyze the DAH power and phase spectra of the corresponding (random) spatio-temporal field, u(t, x), and see whether these spectra provide an intuitive characterization of such an OU process.

For that purpose, we first simulate Eq. (129) as follows. Over a grid of mesh size δ x = L / N x , the discrete approximation u j n of u ( n δ t , j δ x ) is obtained by finite differences for which the noise term is approximated at each time step n δ t (via an explicit Euler-Maruyama scheme) by ξ n δ t , where ( ξ j n ) j = 1 N x denotes an Nx-dimensional vector of random variables each drawn independently with respect to 1 j N x and n (except ξ 1 n = ξ N x n ) from the standard normal distribution N ( 0 , 1 ) .

We chose N x = 2 7 , L = 2 π , D = 2 × 10 1 , and δ t = 10 3 for our simulations. The initial data are taken to be

u 0 ( x ) = 10 1 2 L cos ( 2 π x L ) .
(130)

Without any surprise, the integration reveals striking differences in terms of patterns exhibited by the field u(t, x) compared to those obtained from the nonlinear example of Sec. IX A; compare panel (b) with panel (a) in Fig. 7.

2. DAH spectra, DAHMs, and DAHCs

Excluding transient behavior and subsampling every 8 δ t , we store N = 37 500 time instances of the d-dimensional vector X ( t ) = ( X 1 ( t ) , , X d ( t ) ) , in which the Xp's are also obtained by subsampling every 4 δ x the original field ( u j n ) , i.e., here d = 32. The resulting coarse-grained SPDE field is shown in the top panel of Fig. 13. We built up the matrix C (104) whose Hankel blocks are given in (103) with M = 80, and 1 i , j d . Diagonalizing C gives the multidimensional DAH spectra as explained in Sec. VI D. The corresponding DAH power (resp. DAH phase) spectrum for the coarse-grained SPDE is shown in Fig. 11(a) [resp. Fig. 11(b)].

FIG. 11.

Multidimensional DAH power and phase spectra of the SPDE flow. We chose M = 80 to estimate the cross correlations and form the corresponding matrix C given in (104). The DAH power spectrum and DAH phase spectrum are obtained as described in Sec. VI D. Each non-zero frequency is associated with 32 eigenpairs that correspond to the total of d = 32 spatial channels obtained by coarse-graining the original dataset as described in text. (a) DAH power spectrum and (b) DAH phase spectrum.

FIG. 11.

Multidimensional DAH power and phase spectra of the SPDE flow. We chose M = 80 to estimate the cross correlations and form the corresponding matrix C given in (104). The DAH power spectrum and DAH phase spectrum are obtained as described in Sec. VI D. Each non-zero frequency is associated with 32 eigenpairs that correspond to the total of d = 32 spatial channels obtained by coarse-graining the original dataset as described in text. (a) DAH power spectrum and (b) DAH phase spectrum.

Close modal

As it can be observed by comparison with Fig. 8(a), the DAH power spectrum shown in Fig. 11(a) does not exhibit any broadband peak, and is actually composed by curves of DAH eigenvalues that each exhibits a decay from low to high frequencies, reminiscent with that of a scalar OU process. What is remarkable is that the DAH power spectrum provides here a multidimensional picture of what one should expect for an OU process obtained from a standard 1D Langevin equation, showing thus here again the relevance of the quadratic form Q δ , as an energy characterizing the distribution of eigenvalues of C at a given frequency; see (105) and Remark V.1–(ii).

As another noticeable feature, the phase synchronization such as that observed on the DAH phase spectra of the L96 flow and the traveling wave examples is no longer present for this stochastic example, and is instead replaced in Fig. 11(b) by the presence of some (vaguely visible) bands within a diffuse background. Going back to the interpretation of the phase spectrum discussed at the end of Sec. V in relation with (80), the presence of such a diffuse background, even for f approaching the Nyquist frequency fN, can be attributed to the presence of non-vanishing cross-spectral terms ρ p , q ̂ ( f ) leading to a non-negligible R p ( f ) given in (78). This feature is consistent with the intuitive idea that the noise-term in Eq. (130) is white and thus forces evenly the resolved frequency band. It should be emphasized that such a feature has been also observed in other conducted experiments on datasets issued from various stochastic systems driven by a white and additive noise.

Finally, each pair of DAHCs associated with a frequency f (with f 0 ) is—as for the nonlinear example of Sec. IX A—constituted by time series modulated in amplitude, narrowband about the frequency f, that are nearly in phase quadrature; see the panels corresponding to f = 0.013, f = 0.025, and f = 0.038 in the 3rd column of Fig. 12. We are thus here again in a favorable situation to apply to these DAHCs the MSLM modeling framework of Sec. VIII B. Section IX B 3 below examines examine the corresponding modeling skills.

FIG. 12.

DAHMs and DAHCs from the SPDE flow. Same as in Fig. 4 but for the SPDE flow shown in Fig. 7(b). Here, M = 2 M 1 (here M = 80), and d = 32. The DAHMs and DAHCs shown here correspond also, for each frequency f, to the dominant pair of DAH eigenvalues located above f in Fig. 11(a).

FIG. 12.

DAHMs and DAHCs from the SPDE flow. Same as in Fig. 4 but for the SPDE flow shown in Fig. 7(b). Here, M = 2 M 1 (here M = 80), and d = 32. The DAHMs and DAHCs shown here correspond also, for each frequency f, to the dominant pair of DAH eigenvalues located above f in Fig. 11(a).

Close modal

3. Coarse-grained DAH-MSLM modeling

Recall that the original SPDE field has been sampled every 8-th time step, and every 4-th spatial increment δx in space; the resulting coarse-grained field is shown on the top panel of Fig. 13. The frequency band that is modeled is from f = 0 up to f = 6.33 × 10 2 , which represents a fraction of the variance of 86 % , and corresponds to 10 first Fourier frequencies according to the resolution set by M = 2 × 80 1 = 159 .

FIG. 13.

The stochastic heat flow and its simulated DAH-MSLM field. The coarse-grained SPDE field shown in the top panel corresponds to that obtained from sampling units used for the learning stage, i.e., every 8th time steps and 4th spatial mesh size, and after projection onto the pairs of DAHMs corresponding to the first 10 Fourier frequencies; see text for more details.

FIG. 13.

The stochastic heat flow and its simulated DAH-MSLM field. The coarse-grained SPDE field shown in the top panel corresponds to that obtained from sampling units used for the learning stage, i.e., every 8th time steps and 4th spatial mesh size, and after projection onto the pairs of DAHMs corresponding to the first 10 Fourier frequencies; see text for more details.

Close modal

The corresponding DAHCs are here again modeled for each of these frequencies according to an MSLM of the form described in Sec. VIII B, in which the number of layers per frequency is 4, including the main one. The resulting simulated DAH-MSLM field is shown in Fig. 13 and exhibits striking common features with the coarse-grained SPDE field shown in the same figure.

At a statistical level, the simulated DAH-MSLM field (red curve in Fig. 14) captures remarkably well the distribution of energy across the spatial scales when compared with the coarse-grained version of the SPDE field (blue curve in Fig. 14) from which the parameters of the DAH-MSLM have been estimated. Note that the blue curve in Fig. 14 shows this energy distribution for the coarse-grained SPDE field as projected onto the first pairs of DAHMs corresponding to the first 10 Fourier frequencies. The black curve shows in the same figure, the energy distribution for the full coarse-grained SPDE field.

FIG. 14.

Energy spectrum. Distribution of energy across the spatial scales. The blue curve shows this energy distribution for the coarse-grained SPDE field as projected onto the first pairs of DAHMs corresponding to the first 10 Fourier frequencies; see text for more details. The black curve shows this energy distribution for the full coarse-grained SPDE field (i.e., without projection).

FIG. 14.

Energy spectrum. Distribution of energy across the spatial scales. The blue curve shows this energy distribution for the coarse-grained SPDE field as projected onto the first pairs of DAHMs corresponding to the first 10 Fourier frequencies; see text for more details. The black curve shows this energy distribution for the full coarse-grained SPDE field (i.e., without projection).

Close modal

Finally, as in the case of the L96 flow, in order to assess the ability of the simulated DAH-MSLM field to reproduce the time-variability contained in the coarse-grained SPDE field, we computed for each of these fields the ACFs of each channel, followed by an averaging of the resulting individual ACFs over the channels. The result is shown in Fig. 15(b) and as one can see the ability of the simulated DAH-MSLM field in reproducing the (average) decay of correlations is very satisfactory. Noteworthy is the nearly exponential decay of the later, as intuitively expected from an infinite-dimensional OU process.

FIG. 15.

Averaged autocorrelation function (ACF): modeling skills. In each panel the black curve is obtained by computing the ACFs of each channel of the corresponding dataset, followed by an averaging of the resulting individual ACFs over the channels. The red curve is obtained by averaging furthermore over the realizations of a given DAH-MSLM simulation. (a) DAH-MSLM vs Lorenz 96 and (b) DAH-MSLM vs SPDE.

FIG. 15.

Averaged autocorrelation function (ACF): modeling skills. In each panel the black curve is obtained by computing the ACFs of each channel of the corresponding dataset, followed by an averaging of the resulting individual ACFs over the channels. The red curve is obtained by averaging furthermore over the realizations of a given DAH-MSLM simulation. (a) DAH-MSLM vs Lorenz 96 and (b) DAH-MSLM vs SPDE.

Close modal

To summarize, the DAH-MSLM inverse modeling framework provides here again efficient emulators of the observed dynamics for various metrics. The most striking feature of the approach is that either for the linear stochastic example examined here or the nonlinear chaotic one of Sec. IX A, the predictor functions stay the same, namely, the monomials used for the class of MSLMs described in Sec. VIII B. The data-adaptive character of the DAHMs manifested by the capture of the phase information contained79 in the multivariate signal to be modeled constitutes a key element for a successful modeling. It allows indeed for breaking the signal into elementary bricks, here the DAHCs, that fall naturally within the MSLM class of models. In other words, the DAHMs provide a change of basis which transforms systematically the original dataset into time series generated by a same class of SDEs, provided the corresponding DAHCs are modulated in amplitude and narrowband (for f 0 ). Of course many questions arise, as it did for many novel mathematical concepts and tools, when first introduced. Section X formulates some of these.

Thus, by means of rigorous harmonic analysis tools, this article provides a natural framework for the representation of complex time-evolving fields by a combination of narrowband time series, the DAHCs, suitable for their modeling within a universal class of models, namely, the MSLMs. The framework opens up several possible directions for future research. We outline some of these directions below and comment on recent achievements obtained on non-synthetic datasets.

  1. The modeling approach described in this article may benefit from several natural developments. First, terms of the form γ j ( f ) y j ( x j 2 + y j 2 ) and γ j ( f ) x j ( x j 2 + y j 2 ) in the xj- and yj-equations, respectively, could be incorporated within the formulation of an MSLM. The introduction of such nonlinear “twist” terms is in part motivated by the analysis conducted in Ref. 63. There, it has been indeed shown that such terms are responsible of certain interactions between the noise and the nonlinear terms (associated with hypoellipticity) that may be important to resolve. The inclusion of such twist terms may be thus relevant for the proper modeling of a DAHC pair, for certain applications.

    Also, the driving noise of an MSLM may not be required to belong to the class of correlated white noise processes80 such as in Eq. (MSLM). For instance, an MSLM driven by Lévy α-stable processes could be relevant for certain applications such as the modeling of precipitations in climate81,82 and other climate fields.83 

  2. The use of MSLMs may not be limited to the modeling of pairs of DAHCs as extracted by DAH. For instance, other multivariate time series decomposition techniques such as MSSA84 extract analogue time series known as ST-PCs, which when coming in pairs that are nearly in phase quadrature, could also be efficiently modeled by an MSLM, provided the corresponding pairs are modulated in amplitude while narrowband. We mention that a similar approach was adopted in Ref. 66 for the stochastic modeling of decadal variability in ocean gyres. There, however, the nonlinear coupling within a frequency bin was absent while linear couplings across the frequencies (falling within a targeted frequency band) were used to model the ST-PCs from the dataset considered in Ref. 66.

  3. Similarly, coupling across frequencies that go beyond the driving by the same noise realization may also be useful to build up MSLMs, in certain cases. The inclusion into an MSLM of the contribution of another frequency g than a frequency f may be indeed easily achieved by passing the corresponding terms associated with g into (MSLM). Of course, as far as the pairs of DAHCs are sufficiently narrowband, such a coupling may be superfluous, and in such a case, a “same noise realization coupling”—cheaper in terms of coefficients to be estimated85—may be amply sufficient to reach satisfactory modeling skills as shown for the time-evolving fields of Sec. IX.

  4. Unlike the synthetic datasets of Sec. IX, the DAH-MSLM approach has been applied to various geophysical multivariate datasets exhibiting complex patterns, from “real-world applications.” In that respect, the DAH-MSLM modeling of Arctic sea ice concentration has been successfully addressed in Ref. 86, providing, in particular, stable inverse models for much longer time intervals than their training period (1979–2014), beyond the end of the 21st century. Interestingly, these models exhibit an interdecadal variability consistent with proxy historical records, while reproducing accurately the (nonlinear) trend contained in the observations; see Ref. 86 (Fig. 13).

Prediction experiments using the DAH-MSLM framework have also been recently conducted in Ref. 87 for the Arctic sea ice extent (SIE). Here, the forecasts, as it would be the case for any other dataset, are made for the DAHCs using the best MSLMs learned during the training period. As a result, retrospective forecasts show that the resulting DAH-MSLM modeling is quite skillful in predicting the September SIE; i.e., for the most challenging SIE value to forecast. Noteworthy are the real-time forecasts conducted in 2016 which have shown very competitive skills as documented here: https://www.arcus.org/sipn/sea-ice-outlook/2016/post-season; see Ref. 87 for more details.

Finally, we mention that the DAH-MSLM framework has been shown to be useful for the data-driven study of the coupling between the solar wind and the magnetosphere as discussed in Ref. 88. This example shows, in particular, the relevance of the DAH-MSLM framework to provide inverse models of externally forced systems, another important theme for geophysical applications.89,90

The authors would like to thank the anonymous reviewers for their useful comments that helped improve the presentation of the manuscript.

The authors also thank Georg Gottwald for the inspiring discussions at the start of this project. The authors are grateful to Henk Dijstra, David Neelin, and Alexis Tantet for fruitful and inspiring collaborations about the stochastic Hopf bifurcation. Constructive comments by Pavel Berloff, Sergey Kravtsov, Honghu Liu, and Eric Simonnet in early stages of this project are also most appreciated. The authors are also grateful to James McWilliams for several inspiring discussions and constructive comments about this work. Finally, Michael Ghil has the full gratitude of the authors for the numerous stimulating discussions over the years as well as for his unflinching support.

This research was supported by Grant No. N00014-16-1-2073 from the Multidisciplinary University Research Initiative (MURI) of the Office of Naval Research, and by the National Science Foundation Grant Nos. OCE-1243175, OCE-1658357, and DMS-1616981.

Let us consider a d-dimensional time series X ( t ) = ( X 1 ( t ) , , X d ( t ) ) from which the cross-correlation coefficients ρ i , j between the i-th and the j-th channels have been estimated. Once the DAHMs have been determined as eigenvectors of the matrix C defined in (104) (see Sec. VI D) and the DAHCs according to (112), one can use them to extract from X ( t ) the signal that corresponds to a certain frequency band [ f s , f e ] . We describe hereafter how to proceed to do so.

First, given a M -dimensional DAHM snippet, E p j ( s ) ( 1 s M ), and its corresponding DAHC one, t ξ j ( t ) , form the reconstructed component (RC), R p j ( t ) , given at time t by

R p j ( t ) = 1 M t s = L t U t ξ j ( t s + 1 ) E p j ( s ) , 1 s M .
(A1)

The values of the normalization factor Mt, as well as of the lower and upper bound of summation Lt and Ut, differ between the central part of the time series and its end points, and are given here as for M-SSA; see Ref. 84 [Eq. (12)]. Forming the corresponding reconstructed multivariate time series R j ( t ) = ( R 1 j ( t ) , , R d j ( t ) ) , and summing over all of the (resolved) j's allows us to recover the original multivariate time series X ( t ) . A summation over a subset of { 1 , , d M } gives a partial reconstruction of X ( t ) .

Unlike for M-SSA84 or PCA, such a partial reconstruction is not necessarily guided by a certain targeted percentage of the variance of X ( t ) to be captured, but can also be performed over a targeted frequency band. Due to its harmonic flavor, the DAH approach allows us indeed to consider a notion of harmonic reconstructed component (HRC) which consists of—for a given channel p—summing the RCs corresponding to a given frequency f , namely,

R k ( t ; f ) : = j J ( f ) R p j ( t ) ,
(A2)

where J ( f ) corresponds to the set of indices given by (108).

In applications, the DAH approach thus offers a way to determine the contribution within a component X p ( t ) [of the multivariate signal X ( t ) ] of a particular frequency band I f : = [ f s , f e ] by simply adding up the HRCs over the frequencies f lying within If.

1.
Y.
Bengio
,
A.
Courville
, and
P.
Vincent
,
IEEE Trans. Pattern Anal. Mach. Intell.
35
,
1798
(
2013
).
2.
Y.
LeCun
,
Y.
Bengio
, and
G.
Hinton
,
Nature
521
,
436
(
2015
).
3.
I.
Jolliffe
,
Principal Component Analysis
(
Wiley Online Library
,
2002
).
4.
M. E.
Tipping
and
C. M.
Bishop
,
J. R. Stat. Soc.: Ser. B
61
,
611
(
1999
).
5.
B.
Schölkopf
,
A.
Smola
, and
K.-R.
Müller
,
Neural Comput.
10
,
1299
(
1998
).
6.
D.
Mukhin
,
A.
Gavrilov
,
A.
Feigin
,
E.
Loskutov
, and
J.
Kurths
,
Sci. Rep.
5
,
15510
(
2015
).
7.
A.
Gavrilov
,
D.
Mukhin
,
E.
Loskutov
,
E.
Volodin
,
A.
Feigin
, and
J.
Kurths
,
Chaos
26
,
123101
(
2016
).
8.
R. R.
Coifman
and
S.
Lafon
,
Appl. Comput. Harmonic Anal.
21
,
5
(
2006
).
9.
R. R.
Coifman
,
I. G.
Kevrekidis
,
S.
Lafon
,
M.
Maggioni
, and
B.
Nadler
,
Multiscale Model. Simul.
7
,
842
(
2008
).
10.
D.
Giannakis
and
A. J.
Majda
,
Stat. Anal. Data Min.: ASA Data Sci. J.
6
,
180
(
2013
).
11.
G.
Froyland
,
G. A.
Gottwald
, and
A.
Hammerlindl
,
SIAM J. Appl. Dyn. Syst.
13
,
1816
(
2014
).
12.
P. J.
Schmid
,
J. Fluid Mech.
656
,
5
(
2010
).
13.
M.
Budišić
,
R.
Mohr
, and
I.
Mezić
,
Chaos
22
,
047510
(
2012
).
14.
J. H.
Tu
,
C. W.
Rowley
,
D. M.
Luchtenburg
,
S. L.
Brunton
, and
J. N.
Kutz
,
J. Comput. Dyn.
1
,
391
(
2014
).
15.
M. O.
Williams
,
I. G.
Kevrekidis
, and
C. W.
Rowley
,
J. Nonlinear Sci.
25
,
1307
(
2015
).
16.
That has produced the dataset.
17.
A. J.
Chorin
,
O. H.
Hald
, and
R.
Kupferman
,
Physica D
166
,
239
(
2002
).
18.
D.
Givon
,
R.
Kupferman
, and
A.
Stuart
,
Nonlinearity
17
,
R55
(
2004
).
19.
A. J.
Chorin
,
O. H.
Hald
, and
R.
Kupferman
,
J. Sci. Comput.
28
,
245
(
2006
).
20.
O. H.
Hald
and
P.
Stinis
,
Proc. Natl. Acad. Sci. U. S. A.
104
,
6527
(
2007
).
21.
J.
Wouters
and
V.
Lucarini
,
J. Stat. Mech.: Theory Exp.
2012
,
P03003
.
22.
J.
Wouters
and
V.
Lucarini
,
J. Stat. Phys.
151
,
850
(
2013
).
23.
M. D.
Chekroun
,
H.
Liu
, and
S.
Wang
,
Approximation of Stochastic Invariant Manifolds: Stochastic Manifolds for Nonlinear SPDEs I
, Springer Briefs in Mathematics (
Springer
,
New York
,
2015
).
24.
M. D.
Chekroun
,
H.
Liu
, and
S.
Wang
,
Stochastic Parameterizing Manifolds and Non-Markovian Reduced Equations: Stochastic Manifolds for Nonlinear SPDEs II
, Springer Briefs in Mathematics (
Springer
,
New York
,
2015
).
25.
A. J.
Chorin
and
F.
Lu
,
Proc. Natl. Acad. Sci. U. S. A.
112
,
9804
(
2015
).
26.
M.
Khashei
and
M.
Bijari
,
Expert Syst. Appl.
37
,
479
(
2010
).
27.
K. L.
Hsu
,
H. V.
Gupta
, and
S.
Sorooshian
,
Water Resour. Res.
31
,
2517
, doi: (
1995
).
28.
D.
Mukhin
,
D.
Kondrashov
,
E.
Loskutov
,
A.
Gavrilov
,
A.
Feigin
, and
M.
Ghil
,
J. Clim.
28
,
1962
(
2015
).
30.
C.
Penland
and
P. D.
Sardeshmukh
,
J. Clim.
8
,
1999
(
1995
).
31.
S.
Kravtsov
,
D.
Kondrashov
, and
M.
Ghil
,
J. Clim.
18
,
4404
(
2005
).
32.
S.
Kravtsov
,
D.
Kondrashov
, and
M.
Ghil
, in
Stochastic Physics and Climate Modeling
, edited by
T. N.
Palmer
and
P.
Williams
(
Cambridge University Press
,
2009
), pp.
35
72
.
33.
A. J.
Majda
and
J.
Harlim
,
Nonlinearity
26
,
201
(
2013
).
34.
D.
Kondrashov
,
M. D.
Chekroun
, and
M.
Ghil
,
Physica D.
297
,
33
(
2015
).
35.
Polynomial functions, sigmoid functions, etc…
36.
A.
Zakharova
,
S.
Loos
,
J.
Siebert
,
A.
Gjurchinovski
,
J. C.
Claussen
, and
E.
Schöll
, in
Control of Self-Organizing Nonlinear Systems
, edited by
P. H. E.
Schöll
,
S. H. L.
Klapp
(
Springer
,
Berlin/Heidelberg
,
2016
), pp.
35
72
.
37.
K.-J.
Engel
and
R.
Nagel
,
One-Parameter Semigroups for Linear Evolution Equations
, Graduate Texts in Mathematics (
Springer-Verlag
,
New York
,
2000
), Vol.
194
, pp.
xxii
x+586
.
38.
J.
Van Neerven
,
The Asymptotic Behaviour of Semigroups of Linear Operators
(
Birkhäuser
,
2012
), Vol.
88
.
39.
K.-J.
Engel
and
R.
Nagel
,
A Short Course on Operator Semigroups
(
Springer Science & Business Media
,
2006
).
40.
H.
Brézis
,
Opérateurs Maximaux Monotones Et Semi-Groupes De Contractions Dans Les Espaces De Hilbert
(
Elsevier
,
1973
), Vol.
5
.
41.
M. D.
Chekroun
,
M.
Ghil
,
H.
Liu
, and
S.
Wang
,
Discrete Contin. Dyn. Syst., A
36
,
4133
(
2016
).
42.
A.
Pazy
,
Semigroups of Linear Operators and Applications to Partial Differential Equations
(
Spring-Verlag
,
New York
,
1983
).
43.
H.
Brézis
,
Functional Analysis, Sobolev Spaces and Partial Differential Equations
(
Springer
,
2010
).
44.
P.
Lax
,
Functional Analysis
(
Wiley
,
2002
).
45.
H.
Weinberger
,
Variational Methods for Eigenvalue Approximation
(
SIAM
,
1974
), Vol.
15
.
46.
M.
Embry
and
A.
Lambert
,
J. Math.
7
,
333
(
1977
).
47.
Another name for the eigenvalue problem L ϕ [ Ψ ] = λ Ψ .
48.
That we adopt to lie in [ 0 , 2 π ) in this article.
49.
Except at the zero frequency where the singular values provide d eigenvalues, not paired.
50.
L.
Autonne
,
Ann. Univ. Lyon, Nouvelle Sér. I
38
,
1
(
1915
).
51.
T.
Takagi
,
Jpn. J. Math.: Trans. Abstr.
(
The Mathematical Society of Japan
,
1924
), Vol.
1
, pp.
83
93
.
52.
Z.
Drmac
,
Electron. Trans. Numer. Anal.
44
,
593
(
2015
).
53.
J.
Danciger
,
Linear Algebra Appl.
412
,
22
(
2006
).
54.
J.-P.
Eckmann
and
D.
Ruelle
,
Rev. Mod. Phys.
57
,
617
(
1985
).
55.
P.
Collet
and
J. P.
Eckmann
,
Concepts and Results in Chaotic Dynamics: A Short Course: A Short Course
(
Springer
,
2007
).
56.
M. D.
Chekroun
,
A.
Tantet
,
H.
Dijkstra
, and
J. D.
Neelin
, “
Mixing Spectrum in Reduced Phase Spaces of Stochastic Differential Equations. Part I: Theory
,”
Physica D
(submitted).
57.
L.-S.
Young
,
J. Stat. Phys.
108
,
733
(
2002
).
58.
M. D.
Chekroun
,
J. D.
Neelin
,
D.
Kondrashov
,
J. C.
McWilliams
, and
M.
Ghil
,
Proc. Natl. Acad. Sci U.S.A.
111
,
1684
(
2014
).
59.
J. T.
Zhou
,
Math. Appl. (Wuhan)
9
,
53
(
1996
).
60.
A.
Bose
and
J.
Mitra
,
Stat. Probab. Lett.
60
,
111
(
2002
).
61.
H.
Karner
,
J.
Schneid
, and
C.
Ueberhuber
,
Linear Algebra Appl.
367
,
301
(
2003
).
62.
C. R.
Vogel
,
Computational Methods for Inverse Problems
(
SIAM
,
2002
), Vol.
23
.
63.
A.
Tantet
,
M. D.
Chekroun
,
H.
Dijkstra
, and
J. D.
Neelin
, preprint arXiv:1705.07573 (
2017
).
64.
A.
Selivanov
,
J.
Lehnert
,
T.
Dahms
,
P.
Hövel
,
A.
Fradkov
, and
E.
Schöll
,
Phys. Rev. E
85
,
016201
(
2012
).
65.
The latter point has been empirically observed on various dataset (not shown); see also Sec. VIII C.
66.
D.
Kondrashov
and
P.
Berloff
,
Geophys. Res. Lett.
42
,
1543
, doi: (
2015
).
67.
By a standard EOF analysis, for instance; see e.g. Ref. 66 for such a preprocessing.
68.
E.
Lorenz
, “
Predictability
,” in
Proceedings of the 1995, ECMWF Seminar
(
1996
), p.
1
.
69.
Here we consider the simplified version of the full L96 model given in [Ref. 68, Eqs. (3.2) and (3.3)] that does not include coupling between slow and fast variables.
70.
Since X ( t ) is aperiodic and thus not in the domain D(A); see Sec. VIII A.
71.
R.
Dalang
and
L.
Quer-Sardanyons
,
Expositiones Mathematicae
29
,
67
(
2011
).
72.
J.
Walsh
,
École D'Été De Probabilités De Saint Flour XIV-1984
(
Springer
,
1986
), pp.
265
439
.
73.
G.
Da Prato
and
J.
Zabczyk
, in
Stochastic Equations in Infinite Dimensions, Encycloped Ed.
, edited by
G.-C.
Rota
(
Cambridge University Press
,
Cambridge
,
1992
), Vol.
44
, pp.
xviii+454
.
74.
In which the solutions are sought as a random function of time, with values in an appropriate infinite dimensional function space of the space variable.73
75.
G.
Da Prato
and
J.
Zabczyk
,
Ergodicity for Infinite Dimensional Systems
(
Cambridge University Press
,
1996
),
Vol. 229
.
76.
M.
Gradinaru
,
I.
Nourdin
, and
S.
Tindel
,
J. Funct. Anal.
228
,
114
(
2005
).
77.
J.
Swanson
,
Ann. Probab.
35
,
2122
(
2007
).
78.
M. D.
Chekroun
,
E.
Park
, and
R.
Temam
,
J. Differ. Equations
260
,
2926
(
2016
).
79.
As discussed at the end of Sec. V.
80.
Especially when the stopping criterium of Ref. 34 (Appendix A) is not satisfied.
81.
S. N.
Stechmann
and
J. D.
Neelin
,
J. Atmos. Sci.
68
,
2955
(
2011
).
82.
S. N.
Stechmann
and
J. D.
Neelin
,
J. Atmos. Sci.
71
,
3269
(
2014
).
83.
C.
Penland
and
B. D.
Ewald
,
Philos. Trans. R. Soc. A
366
,
2455
(
2008
).
84.
M.
Ghil
,
M. R.
Allen
,
M. D.
Dettinger
,
K.
Ide
,
D.
Kondrashov
,
M. E.
Mann
,
A. W.
Robertson
,
A.
Saunders
,
Y.
Tian
,
F.
Varadi
, and
P.
Yiou
,
Rev. Geophys.
40
,
1003
, doi: (
2002
).
85.
And also naturally adapted for obtaining simulation through parallelization across the frequencies; see Sec. VIII C.
86.
D.
Kondrashov
,
M. D.
Chekroun
,
X.
Yuan
, and
M.
Ghil
, “
Data-adaptive harmonic decomposition and stochastic modeling of arctic sea ice
,” in
Advances in Nonlinear Geosciences
, edited by
A.
Tsonis
(
Springer
, in press).
87.
D.
Kondrashov
,
M. D.
Chekroun
, and
M.
Ghil
, “
Data-adaptive Harmonic Decomposition and Prediction of Arctic Sea Ice Extent
,”
Dynamics and Statistics of the Climate System
(submitted).
88.
D.
Kondrashov
and
M. D.
Chekroun
,
J. Atmos. Sol.-Terrestrial Phys
(submitted).
89.
M. D.
Chekroun
,
D.
Kondrashov
, and
M.
Ghil
,
Proc. Natl. Acad. Sci. U.S.A.
108
,
11766
(
2011
).
90.
M. D.
Chekroun
,
E.
Simonnet
, and
M.
Ghil
,
Physica D
240
,
1685
(
2011
).