$\alpha$-divergence improves the entropy production estimation via machine learning

How irreversible does a process look? One may pose this question for two distinct reasons. First, whether a biological process requires energy dissipation is often a subject of much debate [1, 2]. To resolve this issue, it is useful to note that irreversibility suggests energy dissipation. Various hallmarks of irreversibility, such as the breaking of the fluctuation-dissipation theorem [3] and the presence of nonequilibrium probability currents in the phase space [4, 5], have been used to determine whether energy is dissipated. Second, whether a nonequilibrium system allows for an effective equilibrium description is an important issue. For instance, in active matter, despite the energy dissipation at the microscopic level, it has been argued that the large-scale phenomena allow for an effective equilibrium description [6, 7, 8, 9, 10]. If we can quantify the irreversibility of an empirical process at various levels of coarse-graining [11, 12], it will provide us with helpful clues as to whether we should look for an effective equilibrium theory for the process.

Based on the framework of stochastic thermodynamics, modern thermodynamics assigns entropy production (EP) to each stochastic trajectory based on its irreversibility [13]. Thus, empirically measuring the irreversibility of a process is closely tied to the problem of estimating EP from sampled trajectories [14, 15, 16, 17, 18, 19, 20, 21]. A straightforward approach to the problem is to evaluate the relevant transition probabilities by directly counting the number of trajectory segments, which is called the plug-in method [14, 15]. The method, readily applicable to discrete systems, can also be applied to continuous systems through the use of kernel functions [16]. However, while this method is simple and intuitive, it requires a huge ensemble of lengthy trajectories for accurate estimations (curse of dimensionality). More recent studies proposed methods based on universal lower bounds of the average EP, such as the thermodynamic uncertainty relations [16, 17, 18, 19] and the entropic bound [20]. While these methods do not suffer from the curse of dimensionality and are applicable even to non-stationary processes [19, 20], their accuracy is impaired when the underlying bounds are not tight. Moreover, these methods are applicable only to the estimation of the average EP, not the EP of each trajectory.

Meanwhile, with the advent of machine learning techniques in physics, a novel method for EP estimation using artificial neural networks has been developed [21]. This method, called the Neural Estimator for Entropy Production (NEEP), minimizes the loss function based on a variational representation of the Kullback-Leibler (KL) divergence. Without any presupposed discretization of the phase space and using the rich expressivity of neural networks, the NEEP suffers far less from the complications of the sampling issues and is applicable to a diverse range of stochastic processes [19].

Still, the NEEP has its limits. Its accuracy deteriorates when the nonequilibrium driving is strong or when the dynamics slows down so that the phase space is poorly sampled. In this study, we show that the NEEP can be significantly improved by changing the loss function. Toward this purpose, we propose the $\alpha$-NEEP, which generalizes the NEEP. Instead of the KL divergence, the $\alpha$-NEEP utilizes the $\alpha$-divergence, which has been mainly used in the machine learning community [22, 23, 24, 25]. We demonstrate that the $\alpha$-NEEP with nonzero values of $\alpha$ shows much more robust performance for a broader range of nonequilibrium driving and sampling quality, with $\alpha=-0.5$ showing the optimal performance overall. This is corroborated by an analytically tractable simplification of the $\alpha$-NEEP that shows the optimality of $\alpha=-0.5$.

The rest of this paper is organized as follows. After reviewing the original NEEP and its limitations (Sec. II), we introduce the $\alpha$-NEEP (Sec. III) and demonstrate its enhanced performance for three different examples of nonequilibrium systems (Sec. IV). Then we investigate the rationale behind the observed results using a simplified model describing how the $\alpha$-NEEP works (Sec. V). Finally, we sum up the results and discuss their implications (Sec. VI).

We first give a brief overview of how the original NEEP [21] estimates EP at the trajectory level. Suppose our goal is to estimate EP of a Markov process in discretized time, $\mathbf{x}_{t}$, in a $d$-dimensional space. For every ordered pair of states, denoted by $\mathbf{x}\equiv(\mathbf{x}_{t},\mathbf{x}_{t+1})$, there is EP associated with the transition between them, which is given by the ratio between the forward and the backward path probabilities

