The entropy of a strictly positive probability vector \( \mathbf{u}\) is defined as:

A uniform probability vector represents a thermodynamic system where all microstates have the same probability of manifesting, so little information is known with regards to the actual state of the system. This means that the system has high entropy, in fact, for the uniform distribution, maximum entropy. The opposite happens when the probability vector concentrates in a single microstate, yielding maximum information about the system, corresponding to minimum entropy.

The entropy of the perturbed system, characterized by the perturbed invariant density \( \mathbf{v}\) — solving Eq. (2)— is to linear order in \( \varepsilon\) , \( \mathbf{u}+\varepsilon\mathbf{v}_1\) . Then, we can expand the entropy functional as:

where we have used the fact that \( \mathbf{u}+\varepsilon\mathbf{v}_1\) must add to unity. More precisely:

Thus, to optimize the entropy of the system, we search for a perturbation matrix \( \mathbf{P}\) that maximizes the functional:

where \( \mathrm{f}_i=\log(u_i)\) and \( \mathbf{P}\) satisfies the constrains (C1)-(C2)-(C3) listed below Eq. (12). This is consistent with the fact that the entropy \( H(\mathbf{u})\) can be interpreted as the expectation value under the distribution \( \mathbf{u}\) of the observable \( \log(\mathbf{u})\) , where the logarithm acts componentwise.

This optimization problem was solved in [25, Section 4] when optimizing the linear response of the expectation value of a generic observable. There the authors explicitly solve the Lagrange multiplier problem to find that the matrix \( \mathbf{P}\) that maximizes Eq. (16) is given component-wise as:

where \( m_j\) is the number of non-zero elements in the \( j\) th column of \( \mathbf{M}\) and \( \nu\) is a normalizing factor modulo a sign dependant on the Hessian matrix of the Lagrangian function. A similar derivation with another constraint is done later in the text in Section 3.3, so more precise details are given then.

While optimizing entropy in absolute terms can be useful in certain cases, it is also important to determine the perturbation that creates the most information relative to the unperturbed measure as discussed in [26]. In other terms, we want to find the perturbation matrix \( \mathbf{P}\) that maximizes the Kullback-Liebler (KL) divergence. The KL-divergence of the probability vector \( \mathbf{v}\) with respect to \( \mathbf{u}\) is given by:

Using the expression derived in Eq.(9) for the response \( \mathbf{v}_1\) , we evaluate the KL-divergence between the leading order perturbed density \( \mathbf{u} + \varepsilon \mathbf{v}_1\) and the unperturbed invariant vector \( \mathbf{u}\) yielding— c.f. [26, Section 2.3]—:

The last equation of the previous display can be recast as a matrix product. Indeed, define \( \mathbf{D} = (D_{ij})\) as a diagonal \( N\times N \) matrix such that \( D_{ii} = u^{1/2}_i>0\) , for \( i=1,…,N\) . Then, we have

which reduces to a weighted 2-norm with respect to the vector \( \mathbf{u}^{-1} = (1/u_i)\) .

The target is to obtain the matrix \( \mathbf{P}\) such that the quantity in Eq. (20) is optimized with the constraints (C1)-(C2)-(C3) above. An entry-wise solution similar to Eq. (17) cannot be obtained since the 2-norm does not respect the sparsity of the matrix \( \mathbf{M}\) . Nevertheless, the problem can be vectorized so that the solutions can be determined by the singular value decomposition. Indeed, the matrix inside the 2-norm in Eq. (20) can be writen as:

where \( \mathbf{A}=\mathbf{u}^{\top}\otimes\left[ \mathbf{D}^{-1}\mathbf{G}\right]\) with “\( \otimes\) " denoting the Kronecker outer product and \( \widehat{\mathbf{P}}\) is a \( N^2\) -dimensional column vector resulting from the vectorization of the matrix \( \mathbf{P}\) . If there were no further constraints to \( \mathbf{P}\) , optimizing the 2-norm of \( \mathbf{D}^{-1}\mathbf{G}\mathbf{P}\mathbf{u}\) is obtained from the singular value decomposition of \( \mathbf{A}\) . In this case, the singular vector associated with the largest singular value is, upon reordering, the desired \( \mathbf{P}\) . With the constraints (C2)-(C3), the authors in [25] show that the problem is equivalent to projecting onto the nullspace of a specific matrix, with an orthonormal basis arranged in the columns of a matrix \( \mathbf{B}\) . We omit here the explicit form of \( \mathbf{B}\) , we refer the reader to [25]— c.f. see Appendix B—. The maximization problem is, then, solved by the SVD of the matrix \( \mathbf{A}\mathbf{B}\) .

Entropy production in Markov matrices quantifies the irreversibility of stochastic processes described by these matrices. Entropy production thus measures how far a system deviates from detailed balance, reflecting the arrow of time and the dissipation of energy. It connects microscopic stochastic dynamics with macroscopic thermodynamic behavior, providing insight of the second law of thermodynamics at a probabilistic level [39].

The entropy production of a (column-normalized) Markov matrix \( \mathbf{M}\) is defined here as [39]:

where \( s(\mathbf{M})\geq 0\) and \( s(\mathbf{M})=0\) when the matrix obeys the detailed balance property so that \( u_iM_{ji}/u_jM_{ij}\) for every \( i,j=1,…,N\) , making the logarithm vanish. The question we want to answer is: what is the matrix \( \mathbf{P}\) such that the Markov matrix \( \mathbf{M} + \varepsilon\mathbf{P}\) is closer to detailed-balance and has \( \mathbf{u}\) as an invariant probability vector? The direct candidate \( \mathbf{P}_r\) is the following:

