A Bayesian nonparametric approach to the approximation of the global stable manifold

TL;DR

This paper introduces a Bayesian nonparametric method using MCMC to reconstruct the global stable manifold from noisy time-series data, applicable to noninvertible maps, demonstrated on polynomial maps with synthetic data.

Contribution

It presents the SM-GSBR model, a novel Bayesian approach for estimating the global stable manifold and underlying dynamics from observed data, including noninvertible maps.

Findings

Successfully reconstructs stable manifolds from synthetic data.

Works for noninvertible maps without modifications.

Provides a new tool for analyzing noisy dynamical systems.

Abstract

We propose a Bayesian nonparametric model based on Markov Chain Monte Carlo (MCMC) methods for unveiling the structure of the invariant global stable manifold from observed time-series data. The underlying unknown dynamical process is possibly contaminated by additive noise. We introduce the Stable Manifold Geometric Stick Breaking Reconstruction (SM-GSBR) model with which we reconstruct the unknown dynamic equations and in parallel we estimate the global structure of the perturbed stable manifold. Our method works for noninvertible maps without modifications. The stable manifold estimation procedure is demonstrated specifically in the case of polynomial maps. Simulations based on synthetic time series are presented.

Figures

Equations

$$W^s(y)=\left\{x\in {\mathbb{R}}^2:g^n(x)\to y,n\to \infty \right\}.$$

$$W^u(y)=\left\{x\in {\mathbb{R}}^2:g^{-n}(x)\to y,n\to \infty \right\}.$$

$$(x,y)\mapsto g(\vartheta ,x,y)=({\theta }_1+{\theta }_2{x}^2+{\theta }_{3}y,x),$$

$$(x,y)\mapsto g(\vartheta ,x,y)=({\theta }_1+{\theta }_{2}x+{\theta }_3{x}^3+{\theta }_{4}y,x).$$

$$X_i=\mathcal{T}(\vartheta ,X_{i-d:i-1},e_i)=g(\vartheta ,X_{i-d:i-1})+e_i,i\ge 1,$$

$$e_i\stackrel{i.i.d.}{\sim }f(x|\mathbb{G})={\int }_{\tau >0}\mathcal{N}(x|0,{\tau }^{-1})\mathbb{G}(d\tau ),$$

$${\mathbb{G}}_{\mathrm{DP}}=\sum _{k=1}^{\infty }{w}_k{\delta }_{{\tau }_k}\sim \mathcal{DP}(c,G_0).$$

$$\mathbb{G}=\sum _{k=1}^{\infty }\mathbb{E}(w_k){\delta }_{{\tau }_k}\sim \mathcal{G}\mathcal{SB}(\alpha ,\beta ,G_0).$$

$$e_i\stackrel{i.i.d.}{\sim }f(x|\lambda ,{\tau }_{1:\infty })=\lambda \sum _{k=1}^{\infty }(1-\lambda {)}^{k-1}\mathcal{N}(x|0,{\tau }_j^{-1}).$$

$$\Pi (\lambda ,{\tau }_{1:\infty },\vartheta ,x_{-T-d+1:-T},x_{-T+1:0}|x_{1:n}).$$

$$\Pi (d_i,N_i|\lambda )\propto (1-\lambda {)}^{N_i-1}\mathcal{I}(d_i\le N_i),$$

$$\Pi (x_i,d_i,N_i|\lambda ,{\tau }_{1:\infty })={\lambda }^2(1-\lambda {)}^{N_i-1}\mathcal{N}(x_i|0,{\tau }_{d_i}^{-1}).$$

$$\Pi (x_i,d_i,N_i|\lambda ,{\tau }_{1:\infty },\vartheta ,x_{i-d:i-1})={\lambda }^2(1-\lambda {)}^{N_i-1}$$

$$\times \mathcal{N}(x_i|g(\vartheta ,x_{i-d:i-1}),{\tau }_{d_i}^{-1})\mathrm{when}{d}_i\le N_i.$$

$$\Pi (x_{-T-d+1:-T})\propto \prod _{1\le i\le d}\mathcal{I}(x_{-T-d+i}\in \mathcal{X}).$$

$$\Pi (x_{-T+1:0}|{\mathcal{R}}_T)\propto \Pi (x_{-T+1:0})\prod _{1\le i\le T}\mathcal{I}(x_{-T+i}\in \mathcal{X}).$$

$$\Pi (\lambda )\Pi ({\tau }_{1:\infty })\Pi (\vartheta )\prod _{-T+1\le j\le d}\mathcal{I}(x_{j-d}\in \mathcal{X})\prod _{\begin{array}{c}-T+1\le j\le n\\ d_j\le N_j\end{array}}{\lambda }^2(1-\lambda {)}^{N_j-1}\mathcal{N}(x_j|g(\vartheta ,x_{j-d:j-1}),{\tau }_{d_j}^{-1}).$$

$$g(\vartheta ,x_{i-2},x_{i-1})={\vartheta }_0+{\vartheta }_1{x}_{i-1}+{\vartheta }_2{x}_{i-2}+{\vartheta }_3{x}_{i-1}{x}_{i-2}+{\vartheta }_4{x}_{i-1}^2+{\vartheta }_5{x}_{i-2}^2.$$

$${\tau }_{1:N^{\ast }},d_{-T+1:n},N_{-T+1:n}$$

$$X_{-T-d-1:-T},X_{-T+1:0},\vartheta ,\lambda ,e_{n+1}.$$

$$f(x_{-T-d+j}|\cdots )\propto \mathcal{I}\left(x_{-T-d+j}\in \mathcal{X}\right)\exp \left\{-\frac{1}{2}\sum _{1\le l\le j}{\tau }_{d_{-T+l}}{\left[x_{-T+l}-g(\vartheta ,x_{-T+l-d:-T+l-1})\right]}^2\right\},1\le j\le d,$$