$\displaystyle\Delta S(\mathbf{x})=\log\frac{p(\mathbf{x})}{p(\tilde{\mathbf{x}})},$

(1)

where $\tilde{\mathbf{x}}\equiv(\mathbf{x}_{t+1},\mathbf{x}_{t})$. Note that, throughout this study, we use the unit system in which the Boltzmann constant can be set to unity ($k_{\mathrm{B}}=1$). Then it follows that the ensemble average of this EP is equivalent to the KL divergence, which satisfies the inequality

	$\displaystyle\langle\Delta S(\mathbf{x})\rangle$	$\displaystyle=D_{\rm{KL}}[p(\mathbf{x}),p(\tilde{\mathbf{x}})]$
		$\displaystyle\geq\left\langle\log r(\mathbf{x})-\frac{p(\tilde{\mathbf{x}})}{p(\mathbf{x})}r(\mathbf{x})+1\right\rangle_{p(\mathbf{x})}$		(2)

for any positive function $r(\mathbf{x})$, given that $\langle\cdot\rangle_{p(\mathbf{x})}$ denotes the average with respect to the distribution $p(\mathbf{x})$. This inequality can be proven as follows: since $\log$ is a concave function, the line tangent to any point never falls below the function. Thus, $\frac{1}{r_{0}}(r-r_{0})+\log r_{0}\geq\log r$ for any $r$ and $r_{0}$. By putting $r_{0}=r_{0}(\mathbf{x})=p(\mathbf{x})/p(\tilde{\mathbf{x}})$ and taking the average with respect to $p(\mathbf{x})$, we get the inequality. In this derivation, we immediately note that the equality condition is satisfied if and only if $r(\mathbf{x})=r_{0}(\mathbf{x})=p(\mathbf{x})/p(\tilde{\mathbf{x}})$. Hence, by varying $r(\mathbf{x})$ to maximize the right-hand side of Eq. (II), we accurately estimate the average EP $\langle\Delta S(\mathbf{x})\rangle$. For this reason, Eq. (II) is called the variational representation of the KL divergence. Moreover, as a byproduct, we also obtain the function $r_{0}(\mathbf{x})$, which yields an accurate estimate for trajectory-level EP by $\Delta S(\mathbf{x})=\log r_{0}(\mathbf{x})$.

Kim et al. [21] used these properties to construct the loss function of the NEEP. More specifically, they introduce $s_{\theta}(\mathbf{x})$, an estimator for trajectory-level EP parametrized by $\theta$, and put $r(\mathbf{x})=e^{s_{\theta}(\mathbf{x})}$. Then, Eq. (II) can be rewritten as

$\displaystyle\langle\Delta S(\mathbf{x})\rangle\geq\left\langle s_{\theta}(\mathbf{x})-e^{s_{\theta}(\tilde{\mathbf{x}})}+1\right\rangle_{p(\mathbf{x})},$

(3)

where $\langle r(\mathbf{x})\,p(\tilde{\mathbf{x}})/p(\mathbf{x})\rangle_{p(\mathbf{x})}=\langle r(\mathbf{x})\rangle_{p(\tilde{\mathbf{x}})}=\langle r(\tilde{\mathbf{x}})\rangle_{p(\mathbf{x})}$ has been used based on the one-to-one correspondence between $\mathbf{x}$ and $\tilde{\mathbf{x}}$. Furthermore, since EP is odd under time reversal, i.e., $\Delta S(\mathbf{x})=-\Delta S(\tilde{\mathbf{x}})$, it is natural to impose the same condition on $s_{\theta}$. This leads to the inequality

$\displaystyle\langle\Delta S(\mathbf{x})\rangle\geq\left\langle s_{\theta}(\mathbf{x})-e^{-s_{\theta}(\mathbf{x})}+1\right\rangle_{p(\mathbf{x})},$

(4)

which motivates the loss function

$\displaystyle\mathcal{L}(\theta)=\left\langle-s_{\theta}(\mathbf{x})+e^{-s_{\theta}(\mathbf{x})}-1\right\rangle_{p(\mathbf{x})},$