where the matrix \( \frac{1}{2}\left( \mathbf{M}+\mathbf{D}\mathbf{M}^{\top}\mathbf{D}^{-1} \right)\) is the so-called additive reversibilization of \( \mathbf{M}\) and satisfies detailed balance by construction [40]. Furthermore, we have that \( (\mathbf{M}+\mathbf{P}_r)\mathbf{u} = \mathbf{u}\) and \( s(\mathbf{M}+\mathbf{P}_r)=0\) . The issue is that we aim at minimizing the distance to reversibility by small perturbations to the original chain \( \mathbf{M}\) , and \( \mathbf{P}_r\) might provoke significant changes in the chain as quantified by a large amplitude \( \|\mathbf{P}_r\|_F\gg \varepsilon\) . To solve such problem, we expand \( s(\mathbf{M}+\varepsilon \mathbf{P})\) perturbatively in terms of \( \varepsilon\) to obtain the leading order contribution of \( \mathbf{P}\) in the entropy production.

Because we aim at leaving the invariant measure unperturbed, \( \mathbf{u}\) must satisfy \( \left(\mathbf{M} + \varepsilon\mathbf{P} \right)\mathbf{u} = \mathbf{u}\) , and it follows that \( \mathbf{P}\mathbf{u} = \mathbf{0}\) . Hence, the matrix \( \mathbf{P}\) in the minimization problem should satisfy another constraint on top of (C1)-(C2)-(C3):

To obtain the leading order contribution of \( \mathbf{P}\) in powers of \( \varepsilon\) , we take the derivative of \( s(\mathbf{M}+\varepsilon\mathbf{P})\) with respect to \( \varepsilon\) evaluated at \( \varepsilon=0\) reads as:

where we have used the fact that \( \sum_{i=1}P_{ji}u_i =0 \) for every \( j=1,…,N\) and the matrix \( \mathbf{C} = (C_{ij})\) is a matrix whose entries are a function of those of \( \mathbf{M}\) and \( \mathbf{u}\) :

This way, the problem has the same linear and entry-wise structure as Eq. (16). So we consider the problem of minimizing the linear functional of Eq. (24) subject to the constraints (C1)-(C2)-(C3)-(C4). Because of the linearity and entry-wise structure of the problem, the sparsity constraint (C3) is automatically satisfied if the correct Lagrange multiplier equation is solved for the relevant \( P_{ij}\) , namely, where \( M_{ij}>0\) . To minimize Eq. (24), we employ the Lagrange multipliers in the same way as [25, Section 4], with the additional (C4) constraint. The corresponding Lagrangian function reads as:

where \( \mathbf{q}\) , \( \mathbf{r}\) and \( \nu\) are the Lagrangian multipliers. To maximize this Lagrangian over the parameters, we need to compute the gradient with respect to each entry of \( \mathbf{P}\) . Because \( \mathbf{P}\) must have the same sparsity mask as \( \mathbf{M}\) , one can impose the same sparsity mask to \( \mathbf{C}\) without changing the value of \( \sum_{i=1}^N\sum_{j=1}^NC_{ij}P_{ij}\) . We then compute the derivative of \( L\) with respect to \( P_{ij}\) and equalize it to zero:

from where we obtain:

Define \( Z = { (i,j)\in \mathbb{N}^{2} : M_{ij}>0 }\) , and \( Z_jĉ={ i \in \mathbb{N}: (i,j) \in Z }\) . From the previous display, we find:

This way we get a formula for \( r_j\) . We now plug the previous display into Eq. (28) to find the expression for \( P_{ij}\) :

we are left with finding \( q_i\) and \( \nu\) .

We now multiply Eq. (26) by \( u_j\) ,

and then we sum over \( j=1,…,N\) , and apply the condition \( \mathbf{P}\mathbf{u} = \mathbf{0}\) , for each value of \( k=1,…,N\) :

Next, we plug in the value of \( r_j\) obtained in Eq. (29):

Now we rearrange the previous equation into a linear equation on \( q_{k}\) as:

where \( \alpha_k\) and \( \beta_j\) are defined as:

Equation (34) can be recast in the following form:

where \( \sum_k \gamma_k = 1\) . Notice that the associated matrix is invertible. Indeed, Eq. (36) can be written as:

where \( \Xi\) is a diagonal matrix with \( \Xi_{ii}=\xi_i>0\) , \( \boldsymbol{\alpha} = (\alpha_i)\) and \( \boldsymbol{\gamma}=(\gamma_i)\) . Since \( \Xi\) is invertible and \( \mathbf{1}\boldsymbol{\gamma}^{\top}\) rank-one, the inverse of \( \Xi + \mathbf{1}\boldsymbol{\gamma}^{\top}\) is given by the Sherman-Morrison formula:

Finally \( \nu\) is chosen so that the Frobenius norm constraint (C1) is satisfied and the sign is chosen negative so that the Hessian matrix of the Lagrangian function is negative; see [25, Eq. (57)]. This guarantees that the matrix \( \mathbf{P}\) in Eq. (30) is a minimizer.

