Machine learning the thermodynamic arrow of time
Alireza Seif Mohammad Hafezi Christopher Jarzynski






Department of Physics, University of Maryland, College Park, MD 20742 
Joint Quantum Institute, NIST/University of Maryland, College Park, MD 20742



Department of Physics, University of Maryland, College Park, MD 20742 
Joint Quantum Institute, NIST/University of Maryland, College Park, MD 20742
Department of Electrical and Computer Engineering, University of Maryland, College Park, Maryland 20742, USA


Department of Physics, University of Maryland, College Park, MD 20742 
Department of Chemistry and Biochemistry, and 
Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742, USA



             


The mechanism by which thermodynamics sets the direction of time's arrow has long fascinated scientists. Here, we show that a machine learning algorithm can learn to discern the direction of time's arrow when provided with a system's microscopic trajectory as input. The performance of our algorithm matches fundamental bounds predicted by nonequilibrium statistical mechanics.  Examination of the algorithm's decision-making process reveals that it discovers the underlying thermodynamic mechanism and the relevant physical observables. Our results indicate that machine learning techniques can be used to study systems out of equilibrium, and ultimately to uncover physical principles.



Machine learning the thermodynamic arrow of time
    Christopher Jarzynski
    December 30, 2023
================================================







§ INTRODUCTION

While the microscopic dynamics of physical systems are time reversible, the macroscopic world clearly does not share this symmetry. 
If we are shown a video of a macroscopic process, it is often easy to guess whether the movie is played in the correct or in time-reversed order, as in the latter case the observed sequence of events is utterly implausible.
In 1927, Sir Arthur Eddington coined the phrase “time's arrow" to express this asymmetry in the flow of events, arguing that it traces back to the second law of thermodynamics <cit.>.


In recent decades there has been increased interest in the out-of-equilibrium physics of microscopic systems, leading to a deepened understanding of non-equilibrium fluctuations and their relation to the second law <cit.>.
In particular, it has become appreciated that fluctuations lead to an effective “blurring” of time's arrow at the nanoscale, and that our ability to discern its direction can be quantified in a system-independent manner <cit.>.



Simultaneously, there have been significant advances in the ability of machine learning (ML) and artificial intelligence (AI) algorithms to tackle practical problems and to automate useful tasks. These include image and video classification <cit.>, medical diagnosis <cit.>, playing games <cit.>, driving cars <cit.>, and most recently analyzing scientific problems such as protein folding <cit.>. ML methods have also emerged as exciting tools to study problems in statistical and condensed matter physics, such as classifying phases of matter, detecting order parameters, and generating configurations of a system from observed data <cit.>. These studies extend and further motivate the use of computer algorithms to learn physics from big data <cit.>. 

In the present work we ask whether a machine can learn to accurately guess the direction of time's arrow from microscopic data, and if so, whether it does so by effectively discovering the underlying thermodynamics, identifying relevant quantities such as work and entropy production. We approach this problem within the framework of nonequilibrium statistical mechanics, numerically generating microscopic trajectories of irreversible physical processes.  
In many of the examples we consider, the system is small and the direction of time's arrow is blurred, in the sense that both a given trajectory and its time-reversed image represent plausible sequences of events.
In these cases the algorithm in principle cannot be perfectly accurate, and it becomes interesting to ask whether it is able to assess its own likelihood to guess the direction of time's arrow correctly.


We find that the machine not only correctly classifies the direction of time's arrow but also approximates the likelihood in the uncertain cases. Moreover, the machine can generate representative trajectories for forward and backward time directions correctly, i.e., it learns what a forward/backward trajectory should look like. We also design a neural network that can detect the underlying process and classify the direction of time's arrow at the same time. 

Finally, we look inside the machine's decision-making process and find that it correctly identifies dissipated work as the key quantity for optimally guessing the direction of time's arrow.

We first introduce the relevant physical laws governing microscopic, non-equilibrium fluctuations. We then briefly review the ML techniques that we will use. Finally, we apply our methods to various model physical examples and we study the ability of ML techniques to learn and quantify the direction of time's arrow.




§ THERMODYNAMICS AND THE ARROW OF TIME


When small systems undergo thermodynamic processes, fluctuations are non-negligible and the second law is expressed in terms of averages.
Thus the Clausius inequality relating the work W performed on a system to the net change in its free energy, Δ F, takes the form 

    ⟨ W ⟩≥Δ F,

where the angular bracket denotes an average over many repetitions of the process.
Moreover, these non-equilibrium fluctuations satisfy strong constraints that allow us to rewrite such inequalities in terms of stronger equalities <cit.>, and to quantify the direction of time's arrow as a problem in statistical inference <cit.>. To frame this problem, let us first specify the class of processes we will study, and introduce notation.

