{'fraction_non_alphanumeric': 0.03745886181266871, 'fraction_numerical': 0.012757460340938505, 'mean_word_length': 4.491167512690355, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ergodic Sensitivity Analysis of One-Dimensional Chaotic Maps

Sensitivity analysis in chaotic dynamical systems is a challenging task from a computational point of view. In this work, we present a numerical investigation of a novel approach, known as the space-split sensitivity or S3 algorithm. The S3 algorithm is an ergodic-averaging method to differentiate statistics in ergodic, chaotic systems, rigorously based on the theory of hyperbolic dynamics. We illustrate S3 on one-dimensional chaotic maps, revealing its computational advantage over naive finite difference computations of the same statistical response. In addition, we provide an intuitive explanation of the key components of the S3 algorithm, including the density gradient function.


Introduction
Sensitivity analysis is a discipline that studies the response of outputs of a certain model to changes in input parameters. It involves computing the derivatives of output quantities of interest with respect to specified parameters. This mathematical tool is essential in many engineering and scientific applications, as it enables optimal design of structures [11,28] and fluid-thermal systems [2], analysis of heterogeneous flows [22], supply chain management [26], estimate errors and uncertainties in measurement, modeling and numerical computations [3,24].
For example, consider an optimal design problem in structural mechanics involving an elastic truss under loading. In such problems, the ultimate goal might be to approximate derivatives of constrained functions (e.g. resulting displacements) with respect to design parameters (e.g. bar cross-sections). Such derivatives can be approximated using analytical methods, as well as finite differences [13]. A classic example from fluid mechanics is a turbulent flow past a rigid object, in which the sensitivity of the drag (resistance) forces with respect to the Reynolds number and other flow parameters [4,19], is of interest. Particularly, aerospace engineers use the computed sensitivity in the design of airfoils [10]. Both mechanical phenomena described above are governed by strongly nonlinear dynamical systems, however the latter features an extra difficulty, namely the chaotic behavior.
Computing such sensitivities in chaotic dynamical systems is a challenging task. The primary issue is the so-called butterfly effect, which is a large sensitivity of the system to initial conditions. This concept is associated with the classical study of Edward Lorenz on climate prediction [16]. Quantitatively, it means that any two points initially separated by an infinitesimal distance diverge at an exponential rate. This implies the prediction of far-future states in chaotic phenomena is hardly possible. We observe this phenomenon in daily weather forecasts, as the predictions of several weeks forward tend to be highly inaccurate. However, we are sometimes interested in predicting the response of long-time averaged behavior, to perturbations [4,19,7].
2 Parameterized one-dimensional chaotic maps and their statistical dependence on parameter In this section, we introduce two families of perturbed one-dimensional chaotic systems that can generally be expressed as x k+1 = ϕ(x k ; s), where x init is a given initial condition, while s denotes a scalar parameter. Let J be a scalar observable. Our quantity of interest is the infinite-time average or ergodic average of J, In particular, we focus on the relationship between J and the parameter s for ϕ. In addition, we review the concepts of Lyapunov exponents and ergodicity through numerical illustrations on the two maps.

Perturbations of the sawtooth map and their Lyapunov exponents
We consider as our first example, perturbations of the sawtooth map, also known as the dyadic transformation, defined in the following way: It is a periodic map that maps [0, 1) to itself. Figure 1 illustrates the sawtooth map for different values of the parameter s. A natural question that arises is whether the chosen is map actually chaotic. Roughly speaking, a chaotic map shows high sensitivity to initial conditions. For example, consider s = 0, and two phase points x and x + δx. Under one iteration of the map, these two points are now separated by a distance of 2δx. Thus, in the limit δx → 0, a trajectory that is infinitesimally separated from x at n = 0 moves away from the trajectory of x at an exponential rate of log 2 ≈ 0.693. This exponential growth of perturbations to the state is the signature of chaotic systems and is measured by the rate of asymptotic growth, known as the Lyapunov exponent (LE) and denoted by λ. More rigorously, the Lyapunov exponent is defined by which clearly indicates that the infinite-time averaged rate of growth converges to a constant. We say that a map is chaotic when its LE is positive. Formula 4 requires computing the derivative of the map at points along a trajectory. Note that the value of the Lyapunov exponent does not depend on initial condition x init , nor on the step k. Figure 2 shows that λ > 0 for all s ∈ (− 1 2π , 1 π ] meaning that the sawtooth map is chaotic in this regime. This can be easily justified by the observation that ∂ϕ ∂x ≥ 1 for all x ∈ [0, 1], when s is in this regime. In case of the classical Bernoulli shift, i.e. when s = 0, repetitive the sawtooth map always appears to converge to a fixed point, after some iterations, when simulated numerically. This is because all machinerepresentable numbers with a fixed number of digits after the binary point, are dyadic-rational numbers, which converge to the fixed point 0, under this map, because the sawtooth map at s = 0 is simply a leftshift operation on binary digits. More details about this problem and possible remedies can be found in Appendix A.1. Note also that if s = 0 and x init is rational, the forward orbit of x init would either converge to a fixed point or be periodic, containing a finite number of distinct values within the interval [0, 1). For example, if x init = 0.1, then all future states belong to a four-element set, {0.2, 0.4, 0.6, 0.8}, and x k = x k+4 for all k > 0. This is an example of an unstable periodic orbit; in this paper, we are interested in chaotic orbits, which are aperiodic and unstable to perturbations.

A family of Cusp maps and their Lyapunov exponents
Another example of a chaotic map is the cusp map ϕ : [0, 1] → [0, h] defined as follows, The above function produces a spade-shaped graph, as shown in Figure 3. The cusp map is a two-parameter map with s = {h, γ}, where h is the height, while γ is a parameter that determines the sharpness of the tip. We use the definition Eq. 4 to compute the LE of the cusp map at different values of h and γ. From the positivity of the Lyapunov exponent shown in Figure 4, we see that the cusp map is always chaotic if γ ∈ [0, 1] and h ≥ 0.6.  Historically, the cusp map has been used as a one-dimensional representation of the three-dimensional Lorenz'63 system [16], a set of ordinary differential equations used as a model for atmospheric convection. Specifically, the iterates of the cusp map are local maxima of the third coordinate of the Lorenz'63 system [17].

Ergodic probability distributions
The long-time average of the objective function, J , was calculated using 100 million iterates of the map, with the initial condition chosen uniformly, at random between (0,1), in the following way: where x i+1 = ϕ(x i ). We choose a sufficiently large N to ensure that the right hand side converges to a fixed value, within numerical precision. Figures 5 and 6 illustrate examples of the mean statistics (i.e. long-time averages) and their dependence on the map parameters for the sawtooth and cusp map, respectively.  In the computation of the long-time averages of the objective function, we used the concept of ergodicity. This property guarantees that long-time averages do not depend on the initial conditions. That is, the time average of the objective function (right hand side of Eq. 6) converges, as N → ∞, to a value independent of the initial condition x 0 , for almost every x 0 chosen uniformly between (0, 1). This limit equals the expected value of the same objective function over an ensemble of initial conditions distributed according to an ergodic, invariant probability distribution ρ. This probability distribution ρ is invariant under ϕ, in the sense that for any open interval A ⊂ (0, 1), ρ(A) = ρ(ϕ −1 (A)). In addition, ρ is defined by the fact that expectations with respect to ρ are the same as infinite-time averages starting from a point uniformly distributed in the unit interval. Such a probability distribution ρ, is known as the SRB distribution [27], and only sometimes coincides with the uniform distribution (it does e.g. for the sawtooth map at s = 0). The above description can be mathematically rephrased as follows, for almost every x 0 uniformly distributed in (0, 1), Thus, in ergodic systems, there exist two alternative ways of computing the long-time average, either through the averaging of the time-series or ensemble averaging. The latter requires prior computation of the probability distribution, which will be explained and illustrated in the next sections. Using these preliminary concepts, we will review the space-split sensitivity method to compute the derivative of J with respect to the map parameter.

Split Space Sensitivity analysis of one-dimensional maps
In [21], Ruelle rigorously derived a formula for the derivative of the quantity of interest, J , with respect to the map parameter s. This expression is an ensemble average (or expectation) with respect to ρ, which can be simplified for one-dimensional maps ϕ : where reflects the parameter perturbation of the map, while U refers to the unit interval [0, 1). A direct evaluation of Eq. 8 is computationally cumbersome for the following reason. Notice that the integrand of the right hand side of Eq. 8 involves a derivative of the composite function that can be expanded using the chain rule to the form As discussed earlier, for a large k, the product of the derivatives exponentially grows with k. However, Ruelle's series converges due to cancellations of these large quantities, upon taking an ensemble average. This problem makes the direct evaluation of Ruelle's formulation computationally impractical since a large number of trajectories are needed for these cancellations. More precisely, since for a large k, at almost every x, we need to increase the number of trajectories by a factor of O(e 2λk ) in order to reduce the mean-squared error in a linear fashion. For example, consider the sawtooth map with s ∈ [− 1 2π , 1 π ]. In this case, (∂ϕ/∂x) ∈ [1,4]. One can easily verify that even for moderate values of k, an overflow error is encountered. Another challenge is that the evaluation of the SRB distribution requires expensive computation of map probability densities [5]. In a recent study [8], Ruelle's formula has been reformulated to a different ensemble average, known as the S3 formula. There, the latter formula has been derived for maps of arbitrary dimension, and is based on splitting the total sensitivity into that due to stable and unstable perturbations. Note the notion of splitting the perturbation space is irrelevant for 1D maps, and the one-dimensional perturbation is, by defintion of chaos, unstable. Therefore we will skip some aspects of the original derivation, and note that our derivation represents only the unstable component of sensitivity in [8] specialized to 1D.
The S3 formula, corresponding to equations 8-9, can be expressed as follows: where and, For one-dimensional maps, the derivation is simple, as it requires integrating 8 by parts and the fact that the integral of df /dx at the boundary of U vanishes; see Appendix A.2 for the full derivation. We observe that both J and df /dx have their analytical forms. However, the function g(x), which will be referred to as density gradient, does not have a closed-form expression, since the SRB distribution ρ, is unknown. The density gradient g represents the variation in phase space, of the logarithm of ρ(x), In the next section, we focus on further interpreting g(x), its computation and verification on the 1D maps introduced in Section 2.

Computation of density gradient
In this section, we focus on the density gradient function, denoted by g(x). First, we present a computable, iterative scheme for g(x). Moreover, we provide an intuitive explanation for g and visualize it on the maps introduced in Section 2.
Based on the S3 formula (Eq. 12), we can conclude the following recursive relation holds. The full derivation of Eq. 16 is included in Appendix A.3. This recursive procedure can be used to approximate g(x) along a trajectory in the asymptotic sense, which means that we need a sufficiently large number of iterations to obtain an accurate approximation of g(x) [8]. In practice, we generate a sufficiently long trajectory, compute first and second derivatives of the map evaluated along the trajectory, and apply Eq. 16. We arbitrarily set g(x init ) = 0, to start the recursive procedure, and obtain g(ϕ(x init )). The recursion is continued by setting x = ϕ(x init ), and so on. For a sufficiently large K, the true value of g(ϕ K (x init )) is approached, for almost every initial condition x init .

Interpretation of the density gradient iterative formula
To intuitively understand the density gradient formula (Eq. 16), we isolate the effects of each term in Eq. 16. In order to do this, we consider a small interval around an iterate x k and examine two cases: 1. the map is a straight line on this interval and 2. the map has a constant curvature on this interval. These two cases are graphically shown on the left (numbered as 1) and right hand sides (numbered as 2) of Figure 7. The x-axis represents an interval around an iterate x k and the y-axis, an interval around x k+1 = ϕ(x k ). The density ρ, around each interval, is shown adjacent to the axes, as a colormap. The colors reflect the distribution of ρ on a logarithmic scale. Figure 7: A graphical representation of two different scenarios in one-dimensional maps, to intuitively understand the derivation of the quantity g. The bold lines illustrate the map, while shaded bars adjacent to each axis represents the corresponding density distribution on that axis. The region around H corresponds to a high value of density, while the region around L to low values. The slope of the line is indicated as t.
1. Consider a small region of (x − , x + ) where the map ϕ(x) has zero second derivative, i.e., the first derivative dϕ/dx is constant. As shown in Figure 7(1), let us assume that the density on the left side of the region, ρ(x − ), is higher than the density on the right side, ρ(x + ). Due probability mass conservation, the mapped density can be calculated using the following equation, Since we consider case dϕ/dx > 0, we drop the absolute value. On this interval where the map is a straight line, this statement says that the density around ϕ(x) is a constant multiple of the density around x. On the logarithmic scale, the density around ϕ(x) is shifted by a constant, when compared to the density around x since, This relationship is graphically depicted in Figure 7 (1), where the regions marked H and L, corresponding to higher and lower densities, are shifted to the left. Notice that Eq. 18 implies that the difference, on the logarithmic scale, between the higher and lower densities on the y-axis equals the difference between the higher and lower densities on the x-axis. However, the small interval is stretched by a factor of dϕ/dx under one iteration of the map. Thus, the derivative of the logarithm of density decreases by a factor of dϕ/dx. Mathematically, we can see this by differentiating both sides of the equation with respect to x (and using that dϕ/dx is constant), From the definition of g, this reduces to, which is confirmed by our formula, Eq. 16, by setting d 2 ϕ/dx 2 = 0.
2. To isolate the effect of curvature of the map on g, we consider a curved map and a constant density region. Thus, by definition, g(x) = 0 in the interval considered. We now describe that g • ϕ(x) becomes non-zero on this interval due to the curvature of the map. Note that Eq. 18 still applies, since it is a restatement of probability mass conservation. This means that ρ is reduced by a factor equal to the slope of the map at every point. This is graphically depicted in Figure 7, in which we have assumed that dϕ/dx is an increasing function that crosses the value 1 at the point indicated using dashed lines.
To the left of this point, the density is therefore increased (shown as H) around ϕ(x) and to the right, the density is decreased (shown as L), when compared to its uniform value around x. Note also the larger the first derivative of the map, the lower the density on the y-axis. Again, by taking the derivative of Eq. 18 with respect to x, and using the definition of g, we obtain As mentioned in Case 1, the first derivative (dϕ/dx)(x), gives the factor by which a length dx around x is stretched (or compressed) by ϕ. The second derivative gives us the change of this stretching (or compression) as a function of x. Thus, the effect of a non-zero second derivative is felt by the derivative of the density, and can again be derived from measure preservation or probability mass conservation.

Numerical examples of density gradients
In the second part of this section, we show numerical results of the density gradient procedure in two examples, the sawtooth and cusp maps, which were introduced in Eq. 3 and Eq. 5 respectively. Figure  8 shows the stationary probability densities of the sawtooth map at different values of s. We observe that all curves appear differentiable , however their derivatives are large, near the interval boundaries, when s is close to −1/(2π) or 1/π. In Figure 9 we show the distribution of the (averaged) density gradient function, g(x), computed using Eq. 16, at different values of s, and compare it against its finite difference approximation: (log(ρ(x + )) − log(ρ(x − )))/ (2 ). Note that the expected value of the density gradient is always zero since We also repeat a similar experiment for the cusp map, whose results are presented in Figures 10-11. We observe a behavior similar to the sawtooth map. Similar to the sawtooth density, the densities computed for the cusp map appear to be differentiable over a range of the parameter γ. However, as γ gets close to 1, the density ρ(x) acquires large derivatives at the boundaries of the interval. The boundedness of dρ/dx is needed for the computation of g(x) to be well-conditioned. Figure 10: The plot shows the empirically estimated stationary probability distributions achieved by the cusp map (Eq. 5), at h = 1 and the indicated value of γ. All curves were generated in the same fashion as for the sawtooth case. Figure 11: The plot compares g(x) (solid lines) against the derivative of the empirically estimated stationary probability distributions (dots) achieved by the cusp map (Eq. 5), at h = 1 and indicated value of γ. All curves were generated in the same fashion as for the sawtooth case.

Spaces-split sensitivity as a sum of time-correlations
The evaluation of Eq. (12) is the main focus of this paper. In practice, expectations with respect to ρ, or ensemble averages, are computed by time-averaging on a single typical trajectory. As mentioned earlier, a time average converges to the ensemble average of the function, as the length of the trajectory approaches infinity. Thus, Eq. 12 can be written as follows, replacing the ensemble averages with ergodic averages where x i = ϕ i (x init ) is the point at time i, along a trajectory starting at a typical point x init . Using the definition of g, and taking a long trajectory,

Numerical examples of sensitivities computed using S3
To numerically verify Eq. 24, we consider a set of objective functions, each of which is an indicator function denoted by δ c , and defined such that its value is a constant 1 in a small interval around c and zero everywhere else on the unit interval. With this particular choice, Eq. 24 gives us the gradient of the probability density, Thus, by varying the constant c in the interval [0, 1), and using the density gradient computed using Eq. 15, one can compute dρ/ds over the unit interval by using Eq. 24. This can be compared with the finite difference approximation of dρ/ds generated using slightly perturbed values of s and approximating the density empirically. This particular choice of J(x) exhibits yet another advantage of the S3 method over Ruelle's formula (Eq. 8). The former is also applicable to objective functions that have non-differentiable points, since unlike a direct evaluation of Ruelle's formula, the derivative of J is not used. Figure 12 shows numerical results for the cusp map, in which the density gradient is computed using the space-split formula (Eq. 24) and compared with the central difference derivative. We observe that only a few terms of the series are required to produce accurate sensitivities.   Figure 13 clearly indicates that the consecutive terms of the series in Eq. 23 exponentially decay in norm. We repeat a similar experiment for the sawtooth map (see Figures 14 and 15). In this case, we only need three terms of Eq. 24 to obtain a result that is indistinguishable from its finite difference approximation. The consecutive terms of Eq. 24 also decay exponentially in norm.  Note that each term of Eq. 12 is in the form of a lag-k time correlation between ∇ ρ f and J. We use the term "lag-k" as ∇ ρ f and the objective function are evaluated at two different states that are k steps apart. In mixing systems, the lag-k time correlations converges to zero as k → ∞. Moreover, for a family of dynamical systems known as Axiom A, the rate of decay of time correlations is proven to be exponential [27]. In the case of one-dimensional maps, Axiom A systems are the ones in which the derivative of the map is different than 1 everywhere. All the map examples we consider in this paper satisfy this requirement. This guarantees that only a small number of time correlation terms are needed to secure high accuracy of the sensitivity approximation.

Computational performance of S3
Finally, we compare the space-split sensitivity and classical finite difference method in terms of computational efficiency. We observe in Figures 16-17, generated for the sawtooth and cusp map, respectively, that the S3 method clearly outperforms its competitor, as it requires a few orders of magnitude lesser samples to generate a result with a similar relative error. This is a very promising observation in the context of analysing higherdimensional systems, since the large cost of generating very long trajectories can make such computations infeasible. Note in the case of both the S3 and finite difference methods, the error is upper-bounded as follows [8], where N denotes the number of samples, while C is some positive number. This means we observe a convergence rate of a typical Monte Carlo simulation in both methods. However, the factor C is substantially larger in case of finite differencing. Moreover decreasing the step size (indicated as δs) in the finite difference calculation, worsens the accuracy, due the dominance of statistical noise. Figure 16: Relative error of the space-split and finite difference methods as a function of the trajectory length. We compute the parametric derivative of density of the sawtooth map at s = 0.1 on the left boundary (x = 0). For the S3 computation (curve marked with ( )), we consider only first three terms of Eq. 12, which corresponds to the line marked with (•) in Figure 14. For the finite difference approximation, we calculate densities at s = 0.105, 0.095 (curve marked with (•)) and s = 0.1005, 0.0995 (curve marked with ( )). We also computed the S3 approximation using 125,829,120,000 samples and 9 terms of Eq. 12, which serves as a reference value. The dashed lines are proportional to the inverse of the square root of the number of samples. For the S3 computation (curve marked with ( )), we consider first ten terms of Eq. 12. For the finite difference approximation, we calculate densities at γ = 0.505, 0.495 (curve marked with ( )) and γ = 0.5005, 0.4995 (curve marked with ( )). We also computed the S3 approximation using 125,829,120,000 samples and 17 terms of Eq. 12, which serves as a reference value. The dashed lines are proportional to the inverse of the square root of the number of samples.

Conclusions and future work
We demonstrate a new method to compute the statistical linear response of chaotic systems, to changes in input parameters. This method, known as space-split sensitivity or S3, is used to compute the derivatives with respect to parameters of the long-time average of an objective function. In the S3 method, a quantity called density gradient, defined as the derivative of the log density with respect to the state, is obtained using a computationally efficient ergodic averaging scheme. An intuitive explanation of this iterative ergodic averaging scheme, based on probability mass conservation, is discussed in this paper. The density gradient plays a key role in the computation of linear response. Specifically, the sum of time correlations between the density gradient and the objective function partially determines the derivative of the mean statistic of the objective function with respect to the parameter. The computational efficiency of the S3 formula when compared to finite difference, which requires several orders of magnitude more samples, stems precisely from this new formula to efficiently estimate the density gradient.
In this work, we restrict ourselves to expanding maps in 1D, which are simple examples of chaotic systems. These examples nevertheless give rich insight into chaotic linear response, and specifically into the behavior of the density gradient. Our study shows that in same cases the derivative of the density gradient might be very large, which corresponds to heavy tailedness of the density gradient distribution. This phenomenon, as well as its implication for analysis of higher-dimensional maps, is the main topic of our future work.

A.3 Derivation of the iterative procedure for g in 1D maps
The purpose of this section is to derive the iterative procedure to calculate the density gradient g. We use the same notational convention as in Appendix A.2. Let us consider a function h that is integrable in U = (0, 1) and vanishes at U L = 0 and U R = 1. Using the definition g = (1/ρ)(∂ρ/∂x), and integrating by parts, we obtain The key property used in this derivation is the density preservation of ϕ. We say that the map ϕ is densitypreserving with respect to the density ρ, if for any scalar observable f , U f ρ dx = U f • ϕ k ρ dx holds for any integer k. This implies the left-hand side of Eq. 1 can be expressed as