ergodic and foliated kernel-differentiation method for linear responses of random systems

The averaged statistic of a dynamical system is of central interest in applied sciences. The average can be taken in two ways. When the system is random, we can average with respect to the randomness. On the other hand, if the system runs for a long-time, we can average with respect to time. The long-time average can also be defined for deterministic systems, and the existence of the limit measure is called the physical measure or SRB measure [43, 38, 6]. If the system is both random and long time, the existence of the limit measure (which we shall also call the physical measure) is discussed in textbooks such as [11].

We are interested in the the linear response, which is the derivative of the averaged observable with respect to the parameters of the system. The linear response is a fundamental tool for many purposes. For example, in aerospace design, we want to know the how small changes in the geometry would affect the average lift of an aircraft, which was answered for non-chaotic fluids [23] but only partially answered for not-so-chaotic fluids [27]. In climate change, we want to ask what is the temperature change caused by a small change of $CO_{2}$ level [16, 13]. In machine learning, we want to extend the conventional backpropagation method to cases with gradient explosion; current practices have to avoid gradient explosion but that is not always achievable [32, 22].

There are three basic methods for expressing and computing the linear response: the path perturbation method, the divergence method, and the kernel-differentiation method. The path perturbation method averages the path perturbation over a lot of orbits; this is also known as the ensemble method or stochastic gradient method (see for example [25, 12, 26]). However, when the system is chaotic, that is, when the pathwise perturbation grows exponentially fast with respect to time, the method becomes too expensive since it requires many samples to average a huge integrand. Proofs of such formula for the physical measure of hyperbolic systems were given in for example [39, 9, 24].

The divergence method is also known as the transfer operator method, since the perturbation of the measure transfer operator is some divergence. Traditionally, the divergence method is functional, and the formula is not pointwisely defined when the invariant measure is singular, which is common for systems with a contracting direction (see for example [19, 5]). This difficulty can not be circumvented by transformations such as Livisic theorem in [19]. Hence, the traditional divergence formulas can not be computed by a Monte-Carlo approach, so the computation is cursed by dimensionality (see for example [15, 42, 20, 2, 35]).

For chaotic deterministic systems, to have a formula without exponentially growing or distributive terms, we need to blend the two formulas, but the obstruction was that the split lacks smoothness. The ‘fast response’ linear response formulas solve this difficulty for deterministic hyperbolic systems. It is composed of two parts, the adjoint shadowing lemma [29] and the equivariant divergence formula [30]. The fast response formula is a pointwise defined function and has no exponentially growing terms. The continuous-time versions of fast response formulas can be found in [29, 31]. The fast response formula is the ergodic theorem for linear response: it computes $2u+2$ $M$-dimensional vectors on each point of an orbit, and then the average of some pointwise functions of these vectors. Here $M$ is the phase space dimension, and $u$ is the number of unstable directions.

The kernel-differentiation formula works only for random system. This method is more popularly known as the likelihood ratio method or the Monte-Carlo gradient method in probability context (see [37, 36, 17, 34]), the idea (dividing and multiplying the local density) is also used in Langevin sampling and diffusion models. However, such names are no longer distinguishing in the larger context: for example, the likelihood ratio idea is also useful in the divergence method. We apologize for the inconvenience in naming. Proofs of the kernel-differentiation method was given in [21, 3, 1]. This method is more powerful than the previous two for random dynamics, since annealed linear response can exist when quenched linear response does not [10]. However, as we shall see, inevitably, this method is expensive when the randomness is small, so it is hard to use this approach as a very good approximation to linear responses of deterministic systems.

In terms of numerics, the kernel-differentiation method naturally allows Monte-Carlo type sampling, so is efficient in high-dimensions. The formula does not involve Jacobian matrices hence does not incur propagations of vectors or covectors (parameter gradients), so it is not hindered by undesirable features of the Jacobian of the deterministic part, such as singularity or non-hyperbolicity. The work in [33] solves a Poisson equation for a term in the formula, which reduces the numerical error, but can be expensive for high-dimensions.

Another issue in previous likelihood ratio method is that the scale of the estimator is proportional to the orbit length. This leads to two undesirable consequences:

• •

  The requested number of sample orbits increases with orbit length.