(5)

so that the minimization of $\mathcal{L}(\theta)$ ensures the accurate EP estimation $s_{\theta}(\mathbf{x})=\Delta S(\mathbf{x})$.

It is notable that $\mathcal{L}(\theta)$ defined above is a convex functional of $s_{\theta}$. Thus, as long as the $\theta$-dependence of $s_{\theta}$ is well behaved, any gradient-descent algorithm can reach the global minimum of $\mathcal{L}(\theta)$ without getting trapped in a local minimum. In this regard, the rugged loss function landscape is not a major issue of the NEEP.

However, the performance of the NEEP strongly depends on how well $p(\mathbf{x})$ is sampled. Since the second term of $\mathcal{L}(\theta)$ depends exponentially on $s_{\theta}(\mathbf{x})$, rare transitions with minute $p(\mathbf{x})$ can make nonnegligible contributions to $\mathcal{L}(\theta)$ when $e^{-s_{\theta}(\mathbf{x})}$ is extremely large. Since the frequency of rare events is subject to considerable sampling noise, the performance of the original NEEP deteriorates in the presence of a strong nonequilibrium driving which induces rare transitions with large negative EP. In the following section, we propose a loss function that remedies this weakness of the NEEP.

Here we formulate a generalization of the NEEP loss function with the goal of mitigating its strong sampling-noise dependence. We note that the loss function needs not be an estimator of average EP $\langle\Delta S(\mathbf{x})\rangle$, for our goal is to estimate $\Delta S(\mathbf{x})$ at the level of each trajectory. Thus, while the original NEEP uses the variational representation of the KL divergence corresponding to $\langle\Delta S(\mathbf{x})\rangle$, we propose a different approach based on the variational representation of the $\alpha$-divergence, which quantifies the difference between a pair of probability distributions $p(\mathbf{x})$ and $q(\mathbf{x})$ as

$\displaystyle D_{\alpha}[p\!:\!q]\equiv\left\langle\frac{[p(\mathbf{x})/q(\mathbf{x})]^{\alpha}-1}{\alpha(1+\alpha)}\right\rangle_{p(\mathbf{x})}.$

(6)

Since this reduces to the KL divergence in the limit $\alpha\to 0$, our approach generalizes the NEEP by introducing an extra parameter $\alpha$. To emphasize this aspect, we term our method the $\alpha$-NEEP.

The goal of the $\alpha$-NEEP is to find $r(\mathbf{x})$ that minimizes the loss function

$\displaystyle\mathcal{L}_{\alpha}[r]\equiv\left\langle-\frac{r(\mathbf{x})^{\alpha}-1}{\alpha}+\frac{q(\mathbf{x})}{p(\mathbf{x})}\frac{r(\mathbf{x})^{1+\alpha}-1}{1+\alpha}\right\rangle_{p(\mathbf{x})},$

(7)

where $p(\mathbf{x})$ and $q(\mathbf{x})$ are probability density functions, and $\alpha$ is a real number other than $0$ and $-1$. See Appendix B for discussions of these two exceptional cases. It can be rigorously shown (see Appendix B) that $\mathcal{L}_{\alpha}[r]$ satisfies the inequality

$\displaystyle\mathcal{L}_{\alpha}[r]\geq-D_{\alpha}[p\!:\!q]$

(8)

where the equality is achieved if and only if $r(\mathbf{x})=p(\mathbf{x})/q(\mathbf{x})$ for all $\mathbf{x}$. In other words, by minimizing $\mathcal{L}_{\alpha}[r]$ to find $D_{\alpha}[p\!:\!q]$, we also obtain an estimate for the ratio $p(\mathbf{x})/q(\mathbf{x})$. We note that the properties of $\mathcal{L}_{\alpha}[r]$ used here are also valid for a much more general class of loss functions, as discussed in [22, 23] (also see Appendix B).

Based on Eq. (8), we can construct a loss function

$\displaystyle\mathcal{L}_{\alpha}(\theta)=\left\langle-\frac{e^{\alpha s_{\theta}(\mathbf{x})}-1}{\alpha}+\frac{e^{-(1+\alpha)s_{\theta}(\mathbf{x})}-1}{1+\alpha}\right\rangle.$