An alternative, yet equivalent, approach can be taken by constructing a suitable vector space that contains the matrices satisfying constraints (C2) and (C4). For the particular problem of entropy production minimization, we run both methods on randomly initialized sparse and diagonal-dominant Markov matrices to illustrate such equivalence and also show that the methods, indeed, minimize the entropy production functional; see Appendix B for details.

Once the optimal perturbation matrix \( \mathbf{P}\) is found, it is of fundamental importance to be able to reconstruct the corresponding perturbation field at the SDE level. In other words, we want to find the vector field such that adding it to Eq. (39) corresponds to the perturbed Markov chain system. Analytical solutions to this problem are expected to be untractable, so here a numerical algorithm is described. First of all, we want to investigate whether we can reconstruct the vector field \( \mathbf{F}\) in Eq. (39) from the Markov matrix discretization, \( \mathbf{M}\) , of \( e^{\tau\mathcal{L}}\) .

The generator \( \mathcal{L}\) of the semigroup is approximated by the logarithm of the Markov matrix normalized by the transition time \( \tau\) — as follows from the Spectral Mapping Theorem [41]—:

Because \( \mathbf{M}\) is column-stochastic, the columns of \( \mathbf{L}\) must add up to zero. The entries of the matrix \( \mathbf{L}\) represent the rate of probability fluxes out of each box in the discretization. Hence, the mean velocity is the sum over all these possible transitions weighted by their rates, given by \( \mathbf{L}\) , and displacement vectors. With this ansatz, we get the approximation formula for the vector field:

where \( \left[\mathbf{F}(\mathbf{c}_i) \right]_k\) is the \( k\) -th component of the vectorfield \( \mathbf{F}\) evaluated at box-center \( \mathbf{c}_i\) .

To illustrate the validity of this method we consider a non-symmetric one-dimensional double well potential equation:

where \( \alpha=-0.1\) , \( \sigma>0\) and \( W_t\) is a one dimensional Wiener process. The equation is integrated for \( 10^6\) time units with a time step of \( \mathrm{d}t=10^{-2}\) time units using an Euler-Maruyama scheme [49]. This gives a time series \( {x_k}_{k=1}^{K}\) of \( K=10^8\) elements with \( x_1=0\) . Then, the interval \( [-1.75,1.75]\) is subdivided into \( N=1024\) boxes of the same size. Such an interval contains most of the trajectory and avoids the least populated regions which can incur in numerical errors when applying Eq. (43). By taking \( \tau = 0.1\) , we estimate the Markov matrix from Eq. (41):

where “\( \#\) ” denotes the counting measure. One can now apply the matrix logarithm in Eq. (42) to obtain the flux matrix \( \mathbf{L}\) . The formula Eq. (43) is systematically applied to the centers of the \( N\) boxes covering the domain to obtain a drift reconstruction vector in \( \mathbb{R}^{N}\) . This reconstructed vector is ploted with the magenta squares in Fig. 1(A) against the the actual drift function \( \mathbf{F}(x)\) plotted in the black curve.

For illustration, we now enquire what is the perturbation \( \mathbf{P}\) of the matrix in Eq. (45) such that the KL-divergence is maximized. Intuitively, in order to amplify the KL-divergence, one expects that the forcing applied to Eq. (44) should equalize the wells that originally were unbalanced due to the term \( \alpha\) . To visualize the forcing, we need Eq. (43). First of all, we apply the algorithm of Section 3.2 to obtain \( \mathbf{P}\) . The matrix \( \mathbf{P}\) is by construction a probability rate flux, just as \( \mathbf{L}\) is in Eq. (42). Therefore, the drift reconstruction algorithm of Eq. (43) is readily applicable to obtain a vector \( \mathbf{g}\) that approximates the perturbation of the original drift \( \mathbf{F}(x)\) . This vector \( \mathbf{g}\) is plotted in red in Fig. 1(A) and, indeed, the bulk of the perturbation is on the side of the deepest well. The invariant measure of the unperturbed system is shown in the black curve of Fig. 1(B) and the perturbed invariant measure \( \mathbf{v}\) is shown in red, where \( \varepsilon = 0.1\) . The leading order correction of \( \mathbf{u}\) due to the perturbation \( \mathbf{P}\) is shown with the square dots. Additionally, we computed \( -\varepsilon \mathbf{G}\mathcal{D}(\mathbf{u}:\mathbf{g})\) , where \( \mathcal{D}\) is a centered-finite-difference derivative matrix. This resulting vector is the discretization of the leading order correction to the solutions of Eq. (40) due to an external perturbation. This vector is shown with the empty triangles. We observe that the applied perturbation falls into the linear response regime due to the accurate match between \( \mathbf{v}\) and the linear order corrections. Furthermore, we illustrate the numerical tractability of Eq. (43), which in spite of being sensitive to unsmoothness— see spikes in the vector \( \mathbf{g}\) of Fig. (1)(A)— the response preserves the accuracy and reveals the underlying differential operator structure of \( \mathbf{P}\) .

Pan/Zoom

In this section we analyse a stochastically perturbed Lanford map, which has also been investigated in the work of [25] using the Markov matrix framework. Here, instead, we employ the algorithms presented in Section 3 to analyse its response to perturbations that maximize, 1) the linear response of the entropy functional and, 2) the KL-divergence of the perturbed measure with respect to the unpertubed state. For each case, the linearity of the response and its convergence to its asymptotic value is investigated.

The stochastic Lanford map \( \mathfrak{L}:[0,1]\longrightarrow[0,1]\) is defined as:

where \( \xi\) is uniformly distributed on the interval \( [-0.1,0.1]\) . The Markov discretization is done as follows. First, the interval \( [0,1]\) is discretized into \( 1024\) boxes \( { B_i }_{i=1}^{1024}\) . Secondly, Eq. (46) is integrated for \( K= 10^8\) timesteps producing a timeseries \( {x_k}_{k=1}^{K}\) , where \( x_k\) is in \( [0,1]\) and \( x_1 = 0\) . Third, the Markov matrix \( \mathbf{M} = (M_{ij})\) describing the pushforward map— the discrete-time version of Eq. (41)— is constructed by:

The resulting matrix \( \mathbf{M}\) has a unique stationary probability vector \( \mathbf{u}\) and it is plotted in the black line of Fig. 2(A).

The linear response of the entropy function is maximized using formula (17) to obtain a perturbation matrix \( \mathbf{P}_1\) . To maximize the KL-divergence from \( \mathbf{u}\) , the functional Eq. (20) is maximized by a matrix \( \mathbf{P}_2\) . Then, the perturbation size of \( \varepsilon = 0.1\) is chosen so that \( \mathbf{M} + \varepsilon \mathbf{P}_i\) is a (perturbed) Markov matrix with perturbed invariant measure \( \mathbf{v}^{(i)}\) satisfying:

for \( i=1,2\) . The respective perturbed steady states are plotted in blue and red in Fig. 2(A). At time \( t=0\) the perturbation \( \varepsilon\mathbf{P}_i\) — for \( i=1,2\) — is applied to \( \mathbf{u}\) and the perturbed system evolves in time towards its new steady state

where \( \mathbf{v}_1^{(i)}(t)= \mathbf{G}(t) \mathbf{P}_i \mathbf{u}\) , and where \( \mathbf{G}(t)\) is defined in Eq. (10). These transient states are plotted for \( \mathbf{P}_1\) and \( \mathbf{P}_2\) , respectively, in Fig. 2(B) and (C). As \( t\) tends to infinity the asymptotic linear response converges to \( \mathbf{v}_{1}^{(i)}=\mathbf{G}\mathbf{P}_i \mathbf{u}\) , for \( i=1,2\) . In this particular case study, after three iterations of the transient linear response we observed convergence and the comparison to the invariant vector computed from Eq. (48) is shown as well in blue and red, in Fig. 2(B) and (C), respectively, were we observe that the linear response correction provides a faithful representation of the perturbed state.

Pan/Zoom

We now extend the study to continuous time stochastic differential equation. In this section we consider a paradigmatic set of equations representing equilibrium dynamics of a particle in a potential with two different competing attractors embedded in thermal (Gaussian) noise, namely, a double well potential.

The double well potential in two dimensions is a simple system displaying non-gaussian statistics and multimodality. We consider the following SDE on \( \mathbb{R}^2\) :

where \( \sigma>0\) and \( \mathbf{W}_t=(W_t^{(1)},W_t^{(2)})^{\top}\) , \( i=1,2\) , denote two dimensional independent standard Wiener processes. It is easy to see that the equation is driven by the potential \( -V(x,y) = x^2/2 - x^4/4 + \alpha x - y^2/2\) . The coefficient \( \alpha\) regulates the asymmetry between the two wells, so that for larger \( \alpha\) the deeper one of the wells is. Here, we select the coefficients of \( \alpha = -0.1\) and \( \sigma = 0.4\) .

In order to approximate the Fokker-Planck semigroup, we integrated Eq. (50) for \( 10^6\) time units with a timestep of \( 10^{-2}\) starting from the origin. This provides a time series \( {(x_k,y_k) }_{k=1}^{K}\) , where \( K = 10^8\) . The selected spatial domain is \( D =[-2,2]\times [-1.5,1.5]\) and it is subdivided into \( 2^{12}\) boxes \( B_i\) , where each axis is divided into \( 2^6\) segments. For a transition time of \( \tau = 25\cdot \mathrm{d} t\) , a Markov matrix \( \mathbf{M}\) is constructed following the two-dimensional analogue of Eq. (45). The resulting invariant measure is given by the eigenvector \( \mathbf{u}\) that solves Eq. (1), namely \( \mathbf{M} \mathbf{u} = \mathbf{u}\) , and it is plotted in Fig. 3(A).

Here we want to find the perturbation \( \mathbf{P}_1\) to the Markov matrix \( \mathbf{M}\) that maximizes the KL divergence. Following the strategy derived in Section 3.2 we have

A value of \( \varepsilon = 0.05\) is chosen so that the matrix \( \mathbf{M} + \varepsilon\mathbf{P}_1\) falls into the linear response regime with respect to \( \varepsilon\) . In Fig. 3(B) we show the perturbation \( \mathbf{P}_1\) applied to \( \mathbf{u}\) , and we observe how the perturbation acts mainly around the deepest well of the potential. The projected response on the \( x\) -axis is shown in Fig. 3(C). The transient values of the linear response are shown in grey scale— the black curve being the unperturbed measure— and the red curve the perturbed invariant measure. Indeed, as expected, the perturbed invariant measure that maximizes KL-divergence with respect to \( \mathbf{u}\) shifts all the probability mass onto the shallowest well.

Secondly, we consider the maximization of the expectation value of an observable as in [20] via (17). We consider the observable