Consider a system in contact with a thermal reservoir at temperature β^-1. The system's Hamiltonian ℋ_λ(𝐱) depends on both the system's microstate  𝐱, and on an externally controllable parameter λ. An external agent performs work by manipulating this parameter. Now imagine that the system begins in equilibrium with the reservoir, and the agent then varies the parameter according to a schedule λ_F(t) from λ_F(0)=A to λ_F(τ)=B. We refer to this as the forward process. The trajectory describing the system's evolution can be pictured as a movie, and is denoted by {𝐱_A→ B(t)}, where the time interval 0≤ t≤τ is implied. We also imagine the reverse process, in which the system starts in an equilibrium state at λ=B, and the agent varies the parameter from B to A according to λ_R(t) = λ_F(τ -t). The trajectory (movie) for this process is denoted by {𝐱_B→ A(t)}. Finally, consider the time reversal of this trajectory, 𝐱̅_B→ A(t)= 𝐱^*_B→ A(τ-t), where the `*' implies negation of momentum coordinates. This time-reversed trajectory corresponds to a movie of the reverse process, played backward in time; the same trajectory may have been achieved during a realization of the forward process, see Fig. <ref>(a). 


Throughout this paper, we will use the term forward trajectory to refer to a trajectory generated during the forward process, i.e. {𝐱_A→ B(t)}, and we will use the term backward trajectory to denote a trajectory generated during the reverse process, but run backward in time, i.e. 𝐱̅_B→ A(t), depicted by the dashed blue line in Fig. <ref>(a).

Guessing the direction of time's arrow can be cast as a game in which a player is shown either a forward or a backward trajectory – thus in either case the player “sees” the parameter being varied from A to B.
The player must then guess which process, forward or reverse, was actually used to generate the trajectory <cit.>. The player's score, or accuracy, is the ratio of correct predictions to the total number of samples.

In order to optimize the likelihood of guessing correctly, it suffices for the player to know the sign of the quantity W - Δ F, where W is the work performed on the system and Δ F = F_B - F_A is the free energy difference between its initial and final states, as depicted in the movie. Specifically, let P(F|{𝐱(t)}) denote the likelihood that a given trajectory, {𝐱(t)}, is obtained by performing the forward process, and let P(R|{𝐱(t)}) denote the likelihood that the trajectory is the time reversal of a realization of the reverse process. Note that P(F|{𝐱(t)})+P(R|{𝐱(t)})=1.  In addition, assume that the game is unbiased, e.g. the choice of performing the forward or reverse process in the first place was decided by flipping a fair coin.
Then the likelihood that the trajectory was generated during the forward process is given by <cit.>

    P(F|{𝐱(t)}) = 1/1+e^-β(W-Δ F),

which is greater than (less than) 50% when W - Δ F is positive (negative).
Here, the work performed by the external agent is  

    W = ∫_0^τ dt λ̇∂ℋ_λ(𝐱)/∂λ,

and the change in free energy is given by 

    Δ F = -1/βlog(Z_B,β/Z_A,β),

where 

    Z_λ,β = ∫ d𝐱exp[-βℋ_λ(𝐱)]

is the partition function. 


In macroscopic systems, the values of work performed on the system corresponding to forward trajectories, W_F, and for backward trajectories, -W_R, are sharply peaked around their mean values, Fig. <ref>(b), and the sign of W-Δ F is a reliable indicator of the direction of time's arrow. (Here, W_R is the work performed during a given realization of the reverse process, therefore for the corresponding backward trajectory the work value is -W_R.)
However, for microscopic systems these distributions can overlap significantly, as in Fig. <ref>(c).
Eq. <ref> shows that the player optimizes the chance of success simply by guessing “forward” whenever W > Δ F, and “reverse” otherwise, without accounting for any further details of the trajectory.
Note that if W - Δ F≫ k_BT then determining the arrow of time is easy, but when W - Δ F≲ k_BT the problem becomes more difficult – in effect, time's arrow is blurred.



§ NEURAL NETWORKS
 


We wish to train a computer program to infer the direction of time's arrow from a movie of the system's trajectory.
To do so, we first simulate a number of trajectories from the forward and the reverse processes, and we “time-reverse” the latter so that each trajectory is chronologically ordered with λ varying from A to B. We attach a label y=0 (reverse) or y=1 (forward) indicating which process was used to generate that trajectory.  We then provide the machine with this collection of labelled trajectories, which will serve as the training data.
A priori, any one of the trajectories could have been generated from either the forward or the reverse process, and the training stage now consists of using a neural network (NN) to construct a model of the function P(F|𝐱(t)), which gives the likelihood the trajectory was generated by the forward process.
Although this function is known analytically, Eq. <ref>, the machine is not provided with this information.
We now sketch how the training is accomplished.

Since each (numerically generated) trajectory consists of a discretized time series of microstates, we represent the trajectory as a matrix 𝐗 whose rows correspond to different times, and whose columns correspond to phase space coordinates.
The training stage amounts to designing a function that maps any such matrix 𝐗 onto a real number p between 0 and 1, whose value is the machine's best estimate of the likelihood that the trajectory was generated by the forward process.

