PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
 
Commun Nonlinear Sci Numer Simul. 2018 January; 54: 248–266.
PMCID: PMC5589148

A matrix-based approach to solving the inverse Frobenius–Perron problem using sequences of density functions of stochastically perturbed dynamical systems

Abstract

The paper introduces a matrix-based approach to estimate the unique one-dimensional discrete-time dynamical system that generated a given sequence of probability density functions whilst subjected to an additive stochastic perturbation with known density.

Keywords: Nonlinear systems, Chaotic maps, Probability density functions, Inverse Frobenius–Perron problem

1. Introduction

A dynamical system, whose evolution is completely dictated by deterministic equations, can under certain conditions exhibit chaotic behavior and generate a density of states [1]. Chaotic behavior has been observed in many real-world systems including biological, physical and economic systems [2], [3], [4]. The simplest dynamical systems that exhibit chaos are one-dimensional maps. Such one-dimensional discrete dynamical systems are used to describe the evolution of many real-world systems including olfactory systems [5], electrical circuits [6], communication networks [7], rotary drills [8], chemical reactions [9] and the heart [10].

An important challenge is to develop such model from experimental observations [11], [12], [13]. Conventional approaches [14], [15], [16] rely on time series data. However, often it is not possible to measure point trajectories. For example, particle image velocimetry, a technique used to generate instantaneous velocity measurements in fluids, identifies individual tracer particles in consecutive images captured at high speeds but cannot resolve their individual orbits [17]. In such cases, in the absence of individual point trajectories, it is desirable to determine the underlying dynamical system that generated the observed density functions. Given a non-singular transformation, the evolution of an initial density function under the action of the transformation is described by the Frobenius–Perron operator associated with the transformation [18]. The fixed point of such an operator represents the invariant density under the transformation. The problem of inferring a dynamical system, whose invariant density function is given, is known as the Inverse Frobenius–Perron problem [19].

In general, solving the inverse problem involves deriving a finite-dimensional representation of the operator, which is then used to construct the dynamical system. Ulam conjectured [20] that a general infinite-dimensional Frobenius–Perron operator can be approximated by a finite rank Markov operator. For one-dimensional transformations Li [21] has shown that given a sequence of piecewise constant approximations Pn of the Frobenius–Perron operator P, the corresponding sequence of fixed points fn of Pn converge to the invariant density (i.e. fixed point) of the operator, thus proving Ulam's conjecture. In this context, the problem of determining the dynamical system that corresponds to the finite dimensional approximation of the Frobenius–Perron operator is also known as the inverse Ulam problem [22], [23].

A numerical algorithm to determine a one-dimensional transformation given the invariant density function was proposed in [19]. The algorithm however does not provide an explicit relationship between the invariant density of the one-dimensional map and the map itself. In [24] a graph-theoretic approach is introduced to construct piecewise-linear transformations that possess piecewise-constant invariant density functions that have value 0 in all relative minima points.

A generalization of these methods is presented in [25] which introduces a relationship between any arbitrary piecewise density function and a semi-Markov piecewise linear transformation defined over a partition of the interval of interest. This forms the basis for a matrix-based method to reconstruct a 3-band transformation, a special class of semi-Markov transformations, which has a given piecewise-constant density function as invariant density. The inverse problem was studied in [26] for a class of symmetric maps that have invariant symmetric Beta density functions and the unique solution can be achieved under given symmetry constraints. This method was generalized in [27] which considers a broader class of continuous unimodal maps for which each branch of the map covers the complete interval and the invariant densities are a asymmetric beta functions. Given arbitrary invariant densities similar approaches were proposed for identifying the maps with specified forms: two types of one-dimensional symmetric maps [28], smooth chaotic map with closed form [29], [30], multi-branches complete chaotic map [31]. Problems of synthesizing one-dimensional maps with prescribed invariant density function or autocorrelation function were considered in [32], [33]. Using positive matrix theory an approach to synthesizing chaotic maps with arbitrary piecewise constant invariant densities and arbitrary mixing properties was developed in [34]. This method was further extended to synthesizing dynamical systems with desired statistical properties [35], developing communication networks [36] and designing randomly switched chaotic maps and two-dimension chaotic maps used for image generation [37]. The global, open-loop strategy to control chaos proposed in [22], [23] is formulated as an inverse Frobenius–Perron problem. The aim is to perturb the original dynamical system to achieve the desired invariant measure. This reduces to the problem of finding a perturbation of the original Frobenius–Perron matrix to achieve the target invariant density function and then solving the inverse Ulam problem to determine the perturbed dynamical system.

In general, the solution to the inverse Frobenius–Perron problem is not unique. Different transformations exhibiting strikingly different dynamics may share the same invariant density functions. Additional limiting assumptions or constraints are required to ensure uniqueness of the solution [26], [27], [28], [29], [30], [31], [32], [33], [34]. In [38] a new method is proposed to solve the inverse Frobenius–Perron problem based not only on the invariant density function but also on sequences of probability density functions generated by the transformation, which ensures uniqueness of the solution. The method has been shown to be quite robust to noise. For small levels of noise this is indeed expected in light of the convergence results for noise perturbed systems established by Bollt et al. [39]. However, the accuracy of the reconstruction starts to deteriorate significantly above a certain level of noise. For large levels of noise, the approximation errors can be drastically reduced by taking into account the density function of the noise, which is often known a priori or can be estimated.

This paper proposes a new method to estimate the piecewise linear and expanding semi-Markov transformation that generated a temporal sequences of probability density functions whilst subjected to constantly applied stochastic perturbations. The method is extended to more general nonlinear transformations that can be approximated arbitrarily close by piecewise linear functions.

To differentiate with the normal deterministic inverse problem, we call this the inverse stochastic Frobenius–Perron problem. The emphasis here is on recovering the unknown transformation that generated the sequence of densities rather than one of the many possible transformation that share the same invariant density function.