This observable is a plateau-type of function that concentrates around the \( (0,0)\) origin in the phase-space precisely where the invariant measure has a local minimum half-way between the two wells. The idea for choosing this observable is to examine to which extent the leading order responose will push probability towards the origin in detriment of the two wells containing most of the probability. We let now \( \mathbf{P}_2\) denote the optimal pertubation for (17). One can guess that in order to maximize the expectation value of \( \Psi\) , the perturbation \( \mathbf{P}_2\) should transfer probability towards the origin. A value of \( \varepsilon = 0.1\) is chosen. The perturbed invariant measure \( \mathbf{v}^{(2)}\) — which solves Eq. (2)— is shown in Fig. 4(A), where the wells of the initially tilted potential are equilized. In Fig. 4(B) we show the perturbation \( \mathbf{P}_2\) applied to \( \mathbf{u}\) . While the big-picture features of Fig. 4(B) resemble those of Fig. 3(B), we must highlight that, as expected, the trasfer of probability is positive towards the origin of phase space in each of the wells, as opposed to Fig. 3(B), where the right-hand-side well experiences a shift to larger positive values of \( x\) .

The projected responses along the \( x\) -axis are shown in Fig. 4(C). The unperturbed and perturbed invariant measures are shown in black and red. Observe that the probability around \( x=0\) has increased. The grey curves represent the transient linear responses at different times indicated in the legends. At time \( t=10\) , the linear response is virtually equal to its asymptotic value, yet it does not match exactly the perturbed invariant vector in red. This suggests that the linear approximation does not fully capture the effects of the perturbation.

Pan/Zoom

Pan/Zoom

Adding non-conservative forces into the system drives it out of equilibrium yielding a nonvanishing net entropy production. In this section, we add a perturbation to the double well system (50) considered in the previous section with \( \alpha=0\) . The system results from substituting \( \mathbf{F} \to \mathbf{F} + \mathbf{R}\) , where \( \mathbf{R}\) is a rotational vector field such that the probability flux associated with it is zero, namely

It is classic to see (e.g. [27, Sec. 4.8]) that this forcing does not alter the invariant density \( \rho_0(x,y)\) of the unperturbed system. There are many ways to choose the rotation \( \mathbf{R}\) , we take it as:

for every \( (x,y)\) in \( \mathbb{R}^2\) and we represent it in Fig. 5(A) using a quiver-plot.

The resulting system \( \mathrm{d} \mathbf{x} = \left(\mathbf{F}\left( \mathbf{x}\right) + \mathbf{R}(\mathbf{x})\right)\, \mathrm{d} t + \mathrm{d}\mathbf{W}_t\) is integrated for \( 10^6\) time units with a timestep of \( 10^{-2}\) . This provides a time series \( {(x_k,y_k) }_{k=1}^{K}\) with \( K = 10^8\) data-points starting at the origin. As spatial domain we consider \( [-2,2]\times [-1.5,1.5]\) and it is subdivided in the same way as in the previous example. A transition Markov matrix is obtained in the same way for a transition time of \( \tau = 25\cdot \mathrm{d} t\) .

We expect that a perturbation that minimizes the entropy production will try to go against the rotational effect of the vector field \( \mathbf{R}(x,y)\) and we test this using the algorithm described in Section 3.3 and then reconstructing the drift via the argument in Section 4.1. The optimal perturbation \( \mathbf{P}_s\) is constructed with the explicit formula (30) derived. After that, we apply Eq. (43) to reconstruct the drift associated with \( \mathbf{P}_s\) . The vector field associated with \( \mathbf{P}_s\) is shown in Fig. 5(B) using a quiver-plot. We observe that there are two main rotational features around the modes of the distribution, just as in panel (A). Furthermore, we note that these rotational features gyrate clockwise, opposite, as expected, to the applied rotation \( \mathbf{R}\) . We observe that the drift associated to the optimal perturbation is not trivially reversing the rotation, but it preserves features of the deterministic field of the double well system.

It remains to check how \( \mathbf{P}_s\) differs from the perturbation matrix \( \mathbf{P}_r\) defined in Eq. (23) using the additive reversibilization of the Markov matrices. First of all, we show in Appendix A that the drift reconstruction algorithm of Eq. (43) applied to \( \mathbf{P}_r\) , indeed, provides the “reversed” rotational vector field \( \mathbf{R}\) . The difference is that if one imposes that the norm of the perturbation is small, say \( \varepsilon\) , the performance of \( \mathbf{P}_r\) is not as good as \( \mathbf{P}_s\) in reducing the entropy production. Indeed, this is shown in Fig. 5(C), where we plot the entropy productions \( s(\mathbf{M}+\varepsilon\mathbf{P}_{s})\) and \( s(\mathbf{M}+\varepsilon\mathbf{P}_{s})\) as a function of \( \varepsilon\) and where \( \mathbf{P}_r\) has been normalized to the unit Frobenius norm, and \( \mathbf{P}_s\) is normalized by construction. When \( \varepsilon \gtrsim 0.3\) , \( \mathbf{P}_r\) starts to outperform \( \mathbf{P}_s\) . This is expected since the reversibilized matrix obeys detailed balance by construction. Indeed, when the matrix \( \mathbf{P}_r\) is not normalized— \( \|\mathbf{P}_r\|_F\approx 4\) — we have \( s(\mathbf{M}+\mathbf{P}_{s})=0\) .

Pan/Zoom