$$f(x_{-T+j}|\cdots )\propto \mathcal{I}\left(x_{-T+j}\in \mathcal{X}\right)\exp \left\{-\frac{1}{2}\sum _{0\le l\le d}{\tau }_{d_{-T+j+l}}{\left[x_{-T+j+l}-g(\vartheta ,x_{-T+j+l-d:-T+j+l-1})\right]}^2\right\},1\le j\le T,$$

$$\mathrm{starting}\text{leavevmode}seq.:\text{leavevmode}f(x_{-T-d+1:-T}|x_{-T+1:-T+d},\cdots )$$

$$\mathrm{intermediate}\text{leavevmode}seq.:\text{leavevmode}f(x_{-T+1:0}|x_{-T-d-1:-T},x_{1:d},\cdots ).$$

$$\lambda \sim \mathcal{B}e(\cdot |0.5,0.5),{\tau }_{1:\infty }\sim \prod _{i=1}^{\infty }\mathcal{G}(\cdot |10^{-3},10^{-3}).$$

$$y_j=g(y_{j-2},y_{j-1})=1.38-y_{j-1}^2+0.27y_{j-2},$$

$$x_i^{(j)}=g(x_{i-2}^{(j)},x_{i-1}^{(j)})+e_i^{(j)},$$

$$e_i^{(j)}\stackrel{i.i.d.}{\sim }f=\frac{9}{10}\mathcal{N}(0,10^{-7})+\frac{1}{10}\mathcal{N}(0,10^{-3}).$$

$${\hat{W}}^s(T;r,n)=\bigcup _{1\le j\le r}\mathrm{supp}\left(X_{-T-d+1:-T}|x_{1:n}^{(j)}\right).$$

$${\hat{W}}_n^s(T;k)=\bigcup _{1\le j\le k}\mathrm{supp}\left(X_{-T-d+1:-T}|x_{j:n-k+j}\right).$$

$$x_i=1-1.4x_{i-1}^2+0.3x_{i-2},$$

$${\hat{W}}_{2000}^s=\bigcup _{1\le j\le 500}\mathrm{supp}(X_{-4},X_{-3}|x_{j:j+1500}).$$

$$e_i\stackrel{i.i.d.}{\sim }f=\frac{3}{4}\mathcal{N}(0,{\sigma }_1^2)+\frac{1}{4}\mathcal{N}(0,{\sigma }_2^2),$$

Full text

Keywords: Bayesian nonparametric inference; Chaotic dynamical systems; Stable manifold; Random dynamical systems

1 Introduction

In the recent literature, there has been a growing interest for the construction and analysis of models that include nonlinear deterministic components and are influenced by a variety of noise processes. There exist two main distinct fields of research associated with random perturbations in dynamical systems. The first approach regards the analysis of Random Dynamical Systems, by defining direct stochastic generalizations of the central notions appearing in the theory of Dynamical Systems, such as random attractors [4] and random invariant manifolds [27]. The second approach, which we will follow in this paper, relates to the behavior of the perturbed system to its zero-noise limit, that is having as a reference system the underlying nonlinear deterministic part.

Invariant manifolds are essential to the theory of dynamical systems, as the behavior of any dynamical system is related to the underlying geometrical structure of the state space, more specifically the organization of the invariant stable and unstable subspaces. There exist several methods in the literature suitable for the approximation of invariant manifolds of deterministic and random dynamical systems.

The need for a development of generic methods suitably approximating the invariant manifolds of nonlinear dynamical systems is highlighted by the fact that in general it is not possible to describe the invariant set of points via closed-form analytic expressions. A special case of a power series expansion approximation for the stable and unstable manifolds of a Hénon map saddle fixed point is presented in Ref. [9]. The most widely used method for the approximation of the global stable manifold of any invertible map is presented in Ref. [34], while the authors also propose a method suitable for noninvertible maps in [20]. Another method for the approximation of the global stable manifold of invertible maps is introduced in [22] and is extended to the Search-Circle algorithm [7] for the noninvertible case. A modified version of the Search-Circle, which is faster and is based on the same concept is presented in [25]. For the numerical computation of higher-dimensional invariant manifolds, we refer to [13, 21, 14]. For a review of numerical methods for two-dimensional manifolds we refer to [23]. As a restriction to the above methods, can be considered the requirement of the complete knowledge of the dynamical system. A method suitable for the approximation of the invariant manifolds given only experimental data, and when there is no available model, is proposed in Ref. [32]. Moreover, the “parameterization method” originally introduced by Cabré et. al [1, 2], forms the basis of a series of results relevant to the establishment of the existence and the computation [33] of invariant manifolds.

In this work, we will use a Bayesian nonparametric approach to reconstruct the unknown random dynamical system and simultaneously jointly estimate variables of past unobserved-observations based solely on time series data. More specifically, we propose an extension of the GSBR sampler introduced in Ref. [26], in order to provide a MCMC based stochastic approximation of the global stable manifold based on observed noisy time-series. We introduce the Stable Manifold Geometric Stick Breaking model, from now on referred to as the SM-GSBR model. Our intention is to compute the support of the marginal posterior densities of certain sets of past unobserved-observations by performing predictions in reversed time. We emphasize on the diffusive support of marginal posterior densities of far-off past random variables, modeling initial conditions. We demonstrate that such supports, contain part of the associated noisy stable manifold. Finally, when the SM-GSBR sampler is applied on proper sliding window time-series, the global noisy stable manifold reveals itself naturally as the union of the supports of the posterior marginal distributions of random variables describing remote in the past initial conditions.