To this end, we formulate the matrix representation of the transfer operator associated with the stochastically perturbed system in terms of the Frobenius–Perron matrix associated with the unperturbed system that we aim to estimate. This representations forms the basis for a the proposed algorithm to estimate the ‘unperturbed’ Frobenius–Perron matrix from sequences of probability density functions generated by the unknown, stochastically perturbed dynamical system, under the assumption that the density function of the perturbation is known. For general nonlinear transformations, we present a practical method to solve the inverse Ulam problem, which allows determining the sign of the derivative for each interval of the partition.

Whilst the sign of the derivative is not important when the goal is to determine a transformation that has a given invariant density function [22], this step is crucial if the aim is to reconstruct/approximate the true dynamical system that generated the data. We demonstrate that the proposed approach can reconstruct the underlying dynamical system that is subject to stochastic perturbations.

This paper is organized as follows. Section 2 introduces the inverse stochastic Frobenius–Perron problem. A matrix approximation of the transfer operator associated with the stochastically perturbed transformation is derived in Section 3. Section 4 introduces a methodology for reconstructing piecewise-linear semi-Markov transformations subject to stochastic perturbations, from sequences of density functions. The approach is extended in Section 5 to general nonlinear maps. Section 6 presents two numerical simulation examples that demonstrate the significant improvement in reconstruction accuracy achieved by the proposed algorithm that incorporates a priori knowledge of the noise in the reconstruction of the unknown transformation. Conclusions are given in Section 7.

2. Description of the inverse problem

Let R=[0,b], B be a Borel σ-algebra of subsets in R, and μ denote the normalized Lebesgue measure on R. Let S: RRbe a measurable, non-singular transformation, that is, μ(S1(A))B for any AB and μ(S1(A))=0 for all AB with μ(A)=0. If xn is a random variable on R having the probability density function fnD(R,B,μ), D={fL1(R,B,μ):f0,f1=1}, such that

Prob{xnA}=Afndμ.
(1)

It follows that xn+1 given by

xn+1=S(xn),
(2)

is distributed according to the probability density function fn+1=PSfn, where PS: L1(R) → L1(R), defined by

APSfndμ=S1(A)fndμ,
(3)

is the Frobenius–Perron operator [1] associated with the unperturbed transformation S.

If A = [a, x],PS can be written explicitly as

fn+1(x)=PSfn(x)=ddxS1([a,x])fndμ,
(4)

Let ={R1,R2,,RN} be a partition of R into intervals, and int(Ri)int(Rj)= if ij. Assuming that S is piecewise monotonic and expanding [18],

PSfn(x)=i=1Nfn(Si1(x))|S(Si1(x))|χS(Ri)(x),
(5)

where Si is the monotonic restriction of S on the interval Ri.

A more complicated situation arises when the dynamical system is subjected to an additive random perturbation [1] such that

xn+1=S¯(xn,ξn)=S(xn)+ξn,
(6)

where S: RR is a given transformation and ξ1,ξ2,,ξn are independent random variables. The ‘stochastic’ Frobenius–Perron operator P¯ corresponding to the perturbed dynamical systems is defined by [1], [39]

P¯f(x)=Rτ(x,y)f(y)dy,
(7)

where τ(x,y)=g(xS(y))is a stochastic kernel, satisfying τ(x, y) > 0, and Rτ(x,y)=1.

Here, we consider the deterministic system with constantly applied stochastic perturbations

xn+1=S(xn)+ξn(modb),
(8)

where S: [0, b] → [0, b] is a piecewise monotonic and expanding transformation, ξn is an independent random variable with a probability density function g that has compact support on [ɛ,ɛ], i.e. ξn is bounded in [ɛ,ɛ], epsilonb.

For an arbitrary Borel set B[subset or is implied by][0, b], the probability of xn+1 falling into B is given by

Prob{xn+1B}=R[ɛ,ɛ]fn(x)g(ξ)dxdξ,
(9)

for xn+1=S¯(xn,ξn), where S¯(xn,ξn)=S(xn)+ξn(modb).

Let y=S(x)+ξ(modb). It follows that

y=S(x)+ξbχ(b,b+ɛ](S(x)+ξ)+bχ[ɛ,0)(S(x)+ξ),
(10)

and

ξ=yS(x)+bχ(b,ɛb](yS(x))bχ[bɛ,b)(yS(x)),
(11)

then, (9) is rewritten as

Prob{xn+1B}=Bfn+1(x)dx=BRfn(x)g(yS(x)+bχ(b,ɛb](yS(x))bχ[bɛ,b)(yS(x)))dydx,
(12)

where fn+1 is the probability density function of xn+1. The stochastic Frobenius–Perron operator P¯:L1(R)L1(R), associated with the perturbed transformation (8), is then defined by

P¯fn(x)=fn+1(x)=Rfn(z)g(xS(z)+bχ(b,ɛb](xS(z))bχ[bɛ,b)(xS(z)))dz.
(13)

It is easy to see that for any ξn there are N1N disjoint intervals {Iξni}i=1N1 such that for xnIξni, i=1,,N1, S(xn) is monotonic and S(xn)+ξn[b,b+ɛ] or S(xn)+ξn[ɛ,0] (i.e. maps outside the interval [0, b]), and N2N disjoint intervals {I¯ξni}i=1N2 such that for xnI¯ξnj, j=1,,N2, S(xn) is monotonic and S(xn)+ξn[0,b].

We have I¯ξnjI¯ξni=, ∀iN1, jN2, and (i=1N1Iξni)(j=1N2I¯ξnj)=[0,b]. For each integers iN1 and jN2, there exist unique integers α(i) ≠ β(j) ≤ N such that IξniRα(i) and I¯ξnjRβ(j). From (13) it follows that

P¯fn(x)=i=1N1Iξnifn(z)g(xS(z)+bχ(b,ɛb](xS(z))bχ[bɛ,b)(xS(z)))dz+j=1N2I¯ξnjfn(z)g(xS(z))dz.
(14)

Substituting y=S(z) gives

