Estimating Phase from Observed Trajectories Using the Temporal 1-Form

This section provides details of the method we used to generate randomized test systems with adjustible complexity and noise characteristics. We provide these details in the hope that other investigators may use them to evaluate the performance of algorithms which aim to analyze the structure of nonlinear oscillators.

Very few oscillators admit a closed-form representation of their phase as a function of state, making validation of our algorithm a significant challenge all on its own. To produce a well specified class of tests for our algorithm, we produced sample paths arising from a Stratonovich stochastic differential equation of the form

and requiring that the deterministic dynamics of $\mathrm{d}{x}(t)=f(x(t))\mathrm{d}{t}$ are known to be a (hyperbolic) limit cycle oscillator, with a known phase coordinate.

We obtained such a deterministic oscillator by starting with a randomly chosen linearized (Floquet) structure converging on the trajectory $x(t)=[t,0,\ldots]$, then winding this affine system around the unit circle with $x_{1}$ as the angle, and finally pushing forward the dynamics through a composition of randomly generated diffeomorphisms to generate $f(\cdot)$ of eqn. (12). We chose diffeomorphisms such that their tangent maps were computable and invertible in closed form. We used a second composition of random diffeomorphisms to generate a function $g(\cdot)$, and the inverse of the tangent map of $g(\cdot)$ was used for $G(\cdot)$. By using the inverse diffeomorphisms, we directly computed the phase of each trajectory point with respect to the deterministic dynamics, to use as a baseline for comparing with other forms of phase estimation.

Differential forms $\Omega$ over $\mathbf{X}$ are an exterior algebra generated using two operators: the “exterior product” $\wedge$—the universal skew symmetric product, and the “exterior derivative” ${\mathbf{d}}$. Elements of rank 0 correspond to (sufficiently smooth) scalar valued functions $\mathbf{X}\to\mathbb{C}$. When $\alpha\in\Omega_{r}$ and $\beta\in\Omega_{s}$, ${\mathbf{d}}\alpha\in\Omega_{r+1}$ and $\alpha\wedge\beta\in\Omega_{r+s}$.

These operations are familiar to many in their 3-dimensional special cases, where ranks $\Omega_{0}\ldots\Omega_{3}$ correspond to scalar functions, vectors, directed areas, and directed volumes. Here the exterior derivative ${\mathbf{d}}$ functions as the gradient, curl, and divergence. The exterior product $\wedge$ functions as a scalar-vector product when one argument is rank 0, as a cross product (resulting in a directed area) when the arguments are rank 1, and as a box-product when applied to rank 2 and rank 1 arguments.

Formally defined, $\mathcal{C}^{k}$ differential 1-forms $\alpha$ on $\mathbf{X}$ are $\mathcal{C}^{k}$ sections of its cotangent bundle $\mathsf{T}^{*}\mathbf{X}$: $\alpha(x)\in\mathsf{T}^{*}_{x}\mathbf{X}$. If $v$ is a $\mathcal{C}^{k}$ vector field, then $\langle\alpha,v\rangle$ is a $\mathcal{C}^{k}$ scalar function. If $\alpha$ is a $\mathcal{C}^{1}$ curve in $\mathbf{X}$, the line integral $\int_{\alpha}\alpha$ is well defined and independent of coordinate system.

A 1-form $\alpha$ is “closed” if its exterior derivative vanishes. It is “exact” if it is the exterior derivative of a scalar function. The identity ${\mathbf{d}}\circ{\mathbf{d}}=0$ implies that all exact forms are closed. On contractible regions of $\mathbf{X}$ the converse is true, i.e. all closed forms are exact (a.k.a. the Poincaré Lemma). The quotient space of closed to exact $k$-forms is known as the $k$-th de-Rham cohomology group. We use some properties of this cohomology to construct our algorithm.

If a closed curve $\gamma_{1}$ can be smoothly deformed into a curve $\gamma_{2}$, $\int_{\gamma_{1}}\alpha$ of a closed 1-form $\alpha$ equals $\int_{\gamma_{2}}\alpha$. In other words, the integral of a closed 1-form depends only on the homotopy class of the underlying closed curve.