In this work, we consider two types of classifiers: (i) logistic regression (LR), and (ii) convolutional neural network (CNN). The input to LR is a vectorized trajectory 𝐚 = vec(𝐗), and the output is p = g(Ω^⊺𝐚 + b), where Ω is a vector of weights, b is the bias, and g(z) = 1/(1+exp(-z)) is the logistic sigmoid function, see the top panel of Fig. <ref>(d). The CNN can compute more complicated functions than the LR <cit.>. The input to our CNN is a trajectory matrix 𝐗, and the output is again a value p. The CNN has convolutional layers that extract useful information by taking advantage of the temporal and spatial structure of the data, see the bottom panel of Fig. <ref>(d). For details of the CNN architecture see Appendix <ref>.

To train the network, we determine optimal values of parameters (such as the weights and biases in LR) by minimizing the cross-entropy<cit.> 

    C=-1/N_ samp∑_m [ y_m log(p_m) + (1-y_m)log(1-p_m) ]   .

Here the sum is carried over the N_ samp training samples, y_m ∈{0,1} is the label attached to the m'th trajectory, indicating which process was actually used to generate the trajectory, and p_m is the output of the network for that trajectory.

Throughout this work, we always split a given data-set into three parts. We use 60%, 20% and 20% of the data for training, validation, and testing the model, respectively. The validation set is used to tune the architecture and hyperparameters of the model, while the test data is used for unbiased evaluation of the  the final model's accuracy. We use Adam optimizer with parameters suggested in the original paper for the training <cit.>. We assess the performance of the network by testing it over a balanced set of trajectories, i.e. half forward and half backward. If  p_m≥0.5 then the algorithm guesses that the trajectory was generated from the forward process, otherwise it guesses the reverse process.
As a figure of merit, we consider the accuracy, i.e. the ratio of correct guesses to total number of samples. 
The best score that an algorithm can achieve, in the limit of a very large test set, is obtained if the output of the network agrees with the theoretical likelihood  (<ref>), in other words if the algorithm “learns” a result from nonequilibrium statistical physics.

For additional considerations in training NNs see Appendix <ref> 



§ CASE STUDIES

We apply the neural network machinery to detect the direction of the time's arrow and assess the NN's accuracy. We also look at the output of the network and compare it with the theoretical optimal result of Eq. (<ref>). Interestingly, the networks not only learn to guess the direction of the time's arrow but also learn to closely reproduce the likelihood function (see Appendix <ref> for a discussion of the sensitivity of the results to the choice of the activation functions).

We first consider a single Brownian particle in a moving potential. This problem is simple and has an analytical solution. 

We then move on to the more complicated problem of a spin chain with nearest-neighbour coupling in a magnetic field. We consider two scenarios involving the spin chain. First, the coupling is assumed to be constant, and the magnetic field is varied in time. Next, the magnetic field is constant the coupling is changed through time. We refer to the former as the  protocol, and the latter as the  protocol.  

For details of the numerical calculations used to generate the trajectories see Appendix <ref>. For details about the NNs and the number of samples used see Appendix <ref>. Table <ref> summarizes the accuracy of the algorithms studied in this and the following sections.



 §.§ Brownian particle in a moving potential 

Consider an overdamped Brownian particle at temperature β^-1 in a harmonic potential (see Fig. <ref>(a)), evolving according to 

    ẋ = -k/γ (x-λ) + ξ(t),

where k denotes the strength of the potential, λ is the position of the center of the potential, and γ is the damping rate. The noise term, ξ(t) satisfies ⟨ξ(t)ξ(t')⟩ = 2 (βγ)^-1δ(t-t'). 
In the forward protocol, the value of λ is changed from A to B at a fixed rate λ̇=u. Hence the reverse protocol changes λ from B to A with λ̇=-u.  

If the potential is moved rapidly, then in most of the forward trajectories the particle lags substantially behind the potential, whereas in most of the backward trajectories the particle leads the potential. In these cases, the direction of time's arrow is clear. However, if the potential is moved slowly, then the particle stays near the center of the potential, the processes approach the reversible limit, and it becomes difficult to determine the direction of time's arrow.  The theoretical likelihood   (<ref>) is determined by the work W performed and the free energy change, Δ F. Note that in this protocol Δ F = 0. For each trajectory, we calculate W by integrating  

    Ẇ = -k u (x-ut),

We generate samples of the forward and backward trajectories by numerically integrating the stochastic differential equation Eq. (<ref>) (see Fig. <ref>(b) and (c)). We then train a classifier to predict the label for a given trajectory, as described earlier. In Fig. <ref>(e) we compare the accuracy and the output of a LR classifier (grey circles) with the theoretical likelihood (solid curve) obtained from  Eqs. (<ref>) and (<ref>), see Table <ref> and Fig. <ref>(e). 