The Lorenz 63 (L63) system is a low-dimensional truncation of the equations of atmospheric convection [28], whose dynamical properties have been extensively reported in the literature; we refer to the survey [50]. In its chaotic regime, while not displaying uniformly hyperbolic structure, it has been demonstrated on numerical grounds the existence of linear response, in the sense that small parameter forcings equate to proportional changes in the statistics [51], also under the perspective of Markov matrix expansions [15, 16]. The L63 system writes as follows:

where \( \mathbf{x} = (x,y,z)^\top\) and the model parameters are set to be \( s=10\) , \( r=28\) and \( b=8/3\) , for which the system displays chaotic behavior [52].

To estimate the semigroup \( { e^{t\mathcal{L}} }_{t\geq 0}\) , the system Eq. (54) is integrated for \( 10^5\) time units with a time step of \( \mathrm{d}t = 10^{-3}\) time units using a fourth order Runge-Kutta algorithm, after removing a transient of \( 10^4\) time units. The initial condition was taken to be \( x=y=z=1\) . A total of \( K = 10^8\) points in phase-space are recorded. The trajectory \( {(x_k,y_k,z_k)}_{k=0}^{K}\) shadows an attractor contained in a compact subset of \( \mathbb{R}^3\) [50], which lies in the domain \( \mathcal{D} = [- 20,20]\times [-30,30]\times [0,50]\) . The \( x\) -\( z\) domain, \( [-20,20]\times [0,50]\) , is then divided into \( N = 2^{12}\) equally sized, non-intersecting boxes \( {B_i}_{i=1}^N\) of dimensions \( 40/2^6 \times 50/2^6\) . Then, the transition time \( \tau=0.1\) is selected and the Markov matrix associted to (54) is computed using Eq. (47) as before.

As a first optimization example, we maximized the entropy functional using Eq. (17) to find a perturbation matrix \( \mathbf{P}_1\) . Secondly, we obtained the perturbation matrix \( \mathbf{P}_2\) such that the resulting invariant vector maximizes the KL-divergence with respect to the unperturbed vector \( \mathbf{u}\) ; see Eq. (20). Then, the perturbation size of \( \varepsilon = 0.05\) is chosen so that \( \mathbf{M} + \varepsilon \mathbf{P}_i\) is a (perturbed) Markov matrix with perturbed invariant measure \( \mathbf{v}^{(i)}\) satisfying:

for \( i=1,2\) , where we recall that the (asymptotic) linear response is given by \( \mathbf{v}_1^{(i)}=\mathbf{G}\mathbf{P}_i\mathbf{u}\) for \( i=1,2\) .

Figure 6 summarizes the results for the first perturbation matrix \( \mathbf{P}_1\) . Figure 6(A) shows the (perturbed) invariant measure of \( \mathbf{M}+\varepsilon\mathbf{P}_1\) . Compare to Figure 7(A), which contains the unperturbed invariant vector. Figure 6(B) contains the linear response of the system, revealing that the leading order correction of the invariant vector due to the perturbations is symmetric about the \( x=0\) axis. While the \( x\) -marginal experiences a contraction which seemingly reduces entropy, the \( z\) -axis marginal increases in breadth accounting for the full increase in entropy of the system due to the perturbation \( \mathbf{P}_1\) . This balance is better visualized in Fig. 6(C), where we computed the marginal distributions along the \( z\) -axis. The black and red lines show the unperturbed and perturbed invariant vectors, respectively. The grey line is the linear response correction, which in fact matches perfectly with the perturbed invariant vector, implying that we are within the linear response regime. We observe how, in order to amplify the entropy of the system, the perturbed distribution has to spread along the \( z\) -axis. This results align with the fact that in order to increase the entropy of a distribution, the latter has to resemble a uniform distribution as much as possible.

For the second perturbation problem considered, we show in Fig. 7(A) the invariant vector of the matrix \( \mathbf{M}\) . Figure 7(B) contains the linear response of the system, revealing this time that the leading order correction of the invariant vector due to the perturbations is asymmetric about the \( x=0\) axis. This asymetry is better seen in Fig. 7(C), where we computed the marginal distributions along the \( x\) -axis. The black and red lines show the unperturbed and perturbed invariant vectors, respectively. The grey line is the linear response correction, which coincides perfectly with the perturbed invariant vector, implying that we are within the linear response regime. This time, the perturbation seems to tilt the weight of the distribution towards one of the lobes of the attractor. In comparison with the results for the double well potential in Fig. 3(C), we observe an oppositve behaviour. Whereas for the double well potential the perturbation— which maximizes relative entropy— tends to equalize the competing wells, for the L63 system we observe that the perturbation that creates more information with respect to the unforced system, actually makes one of the lobes (or wings) of the attractor more probable. One of the possible explanations is that for the L63 system, while the lobes represent regimes of the system they are not competing attractors as it happens with the double well potential.

Pan/Zoom

Pan/Zoom

The role of invariant periodic solutions in shaping dynamical properties of chaotic and turbulent systems has been the subject of intense research [53, 54]. On the one hand, the subcritical transition to turbulence can be recast in terms of a sequence of bifurcations of a non trivial originating invariant solution of the system [55, 56]. On the other hand, the chaotic dynamical behaviour of turbulent systems can be described in terms of the underlying unstable periodic orbits which coexist with the chaotic attractor [57, 58, 59]. Of interest in these findings is the possibility of providing a low order, effective stochastic description of chaotic systems in terms of (a subset of) the underlying unstable periodic orbits [30, 29]. Here, as a last numerical example, we investigate the response properties of the reduced Markov Chain description of a turbulent flow provided in [29]. We consider a \( 3\) D Kolmogorov flow described by the following forced incompressible Navier Stokes equations

