Response Operators for Markov Processes in a Finite State Space: Radius of Convergence and Link to the Response Theory for Axiom A Systems

Using straightforward linear algebra we derive response operators describing the impact of small perturbations to finite state Markov processes. The results can be used for studying empirically constructed—e.g. from observations or through coarse graining of model simulations—finite state approximation of statistical mechanical systems. Recent results concerning the convergence of the statistical properties of finite state Markov approximation of the full asymptotic dynamics on the SRB measure in the limit of finer and finer partitions of the phase space are suggestive of some degree of robustness of the obtained results in the case of Axiom A system. Our findings give closed formulas for the linear and nonlinear response theory at all orders of perturbation and provide matrix expressions that can be directly implemented in any coding language, plus providing bounds on the radius of convergence of the perturbative theory. In particular, we relate the convergence of the response theory to the rate of mixing of the unperturbed system. One can use the formulas derived for finite state Markov processes to recover previous findings obtained on the response of continuous time Axiom A dynamical systems to perturbations, by considering the generator of time evolution for the measure and for the observables. A very basic, low-tech, and computationally cheap analysis of the response of the Lorenz ’63 model to perturbations provides rather encouraging results regarding the possibility of using the approximate representation given by finite state Markov processes to compute the system’s response.


A Brief Summary of Response Theory
The development of methods for computing the response of a complex system to small perturbations affecting its dynamics is the subject of very active investigation in many fields of science and of technology. Statistical mechanics provides tools for approaching such a problem through so-called response theories, which allow for evaluating the change in the properties of a system through suitably defined operators that factor in the statistical properties of the unperturbed system and the specific nature of the perturbation one wants to study.
One can see a response theory as a virtual experimental setting where one has at hand a given system, various measurement instruments, and a knob controlling the value of a parameter, and knows how to relate the position of the knob with the reading of the instruments. In other terms, response theories provide the basis for understanding the outcome of experiments, and, not by chance, physical sciences have been at the forefront of the theoretical investigation in this direction. The monumental contribution by [1] provided the basis and the explicit formulas needed for studying the impact of very general perturbations to statistical mechanical systems at equilibrium, as described by the canonical ensemble. The Kubo formulas are extremely useful for studying a large class of problems in e.g. transport, optics, and acoustics. A cornerstone of Kubo's theory is the fluctuation-dissipation relation, which enables connecting-within linear approximation-the free fluctuations of a system to its response to perturbations. This property is closely related to the celebrated diffusion law for the brownian motion and has been recently extend to a fully nonlinear case [2]. Despite its obvious relevance, Kubo's approach has been criticized for several reasons: • it is not physically consistent in treating the transition from equilibrium to nonequilibrium dynamics, because it studies the impact on equilibrium systems of perturbations that drive them near (but out of) equilibrium, but does not clarify how a new stationary state is reached and maintained; additionally, it is not suited for studying the response to perturbations of non-equilibrium systems; • it lacks mathematical rigour, as it is not clear which are the systems for which the response formulas apply, and why it should apply at all.
In [3][4][5] it was clarified that it is possible to establish a rigorous response theory for Axiom A [6] continuous or discrete time dynamical systems. One obtains that the invariant SRB measure is smooth with respect to the parameter that controls the strength of the perturbation changing the dynamics of the system fromẋ = F(x) toẋ = F(x) + X(x), in the case of continuous time evolution, and from x k+1 = F(x k ) to x k+1 = F(x k ) + X(x k ), in the discrete case. We continue our discussion taking into consideration the continuous case. We can introduce the unperturbed evolution operator S t 0 = exp(tF·), which moves forward in time any function of phase space O(x) by an interval t according to the unperturbed dynamics, so that O(x(t)) = S t 0 O(x(0)), and its perturbed counterpart S t = exp(t (F+ X)·), which instead describes the evolution in the perturbed system.
We define ρ 0 (dx) and ρ (dx) the invariant measures of the unperturbed and perturbed states, respectively. In particular, one obtains that the expectation value of sufficiently smooth observables O(x) in the perturbed state can be expressed in the form: where [Q] = ν (dx)Q(x) and [Q] 0 = ν 0 (dx)Q(x), while the various terms of the perturbative expansion can be written as: where (•) = X · ∇(•). In particular, the linear term can be written as: All terms δ[O] j can be written as an expectation value on the unperturbed measure of a new observable expressed as a functional of the background vector field F, of the perturbative vector field X, and of the observable O. The somewhat surprising conclusion we draw is that the invariant measure of the system, despite being supported on a strange geometrical set, is differentiable with respect to . Among the many merits of the Ruelle response theory, one can mention that a) it clarifies the mathematical framework needed for developing a response theory, whose main ingredient, roughly speaking, is the robustness deriving from having a uniformly hyperbolic dynamics on the attractor supporting an SRB measure; and b) it works seamlessly, in principle, in equilibrium and non equilibrium statistical mechanical systems, reducing to Kubo's formulas when considering the first scenario, if one assumes that statistical mechanical systems are Axiom A. Non-trivial implications of the nonequilibrium/equilibrium dichotomy regarding the validity of the fluctuation-dissipation relations are discussed in [2,5,7], while the a physical interpretation of the first and second order terms occurring in Ruelle's response formalism is provided in [8].
Of course, at this stage one needs to bridge the gap between mathematical formalism and physical meaningfulness, One manages to bring Ruellle's formalism into the realm of applicability by adopting the chaotic hypothesis [9,10], which basically says that a high-dimensional chaotic physical system can be treated at all practical purposes as if it were Axiom A if we focus on macroscopic observables. The chaotic hypothesis is the generalisation of the ergodic hypothesis, and provides a firm background for translating the mathematical properties of Axiom A systems into physically meaningful statements. Clearly, the chaotic hypothesis applies far from regimes of metastability and far from critical transitions, where entirely different phenomena appear. The chaotic hypothesis might also be practically problematic in the case one treats multiscale systems featuring many near-zero Lyapunov exponents; see discussion in [11].
Taking the point of view of the chaotic hypothesis, one has that, after transients have died out, nonequilibrium systems reach a nonequilibrium steady state (NESS) where the phase space is on the average contracting (with the rate of contraction corresponding, broadly speaking, to the entropy production of the system [12]), so that one can associate to the hyperbolic strange attractor supporting the invariant measure a Hausdorff dimension that is lower that the dimensionality of the phase space and, in general, not integer [6,13].
The last piece of the puzzle one needs to lay in order to sort out the above-mentioned criticisms to Kubo's theory relies on the physical interpretation of how a perturbed equilibrium system reaches a steady state. A convincing point of view on this relies on emphasizing the role of thermostats, which are large physical systems interacting with the system of interest in such a way to extract the excess of heat generated as result of the energy input due to the perturbation. Thermostats are also responsible for making it possible the set-up of stationarity in the case of forced and dissipative non equilibrium systems. An extensive treatment of the role of thermostats in equilibrium and nonequilibrium systems in the context of the chaotic hypothesis is given in [14]. We will not elaborate further on this aspect here.

