Capabilities and Limitations of Time-lagged Autoencoders for Slow Mode Discovery in Dynamical Systems

TL;DR

This paper analyzes the capabilities and limitations of time-lagged autoencoders in discovering slow modes in dynamical systems, providing theoretical bounds and numerical comparisons with SRVs.

Contribution

It offers the first rigorous theoretical analysis of nonlinear TAEs and compares their performance with SRVs in identifying slow modes.

Findings

TAEs learn a mixture of slow and maximum variance modes

TAEs can fail to identify the slowest mode in certain systems

SRVs outperform TAEs in some cases of slow mode discovery

Abstract

Time-lagged autoencoders (TAEs) have been proposed as a deep learning regression-based approach to the discovery of slow modes in dynamical systems. However, a rigorous analysis of nonlinear TAEs remains lacking. In this work, we discuss the capabilities and limitations of TAEs through both theoretical and numerical analyses. Theoretically, we derive bounds for nonlinear TAE performance in slow mode discovery and show that in general TAEs learn a mixture of slow and maximum variance modes. Numerically, we illustrate cases where TAEs can and cannot correctly identify the leading slowest mode in two example systems: a 2D "Washington beltway" potential and the alanine dipeptide molecule in explicit water. We also compare the TAE results with those obtained using state-free reversible VAMPnets (SRVs) as a variational-based neural network approach for slow modes discovery, and show that SRVs can correctly discover slow modes where TAEs fail.

Equations

$$A(z)=\frac{\mathbb{E}\left[\delta z(x_t)\delta z(x_{t+\tau })\right]}{{\sigma }^2(z)},$$

$$d_{\tau }$$

$$=\mathbb{E}\left[{\Vert D(z_t)-x_{t+\tau }\Vert }^2\right]$$

$$=\mathbb{E}\left[{\Vert D(E(x_t))-x_{t+\tau }\Vert }^2\right].$$

$$D(z)=(c_1{x}_{t,1}+c_0,0)(c_1\ne 0),$$

$$d_{\tau }(z\sim x_{t,1})$$

$$=\mathbb{E}\left[{\Vert (c_1{x}_{t,1}+c_0,0)-(x_{t+\tau ,1},x_{t+\tau ,2})\Vert }^2\right]$$

$$\stackrel{\text{mean-free features}}{=}\mathbb{E}\left[c_1^2{x}_{t,1}^2+c_0^2+x_{t+\tau ,1}^2-2x_{t+\tau ,1}{c}_1{x}_{t,1}+x_{t+\tau ,2}^2\right]$$

$$\stackrel{\text{process is stationary}}{=}\mathbb{E}\left[c_1^2{x}_{t,1}^2+c_0^2+x_{t,1}^2-2x_{t+\tau ,1}{c}_1{x}_{t,1}+x_{t,2}^2\right].$$