• •

  We have to restart a new orbit very often. For each orbit, the observable function is only evaluated at the last step, while the earlier steps are somewhat wasted waiting for the orbits to land onto the physical measure.

This paper shows that, for physical measures, we can resolve these issues by invoking the decay of correlation and the ergodic theorem.

Moreover, sometimes we are interested in the case where some directions in the phase space have special importance, such as the level sets in Hamiltonian system. For these cases, the perturbation in the dynamics and the added noise are sometimes aligned to a subspace. This paper extends the kernel-differentiation method to foliated spaces.

Finally, for pedagogical purposes, it is beneficial to derive the kernel-differentiation method from a local or microscopic perspective. This will be somewhat more intuitive, and the above extensions will be more natural.

We consider the random dynamical system

$$X_{n+1}=f(\gamma,X_{n})+Y_{n+1}.$$

Let $h_{n}$ and $p_{n}$ be the probability density for $X_{n}$ and $Y_{n}$. We need the expressions of $p_{n}$; but we do not need expressions of $h_{n}$. We denote the density of the physical measure by $h$, which is the limit measure under the dynamics. Here $\gamma$ is the parameter controlling the dynamics and hence $h_{n}$ and $h$; by default $\gamma=0$. We denote the perturbation $\delta(\cdot):=\partial(\cdot)/\partial\gamma|_{\gamma=0}$, and when $f$ is invertible, define

$${\tilde{f}}:=f_{\gamma}\circ f^{-1}.$$

$\Phi$ is a fixed $C^{2}$ observable function.

As a warm-up, through basic calculus, section 2.2 rederives the kernel-differentiation formula for $\delta\int\Phi(x_{1})h(x_{1})dx_{1}$, the linear response for 1-step. Section 2.3 gives a pictorial explanation for the main lemma. In this simplest setting, we can see more easily how the kernel-differentiation method relates to other methods, and why it can be generalized to foliated situations. Section 2.4 rederives the formula for finitely many steps, which was known as the likelihood ratio method and was given in [37, 36, 17]. Our derivation use the measure-transfer-operator point of view, which is more convenient for our extensions.

The first contribution of this paper, given in section 2.5, is the ergodic kernel-differentiation formula for physical measures. We invoke decay of correlations to constrain the scale of the estimator, and the ergodic theorem to reduce the required data to one long orbit. Hence, we only need to spend one spin-up time to land onto the physical measure, and then all data are fully utilized. This improves the efficiency of the algorithm. More specifically, we prove

{theorem}

[orbitwise kernel-differentiation formula for $\delta h$] Under assumptions 1 and 2, we have

$$\begin{split}\delta\int\Phi(x)h(x)dx\,\overset{\mathrm{a.s.}}{=}\,\lim_{W\rightarrow\infty}\lim_{L\rightarrow\infty}-\frac{1}{L}\sum_{n=1}^{W}\sum_{l=1}^{L}\left(\Phi(X_{n+l})-\Phi_{avg}\right)\,\delta f_{\gamma}X_{l}\cdot\frac{dp}{p}(Y_{1+l})\end{split}$$

almost surely according to the measure obtained by starting the dynamic with $X_{0},Y_{1},Y_{2},\cdots\sim h\times p\times p\times\cdots$. Here $\Phi_{avg}:=\int\Phi h$.

Here $W$ and $L$ are two separate and sequential limits, so we can typically pick $W\ll L$. This makes the estimator smaller than previous likelihood ratio formulas, where $W=L$. Our formula runs on only one orbit; hence, it is faster than previous likelihood ratio methods. Section 2.6 generalizes the result to cases where the randomness depends on locations and parameters.

The second contribution of this paper is that section 3 extends the kernel-differentiation method to the case where the phase space is partitioned by submanifolds, and the added noise and the perturbation are along the submanifolds (hence the noise is ‘singular’ with respect to the Lebesgue measure). More specifically, Let $F=\{F_{\alpha}\}_{\alpha\in A}$ be a co-dimension $c$ foliation on $\mathcal{M}$, which is basically a collection of $(M-c)$-dimensional submanifolds partitioning $\mathcal{M}$, and each submanifold is called a ‘leaf’. Note that we do not require $f$ to respect the foliation, that is, $x$ and $fx$ can be one two different leaves. Let $p(z,x_{n+1})$ be the density of probability of $X_{n+1}$, conditioned on $f(X_{n})=z$, with respect to the Lebesgue measure on $F_{\alpha}(z)$. The structure of section 3 is similar to section 2, and the main theorem is