(9)

Note that this reduces to the loss function of the original NEEP shown in Eq. (5) in the limit $\alpha\to 0$. If $s_{\theta}$ is sufficiently well behaved, the minimization of $\mathcal{L}_{\alpha}(\theta)$ yields the minimizer $\theta^{*}$ which satisfies $\mathcal{L}_{\alpha}(\theta^{*})=-D_{\alpha}[p(\mathbf{x}),p(\tilde{\mathbf{x}})]$ and $\Delta S(\mathbf{x})=s_{\theta^{*}}(\mathbf{x})$. The former is generally not equal to average EP $\langle\Delta S(\mathbf{x})\rangle$ (unless $\alpha\to 0$), but the latter ensures the accurate estimation of trajectory-level EP $\Delta S(\mathbf{x})$.

Comparing Eqs. (5) and (9), one readily observes that the exponential dependence on $s_{\theta}(\mathbf{x})$ can be made much weaker in $\mathcal{L}_{\alpha}(\theta)$ by choosing the value of $\alpha$ between $-1$ and $0$. Since this mitigates the detrimental effects of the sampling error associated with rare trajectories with large negative $s_{\theta}(\mathbf{x})$, one can naturally expect that the performance of the $\alpha$-NEEP is much more robust against strong nonequilibrium driving. This is confirmed in the following sections.

Before proceeding, a few remarks are in order:

1. 1.

   The loss function satisfies $\mathcal{L}_{-(1+\alpha)}(\theta)=\mathcal{L}_{\alpha}(\theta)$, so the $\alpha$-NEEP is symmetric under the exchange $\alpha\leftrightarrow-(1+\alpha)$. For this reason, in the rest of this paper, we focus on the regime $-0.5\leq\alpha\leq 0$ (the regime $\alpha>0$ leads to very poor performance and is left out).

2. 2.

   From the antisymmetry $\Delta S(\mathbf{x})=-\Delta S(\tilde{\mathbf{x}})$, we may set the estimator $s_{\theta}$ to be related to the feedforward neural network (FNN) output $h_{\theta}$ as

   $\displaystyle s_{\theta}(\mathbf{x})=h_{\theta}(\mathbf{x})-h_{\theta}(\tilde{\mathbf{x}})\;,$

   (10)

   so that the neural network focuses on the estimators that satisfy the antisymmetry of EP for more efficient training. The method described so far is schematically illustrated in Fig. 1.

3. 3.

   We emphasize that the minimized $\mathcal{L}_{\alpha}$ is not directly related to average EP. In all cases, we compute the average EP by averaging $s_{\theta}$ over the sampled transitions.

To assess the performance of the $\alpha$-NEEP for various values of $\alpha$, we apply the method to toy models of nonequilibrium systems, namely the two-bead model, the Brownian gyrator, and the driven Brownian particle.

(i) The two-bead model. This model has been used in a number of previous studies as a benchmark for testing EP estimators [4, 16, 18, 21]. The model consists of two one-dimensional (1D) overdamped beads which are connected to each other and to the walls on both sides by identical springs, see Fig. 2(a). The beads are in contact with heat baths at temperatures $T_{\rm{h}}$ and $T_{\rm{c}}$ with $T_{\rm{h}}>T_{\rm{c}}$. Denoting by $x_{\rm{h}}$ ($x_{\rm{c}}$) the bead in contact with the hot (cold) bath, the stochastic equations of motion are given by

	$\displaystyle\gamma\dot{x}_{\rm{h}}$	$\displaystyle=k(-2x_{\rm{h}}+x_{\rm{c}})+\sqrt{2\gamma T_{\rm{h}}}\xi_{\rm{h}}(t)\;,$		(11a)
	$\displaystyle\gamma\dot{x}_{\rm{c}}$	$\displaystyle=k(-2x_{\rm{c}}+x_{\rm{h}})+\sqrt{2\gamma T_{\rm{c}}}\xi_{\rm{c}}(t)\;.$		(11b)