The unique feature of the random variable which is modeling a starting point for an observed time-series, is that it is the only variable with a full conditional density, depending exclusively on observations that are ahead in time. All other variables, have full conditional densities that depend on both future and past observations. This apparent lack of information regarding past observations, makes the support of the associated full conditional to diffuse along the direction of the local stable manifold contained on a neighborhood of the true initial condition. We will see that when we perform predictions in reversed time, the fluctuations of the far-off past variables describing initial conditions is increased, thus spreading the sampled values along the stable direction.

Likewise, conditionally on proper subsets of the observed noisy time-series, we can approximate parts of the global unstable noisy manifold, by means of the diffusive supports of the random variables describing the ending points of forward out-of-sample future unobserved-observations.

We will demonstrate the efficiency of the proposed method using different types of polynomial maps, which are of particular interest not only because of their rich dynamical behavior, but also due to their ability to approximate more complicated maps, as finite degree Taylor approximations of non-polynomial nonlinearities. Furthermore, in essence all invertible polynomial maps, i.e. any nontrivial polynomial diffeomorphism of $\mathbb{R}^{2}$ with constant Jacobian, are conjugate to compositions of generalized Hénon maps [5].

The basic advantage of the proposed method is its ability to provide an adequate stochastic approximation of the stable manifold, under a data driven method. In particular, no knowledge regarding the parameters of the system is required, not even the location of the saddle fixed point, just a general functional representation of the deterministic part. Namely, the procedures of the system identification and of the stochastic approximation are performed in parallel, in a similar fashion as the Dynamic Noise Reduction Replicator model described in Ref. [19]. Our method is parsimonious, due to its flexibility induced by the general functional form of the deterministic part, and the application of a GSB mixture process prior over the additive stochastic component. Our method is applicable under small data sets which are possibly corrupted by (perhaps) non-Gaussian noise. Moreover, an important feature of the SM-GSBR model is its wide applicability. In fact we will demonstrate the applicability of the SM-GSBR algorithm equally to the case of invertible and the case of noninvertible maps.

This work relates to utilizing a Bayesian nonparametric framework, in order to approximate dynamical invariants, based entirely on observed time series data. As we have seen in a series of previous works, one can in principle:

Approximate the quasi-invariant measure, of the underlying random dynamic system, thus implying the existence of a prediction barrier [26]. 2. 2.

Calculate the positions of primary homoclinic tangencies, revealing the non-hyperbolic nature of the underlying perturbed deterministic dynamic system [19]. 3. 3.

Utilize borrowing of strength relationships among the dynamical error pairs of multiple time-series, such that underrepresented data sets benefit in terms of model estimation accuracy [16].

The layout of the paper is as follows. In Section 2, we discuss some preliminary notions. In Section 3, we derive the SM-GSBR model, a Bayesian nonparametric model incorporating a GSB mixture process prior, for making predictions in reversed time, suitable for the stochastic approximation of the global stable manifold from observed time-series data. In Section 4, we resort to simulation. We apply the proposed SM-GSBR model on deterministic and dynamically perturbed versions of the classical Hénon map, the Dual-Hénon map and a noninvertible map with a Hénon-like strange attractor. Finally, in Section 5 are the concluding remarks.

2 Preliminaries

In what follows, we will consider a planar diffeomorphism $g$ of $\mathbb{R}^{2}$ to itself. The global dynamical behavior of the map $g$ is determined by dynamically invariant objects, such as the stable and the unstable manifolds. Let $y$ be a saddle fixed point of the map $g$, i.e. the Jacobian matrix $D_{y}\,g$ has eigenvalues $\lambda_{\rm s}$ and $\lambda_{\rm u}$ such that $|\lambda_{\rm s}|<1$ and $|\lambda_{\rm u}|>1$. The global stable manifold [12, 22] $W^{s}(y)$ of $y$, is the set of points whose orbit tend to $y$ in forward time, namely

[TABLE]

Similarly, the corresponding unstable manifold is defined as the set of points whose orbits tend to $y$ in reversed time, that is

[TABLE]

Since $g$ is invertible, the stable manifold of $g$ is the unstable manifold of $g^{-1}$ and vice versa.

Due to the stretching and folding mechanisms, the two manifolds do not only intersect at the fixed point but also in other locations called homoclinic points. If the intersection $p$ between the invariant manifolds is tangential, we have a point of homoclinic tangency (HT). There are infinitely many HTs, due to the fact that $g(p)$ and $g^{-1}(p)$ are themselves HTs. HTs are important in the study of how noise affects chaotic systems, as they lead to noise amplification. In the symbolic dynamics literature, they have been proposed as a basis for building generating partitions. Furthermore, the existence of HTs is mutually exclusive with the property of hyperbolicity, thus making shadowing impossible. Dynamical perturbations in the vicinity of HTs are driving away the dynamics from the neighborhood of the attractor. Nevertheless, due to the folding effect of the nonlinear map, the noisy orbits are mapped onto the neighborhood of the attractor again. For a detailed description of the above mechanism see Ref. [18].

The Hénon map[17], is defined as

[TABLE]

where $\vartheta$ denotes the dependence of $g$ to control parameters. The Hénon map exhibits the same qualitative dynamical behavior with a wide variety of polynomial maps, namely quadratic maps with a constant Jacobian.

In Fig. 1 we present portions of the stable manifold of the saddle fixed point $(0.631,0.631)$, for the classical parameter values $\vartheta=(1,-1.4,0.3)$. We obtain points of the stable manifold by iterating a set of points in small line segments, aligned along the stable direction of the saddle fixed point [12]. We notice the emergence of homoclinic tangencies, as the map exhibits non-hyperbolic chaotic behavior. The attractor has been conjectured to be the closure of the unstable manifold.