A more complete introduction to the theory of differential forms can be found in Spivak (1971); Guillemin and Pollack (2010); Lee (2013).

This self-contained section serves as an alternative introduction to our paper which contains more details and is aimed at mathematically inclined readers. It includes a description of the Temporal 1-Form assuming some familiarity with differential forms, which are reviewed in sec. A.2.

We consider continuous-time dynamical systems, or flows $\Phi$, on a smooth $n$-dimensional manifold $\mathbf{X}$. That $\Phi$ is a flow means that $\Phi:\mathbf{X}\times\mathbb{R}\to\mathbf{X}$, $\Phi_{0}=\textnormal{id}_{\mathbf{X}}$, and $\Phi_{t}\circ\Phi_{s}=\Phi_{t+s}$. We assume that the vector field $f:=\frac{\partial}{\partial t}\Phi_{t}|_{t=0}$ is well-defined and $\mathcal{C}^{k\geq 1}$, so that $\Phi\in\mathcal{C}^{k}$. The trajectories $t\mapsto\Phi_{t}({x_{0}})\eqqcolon x(t)$ satisfy the ordinary differential equation (ODE) $\dot{x}=f(x)$. A periodic orbit $o$ of period $T>0$ is a trajectory that satisfies $o(t+T)=o(t)$ with $o(t)\neq o(0)$ for all $0<t<T$. We will always assume that $T>0$ is the minimal period of a nonstationary periodic orbit. We focus on dynamical systems possessing only one nonstationary periodic trajectory $o$ and denote by $\Gamma\subset\mathbf{X}$ the image of this $o$.

We further assume that this periodic trajectory is exponentially stable, a property that holds in many practical cases. Exponentially stable limit cycles are always “normally hyperbolic” (Hirsch et al., 1977; Eldering et al., 2018); for brevity we will simply refer to these as “limit cycles”. We refer to the stability basin of a limit cycle and the dynamics within it as an “oscillator”. The asymptotic behavior of any oscillator is described fully by its “asymptotic phase”: the stability basin of an oscillator is partitioned into codimension-$1$ $\mathcal{C}^{k}$ embedded submanifolds (Eldering et al., 2018, pp. 4208–4209) called “isochrons” which are asymptotically in phase with one another Guckenheimer (1975); Winfree (1980).

For systems in which the dimension is low, the noise is small or the data is plentiful, data pairs $(x_{k},{\dot{x}}_{k}),~{}k=1\ldots N$ may be used to estimate the vector field $f$. Numerous investigators proposed such approaches, primarily in the 1990s. Classical papers include Kostelich and Yorke Kostelich and Yorke (1990), who fit arbitrary continuous dynamics to trajectories, and the rich literature on time-delay embeddings, derivatives and principal components Takens (1981); Sauer et al. (1991); Schelter et al. (2006); Gibson et al. (1992). When dynamical equations are known precisely, automatic differentiation and continuation methods can be used to compute the isochrons, a problem equivalent to that of phase estimation Huguet and de-la Llave (2012); Osinga and Moehlis (2010). However, care must be taken if it is desired to extract phase response curves from this, as these quantities are coordinate dependent (Wilshin et al., 2021).

With fewer data points and higher levels of noise, one expects that the full off-cycle dynamics are less amenable to accurate estimation, and that only a few of the slower-decaying “modes” might be observable off of the limit cycle Revzen and Guckenheimer (2012). At the extreme, only phase itself might be detectable, through methods such as we proposed in Revzen & Guckenheimer Revzen and Guckenheimer (2008) using various phase-locking approaches Peters and Schaal (2008); Schaal and Atkeson (2010); Schaal et al. (2005). This motivates the main goal of the present paper: estimate asymptotic phase from empirically observed dynamics, e.g., from an ensemble of noisy measurements of (possibly short and noisy) system trajectory segments.