on a box domain \( \mathcal{D} =[0,L_x] \times [0,L_y] \times [0,L_z]\) with periodic boundary conditions. Here \( \mathbf{x}= (x,y,z)\in \mathcal{D}\) , \( \boldsymbol{v}= (v_x,v_y,v_z) =\boldsymbol{v}(\mathbf{x}) \) is the velocity field of the flow, \( p=p(\mathbf{x})\) is the pressure field and \( \nu\) is the kinematic viscosity. Furthermore, the flow is forced with \( \mathbf{F}(\mathbf{x})= \Omega \sin \left(\frac{2\pi y}{L_y} \right)\mathbf{e}_x\) where \( \mathbf{e}_x\) is the canonical vector in the \( x\) direction and \( \Omega\) is the amplitude of the perturbation. We set the parameters for the flow as \( \nu =0.05\) , \( \Omega=1\) , \( L_x=L_y=2\pi\) and \( L_z = \pi\) .

Pan/Zoom

Pan/Zoom

The above partial differential equation admits a laminar solution \( \boldsymbol{v}_L = (\Omega \nu^{-1} (L_y/2\pi)^2 \sin(2\pi y/L_y),0,0)\) which is linearly stable for all values of \( \nu\) . However, one commonly observes transient (but long lasting) turbulent behaviour if the system is triggered with finite amplitude perturbations. The authors in [29] have shown that it is possible to obtain a reduced order model of the high dimensional turbulent dynamic in terms of the underlying unstable periodic orbits of the system. In particular, they managed to identify \( 17\) unstable periodic orbits \( \gamma_i : \mathbb{R} \to \mathbb{R}^3\) such that \( \gamma_i(t+T_i) = \gamma_i(t)\) \( \forall t\) where \( T_i\) is the period of the \( i-th\) unstable periodic orbit. The shadowing properties of the turbulent trajectory, i.e., the amount of time that the turbulent trajectory spends in the neighborhood of an unstable orbit, have been assessed with topological data analysis techniques by comparing the difference between the shapes of the periodic orbits and the turbulent trajectory. The resulting Markov Chain description corresponds with the identification of the periodic orbits as fundamental states of the system and with modelling the turbulent dynamic as a stochastic hopping between them. The statistics of such hoppings are given in terms of a transition matrix \( \mathbf{M}=(M_{ij})\) describing the transition rates from the \( j-th\) periodic orbit to the \( i-th\) periodic orbit. The transition matrix is obtained by carefully studying the shadowing properties of the turbulent trajectory, see details in [29], and is reported here in Fig. 8.

The validity of this reduced order model is corroborated by the fact that statistics of the turbulent flow can be reconstructed from the Markov Chain description. In particular, given any observable of the system \( h(\mathbf{x})\) , such as the kinetic energy \( E(\mathbf{x}) = |\boldsymbol{v}(\mathbf{x})|^2 / 2\) or the dissipation \( D(\mathbf{x})= \nu |\nabla \times \boldsymbol{v}|^2\) , it is possible to associate to each periodic orbit an average value

This way we construct the coarse-grained observable vector \( \mathbf{h}=(\bar{h}_i)\) in \( \mathbb{R}^{17}\) . The invariant vector \( \mathbf{u}=(u_i)\) of the coarse grained system satisfies \( \mathbf{M}\mathbf{u}=\mathbf{u}\) , and the averages of observables of the turbulent flow can be reconstructed from the Markov matrix picture by considering suitably weighted expectation values

where the weights \( w_i = u_i T_i / (\sum_{j=1}^{17} u_j T_j)\) take into account the average residence time \( \tau_i = u_i T_i\) of the turbulent trajectory on each periodic orbit. We highlight that Eq. (58) is a coarse grained approximation of the expectation value of the actual observable \( h(\mathbf{x})\) evaluated on the turbuent trajectory.

We now consider a situation where a perturbation is applied to the system, such as the transitions between the unperturbed system’s UPOs are described by a perturbed transition matrix \( \mathbf{M}_\varepsilon = \mathbf{M} + \varepsilon \mathbf{P}\) where \( \varepsilon\) is a small parameter. In the long time limit, the perturbed system will reach its statistical steady state characterized by the invariant distribution \( \mathbf{u}+\varepsilon \mathbf{v_1}\) , with \( \mathbf{v}_1\) being given by (9). Expectation values in the perturbed system can be written in a linear approximation as \( \langle \mathbf{h} \rangle_\varepsilon = \langle \mathbf{h} \rangle + \varepsilon \sum_{i} \bar{h}_i \delta w_i\) , where the perturbations to the weights is given by

We seek for the perturbation matrix \( \mathbf{P}\) that maximizes the change in expectation value \( \sum_i \bar{h}_i \delta w_i\) . From the previous expression it is simple to see that the change in the weights can be written as \( \delta w_i = ( \mathbf{A} \mathbf{v}_1)_i\) , where \( \mathbf{A}\) is a suitable matrix, whose expression we omit here. We are then trying to optimize for the quantity \( \sum_i \bar{h}_i \delta w_i = \sum_i \bar{h}_i (\mathbf{A}\mathbf{v}_1)_i= \sum_j(\sum_i \bar{h}_i A_{ij}) v_{1,j} \mathrm{:=} \mathbf{f}^\top \mathbf{v}_1\) , where \( \mathbf{f} = \mathbf{A}^\top \mathbf{h}\) . This shows that this optimisation problem falls into the same class given by  (16), and we can then apply the same formula  (17).