An illustrating example, regarding the relevance of the stable manifold and the basins of attraction in the case of a multistable behavior, is the Dual-Hénon map [31]. The latter map is a special case of the Generalized-Hénon map[5], which is defined as

[TABLE]

When $\vartheta=(0,2,-0.1,0.3)$, we have the particular case of the Dual-Hénon map. In Fig. 2 we present two symmetric, with respect to the origin, isolated strange attractors along with their corresponding basins of attraction. The two basins of attraction are extending to infinity along both directions of the $y$-axis and have a fractal boundary [31].

In Fig. 3 we plot the pair of strange attractors, along with the two global stable manifolds superimposed. The two associated saddle fixed points are situated in $(3.605,3.605)$ and $(-3.605,-3.605)$. The interlaced global stable manifolds form the boundaries between the two different basins. Points situated on the stable manifolds, will eventually converge to the associated saddle fixed points.

3 The SM-GSBR model

Suppose that we have at our disposal, a time-series $x_{1:n}:=(x_{1},\ldots,x_{n})$ of length $n$, which is a realization of the random recurrence relation:

[TABLE]

where $g:\Theta\times\mathbb{X}^{d}\to\mathbb{X}$, for some compact subset $\mathbb{X}$ of $\mathbb{R}$. The states of the system $(X_{i})_{i\geq-d+1}$ and the noise perturbations $(e_{i})_{i\geq 1}$, are all real random variables over some probability space $(\Omega,{\cal F},{\rm P})$. We denote by $\vartheta\in\Theta\subseteq\mathbb{R}^{m}$ any dependence of the deterministic part of the random map on control parameters. We assume, that the deterministic part $g$ of $\cal T$, is nonlinear and continuous in the variable $X_{i-d\,:\,i-1}$.

In addition, we assume that the random variables $e_{i}$, are independent to each other and independent of the states $X_{i-r}$ for all $r<i+d$. We relax the conventional assumption of normality for the additive dynamical errors $e_{i}$, by making the more general assumption, that the $e_{i}$s are are independent and identically distributed (i.i.d.), from some symmetric with respect to zero distribution, with unknown density $f$ with support over the real line, so that ${\cal T}:\Theta\times\mathbb{X}^{d}\times\mathbb{R}\to\mathbb{R}$.

We will not take into account observational noise, so that the time series $x_{1:n}$ is completely specified by the initial point $x_{-d+1:0}$, the control parameters $\vartheta$ and the particular realization $e_{1:n}$ of the noise process. Also, we note that the observed time-series, can be considered as the one-dimensional measurements of the states of a dynamical system of unknown dimension, which can be embedded in a space of proper dimension, using delay coordinates with a suitably chosen time delay.

We propose a Bayesian nonparametric model, the SM-GSBR model, in order to jointly estimate a fixed number of consecutive past observations, by performing prediction in reversed time. Given the time-series $x_{1:n}$ and a backward prediction horizon $T$, we aim to estimate the support of the $T$-shifted to the past initial variable $X_{-T-d+1:-T}$, associated with the sequence of variables $X_{-T+1:n}$.

3.1 The infinite random mixture prior

We will model a-priori the errors in recurrence (3) as a random infinite mixture of zero-mean normal kernels with mixing measure the GSB random measure $\mathbb{G}$, as

[TABLE]

where $\tau$ is the precision (inverse variance).

The random measure $\mathbb{G}$ is closely linked to the Dirichlet random measure [8, 30]

[TABLE]

The Dirac measures $\delta_{\tau_{k}}$ are concentrated on the random precisions $\tau_{k}$, which are drawn independently from the mean parametric distribution $G_{0}$, which is the prior guess of $\mathbb{G}_{\rm DP}$ i.e. $\mathbb{E}(\mathbb{G}_{\rm DP}(A))=G_{0}(A)$, for measurable subsets $A$ of $\mathbb{R}_{+}$. The weights $w_{k}$ satisfy $\sum_{k\geq 1}w_{k}=1$ a.s. and are stick-breaking, that is $w_{1}=v_{1}$ and $w_{k}=v_{k}\prod_{l<k}(1-v_{l})$, and random because $v_{k}\stackrel{{\scriptstyle{\rm i.i.d.}}}{{\sim}}{\cal B}e(1,c)$, a beta density with shapes 1 and $c$. The GSB random measure $\mathbb{G}$ is defined by removing a hierarchy from the random measure $\mathbb{G}_{\rm DP}$[11] as

[TABLE]

Then, letting $\lambda:=\mathbb{E}(v_{k})=(1+c)^{-1}$, it is that $\mathbb{E}(w_{k})=\lambda\,(1-\lambda)^{k-1}$ (the geometric probability of the GSB prior measure), and $\mathbb{E}(\mathbb{G}(A))=G_{0}(A)$. Finally, we randomize the probability-weights of $\mathbb{G}$ by letting $\lambda\sim{\cal B}e(\alpha,\beta)$. In Ref.[26] it is shown that the $\mathbb{G}$-based Bayesian nonparametric framework for dynamical system estimation is efficient, faster and less complicated when compared to Bayesian nonparametric modeling via the Dirichlet process.

Taking into account the previous discussion, letting $\tau_{1:\infty}=(\tau_{j})_{j\geq 1}$ and the fact that the map $\mathbb{G}\mapsto(\lambda,\tau_{1:\infty})$ is one-to-one, the random mixture in (4) can be represented as

[TABLE]

Now, the posterior density $\pi(x_{-T-d+1:-T}\,|\,x_{1:n})$ of the $T$-shifted to the past initial condition variable $X_{-T-d+1:-T}$, given the observed time-series $x_{1:n}$, is a marginal density of the posterior