P¯fn(x)=i=1N1Iξnifn(S1(y))S(S1(y))χS(Iξni)(y)g(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))dy+j=1N2I¯ξnjfn(S1(y))S(S1(y))χS(I¯ξnj)(y)g(xy)dy.
(15)

We have

i=1N1Iξnifn(S1(y))S(S1(y))χS(Iξni)(y)dy+j=1N2I¯ξnjfn(S1(y))S(S1(y))χS(I¯ξnj)(y)dy=k=1NI˜ξnkfn(S1(y))S(S1(y))χS(I˜ξnk)(y)dy.
(16)

where I˜ξnk{Iξn1,,IξnN1,I¯ξn1,,I¯ξnN2}. Since {S(Iξni)}i=1,,N1{S(I¯ξnj)}i=1,,N2=R, the following equality holds

k=1NI˜ξnk[fn(S1(y))S(S1(y))χS(I˜ξnk)(y)]dy=Ri=1N[fn(S1(y))S(S1(y))χS(Ri)(y)]dy.
(17)

It follows that

P¯f(x)=Rg(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))·i=1N[fn(S1(y))S(S1(y))χS(Ri)(y)]dy.
(18)

So that

P¯f(x)=Rg(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))PSfn(y)dy.
(19)

Eq. (19) provides the link between the operator P¯corresponding to the randomly perturbed dynamical system (8) and the Frobenius–Perron operator PS associated to the noise-free system. The inverse problem is formulated as follows.

Let {x0,ij}i,j=1θ,K and {x1,ij}i,j=1θ,K be K sets of initial and final state observations, respectively, such that

x1,ij=S(x0,ij)+ξ0j(modb),j=1,,K
(20)

It is assumed that the measurement system does not allow associating an initial state x0,ij with its image x1,ij under the transformation. The inverse problem considered here is to determine the transformation S in (8) given the noise density function g and probability density functions {f0j}j=1K, {f1j}j=1K associated with the initial states{x0,ij}i,j=1θ,K and final states {x1,ij}i,j=1θ,K , that isf1j=P¯f0j, j=1,,K, whereP¯ is the transfer operator associated with perturbed transformation (8).

3. A matrix representation of the transfer operator P¯

Let S be a piecewise linear and expanding semi-Markov transformation over the N-interval partition, ={R1,R2,,RN}.

Definition 1

A transformation S: RR is said to be semi-Markov with respect to the partition R (or R-semi-Markov) if there exist disjoint intervals Qj(i) so that Ri=k=1p(i)Qk(i), i=1,,N, the restriction of S to Qk(i), denoted S|Qk(i), is monotonic and S(Qk(i)). [25]

The restriction S|Riis a homeomorphism from Ri to a union of intervals of R

k=1p(i)Rr(i,k)=k=1p(i)S(Qk(i)),
(21)

where Rr(i,k)=S(Qk(i)), Qk(i)=[qk1(i),qk(i)], i=1,,N,k=1,,p(i) and p(i) denotes the number of disjoint subintervals Qk(i) corresponding to Ri.

Let fn be a piecewise constant function over the partition R such that fn(x)=i=1NwiχRi(x). Its image under transformation PSfn is also a piecewise constant function over R [18] such that PSfn(x)=i=1NuiχRi(x) and the Frobenius–Perron operator can be represented by a finite-dimensional matrix

PSfn(x)=j=1N(i=1N(wimi,j))χRj(x),
(22)

where M=(mi,j)1i,jN is the Frobenius–Perron matrix induced by S with entries given by

mi,j={|(S|Qj(i))|1,ifS(Qk(i))=Rj;0,otherwise.
(23)

From (22) it follows that

uj=i=1Nwimi,j,
(24)

for j=1,,N.

Let ¯={R¯1,R¯2,,R¯N} be a regular partition of R into N equal sized intervals. By integrating (19) over an interval R¯k¯ gives

R¯kP¯fn(x)dx=R¯kIg(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))·PSfn(y)dydx
(25)

Consider the following approximation

P¯fn(x)=pN(x)+qN(x)=k=1NvkχR¯k(x)+qN(x),
(26)

where pN(x)=P¯Nfn(x) is the orthogonal projection of P¯fn in L1 on the finite-dimensional space spanned by, qN(x) is the orthogonal complement in L1 and

vk=1λ(R¯k)R¯kIg(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))·PSfn(y)dydx,
(27)

where λ(R¯k) is the Lebesgue measure on R¯k. Clearly, qN(x) → 0 as N+.

It follows that

vk=1λ(R¯k)R¯kj=1N[Rjg(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))dy·uj]dx=Nbj=1N[R¯kRjg(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))dydx·uj].
(28)

Let D=(dk,j)1kN;1jN be the N × N matrix with entries given by

dk,j=NbR¯kRjg(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))dydx.
(29)

Substituting (24), (29) in (28) leads to

vfn+1=wfn·M·D=wfn·Q,
(30)

where wfn=[w1,,wN], vfn+1=[v1,,vN] are the coefficient vectors associated with the piecewise constant density functions fn and P¯fn respectively and Q=M·D is the matrix approximation of the operator P¯. Eq. (30) maps a piecewise-constant density function over the N-dimensional partition, which in general is non-uniform, to a piecewise-constant density function over a uniform N-dimensional partition. Eq. (30) is the basis for the new algorithm to reconstruct the transformation S given pairs of successive density functions generated by the stochastically perturbed transformation.

In practice, we can chose a finer N1-interval partition ¯, N1 >> N. For example, we can construct ¯ as a refinement of the partition R such that the cut points of partition R are a subset of cut-points associated with ¯. This leads to an alternative formulation of (30) where both the initial and final densities are defined over the same partition. Given an initial piecewise density function f over the partition¯, the matrix approximation can then be used to compute a sequence of successive iterations by the corresponding finite-dimensional approximation of the stochastic Frobenius–Perron operator.

4. Solving the inverse stochastic Frobenius–Perron problem for piecewise linear semi-Markov transformations

This section presents an approach to solving the inverse stochastic Frobenius–Perron problem, under the assumption that S: [0, b] → [0, b] is a piecewise linear semi-Markov transformation over a partition R,