The seemingly remarkable agreement with the theory can be understood by examining Ẇ (<ref>). Namely, the work W calculated by numerically integrating Ẇ for a given trajectory is linearly related to the sum of the components of that vector. Therefore, LR is well-equipped to calculate this quantity and reproduce the likelihood function. See Appendix <ref> for detailed analysis of the optimal network. 



 §.§ Spin chain - time-dependent field  

Now let us consider a more complicated, many-particle system and a non-linear work protocol. Specifically, we consider a spin chain in a time-dependent magnetic field B(t) and in contact with a thermal reservoir at temperature β^-1, see Fig. <ref>(a), described by a Hamiltonian

    H = ∑_i J σ_i σ_i+1 - B(t) ∑_i σ_i,

where σ_i ∈{-1,+1} is the spin variable at site i, and J is the nearest-neighbour coupling strength.  The dynamics of this system are modeled as a Markov process (See Appendix. <ref>). The Hamiltonian aligns the spin in preferred energy configurations, while thermal fluctuations cause the spins to flip randomly according to a rate related to β. We refer to this example as the 𝖡 protocol.


In the forward process, B(t) changes from a positive value B_0 at t=0 to a negative value -B_0 at t=τ, as shown in Fig. <ref>(b). In the limit where B_0/J ≫β^-1/J ≫ 1, the spins start mostly aligned with the magnetic field. As the field magnitude is lowered, thermal fluctuations become dominant and the spins flip randomly. Eventually, B(t) becomes large and negative, and aligns the spin in the other direction (see the top row of Figs. <ref> (c) and (d)). In this limit it is easy to detect the direction of the time's arrow, as the work distributions have a modest overlap, see Fig. <ref>(e) top row. As the temperature is increased, thermal fluctuations increase the overlap in work distributions, blurring the direction of time's arrow, see the middle and the bottom rows of Figs. <ref>(c), (d), and (e). 

To train the classifier, we generate samples of forward and backward trajectories for three different temperatures using the Metropolis algorithm. The trajectories are matrices with ±1 entries, whose rows and columns correspond to time steps and spin positions, respectively. We are interested in training a single LR classifier that is capable of detecting the direction of time's arrow for different temperatures. Therefore, the information about β is provided through normalizing the elements of the trajectory data with their corresponding temperatures. This matrix is then reshaped as a vector to serve as the input to an LR classifier. We observe that the success of LR in learning both the correct labels and in approximating the likelihood function persists, see Table <ref> and Fig. <ref>(f).  The reason, again, lies in the functional form of W, which can be evaluated by numerically integrating

    Ẇ = -Ḃ(t)∑_i σ_i.

It can be seen that W is proportional to the weighted sum of the elements of the input vector. Note that, in this protocol Δ F =0. Consequently, LR is a perfect model of the likelihood function for all the temperatures, see Appendix <ref>.



 §.§ Spin chain - time-dependent coupling 


In this example, we consider a more complicated version of the spin chain problem, with a ferromagnetic-antiferromagnetic transition. Here, we keep B constant and positive and allow for the time-dependent couplings J(t), see Fig. <ref>(a). We refer to this example as the 𝖩 protocol.  The Hamiltonian is given by

    H = ∑_i J(t) σ_i σ_i+1 - B ∑_i σ_i.



The protocol J(t) is shown in Fig. <ref>(b). In the forward case  J(t) is varied from J_0>0 at t=0 to a -J_0 at t=τ. Note that at low temperatures, where J_0/B≫β^-1/B≫1, the spins start in a state with anti-ferromagnetic ordering. As J(t) grows weaker during the protocol, thermal fluctuations dominate. By the end of the protocol, J(t)=-J_0 and the system settles in a state with ferromagnetic ordering. In this case, the forward and reverse work distributions are distinguishable, and so is the arrow of time, see the top rows of Fig. <ref>(c) - (e). As the temperature increases, so does the overlap between the distributions (the bottom two rows of Fig. <ref>(c) - (e)). 

In this case the work is given by the time integral of

    Ẇ = J̇(t)∑_i σ_i σ_i+1.

We see that W is no longer linearly related to the input, and the LR classifier is incapable of calculating it. Therefore, we use a CNN with periodic boundary condition that can capture more complicated functions. With this CNN, we are able to recover the optimal accuracy again. Note that in this process Δ F≠0, which adds another layer of complexity to the problem. The convolution layer in a CNN has filters that can capture the two-body nearest-neighbor correlations required to calculate the work, without introducing too many parameters. In fact, for a single temperature, we are able to analytically derive the parameters of a CNN that exactly calculate the likelihood function. The performance and the output of the network are shown in Table <ref> and Fig. <ref>(f).   For more details on the optimal network construction and the performance of sub-optimal strategies see Appendix <ref>. 





§ INTERPRETATION AND EXTENSIONS 

In this section, we use three approaches to investigate trained networks and to develop insight into what they have learned. 