[TABLE]

Nevertheless, the prior mixture density over the additive errors in (5) is infinite dimensional, to create a Gibbs sampler with an a.s. finite number of components, having as its stationary distribution the latter posterior, we need to augment (6) with the following $T+n$ pairs of auxiliary variables ${\cal V}_{T}=\{(d_{i},N_{i}):i=-T+1,\ldots,n\}$. $T$ pairs corresponding to the past unobserved-observations in $x_{-T+1:0}$, and $n$ pairs corresponding to the true observations in $x_{1:n}$. Suppose that the pair of variables $(d_{i},N_{i})$ is associated with the observation $x_{i}$ (observed or unknown). Then, the variable $d_{i}$ is the infinite mixture allocation variable for $x_{i}$, with $\Pi(d_{i}\,|\,\lambda)=\lambda(1-\lambda)^{d_{i}-1}$. Moreover, if jointly $d_{i}$ and $N_{i}$ are distributed as

[TABLE]

where ${\cal I}(d_{i}\leq N_{i})$ equals to $1$ whenever $d_{i}\leq N_{i}$ and 0 elsewhere, then and only then, $d_{i}$ conditionally on $N_{i}$, will attain the discrete uniform distribution over the random set $\{1,\ldots,N_{i}\}$, thereby, making the associated (${\cal V}_{T}$-augmented) Gibbs sampler having order $N^{*}=\max_{-T+1\leq i\leq n}N_{i}<\infty$, a.s. Finally, when $d_{i}\leq N_{i}$, because $\Pi(x_{i}\,|\,d_{i},\tau_{1:\infty})={\cal N}(x_{i}\,|\,0,\tau_{d_{i}}^{-1})$, we have

[TABLE]

In fact, conditioning further the latter density with the variables $\vartheta$ and $x_{i-d:i-1}$, yields

[TABLE]

Because the global stable manifold possibly extends to infinity, to make the Gibbs sampler more efficient, we truncate the coordinates of the past unobserved-observation variables $x_{-T-d+1:-T}$ and $x_{-T+1:0}$ to a finite set $\cal X$ centered at the origin. In that fashion:

To the $T$-shifted to the past initial condition variable $x_{-T-d+1:-T}$ we will assign the uniform prior over ${\cal X}^{d}$, i.e. $x_{-T-d+1:-T}\sim{\cal U}({\cal X}^{d})$, equivalently

[TABLE] 2. 2.

To the variables intermediate between $X_{-T-d+1:-T}$ and and the observed time-series $x_{1:n}$, we apply the restriction event ${\cal R}_{T}=\{x_{-T+1:0}\in{\cal X}^{T}\}$. Then, it is clear that

[TABLE]

We have the following proposition:

Proposition $1$. The ${\cal R}_{T}$-restricted and ${\cal V}_{T}$-augmented posterior incorporating the $T$-shifted to the past initial condition variable $x_{-T-d+1:-T}$, under the assumption that a-priori and independently, $x_{-T-d+1:-T}\sim{\cal U}({\cal X}^{d})$, is proportional to:

[TABLE]

*The proof is given in the Appendix A. *

Completing the model:  Finally, we have to assign independent priors over the variables $\tau_{1:\infty}$ and $\vartheta$. Following standard Bayesian modeling, to the precisions $\tau_{i}$ we apply the conjugate gamma priors $\tau_{i}\stackrel{{\scriptstyle{\rm i.i.d.}}}{{\sim}}{\cal G}(b_{1},b_{2})$, $i\geq 1$, with fixed shape and rate hyperparameters, $b_{1}$ and $b_{2}$, respectively. Although there is an infinity of them, only at most $N^{*}<\infty$ of them will be used in the posterior computations. We model the deterministic part $g$ of $\cal T$, with the complete quadratic polynomial in the two variables, given by

[TABLE]

Over $\vartheta$, we place the prior $\Pi(\vartheta)\propto 1$. Such a prior, although improper, is acceptable due to the fact that it leads to a proper posterior.

At each iteration of the Gibbs sampler, we we will sample from the full conditional densities of the variables:

[TABLE]

It can be verified via (11), that the full conditional densities of the individual variables in the $T$-shifted to the past initial condition $x_{-T-d+i:-T}$, and the intermediate part $x_{-T+1:0}$, are given by:

[TABLE]

respectively, where $f(x\,|\cdots)$ denotes the density of $x$ conditional on the rest of the variables.

We remark that, although both previous collections of densities are modeling past out-of-sample observations, the full conditionals associated with the $T$-shifted to the past starting point in (13), exhibit the unique characteristic that they jointly depend purely on future observations, in the sense that:

[TABLE]

All other variables, have full conditionals that depend equally on both future and past observations. This apparent lack of information, makes the support of the associated full conditional diffusive along the direction of the local stable manifold contained on a neighborhood of the true (but unknown) initial point $x_{-d+1:0}$ of the observed time-series $x_{1:n}$.

3.2 The stable manifold stochastic approximation

The proposed method is based on the observation that the posterior marginal distribution of the initial condition vector lies along the stable direction, i.e. the direction of the stable manifold [15]. We use the stable direction, for example in order to find regions of HTs, by iterating a random normalized vector under the inverse Jacobian matrix at a given point[3, 24]. In practice, due to the influence of numerical errors, we use the evaluation of the inverse Jacobian at each orbit point. More details regarding the calculation of stable and unstable directions, can be found in [24].

In all our numerical illustrations, we will attempt a completely noninformative prior specification, namely we will set:

[TABLE]

We remark that $\lambda$ a-priori, follows the arcsine density, coinciding with the Jeffrey’s prior. The prior over the independent inverse variances $\tau_{i}$ is a vague gamma prior, being very close to a scale invariant prior.