={R1,R2,,RN}={[0,a1],(a1,a2],,(aN1,aN]},
(31)

aN=b, which is assumed to be known. In what follows we assume that ¯ is defined as the uniform partition of dimension N of [0, b].

The main steps of the approach are summarized below:

  • Step 1: Given the observations Xt={xt,j}j=1θ, t = 0,…,T, estimate the coordinate vectors wft=[wt,1,,wt,N] and vft+1=[vt+1,1,,vt+1,N], t = 0,…,T−1 corresponding to the piecewise constant density functions ft(x) over R and ft+1(x)=P¯Sft(x) over ¯, respectively. Compute the matrix D defined in (29).
  • Step 2: Estimate M, the matrix representation of the Frobenius–Perron operator PS associated with the deterministic transformation S.
  • Step 3: Construct the piecewise linear semi-Markov transformation over R.

    These steps are described below in more detail.

4.1. Step 1: estimate w and v and compute d

Let f0(x) be an initial density function that is piecewise constant on the partition={R1,,RN}.

f0(x)=i=1Nw0,iχRi(x),
(32)

where the coefficients satisfy i=1Nw0,iλ(Ri)=1. Let X0={x0,j}j=1θ be the set of initial states obtained by sampling f0(x). The states Xt={xt,j}j=1θ at a given sampling time t > 0 are assumed to be generated by applying t times the process defined in (8), where Ω={ξi}i=1θ are generated by sampling g(ξ).

The density function ft(x) on ¯ associated with the states Xt is given by

ft(x)=i=1Nvt,iχR¯i(x),
(33)

where the coefficients vt,i=1λ(R¯i)·θj=1θχR¯i(xt,j)=Nbθj=1θχR¯i(xt,j). In practice the densities ft(x) are estimated directly from observations.

We define the following matrices

W0=[wf0wf1wfT1]=[w0,1w0,2w0,Nw1,1w1,2w2,NwT1,1wT1,2wT1,N],
(34)

and

V=[vf1vf2vfT]=[v1,1v1,2v1,Hv2,1v2,2v2,HvT,1vT,2vT,H].
(35)

The matrix D is obtained by numerical integration of (29).

4.2. Step 2: estimate the Frobenius–Perron matrix m

This is carried out in two stages. Firstly, the coordinate vector corresponding to the piecewise constant densities PSft(x)=i=1Nut,iχRi(x) over the partition R are obtained by solving the following constrained optimization problem

min0{ut,j}t=0,,T1,j=1,,NNVY·DF,
(36)

subject to

i=1Nut,iλ(Ri)=1,fort=0,,T1.
(37)

where

Y=[uf0uf1ufT1]=[u0,1u0,2u0,Nu1,1u1,2u2,NuT1,1uT1,2uT1,N],
(38)

and || · ||F denotes the Frobenius norm.

In the second stage, the matrix representation of the Frobenius–Perron operator associated with the unperturbed transformation S is obtained as a solution to the following constrained optimization problem

min{mi,j}i,j=1N0YW0M||F,
(39)

subject to

j=1Nmi,jλ(Rj)=λ(Ri),fori=1,,N.
(40)

In the following it is shown that the matrices Φ1=DD and Φ2=W0W0 are non-singular, which ensures uniqueness of solutions.

Proposition 1

For a piecewise linear R-semi-Markov transformation S subjected to additive perturbation, where=¯ is an N-dimensional regular partition of [0, b], the matrixΦ1=DD is non-singular.

Proof

For ¯={R¯k=Rj,k=1,,N}, λ(Rj)=λ(R¯k), Rj [set membership] R, Rk¯, for j,k=1,,N,

dk,j=NbRkRjg(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))dydx.
(41)

The matrix D satisfies that

dk,j=dj,k,dj,j=dk,k,j,k{1,,N}
(42)

and dk,j=dk+1,j+1, for 1k,jN1.

Let di=d1,i, for i=1,,N. The matrix D is decomposed into two triangle matrices as follows

D=Dl+Du,
(43)

where

Dl=[d1/2d2d1/20d3d2d1/2dNdN1dN2d1/2],
(44)
Du=[d1/2d2d3dNd1/2d2dN1d1/2dN20d1/2].
(45)

Then, det(Dl)=det(Du)=d1N2N. According to the Minkowski determinant theorem, Dl and Du are both non-negative, thendet(D)det(Dl)+det(Du)>0. Hence, D and Φ1 are non-singular, and this completes the proof.

Theorem 1

Let S: [0, b] → [0, b] be a piecewise linear R-semi-Markov transformation subjected to additive noise,=¯. Then the matrix Q representing the transfer operator associated with the noisy dynamical system has 1 as the eigenvalue of maximum modulus and also has the unique eigenvalue of modulus 1.

Proof

For =¯, the matrix Q=M·D is square. Let Q=(qi,j)1i,jN where

qi,j=k=1N(mi,kdj,k).
(46)

The sum of ith row of Q is given by

j=1Nqi,j=j=1N([mi,1mi,kmi,N][dj1dj,kdj,N])=[mi,1mi,kmi,N][d11++dj,1++dN,1d1,k++dj,k++dN,kd1,N++dj,N++dN,N].
(47)

The column sum of kth column of D is given by

i=1Ndi,k=i=1NNbRiRjg(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))dydx=NbIRjg(xy+bχ(b,ɛb](xy)bχ[bɛ,b)(xy))dydx=Nb·bN=1.
(48)

It follows that

j=1Nqi,j=j=1Nmi,j.
(49)

For a regular partition R, j=1Nmi,j=k=1p(i)Qk(i)λ(Ri)=1, then j=1Nqi,j=1. Thus, Q is row stochastic. Because Q is also a positive matrix, matrix Q has 1 as the eigenvalue of maximum modulus, and the algebraic and geometric multiplicities of this eigenvalue are 1. This concludes the proof.

Remark

