Lyapunov exponent and natural invariant density determination of chaotic maps: An iterative maximum entropy ansatz
LLyapunov exponent and natural invariant densitydetermination of chaotic maps: An iterativemaximum entropy ansatz
Parthapratim Biswas
Department of Physics and Astronomy, The University of Southern Mississippi,MS 39406, USAE-mail: [email protected] ‡ Hironori Shimoyama
Department of Physics and Astronomy, The University of Southern Mississippi,MS 39406, USA
Lawrence R. Mead
Department of Physics and Astronomy, The University of Southern Mississippi,MS 39406, USA
Abstract.
We apply the maximum entropy principle to construct the natural invariantdensity and Lyapunov exponent of one-dimensional chaotic maps. Using a novelfunction reconstruction technique that is based on the solution of Hausdorffmoment problem via maximizing Shannon entropy, we estimate the invariantdensity and the Lyapunov exponent of nonlinear maps in one-dimension from aknowledge of finite number of moments. The accuracy and the stability of thealgorithm are illustrated by comparing our results to a number of nonlinear mapsfor which the exact analytical results are available. Furthermore, we also considera very complex example for which no exact analytical result for invariant densityis available. A comparison of our results to those available in the literature is alsodiscussed.
1. Introduction
The classical moment problem (CMP) is an archetypal example of an inverse problemthat involves reconstruction of a non-negative density distribution from a knowledge of(usually finite) moments [1, 2, 3, 4, 5, 6]. The CMP is an important inverse problemthat has attracted researchers from many diverse fields of science and engineeringranging from geological prospecting, computer tomography, medical imaging totransport in complex inhomogeneous media [7]. Much of the early developments inthe fields such as continued fractions and orthogonal polynomials have been inspiredby this problem [2, 8]. The extent to which an unknown density function can bedetermined depends on the amount of information available in the form of moments ‡ Corresponding author a r X i v : . [ n li n . C D ] O c t provided that the underlying moment problem is solvable. For a finite number ofmoments, it is not possible to obtain the unique solution and one needs to supplementadditional information to construct a suitable solution. The maximum entropyprovides a suitable framework to reconstruct a least biased solution by simultaneouslymaximizing the entropy and satisfying the constraints defined by the moments [9].In this communication we address how the maximum entropy (ME) principle canbe applied to the Hausdorff moment problem [1] in order to estimate the Lyapunovexponent and the associated natural invariant density of a nonlinear dynamical system.In particular, we wish to apply our ME ansatz to a number of nonlinear iterative mapsin one-dimension for which the analytical results in the closed form are available.The problem was studied by Steeb et al [10] via entropy optimization for the tentand the logistic maps using the first few moments (up to 3). Recently, Ding andMead [11, 12] addressed the problem and applied their maximum entropy algorithmbased on power moments to compute the Lyapunov exponents for a number ofchaotic maps. The authors generated Lyapunov exponents using up to the first 12moments, and obtained an accuracy of the order of 1%. In this paper, we addressthe problem using a method based on an iterative construction of maximum entropysolution of the moment problem, and apply it to compute Lyapunov exponents andthe natural invariant densities for a number of one-dimensional chaotic maps. Unlikethe power moment problem that becomes ill-conditioned with the increasing numberof moments, the hallmark of our method is to construct a stable algorithm by resortto moments of Chebyshev polynomials. The resulting algorithm is found to be verystable and accurate, and is capable of generating Lyapunov exponents with an errorless than 1 part in 10 , which is significantly lower than any of the methods reportedearlier [10, 11]. Furthermore, the method can reproduce the natural invariant densityof the chaotic maps that shows point-wise convergence to the exact density functionwhenever available.The rest of the paper is organized as follows. In Section 2, we briefly introducethe Hausdorff moment problem and a discrete maximum entropy ansatz to constructthe least biased solution that satisfies the moment constraints. This is followed bySection 3, where we introduce the natural invariant density as an eigenfunction of thePerron-Frobenius operator associated with the dynamical system represented by theiterative maps [13]. In Section 4, we discuss how the moments of the invariant densityare computed numerically via time evolution of the dynamical variable, which are thenused to construct the Lyapunov exponents and the natural invariant densities of themaps. Finally, in Section 5, we discuss the results of our method and compare ourapproximated results to the exact results and to those available in the literature.
2. Maximum Entropy approach to the Hausdorff moment problem
The classical moment problem for a finite interval [a, b], also known as the Hausdorffmoment problem, can be stated loosely as follows. Consider a set of moments µ i = (cid:90) ba x i ρ ( x ) dx i = 0 , , , . . . , m, i ≤ m (1)of a function ρ ( x ) integrable over the interval with µ i < ∞ ∀ x ∈ [a,b] and ρ ( x )has bounded variation, the problem is to construct the non-negative function ρ ( x )from a knowledge of the moments. The necessary and sufficient conditions for asolution to exist were given by Hausdorff [1]. The moment problem and its variantshave been discussed extensively in the literature [2, 3, 14, 15] at length, and anauthoritative treatment of the problem with applications to many physical systemswas given by Mead and Papanicolaou [4]. For a finite number of moments, the problemis underdetermined and it is not possible to construct the unique solution from themoment sequence unless further assumptions about the function are made. Within themaximum entropy framework, one attempts to find a density ρ A ( x ) that maximizesthe information entropy functional, S [ ρ ] = − (cid:90) ba ρ A ( x ) ln[ ρ A ( x )] dx (2)subject to the moment constraints defined by Eq. (1). The resulting solution is anapproximate density function ρ A ( x ) and can be written as ρ A ( x ) = exp (cid:32) − m (cid:88) i =0 λ i x i (cid:33) . (3)The normalized density function ρ ( x ) is often referred to as probability density bymapping the interval to [0,1] without any loss of generality. For a normalized densitywith µ = 1, the Lagrange multiplier λ is connected to the others via e λ = (cid:90) exp (cid:32) − m (cid:88) i =1 λ i x i (cid:33) A reliable scheme to match the moments numerically for the entropy optimizationproblem (EOP) was discussed by one of us in Ref. [16]. The essential idea behindthe approach was to use a discretized form of entropy functional and the momentconstraints using an accurate quadrature with a view to reduce the original constraintoptimization problem in primal variables to an unconstrained convex optimizationprogram involving dual variables. This guarantees the existence of the uniquesolution [17], which is least biased and satisfies the moment constraints defined byEq. (1). Using a suitable quadrature, the discretized entropy and the momentconstraints can be expressed as respectively, S [ ρ ] = − (cid:90) ρ ( x ) ln[ ρ ( x )] dx ≈ − n (cid:88) j =1 ω j ρ j ln ρ j (4) µ i = (cid:90) x i ρ ( x ) dx ≈ n (cid:88) j =1 ( x j ) i ω j ρ j (5)where ω i ’s are a set of weights associated with the quadrature and ρ j is the value ofthe distribution at x = x j . If ω j and x j are the weight and abscissas of the Gaussian-Legendre quadrature, the Eq. (4) is exact for polynomials of order up to 2 n − n (cid:88) j =1 ω j = 1 , n (cid:88) j =1 ω j ρ j = 1 . (6)The task of our EOP can now be stated as, using ˜ ρ j = ω j ρ j and t ij = ( x j ) i , tooptimize the Lagrangian L (˜ ρ, ˜ λ ) = n (cid:88) j =1 ˜ ρ j ln (cid:18) ˜ ρ j ω j (cid:19) − m (cid:88) i =1 ˜ λ i n (cid:88) j =1 t ij ˜ ρ j − µ i (7)where 0 ≤ ˜ ρ ∈ R n and ˜ λ ∈ R m , respectively are the primal and dual variables ofthe EOP, and the discrete solution is given by functional variation with respect to theunknown density,˜ ρ j = ω j exp (cid:32) m (cid:88) i =1 t ij ˜ λ i − (cid:33) , j = 1 , , . . . n. (8)The Eqs. (4) to (8) can be combined together and a set of nonlinear equations can beconstructed to solve for the Lagrange multipliers ˜ λF i ( ˜ λ ) = n (cid:88) j =1 t ij ω j exp (cid:32) m (cid:88) k =1 t kj ˜ λ k − (cid:33) − µ i = 0 , i = 1 , , . . . , m. The set of nonlinear equations above can be reduced to an unconstrained convexoptimization problem involving the dual variables:min ˜ λ ∈ R m D (˜ λ ) ≡ n (cid:88) j =1 ˜ ρ j exp (cid:32) m (cid:88) i =1 t ij ˜ λ i − (cid:33) − m (cid:88) i =1 µ i ˜ λ i . (9)By iteratively obtaining an estimate of ˜ λ , D (˜ λ ) can be minimized, and theEOP solution ρ (˜ λ ∗ ) can be constructed from Eq. (8). In the equation above, t ij = x ij corresponds to power moments, but the algorithm can be implementedusing Chebyshev polynomials as well. The details of the implementation of the aboveapproach for shifted Chebyshev polynomials was discussed in Ref [16]. The maximumentropy solution in this case is still given by the Eq.(3) except that x i within theexponential term is now replaced by T ∗ i ( x ), where T ∗ i ( x ) is the shifted Chebyshevpolynomials. In the following, we apply the algorithm based on the shifted Chebyshevmoments to construct the invariant density of the maps.
3. Lyapunov exponent and the natural invariant density of chaotic maps
The Lyapunov exponent of an ergodic map can be expressed in terms of the naturalinvariant density of the map:Γ = (cid:90) ρ ( x ) ln | f (cid:48) ( x ) | dx (10)where ρ ( x ) is the invariant density and f (cid:48) ( x ) is the first derivative of the map f ( x ) withrespect to the dynamical variable x . The invariant density of a map can be definedas an eigenfunction of Perron-Frobenius operator associated with the map. Given aniterative map, x n +1 = f ( x n ), one can construct an ensemble of initial iterates { x } defined by a density function ρ ( x ) in some subspace of the phase space and considerthe time evolution of the density in the phase space instead of initial iterates x . Thecorresponding evolution operator L is known as Perron-Frobenius operator, which islinear in nature as each member of the ensemble in the subspace evolves independently.The invariant density can be written as, L ρ ( x ) = ρ ( x ) (11)where ρ ( x ) is a fixed point of the operator L in the function space. In general, theremay exist multiple fixed points but only one has a distinct physical meaning, whichis referred to as the natural invariant density. Following Beck and Schl¨ogl [13], thegeneral form of the operator in one-dimension can be written as, L ρ ( y ) = (cid:88) x ∈ f − ( y ) ρ ( x ) | f (cid:48) ( x ) | . (12)For an one-dimensional map, one can define the Lyapunov exponent as theexponential rate of divergence of two arbitrarily close initial points separated by δx n =0 = | x − x (cid:48) | in the limit n → ∞ , and the exponent can be expressed as theaverage of the time series of the iterative map,Γ = 1 N lim N →∞ N − (cid:88) n =0 ln | f (cid:48) ( x n ) | . (13)For ergodic maps the time average of the Lyapunov exponent can be replaced bythe ensemble average,Γ = (cid:90) dx ρ ( x ) ln | f (cid:48) ( x ) | (14)using the natural invariant density. Equation (14) suggests that the Lyapunovexponent can be obtained from a knowledge of the reconstructed natural invariantdensity from the moments. In the following we consider some nonlinear mapsto illustrate how the normalized invariant density and Lyapunov exponent can becalculated using our discrete entropy optimization procedure.
4. Reconstruction of invariant density as a maximum entropy problem
In the preceding sections, we have discussed how a probability density can beconstructed from a knowledge of the moments (of the density) by maximizingthe information entropy along with the moment constraints. Once the densityis reconstructed, the Lyapunov exponent can calculated from Eq. (14) using thereconstructed density. The calculation of the moments can proceed as follows. Weconsider a dynamical system represented by a nonlinear one-dimensional map, x n +1 = f ( x n ) (15)where n =0, 1, 2, . . . and x ∈ [0 , x n can be expressed as, < x i > = lim t →∞ t t (cid:88) n =0 ( x n ) i Since we are working with the shifted Chebyshev polynomials, the correspondingmoments are, µ i = lim t →∞ t t (cid:88) n =0 T i ∗ ( x n ) (16)where T ∗ i ( x ) are the shifted Chebyshev polynomials and are related to Chebyshevpolynomials via T ∗ i ( x ) = T i (2 x − x ∈ [0 ,
5. Results and Discussions
Let us first consider the case for which the invariant density function and the Lyapunovexponent can be calculated analytically. We begin with the map, f ( x ) = x − x ≤ x ≤ √ − − x x √ − ≤ x ≤ ρ ( x ) = 4 π (1 + x ) (18)and the Lyapunov exponent is given by ln 2, which can be obtained analytically fromEq. (14). The approximated Chebyshev moments for the map f ( x ) can be obtained Table 1.
Test Map f Moments Γ maxent
Percentage error20 0.691577 0.22640 0.692786 0.05560 0.692999 0.02180 0.693061 0.012100 0.693109 0.006Γ exact ln 2 ≈ Ref. [11] ≈ r defined via, f ( x ) = r x (1 − x ) . (19)We consider three representative values of r to illustrate our method in the chaoticand non-chaotic domain. In particular, we choose (a) r = , (b) r =4, and (c) r =3.79285. The analytical densities are known only for the first two cases, and are givenby respectively, ρ chaotic2 ( x ) = 1 π √ x − x ; r = 4 (20) ρ fixed2 ( x ) = δ (cid:18) x − (cid:19) ; r = 52 (21)For the remaining value of r = 3 . x , and constructing a histogram averaging over a number of configurations [13].For the purpose of comparison to our maximum entropy results, we use this numericaldensity here. Table 2.
Logistic Map: Case A f = x (1 − x ) Moments Γ maxent
Percentage error10 − − − − exact − ln 2 ≈ − . Table 3.
Logistic Map: Case B f = 4 x (1 − x ) Moments Γ maxent
Percentage error40 0.690703 0.3560 0.692101 0.1580 0.692850 0.04100 0.693319 0.02Γ exact ln 2 ≈ Ref. [11] r = and r = 4 respectively. The errors associated with Γ arealso listed in the respective tables. The invariant density for r = is a δ -function, andthe exact analytical value of the exponent is given by − ln 2. Since the invariant densityis a δ -function at x = , it is practically impossible to reproduce the density veryaccurately using a finite number of quadrature points. However, our maximum entropyalgorithm produces an impressive result by generating only two non-zero values in theinterval containing the point x = using Gaussian quadrature with 192 points. Theapproximate density for a set of 40 moments is shown in fig. 3. The two non-zero valuesof the density are given by 19.781 and 105.264 within the interval [0.593, 0.601]. It maybe noted that for a normalized density, one can estimate the maximum height of the δ -function to be of the order of (∆ x ) − ≈ .
0, where ∆ x is the interval containingthe point x = point [19]. Furthermore, we have found that the result is almostindependent of the number of moments (beyond the first 20), and the δ -function hasbeen observed to be correctly reproduced with few non-zero values using only as fewas first 10 moments. Table 3 clearly shows that the first 3 digits have been correctlyreproduced using only the first 10 moments. On increasing the number of moments,there is but very little improvement of the accuracy of Lyapunov exponent. For each ofthe moment sets, the density is found to be zero throughout the interval except at few(two for the set 40 and higher) points mentioned above. In absence of any structure inthe density, higher moments do not contribute much to the density reconstruction, andhence it’s more or less independent of the number of moments. Since the contributionto the Lyapunov exponent is coming only from the few (mostly two) non-zero values,and that these values fluctuate with varying moments, an oscillation of these valuescauses a mild oscillation in the Lyapunov exponent.We now consider the case r = 4. The exact density in this case is given by Eq. (20)that has singularities at the end points x = 0 and 1. It is therefore instructive to studythe divergence behavior of the reconstructed density near the end points. In fig.4 wehave plotted the approximate density obtained using the first 90 moments along withthe exact density. The reconstructed density matches excellently within the interval.The divergence behavior near x = 0 is also plotted in the inset. Although there issome deviation from the exact density, the approximate density matches very goodexcept at very small values of x . Such observation is also found to be true near x = 1. The results for the Lyapunov exponent are listed in table 3 for differentnumber of moments. It is remarkable that the exponent has been correctly producedup to 3 decimal points with 100 moments. While the error in this case is largercompared to the cases discussed before for the same number of moments, it is muchsmaller than the result reported earlier [11]. Our numerical investigation suggeststhat the integral converges slowly for Gaussian quadrature in this case owing to thepresence of a logarithmic singularity in the integrand. This requires one to use moreGaussian points to evaluate the integral correctly. However, since the density itselfhas singularities at the end points, attempts to construct the density very close to theend points introduce error in the reconstructed density that affects the integral value.The use of Gauss-Chebyshev quadrature would ameliorate the latter problem, but forthe purpose of generality (and in absence of prior knowledge of the density) we refrainourselves from using the Gauss-Chebyshev quadrature.Finally, we consider a case where analytical results are not available and thedensity consists of several sharp peaks having fine structure in the interval [0,1].As mentioned earlier, the case r = 3 . x and constructing an histogram of the distribution of the iterates in the longtime limit, which is then finally averaged over many configurations. While our MEPansatz produces most of the peaks in the exact density using the first 20 moments,the fine structure of the peaks is missing and so is the location of the peaks. Thereconstructed density can be improved systematically by increasing the number ofmoments, and for 150 moments the density matches very good with the exact density.In fig. 6 we have plotted both the reconstructed density for the first 150 momentsand the numerical density from the histogram method. The result suggests that forsufficient number of moments our algorithm is capable of reproducing density whichis highly irregular, non-differentiable and consists of several sharp peaks.
6. Conclusion
We apply an iterative maximum entropy optimization technique based on Chebyshevmoments to calculate the invariant density and the Lyapunov exponent for a numberof one-dimensional nonlinear maps. The method consists of evaluating approximatemoments of the invariant density from the time evolution of the dynamical variableof the iterative map, and to apply a novel function reconstruction technique viamaximum entropy optimization subject to moment constraints. The computedLyapunov exponents from the approximated natural invariant density are found tobe in excellent agreement with the exact analytical values. We demonstrate thatthe accuracy of the Lyapunov exponent can be systematically improved by increasingthe number of moments used in the (density) reconstruction process. An importantaspect of our method is that it is very stable and accurate, and that it does notrequire the use of extended precision arithmetic for solving the moment problem. Acomparison to the results obtained from power moments suggest that the algorithmbased on Chebyshev polynomials gives more accurate results than the power moments.This can be explained by taking into account the superior minimax property of theChebyshev polynomials and the form of the maximum entropy solution for Chebyshevmoments [20, 21]. Our method is particularly suitable for maps for which exactanalytical expression for invariant density are not available. Since the method candeal with a large number of moments, an accurate invariant density function canbe constructed by studying the convergence behavior with respect to the number ofmoments. The Lyapunov exponent can be obtained from a knowledge of the invariantdensity of the maps. Finally, the method can also be adapted to solve non-lineardifferential and integral equations as discussed in Refs.[22] and [23], which we willaddress in a future communication.PB acknowledges the partial support from the Aubrey Keith Lucas and Ella GinnLucas Endowment in the form of a fellowship under faculty excellence in researchprogram. He also thanks Professor Arun K. Bhattacharya of the University ofBurdwan (India) for several comments and discussions during the course of the work. [1] Hausdorff F, Math Z. , 220-248 (1923)[2] Shohat J A and Tamarkin J D, The Problem of Moments (American Mathematical Society,Providence, RI, 1963)[3] Akheizer, N. I,
The classical moment problem and some related questions in analysis , HafnerPublishing Co., New York 1965 [4] Mead, L.R. and Papanicolaou, N. J. Math. Phys.
25, 8 (1984)[5] A. Tagliani, J. Math. Phys. , 326 (1993)[6] R.V. Abramov, J. of Computational Physics, , 621 (2007)[7] Kirsch, A. An Introduction to Mathematical Theory of Inverse Problems , Springer-Verlag, NewYork[8] R. Haydock in
Solid State Physics: Advances in Research and Applications , Vol 35 editedF.Seitz, Academic Press (1980)[9] Jaynes E. T,
Probability Theory: The Logic of Science (Cambridge University Press, Cambridge,U. K., 2003)[10] Steeb W-H, Solms F and Stoop R 1994
J. Phys. A: Math. Gen. L399[11] Ding J and Mead L R 2002
J. of Math. Phys. , 658 (2007)[13] Beck, C and Schl¨ogl, F Thermodynamics of chaotic systems (Cambridge University Press,Cambridge, MA, 1993)[14] Wimp, J., Proc. Roy. Soc. Edinburgh 82 A, 273 (1989)[15] M. Junk, Math. Models Methods Appl. Sci , 1001 (2000)[16] Bandyopadhyay K, Bhattacharya A K, Biswas P and Drabold D A 2005 Phys. Rev. E least biased as far as entropy of the density isconcerned, and the Hausdorff conditions are satisfied. There is no guarantee that the maximumentropy solution would be close to the exact solution particularly for very few moments.However, the quality of the MEP solution drastically improves with increasing number ofmoments, and numerical experiments for cases where exact solutions are known confirm thatmaximum entropy principle indeed can produce the correct solution.[18] Feigenbaum, M.J. J. Stat. Phys.
19, 25 (1978); May, R.M.
Nature δ ( x − x ) can be approximated by a box function of width ∆ x around x and of height h , giving h ∆ x = 1 to satisfy the normalization condition. The heightof the δ -function obtained in our work is very close to this limiting value.[20] The presence of λ i x i in the maximum entropy solution for power moments makes it very difficultto exploit information from the high order moments for i >
20 in the interval [0, 1] even withextended precision arithmetic. Such a problem does not appear in our formulation of theproblem via Chebyshev moments owing to the bounded nature of Chebyshev polynomials,and consequently provides a way to incorporate systematically the information from the highermoments.[21] Mason, J.C and Handscomb, D.C
Chebyshev Polynomials (Chapman and Hall/CRC, New York,USA 2003)[22] J. Baker-Jarvis and D. Schultz,
Numerical Methods for Partial Differential Equations , , 133-142(1989), John Wiley & Sons Inc,[23] L.R. Mead, J. Math. Phys. , December 1986 Figure 1.
The reconstructed density from the first 20 Chebyshev moments forthe map f ( x ) along with the exact density. Although the general shape of thedensity appears correctly, some oscillations are present in the data in absence ofsufficient information. Figure 2.
The reconstructed density obtained from the first 100 moments for themap f ( x ) along with the exact density. The approximate density now effectivelymatches point-by-point with the exact one with the exception of few points nearthe edges of the interval. Figure 3.
The reconstructed density for the logistic map for r = obtainedfrom the first 40 Chebyshev moments. The density consists of only two non-zerovalues in the interval [0.593, 0.601]. The exact density is a δ -function centered at x = . Figure 4.
The reconstructed density for the map f ( x ) for r = 4 obtained fromthe first 90 moments along with the exact density. The density matches excellentlywithin the interval. The divergence behavior at the left edge near x = 0 is alsoshown in the inset. Figure 5.
The reconstructed density for the logistic map for r = 3 . Figure 6.
The reconstructed density for logistic map for r = 3 ..