First, we use inceptionism techniques <cit.> to learn the network's ideal representative of forward and backward trajectories. 

 Specifically, we use gradient descent on a random input such that the trained networks in Secs. <ref> and <ref> output 1 or 0 corresponding to forward and backward trajectories, respectively. This is in contrast with the previous section where we optimized for the weights and biases of the network. Among the simulated trajectories in the test set, we choose one with p≈0.5 <cit.> – this is a trajectory for which the classifier has difficulty assigning the direction of time's arrow. We project the configurations to discrete values after each step of the gradient descent, and demand that there be at most 1 spin-flip per time step, to ensure that the network `dreams' of physically realizable trajectories. We find that the networks' ideas of the forward and backward trajectories show strong agreement with the true physical picture, see Fig. <ref>(a). 


Secondly, to assign a physical interpretation to the networks' decision-making process, we project the trajectories onto a two-dimensional reduced phase space corresponding to the collective coordinates {x̃^(1)(t)} = {∑_i σ_i(t)} and {x̃^(2)(t)} = {∑_i σ_i(t)σ_i+1(t)} (taking period boundary conditions), representing  magnetization and nearest-neighbour correlations, respectively.
We also replace the value of x̃^(ℓ)(t) within each time window of ten time steps, by the sum of the values within that window.
By thus coarse-graining in both phase space and time, we reduce the noise due to finite size effects and variations over samples.
Next, we use these coarse-grained trajectories to train LR classifiers for both protocols in Secs. <ref> and <ref> (See Table <ref> for the performance of these networks). Finally, we investigate the weights Ω^(ℓ) that the networks assign to the magnetization (ℓ=1) and the nearest neighbor correlations (ℓ=2).  Fig. <ref>(b) reveals that for the  protocol (top row), the network mostly cares about the magnetization, whereas when the  protocol is performed (bottom row), the network bases its decision on the nearest-neighbor correlations. Moreover, the learned values of Ω^(ℓ) agree with our analytical results that reproduces the correct likelihood value (see Appendices <ref> and <ref> for details). These observations suggest that the network learns that the time derivative of the Hamiltonian, and by extension the work (<ref>), is an important feature in guessing the direction of time's arrow. We note that when the process is highly irreversible, the distributions of the forward and reverse work are well-separated. In this case, the network easily determines the arrow of time, but does not learn about the importance of work and bases its decision on other visible differences in the trajectories, see Appendix <ref>. 







Lastly, we ask whether a single algorithm can learn to accurately guess the direction of time's arrow for trajectories generated using multiple protocols, when the identity of the protocol is not specified.
One approach is to take a large neural network and train it on trajectories from both 𝖡 and 𝖩 protocols. We have found that this approach works to a certain degree, but does not reach the accuracy of the individual networks in Secs. <ref> and <ref>. However, by using our knowledge about the structure of the problem we can design an algorithm that attains the optimal performance. Specifically, we use a  mixture of experts (MoE), with an output that is the weighted sum of expert networks <cit.>. When the protocol is not specified, the net forward likelihood is

    P(F|𝐗) = P(F|𝐗,𝖡)P(𝖡|𝐗)+P(F|𝐗,𝖩)P(𝖩|𝐗).

The quantities P(F|𝐗,𝖡) and P(F|𝐗,𝖩) are modeled using neural networks similar to those considered in Secs. <ref> and <ref>, respectively. These networks are referred to as experts.  Additionally, we use a CNN  to model P(𝖡|𝐗) = 1 -P(𝖩|𝐗). This CNN, which is called the gating network, learns the protocol from trajectories. Therefore, we obtain a larger three-headed network by combining the output of the three neural networks as in Eq. (<ref>), as illustrated in Fig. <ref>(a). For the training, we use the pre-trained expert networks for the 𝖡 and 𝖩 protocols, and optimize the cost function (<ref>) over sample trajectories from both protocols. We observe that the performance of this network is similar to that of the individual networks, as the gating network learns to accurately identify the protocol of input trajectories (see Fig. <ref>(b)). Note that the predictions of the gating network are more accurate at lower temperatures. This makes sense as the distribution of the initial state in the two protocols are distinguishable in low temperatures, but become less so as the temperature is increased.  



§ CONCLUSION AND OUTLOOK

Starting with a simple, solvable harmonic oscillator model, then proceeding to more complicated spin systems, we have shown that machine learning algorithms can be trained to discern the direction of time's arrow in irreversible thermodynamic processes.
We have found that neural networks not only learn to guess the direction of time's arrow but also to accurately evaluate the likelihood that the guess is correct, when the direction of the arrow is not entirely clear.
Moreover, we have used various techniques to interpret what the network learns.
In particular, by examining the optimized parameter values that emerge from the training, we have been able to identify which physical quantities the network uses to guess the direction of time's arrow.
In this sense, our study represents a step toward AI driven discovery of physical concepts.