{theorem}

[centralized foliated kernel-differentiation formula for $\delta h_{T}$ of time-inhomogeneous systems, comparable to Remark] Under assumptions 1 and 3 for all $n$, we have

$$\begin{split}\delta\int\Phi(x_{T})d\mu_{T}(x_{T})=\mathbb{E}\left[\left(\Phi_{T}(x_{T})-\Phi_{avg,T}\right)\sum_{m=1}^{T}\left(\delta{\tilde{f}}_{m}(z_{m})\cdot\frac{d_{z}p_{m}(z_{m},x_{m})}{p_{m}(z_{m},x_{m})}\right)\right]\\ :=\int_{x_{0}\in\mathcal{M}}\int_{x_{1}\in F^{1}_{\alpha}(z_{1})}\cdots\int_{x_{T}\in F^{T}_{\alpha}(z_{T})}\left(\Phi_{T}(x_{T})-\Phi_{avg,T}\right)\\ \sum_{m=1}^{T}\left(\delta{\tilde{f}}_{m}(z_{m})\cdot\frac{d_{z}p_{m}(z_{m},x_{m})}{p_{m}(z_{m},x_{m})}\right)p_{T}(z_{T},x_{T})dx_{T}\cdots p_{1}(z_{1},x_{1})dx_{1}d\mu_{0}(x_{0}).\end{split}$$

Here $\Phi_{avg,T}:=\int\Phi_{T}d\mu_{T}$. Note that here $F,{\tilde{f}},f,p$ can be different for different steps.

Section 4 considers numerical realizations. Section 4.1 gives a detailed list for the algorithm. Note that the foliated case uses the same algorithm; the difference is in the set-up of the problem. Section 4.2 illustrates the ergodic version of the algorithm on an example where the deterministic part does not have a linear response, but adding noise and using the kernel-differentiation algorithm give a reasonable reflection of the relation between the long-time averaged observable and the parameter. Section 4.3 illustrates the foliated version of the algorithm on a unstable neural network with 51 layers $\times$ 9 neurons, where the perturbation is along the foliation. We show that adding a low-dimensional noise along the foliation induces a much smaller error than adding noise in all directions. And the foliated kernel-differentiation can then give a good approximation of the true parameter-derivative.

In the discussion section, section 5.1 gives a rough cost-error estimation of the kernel-differentiation algorithm. This helps to set some constants in the algorithm; it also shows that the small-noise limit is still very expensive. Section 5.2 proposes a program which combines the three formulas to obtain a good approximation of linear responses of deterministic systems.

The paper is organized as follows. Section 2 derives the kernel-differentiation formula where the noise is in all directions. Section 3 derives the formula for foliated perturbations and noises. Section 4 gives the procedure list of the algorithm based on the formulas, and illustrates on two numerical examples. Section 5 discusses the relation to deterministic linear responses.

Let $f_{\gamma}$ be a family of $C^{2}$ map on the Euclidean space $\mathbb{R}^{M}$ of dimension $M$ parameterized by $\gamma$. Assume that $\gamma\mapsto f_{\gamma}$ is $C^{1}$ from $\mathbb{R}$ to the space of $C^{2}$ maps. The random dynamical system in this paper is given by

(1)

$$\begin{split}X_{n}=f_{\gamma}(X_{n-1})+Y_{n},\quad\textnormal{where}\quad Y_{n}{\,\overset{\mathrm{iid}}{\sim}\,}p.\end{split}$$

The default value is $\gamma=0$, and we denote $f:=f_{\gamma=0}$. In this section, we consider the case where $p$ is any fixed $C^{2}$ probability density function. The next section considers the case where $p$ is smooth along certain directions but singular along other directions.

Moreover, except for section 1.2, our results also apply to time-inhomogeneous cases, where $f$ and $p$ are different for each step. More specifically, the dynamic is given by

(2)

$$\begin{split}X_{n}=f_{\gamma,n-1}(X_{n-1})+Y_{n},\quad\textnormal{where}\quad Y_{n}\sim p_{n}.\end{split}$$