**The multiple time-series approach: ** We will now set up a basic illustrative numerical example, concerning the support of the T-shifted to the past initial points for the Hénon map. We will consider the sequence of $r=250$ initial points ${\cal C}_{g}^{r}=\{(y_{j},y_{j+1}):1\leq j\leq r\}$, generated by means of the Hénon map

[TABLE]

with starting point $y_{-1:0}=(-0.61,1.37)$. Next, for each initial point $x_{-1:0}^{(j)}\in{\cal C}_{g}^{r}$, $j=1,\ldots,r$, we generate an associated $f$-noisy time-series $x_{1:n}^{(j)}$ of length $n=500$, via

[TABLE]

Now, conditional on the $r$ “observed” $f$-noisy time-series $x_{1:n}^{(j)}$, we sample from the posterior marginal distribution of the T-shifted to the past initial point variable $X_{-T-d+1:-T}$. We restrict the intermediate variables $X_{-1}$ and $X_{0}$ in the posterior (11) over the interval ${\cal X}=(-2,2)$, and we set a-priori $(X_{-3},X_{-2})\sim{\cal U}({\cal X}^{2})$. Finally, we consider the union of the supports of the variable $(X_{-T-d+1:-T}\,|\,x_{1:n}^{(j)})$ for $j=1,\ldots,r$, as an $(r,n)$-approximation to the perturbed global manifold of $g$; so that in the multiple time-series approach the SM-GSBR approximation is given by:

[TABLE]

In Fig. 4(a) we display the set of initial conditions ${\cal C}_{g}^{r}$ in black small circles. Together, we superimpose the set of points in ${\hat{W}}^{s}(T=0;r,n)$ in red. This is the union of the supports of the initial points $(X_{-1}^{(j)},X_{0}^{(j)})$ corresponding to the time-series $x_{1:n}^{(j)}$, for $j=1,\ldots,r$. Together we superimpose the set of points ${\hat{W}}^{s}(T=2;r,n)$ in gray. We can see how the increase of the backward prediction horizon to $T=2$, drives the supports of the posterior distributions along the stable direction. In Fig. 4(b)-(c) we provide enlargements of certain regions of the state space.

We remark that the quality of the initial points estimations $\hat{x}_{-1:0}^{(j)}$, associated with the time series $x_{1:n}^{(j)}$, for $j=1,\ldots,m$, depend not only on the perturbation realizations $e_{1:n}^{(j)}$, but also on the true position of the initial points in the state space.

In Fig. 5 we depict the sequence of points in ${\cal C}_{g}^{s}$ by small black circles. Together, we superimpose the predictive samples of the posterior marginal initial conditions of two special initial points that are in the neighborhood of primary HTs, for $T=0$ and $T=2$, in Fig. 5(a),(c),(e) and Fig. 5(b),(d),(f), respectively. Near the initial point $x_{-1:0}^{(87)}=(0.08,-1.32)$ (in blue) the stable manifold exhibits locally relatively high curvature, while near the initial point $x_{-1:0}^{(150)}=(1.71,-0.02)$ (in red) the curvature of the local stable manifold is close to zero. More specifically, a-posteriori the initial point variables have their support very close to the local stable direction, so that the structure of the support depends on the curvature of the local stable manifold. In Fig. 5(a)-(b) the posterior computations are based on the synthetic time-series produced by the non-Gaussian noise process in (16). In Fig. 5(c)-(d) the posterior computations are based on time series produced with ${\cal C}_{g}^{r}$ with Gaussian noise, of variance equal to the variance of the non-Gaussian noise. Finally, in Fig. 5(e)-(f) we have employed the corresponding deterministic time-series. Clearly, the most fuzzy results qualify the non Gaussian perturbations.

**The sliding window time-series approach: ** In practice, most of the times we have only one available data set $x_{1:n}$, generated by some unknown random dynamical system. In order to reveal the macroscopic structure of the stable manifold, we apply the T-backward prediction algorithm a multiple of times, each time based on a sliding window time-series subset of the observed time-series.

Therefore, after restricting the intermediate variable $X_{-T+1:0}$ in the posterior (11) over ${\cal X}^{T}$ and assigning to the variable $X_{-T-d+1:-T}$ the uniform prior over ${\cal X}^{d}$, conditionally on the $k$ sliding window time-series $x_{j:n-k+j}$ of length $n-k+1$, we sample from the posterior marginal distribution of the T-shifted to the past initial point variable $X_{-T-d+1:-T}$. So, based on a single observed time-series, we define the SM-GSBR approximation to the perturbed stable manifold of $g$ as the union of all the associated supports of the variables $(X_{-T-d+1:-T}\,|\,x_{j:n-k+j})$ for all $j=1,\ldots,k$, namely

[TABLE]

Regarding the choice of the prediction horizon $T$, we have observed that small values of $T=2,3,4$, are adequate.

4 Simulation results

In this section, we will provide numerical illustrations of the SM-GSBR algorithm, in approximating the global stable manifold of the Hénon and the Dual-Hénon maps, using single synthetic time series. Moreover, we will apply the SM-GSBR model on the orbits of a Hénon type noninvertible polynomial map.

In all our numerical experiments, we use the noninformative prior specifications (15) and the model polynomial map in (12). We will set the restriction-prior interval $\cal X$ in posterior (11) on a case-by-case basis. The SM-GSBR Gibbs sampler executions ran for $20\times 10^{4}$ iterations, leaving the first $5\times 10^{4}$ samples as a burn-in period. Thinning has been performed in every 150 iterations as a means for reducing the correlation between sampled values.

4.1 The Hénon map