Machine learning techniques have been applied extensively to the study of equilibrium statistical physics <cit.>.
Our results extend this computational toolkit to out-of-equilibrium phenomena.
While we have focused on the arrow of time, we expect that other important issues and questions in non-equilibrium physics can usefully be studied with these tools.
We anticipate that the techniques considered in this work can be extended to estimate free energy differences, as well as to identify physical quantities that distinguish different regimes of dynamics in out-of-equilibrium quantum phenomena.
In addition, using unsupervised learning techniques such as generative modeling may be especially useful in studying non-equilibrium phenomena <cit.>. Unlike the equilibrium case where the state of the system is given by the Boltzmann distribution, the general form of the non-equilibrium steady-state is not known. Generative models are an ideal candidate to model and learn these distributions.


Moreover, machine learning researchers have shown that ML techniques can be used to detect the playback direction of real-world videos  <cit.>. These studies are concerned with videos of macroscopic objects that are in principle irreversible, and the arrow of time has a clear direction. In such scenarios, there are many indicators that can reveal the true playback direction, and therefore it is hard to quantify the optimal performance. However, in the physical examples the optimal attainable accuracy of the classifier is dictated by the laws of physics. Therefore, problems with large number of phase-space coordinates and with complicated dynamics, such as the  protocol for 2D Ising model, can serve as a standardized benchmark for video classification algorithms.   





§ ACKNOWLEDGEMENTS

AS thanks Evert van Nieuwenburg, Grant Rostkoff, and Ali Izadi Rad for helpful discussions. AS and MH gratefully acknowledge support from ARO-MURI and Physics Frontier Center by National Science Foundation at the JQI and CJ from the National Science Foundation under grant DMR-1506969.



apsrev4-1






§ NEURAL NETWORKS 



 §.§ Convolutions

A convolution layer convolves the input with a number of filters, and then applies a non-linear function to the output of the filters. Each convolution operation with a kernel Ω and bias b, maps an input matrix 𝐗,  to another matrix 𝐙 = Ω*𝐗 given by <cit.>

    Z_j,k= ∑_m,n X_j× s+m,k× s+nΩ_m,n + b

where s specifies the number of steps the filter moves in each direction. It is called the stride of the convolution and is a hyperparameter that is tuned using the cross-validation data. The output of the convolution layer is obtained by applying a non-linear function g element-wise to 𝐙. The convolution layers can be repeated many times, and combined with pooling layers where the dimension of the output is reduced through a procedure such as averaging. At the end, the output of the convolution layer is flattened to form a vector and that vector is fed into a series of fully connected layers to produce the network's output <cit.>.

The CNN's that we consider has four 2×2 filters, with the stride of 1, and with periodic boundary condition. We choose the rectifier, i.e., g(z) = max(0,z), for the activation of these filters. The output of all the filters is then combined to form a single vector. For the CNN classifying the  protocol (Sec. <ref>), this vector is fed into a single neuron with sigmoid activation, whose values determine the direction of time's arrow. For the gating network (Sec. <ref>), this vector is fed into a fully connected layer with 50 hidden neurons and the rectifier activation, followed by the output neuron with the sigmoid activation.  



 §.§ Regularization and sample size

To reduce overfitting it is helpful to include a regularization term.  This will help to reduce the difference between the training error and the test error. We consider L_2 regularization α∑_ℓΩ_ℓ^2, that is adding the square of all the weights in the network to the cost function. The parameter α is a hyper-parameter of the model and is tuned using the cross-validation data. 

Additionally, in training the CNN in Sec. <ref>, we use the dropout technique to reduce overfitting. Dropout refers to deactivating and ignoring certain neurons during the training phase. Specifically, at every training step, a random fraction of p_drop of neurons are deactivated <cit.>. 

We find that the performance of our algorithms does not vary significantly with the choice of hyper-parameters. We choose p_drop=0.25 for the dropout rate of neurons of the convolutional layer in the  network, and p_drop=0.5 for the gating network. The L_2 regularization rates are shown in Table <ref>.


Another important quantity in training the neural networks is the sample size. We use a total of 20000 samples for the Brownian particle. For the spin chain examples ( and  protocols), we use 20000 samples for each temperature. The samples are then split into three sets and are used to train, validate, and test the models.



§ GENERATING THE DATA 

To generate trajectories we closely follow Ref. <cit.>. We consider a discrete set of time steps t∈{0,1,…,τ}. The value of the control parameter and the state of the system at each time step is denoted by λ_t and 𝐱_t, respectively. Note that in dealing with discrete time steps, rather than using u(t), we use the notation u_t for the value of variable u at the time step t. In the forward process, the initial state of the system is drawn from equilibrium with λ = λ_0. The time evolution can be broken into two substeps: 

    
  * With the state of the system fixed, the control parameter is changed λ_t →λ_t+1
    
  * At fixed λ_t+1, the state of the system evolves 𝐱_t →𝐱_t+1