We note that, while there has been work to define generalized notions of asymptotic phase for stochastic oscillators (Schwabedal and Pikovsky, 2013; Thomas and Lindner, 2014; Cao et al., 2020; Engel and Kuehn, 2021) (see Pikovsky (2015); Thomas and Lindner (2015) for a spirited discussion), in this paper we restrict ourselves to estimation of the classical asymptotic phase of a deterministic oscillator using data from a perturbed version of the underlying deterministic system. We also note that there are operator-theoretic methods of recent interest revealing phase as the generator of a family of eigenfunctions of the Koopman operator with purely imaginary eigenvalues Mauroy and Mezić (2012); Kvalheim and Revzen (2021); Kvalheim et al. (2021a). These eigenfunctions and others can be estimated using Fourier/Laplace averages Mauroy and Mezić (2012); Mauroy et al. (2013); Kvalheim and Revzen (2021); Kvalheim et al. (2021a); Brunton et al. (2021) when dynamical equations are known, and from data using Dynamic Mode Decomposition Kutz et al. (2016), at least in certain cases. However, the results of these emerging spectral methods seem to be sensitive to the choice of observables and their nonlinearities (Wu et al., 2021).

Several investigators have pursued the construction of linearized (Floquet) models for dynamics near the limit cycle. We proposed an approach termed “Data Driven Floquet Analysis (DDFA)” Revzen et al. (2011); Revzen (2009) consisting of estimating phase, then constructing affine models conditioned on phase. Tytell Tytell (2013) reported a DDFA approach based on harmonic balance. Wang and Srinivasan Wang and Srinivasan (2012) proposed to construct a Floquet model using “factored Poincaré maps” and derive a phase estimate from this model. Ankarali et. al Ankarali and Cowan (2014) applied “harmonic transfer functions” and achieved good results on a “hybrid” Back et al. (1993); Simić et al. (2000); Westervelt et al. (2003); Ames and Sastry (2005); Goebel et al. (2009); Lerman (2016); Lerman and Schmidt (2020); Clark and Bloch (2021); Kvalheim et al. (2021b) spring-mass hopper model. Regardless of how they are obtained, linear models of the dynamics around the limit cycle can at best only represent the hyperplanes tangent to the isochrons at their intersection with the limit cycle. By contrast, the method we propose here can construct nonlinear isochron approximations on suitable neighborhoods of the limit cycle, as we explain in sec. A.5.

We now begin our theoretical description of the object we call the “Temporal 1-Form”. By redefining $\mathbf{X}$ to be the basin of attraction of $\Gamma$, we henceforth assume that $\Gamma$ is globally asymptotically stable. Asymptotic phase can be viewed as a map $\mathrm{P}:\mathbf{X}\to\Gamma$ which is a retraction ($\mathrm{P}|_{\Gamma}=\textnormal{id}|_{\Gamma}$) and a semiconjugacy ($\forall t\in\mathbb{R}:\mathrm{P}\circ\Phi_{t}=\Phi_{t}\circ\mathrm{P}$). From this it follows that $\mathrm{P}\in\mathcal{C}^{k}$ when $f\in\mathcal{C}^{k}$⁶⁶6 In more detail: fix $x\in\mathbf{X}$, let $U\subsetneq\Gamma$ be a contractible open neighborhood of $\mathrm{P}(x)$, and define $\tau_{x}\in[0,\infty)$ to be the smallest nonnegative number such that $\mathrm{P}\circ\Phi_{-\tau_{x}}({x})=o(0)$. Since the isochron $\mathrm{P}^{-1}(o(0))$ is a $\mathcal{C}^{k\geq 1}$ manifold transverse to the vector field Guckenheimer (1975); Eldering et al. (2018), the implicit function theorem implies that there is a $\mathcal{C}^{k}$ function $\tau\colon\mathrm{P}^{-1}(U)\to\mathbb{R}$ satisfying $\tau(x)=\tau_{x}$ and $\mathrm{P}\circ\Phi_{-\tau(y)}({y})=o(0)$ for all $y\in\mathrm{P}^{-1}(U)$. Thus, $\mathrm{P}|_{\mathrm{P}^{-1}(U)}=o\circ\tau$ is $\mathcal{C}^{k}$. Since $x\in\mathbf{X}$ was arbitrary, it follows that $\mathrm{P}\in\mathcal{C}^{k}$. . Each $p\in\Gamma$ represents the “isochron” $\mathrm{P}^{-1}(p)$, which is the set of initial conditions from which trajectories converge with the one starting at $p$⁷⁷7 The geometry of isochrons can be surprisingly complicated, especially when the vector field $f$ has multiple time scales (Langfield et al., 2015) . In the present context of oscillators, asymptotic phase is more commonly represented as a circle-valued (“phasor”-valued) map $\varphi\colon\mathbf{X}\to\mathrm{S}^{1}$; we explain the relationships between different representations of phase in sec. A.4.