Transfer Operator Approach
One can point out that the formulas above describe the impact of and expressed in terms of expectation values of a generic observable O, whereas one might like to derive directly results for the impacts of the perturbations on the invariant measure.
In [3][4][5] one constructs the response of the system to perturbations by following the changes in the individual trajectories and summing over the possible initial configurations distributed according to the unperturbed invariant measure. A different point of view on response theory focuses on studying the properties of the unperturbed and perturbed transfer operators and of their generators (see [15] for an introduction on these mathematical objects), through the construction of an appropriate framework of suitable (Banach) functional spaces where their actions are well defined, able to carefully treat the fundamental differences between the (smooth) unstable and (singular) stable manifolds of the Axiom A systems [16][17][18][19].
The evolution of the measure driven by the systemẋ = F(x) up to time t ≥ 0 starting from an initial condition at time t = 0 is described by the Perron-Frobenius transfer L t (see, e.g., [15]), so that ρ(x, t) = L t ρ(x, 0). We have that the family of {L t } t≥0 forms a oneparameter semigroup, such that L t+s = L t L s and L 0 = 1. The Perron-Frobenius operator L t is the adjoint of the evolution operator S t = L t , so that S t O, ρ = O, L t ρ , where f, g is the action (computation of the expectation value) of the linear functional g (the probability measure) on the test function f (the observable). We have that L t ν 0 = ν 0 ∀t ≥ 0, meaning that the invariant measure is an eigenvector corresponding to unitary eigenvalue of the Perron-Frobenius operator.
Assuming strong continuity and boundedness of the semigroup given by {L t } t≥0 , we can introduce the unperturbed Liouvillian operator L, which is the generator of the unperturbed Perron-Frobenius operator L t = exp(t L), and write the Liouville evolution equation for ρ(x, t) as follows [20]: One immediately obtains that Lν 0 = 0. In general, the spectrum of L is complex and in a strip of finite width including and below the imaginary axis consists only of isolated eigenvalues of finite multiplicity corresponding to the Ruelle-Pollicott resonances, while below such a strip one finds the essential spectrum, which is responsible for the continuum of the power spectra of integrable observables. Furthermore, the presence of a unique SRB measure comes from the presence of a simple vanishing eigenvalue, while mixing properties result from the absence of any other eigenvalue along the imaginary axis. The relevance of these properties for constructing a response theory are discussed in great detail in [18,19]. In [21] it is argued, using mathematical considerations and examples of geophysical relevance, that the presence of Ruelle-Pollicott resonances having real part close to zero may lead to the presence of rough parameter dependence, as the smoothness of the response if lost. Additionally, in [22], it is shown, along similar lines, that the crisis of a very high-dimensional chaotic attractor near a critical transition-namely, of a climate model in the vicinity of the tipping point responsible for the transition between warm and snowball climate [23][24][25][26]-can be detected and anticipated by looking at spectrum of the transfer operator. We then have that the presence of the perturbation to the dynamics changes the Liouville equation as follows: so that we can introduce the perturbed Perron-Frobenius operator L t = exp(t L ), which pushes forward in time the measure according to the perturbed dynamics: ρ(x, t) = L t ρ(x, 0). Clearly, S t O, ρ = O, L t ρ . Additionally, we have that L t ν = ρ ∀t ≥ 0 and L ν = 0. While this approach is in some sense mathematically more problematic, because it is based on studying a partial differential equation instead of a finite dimensional dynamical system, it seems to provide a more comprehensive set of tools for studying the response of a system and relating it to its unperturbed fluctuations, see, e.g., [16], where Ruelle's formulas are obtained along these lines. See also a comprehensive review given in [18], where the applicability of the response theory beyond the case of Axiom A systems is discussed in detail.. One needs to emphasise that the transfer operator approach is more natural in all the cases when our interest focuses on studying the properties of the response of an ensemble of trajectories (initialised according to the unperturbed invariant measure) rather than on individual orbits of a system.
Note that in some applications there is not an obvious separation between the two approaches. Let's take the problem of constructing climate projections through the use of (extremely complex) numerical climate models, which is one of the core activities summarized in the IPCC reports [27]. Indeed, modelling centers are actively pursuing the preparation of multiple runs starting from an ensemble of initial conditions for a given scenario of forcing in order to estimate more accurately the uncertainties in the projections. Nonetheless, we will not experience an ensemble of realizations of the climatic evolution, but just one.