Here, the second substep is either generated by a stochastic differential equations (Sec. <ref>) or Metropolis algorithm (Secs. <ref>  <ref>). 
The total work performed in this process is 

    W = ∑_t=0^τ-1 [ℋ_λ_t+1(𝐱_t) - ℋ_λ_t(𝐱_t) ]

For producing backward trajectories, the system is initialized in an equilibrium state with λ=λ_τ. The dynamics begin with a change in the system state, followed by a change in λ. In the end, the history of the system state is reversed, and the calculated work is negated to obtain backward trajectories and their corresponding work values. 


§ OPTIMAL NETWORKS 

For some of the examples that we considered, it is possible to derive an analytical expression for the optimal weights and biases of the network. Specifically, we examine the expression that is used to calculate the work W and the change in free energy Δ F. Because the logistic sigmoid activation function, i.e., g(z)=1/(1+exp(-z)), used for classification coincides with the form of the likelihood function (<ref>) in the arrow of time problem, we are able to find the networks parameters {Ω,𝐛} that reproduce the same likelihood function.  To illustrate, consider the LR model with the output p=1/(1+exp(-z)) and z=Ω^⊺𝐚 + b. If we find Ω and b such that z = β(W-Δ F), the output of the network p correctly represent P(F|𝐗). In the following we show that when W (<ref>), is linear in the elements of 𝐱_t and is subsequently linear in elements of 𝐚, we are able to find such optimal Ω and b.  



 §.§ Brownian particle in a moving potential

In this example, the system's state at each time step is described by a scalar x_t, i.e. the position of the particle. We have a total of τ+1 time steps, therefore the input to the NN is a τ+1 dimensional vector. The LR classifier considered here, is parameterized by a τ+1 dimensional weight vector with elements Ω_t for t=0,…,τ and a bias b. Using Eq. (<ref>) we find

    W = ∑_t=0^τ-1δλ(k x_t - kλ_t +1/2 k δλ),

where δλ = λ_t+1 - λ_t is independent of t, because the protocol is linear. Note that Δ F =0 in this example. With the choice of  

    Ω_t    = β k δλ   (t≠τ) 
    
        0     (t=τ)
        ,
    
        b    =  β∑_t=0^τ-1δλ(- kλ_t +1/2 k δλ),

we can see that Ω^⊺𝐚 + b = β(W-Δ F), where (𝐚)_t = x_t.



 §.§ Spin chain -  protocol 

The full trajectory of an n spin system over τ time steps is represented by a τ× n matrix 𝐗. We denote the orientation (up or down) of the ith spin at time t with X_t,i = ± 1. The input to LR classifier, is a vector obtained from rearranging the trajectory matrix 𝐗 to shape it into an τ n × 1 array. By using Eq. (<ref>) we find that 

    W = -∑_t=0^τ-1 (δ B_t ∑_i=1^n X_t,i),

where δ B_t = B_t+1-B_t. Work calculated using Eq. (<ref>) is the discrete time version of W obtained from Eq. (<ref>). In this example, Δ F = 0 again, and we find that the optimal weights and bias are given by

    Ω_t,i   = - βδ B_t      (t≠τ),
    
        b    = 0,

where Ω_τ,i = 0. Note that if the input X_t,i is scaled with β, a single LR classifier is able to reproduce the correct likelihood for different temperatures. 



 §.§ Spin chain -  protocol 

Using the same notation in the previous section, we find that the work, W, is given by

    W = ∑_t=0^τ-1 (δ J_t ∑_i=1^n X_t,iX_t,i+1),

where δ J_t = J_t+1-J_t. Work calculated using Eq. (<ref>) is the discrete time version of W obtained from Eq. (<ref>).  We also use periodic boundary condition, which implies X_n+1 = X_1. Moreover, there is a non-zero change in the free energy, which is given by <cit.>

    Δ F = -1/βlog(ϵ_-(β B, β J_0)^n+ϵ_+(β B,β J_0)^n/ϵ_-(β B,-β J_0)^n+ϵ_+(β B,-β J_0)^n)

where 

    ϵ_±(β B,β J)    = exp(β  J) cosh (β  B)
       ±√(exp (2 β  J) cosh ^2(β  B)-2 sinh (2 β  J)).

We see that it is not possible to have a logistic regression model that calculates W (see Fig. <ref>). 

However, CNN's can, in principle, learn the relevant representation (i.e. nearest-neighbor correlations) from the input data, and learn the corresponding weights to calculate work. Specifically, to show that in principle a CNN can exactly calculate the correct likelihood from the input, we consider a CNN with four 1×2 filters with periodic boundary condition and the rectifier activation g_1, followed by the output layer with sigmoid activation g_2. We set the biases of the convolutional layer to zero, and choose the weights Ω^[1,i] for filters i=1,2,3,4 as follows

    Ω^[1,1]
    Ω^[1,2]
    Ω^[1,3]
    Ω^[1,4]=
         
        (1     1)
    
        (1     -1)
    
        (-1     1)
    
        (-1     -1)
        .