Here $k$ is the spring constant, $\gamma$ the friction coefficient, and $\xi_{\rm{h},\,\rm{c}}$ the Gaussian thermal noise with zero means and $\langle\xi_{\rm{h}}(t)\xi_{\rm{h}}(t^{\prime})\rangle=\langle\xi_{\rm{c}}(t)\xi_{\rm{c}}(t^{\prime})\rangle=\delta(t-t^{\prime})$. For infinitesimal displacements $(dx_{\rm{h}},dx_{\rm{c}})$, the associated EP is given by

$\displaystyle\Delta S=\frac{k}{\gamma}\left[\frac{2x_{\rm{h}}-x_{\rm{c}}}{T_{\rm{h}}}\circ dx_{\rm{h}}+\frac{2x_{\rm{c}}-x_{\rm{h}}}{T_{\rm{c}}}\circ dx_{\rm{c}}\right]+\Delta S_{\text{sys}}\;,$

(12)

where $\circ$ denotes the Stratonovich product and $\Delta S_{\text{sys}}$ the change of the system’s Shannon entropy, namely

$\displaystyle\Delta S_{\text{sys}}=-\ln\frac{p_{\mathrm{s}}(x_{\rm{h}}+dx_{\rm{h}},x_{\rm{c}}+dx_{\rm{c}})}{p_{\mathrm{s}}(x_{\rm{h}},x_{\rm{c}})}$

(13)

for the steady-state distribution $p_{\mathrm{s}}(x_{\rm{h}},x_{\rm{c}})$. Since the system is fully linear, $p_{s}(x_{\rm{h}},x_{\rm{c}})$ can be calculated analytically. Thus the EP of this model can be calculated exactly using Eq. (12) and compared with the $\alpha$-NEEP result.

To see how the predicted EP differs from the true EP, we observe the behavior of the mean square error (MSE) $\langle(s_{\theta}-\Delta S)^{2}\rangle$. In Fig. 2(b), we observe that strengthening the nonequilibrium driving (by increasing $T_{h}$ while keeping $T_{c}=1$) tends to impair the EP estimation. This is because a stronger driving makes the reverse trajectories of typical trajectories rarer, lowering the sample quality. The adverse effects of the nonequilibrium driving are the strongest for the original NEEP ($\alpha=0$), which are mitigated by choosing different values of $\alpha$. Remarkably, choosing $\alpha=-0.5$ leads to the most robust performance against the driving.

As an alternative measure of the estimator’s performance, we also observe the ratio between the predicted average EP $\sigma_{\text{pred}}$ and the exact average EP $\sigma$. The results are shown in Fig. 2(c), which exhibit two different regimes. As $T_{h}$ increases, there is a regime where the estimator overestimates average EP, which is followed by an underestimation regime. A detailed explanation for this behavior will be given in Sec. V using a simplified model. At the moment, we note that $\sigma_{\text{pred}}/\sigma$ tends to deviate away from $1$ most strongly for the original NEEP ($\alpha=0$), while choosing different values of $\alpha$ makes the ratio stay closer to $1$. Again, the optimal value of $\alpha$ seems to be $-0.5$.

(ii) The Brownian gyrator. This simple model of a single-particle heat engine allows us to check the effects of a nonequilibrium driving apart from the temperature difference $T_{\mathrm{h}}>T_{\mathrm{c}}$. The dynamics of the model is governed by

	$\displaystyle\gamma\dot{x}_{\rm{h}}$	$\displaystyle=-\partial_{x_{\rm{h}}}U(x_{\rm{h}},x_{\rm{c}})+\varepsilon x_{\rm{c}}+\sqrt{2\gamma T_{\rm{h}}}\,\xi_{\rm{h}}(t)\;,$		(14a)
	$\displaystyle\gamma\dot{x}_{\rm{c}}$	$\displaystyle=-\partial_{x_{\rm{c}}}U(x_{\rm{h}},x_{\rm{c}})+\delta x_{\rm{h}}+\sqrt{2\gamma T_{\rm{c}}}\,\xi_{\rm{c}}(t)\;,$		(14b)