. From Theorem 1 it follows that f*=P¯Nf* has a unique, no-trivial solution, where P¯N is the N-dimensional approximation of the stochastic Frobenius–Perron operator associated with the noise perturbed piecewise linear R-semi-Markov transformation (8).

Proposition 2

A noise perturbed piecewise linear R-semi-Markov transformation S can be uniquely reconstructed given N linearly independent piecewise constant density functionsf0i defined over a regular partition R, which correspond to initial states, and the piecewise constant densitiesf1i=P¯Nf0i.

Proof

Let

f0i(x)=j=1Nwi,j0χRj(x),i=1,,N
(50)

be a set of initial piecewise constant densities over a partition R and

f1i(x)=j=1Nwi,j1χRj(x),i=1,,N
(51)

be the piecewise constant densities, corresponding to the final states

From (30), we have

W1=W0·M·D,
(52)

where W0=[wf01wf02wf0N]=[w1,10w1,20w1,N0w2,10w2,20w2,N0wN,10wN,20wN,N0], W1=[wf11wf12wf1N]=[w1,11w1,21w1,N1w2,11w2,20w2,N0wN,11wN,21wN,N1].

Because the initial densities f0i, i=1,,N are linearly independent, the matrix W0 is non-singular. Moreover, from Proposition 1, D is non-singular and the Frobenius–Perron matrix is given by

M=W01W1D1.
(53)

4.3. Step 3: construct the piecewise linear semi-Markov transformation over R

It is assumed that each branch of the map, S|Qk(i), is monotonically increasing. The derivative of S|Qk(i) is 1/mi, j, the length of Qk(i) is given by

λ(Qk(i))=qk(i)qk1(i)=mi,jλ(Rj)
(54)

which allows computing iteratively qk(i) for each interval Ri starting with q0(i)=ai1.Then the map is given by

S|Qk(i)(x)=1mi,j(xqk1(i))+aj1,
(55)

for k=1,,p(i), j is the index of image Rj of Qk(i), i.e. S(Qk(i))=Rj, i=1,,N, j=1,,N, where mi, j ≠ 0.

5. Solving the inverse stochastic Frobenius–Perron problem for general nonlinear transformations

This section considers more general nonlinear maps that are not piecewise linear semi-Markov. Starting with Lasota and Yorke [40] who established the existence of invariant measure for piecewise monotonic transformations and Li [21] who proposed a numerical procedure to calculate the invariant density function corresponding to the invariant measure, the problem of approximating the invariant density of a transformation, which is closely linked with the problem of approximating the Frobenius–Perron operator and the transformation itself, has been studied by a number of authors [41], [42], [43], [44], [45]. In [43] Gora and Boyarsky approximated a nonsingular transformation S, that may have infinitely many pieces of monotonicity, by a sequence of piecewise linear functions Sn and shown that the invariant density of the map S can be approximated arbitrarily well by densities that are invariant under finite approximations Sn of S. In general, for continuous nonsingular transformation we have the following result [22], [46]

Theorem 2

Let S: RR be a continuous transformation and let {Sn}n ≥ 1be a sequence of transformations converging to S in the C0 topology. Let μn be a probabilistic measure invariant under Sn, n = 1,2,…. If μ is a weak-* limit point of the sequence {μn}n ≥ 1then μ is S-invariant.

This shows that the invariant densities fn* of successive piecewise linear approximations Sn of S converge in a weak sense to the invariant density f*of the original transformation as n→∞. This means that in practice we can estimate transition matrices to approximate arbitrarily well the Frobenius–Perron operator associated with S. A generalization of this result for dynamical systems subjected to additive noise is presented in [39].

Here, the goal is to estimate from sequences of density functions a piecewise linear semi-Markov approximation S^, defined over a uniform Markov partition ={R1,R2,,RN}of [0, b], λ(Ri)=b/N, of the unknown nonlinear map S: [0, b] → [0, b] subjected to stochastic perturbation. It is assumed that the nonlinear map S has an invariant density and that it can be approximated arbitrarily well by piecewise linear functions. Unlike the control problem studied in [33], here the challenge is to estimate the unknown nonlinear transformation that generated a sequence of density functions rather than one of the many possible perturbations of the original map, which yield a desired invariant density.

The proposed identification scheme for general nonlinear maps is summarized as follows:

Step 1: Given the observationsXn={xn,j}j=1θ, n=0,,T, estimate the coordinate vectors wft and vft+1 of the piecewise constant density ft(x) and ft+1(x) defined over a regular partition =¯ of size N. Compute D using the given probability density function of perturbation.

Step 2: Estimate the matrix Y corresponding to PSft(x)by solving the optimization problem (36). PSft(x) can be used to identify the Frobenius–Perron matrix associated with the unknown map;

Step 3: Identify a trial Frobenius–Perron matrix M^=(m^i,j)1i,jN , j=1Nm^i,j=1 by solving the constrained optimization problem (39). Since the map is continuous, this is then used to further determine the indices of consecutive positive entries of each row and in this way M^ is refined. Let Bi={rsi,rsi+1,,rei} be the set of column indices corresponding to consecutive positive entries of ith row and rmiBisatisfying

m^i,rmi=max{m^i,j}j=1N.
(56)

Therefore, the piecewise linear R-semi-Markov map associated with the refined Frobenius–Perron matrix M should satisfy that Rr(i,k)=S(Qk(i)), where k=1p(i)Rr(i,k) is the image of the interval Ri,i=1,,N, p(i)=reirsi+1 and r(i,k)Bi is the column index of a positive entry on the ith row of M satisfying

r(i,k+1)=r(i,k)+1.
(57)

for i=1,...,N,k=1,...,p(i)1.

Step 4: Solve the following optimization problem to determine the Frobenius–Perron matrix M

min{mi,j}0YW0M||F,
(58)

subject to

k=1p(i)mi,r(i,1)+k1=1.
(59)

and