Each filter is only activated for one of the possible configuration of two neighboring spins. Specifically, given the τ× n input 𝐗, the output of each filter g_1(Ω^[1,j] * 𝐗) is a τ× n matrix  𝐗̃^(j) such that 

    X̃^(1)_t,i   = 1    if    (X_t,i X_t,i+1) = (1,1),  
    X̃^(2)_t,i   = 1    if    (X_t,i X_t,i+1) = (1,-1),  
    X̃^(3)_t,i   = 1    if    (X_t,i X_t,i+1) = (-1,1),  
    X̃^(4)_t,i   = 1    if    (X_t,i X_t,i+1) = (-1,-1),

and X̃^(j)_t,i = 0 otherwise. We can now rewrite the output of the network as 

    g_2(b_2 + ∑_j (∑_t,iΩ^[2,j]_t,iX̃^(j)_t,i))),

where Ω^[2,j] contains the weights of the output layer corresponding to 𝐗̃^(j). 
The optimal values of these weights are given by 

    Ω^[2,1]_t,i=Ω^[2,4]_t,i   = - βδ J_t      (t≠τ),
    Ω^[2,2]_t,i=Ω^[2,3]_t,i   = + βδ J_t      (t≠τ),
    
        b    = -βΔ F,

where Ω^[2,j]_τ,i = 0. 

Therefore, a CNN with four 1 × 2 filters is sufficient to capture the likelihood at a single temperature of β^-1. We find that in practice, a CNN with such an architecture is likely to get stuck at local minima, and finding the optimal parameters shown above greatly depends on the initial weights of the network. However, we observe that a CNN with four 2 × 2 filters can achieve a close to optimal performance more easily. 





§ COARSE-GRAINED FEATURES 


To reduce the parameters of the neural network and simplify the task of learning we pre-calculate a set of features for the network. Specifically, for the two protocols concerning the spin chain in a magnetic field, the coarse-grained features are 

    x̃^(1)_s   = ∑_t=m s^m (s+1)-1∑_i=1^n X_t,i,
    x̃^(2)_s   = ∑_t=m s^m (s+1)-1∑_i=1^n X_t,iX_t,i+1,

where m is an integer and s is the scaled time. Using this feature map, LR classifier can calculate W for both  and  protocols (See Fig. <ref>). The input to the network is a 2τ/m dimensional vector [ 𝐱̃^(1); 𝐱̃^(2) ]. We denote the weights corresponding to 𝐱̃^(ℓ) by τ/m dimensional vectors Ω^(ℓ) for ℓ = 1,2. 
In this case, we approximate the optimal weights and bias of the networks by their average over the coarse-grained time window. For the  protocol we find 

    Ω^(1)_s    = β/m∑_t=m s^m (s+1)-1δ B_t,
    Ω^(2)_s    = 0,
    
        b    = 0.

Similarly, the weights and bias for the  protocol are given by 

    Ω^(1)_s    = 0,
    Ω^(2)_s    = β/m∑_t=m s^m (s+1)-1δ J_t,
    
        b    = -βΔ F.

In both cases we can see that β(W-Δ F)≈(Ω^(1))^⊺𝐱̃^(1)+(Ω^(2))^⊺𝐱̃^(2), where the approximation comes from coarse-graining. 

Note that in this case, even though LR classifier can calculate β W if the input is scaled with β, it is not possible to train the network over different temperatures. This is because Δ F ≠ 0, and a simple bias cannot capture multiple values of βΔ F. Therefore, we only consider a single temperature for the  protocol in studying the optimal networks with coarse-grained features.  



§ ALTERNATIVE ACTIVATION FUNCTIONS 

The logistic function that appears in the theoretically calculated likelihood in the time's arrow problem, is similar to the activation function that is commonly used for classification in machine learning. To assess the general ability of the networks in approximating the likelihoods we try different activation functions. Specifically, we choose 
g(z) = exp(z^2) and g(z)=sin(z)^2 as the activation of the last layer of the neural network so that the output is always between 0 and 1, and can be interpreted as probabilities.  We also add a hidden layer to give the network the ability to calculate complex functions.
We compare the network's output with the theoretical likelihoods for the spin chain under the  protocol with coarse-grained feature discussed in Sec. <ref>. We only train the network at a single temperature and observe that the network can still approximate the likelihood function as shown in Fig. <ref>. 




§ HIGHLY IRREVERSIBLE PROCESSES 

 When the process is highly irreversible, the arrow of time has a clear direction. In such cases, we observe that while LR has 100% accuracy, it does not learn work (obtained by inverting the sigmoid function in the output). This is because, the events that enable the network to learn work are extremely rare, and are usually absent in the training data. However, there are other evident differences that can show the direction of time's arrow. We show an example of such a process for the  protocol at low temperatures. We observe that the orientation of spins undergoes a sharp transition as the magnetic field changes sign. However, the time that this transition occurs is different in the forward and backward trajectories. The classifier makes a decision based on the spin configuration at this particular time.