where $U(x_{\rm{h}},x_{\rm{c}})=\frac{1}{2}k(x_{\rm{h}}^{2}+x_{\rm{c}}^{2})$ is the harmonic potential, and $(\varepsilon x_{\mathrm{c}},\delta x_{\mathrm{h}})$ is a nonconservative force that drives the system out of equilibrium and enables work extraction. See Fig. 3(a) for an illustration of this system. For infinitesimal displacements $(dx_{\rm{h}},dx_{\rm{c}})$, the associated EP is given by

$\displaystyle\Delta S=-\frac{Q_{\rm{h}}}{T_{\rm{h}}}-\frac{Q_{\rm{c}}}{T_{\rm{c}}}+\Delta S_{\text{sys}},$

(15)

where

	$\displaystyle Q_{\rm{h}}$	$\displaystyle=(\partial_{x_{\rm{h}}}U-\epsilon x_{c})\circ dx_{\rm{h}},$		(16a)
	$\displaystyle Q_{\rm{c}}$	$\displaystyle=(\partial_{x_{\rm{c}}}U-\delta x_{h})\circ dx_{\rm{c}},$		(16b)

and $\Delta S_{\text{sys}}$ the change of the system entropy. Again, the system is fully linear and the steady-state distribution can be calculated analytically, allowing exact calculations of EP at the trajectory level.

Setting $T_{\mathrm{h}}/T_{\mathrm{c}}=10$ and $\varepsilon=-\delta$, we vary the magnitude of $\varepsilon$ to assess the robustness of the $\alpha$-NEEP in terms of the MSE and the ratio $\sigma_{\text{pred}}/\sigma$, as shown in Figs. 3(b) and (c), respectively. The results are qualitatively similar to the case of the two-bead model: as the nonconservative driving gets stronger, the performance of the original NEEP ($\alpha=0$) deteriorates the most, while other values of $\alpha$ yield more robust results. Again, $\alpha=-0.5$ seems to be the optimal choice.

(iii) The driven Brownian particle. While the two examples given above were both linear systems, we also consider a nonlinear system featuring a 1D overdamped Brownian particle in a periodic potential $U(x)=A\sin x$ driven by a constant force $f$. The motion of the particle is described by the Langevin equation

$\displaystyle\gamma\dot{x}=f-U^{\prime}(x)+\sqrt{2\gamma T}\xi(t)\;,$

(17)

where $\xi(t)$ is a Gaussian white noise with unit variance. See Fig. 4(a) for an illustration of the model. For sufficiently large $A$, this model can approximate the behaviors of the Markov jump process on a discrete chain. For this model, the EP associated with the infinitesimal displacement $dx$ is given by

$\displaystyle\Delta S=\frac{-f\,dx+U(x+dx)-U(x)}{T}+\Delta S_{\text{sys}}\;,$

(18)

where $\Delta S_{\text{sys}}=-\ln p_{\mathrm{s}}(x+dx)/p_{\mathrm{s}}(x)$ again denotes the Shannon entropy change for the steady-state distribution $p_{\mathrm{s}}(x)$. Since the system is 1D, it is straightforward to obtain $p_{\mathrm{s}}(x)$ by numerical integration. Thus, the EP of this model can also be calculated exactly and compared to the $\alpha$-NEEP result.

Fixing $f=32$, the performance of the $\alpha$-NEEP for this model is shown in Figs. 4(b) and (c) in terms of the MSE and the ratio $\sigma_{\mathrm{pred}}/\sigma$, respectively. Due to the presence of a strong background driving ($f=32$), there are already considerable differences among different methods at $A=0$. But it is worth noting that increasing the amplitude $A$ of the periodic potential $U(x)$ clearly increases the MSE and makes $\sigma_{\mathrm{pred}}/\sigma$ deviate farther away from $1$ for the original NEEP ($\alpha=0$). This may be the consequence of rarer movements (jumps from one potential well to the next) across the system as the potential well gets deeper, which means rare trajectories are even more poorly sampled. The $\alpha$-NEEPs with nonzero values of $\alpha$ are much more robust against the increase of $A$, with $\alpha=-0.5$ showing the best performance overall.