Computing the Response
The analysis of high-dimensional complex system in terms of direct numerical simulation and of time series analysis suffers from the (almost) ubiquitous curse of dimensionality, which makes it hard to represent correctly the details of the dynamics because computational complexity explodes with the number of degrees of freedom. The construction of efficient and accurate algorithms for studying the response of a complex system to perturbations faces serious difficulties. Let's focus now on the linear case. Some previous studies have emphasised the need for treating separately the contributions to the response coming from short and long-time delayed contributions in Eq. 3, and have underlined the need for reducing the complexity of the invariant measure by adding in the background state some stochastic forcing, able to smooth out the singularity of the SRB measure [28,29].
A promising way to deal with the actual computation of the scalar product in Eq. 3 is to use as time-dependent basis the covariant Lyapunov vectors [30,31], which automatically separate the contributions to the response coming from the unstable, neutral, and stable directions. This clarifies that the convergence of the formula given in Eq. 3 comes from the two distinct facts that (a) perturbations along the stable directions naturally decay, and (b) perturbations along the unstable directions grow in size, but are dominated by the loss of correlation due to mixing.
Recently, algorithms based upon adjoint methods have shown a good degree of accuracy and seem promising, even if scaling them up to high-dimensional systems has not been attempted yet [32,33]. A different approach to the problem has been proposed in [7,[34][35][36], where, instead of trying to computing ab initio and directly the response given in Eq. 3, the authors construct it a posteriori, probing the system with some test forcings and using the formal properties of the theory to be able to predict the response for new patterns of forcings. One can say that by studying the differential response to similar yet differently modulated perturbations, it is possible to derive the overall response properties of the system.

This Paper
Any numerical representation of a continuum system builds upon the need of discretizing the phase space and, in the case of time-continuous system, of time.
In this case, we partition the phase space of the system in say N states φ 1 , . . . , φ N . In many cases, the states are constructed by discretizing the phase space in a grid of boxes, which provide a (Galerkin) basis of orthogonal functions. We then construct an initial ensemble as defined by the occupancy where 1(A) is the characteristic function in the set A, and we want to approximate the evolution of such occupancies change with time, considering discrete time steps t, so that, to a good approximation the occupancy at time k t is Moreover, in such a discrete representation, we have that the value of an observable O in the state φ i is given by its average Let's emphasize that when analyzing virtually any sort of complex system, almost invariably one proposes a natural spatial and temporal cut-off, so that one on not in fact interested in really being able to compute the response of any possible observable defined at any possible spatial and temporal resolution, whereas meso-or macroscopic properties are relevant. Going again to the useful example of climate science, it is commonly regarded as a good and useful question to learn about the change in the surface temperature in response to climate forcing on a spatial scale corresponding to say a continent or a fraction thereof, and on a temporal scale of say one year. Nobody would find useful nor intelligent to study the surface temperature response over extremely small temporal and spatial scales. Empirically, using long numerical integrations and defining the set of finite states φ i , i = 1, . . . , N , we can construct the stochastic matrix M i, j describing the probability of performing a transition from state φ i to state φ j in a period of time t. The same operation can in principle be performed using experimental and observational data. A fundamental issue at the core of such procedure is whether for some dynamical systems in the limit of finer and finer partitions covering the phase space (actually, the attractor of the system) with N → ∞ one reconstructs the actual invariant measure of the original system. See in [37] a comprehensive discussion of such an issue, the so-called Ulam conjecture, and in [38] some extremely promising applications of finite state Markov processes for studying severely reduced representations of complex systems.
Following the idea that the performing the discretization of the phase space amounts to adding a stochastic perturbation of the original dynamical systems, with intensity going to zero with the scale of the actual partitions, and exploiting the fact that the SRB measure can be constructed as zero-noise limit (with measure that is absolutely continuous with respect to Lebesgue) of the physical measure, in [39,40] it has been proposed that the Ulam conjecture applies in the case of Axiom A systems, which are endowed with an SRB measure. The convergence in the case of Anosov diffeormorphism has indeed been proved provided one adds some noise of asymptotically vanishing intensity (through stronger than the noise induced by the partition itself) to the underlying dynamics [41]. Somehow this is not so surprising because by adding noise one introduces a cutoff below which partitions do indeed work. At any practical level, these results suggest that in the case of Axiom A system constructing finite state Markov processes using Ulam partitions can do a pretty good job in simulating the true dynamics, if one consider reasonably well-behaved, smooth observables as test functions. Nonetheless, one has to note that different choices for the partitions can lead to very different rates of convergence [37]. See also the discussion and the numerical examples presented in [42].
Apart from the Ulam method, one can follow a mathematically more elegant but practically much harder way to construct finer and finer partitions. As well known, Axiom A systems possess Markov partitions, i.e. well-defined, metric independent, finite resolution representations of the phase space that refine themselves with the dynamics [6,14]. Such Markov partitions can be used to construct in the limit the actual SRB measure of the system, and, additionally, following [43], they provide a natural way to build finite Markov chains whose properties converge in the limit to those of the Perron-Frobenius operator of the system.
Having a response formulas in the finite case has direct relevance for finite Markov chains and for interpreting the results of reduced models. Another good reason to construct a response theory in a finite state space has to do with the fact that the response operators for Axiom A systems introduced by Ruelle can be written as expectation value of certain observables on the unperturbed SRB measure. Therefore, given what said above, one can hope to have convergence of the finite state reconstructed response operators to the corresponding true response operator in the limit of infinitely fine partitions of the dynamics. Actually, providing explicit formulas for the response operator for a finite state partition of a system the response operator and taking the limit for (suitably defined) finer and finer partitions could be interpreted as a rigorous way for constructing the actual response on the asymptotic SRB measure. One needs to note-see discussion in Sects. 2.1 and 3-that special attention has to be paid when studying the convergence of such operators.
In what follows, we present the derivation of the response formulas at all orders of perturbations (as well as the full nonlinear versions) for finite state spaces of arbitrary size N . All expressions are given in terms of the transitions matrix of the unperturbed system, to its corrections to the perturbation, and of the parameter controlling the strength of the perturbation. The interest we see in the calculations we present below is mostly three-fold: • our results are obtained using basic linear algebra operations in finite dimensional spaces, which can used to interpret more complex operators acting on infinite dimensional spaces. It is also possible to use the finite dimensional expressions to derive, e.g., the the actual response operators for continuous time Axiom A dynamical systems; • we are able to derive an explicit expression for the a lower bound to for the radius of convergence of the perturbative theory, and relate it with the mixing properties of the unperturbed system. We also find a (very tentative) expression for such a lower bound in the case of continuous time case Axiom A dynamical systems; • our formulas can be translated into one-line commands in now widely available software tools like R, Octave, or MATLAB . This might greatly facilitate the actual implementation of response operators. In particular, we can say that our results provide a direct translation of the response theory into a readily implementable algorithms.
The paper is organised as follows. In Sect. 2, we introduce some notation and provide basic properties of ergodic finite state Markov chains, which can be taken as mathematical model on its own or as finite precision representation of ergodic (in this case, Axiom A) systems. We also show how it is possible to find an exact expression for the impact of a perturbation on the invariant measure of the Markov process and we study the radius of convergence of the perturbative expansion. In Sect. 3 we rephrase our results in terms of observables, by constructing straightforward adjoint operators in finite dimensions. In Sect. 4 we show how our findings agree with the response theory for continuous time systems when we suitably translate the matrix operations into operators. In Sect. 5 we present a simple yet instructive investigation of the response of the Lorenz '63 system [44] to perturbations using Ulamlike partitions and the formalism developed here. In Sect. 6 we recapitulate and discuss our results.

Response Operators for Finite-State Markov Processes
Let's consider an ergodic Markov process with a finite number of states defined by the Ncomponent vector u. We consider the infinite Markov chain generated as u 0 , Mu 0 ,. . . M n u 0 , . . . where u 0 is the initial ensemble of states, and M i, j ∈ R N ×N is the stochastic transition matrix determining the probability of reaching the state i at step n if at step n − 1 we are in the state j. The process is taken to be stationary, so that M does not change with n. We The invariant measure is obtained by solving the eigenvalue problem and selecting the unique solution with eigenvalue λ = 1. The corresponding (column) eigenvector u 1 1 is the invariant measure of the system. We also remind that where {λ j , u j } j = 1, . . . , N are the pairs of eigenvalues and eigenvectors of M, where λ 1 = 1, |λ j | < 1 if j > 1, and z can be expressed as z = N j=1 α j u j . Our goal is to find a formula for expressing the change in the invariant measure resulting from perturbing the transition matrix M → M + m.
We note that in order to preserve the Markov property of the system, m obeys the following constraint: N j=1 m i, j = 0, so that N j=1 M i, j + m i, j = 0. Moreover, an additional constraint on m comes from the fact that all elements of M + m have to be positive. We define and clearly, − ≤ 0 ≤ + , and the perturbed matrix is a stochastic matrix ∀ ∈ [ − , + ]. In order to have some room for studying the impacts of perturbations, we require that + − − > 0.
with unitary eigenvalue. We define v 1 as the invariant measure of the perturbed system. Our goal is to express it as a function of M, m, and u. This amounts to constructing a response theory. We first present the results of the explicit calculation, and then discuss issues of well-posedness of the problem and convergence of the procedure in Sect. 2.1. Let's express v 1 = u 1 + ∞ n=1 n w n , so that we obtain: Note that the first eigenvalue is not changed by the perturbation M → M+ m, because also M + m is a stochastic matrix. Using the definition of u 1 we obtain a system of concatenated equations (1 − M)w n = mw n−1 , ∀n ∈ N, n > 1.
We obtain Given the recursive structure, we immediately derive the general formula: where n = n 1 . Concluding, we have that: (1 − M) −1 m u 1 (18) which provides the formula we have been looking for. We note that the term responsible for the nth order of perturbation to the measure can be expressed as Using the matrix identity we can also formally express the previous result as: Using again the matrix identity (1 − M) −1 = ∞ k=0 M k , the previous expression can be rewritten as: or

Well-Posedness and Convergence
In the previous equations, we have used somewhat carelessly the expression (1 − M) −1 . Unfortunately, the matrix 1 − M is not invertible, because all of its columns sum up to zero, or, alternatively, because we know that 1 is an eigenvalue of M. Nonetheless, the expression makes sense if we apply it to a vector belonging to span{u 2 , . . . , u n }. We now want to prove that: eigenvectors (u 1 , u 2 , . . . , u N ), and corresponding eigenvalues (λ 1 = 1, λ 2 , . . . , λ N ), 1 > |λ 2 | ≥ . . . |λ N |, and m is a matrix matrix R N → R n such that n i=1 m i, j = 0, then mz ∈ span{u 2 , . . . , u n } ∀z ∈ R n . Proof Let's consider the vector y = mz. Its ith component can be written as Remark One needs note that finite numerical precision might cause troubles, so that one should be careful in eliminating any component along u 1 at each before applying ∞ j=1 M j . Note that we must use ∞ j=1 M j expression for (1 − M) −1 in any code, because otherwise any software would give us automatically a NaN as error message.
Remark We wish to underline another method for avoiding the NaN problem mentioned above. Following [45], we introduce the fundamental matrix of the Markov chain as Z = (1 − M + M ∞ ) −1 , where M ∞ is the limit matrix whose columns are all equal to u. One can show that Z exists as the operation of inverse is well defined given the spectral properties of M − M ∞ [39]. One can show that M ∞ mz = 0 ∀z ∈ R N . Therfore, in all the previous Eqs.

16-22 we can substitute
Let's consider the problem of convergence of the expression in Eq. 18. We want to make sure that the L 1 norm of ∞ n=1 n w n does not diverge, and use this to find a bound for the value of . A simple way to approach this problem is to study the ratio of the L 1 norm of two consecutive terms in the previous series. Using Eqs. 15-16, we have: where we use the submultiplicative property of the norm and we introduce a modified definition of the L 1 norm taking into account that the vector mv ∈ span{u 2 , . . . , u N }∀v ∈ R N : Using expression 25, we have that the perturbative expression converges if The previous expression provides an explicit bound for our calculations. We note that max is finite because of the restriction imposed in the definition of the norm || • || * . Such a bound ensures also the invertibility of (1 − 1 ) −1 . From the previous result, we find the following bound for the first order correction to the invariant measure: so that ||m|| 1 /(1 − ||M|| * ) can be though as a bound to the first order sensitivity of the measure to perturbations.
Using expression 24, we can derive a more generous bound for : while ||m|| 1 ||(1−M) −1 || * 1 provides an additional (stricter) bound to the first order sensitivity. Note that in all the previous expressions we can substitute ||(1 − M) −1 || * 1 with ||Z|| 1 . At this point, we wish to refer to previous results (see, e.g., [46]) providing bounds for the L 1 norm of the difference between the perturbed and unperturbed invariant measure: where τ M (1) is the so-called ergodicity coefficient [47] defined as: ||M(e i − e j )|| 1 with e i indicating the unit vector having 1 at the ith entry and zero elsewhere. We remind that τ 1 (M) is larger than any subdominant eigenvalue of M, and 1/(1 − τ M (1)) can be taken as a definition of conditioning number of M [48]. Clearly if τ M (1) is close to 1, the bound given in Eq. 28 diverges. Note that 1/(1 − τ M (1)) is the bound to non-perturbative sensitivity mirroring the bound to the perturbative, linearized sensitivity given previously as 1/(1 − ||M|| * 1 ). See also additional results presented in [49]. The sensitivity of the unperturbed measure to perturbations given in Eq. 28 can also be cast in terms ρ M , the smallest possible value for the constant controlling the rate of convergence of iterates Me i , M 2 e i , . . ., M n e i to u 1 , so that ∀n ∈ N + , ∀i ∈ 1, . . . N we have that ||M n e i − u 1 || 1 ≤ Cρ n M , C ≥ 1 [46,48]. The sensitivity diverges as ρ M approaches 1, i.e. when the unperturbed matrix has slow properties of convergence.
While the quantities ||M|| * 1 , τ M (1), and ρ M are indeed different, they all point to the fact that if the mixing rate of the unperturbed matrix M is slow-so that such quantities are close to 1 (so that ||(1 − M) −1 || * 1 and ||Z|| 1 are very large)-then the sensitivity of the measure to perturbations is high. See in [21] a discussion of the link between slow mixing of a system and the presence of rough parameter dependence in its response to perturbations, with some examples of applications in a geophysical context.
Bringing together the results presented in Eqs. 9, 10 and in Eq. 27, we conclude that Eqs.

Response Theory for Observables
Let's now look at the problem in terms of impact of the perturbation m on the expectation value of observables. Observables live in the dual space of the densities, and, given our convention, they are row vectors. They are approximated as having a constant value within each cell of the chosen partition of the phase space. The expectation value of the observable π with respect to a measure w can be written as π, w , where •, • denotes the scalar product. By definition, we have that π, Aw = A π, w , where A indicates the transpose (and adjoint, because we are studying real functions) of A. Let's look at the change in the expectation value of the observable π as a result of M → M + m. We can write: = π, u 1 + ∞ n=1 n n π, u 1 (30) where [π] 0 = π, u 1 is the expectation value of π in the unperturbed system, [π] = π, v 1 is the expectation value of π in the perturbed system, δ[π] n is the nth order perturbation, which can be expressed as Moreover, n is the nth order adjoint response operator, acting on the observables, which can be written as: We can also wrote Eq. 31 as: where the last two expressions provide the nonperturbative formulas.
Remark Equations 22 and 31 provide at all orders the response formulas for the discrete Markov process studied here. If we are constructing empirically the discrete phase space, we expect that different choices of the partitions, corresponding to different approximate representations of the full dynamics, will deliver different results in terms of response. Hence, our results can be model dependent, which is reasonable, as we are starting from a subjective choice on the way we approximate the phase space. In fact, one can empirically test the robustness of the obtained results against a set of given criteria by comparing whether the perturbations to a certain set of relevant observables weakly depend on the specific partition used. We present a very preliminary (and encouraging) numerical study performed on the Lorenz '63 model [44] later in Sect. 5. Moreover, as discussed in Sect. 1.4, if we construct finer and finer partitions of for studying the response of systems whose unperturbed dynamics features an SRB invariant measure (most notably in the case of Axiom A systems), and indeed if we follow the self-refining Markov partitions of the dynamics, our results should converge to the exact response theory built upon the true SRB measure.
One needs to note that Eq. 27 gives an estimate of the largest possible value of for a given partition, but we are are not sure whether the minimum over all the finer and finer partitions of * max is positive-this corresponds to imposing the uniform-in N -bound on the norm of ||(1 − M) −1 || * 1 or ||Z|| 1 . In [39] it is shown that L 1 convergence of the finite state measure constructed using the Ulam method to the actual SRB measure is realized when ||Z|| 1 grows asymptotically not faster than log N , where N is the number of states. The requirement we seem to have here for applying response theory here is unavoidably stricter because computing the response entails considering the expectation value of not necessarily well behaved observables, constructed through nontrivial operations of differentiation of the actual observables of which we want to study the sensitivity to perturbations, see Eq. 2 and [3][4][5]. This essential difficulty is exactly what motivates the point of view discussed in [18,50], where a delicate analysis of the relationship between tangent space of the unperturbed dynamics, the perturbation flow, and of the observable allow to set up a robust framework for the response theory.
Similarly, in our case, making the response theory work at practical level means having/choosing m and u in such a way that ||(1 − M) −1 || * 1 or ||Z|| 1 grossly overestimates in terms of norm the effect of applying (1 − M) −1 or equivalently Z in, e.g., Eq. 22. Additionally, a suitable choice of the observable π can help avoiding potential singularities in Eq. 36. In other terms, response theory can work much more easily once we get rid of or cure pathological cases.

Towards Continuous Time Dynamical Systems
We want to rephrase the previous results in the context of continuous time dynamical systems and derive some formulas previously presented in the literature concerning Axiom A systems. We coonsider a time continuous dynamical system of the formẋ = F(x) and study its response to the perturbation F(x) → F(x)+ X(x). Correspondingly, as a result of the perturbation, the original invariant measure ν 0 (dx) is changed into ν (dx). The Liouville equation describing the evolution of a given initial density of states ρ(x) for the unperturbed system can be written as considering two instants of time separated by a small time interval dt, we have: We understand that M is in this context the unperturbed Perron-Frobenius operator L dt pushing forward the measure ρ from t to t + dt. When looking at the perturbed flow we have: where m = dtX X = −∇ · (X(x)•) (40) In this case, starting from Eq. 23, and considering that no normalization is applied to the perturbation operator, it is possible to propose a definition of * max for the continuous time dynamics taking inspiration from Eq. 27: * such that the perturbative expansion converges if ≤ * max , where || • || B describes the norm of the operator in the appropriate Banach space B it belongs to, while || • || * B is such that the computation of the norm excludes the SRB measure. Note that * max is finite if both ||X || B || and ||F −1 || * B are finite. This expression is admittedly tentative. As mentioned before, the problem of selecting appropriate functional spaces for constructing the response theory for Axiom A systems along the lines of studying the perturbations to the transfer operator requires a careful construction of suitable Banach spaces and of the related metrics [16,18,19] and is beyond the scope of this paper. 2

Linear Response
We now want to derive the Ruelle response formulas for computing the linear correction to the invariant measure resulting from the perturbation. We write where n indicates the order of perturbation. Let's first go back to the first order term in Eq. 15: Each term of the form M k pushes forward up to time t k = k × dt what is positioned to its right. Summing over k in, in fact, amounts to looking forward in time. If we insert the definition of m given above, we get the integrating factor dt, so that we obtain the following expression: where the evolution takes place according to the unperturbed system, and we have used the invariance of ν(dx) with respect to such an evolution law. By going into the dual space of the observables, we have that the change in the value of an observable O(x) from time t to time t + dt in the unperturbed system can be written as: where the operator M = 1 + dtF = 1 + dtF(x) · ∇(•). Along the same lines, one derives that the perturbation operator m acting on the observable can be written as m = dtX = dtX(x) · ∇(•). Furthermore, we introduce the following expansion for the expectation value of O(x): where [O] is the expectation value in the perturbed system, [O] 0 is the unperturbed expectation value, and the corrections are included in the summation. Applying this expression to the first order term in Eq. 31-33: we get: which is exactly the original version of Ruelle's linear response formula given in Eq. 3. One needs to note that what in Ruelle's formulation is causality (time integration in the response starts from 0), in the context of the Markov matrices formalism followed here comes from the algebraic expansion of (1 − M) −1 . The issues of convergence mentioned in the original paper by Ruelle can be translated in the rate of mixing of the system as determined by the properties of M discussed in Sect. 2.1.

Higher Order Terms
We can repeat the same construction to derive the higher order perturbation terms in the case of the continuous time dynamical systems. Inserting in Eqs. 15, 16 the expression 38 for M and expression 40 for m, we obtain for the second order the following expression for the perturbation to the invariant density: while the expression for the nth order correction reads like Considering the adjoint problem and computing the higher order corrections to the expectation value of the observable O, we derive the general response formula proposed by Ruelle as reported in Eq. 2.

A Very Basic Numerical Experiment
In order to make a (very) preliminary assessment of the potential of some of the ideas presented in this paper, we have focused on investigating some properties of the celebrated Lorenz '63 system [44]:ẋ where we have chosen the standard value for the parameters σ = 10, ρ = 28, and β = 8/3. We remark that such a system is not an Axiom A, but instead a singular hyperbolic system [52], which possesses a chaotic attractor and an invariant SRB measure [53]. In a previous publication [34], we have performed an analysis of the linear and nonlinear response of the Lorenz '63 to perturbations, extending a previous investigation by Reick [54], which makes us confident that response theory can be safely applied at all practical purposes also in this case. We consider the special case of time-indepedent perturbations to the dynamics resulting from substituting ρ → ρ + in Eq. 53, so that the perturbation flow can be written as We have then identified a 3-dimensional box B containing the attractor, defined as B = {(x, y, z) ∈ R 3 |x ∈ [−20, 20], y ∈ [− 30,30], z ∈ [−0, 50]}, and subdivided it, á la Ulam, in smaller boxes of identical size using a regularly spaced cartesian grid. We have considered partitions obtained using small boxes with linear dimension given by dx = 2 × j, dy = 3 × j, and dz = 2.5 × j, along the three directions, with j = 1, 2, 4, see Fig. 1. This amounts to partitioning B into 8000/j 3 smaller boxes. Note that our construction delivers a much lower resolution with respect to what used in, e.g., [55].
We run the model with standard values of the parameters choosing as initial condition [1 1 1] (in fact, given the global attractivity and ergodicity of the Lorenz attractor, any initial condition can be chosen), and, after discarding a transient of 1000 time units, which brings us safely into the asymptotic regime, we run the model for 50,000 time units with a simple Runge-Kutta 4th order adaptive scheme and obtain the output with time step of 0.001 time units. This takes less than 10 minutes in a today's commercial laptop with stan- Fig. 1 Attractor of the Lorenz '63 system with indication of the cartesian grids used for constructing the partitions of its phase space. See text The first row refers to the integration of the Lorenz model given in Eq, 53. The other rows refer to the empirical discrete Markov chain constructed using boxes of different sizes. N j B refers to the number of states The linear response of the observables defined in Eq. 32 has been obtained using Eq. 48. The derivative with respect to is estimated using finite differences with = 0.1. See text dard specifics using MATLAB . We present results at such a low level of sophistication in order to clarify that the appracch proposed here is rather robust and of relatively simple implementation.
As the box-counting dimension or capacity of the attractor of the model given in Eq. 53 is d 0 ∼ 2.05, we expect that the number of boxes B j k , k = 1, . . . , N j B needed to cover the attractor decreases N j B ∝ 1/j d 0 . We obtain a slightly lower exponent ∼1.9, which is perfectly acceptable as we are far from the asymptotic regime where the scaling given by d 0 is realized.
For each value of j, the boxes B j k define the discrete states φ j k , k = 1, . . . , N j B . By counting the number of times the trajectory is included in each state φ j k and normalizing we derive experimentally the asymptotic normalized occupanciesū j k . Instead, by tracking the transitions between the various discrete states, we construct the estimate of the stochastic transition matrix M j p,q describing the probability that the state φ j q makes a transition to the state φ j p in one time step. By finding the eigenvector corresponding to the unique unitary eigenvalue of M j p,q , we find the invariant measure, which agrees up to very high precision with the empirical occupancy rateū k computed from the trajectory. As a first step, we evaluate the expectation values of four meaningful observables given by x 2 , y 2 , z 2 , and z, as obtained from the time integration of the Lorenz model and from its discrete representation in terms of Markov chain. Table 1 shows that the agreement is rather good even when extremely coarse resolution is used.
We then show how to compute the response of the system to the perturbation due to the introduction of the vector field X(x). We keep in mind that when continuous time dynamics is considered, there is a very simple linear relation between the perturbation flow and the corresponding perturbation to the Perron-Frobenius operator, see Eqs. 38-40. Therefore, we repeat the the steps described above for the −perturbed flow (we choose = 0.1 in order to be on the safe side in terms of convergence), compute the new stochastic transition matrices M One needs to note that because of the non-infinite integration time considered, of the non-infinitesimal perturbation applied, and of the somewhat arbitrary choice of the boxes, it can happen that the original and perturbed flow may be characterized by a different number of discrete states. We have observed such a difference only in the case j = 1, involving one single extra state for the perturbed flow, with normalized relative occupancy (≤10 −6 ). This problem can be easily sorted out by imposing a cutoff and removing from the the discrete description all states with very low.
As discussed above, one needs to test accurately the well-posedness and convergence of the expansion in order to be sure to obtain meaningful results. This is not our goal at this stage given such a preliminary numerical test of our results. Therefore, we limit ourselves to the less ambitious yet interesting goal of computing the linear response defined in Eq. 32 for the observables indicated above, using Eq. 48. The results are reported in Table 1 and seem very encouraging. We have that the estimates of the response are very stable with respect to changes in the resolution of the boxes, and agree to a high degree of precision with the results one obtains by empirically evaluating the sensitivity of the observables with respect to the introduction of the perturbation flow using two integrations, as well, in the case of the z observable, with what reported in [34]. We note that the results are virtually unchanged if one uses instead of the high resolution time series with time step of 0.001 time units sparser observations corresponding to, e.g. a time step of 0.01 time units. Obviously, using a time resolution lower by a factor of s with respect to what considered here, one derives by tracking the transitions a stochastic transition matrix corresponding to the sth power of the one obtained at higher resolution. This does not affect the results as long as the sampling is much higher than the characteristic time scale of the system, which can be approximated in ∼1/λ 1 ∼ 1.1 time units, where λ 1 is the positive Lyapunov exponent of the system. On much longer time scales, instead, the stochastic matrix is quasi-degenerate, with all columns almost equal to the invariant measure

Conclusions
Taking the point of view of finite state Markov systems, we have been able to construct a perturbation theory for studying the impact of small perturbations to the background dynamics. While previous approaches focus on the constructing a theory able to account for the effect of adding small perturbations to the baseline flow, we focus on computing the change in the invariant measure and for the change in the expectation values of general observables (one problem being the adjoint of the other) occurring when the Markov transition matrix M → M + m.
The perturbation term m has to be such that all the columns of the new stochastic matrix sum up to 1 and all entries are positive. All of our findings are obtained with rather simple linear algebra manipulations and using basic properties of the stochastic matrices. We can express the response as a perturbation series or, after suitable resummation, using compact exact formulas. We are also able to assess the convergence properties of the response theory by defining a value * max such that if | | ≤ * max the perturbative expansion converges. We have that the stronger is the mixing of the unperturbed system, the larger is the value of max . These findings match well with previous results providing upper bounds to the sensitivity of stochastic matrices to perturbations.
Our results provide a direct algorithmic method for studying the response to perturbations for finite state Markov processes and have the advantage of allowing for an immediate and practical change of point of view between response theory seen in terms of changes of the invariant measure or in terms of changes in the expectation values of observables, by simply computing the transpose of the resulting finite dimensional linear operators. Our findings give closed formulas for the linear and nonlinear response theory at all orders of perturbations through explicit matrix expressions that can be directly implemented in any coding language.
We can use our formulas to study the response to perturbations of finite state Markov processes constructed in order to have a simplified and treatable picture of a complex system. Given two different state spaces constructed using different finite partitions covering the attractor of the system, we cannot expect to obtain the same results for the change in the expectation value of a given observables. The results might indeed be model dependent, but this is the obvious price one has to pay because of the subjective choice of the reduced state space. An assessment of the robustness of the obtained results is key to applying our methods in the context of reduced models. Nonetheless, the extremely unsophisticated numerical study reported here on the Lorenz'63 model is quite encouraging at this regard, even if test should be made on much higher dimensional models.
If the underlying dynamics is Axiom A (or Axiom A equivalent, as in the cases where the chaotic hypothesis applies), one can impose conditions such that the response operators constructed using finer and finer partitions converge to to the actual corresponding response operators constructed on the SRB measure. Having in mind the Ulam method, the conditions are stricter than what needed in order to have convergence of the unperturbed measure, the basic reason being that Ruelle response operators correspond to nontrivial observables. One expects better convergence if the self-refining Markov partitions of the system are considered when constructing the finite state approximations.
Our results can be thought as intermediate steps at finite precision leading to the correct response formulas in the limit. One needs to add as a caveat that going from finite state to functional spaces is far from trivial and requires a high degree of mathematical precision, which is beyond the scopes of this paper. Nonetheless, the finite construction proposed here seems to somehow point at why some important mathematical issues emerge when the Perron-Frobenius operator formalism is considered in a continuum setting. In particular, the need for selecting suitable norms for vectors and linear operators in finite dimension points to the complex requirements in terms of functional spaces described in e.g. [51].
Interestingly, we can use the formulas obtained for finite state Markov processes to study the impact of perturbations to continuous time dynamical systems, after making a suitable identification between the considered transition matrices and the evolution operators for measures and observables. This operation is straightforward because there is a simple linear exact relation between the perturbation in the vector flow of the dynamical system and the perturbation in the Perron-Frobenius operator when infinitesimal time intervals are considered. As a result, we are able to derive in a very simple way previous formulas obtained studying the perturbations to the transfer operator as well as the original expressions proposed by Ruelle for the linear and higher order perturbations in the expectation values of observables. Using the results obtained in the finite state case, we propose a formula for the radius of expansion of the perturbative theory.
One can envision that in the case the underlying dynamics is discrete, there is not such a one-to-one correspondence between perturbations to the vector field and perturbations to the Markov transition matrix. This can be easily checked when constructing the perturbed Perron-Frobenius operator resulting from adding a correction to the vector field, which results into changes in the Perron-Frobenius operator at all orders in . Therefore, the perturbative expansion is different in the two cases. Agreement is instead found in the limit → 0, or, more practically, when we retain only the linear terms in perturbative expansion, i.e. when aiming only at the linear response function.
Future investigations will try, on the one side, to have a sharper mathematical look at the problem of going from finite to infinitely small partitions of the phase space, and, on the other side, to delve in the numerical study of the effectiveness and efficiency of the proposed tools. Apart from testing the results on specific finite state Markov systems, we will test how robust the proposed methods are when studying finite state Markov processes that have been empirically constructed from time series of observations or of numerical simulations of highdimensional complex systems. One may be led to hoping that it could be possible to have an accurate representation of the response of a high dimensional system to perturbations by constructing a smart finite state model well suited to studying specific observables of interest. Of course, in order to deal with the curse of dimensionality, one would like to be able to go beyond the Ulam method and deal with finite partition of reduced phase spaces where projection is applied on many or even most dimensions.
Our formulas may address the now long-standing problem of constructing suitable algorithms for studying the response of chaotic systems to perturbations. It is extremely hard to construct an algorithm for computing the (linear) response theory directly on the flow, because serious problems emerge when considering the contributions coming from the unstable directions in the tangent space. This might have great relevance for studying problems, like climate dynamics, where a direct construction of the response operator is especially challenging and slightly indirect methods have to be used [35] and a lot of effort has been devoted to defining the so-called atmospheric regimes and predicting their response to forcings [56].