**The zero noise limit case: ** In this case, the observed time-series $x_{1:n}$ is deterministic of size $n=2000$, coming from the classical Hénon map

[TABLE]

with initial point at $(-1,0.5)$. In fact, we consider the deterministically generated time-series $x_{1:n}$, as being additively corrupted by a close to zero intensity, unknown dynamical noise process $f^{*}$. For example, we could assume the existence of such a zero noise limit component $e_{i}\sim f^{*}$ as coming from the errors produced at each iteration of the map, caused by the finite precision of the floating point computer arithmetic.

In this case, we have fixed the time reversed prediction horizon to $T=3$ with a $k=500$ sliding window time-series. We restrict the intermediate variables $X_{-2},X_{-1}$ and $X_{0}$ over the interval ${\cal X}=(-3,3)$, and we set a-priori $(X_{-4},X_{-3})\sim{\cal U}({\cal X}^{2})$. The resulting SM-GSBR approximation of the stable manifold in (17) is given by

[TABLE]

In Fig. 6 we present the $\mathbb{R}^{2}$-embedded data points lying close to the Hénon attractor in black, superimposed with the sampled values in the set ${\hat{W}}_{2000}^{s}$ in red.

It is evident that the union of the supports of the joint posterior marginals forms a stochastic approximation of the associated global stable manifold of the deterministic map, as shown in Fig. 1.

In order to further investigate the quality of the stochastic approximation, we show in Fig. 7 data points and sampled values from the joint posterior marginals over the rectangle $\left[0.65,1.25\right]\times\left[-0.5,0.5\right]$.

**Non-Gaussian noise: ** Here we generate a time series $x_{1:n}$ of length $n=1000$, using the Hénon map in (18), with the same initial point, yet with an additive stochastic component which is given by

[TABLE]

with $\sigma_{1}^{2}=10^{-6}$ and $\sigma_{2}^{2}=10^{-3}$.

We set $T=3$, the number of sliding time-series to $k=250$, and the restriction-prior interval to ${\cal X}=(-6,2)$. We denote the SM-GSBR approximation of the global stable manifold by ${\hat{W}}_{1000}^{s}=\cup_{1\leq j\leq 250}\,{\rm supp}\left(X_{-4},X_{-3}\,|\,x_{j:j+750}\right).$

We note that due to the presence of corrupting noise, the approximation of the global stable manifold is subjected to a blurring effect. However, with the proposed SM-GSBR method, we are still able to gain insight about the qualitative characteristics of the dynamical behavior of the underlying system responsible for the observed data set.

In such cases, instead of presenting the approximation to the stable manifold as the union of supports of posterior marginals, it will be more instructive to present the approximation, with respect to higher posterior density regions (HPDRs) supported over the truncating set $\cal X$. More specifically, we apply a high analysis grid over $\cal X$, and we color-code regions according to the frequency of the predictive samples in ${\hat{W}}_{1000}^{s}$.

In Fig. 8 we present the $\mathbb{R}^{2}$-embedded noisy time-series $x_{1:n}$ of length $n=1000$ in black. Together, we superimpose the HPDR set based on the SM-GSBR approximation.

4.2 The Dual-Hénon map

**The zero noise limit case: ** In this subsection we will consider the Dual-Hénon map

[TABLE]

which is a special case of the Generalized Hénon map in (2). These is the case in which, two symmetric with respect to the origin, isolated Hénon-like attractors ${\cal A}^{+}$ and ${\cal A}^{-}$ coexist (see Fig. 5.3), located in the first and third quadrant, respectively.

We generate two deterministic orbits $x_{1:n}^{+}$ and $x_{1:n}^{-}$ of length $n=1000$, by iterating the initial points $p^{+}=(1,0.5)$ and $p^{-}=(0.5,1)$, lying on the basins of attraction of ${\cal A}^{+}$ and ${\cal A}^{-}$, respectively.

We use the two orbits as observed data sets, and to each orbit, we apply the SM-GSBR sampler with reversed time horizon $T=2$, with $k=500$ sliding time-series. We set the restriction-prior interval to ${\cal X}=(-5,5)$, so that a-priori $(x_{-3},x_{-2})\sim{\cal U}({\cal X}^{2})$, and the variable $(X_{-1},X_{0})$ is truncated over the square ${\cal X}^{2}$. The two corresponding SM-GSBR approximations are given by

[TABLE]

In Fig.9 we present the $\mathbb{R}^{2}$-embedded data sets $x_{1:n}^{+}$ and $x_{1:n}^{-}$ (in black), lying close to the coexisting strange attractors ${\cal A}^{+}$ and ${\cal A}^{-}$, respectively. Together we superimpose the SM-GSBR approximations ${\hat{W}}_{1000}^{s\,+}$ and ${\hat{W}}_{1000}^{s\,-}$ of the two stable manifolds emanating from the two coexisting saddle fixed points (see Fig. 3), in red and gray, respectively.

We remark, that the two bivariate densities $\Pi(x_{-3},x_{-2}\,|\,x_{1,n}^{\pm})$, for computational simplicity have common support $\cal X$. This is the reason why a very sparse cloud of points (in red) belonging to ${\hat{W}}_{1000}^{s\,+}$ can be observed inside the basin of attraction of ${\cal A}^{-}$ and vise-versa.

At this point, we would like to emphasize the fact that a-priori, there was no information on the control parameters of the map or the location of the saddle fixed points. Our stochastic approximation, given the general functional representation of the deterministic part is completely driven by the observed time-series.

**The case of an impulsive noise process: ** Here we mostly elaborate on the efficiency of the proposed SM-GSBR model. We generate a short time-series $x_{1:n}^{+}$ of length $n=500$, via the Dual-Hénon map in (19), yet influenced by the additive impulsive stochastic component