The results shown thus far clearly indicate that, by choosing a nonzero value of $\alpha$, the $\alpha$-NEEP can exhibit a much more robust performance against the adverse effects of the nonequilibrium driving. Moreover, $\alpha=-0.5$ seems to exhibit the best performance in many cases. To gain more intuition into these results, we simplify the EP estimation problem to the density-ratio estimation problem for a 1D random variable. To be specific, we estimate the log ratio $s(x)\equiv\log p(x)/p(-x)$ given samples drawn from the distribution $p(x)$. It is intuitively clear that this problem is structurally equivalent to EP estimation.

For further simplification, we set

$\displaystyle p(x)=\begin{cases}\mathcal{N}\exp\left[-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right]&\text{for $|x-\mu|\leq k\sigma$,}\\ 0&\text{otherwise.}\end{cases}$

(19)

Here $\mathcal{N}$ is a suitable normalization factor, $\mu$ the positive mean of the distribution, $\sigma$ the width of the distribution, and $k$ a positive number truncating the tails of the distribution. While $k=\infty$ corresponds to the perfect sampling of a Gaussian distribution, a finite $k$ corresponds to the case where the tails of the distribution are poorly sampled.

For $k=\infty$, the correct answer to the problem is a linear function $s(x)=\theta_{0}x$, where $\theta_{0}\equiv 2\mu/\sigma^{2}$. Thus, for further simplicity, we focus on the one-parameter model $s_{\theta}(x)=\theta x$, which estimates $s(x)$ using only a single parameter $\theta$. For this problem, the suitable loss function is obtained as an analog of Eq. (9):

$\displaystyle\mathcal{L}_{\alpha}(\theta)=\left\langle-\frac{(1+\alpha)e^{\alpha\theta x}-1}{\alpha}+e^{-(1+\alpha)\theta x}\right\rangle_{p(x)}\;.$

(20)

If $k$ is large but finite, the minimum of this loss function shifts to $\theta_{0}+\Delta\theta$, where $\Delta\theta$ can be expanded to the leading orders in $1/k$:

$\displaystyle\Delta\theta\sim\exp\left[-\frac{k^{2}}{2}+\frac{2\mu}{\sigma}\left(\left|\alpha+\frac{1}{2}\right|+\frac{1}{2}\right)k\right]\;.$

(21)

This clearly shows that $\alpha=-0.5$ gives the least shift $\Delta\theta$, as also illustrated by various results shown in Fig. 5.

In Fig. 5(a), we show that the shift of the minimum $\Delta\theta$ tends to increase as the tail sampling becomes poorer (i.e., $k$ decreases). The landscapes of the loss function $\mathcal{L}_{\alpha}(\theta)$, shown in the inset of Fig. 5(a), also confirm this observation. The increase of the error with the potential depth $A$ in Figs. 1(d) and 2(b) may primarily be due to the same effect.

In Fig. 5(b), we plot the ratio between the estimated minimum $\theta^{*}$ and the true minimum $\theta_{0}$ as a function of the mean $\mu$, which is an analog of the nonequilibrium driving. We note that here $\theta^{*}$ is the lowest value of $\theta$ at which the slope of the loss function becomes less then $10^{-3}$. We observe that an overestimation regime ($\theta^{*}/\theta_{0}>1$) crosses over to an underestimation regime ($\theta^{*}/\theta_{0}<1$) as $\mu$ grows. This is in striking agreement with the trends shown in Fig. 2(a). The reason why $\theta^{*}$ underestimates $\theta_{0}$ for large $\mu$ can be understood by the flattened loss function landscapes shown in the inset of Fig. 5(b). In this regime, the dynamics of $\theta$ (starting from $\theta$ = 0) slows down, ending up at a value (filled diamonds) even lower than $\theta_{0}$ (empty diamonds). This effect is due to the samples with $x<0$ vanishing when $\mu$ is too large. We expect that a similar mechanism might be at play behind the observed behavior of $\sigma_{\mathrm{pred}}/\sigma$ shown in Fig. 2(a). If we had used a broader range of nonequilibrium driving, the same behaviors might have been observed for other models as well, although this remains to be checked.

