Inferring entropy production in many-body systems using nonequilibrium maximum entropy

TL;DR

The paper tackles the challenge of inferring entropy production (EP) in high-dimensional nonequilibrium stochastic systems where reconstructing full trajectory probabilities is intractable. It introduces a nonequilibrium Maximum Entropy (MaxEnt) variational principle that yields a convex dual optimization problem for a lower bound on EP, expressed as ${\5ig\} = \max_{\bm{\theta}} \left( {\bm{\theta}}^\top \langle \bm{g} \rangle_p - \ln \langle e^{ {\bm{\theta}}^\top \bm{g}} \rangle_{\tilde{p}} \right)$, with the true EP given by ${\Sigma} = D(p \Vert \tilde{p})$. This framework provides a trajectory-level EP estimator ${\sigma_{\bm{\theta}^*}}(\bm{x}) = {\bm{\theta}^{*}}^\top \bm{g}(\bm{x}) - \ln \langle e^{ {\bm{\theta}^{*}}^\top \bm{g}} \rangle_{\tilde{p}}$, and yields ${\Sigma}_{\bm{g}} \le {\Sigma}$ with a Pythagorean decomposition ${\Sigma} = {\Sigma}_{\bm{g}} + {\Sigma}^{\perp}_{\bm{g}}$, where ${\Sigma}^{\perp}_{\bm{g}} = D(p \Vert q^*)$. The method further supports a hierarchical decomposition by interaction order and extends to multipartite observables, enabling independent subproblems and improved scalability. Importantly, only samples of trajectory observables (e.g., spatiotemporal correlations) are required, avoiding reconstruction of high-dimensional distributions, and the approach connects to thermodynamic uncertainty relations. Demonstrations on a disordered nonequilibrium spin model with 1000 spins and on Neuropixels spike-train data show accurate EP bounds and interpretable coupling-inference, with a public code release to facilitate application.

Abstract

We propose a method for inferring entropy production (EP) in high-dimensional stochastic systems, including many-body systems and non-Markovian systems with long memory. Standard techniques for estimating EP become intractable in such systems due to computational and statistical limitations. We infer trajectory-level EP and lower bounds on average EP by exploiting a nonequilibrium analogue of the Maximum Entropy principle, along with convex duality. Our approach uses only samples of trajectory observables, such as spatiotemporal correlations. It does not require reconstruction of high-dimensional probability distributions or rate matrices, nor impose any special assumptions such as discrete states or multipartite dynamics. In addition, it may be used to compute a hierarchical decomposition of EP, reflecting contributions from different interaction orders, and it has an intuitive physical interpretation as a "thermodynamic uncertainty relation." We demonstrate its numerical performance on a disordered nonequilibrium spin model with 1000 spins and a large neural spike-train dataset.

Figures (4)

• Figure 1: Disordered nonequilibrium spin model with 1000 spins. (a) Steady-state EP estimates for different inverse temperatures $\beta$. (b) Asymmetry of inferred parameters, shown against the true coupling asymmetries in the model for $\beta=2.5$ ($R^2=0.9831$). Estimates are based on $10^9$ state transitions sampled by Monte Carlo.
• Figure 2: (a) Estimated EP per expected number of spikes $R$, in the Neuropixels Visual Behavior dataset for three conditions. Error bars indicate standard error of the mean. (b) A sample of inferred coupling coefficients $\theta^*_{ij}$ grouped by visual area. Here we select 200 neurons with highest firing rate from an active trial. To improve visualization, lower triangle shows $\theta_{ji}^*\equiv-\theta_{ij}^*$ for $i<j$.
• Figure 3: Computation time for ${\Sigma}_{\bm{g}}$ using the regular (non-multipartite) optimization \ref{['eq:dual']} versus the multipartite decomposition \ref{['eq:lb']}. (a) Runtime versus number of spins at fixed number of samples $2\times 10^4$. (b) Runtime versus number of samples per spin (system with $40$ spins). Shaded bands: standard deviations over 1000 trials. The end of the "regular" curve indicates the point where the GPU runs out of memory. Hardware: Intel Core i9-12900KF CPU, NVIDIA GeForce RTX 3050 GPU with 8 GB VRAM.
• Figure 4: Comparison between our EP bound ${\Sigma}_{\bm{g}}$\ref{['eq:dual']} and the variational bound $\Sigma_{\bm{g}}^{\text{KO}}$\ref{['eq:fdivdual']} on a 3-state unicyclic system. (a) For the non-antisymmetric observable $g({{\bm{x}}})=1 + \delta_{x_{0}+1,x_{1}}-\delta_{x_{0}-1,x_{1}}$. (b) For the antisymmetric observable $g({{\bm{x}}})=\delta_{x_{0}+1,x_{1}}-\delta_{x_{0}-1,x_{1}} - 1.2 (\delta_{x_{1},1} \delta_{x_{0},0}-\delta_{x_{0},1} \delta_{x_{1},0})$. In both cases, ${\Sigma}_{\bm{g}}$ diverges in the irreversible limit $\lambda\to \infty$, while $\Sigma_{\bm{g}}^{\text{KO}}$ saturates at a finite value.