Here $X_{n}\in\mathbb{R}^{M_{n}}$ are not necessarily of the same dimension. We shall exhibit the time dependence in Remark and also in the numerics section. On the other hand, for the infinite time case in section 1.2, we want to sample by only one orbit; this requires that $f$ and $p$ be repetitive among steps.

We define $h_{T}$ as the density of pushing-forward the initial measure $h_{0}$ for $T$ steps.

$$\begin{split}h_{T}:=(L_{p}L_{f_{\gamma}})^{T}h_{0}.\end{split}$$

Here $h_{T}=h_{T,\gamma}$ depends on $\gamma$ and also $h_{0}$; we sometimes omit the subscript of $\gamma$. Here $L_{f}$ is the measure transfer operator of $f$, which are defined by the integral equality

(3)

$$\begin{split}\int\Phi(x)\,(L_{f}h)(x)dx:=\int\Phi(fx)h(x)dx.\end{split}$$

Here $\Phi\in C^{2}(\mathbb{R}^{M})$ is any fixed observable function. When $f$ is bijective, there is an equivalent pointwise definition, $(L_{f}h)(x):=\frac{h}{|Df|}(f^{-1}x)$.

For density $q$, $L_{p}q$ is pointwisely defined by convolution with density $p$:

(4)

$$\begin{split}L_{p}q(x)=\int q(x-y)p(y)dy=\int q(z)p(x-z)dz.\end{split}$$

We shall also use the integral equality

(5)

$$\begin{split}\int\Phi(x)L_{p}q(x)dx=\int p(y)\left(\int\Phi(x)q(x-y)dx\right)dy\\ \,\overset{z=x-y}{=}\,\iint\Phi(y+z)q(z)p(y)dzdy.\end{split}$$

If we want to compute the above integration by Monte-Carlo method, then we should generate i.i.d. $Y_{l}=y_{l}$ and $Z_{l}=z_{l}$ according to density $p$ and $q$, then compute the average of $\Phi(y_{l}+z_{l})$.

To consider differentiability of $h$ with respect to $\gamma$, we make the following simplifying assumptions on differentiability of the basic quantities.

The linear response formula for finite time $T$ is an expression of $\delta h_{\gamma,T}$ by $\delta f_{\gamma}$, and

$$\begin{split}\delta(\cdot):=\left.\frac{\partial(\cdot)}{\partial\gamma}\right|_{\gamma=0}.\end{split}$$

Here $\delta$ may as well be regarded as small perturbations. Note that by our notation,

$$\begin{split}\delta f_{\gamma}(x)\in T_{fx}\mathbb{R}^{M};\end{split}$$

that is, $\delta f_{\gamma}(x)$ is a vector at $fx$. The linear response formula for finite-time is given by the Leibniz rule,

(6)

$$\begin{split}\delta h_{T}=\delta(L_{p}L_{f_{\gamma}})^{T}h_{0}=\sum_{n=0}^{T-1}(L_{p}L_{f})^{n}\,\delta(L_{p}L_{f_{\gamma}})\,(L_{p}L_{f})^{T-1-n}h_{0}.\end{split}$$

We shall give an expression of the perturbative transfer operator $\delta(L_{p}L_{f_{\gamma}})$ later.

For the case where $T\rightarrow\infty$, we define the density $h$ of the (ergodic) physical measure $\mu$, which is the unique limit of $h_{T}$ in the weak$*$ topology for any smooth initial density $h_{0}$:

$$\begin{split}\lim_{T\in\mathbb{N},\,T\rightarrow\infty}(L_{p}L_{f_{\gamma}})^{T}h_{0}=:h_{\gamma},\end{split}$$

note that $h$ depends on $\gamma$ but not on $h_{0}$. The physical measure is also called the SRB measure or the stationary measure. We also define the correlation function, $\mathbb{E}[\varphi(X_{n})\psi(X_{0},X_{1})]$, where $X_{0}$ is distributed according to the physical measure $\mu$, and $X_{n}$ is distributed according to the dynamics in equation 1. For this section,

$$\mathbb{E}_{\gamma}[\varphi(X_{n})\psi(X_{0},X_{1})]:=\int\varphi(x_{n})\psi(x_{0},x_{1})h_{\gamma}(x_{0})dx_{0}p(y_{1})dy_{1}\cdots p(y_{n})dy_{n}.$$

We shall make the following assumptions for infinitely many steps:

Assumption 2 can be proved from more basic assumptions, which are typically more forgiving than the hyperbolicity assumption for deterministic cases. For example, under assumption 1, the proof of unique existence of the physical measure, where the system has additive noise with smooth density, can be found in probability textbooks such as [11]. [14] proves the decay of correlations when, roughly speaking, $f_{\gamma}$ has bounded variance and is mixing and that the transfer operator is bounded. The technical assumptions and rigorous proof for linear responses can also be found in [21]. Our main purpose is to derive the ergodic kernel-differentiation formula with the detailed expression in section 2.2. The formula seems to be the same in most cases where linear response exists.

We can give a pointwise formula for $\delta(L_{p}L_{f_{\gamma}}q)(x)$, which seems to be a long-known wisdom. An intuitive explanation of the following lemma is given in section 2.3. Then we give the integral version which can be computed by Monte Carlo algorithms. An intuition of the lemmas are given in section 2.3.

{lemma}

[local one-step kernel-differentiation formula] Under assumption 1, then for any density $h_{0}$

$$\begin{split}\delta(L_{p}L_{f_{\gamma}}h_{0})(x_{1})=\int-\delta f_{\gamma}(x_{0})\cdot dp(x_{1}-fx_{0})\,h_{0}(x_{0})dx_{0}.\end{split}$$

Here $\delta:=\partial/\partial\gamma|_{\gamma=0}$, and $\delta f_{\gamma}(x_{0})\cdot dp(x_{1}-fx_{0})$ is the derivative of the function $p(\cdot)$ at $x_{1}-fx_{0}$ in the direction $\delta f_{\gamma}(x_{0})$, which is a vector at $fx_{0}$.

{lemma}

[integrated one-step kernel-differentiation formula] Under assumption 1, then

$$\begin{split}\int\Phi(x_{1})\delta(L_{p}L_{f_{\gamma}}h_{0})(x_{1})dx_{1}=\iint-\Phi(x_{1})\delta f_{\gamma}(x_{0})\cdot\frac{dp(y_{1})}{p(y_{1})}\,h_{0}(x_{0})dx_{0}\,p(y_{1})dy_{1}.\end{split}$$

Here $x_{1}$ in the right expression is a function of the dummy variables $x_{0}$ and $y_{1}$, that is, $x_{1}=y_{1}+fx_{0}$.

We give an intuitive explanation to section 2.2 and Remark. This will help us to generalize the method to the foliated situation in section 3. First, we shall adopt a more intuitive but restrictive notation. Define ${\tilde{f}}(x,\gamma)$ such that

$$\begin{split}{\tilde{f}}(f(x),\gamma)=f_{\gamma}(x).\end{split}$$

In other words, we write $f_{\gamma}$ as appending a small perturbative map ${\tilde{f}}$ to $f:=f_{\gamma=0}$. Here ${\tilde{f}}$ is the identity map when $\gamma=0$. Note that ${\tilde{f}}$ can be defined only if $f$ satisfies

$$\begin{split}f_{\gamma}(x)=f_{\gamma}(x^{\prime})\quad\textnormal{whenever}\quad f(x)=f(x^{\prime}).\end{split}$$

For example, when $f$ is bijective then we can well-define ${\tilde{f}}$. Hence, this new notation is more restrictive than what we used in other parts of this section, but it lets us see more clearly what happens during the perturbation.

With the new notation, the dynamics can be written now as

$$\begin{split}X_{n}={\tilde{f}}f(X_{n-1})+Y_{n},\quad\textnormal{where}\quad Y_{n}{\,\overset{\mathrm{iid}}{\sim}\,}p.\end{split}$$

By this notation, only $f$ changes time step, but ${\tilde{f}}$ and adding $Y$ do not change time. In this subsection we shall start from after having applied the map $f$, and only look at the effect of applying ${\tilde{f}}$ and adding noise: this is enough to account for the essentials of section 2.2 and Remark. Roughly speaking, the $q_{n}$ we use in the following is in fact $q_{n}:=L_{f}h_{n-1}$; $q_{n}$ is the density of $z_{n}:=fx_{n-1}$, so $z_{n}+y_{n}=x_{n}$, and we omit the subscript for time $0$.

With the new notation, section 2.2 is essentially equivalent to

(7)