The Temporal 1-Form appears naturally as a consequence of the existence of $\mathrm{P}$. Noting that $f$ is nowhere zero on the $1$-dimensional $\mathcal{C}^{k}$ manifold $\Gamma$ and defining the angular frequency $\omega\coloneqq 2\pi/T$, it follows that there exists a unique $\mathcal{C}^{k-1}$ differential 1-form $\alpha$ on $\Gamma$ satisfying $\langle\alpha,f\rangle=\omega$ identically on $\Gamma$⁸⁸8In a Euclidean space, $\alpha$ can be expressed as $\alpha:=(f^{\mathsf{T}}f)^{-1}\omega f^{\mathsf{T}}$. The “Temporal 1-Form ” $d\phi:\,\mathbf{X}\to\mathsf{T}^{*}\mathbf{X}$ is defined to be the pullback $\mathrm{P}^{*}\alpha$ of $\alpha$, which means that $\langle d\phi(x),v\rangle:=\langle\alpha(\mathrm{P}(x)),\mathsf{D}\mathrm{P}(x)\cdot v\rangle$ for any $v\in\mathsf{T}_{x}\mathbf{X}$. Since $\mathrm{P}\in\mathcal{C}^{k}$, it follows that $d\phi\in\mathcal{C}^{k-1}$. It is easy to see that the Temporal 1-Form satisfies two properties (further discussed in sec. A.4): (1) $\langle d\phi,f\rangle=\omega$ everywhere, and (2) $d\phi$ is closed.

By “closed” we mean that each $x\in\mathbf{X}$ possesses an open neighborhood $U_{x}$ on which $d\phi$ is the exterior derivative of a $\mathcal{C}^{k}$ function. This definition permits us to discuss continuous closed forms, needed for the case $f\in\mathcal{C}^{1}$; it reduces to the usual definition for $\mathcal{C}^{1}$ forms (cf. Farber et al. (2004)).

Somewhat less obvious is the fact that the Temporal 1-Form is the unique continuous closed $1$-form on $\mathbf{X}$ satisfying the preceding two properties; we show that this is the case in sec. A.4. In that section we also relate $\mathrm{P}$, $d\phi$, and circle-valued phases $\varphi\colon\mathbf{X}\to\mathrm{S}^{1}$. By a “circle-valued (asymptotic) phase” $\varphi\colon\mathbf{X}\to\mathrm{S}^{1}\subset\mathbb{C}$ we mean a continuous map satisfying $\varphi\circ\Phi_{t}=e^{i\omega t}\varphi$ for all $t\in\mathbb{R}$. To summarize Theorem 1 of sec. A.4, some of these relationships are as follows: $\varphi\circ\mathrm{P}=\varphi$, $\varphi$ is unique modulo rotations of $\mathrm{S}^{1}$, and $d\phi=\varphi^{*}(d\theta)$ is the pullback of the standard angular form on $\mathrm{S}^{1}$ via any circle-valued phase. Moreover, any choice of basepoint $x_{0}\in\mathbf{X}$ uniquely determines a circle-valued phase by integrating $d\phi$ along continuous curves from $x_{0}$.