The one-parameter model also allows us to examine the effects of the finite minibatch size $M$. While the ideal loss function is given in Eq. (20), the loss function used in the actual training looks like

$\displaystyle\mathcal{L}_{\alpha}(\theta;M)=\frac{1}{M}\sum_{i}{\left[-\frac{(1+\alpha)e^{\alpha\theta X_{i}}-1}{\alpha}+e^{-(1+\alpha)\theta X_{i}}\right]}\;,$

(22)

where $X_{1},\,\ldots,\,X_{M}$ are i.i.d. Gaussian random variables of mean $\mu$ and variance $\sigma^{2}$. When $M$ is large and finite, using the central limit theorem (CLT), the gradient of this loss function can be approximated as [26, 27]

$\displaystyle\frac{\partial\mathcal{L}_{\alpha}}{\partial\theta}\bigg{|}_{\theta=\theta_{t}}=\bar{K}(\theta_{t}-\theta_{0})+\sqrt{\frac{\Lambda}{M}}N_{t}+o(M^{-1/2})\;,$

(23)

where $\theta_{0}=\text{argmin}(\langle\mathcal{L}_{\alpha}(\theta)\rangle)$, $\bar{K}=\partial^{2}_{\theta}\langle\mathcal{L}_{\alpha}\rangle\big{|}_{\theta=\theta_{0}}$, and $\Lambda=\text{Var}\big{[}\partial_{\theta}\mathcal{L}_{\alpha}\big{|}_{\theta=\theta_{0},\,M=1}\big{]}$. When the stochastic gradient descent $\theta_{t+1}=\theta_{t}-\lambda(\partial\mathcal{L}_{\alpha}/\partial\theta)\big{|}_{\theta=\theta_{t}}$ reaches the steady state, the MSE of $\theta$ is given by

$\displaystyle\langle\left(\theta-\theta_{0}\right)^{2}\rangle_{s}=\frac{\lambda\Lambda}{M\bar{K}(2-\lambda\bar{K})}+o(M^{-1})\;.$

(24)

This leading-order behavior is shown in Fig. 5(c) for various values of $\mu$. For all cases, the MSE of $\theta$ is minimized at $\alpha=-0.5$, which is consistent with the smallest error bars observed at $\alpha=-0.5$ in Figs. 1 and 2. Hence, $\alpha=-0.5$ yields the most consistent EP estimator.

Direct measurements of the loss function gradient at the minimum also confirm the above result. As shown in Fig. 5(d), the gradient $\partial_{\theta}\mathcal{L}_{\alpha}$ is far more broadly distributed for $\alpha=0$ than for $\alpha=-0.5$. Moreover, due to the subleading effects (beyond the CLT) of finite $M$, the gradient for $\alpha=0$ features a large skewness. These show that the training dynamics for the original NEEP ($\alpha=0$) tends to be far more volatile and unstable than for the $\alpha$-NEEP with $\alpha=-0.5$.

We proposed the $\alpha$-NEEP, a generalization of the NEEP for estimating steady-state EP at the trajectory level. By choosing a value of $\alpha$ between $-1$ and $0$, the $\alpha$-NEEP weakens the exponential dependence of the loss function on the EP estimator, effectively mitigating the adverse effects induced by poor sampling of transitions associated with large negative EP in the presence of strong nonequilibrium driving and/or deep potential wells. We also observed that $\alpha=-0.5$ tends to exhibit the optimal performance, which can be understood via a simplification of the original EP estimation problem, whose loss function landscape and relaxation properties are analytically tractable. The $\alpha$-NEEP thus provides a powerful method for estimating the EP for much broader range of the nonequilibrium driving force and the time scale of dynamics. Identification of even better loss functions and optimization of other hyperparameters (network size, number of iterations, etc.) are left as future works. It would also be interesting to apply the $\alpha$-NEEP to estimations of the EP of the Brownian movies [28] and stochastic systems with odd-parity variables [29], which have been studied using the original NEEP method.

Acknowledgments. — This work was supported by the POSCO Science Fellowship of the POSCO TJ Park Foundation. E.K. and Y.B. also thank Junghyo Jo and Sangyun Lee for helpful comments.