$$\begin{split}\delta(L_{p}L_{{\tilde{f}}}q)(x)=\int-\delta{\tilde{f}}z\cdot dp(x-z)\,q(z)dz.\end{split}$$

We explain this intuitively by figure 1. Let ${\tilde{z}}:={\tilde{f}}z$ be distributed according to ${\tilde{q}}$, then $L_{p}{\tilde{q}}$ is obtained by first attach a density $p$ to each ${\tilde{z}}$, then integrate over all ${\tilde{z}}$. $L_{p}L_{\tilde{f}}q$ is obtained by first move $z$ to ${\tilde{z}}$ and then perform the same procedure. Hence, $\delta L_{p}L_{\tilde{f}}q(x)$ is first computing $\delta p_{{\tilde{f}}z}(x)$ for each $z$, then integrate over $z$. Let

$$\begin{split}p_{{\tilde{f}}z}(x):=p(x-{\tilde{f}}z)\end{split}$$

be the density of the noise centered at ${\tilde{f}}z$. So

$$\begin{split}\delta p_{{\tilde{f}}z}(x)=-dp_{{\tilde{f}}z}(x)\cdot\delta{\tilde{f}}z(x)=-dp(x-{\tilde{f}}z)\cdot\delta{\tilde{f}}z.\end{split}$$

Here $\delta{\tilde{f}}z(x)$ in the middle expression is the horizontal shift of $p_{{\tilde{f}}z}$’s level set previously located at $x$. Since the entire Gaussian distribution is parallelly moved on the Euclidean space $\mathbb{R}^{M}$, so $\delta{\tilde{f}}z(x)=\delta{\tilde{f}}z$ is constant for all $x$. Then we can integrate over $z$ to get equation 7.

Using ${\tilde{f}}$ notation, Remark is equivalent to

$$\begin{split}\delta\int\Phi(x)(L_{p}L_{\tilde{f}}q)(x)dx=\iint-\Phi(x)\delta{\tilde{f}}(z)\cdot\frac{dp(y)}{p(y)}\,q(z)dz\,p(y)dy,\end{split}$$

where $x$ on the right is a function $x=z+y$. Intuitively, this says that the left side equals to first compute $\delta\int\Phi(z+y)p_{{\tilde{f}}z}(y)dy$ for each $z$, then integrate over $z$. The integration on the left uses $x$ as dummy variable, which is convenient for the transfer operator formula above. But it does not involve a density for $x$, so is not immediately ready for Monte-Carlo. The right integration is over $q(z)$ and $p(y)$, which are easy to sample.

It is important that we differentiate only $p$ but not $f$. In fact, our core intuitions are completely within one time step, and do not even involve $f$. Hence, we can easily generalize to cases where $f$ is bad, for example, when $f$ is not bijective, or when $f$ is not hyperbolic: these are all difficult for deterministic systems.

We use Remark to get integrated formulas for $h_{T}:=(L_{p}L_{f_{\gamma}})^{T}h_{0}$ and its perturbation $\delta h_{T}$, which can be sampled by Monte-Carlo type algorithms. Note that here $h_{0}$ is fixed and does not depend on $\gamma$. This was previously known as the likelihood ratio method or the Monte-Carlo gradient method and was given in [37, 36, 17]. The magnitude of the integrand grows as $O(\sqrt{T})$ after the centralization, where $T$ is the number of time steps.

{lemma}

Under assumption 1, for any $C^{2}$ (not necessarily positive) densities $h_{0}$, any $n\in\mathbb{N}$,

$$\begin{split}\int\Phi(x_{n})((L_{p}L_{f})^{n}h_{0})(x_{n})dx_{n}=\int\Phi(x_{n})h_{0}(x_{0})dx_{0}p(y_{1})dy_{1}\cdots p(y_{n})dy_{n}.\end{split}$$

Here the $x_{n}$ on the left is a dummy variable; whereas $x_{n}$ on the right is recursively defined by $x_{m}=f(x_{m-1})+y_{m}$, so $x_{n}$ is a function of the dummy variables $x_{0},y_{1},\cdots,y_{n}$.

{theorem}

[kernel-differentiation formula for $\delta h_{T}$] Under assumption 1, then