Pan/Zoom

Pan/Zoom

We consider the optimal perturbation \( \mathbf{P}\) that maximizes the kinetic energy of the system. We remark that, in a linear response approach, the perturbation that instead minimizes the energy is simply \( - \mathbf{P}\) . Panel (a) in Fig. 11 shows the result of such perturbation. Here, we have checked that linearity, and thus the validity of (59), holds for \( \varepsilon \in (0,0.05]\) . The perturbation strongly penalizes the weight \( w_6\) which corresponds to a state with a high occupancy (\( u_6\) corresponds to the maximum of the invariant density) but with quite low energy compared to the other states. For a similar reason, state \( 11\) gets less populated than \( 10\) because, despite having \( w_{10} \approx w_{11}\) , the energy of the two states differs substantially. In general, the perturbation favours states with a longer period, which correspond to by and large high energy states. The perturbation acts in a non trivial way, as, for example, mass leaks from state \( 17\) , corresponding to the longest and quite energetic periodic orbit (in fact with the second largest weight), to favour transitions to other states in order to maximize the energy of the system.

We remark that, as opposed to the previous section where an Ulam approach has been taken, this analysis does not allow for an immediate reconstruction of the optimal perturbation field to be applied to (56). However, we provide evidence that the above perturbation of the Markov matrix corresponds to forcings that push the system towards a specific bundle of periodic orbits. In particular, we concentrate on the spatially dependent velocity profile \( U(y) = \frac{1}{L_x L_z}\int v_x(x,y,z) \mathrm{d}x\mathrm{d}z\) of solutions of (56) and investigate the effect of the perturbation on the average velocity profile obtained using (58). Panel (b) of Figure 11 (bottom) shows that the velocity profile, in a linear regime, is mostly left unaffected except for stretching the peak to higher values. This is in agreement with the velocity profiles of the underlying periodic orbits that are favoured by the perturbation. In fact, the periodic orbits with \( \delta w_i > 0\) correspond to velocity profiles that are higher than the average one. There is one exception to the above: the periodic orbit \( i=15\) has a velocity profile below the average one. We remark that the weight of this state is very small due to a very low \( u_{15} \approx 0.005\) invariant density.

  The authors thank Gökhan Yalnız for sharing the data regarding the unstable periodic orbits for the Kolmogorov flow. MSG acknowledges the partial support provided by the Horizon Europe project Past2Future (Grant No. 101184070). N.Z. has been supported by the Wallenberg Initiative on Networks and Quantum Information (WINQ). GC is a member of GNAMPA of the Italian Istituto Nazionale di Alta Matematica (INdAM).

In the main text we provided the solution to the Lagrange multiplier equation for the optimization of Eq. (60). We call this Method 1. In this appendix, we also show that the solutions to this optimization problem can be viewed as the orthogonal projection of the data matrix onto the kernel of a suitable functional— described in details in the previous section—. We shall call this Method 2. In principle, both methods should provide the same results, but the radically different implementation and conceptual basis is prone to numerical disagreements. Here we show that both methods give the same results. Method 1 and Method 2 are applied to the entropy production minimization problem of Eq. (24). This, in addition, confirms that the methods indeed minimize the entropy production.

We randomly initialized \( K = 10^4\) matrices of size \( 50\times 50\) with a diagonal dominance of \( 100\%\) — the diagonal element is the sum of the non-diagonal entries in its column— and a sparsity of \( 50\%\) . The entries are randomly sampled from a uniform distribution in \( [0,1]\) and normalized to obtain Markov matrices. Non-mixing matrices— for instance, if a diagonal element was equal to unity— are discarded and randomly replacing them so that the ensemble consists of \( K=10^4\) members. The factors of diagonal dominance and sparsity control the spectral gap— defined as \( 1-|\lambda_1|\) , where \( \lambda_1\) is the subdominant eigenvalue—. We chose the values so that they resemble the sparsity and diagonal dominance of the matrices considered earlier in this study. A mean spectral gap of \( 0.31\pm 0.03\) was obtained.

Then Method 1 and Method 2 are applied to the matrices to obtain perturbation matrices \( \mathbf{P}_i^{(k)}\) , for \( k=1,…,K\) and \( i=1,2\) — the subscript \( i\) refers to the method used—. Then, for a value of \( \varepsilon = 10^{-3}\) , we computed the entropy production \( s\left( \mathbf{M}^{(k)} + \varepsilon\mathbf{P}_{i}^{(k)}\right)\) and computed the relative difference with respect to \( s\left( \mathbf{M}^{(k)} \right)\) for each random initialization and method. Then we computed the (normalized) histogram shown in Fig. 15(A). Indeed, this shows that both methods actually minimize the entropy production consistently. Secondly, we computed the histogram of relative difference between methods for each random initialization of the matrix \( \mathbf{M}^{(k)}\) . The normalized histogram in Fig. 15(B) shows that both methods provide different results, but that their difference is small and here we attribute it to numerical issues.

Pan/Zoom