$$\{\begin{array}{ll}{c}_0& =0,\\ c_1& =\mathbb{E}\left[x_{t,1}{x}_{t+\tau ,1}\right]/\mathbb{E}\left[x_{t,1}^2\right]=A(x_{t,1}),\end{array}$$

$$d_{\tau }(z\sim x_{t,1})={\sigma }^2(x_{t,1})(1-A^2(x_{t,1}))+{\sigma }^2(x_{t,2}),$$

$$d_{\tau }(z\sim x_{t,1})$$

$$=$$

$$d_{\tau }(z\sim x_{t,2})$$

$$=$$

$${\sigma }^2(x_{t,1})A^2(x_{t,1})<{\sigma }^2(x_{t,2})A^2(x_{t,2}),$$

$$d_{\tau }(z\sim x_{t,1})>d_{\tau }(z\sim x_{t,2}),$$

$$d_{\tau }(z\sim (b_1{x}_{t,1}+b_2{x}_{t,2}))$$

$$=$$

$$d_{\tau }(z\sim b_1{x}_{t,1}+b_2{x}_{t,2})=\{\begin{array}{ll}{\sigma }^2(x)-{\sigma }^2(x_{t,1})A^2(x_{t,1}),& \text{for }{b}_2=0\\ {\sigma }^2(x)-{\sigma }^2(x_{t,1})A^2(x_{t,1})-\frac{{\sigma }^2(x_{t,2})({\sigma }^2(x_{t,2})A^2(x_{t,2})-{\sigma }^2(x_{t,1})A^2(x_{t,1}))}{(1/b_2^2-1){\sigma }^2(x_{t,1})+{\sigma }^2(x_{t,2})},& \text{for }0<b_2<1\\ {\sigma }^2(x)-{\sigma }^2(x_{t,2})A^2(x_{t,2}),& \text{for }{b}_2=1\end{array}$$

$${\sigma }^2(\bar{x}(z))={\sigma }^2({\mathbb{E}}_x\left[x|z\right])$$

$$d_{\tau }(z)\approx {\sigma }^2(\bar{x}(z))(1-G^2(z))+\left({\sigma }^2(x)-{\sigma }^2(\bar{x}(z))\right)$$

$$\frac{V(r,\theta )}{k_{B}T}=\{\begin{array}{l}0\text{, if }{r}_1-dr<r<r_1+dr\\ 0\text{, if }{r}_2-dr<r<r_2+dr\\ 4\text{, if }{r}_1+dr<r<r_2-dr\\ 1.25+7.5\left(\left|r-r_1\right|+\left|r-r_2\right|\right)\text{, otherwise}\end{array}$$

$$p_{ij}=\{\begin{array}{l}{C}_i{e}^{-(V_j-V_i)/(k_{B}T)},\text{if i,j are neighbors or i=j}\\ 0,\text{otherwise}\end{array}$$

$$d_{SRV}=-A(z)=-A(E(x_t)).$$

$$x_1=$$

$$x_2=$$

$$r=$$

$$x_1=$$

$$x_2=$$

$$r=$$

$$d_{\tau }^M=\frac{\mathbb{E}\left[{\Vert E(x_t)-E(x_{t+\tau })\Vert }^2\right]}{{\sigma }^2(E(x_t))},$$

$$\mathbb{E}\left[{\Vert E(x_t)-E(x_{t+\tau })\Vert }^2\right]$$

$$=$$

$$d_{\tau }^M=2-2A(E(x_t))$$

$$d_{VDE}=\lambda \left(\mathbb{E}\left[{\Vert D(z_t)-x_{t+\tau }\Vert }^2\right]+L_{KL}\right)-(1-\lambda )A(z),$$

Full text

Introduction

Estimation of the slow (i.e., maximally autocorrelated) collective modes of a dynamical system from trajectory data is an important topic in dynamical systems theory in understanding, predicting, and controlling long-time system evolution Andrew et al. (2013); Mardt et al. (2018); Pathak et al. (2018); Ye et al. (2015); Giannakis and Majda (2012); Korda and Mezić (2018a); Sharma et al. (2016); Korda and Mezić (2018b); Koopman (1931); Mezić (2005). In the context of molecular dynamics, identification of the leading slow modes is of great value in illuminating conformational mechanisms, constructing long-time kinetic models, and guiding enhanced sampling techniques Noé and Nuske (2013); Nüske et al. (2014); Schütte et al. (2001); Prinz et al. (2011); Schwantes and Pande (2015, 2013); Pande et al. (2010); Chodera and Noé (2014, 2014); Pérez-Hernández et al. (2013); Schwantes and Pande (2015); Harrigan and Pande (2017); Mardt et al. (2018); Hernández et al. (2018); Sultan et al. (2018); Wehmeyer and Noé (2018). Many machine learning models have been applied to learn slow modes from molecular trajectory data, some falling into the category of more traditional techniques, including time-lagged independent component analysis (TICA) Pérez-Hernández et al. (2013); Noé and Nuske (2013); Nüske et al. (2014); Noé and Clementi (2015); Noé et al. (2016); Pérez-Hernández and Noé (2016); Schwantes and Pande (2013); Klus et al. (2018); Husic and Pande (2018), kernel TICA Schwantes and Pande (2015); Harrigan and Pande (2017), Markov state models (MSMs) Prinz et al. (2011); Pande et al. (2010); Schwantes and Pande (2013); Husic and Pande (2018); Trendelkamp-Schroer et al. (2015); Sultan and Pande (2017); Mittal and Shukla (2018); Harrigan et al. (2017); Wehmeyer et al. (2018); Scherer et al. (2018), while others employ more recently-developed deep learning models, including time-lagged autoencoders (TAEs) Wehmeyer and Noé (2018), variational dynamics encoders (VDEs) Wayment-Steele and Pande (2018); Hernández et al. (2017); Sultan et al. (2018), variational approach for Markov processes nets (VAMPnets) Mardt et al. (2018), and state-free reversible VAMPnets (SRVs) Chen et al. (2019). These approaches all employ variants of deep neural networks, but differ in the details of their architecture and loss function: TAEs are a regression approach that minimize time-lagged reconstruction loss, VAMPnets and SRVs are variational approaches that maximize autocorrelation of the slow modes, and VDEs can be conceived as a mixture of the regression and variational approaches. Although there are existing theoretical guarantees of slow modes discovery for variational approaches Chen et al. (2019); Noé and Nuske (2013); Nüske et al. (2014); Schütte et al. (2001), similar theoretical guarantees for regression approaches is currently limited to linear cases Wehmeyer and Noé (2018). Specifically, linear TAEs are known to be equivalent to time-lagged canonical correlation analysis, and closely related to TICA and kinetic maps Wehmeyer and Noé (2018); Wu and Noé (2017); Hotelling (1936); Noé and Clementi (2015); Pérez-Hernández et al. (2013). In this work, we aim to fill this gap by presenting a theoretical and numerical analysis of the capabilities and limitations of TAEs as a nonlinear regression approach for slow mode discovery.

Results and Discussion

Consider a trajectory of a dynamical system $\{x_{t}\}$ where $x_{t}$ is a system configuration, or a derived featurization of the configuration, at time $t$. We define the slowest mode for a given lag time $\tau$ as the functional mapping $z(\cdot)$ that maximizes autocorrelation $A(z)$ for a lag time $\tau$,

[TABLE]

where $\delta z(x_{t})=z(x_{t})-\mathop{\mathbb{E}}\left[z(x_{t})\right]$ is the mean-free slow mode and $\sigma^{2}(z)$ is its variance. This definition is closely related to the dominant eigenfunction of the transfer operator Chen et al. (2019); Noé and Nuske (2013); Nüske et al. (2014); Schütte et al. (2001). A TAE seeks to estimate the slowest mode by training a time-lagged autoencoder to encode a system configuration $x_{t}$ at time $t$ into a low-dimensional latent space $z_{t}=E(x_{t})$, and then decode this latent space embedding to reconstruct the system configuration $x_{t+\tau}=D(z_{t})=D(E(x_{t}))$ at time $(t+\tau)$. The operational principle is that minimizing the reconstruction loss $\mathop{\mathbb{E}}\left[\left\lVert D(E(x_{t}))-x_{t+\tau}\right\rVert^{2}\right]$ at a lag time $\tau$ promotes discovery of slow modes $z_{t}=E(x_{t})$ within the latent space. In the following sections, we define under what conditions TAEs are able to correctly learn the slowest mode and when they will fail to do so. To simplify our discussion, we restrict our analysis to recovery of the leading slowest mode of the system. A schematic of a fully-connected feedforward TAE with a 1D latent space is presented in Fig. 1.

.1 Linear time-lagged autoencoders (TAEs) can learn the slowest mode by employing whitened features

Consider a sufficiently long 2D trajectory $\{x_{t}\}=\{(x_{t,1},x_{t,2})\}$ for a stationary process such that the mean and variance do not change over time and the two components $x_{t,1}$ and $x_{t,2}$ are mean-free and mutually independent. Let the autocorrelation and the variance for component $i\ (i=1,2)$ be $A(x_{t,i})$ and $\sigma^{2}(x_{t,i})$, respectively. Let $z_{t}=E(x_{t})$ be the latent variable, where $E$ is the encoder mapping for TAE, and $\tilde{x}_{t+\tau}=D(z_{t})$ be the reconstructed time-lagged output, where $D$ is the decoder. The TAE seeks to find the encoding and decoding functional mappings $E$ and $D$ to minimize the time-lagged reconstruction loss,

[TABLE]

Let us assume that $A(x_{t,1})>A(x_{t,2})$, which defines $x_{t,1}$ to be a slower component than $x_{t,2}$. If $z\sim x_{t,1}$ (denoting that $z$ is a bijection of $x_{t,1}$, which in the linear case is a non-trivial linear transformation), the reconstructed output $D(z)$ should be a linear function of $x_{t,1}$ given by,

[TABLE]

where $c_{1}$ and $c_{0}$ are constants, and the second component is 0 since $z$ does not contain information about $x_{t,2}$. The corresponding time-lagged reconstruction loss is given by,

[TABLE]

The optimal reconstruction coefficients are given by minimization of Eq. .1 with respect to $c_{0}$ and $c_{1}$,

[TABLE]

and the corresponding minimal loss is given by,

[TABLE]

where we employed the substitution $\mathop{\mathbb{E}}\left[x_{t,1}x_{t+\tau,1}\right]=A(x_{t,1})\mathop{\mathbb{E}}\left[x_{t,1}^{2}\right]=A(x_{t,1})\sigma^{2}({x_{t,1}^{2}})$. We identify the first term $\sigma^{2}(x_{t,1})(1-A^{2}(x_{t,1}))$ as the “propagation loss”, which increases as autocorrelation decreases and therefore generally increases with lag time $\tau$. We identify the second term $\sigma^{2}(x_{t,2})$ as the “irreducible capacity loss”, which reflects the fact that the one-dimensional latent variable $z$ does not contain any information about component $x_{t,2}$, and is independent of lag time. Eq. 6 can be rearranged as,

[TABLE]

where $\sigma^{2}(x)$ is the time-independent total variance of configurations.

If we now consider the case that $z\sim x_{t,2}$ by an analogous analysis the loss is given by,

[TABLE]

From Eq. .1 and Eq. .1, we see that in the time-lagged reconstruction loss there are two contributing factors: variance and autocorrelation. By construction $A(x_{t,1})>A(x_{t,2})$, so it is the objective of the TAE to learn $z\sim x_{t,1}$ as the slowest mode. However, if $\sigma^{2}(x_{t,2})$ is sufficiently large compared to $\sigma^{2}(x_{t,1})$ such that,

[TABLE]

then,

[TABLE]

and the TAE loss is minimized by learning $z\sim x_{t,2}$.

We also consider whether it is possible that the optimal linear mode is actually a mixed linear mode where $z\sim(b_{1}x_{t,1}+b_{2}x_{t,2})$ and $b_{1}^{2}+b_{2}^{2}=1$. It can be shown by an analogous analysis, and recalling that the two components are independent and mean free, that the minimal loss is,

[TABLE]

and that this expression reduces to Eq. .1 for $b_{2}=0$ and Eq. .1 for $b_{1}=0$ as expected.

Using $b_{1}^{2}+b_{2}^{2}=1$ to eliminate $b_{1}$, the minimal loss can be simplified to,

[TABLE]

which is a monotonic function with respect to $b_{2}$. Recalling that all variances and squared autocorrelations constrained to non-negative values, if $\sigma^{2}(x_{t,1})A^{2}(x_{t,1})<\sigma^{2}(x_{t,2})A^{2}(x_{t,2})$ then the loss is minimized for $b_{2}=1$ and $z\sim x_{t,2}$, whereas if $\sigma^{2}(x_{t,1})A^{2}(x_{t,1})>\sigma^{2}(x_{t,2})A^{2}(x_{t,2})$ then the loss is minimized for $b_{2}=0$ and $z\sim x_{t,1}$. Accordingly, except for the case $\sigma^{2}(x_{t,1})A^{2}(x_{t,1})=\sigma^{2}(x_{t,2})A^{2}(x_{t,2})$, the loss is globally minimized by learning one of the pure component modes, the optimum is not a mixture of the two modes, and the analysis presented for the two pure modes is sufficient and complete for the case of a linear TAE with independent components.

This theoretical development demonstrates that a linear TAE is not guaranteed to find the slowest mode in the event that the associated variance of a faster mode is sufficiently large such that its erroneous identification leads to a smaller reconstruction loss. A straightforward solution to this issue is to apply a whitening transformation Wehmeyer and Noé (2018) to the input data such that $\sigma^{2}(x_{t,1})=\sigma^{2}(x_{t,2})=1$, and any linear combination $b_{1}x_{t,1}+b_{2}x_{t,2}$ with $b_{1}^{2}+b_{2}^{2}=1$ has unit variance. By eliminating the variance of the learned mode as a discriminating feature within the loss functions Eq. .1 and Eq. .1, learning is performed exclusively on the basis of autocorrelation, and the linear TAE can correctly learn the slowest mode by minimizing the reconstruction loss. A rigorous proof that linear TAEs employing whitened features is equivalent to TICA for reversible process and can correctly identify the slowest mode is presented in Ref. Wehmeyer and Noé (2018).

.2 Nonlinear TAEs cannot equalize the variance explained within the input features and so cannot be assured to learn the slowest mode

We now proceed to perform a similar analysis for nonlinear TAEs. To aid in our discussion, we define $\bar{x}(z)=\mathop{\mathbb{E}_{x}}\left[x|z\right]$ and introduce the “variance explained” $\sigma^{2}(\bar{x}(z))$ by the latent variable $z$ as,

[TABLE]

which measures how much variance in the feature space can be explained by $z$. The idea of this concept is as follows. In a standard (non-time-lagged) autoencoder, the optimal reconstruction $\tilde{x}=D(z)$ is $D(z)=\mathop{\mathbb{E}_{x}}\left[x|z\right]$ and the variance of the outputs is $\sigma^{2}(\mathop{\mathbb{E}_{x}}\left[x|z\right])$, which explains part of the variance in the feature space while leaving $\left(\sigma^{2}(x)-\sigma^{2}(\mathop{\mathbb{E}_{x}}\left[x|z\right])\right)$ unexplained. In analogy to the linear case, it may be shown that for a long trajectory generated by a stationary process with finite states, the optimal loss for nonlinear TAE is approximately

[TABLE]

where $G(z)$ is the nonlinear generalization of autocorrelation $A(z)$ and $\sigma^{2}(x)$ is the total variance of the input features. Here $\sigma^{2}(\bar{x}(z))(1-G^{2}(z))$ is the “propagation loss” and $\left(\sigma^{2}(x)-\sigma^{2}(\bar{x}(z))\right)$ is the “irreducible capacity loss” for the nonlinear case. The details of the proof and related concepts can be found in Section .1 of the Appendix.

From Eq. 14, we see that both variance explained and the autocorrelation contribute to the TAE loss as in the linear case. If a faster mode has much larger variance explained than the true slowest mode, it is possible that nonlinear TAE would learn this faster mode in the latent encoding to minimize the loss, and in this case the nonlinear TAE fails to learn the correct slowest mode.

We demonstrate this idea in the context of “Washington beltway” potential, which consists two circular potential valleys at 0 $k_{B}T$, separated by a circular barrier of 4 $k_{B}T$ (Fig. 2). The expression for the potential is given by

[TABLE]

where $r_{1}=0.7,r_{2}=0.9,dr=0.05$. To simulate a particle moving in this potential, we conduct a Markov state model (MSM) simulation to generate a trajectory following the procedure detailed in Ref. Chen et al. (2019). Specifically, we split $(r,\theta)\in[0.6,1]\times[0,2\pi]$ into 20-by-200 evenly spaced bins, and run a 5,000,000 step MSM simulation according to transition probabilities from bin $i$ to bin $j$,

[TABLE]

where $C_{i}$ is the normalization factor that ensures $\sum_{j}p_{ij}=1$.

Given the analytical expression for the potential, the slow modes are analytically calculable by computing eigenvectors of the MSM transition matrix. The first six leading slow modes with their implied timescales are presented in Fig. 3. As expected, the slowest mode corresponds to transitions over the potential barrier radial direction ($r$-direction). Transitions around the circular potential in the polar direction ($\theta$-direction) appear as a degenerate pair within the second and third slowest modes with an order of magnitude faster timescale. Higher-order modes correspond to mixed $r$-$\theta$ transitions and appear as degenerate pairs due to the circular symmetry in $\theta$.

We now analyze the MSM trajectory using a nonlinear TAE with a [2-50-50-1-50-50-2] architecture and tanh activation functions for encoding and decoding layers and linear activation functions for input/output/latent layers. We supply to the input and output layers $(x(t),y(t))$ and $(x(t+\tau),y(t+\tau))$ pairs, respectively, employing a lag time of $\tau$ = 3000 steps. The data are naturally mean-free due to the topology of the potential landscape and are pre-whitened to normalize their variance. In Fig. 4a, we show the contour plot of the latent variable for the system. Clearly, the TAE did not correctly identify theoretical slowest mode, which should be along $r$ direction. Instead, since $\theta$-direction has much larger variance explained than $r$-direction, the TAE is biased towards identifying $\theta$ as slower mode in order to minimize time-lagged reconstruction loss. Therefore the learned latent variable is actually the mixture of $r$ and $\theta$ with most contribution coming from $\theta$.

To be quantitative, we compare the TAE loss employing $r$ and $\theta$ as the latent variable. Due to symmetry, when $z\sim r$ (denoting that $z$ is a bijection of $r$), the average configuration given $r$ is $\bar{x}(z\sim r)=0$, so $\sigma^{2}(\bar{x}(z\sim r))=0$ and $d_{\tau}(z\sim r)=\sigma^{2}(x)$ from Eq. 14. Conversely if $z\sim\theta$, the average configuration given $\theta$ is $\bar{x}(z\sim\theta)=(\mathop{\mathbb{E}}\left[r\right]\cos\theta,\mathop{\mathbb{E}}\left[r\right]\sin\theta)$, so $\sigma^{2}(\bar{x}(z\sim\theta))=\mathop{\mathbb{E}}\left[r\right]^{2}$ and $d_{\tau}(z\sim\theta)=\sigma^{2}(x)-\mathop{\mathbb{E}}\left[r\right]^{2}G^{2}(z\sim\theta)$. For $\mathop{\mathbb{E}}\left[r\right]^{2}G^{2}(z\sim\theta)>0$, learning $z\sim\theta$ as the slow mode results in a lower time-lagged reconstruction loss than learning the true slow mode $z\sim r$. We numerically estimate from the simulation trajectory $G(z\sim\theta)\approx 0.535,\sigma^{2}(x)=2$, and $\mathop{\mathbb{E}}\left[r\right]\approx 1.4$, from which we compute $d_{\tau}(z\sim r)=2$ and then the estimated loss $d_{\tau}(z\sim\theta)\approx 1.44$. The actual training loss is computed to be $d_{\tau}\approx 1.48$, showing that the nonlinear TAE approximately learns $\theta$ as the slowest mode. As a side note, if we use an optimal encoding corresponding to a bijective encoding of the feature space (see last part of Section .1 of the Appendix for details), the minimal possible loss is $d_{\tau}\approx 1.42$, which implies that $\theta$ is very close to the optimal encoding.

Can we somehow transform the input features $(x,y)$ to equalize the variance explained? This operation would serve the same purpose as the whitening transformation in the linear case to eliminate the variance explained as a discriminating factor in the time-lagged reconstruction loss and force the TAE to identify the slow mode based on (generalized) autocorrelation alone. We first note that such a transform is not always possible even if we know the theoretical slow modes. For instance, in our Washington beltway potential example, we cannot equalize $\sigma^{2}(\bar{x}(z\sim r))$ and $\sigma^{2}(\bar{x}(z\sim\theta))$ since the intrinsic circular symmetry assures that $\sigma^{2}(\bar{x}(z\sim r))=0$. What about other systems with slow modes that are possible to be theoretically whitened? We note that for all but the most trivial systems we do not know the nonlinear slow modes a priori and so there is no way to equalize the variance explained. Even if we are able to identify putative slow modes using nonlinear dimensionality reduction techniques, if these modes are identified on the basis of anything other than slowness (e.g., variance) then equalization of the variance explained can still lead to incorrect results. If they are identified only on the basis of slowness then we have already obtained the slowest mode and there is no need to apply nonlinear TAE. In short, there is no general procedure to correctly equalize the variance explained within the input features and therefore no way to guarantee that nonlinear TAEs will correctly discover the slowest mode.

.3 State-free reversible VAMPnets (SRVs) correctly identify nonlinear slow modes since the explained variance does not appear in the loss function

We now demonstrate that in contrast to TAEs, state-free reversible VAMPnets (SRVs) can correctly identify the slowest mode using the same input features. Given a trajectory $\{x_{t}\}$, SRVs employ an artificial neural network to simultaneously discover the optimal nonlinear featurization of the input data $E(x_{t})=\{E_{j}(x_{t})\}_{j=1}^{d}$, where $E_{j}(x_{t})$ is the $j^{\text{th}}$ component of the encoder output, and the linear combination of these features $z_{i}=\sum_{j=1}^{d}s_{ij}E_{j}(x_{t})$ that maximizes the squared sum of autocorrelation of the slow modes. Full details of SRVs are presented in Ref. Chen et al. (2019). A schematic of an SRV with $d$=1 corresponding to a 1D latent space embedding is presented in Fig. 5. The corresponding loss function for the 1D SRV is,

[TABLE]

We apply a SRV model with [2-50-50-1] architecture and tanh activation functions for all layers except input/output layers (an analogous design to the TAE above) to the whitened MSM trajectory on the Washington beltway potential. In Fig. 4b, we show the contour plot of the learned SRV mode that correctly identifies transitions in $r$ as the slowest mode of the system. The reason why the SRV correctly identifies slowest mode is that the loss function is exactly equivalent to maximizing the autocorrelation of the learned mode with no contribution from variance explained. Accordingly, we generally recommend to use a SRV rather than a TAE for estimation of the slowest mode. Moreover, we note that when it comes to higher-order slow modes discovery through multidimensional latent variables, it is not recommended to use TAEs since there is no constraint that different components of the latent variable are orthogonal to each other, and it is therefore possible that each component becomes mixture of many slow modes. In SRVs, however, the orthogonality constraints are satisfied naturally in the variational optimization procedure and it can discover a hierarchy of slow orthogonal modes Chen et al. (2019).

.4 The “slow” mode learned by TAEs can be controlled by feature engineering

Following the ideas developed above, we now show that we can design features of a system to mislead nonlinear TAEs to learn a desired mode as the slowest mode simply by assuring that the target mode has much larger variance explained than the true slowest mode. To show this, we employ molecular dynamics simulations of alanine dipeptide as 22-atom peptide that serves as the “fruit fly” for testing new numerical methods in molecular systems (Fig. 6).

It is well known from extensive prior study that the slowest mode for alanine dipeptide corresponds to transitions in the backbone dihedral angle $\phi$ Chen et al. (2019); Trendelkamp-Schroer et al. (2015). Is it possible to design a feature set such that a TAE is misled into learning the backbone dihedral angle $\psi$ as the slowest mode? To do so, consider 2D features $(x_{1},x_{2})$ given by,

[TABLE]

where $\Delta r<<r_{0}$. The idea is to map $(\phi,\psi)$ to a ring such that $\psi$ encodes polar angle direction while $\phi$ encodes radial direction such that $\psi$ has much larger variance explained than $\phi$. Scatter plots of the engineered feature space colored by two dihedral angles $\phi$ and $\psi$ are presented in Fig. 7.

We run molecular dynamics simulation for alanine dipeptide in explicit solvent to generate a trajectory of 2,000 ns saving frames every 2 ps to generate a 1,000,000-frame trajectory. Then we use a TAE with [2-50-50-1-50-50-2] architecture and tanh activation functions for encoding and decoding layers and linear activation functions for input/output/latent layers to learn over the mean-free whitened trajectory data with features described above. We show the Ramachandran plot colored by the slowest mode discovered by TAE in Fig. 8a, which verifies that we successfully misled the TAE to learn $\psi$ as the slowest mode. Again, the SRV employing an analogous architecture has no difficulty in correctly identifying the slowest mode (Fig. 8b), which is consistent with the ground truth results of a state-of-the-art Markov state model trained in Ref. Chen et al. (2019) (Fig. 8c).

Similarly if we swap $\phi$ and $\psi$ and define our features according to,

[TABLE]

then the TAE learns $\phi$ as the slowest mode (Fig. 9a), which is closer to the ground truth MSM results (Fig. 9c), and again the SRV has no difficulty identifying slowest mode with these features (Fig. 9b)

.5 The 1D TAE structure can be modified to explicitly equalize the variance explained within the latent space to render it equivalent to a 1D SRV

To resolve the variance explained issue for nonlinear TAE, we can modify its structure. Instead of employing an autoencoding architecture to minimize the time-lagged reconstruction loss, we can instead employ an encoder-only architecture and optimize the time-lagged loss for the “whitened” encoder output. This modification in architecture replaces the standard TAE loss given by Eq. .1 with the modified TAE loss given by,

[TABLE]

where $E(x_{t})$ is the 1D encoder output for input $x_{t}$ at time $t$, the variance $\sigma^{2}(E(x_{t}))$ in the denominator is used to correctly “whiten” the learned slow mode such that the variance (or “variance explained”) does not play a role in learning of the slowest mode.

How does the modified TAE architecture and loss function relate to the 1D SRV? Following a similar derivation to Eq. .1, we have,

[TABLE]

where $A(E(x_{t}))$ is the autocorrelation for $E(x_{t})$. Therefore Eq. 20 becomes

[TABLE]

which, up to a trivial affine transformation, is equivalent to the SRV loss given by Eq. 17. Accordingly, the modified nonlinear TAE with loss given by Eq. 22 is equivalent to a 1D SRV and can correctly identify slow modes without the misleading influence of variance explained.

.6 Variational dynamics encoders (VDEs) learn mixtures of the slowest mode and the maximum variance mode

Variational dynamics encoders (VDEs) employ a variational autoencoder architecture with a loss function comprising both reconstruction loss for inputs/outputs $x$ and the autocorrelation loss for the 1D latent variable $z$. The loss function of VDE can be written as Hernández et al. (2018),

[TABLE]

where $D(z_{t})$ is the reconstructed output, $\mathop{\mathbb{E}}\left[\left\lVert D(z_{t})-x_{t+\tau}\right\rVert^{2}\right]$ is the time-lagged reconstruction loss, $(-A(z))$ is the autocorrelation loss for $z$, and $\lambda$ is a linear mixing parameter. If we ignore the Kullback-Leibler divergence term $L_{KL}$, which measures the similarity of the encoded probability distribution of $z$ to a Gaussian distribution, can be considered a form of regularization, and goes to zero in the case of well-trained variational autoencoders Doersch (2016), then the loss is exactly the mixture of 1D SRV loss and 1D TAE loss. Since SRVs learn the slowest mode, and TAEs learn a mixture of the slowest mode and the maximum variance mode, in general VDEs also learn a mixture of the slowest mode and the maximum variance mode. As was the case for TAEs, in any applications where we aim to find the slow modes it is not recommended to include terms associated with the explained variance (here within the reconstruction loss) and it is therefore not recommended to use VDEs for this goal.

Methods

The TAE and SRV neural networks were constructed in Python using the Keras Chollet deep learning libraries and training performed on an NVIDIA GeForce GTX 1080 GPU card. The MSM simulations of particle motion over the “Washington beltway” potential were conducted in Python. Simulations of alanine dipeptide in water were conducted using the OpenMM 7.3 simulation suite Eastman et al. (2012, 2017) employing the Amber99sb-ILDN forcefield for the biomolecule Lindorff-Larsen et al. (2010) and TIP3P forcefield for water Jorgensen et al. (1983). The temperature and pressure were maintained at $T$ = 300 K and $P$ = 1 bar using Andersen thermostat Andersen (1980) and $P$ = 1 atm using a Monte-Carlo barostatChow and Ferguson (1995); Åqvist et al. (2004). Lennard-Jones interactions were smoothly switched to zero at a cuttoff of 1.4 nm cutoff, and Coulombic interactions were treated by particle-mesh Ewald Essmann et al. (1995) with a real space cutoff of 1.4 nm and a reciprocal space grid spacing of 0.12 nm.

Conclusions

In this work, we present a theoretical analysis of the capabilities and limitations of TAEs in slow mode discovery. We show that linear TAEs correctly learn the slowest linear mode with whitened features, while nonlinear TAEs cannot be assured of discovering the slowest nonlinear mode since it is not in general possible to perform an equivalent nonlinear “whitening” of the input features to equalize their variance explained. We prove the theoretical bounds for time-lagged reconstruction loss for nonlinear TAE and demonstrate that a faster nonlinear mode could be erroneously identified as the slowest mode if it has significantly higher variance explained. We validate our theoretical analysis in applications to a 2D “Washington beltway” potential and in molecular simulations of alanine dipeptide. We also show how the 1D nonlinear TAE network structure can be modified to become equivalent to the 1D SRV to remove the misleading influence of variance explained and permit variable discrimination exclusively on the basis of autocorrelation. We also show that 1D VDEs mix the loss functions of TAEs and SRVs and therefore also suffer from the misleading influence of variance explained. Accordingly, SRVs serve as a more appropriate tool than TAEs or VDEs for the discovery of slow modes. If variance explained of the modes is also of interest, then TAEs or VDEs can prove useful.

Acknowledgments

This material is based upon work supported by the National Science Foundation under Grant No. CHE-1841805. H.S. acknowledges support from the Molecular Software Sciences Institute (MolSSI) Software Fellows program (NSF grant ACI-1547580)Krylov et al. (2018); Wilkins-Diehr and Crawford (2018).

Appendix

.1 Theoretical bounds for time-lagged autoencoder reconstruction loss

We present a derivation of the theoretical bounds on the time-lagged autoencoder (TAE) reconstruction loss for a stationary process with finite states. We show that this leads to the approximate expression Eq. 14 in the main text, and that the approximation can be made arbitrarily accurate by employing sufficiently large neural network architectures. Let $x_{t}$ be the mean-free featurization of the system at time $t$ in a finite trajectory that therefore comprises a finite number of states. The process is stationary and the trajectory is long enough such that for any $\tau\geq 0$,

[TABLE]

where $\sigma^{2}(x)$ is the total variance of featurized configurations.

Let $z_{t}=z(x_{t})$ be the encoded embedding of $x_{t}$, where $z(\cdot)$ represents the encoding mapping. Let $f_{t}=f(z_{t})$ be the decoded output of $z_{t}$, where $f(\cdot)$ represents the decoding mapping. The time-lagged reconstruction loss $d_{\tau}(z)$ can be written as,

[TABLE]

We define the expected evolution $\tilde{D}_{\tau}(z)$ after lag time $\tau$ given latent space encoding $z$ as,

[TABLE]

which can be viewed as the “optimal reconstructed output” given $z$. Note that when $\tau=0$, $\bar{x}(z)=\tilde{D}_{0}(z)$ denotes the average system featurization corresponding to a particular latent-space encoding $z$.

To quantify how much variance information of $x$ is captured by $z$, we define the variance explained,

[TABLE]

which measures the variance of the feature average given the encoding $z$ and is consistent with our definition Eq. 13 in the main text.

Now Eq. 25 becomes,

[TABLE]

Since we have finite number of possible $x_{t}$ values, the number of $z_{t}$ values should also be finite. Therefore we can split the summation over the N-frame trajectory into two parts (sum-splitting trick): first summing over all frames with the same $z$ value, then summing over all $z$ values. Considering the third term in Eq. .1,

[TABLE]

Therefore Eq. .1 becomes,

[TABLE]

We can apply the same sum-splitting idea to the third term of Eq. .1, where $N$ is total number of frames, and $N(z)$ is the number of frames with encoding equal to $z$,

[TABLE]

Therefore Eq. .1 becomes,

[TABLE]

If we define the “generalized autocorrelation” $G(z)$ as,

[TABLE]

then the lower bound of the TAE reconstruction loss is given by,

[TABLE]

Now we consider how good our lower bound is. Due to the universal approximation theorem Hassoun (1995); Chen and Chen (1995), for any $\epsilon>0$, there exists a finite size decoder neural network $f$ such that the following inequality holds,

[TABLE]

which means $f$ can be made arbitrarily close to the “optimal reconstructed output” as given in Eq. 26 by employing a sufficiently large neural network.

Therefore Eq. .1 becomes

[TABLE]

This indicates that the lower bound is actually quite tight, and it is a relatively good approximation of the TAE loss, yielding,

[TABLE]

corresponding to Eq. 14 in the main text.

What are the optimal encoding $z$ and the corresponding minimal possible loss for the TAE? From Eq. .1 we see that if $f(z_{t})=\tilde{D}_{\tau}(z_{t})$, the optimal encoding $z$ should maximize $\mathop{\mathbb{E}_{z}}\left[\tilde{D}_{\tau}(z)^{2}\right]$. We define a finite set $S^{z}(z^{\prime})=\{x|z(x)=z^{\prime}\}$ that includes all configurations mapped to the same encoding value $z^{\prime}$ under encoding mapping $z(\cdot)$. If $z^{(1)}$ is a bijective encoding of $x$, then $\tilde{D}_{\tau}(z^{(1)})$ is,

[TABLE]

where $N(x(z^{(1)}))$ is the number of frames with configuration $x_{t}=x(z^{(1)})$ and $x(z^{(1)})$ is the configuration corresponding to $z^{(1)}$. For any encoding $z$, we have,

[TABLE]

where in the fourth line, we use Jensen’s inequality, considering that,

[TABLE]

Therefore,

[TABLE]

This means that the optimal encoding that minimizes time-lagged reconstruction loss should be a bijective encoding of configurations.

.2 Interpretation of $G(z)$ as generalized autocorrelation

We now demonstrate why $G(z)$ can be interpreted as a nonlinear generalization of the linear autocorrelation $A(z)$ that appears in the loss function for linear TAEs (Eq. 1). Consider a linear TAE applied on trajectory $\{x_{t}\}$ with independent components. If $z_{t}$ is a linear transformation of the $k^{\text{th}}$ component of input features $x_{t,k}$, then $f(z_{t})$ is also a linear transformation of $x_{t,k}$. Without loss of generality and considering the independence among different components and that $z_{t}$ only includes information of $k^{\text{th}}$ component, all components of optimal output $f$ except $k^{\text{th}}$ component should be equal to 0, therefore we can write the encoding and decoding mapping as,

[TABLE]

Similarly the $s^{\text{th}}$ component of $\bar{x}(z_{t})$ and $\tilde{D}_{\tau}(z_{t})$ should be given by,

[TABLE]

So all terms in $\min_{f}\mathop{\mathbb{E}}\left[\left\lVert\tilde{D}_{\tau}(z_{t})-f(z_{t})\right\rVert^{2}\right]+\left(\sigma^{2}(\bar{x}(z))-\mathop{\mathbb{E}_{z}}\left[\tilde{D}_{\tau}(z)^{2}\right]\right)$ of Eq. 33 have non-zero values only in their $k^{\text{th}}$ component.

Now Eq. 33 becomes,

[TABLE]

where in the fourth line we use sum-splitting trick to simplify the second term in the numerator of the fraction,

[TABLE]

where $N(x_{t,k})$ is number of frames with $k^{\text{th}}$ component equal to $x_{t,k}$, $\sum_{x_{t,k}}$ denotes summation over all possible $x_{t,k}$ values, and $\sum_{x_{t^{\prime},k}=x_{t,k}}$ denotes summation over all frames with $k^{\text{th}}$ component equal to $x_{t,k}$.

The optimal $c_{0}$ and $c_{1}$ follow from minimization with respect to these parameters and satisfy,

[TABLE]

where $A(x_{t,k})$ is the traditional autocorrelation of component $x_{t,k}$.

So Eq. .2 becomes,

[TABLE]

We see that in the linear case with independent components, the generalized autocorrelation $G(z)$ reduces to the standard linear autocorrelation $A(z)$.

In the nonlinear case, the minimal possible loss in Eq. .1 is obtained if $f(z_{t})=\tilde{D}_{\tau}(z_{t})$, wherein $G(z)$ becomes,

[TABLE]

We note the following two attractive properties of $G(z)$. First, $G(h(z))=G(z)$ if $h$ is a bijection. This is important for nonlinear TAEs since under the transformation $z\to h(z),f\to f\circ h^{-1}$ the TAE loss is invariant and $G(\cdot)$ is invariant, whereas the autocorrelation $A(\cdot)$ is not invariant. Second, $G(z)$ can be naturally applied to multi-dimensional $z$ where it is difficult to define the traditional autocorrelation $A(z)$.