$$\begin{split}\delta\int\Phi(x_{T})h_{T}(x_{T})dx_{T}=-\int\Phi(x_{T})\sum_{m=0}^{T-1}\left(\delta f_{\gamma}x_{m}\cdot\frac{dp}{p}(y_{m+1})\right)h_{0}(x_{0})dx_{0}(pdy)_{1\sim T}.\end{split}$$

Here $dp$ is the differential of $p$, $(pdy)_{1\sim T}:=p(y_{1})dy_{1}\cdots p(y_{T})dy_{T}$ is the independent distribution of $Y_{n}$’s, and $x_{m}$ is a function of dummy variables $x_{0},y_{1},\cdots,y_{m}$.

Note that subtracting any constant from $\Phi$ and does not change the linear response. We prove a more detailed statement in a Remark. Hence, we can centralize $\Phi$, i.e. replacing $\Phi(\cdot)$ by $\Phi(\cdot)-\Phi_{avg,T}$, where the constant

$$\begin{split}\Phi_{avg,T}:=\int\Phi(x_{T})h_{T}(x_{T})dx_{T}.\end{split}$$

We sometimes centralize by subtracting $\Phi_{avg}:=\int\Phi hdx$. The centralization reduces the amplitude of the integrand, so the Monte-Carlo method converges faster: this is also a known result and was discussed in for example [18].

{lemma}

[free centralization] Under assumption 1, for any $m\geq 0$, we have

$$\begin{split}\int\delta f_{\gamma}x_{m}\cdot\frac{dp}{p}(y_{m+1})h_{0}(x_{0})dx_{0}(pdy)_{1\sim m+1}=0.\end{split}$$

For the convenience of computer coding, we explicitly rewrite this theorem into centralized and time-inhomogeneous form, where $\Phi$, $f$, and $p$ are different for each step. This is the setting for many applications, such as finitely-deep neural networks. The following proposition can be proved similarly to our previous proofs. Note that we can reuse a lot of data should we also care about the perturbation of the averaged $\Phi_{n}$ for other layers $1\leq n\leq T$.

In the following proposition, let $\Phi_{T}:\mathbb{R}^{M_{T}}\rightarrow\mathbb{R}$ be the observable function defined on the last layer of dynamics. Let $h_{T}$ be the pushforward measure given by the dynamics in equation 2, that is, defined recursively by $h_{n}:=L_{p_{n}}L_{f_{\gamma,n-1}}h_{n-1}$; Let $(pdy)_{1\sim T}:=p_{1}(y_{1})dy_{1}\cdots p_{T}(y_{T})dy_{T}$ be the independent but not necessarily identical distribution of $Y_{n}$’s.

{proposition}

[centralized kernel-differentiation formula for $\delta h_{T}$ of time-inhomogeneous systems] If $f_{\gamma,n}$ and $p_{n}$ satisfy assumption 1 for all $n$, then

$$\begin{split}&\delta\int\Phi_{T}(x_{T})h_{T}(x_{T})dx_{T}\\ =&-\int\left(\Phi_{T}(x_{T})-\Phi_{avg,T}\right)\sum_{m=0}^{T-1}\left(\delta f_{\gamma,m}x_{m}\cdot\frac{dp}{p}(y_{m+1})\right)h_{0}(x_{0})dx_{0}(pdy)_{1\sim T}.\end{split}$$

Here $\Phi_{avg,T}:=\int\Phi_{T}h_{T}$.

For the perturbation of physical measures, the kernel-differentiation formula for $\delta h$ can further take the form of long-time average on an orbit. This reduces the magnitude of the integrand so the convergence is faster. Here we can apply the ergodic theorem, forget details of initial distribution, and sample $h$ by points from one long orbit. Note that this subsection requires that the dynamics is time-homogeneous ($f$ and $p$ are the same for each step) to invoke ergodic theorems.

{lemma}

[kernel-differentiation formula for physical measures] Under assumptions 1 and 2,

$$\begin{split}\delta\int\Phi(x)h(x)dx=\int\Phi(x)\delta h(x)dx\\ =-\lim_{W\rightarrow\infty}\sum_{n=1}^{W}\int\Phi(x_{n})\left(\delta f_{\gamma}x_{0}\cdot\frac{dp}{p}(y_{1})\right)h_{\gamma=0}(x_{0})dx_{0}(pdy)_{1\sim n}\end{split}$$