{mi,r(i,k)>0,i=1,,N;mi,j=0,jr(i,k),
(60)

for k=1,...,p(i).

Step 5. Determine the monotonicity of each branch S|Qk(i). Let Ri=[ar(i,1)1,ar(i,p(i))] be the image of the interval Ri under the semi-Markov transformation S^ associated with the identified Frobenius–Perron matrix M. Denote ar(i,1)1 as the starting point of Rr(i, 1) mapped from the subinterval Q1(i), and ar(i, p(i)) as the end point of Rr(i, p(i)), the image of the subinterval Qp(i)(i). Let c¯i be the midpoint of the image Ri. The sign γ(i)of {S^(x)|Qk(i)}k=1p(i) is given by

γ(i)={1,ifc¯ic¯i1<0;1,ifc¯ic¯i10;γ(i1),ifc¯i=c¯i1,
(61)

for i=2,,N and γ(1)=γ(2).

Step 5. Construct semi-Markov map based on the Frobenius–Perron matrix M and the monotonicity of each branch. Given that the derivative of S|Qk(i) is 1/mi, j, the end point qk(i)of subinterval Qk(i) within Ri is given by

qk(i)={ai1+j=1kmi,r(i,j)λ(Rr(i,j)),ifγ(i)=+1;ai1+j=1kmi,r(i,p(i)k+1)λ(Rr(i,p(i)k+1)),ifγ(i)=1.
(62)

where k=1,...,p(i)1 and qp(i)(i)=ai. The piecewise linear semi-Markov transformation S^ on each subinterval Qj(i) is given by

S^|Qj(i)(x)={1mi,j(xqk1(i))+aj1,ifγ(i)=+1;1mi,j(xqk1(i))+aj,ifγ(i)=1.
(63)

for mi, j ≠ 0, i=1,,N,j=1,,N,k=1,...,p(i)1. A smooth nonlinear map is obtained by fitting a polynomial smoothing spline. Fig. 1 shows the construction of monotonically increasing and decreasing piecewise linear semi-Markov transformation and the resulting smooth continuous map.

Fig 1
Schematic diagram of construction of piecewise linear semi-Markov transformations for monotonically increasing (a) and decreasing (b) cases. The indices of positive entries of the Frobenius–Perron matrix are determined from a trial matrix. The ...

6. Numerical simulation studies

The proposed algorithms are demonstrated using simulated data generated by two chaotic maps.

6.1. Example A

Consider the noise perturbed dynamical system

xn+1=S(xn)+ξn,(mod1)
(64)

where {ξn} is white noise that follows a Gaussian distribution truncated to the range [ɛ,ɛ] where ɛ=0.02. The piecewise linear and expanding transformation S: [0, 1] → [0, 1] is defined by

S|Ri(x)=αi,jx+βi,j,
(65)

for i=1,,4, j=1,,4, over the partition ={Ri}i=14={[0,0.3],(0.3, 0.4], (0.4, 0.8], (0.8, 1]} where

(αi,j)1i,j4=[3.330.8313.333.3315.003.3310.0020.003.752.503.331.257.501.256.6710.00],(βi,j)1i,j4=[00.232.4004.500.773.1071.500.901.3306.000.755.739.00].

The graph of S is shown in Fig. 2.

Fig 2
Example A: Original piecewise linear transformation S.

A set of initial densities f0i, i=1,,4 over the partition R is shown in Fig. 3. These are used to generate the set of initial states X0i={x0,ji}j=15×103i=1,,4. The corresponding final states X1i={x1,ji}j=15×103, obtained by applying the map (64), were used to estimate piecewise constant densities f1i, i=1,,4 over a uniform partition {[0, 0.25], (0.25, 0.5], (0.5, 0.75], (0.75, 1]} of [0, 1].

Fig 3
Example A: The density functions f0i(x), f1i(x) and Psf0i(x).

For this partition, the matrix D calculated using (29)

D=[0.9921000.00800.20000.40000.4000000100.007900.20000.7920],
(66)

is non-singular.

The two stage approach detailed in Section 4, was used to estimate the matrix representation of the Frobenius–Perron operator associated with the deterministic transformation

M=[0.30661.20940.07370.28800.06850.30060.09890.04910.27260.40910.29700.79260.13630.82470.14190.0994].
(67)

Specifically, the constrained optimization problems (36), (37) and (39), (40) were solved using the lsqlin function in the Matlab Optimization Toolbox.

The reconstructed piecewise linear map S^ is shown in Fig. 4. The estimated coefficients of the identified piecewise linear semi-Markov transformation S^|Ri(x)=α^i,jx+β^i,j are (α^i,j)1i,j4=[3.260.8313.573.4714.593.3310.1120.373.672.443.371.267.341.217.0510.06],(β^i,j)1i,j4=[00.222.490.044.380.773.157.151.470.881.360.015.870.726.119.06].

Fig 4
Example A: The identified transformation S^ of the underlying noisy system.

The performance of the reconstruction algorithm is evaluated by computing the relative error

δS(x)=100|S(x)S^(x)S(x)|,
(68)

for xX={0.01,0.02,,0.99}, which is shown in Fig. 5. The performance of the new reconstruction algorithm for different levels of noise was compared with that of a previous algorithm [38] that does not incorporate knowledge of the noise density. Table 1 shows for comparison the mean absolute percentage error (MAPE)

MAPE=100θδSi=1θδS|S(xi)S^(xi)S(xi)|,
(69)

for the two reconstruction approaches, where {xi}i=1θδS=X. The results, clearly demonstrate the advantages of the new algorithm and in particular its robustness even for in the presence of very significant noise levels.

Fig 5
Example A: Relative error between the original map S and the identified map S^ evaluated for 99 uniformly spaced points.
Table 1
Example A: Performance comparison between the new algorithm (A1) and the algorithm (A2) in [38] for different levels of noise.

6.2. Example B

This example demonstrates the proposed algorithm to reconstruct a nonlinear continuous map. Specifically, we consider the logistic map defined by

xn+1=4xn(1xn)+ξn(mod1)
(70)

where ξ is white noise following a Gaussian distribution function N(0,(5×103)2) truncated to the interval [ɛ,ɛ] with ɛ=0.02. The aim is to infer a piecewise linear semi-Markov defined over a uniform partition R with N=40 intervals, which approximates the logistic map S. The initial states X0i={x0,ji}j=1θ, i = 1,…,100, θ=5×103were generated by sampling a set of initial densities f0,ij(i)shown in Fig. 6 (see Appendix for more details). The corresponding final densities f1,ij(i) over the uniform partition R were estimated from X1i={xi,j}j=1θ, i = 1,…,100, the images under the noise perturbed transformation of the initial states X0i={x0,ji}j=1θ .

Fig 6
Example B: Examples of initial densities (gray lines) f0,ij(i) and the corresponding final densities after one iteration (black lines) f1,ij(i).

The approximate piecewise linear semi-Markov map, identified for ε = 0.02 using the algorithm in Section 5, is shown in Fig. 7.

Fig 7
Example B: Reconstructed piecewise linear semi-Markov map Ŝ over the uniform partition ={Ri}i=140.

The smoothed map obtained with the smoothing parameter 0.999 is shown in Fig. 8, and the relative error calculated on the uniformly spaces points is shown in Fig. 9.

Fig 8
Example B: Identified smooth map S¯ resulted from piecewise linear semi-Markov map in Fig. 7 with smoothing parameter 0.999.
Fig 9
Example B: Relative error between the original map S and the identified smooth map S¯ in Fig. 8 evaluated for 99 uniformly spaced points.

The root mean square error (RMSE) between the predicted density functions using original and identified maps calculated by

RMSE=1Ni=1N(viv^i)2,
(71)

where v^i is the coefficient of predicted density function, is given in Table 2. As can be seen in Fig. 9 and Table 2, the approximation error has been very low. With the increase of the interval number, the reconstructed map is more close to the original one, and the stabilized distribution converges to the invariant density of the system, as proven in [22].

Table 2
Example B: RMSE between the predicted density functions fn using the identified map and those generated by the original noisy system.

As in the previous example, the performance of the new reconstruction method was compared with that of a previous method [38] for different levels of noise. The results are summarized in Table 3. As it can be seen, the new algorithm performs significantly and consistently better. This is also illustrated in Fig. 10 in which the maps reconstructed using the two algorithms for epsilon = 0.15 and epsilon = 0.50 are shown side-by-side. It is worth noting however, that one of the advantages of the original algorithm in [38] is that it includes an additional step to optimize the partition.

Fig 10
Example B: Reconstructed maps for noise level epsilon = 0.15 using (a) the new algorithm and (b) the algorithm in [38] and the reconstructed maps for noise epsilon = 0.50 using (c) the new algorithm and (d) the algorithm ...
Table 3
Example B: Performance comparison between the new algorithm (A1) and the algorithm (A2) in [38] for different levels of noise.

7. Conclusions

This paper introduced a new method for reconstructing/approximating an unknown one-dimensional chaotic map that is perturbed by additive noise, from sequences of density functions. The emphasis here is on recovering the true transformation that generated the data rather than one of the many possible transformations that share the same invariant density functions.

By incorporating knowledge of the noise distribution, the new estimation method achieves dramatically better accuracy (i.e. over tenfold error reduction in some cases) for high levels of noise compared with a previous method that does not account for the noise density.

As highlighted in [38], it would be of interest to develop similar reconstruction approaches for higher-dimensional systems. The main challenge is to construct the transformation given the matrix-representation of the Frobenius–Perron operator [22].

Acknowledgments

X. N. gratefully acknowledges the support from the Department of Automatic Control and Systems Engineering at the University of Sheffield and China Scholarship Council. D. C. gratefully acknowledges the support from MRC (G0802627), BBSRC (BB/M025527/1) and the Human Frontier Science Program grant no. RGP0001/2012.

Appendix: Initial states for Example B in Section 6

The 100 sets of initial states used in the example are obtained by sampling the following density functions.

f0,1β1(x,β1)=710·x29(1x)β11B(30,β1)+310·xβ11(1x)29B(β1,30),β1=1,2,,30;f0,2β2(x,β2)=xβ21(1x)29B(β2,30),β2=1,2,,25;f0,3β3(x,β3)=x29(1x)β31B(30,β3),β3=1,2,,25;f0,4β4(x,β4)=12·x39(1x)β4+19B(40,β4)+12·x39(1x)β4+19B(40,β4),β4=1,2,,10;f0,5β5(x,β5)=12·xβ5+19(1x)39B(β5,40)+12·xβ5+19(1x)39B(β5,40),β5=1,2,,10;

where B( ·, ·) is beta function.

References

1. Lasota A, Mackey MC. 2nd ed. Springer-Verlag; New York: 1994. Chaos, fractals, and noise: stochastic aspects of dynamics.
2. Strogatz SH. Westview Press; 2014. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering.
3. Skinner JE. Low-dimensional chaos in biological systems. Nat Biotechnol. 1994;12:596–600. [PubMed]
4. Swishchuk A, Islam S. Taylor & Francis; 2013. Random dynamical systems in finance.
5. Lozowski AG, Lysetskiy M, Zurada JM. Signal Processing with temporal sequences in olfactory systems. IEEE Trans Neural Netw. 2004;15:1268–1275. [PubMed]
6. van Wyk MA, Ding J. Stochastic analysis of electrical circuits. In: Chen G, Ueta T, editors. Chaos in Circuits and Systems. World Scientific; 2002. pp. 215–236. edited by. edited by.
7. Huijberts H, Nijmeijer H, Willems R. System identification in communication with chaotic systems. IEEE Trans Circuits Syst I Fundam Theory Appl. 2000;47:800–808.
8. Lasota A, Rusek P. An application of ergodic theory to the determination of the efficiency of cogged drilling bits. Arch Górnictwa. 1974;19:281–295.
9. Simoyi RH, Wolf A, Swinney HL. One-dimensional dynamics in a multicomponent chemical reaction. Phys Rev Lett. 1982;49:245–248.
10. Guevara MR, Glass L. Phase locking, period doubling bifurcations and chaos in a mathematical model of a periodically driven oscillator: a theory for the entrainment of biological oscillators and the generation of cardiac dysrhythmias. J Math Biol. 1982;14:1–23. [PubMed]
11. Perry J, Smith RH, Woiwod IP, Morse DR. Springer Science & Business Media; 2012. Chaos in real data: the analysis of non-linear dynamics from short ecological time series.
12. Lai Y-C, Grebogi C, Kurths J. Modeling of deterministic chaotic systems. Phys Rev E. 1999;59:2907.
13. Skiadas CH, Skiadas C. CRC Press; 2008. Chaotic modelling and simulation: analysis of chaotic models, attractors and forms.
14. Maguire LP, Roche B, McGinnity TM, McDaid L. Predicting a chaotic time series using a fuzzy neural network. Inf Sci. 1998;112:125–136.
15. Han M, Xi J, Xu S, Yin F-L. Prediction of chaotic time series based on the recurrent predictor neural network. IEEE Trans Signal Process. 2004;52:3409–3416.
16. Príncipe J, Kuo J-M. Dynamic modelling of chaotic time series with neural networks. Adv Neural Inf Process Syst. 1995:311–318.
17. Lueptow R, Akonur A, Shinbrot T. PIV for granular flows. Exp Fluids. 2000;28:183–186.
18. Boyarsky A, Góra P. Laws of chaos: invariant measures and dynamical systems in one dimension. 1997
19. Ershov SV, Malinetskii GG. The solution of the inverse problem for the Perron–Frobenius equation. USSR Comput Math Math Phys. 1988;28:136–141.
20. Ulam SM. Interscience; New York: 1960. A collection of mathematical problems: interscience tracts in pure and applied mathematics.
21. Li T-Y. Finite approximation for the Frobenius–Perron operator: a solution to Ulam's conjecture. J Approx Theory. 1976;17:177–186.
22. Bollt EM. Controlling chaos and the inverse Frobenius–Perron problem: global stabilization of arbitrary invariant measures. Int J Bifurc Chaos. 2000;10:1033–1050.
23. Bollt EM, Santitissadeekorn N. SIAM; 2013. Applied and computational measurable dynamics.
24. Friedman N, Boyarsky A. Construction of ergodic transformations. Adv Math. 1982;45:213–254.
25. Góra P, Boyarsky A. A matrix solution to the inverse Perron–Frobenius problem. Proc Am Math Soc. 1993;118:409–414.
26. Diakonos FK, Schmelcher P. On the construction of one-dimensional iterative maps from the invariant density: the dynamical route to the beta distribution. Phys Lett A. 1996;211:199–203.
27. Pingel D, Schmelcher P, Diakonos FK. Theory and examples of the inverse Frobenius–Perron problem for complete chaotic maps. Chaos. 1999;9:357–366. [PubMed]
28. Koga S. The inverse problem of Flobenius–Perron equations in 1D difference systems―1D map idealization. Progr Theor Phys. 1991;86:991–1002.
29. Huang W. Constructing chaotic transformations with closed functional forms. Discret Dyn Nat Soc. 2006;2006
30. Huang W. Proceedings of the 2009 International Conference on Topics on Chaotic Systems: Selected Papers from CHAOS. 2009. On the complete chaotic maps that preserve prescribed absolutely continuous invariant densities; pp. 166–173.
31. Huang W. Constructing multi-branches complete chaotic maps that preserve specified invariant density. Discret Dyn Nat Soc. 2009:14. (2009)
32. Baranovsky A, Daems D. Design of one-dimensional chaotic maps with prescribed statistical properties. Int J Bifurc Chaos. 1995;5:1585–1598.
33. Diakonos FK, Pingel D, Schmelcher P. A stochastic approach to the construction of one-dimensional chaotic maps with prescribed statistical properties. Phys Lett A. 1999;264:162–170.
34. Rogers A, Shorten R, Heffernan DM. Synthesizing chaotic maps with prescribed invariant densities. Phys Lett A. 2004;330:435–441.
35. Rogers A, Shorten R, Heffernan DM. A novel matrix approach for controlling the invariant densities of chaotic maps. Chaos Solitons Fractals. 2008;35:161–175.
36. Berman A, Shorten R, Leith D. Positive matrices associated with synchronised communication networks. Linear Algebra Appl. 2004;393:47–54.
37. Rogers A, Shorten R, Heffernan DM, Naughton D. Synthesis of piecewise-linear chaotic maps: invariant densities, autocorrelations, and switching. Int J Bifurc Chaos. 2008;18:2169–2189.
38. Nie X, Coca D. Reconstruction of one-dimensional chaotic maps from sequences of probability density functions. Nonlinear Dyn. 2015;80:1373–1390.
39. Bollt E, Góra P, Ostruszka A, Zyczkowski K. Basis Markov partitions and transition matrices for stochastic systems. SIAM J Appl Dyn Syst. 2008;7:341–360.
40. Lasota A, Yorke JA. On the existence of invariant measures for piecewise monotonic transformations. Trans Am Math Soc. 1973;186:481–488.
41. Ding J, Li TY. Markov finite approximation of Frobenius–Perron operator. Nonlinear Anal Theory Methods Appl. 1991;17:759–772.
42. Ding J, Li TT. Projection solutions of Frobenius–Perron operator equations. Int J Math Math Sci. 1993;16:465–484.
43. Góra P, Boyarsky A. Approximating the invariant densities of transformations with infinitely many pieces on the interval. Proc Am Math Soc. 1989;105:922–928.
44. Hunt FY, Miller WM. On the approximation of invariant measures. J Stat Phys. 1992;66:535–548.
45. Kohda T, Murao K. Piecewise polynomial galerkin approximation to invariant densities of one‐dimensional difference equations. Electron Commun Japan (Part I: Commun) 1982;65:1–11.
46. Góra P, Boyarsky A. An algorithm to control chaotic behavior in one-dimensional maps. Comput Math Appl. 1996;31:13–22.