[TABLE]

with $\sigma_{1}^{2}=10^{-7}$ and $\sigma_{2}^{2}=10^{-2}$. The initial condition of the noisy orbit $x_{1:n}^{+}$ is given by $p^{+}=(1,0.5)$.

We set $T=2$, with the number of sliding time-series decreased to $k=150$. At the same time we expand the restriction-prior interval to ${\cal X}=(-10,10)$. We denote the associated SM-GSBR approximation by

[TABLE]

In Fig. 10 we present the $\mathbb{R}^{2}$-embedded noisy time-series $x_{1:n}^{+}$ of length $n=500$ in black. Together, we superimpose the HPDI set corresponding to the SM-GSBR approximation ${\hat{W}}_{500}^{s\,+}$.

We remark that although the noise processes was non-Gaussian and the observed time-series was short, we were able to produce an adequate approximation to the stable manifold of the saddle fixed point located on the first quadrant.

4.3 A noninvertible 2-d quadratic map

In this subsection we illustrate the performance of the SM-GSBR model based on nonlinear time-series data coming from a noninvertible polynomial map. For the generation of the synthetic time-series, we make use of the noninvertible quadratic map, given by:

[TABLE]

For a detailed analysis of the dynamical behavior of two–dimensional noninvertible quadratic maps we refer to [6].

Since there exist orbit points with more than one preimages, we extend the definition of the stable and unstable manifolds [29] to the following:

[TABLE]

and,

[TABLE]

with $y$ being a fixed point for the map $g$. Whenever the invertibility of $g$ is not guaranteed, usually ${\cal S}^{s}$ and ${\cal S}^{u}$ are being referred to as the stable and unstable invariant sets, respectively. The global invariant sets are no longer guaranteed to be manifolds and due to the multiple preimages, they can have intersections, being non–smooth or being totally disconnected [10]. The approximation of the stable manifold in such case is useful, because it gives us information on the complicated basins of attraction that the multiple preimages created [7]. The simplest case of noninvertible quadratic maps are of the $(Z_{0}-Z_{2})$ type, meaning that there exist points in the state space with only two or zero preimages, defining the two mutually exclusive basins of attraction.

**The zero noise limit case: ** We generate the deterministic orbit $x_{1:n}$ of length 2000, via the map in (21), with initial point at $(0.5,1.5)$. We aim to approximate the global stable set of the noninvertible map in (21).

We set $T=3$ with $k=500$ sliding time-series and the restriction-prior interval to ${\cal X}=(-2,2)$. We denote the SM-GSBR approximation of the stable set over $\cal X$, by

[TABLE]

In Fig. 11 we present the deterministic time-series $x_{1:n}$ in black, and the union of the predictive sample in ${\hat{W}}_{2000}^{s}$ in red.

**The case of an impulsive noise process: ** We generate a time-series $x_{1:n}$ of length $n=1000$, via the noninvertible map in (21), influenced by the additive impulsive stochastic component in (20).

We set $T=2$ with $k=500$ sliding time-series, and we set the restriction-prior interval to ${\cal X}=(-4,4)$. We denote the SM-GSBR approximation of the stable set over $\cal X$, by

[TABLE]

In Fig. 12 we present the $\mathbb{R}^{2}$-embedded noisy time-series $x_{1:n}$, together with the approximation of the stable set based on the noisy time-series $x_{1:n}$, via an HPDR color map. We remark that the qualitative characteristics persist even when the approximation is based to short time-series contaminated by dynamical noise.

5 Discussion

We have presented a novel approach for the stochastic approximation of the global stable manifold, based on MCMC methods. Our approach is Bayesian, modeling the global stable manifold as the union of supports of T-shifted to the past initial condition variables given appropriate sliding window subsets of the time-series observations. By proper applications of the proposed SM-GSBR model, we were able to approximate the perturbed global stable manifold, and gain insight about the qualitative dynamical behavior in the zero noise limit, having as reference the associated deterministic system.

Our proposed SM-GSBR sampler is flexible because it requites no prior knowledge on the type of the noise process perturbing the observed time-series. The reason for that, is that over the additive errors, we apply a Geometric Stick Breaking mixture process prior, which is supported over the space of symmetric zero-mean densities. In fact the SM-GSBR sampler requires no prior information except for the model space required for the functional representation of the unknown deterministic part $g$ of the underlying map. Because we do not use any representation of the preimages of $g$ the SM-GSBR sampler is directly applicable to noninvertible maps.

A more versatile version of the SM-GSBR sampler, would involve relaxing the assumption of the functional representation of the unknown deterministic part $g$. For example by applying over $g$ a Gaussian process [28] prior.

Appendix A Proof of Proposition 1

To augment the posterior in (6), we will group the ${\cal V}_{T}$ auxiliary variables into $d_{-T+1:0},N_{-T+1:0}$ and $d_{1:n},N_{1:n}$. We will also group the variables $\lambda,\tau_{1:\infty}$ and $\vartheta$, into the variable $\nu=(\lambda,\tau_{1:\infty},\vartheta)$. The variables $x_{-T+1:0},\,d_{-T+1:0}$ and $N_{-T+1:0}$ to the variable $y_{-T+1:0}=(x_{-T+1:0},d_{-T+1:0},N_{-T+1:0})$, and the variables $x_{1:n},\,d_{1:n}$ and $N_{1:n}$ to the variable $y_{1:n}=(x_{1:n},d_{1:n},N_{1:n})$. To truncate the intermediate variables between the $T$-shifted to the past initial condition and the observed time-series, we condition the posterior density on the event ${\cal R}_{T}$. Then, using Bayes’ theorem, the posterior becomes

[TABLE]

Finally, the desired result comes from equation (8). $\square$