While we are interested in asymptotic phase as defined for deterministic oscillators, our primary goal is to compute asymptotic phase from real-world data which we assume may be subject to system and measurement noise. In sec. A.5 we derive fairly general performance guarantees applicable to our algorithm and potentially others. In that section, we define “unobserved” and “observed” empirical cost functions $J_{N}$ and $\widehat{J_{N}}$ mapping $1$-forms to nonnegative numbers, where $N$ is the number of observed state-velocity data pairs. Intuitively, we would like to minimize the former cost (for which $J_{N}(d\phi)=0$), but only the latter cost is observable, so our algorithms attempts to minimize the latter.

Assuming that our algorithm outputs estimates $\widehat{d\phi}_{N}$ satisfying $\widehat{d\phi}_{N}\leq d\phi+\epsilon$ for some $\epsilon\geq 0$, that $\mathbf{X}$ is an open subset of $\mathbb{R}^{n}$, and that state measurements belong to some fixed compact set $K\subset\mathbf{X}$, we establish in Theorem 2 the following inequality holding with probability $1$:

Here $\|\,\cdot\,\|_{\mathcal{C}^{0}(K)}$ denotes the supremum norm. Assuming $\epsilon=0$ for simplicity, it follows in particular that, with $f$ fixed, the performance of the estimator with respect to $J_{N}$ is within $\mathcal{O}(\kappa)$ of optimality in the limit of large data. Moreover, the full statement of Theorem 2 also includes explicit convergence rates (see Remark 3). See the remarks following Theorem 2 for further discussion.

However, as discussed in sec. A.5, bounding the performance with respect to $J_{N}$ as in Theorem 2 does not say anything about closeness of the estimates $\widehat{d\phi}_{N}$ themselves to $d\phi$. Thus, even if $J_{N}(\widehat{d\phi}_{N})$ is small, the estimated isochrons corresponding to $\widehat{d\phi}_{N}$ need not resemble the true isochrons. This motivates Theorem 3 which shows, roughly speaking, that under some additional assumptions the estimated isochrons converge $\mathcal{C}^{1}$-uniformly to the true isochrons on forward invariant open subsets of $K$. Moreover, the full statement of Theorem 1 also includes explicit convergence rates.

In this section we explain the relationship (mentioned in sec. A.3) between the asymptotic phase map $\mathrm{P}\colon\mathbf{X}\to\Gamma$, circle-valued phase maps $\varphi\colon\mathbf{X}\to\mathrm{S}^{1}$, and the Temporal 1-Form $d\phi$. We continue under the assumptions and notation of sec. A.3; in particular, the state space $\mathbf{X}$ is also the basin of $\Gamma$.

Given an open neighborhood $U\subset\mathbf{X}$ of $\Gamma$ which is forward invariant for the flow of $f$, we extend the definition of circle-valued phase (sec. A.3) to include continuous maps $\varphi\colon U\to\mathrm{S}^{1}\subset\mathbb{C}$ satisfying $\varphi\circ\Phi_{t}=e^{i\omega t}\varphi$ for all $t\geq 0$. (We use this generality in the proof of Theorem 3.) For the following lemmas and theorem, $d\theta$ denotes the standard angular form on the circle $\mathrm{S}^{1}$ and $i=\sqrt{-1}$.

The following result summarizes some relationships between different representations of asymptotic phase.

In this section we prove two theorems relating the performance of the true Temporal 1-Form to that of estimates computed in the presence of noise, to which data from real-world systems and measurements are subjected. These are extensions of corresponding results in Kvalheim (2018, Sec. 3.4).

Throughout the remainder of this section, we assume for simplicity that the state space $\mathbf{X}$ (also assumed to be the basin of attraction of the limit cycle) is an open subset of $\mathbb{R}^{n}$. We also assume for the remainder of this section that data takes values in a fixed compact set $K\subset\mathbf{X}$.

In what follows we let $\|\,\cdot\,\|$ be the Euclidean norm, and we denote induced operator norms by the same symbol. Given a continuous function $\beta$ on a set $A\subset\mathbf{X}$, we define

We denote by $\Omega^{1}_{0}(A)$ the set of continuous $1$-forms over $A$, i.e., continuous sections of $\mathsf{T}^{*}\mathbf{X}|_{A}$.