Inexact iterative numerical linear algebra for neural network-based
  spectral estimation and rare-event prediction
Understanding dynamics in complex systems is challenging because there are
many degrees of freedom, and those that are most important for describing
events of interest are often not obvious. The leading eigenfunctions of the
transition operator are useful for visualization, and they can provide an
efficient basis for computing statistics such as the likelihood and average
time of events (predictions). Here we develop inexact iterative linear algebra
methods for computing these eigenfunctions (spectral estimation) and making
predictions from a data set of short trajectories sampled at finite intervals.
We demonstrate the methods on a low-dimensional model that facilitates
visualization and a high-dimensional model of a biomolecular system.
Implications for the prediction problem in reinforcement learning are
discussed.
John Strahan, Spencer C. Guo, Chatipat Lorpaiboon, Aaron R. Dinner, Jonathan Weare
Inexact iterative numerical linear algebra for neural network-based spectral estimation and rare-event prediction

###### Abstract

Understanding dynamics in complex systems is challenging because there are many degrees of freedom, and those that are most important for describing events of interest are often not obvious. The leading eigenfunctions of the transition operator are useful for visualization, and they can provide an efficient basis for computing statistics such as the likelihood and average time of events (predictions). Here we develop inexact iterative linear algebra methods for computing these eigenfunctions (spectral estimation) and making predictions from a data set of short trajectories sampled at finite intervals. We demonstrate the methods on a low-dimensional model that facilitates visualization and a high-dimensional model of a biomolecular system. Implications for the prediction problem in reinforcement learning are discussed.

## 1 Introduction

Modern observational, experimental, and computational approaches often yield high-dimensional time series data (trajectories) for complex systems. In principle, these trajectories contain rich information about dynamics and, in particular, the infrequent events that are often most consequential. In practice, however, high-dimensional trajectory data are often difficult to parse for useful insight. The need for more efficient statistical analysis tools for trajectory data is critical, especially when the goal is to understand rare-events that may not be well represented in the data.

We consider dynamics that can be treated as Markov processes. A common starting point for statistical analyses of Markov processes is the transition operator, which describes the evolution of function expectations. The eigenfunctions of the transition operator characterize the most slowly decorrelating features (modes) of the system [1, 2, 3, 4, 5]. These can be used for dimensionality reduction to obtain a qualitative understanding of the dynamics [6, 7], or they can be used as the starting point for further computations [8, 9, 10]. Similarly, prediction functions, which provide information about the likelihood and timing of future events as a function of the current state, are defined through linear equations of the transition operator[11, 10].

A straightforward numerical approach to obtaining these functions is to convert the transition operator to a matrix by projecting onto a finite basis for Galerkin approximation[1, 2, 10, 11, 12, 13, 14, 15]. The performance of such a linear approximation depends on the choice of basis [11, 10, 15], and previous work often resorts to a set of indicator functions on a partition of the state space (resulting in a Markov state model or MSM [14]) for lack of a better choice. While Galerkin approximation has yielded many insights [16, 17], the limited expressivity of the basis expansion has stimulated interest in (nonlinear) alternatives.

In particular, artificial neural networks can be harnessed to learn eigenfunctions of the transition operator and prediction functions from data [18, 19, 20, 21, 22, 23, 24, 25]. However, existing approaches based on neural networks suffer from various drawbacks. As discussed in Ref. [5], their performance can often be very sensitive to hyperparameters, requiring extensive tuning and varying with random initialization. Many use loss functions that are estimated against the stationary distribution [25, 26, 27, 28, 29, 30], so that metastable states contribute most heavily, which negatively impacts performance [24, 30]. Assumptions about the dynamics (e.g., microscopic reversibility) limit applicability. In Ref. [24] we introduced an approach that overcomes the issues above, but it uses multiple trajectories from each initial condition; this limits the approach to analysis of simulations and moreover requires specially prepared data sets.

The need to compute prediction functions from observed trajectory data also arises in reinforcement
learning. There the goal is to optimize an expected future reward (the prediction function) over a policy (a Markov process). For a fixed Markov process, the prediction problem in reinforcement learning is often solved by temporal difference (TD) methods, which allow the use of arbitrary ensembles of trajectories without knowledge of the details of the underlying dynamics[31]. TD methods have a close relationship with an inexact form of power iteration, which, as we describe, can perform poorly on rare-event related problems.

Motivated by this relationship, as well as by an inexact power iteration scheme previously proposed for approximating the stationary probability distribution of a Markov process using trajectory data [32], we propose a computational framework for spectral estimation and rare-event prediction based on inexact iterative numerical linear algebra. Our framework includes an inexact Richardson iteration for the prediction problem, as well as an extension to inexact subspace iteration for the prediction and spectral estimation problems. The theoretical properties of exact subspace iteration suggest that eigenfunctions outside the span of the approximation will contribute significantly to the error of our inexact iterative schemes[33]. Consistent with this prediction, we demonstrate that learning additional eigenvalues and eigenfunctions simultaneously through inexact subspace iteration accelerates convergence dramatically relative to inexact Richardson and power iteration in the context of rare events. While we assume the dynamics can be modeled by a Markov process, we do not require knowledge of their form or a specific underlying model. The method shares a number of further advantages with the approach discussed in Ref. [24] without the need for multiple trajectories from each initial condition in the data set. This opens the door to treating a wide range of observational, experimental, and computational data sets.

The remainder of the paper is organized as follows. In Section 2, we describe the quantities that we seek to compute in terms of linear operators. In Sections 3 and 4, we introduce an inexact subspace iteration algorithm that we use to solve for these quantities. Section 5 illustrates how the loss function can be tailored to the known properties of the desired quantity. Section 6 summarizes the two test systems that we use to illustrate our methods: a two-dimensional potential, for which we can compute accurate reference solutions, and a molecular example that is high-dimensional but still sufficiently tractable that statistics for comparison can be computed from long trajectories. In Section 7, we explain the details of the invariant subspace iteration and then demonstrate its application to our two examples. Lastly, Section 8 details how the subspace iteration can be modified to compute prediction functions and compares the effect of different loss functions, as well as the convergence properties of power iteration and subspace iteration. We conclude with implications for reinforcement learning.

## 2 Spectral estimation and the prediction problem

We have two primary applications in mind in this article. First, we would like to estimate the _dominant eigenfunctions and eigenvalues_ of the transition operator \(\mathcal{T}^{t}\) for a Markov process \(X^{t}\in\mathbb{R}^{d}\), defined as

\[\mathcal{T}^{t}f(x)=\mathbb{E}_{x}\left[f(X^{t})\right], \tag{1}\]

where \(f\) is an arbitrary real-valued function and the subscript \(x\) indicates the initial condition \(X^{0}=x\). The transition operator encodes how expectations of functions evolve in time. The right eigenfunctions of \(\mathcal{T}^{t}\) with the largest eigenvalues characterize the most slowly decorrelating features (modes) of the Markov process [1, 2, 4, 5].

Our second application is to compute _prediction functions_ of the general form

\[u(x)=\mathbb{E}_{x}\left[\Psi(X^{T})+\sum_{t=0}^{T-1}\Gamma(X^{t})\right], \tag{2}\]

where \(T\) is the first time \(X^{t}\notin D\) from some domain \(D\), and \(\Psi\) and \(\Gamma\) are functions associated with the escape event (rewards in reinforcement learning). Prototypical examples of prediction functions that appear in our numerical results are the mean first passage time (MFPT) from each \(x\):

\[m(x)=\mathbb{E}_{x}\left[T\right] \tag{3}\]

and the committor:

\[q(x)=\mathbb{E}_{x}\left[\mathbbm{1}_{B}(X^{T})\right]=\mathbb{P}_{x}\left[X ^{T}\in B\right], \tag{4}\]

where \(A\) and \(B\) are disjoint sets ("reactant" and "product" states) and \(D=(A\cup B)^{\mathsf{c}}\). The MFPT is important for estimating rates of transitions between regions of state space, while the committor can serve as a reaction coordinate[29, 34, 35, 36] and
as a key ingredient in transition path theory statistics [24; 37; 38]. For any \(\tau>0\), the prediction function \(u(x)\) satisfies the linear equation

\[\left(\mathcal{I}-\mathcal{S}^{\tau}\right)u(x)=\mathbb{E}_{x}\left[\sum_{t=0}^ {\left(\tau\wedge T\right)-1}\Gamma(X^{t})\right] \tag{5}\]

for \(x\in D\), with boundary condition

\[u(x)=\Psi(x) \tag{6}\]

for \(x\notin D\). In (5), \(\mathcal{I}\) is the identity operator and

\[\mathcal{S}^{t}f(x)=\mathbb{E}_{x}\left[f(X^{t\wedge T})\right] \tag{7}\]

is the "stopped" transition operator [10]. We use the notation \(\tau\wedge T=\min\{\tau,T\}\).

Our specific goal is to solve both the eigenproblem and the prediction problem for \(X^{t}\) in high dimensions and without direct access to a model governing its evolution. Instead, we have access to trajectories of \(X^{t}\) of a fixed, finite duration \(\tau\). A natural and generally applicable approach to finding an approximate solution to the prediction problem is to attempt to minimize the residual of (5) over parameters \(\theta\) of a neural network \(u_{\theta}(x)\) representing \(u(x)\). For example, we recently suggested an algorithm that minimizes the residual norm

\[\frac{1}{2}\Big{\|}\left(\mathcal{I}-\mathcal{S}^{\tau}\right)u_{\theta}-r \Big{\|}_{\mu}^{2}, \tag{8}\]

where \(r(x)\) is the right-hand side of (5) and \(\mu\) is an arbitrary distribution of initial conditions \(X^{0}\) (boundary conditions were enforced by an additional term) [24]. The gradient of the residual norm in (8) with respect to neural-network parameters \(\theta\) can be written

\[\left\langle\left(\mathcal{I}-\mathcal{S}^{\tau}\right)u_{\theta}-r,\nabla_{ \theta}u_{\theta}\right\rangle_{\mu}-\left\langle\left(\mathcal{I}-\mathcal{S }^{\tau}\right)u_{\theta}-r,\mathcal{S}^{\tau}\nabla_{\theta}u_{\theta}\right \rangle_{\mu} \tag{9}\]

where \(\left\langle f,g\right\rangle_{\mu}=\int f(x)g(x)\,\mu(dx).\) Given a data set of initial conditions \(\left\{X_{j}^{0}\right\}_{j=1}^{n}\) drawn from \(\mu\) and a single sample trajectory \(\left\{X_{j}^{t}\right\}_{t=0}^{\tau}\) of \(X^{t}\) for each \(X_{j}^{0}\), the first term in the gradient (9) can be approximated without bias as

\[\left\langle\left(\mathcal{I}-\mathcal{S}^{\tau}\right)u_{\theta}-r,\nabla_{ \theta}u_{\theta}\right\rangle_{\mu}\approx\frac{1}{n}\sum_{j=1}^{n}\left(u_{ \theta}(X_{j}^{0})-u_{\theta}(X_{j}^{\tau\wedge T_{j}})-\sum_{t=0}^{(\tau \wedge T_{j})-1}\Gamma(X_{j}^{t})\right)\nabla_{\theta}u_{\theta}(X_{j}^{0}) \tag{10}\]

where \(T_{j}\) is the first time \(X_{j}^{t}\notin D\).

In contrast, unbiased estimation of the second term in (9) requires access to two independent trajectories of \(X^{t}\) for each sample initial condition since it is quadratic in \(\mathcal{S}^{\tau}\)[24; 31]. In the reinforcement learning community, TD methods were developed for the purpose of minimizing a loss of a very similar form to (8), and they perform a "semigradient" descent by following only the first term in (9) [31]. As in the semigradient approximation, the algorithms proposed in this paper only require access to inner products of the form \(\left\langle f,\mathcal{A}g\right\rangle_{\mu}\) for an operator \(\mathcal{A}\) related to \(\mathcal{T}^{\tau}\) or \(\mathcal{S}^{\tau}\), and they avoid terms that are nonlinear in \(\mathcal{A}\). As we explain, such inner products can be estimated using trajectory data. However, rather than attempting to minimize the residual directly by an approximate gradient descent, we view the eigenproblem and prediction problems through the lens of iterative numerical linear algebra schemes.

## 3 An Inexact Power Iteration

To motivate the iterative numerical linear algebra point of view, observe that the first term in (9) is the exact gradient with respect to \(\theta^{\prime}\) of the loss

\[\frac{1}{2}\big{\|}u_{\theta^{\prime}}-\mathcal{S}^{\tau}u_{\theta}-r\big{\|}_ {\mu}^{2}, \tag{11}\]

evaluated at \(\theta^{\prime}=\theta\). The result of many steps of gradient descent (later, "inner iterations") on this loss with \(\theta\) held fixed can then be viewed as an inexact Richardson iteration for (5), resulting in a sequence

\[u_{\theta^{\prime+1}}\approx\mathcal{S}^{\tau}u_{\theta^{\prime}}+r, \tag{12}\]
where, for each \(s\), \(u_{\theta^{s}}\) is a parametrized neural-network approximation of the solution to (5). To unify our discussion of the prediction and spectral estimation problems, it is helpful to observe that (12) can be recast as an equivalent inexact power iteration:

\[\bar{u}_{\theta^{s+1}}\approx\mathcal{A}\bar{u}_{\theta^{s}} \tag{13}\]

where

\[\bar{u}_{\theta}=\left(\begin{array}{c}1\\ u_{\theta}\end{array}\right)\quad\text{and}\quad\mathcal{A}=\left[\begin{array} []{cc}1&0\\ r&\mathcal{S}^{\tau}\end{array}\right]. \tag{14}\]

Note that \(\left(1,u\right)^{\top}\) is the dominant eigenfunction of \(\mathcal{A}\) and has eigenvalue equal to \(1\).

Ref. [32] introduced an inexact power iteration to compute the stationary probability measure of \(\mathcal{T}^{\tau}\), i.e., its dominant left eigenfunction. As those authors note, an inexact power update such as (13) can be enforced using a variety of loss functions. In our setting, the \(L^{2}_{\mu}\) norm in (11) can be replaced by any other measure of the difference between \(u_{\theta^{\prime}}\) and \(\mathcal{S}^{\tau}u_{\theta}+r\), as long as an unbiased estimator of its gradient with respect to \(\theta^{\prime}\) is available. This flexibility is discussed in more detail in Section 5, and we exploit it in applications presented later in this article. For now, we focus instead on another important implication of this viewpoint: the flexibility in the form of the iteration itself.

We will see that when the spectral gap of \(\mathcal{A}\), the difference between its largest and second largest eigenvalues, is small, the inexact power iteration (or the equivalent Richardson iteration) described above fails to reach an accurate solution. The largest eigenvalue of \(\mathcal{A}\) in (14) is \(1\) and its next largest eigenvalue is the dominant eigenvalue of \(\mathcal{S}^{\tau}\). For rare-event problems the difference between these two eigenvalues is expected to be very small. Indeed, when \(X^{0}\) is drawn from the quasi-stationary distribution of \(X^{t}\) in \(D\), the logarithm of the largest eigenvalue of \(\mathcal{S}^{\tau}\) is exactly \(-\tau/\operatorname{\mathbb{E}}[T]\)[39]. Thus, when the mean exit time from \(D\) is large, we can expect the spectral gap of \(\mathcal{A}\) in (14) to be very small. In this case, classical convergence results tell us that exact power iteration converges slowly, with the largest contributions to the error coming from the eigenfunctions of \(\mathcal{S}^{\tau}\) with largest magnitude eigenvalues [33]. Iterative schemes that approximate multiple dominant eigenfunctions simultaneously, such as subspace iteration and Krylov methods, can converge much more rapidly[33]. In the next section, we introduce an inexact form of subspace iteration to address this shortcoming. Beyond dramatically improving performance on the prediction problem for rare-events when applied with \(\mathcal{A}\) in (14), it can also be applied with \(\mathcal{A}=\mathcal{T}^{\tau}\) to compute multiple dominant eigenfunctions and values of \(\mathcal{T}^{\tau}\) itself.

## 4 An Inexact Subspace Iteration

Our inexact subspace iteration for the \(k\) dominant eigenfunctions/values of \(\mathcal{A}\) proceeds as follows. Let \(\{\varphi_{\theta^{s}}^{a}\}_{a=1}^{k}\) be a sequence of \(k\) basis functions parametrized by \(\theta^{s}\) (these can be scalar or vector valued functions depending on the form of \(\mathcal{A}\)). We can represent this basis as the vector valued function

\[U_{\theta^{s}}=\left(\varphi_{\theta^{s}}^{1},\varphi_{\theta^{s}}^{2},\ldots, \varphi_{\theta^{s}}^{k}\right). \tag{15}\]

Then, we can obtain a new set of \(k\) basis functions by approximately applying \(\mathcal{A}\) to each of the components of \(U_{\theta^{s}}\):

\[U_{\theta^{s+1}}K^{s+1}\approx\mathcal{A}\,U_{\theta^{s}} \tag{16}\]

where \(K^{s+1}\) is an invertible, upper-triangular \(k\times k\) matrix that does not change the span of the components of \(U_{\theta^{s+1}}\) but is included to facilitate training. One way the approximate application of \(\mathcal{A}\) can be accomplished is by minimizing

\[\frac{1}{2}\sum_{a=1}^{k}\!\left\|\sum_{b=1}^{k}\varphi_{\theta}^{b}K_{ba}- \mathcal{A}\,\varphi_{\theta^{s}}^{a}\right\|_{\mu}^{2} \tag{17}\]

over \(\theta\) and \(K\) with \(\theta^{s}\) held fixed.

The eigenvalues and eigenfunctions of \(\mathcal{A}\) are then approximated by solving the finite-dimensional generalized eigenproblem

\[C^{\tau}W=C^{0}W\Lambda \tag{18}\]

where

\[C^{\tau}_{ab} =\langle\varphi_{\theta^{s}}^{a},\mathcal{A}\varphi_{\theta^{s}}^{ b}\rangle_{\mu} \tag{19}\] \[C^{0}_{ab} =\langle\varphi_{\theta^{s}}^{a},\varphi_{\theta^{s}}^{b}\rangle_{ \mu}, \tag{20}\]

each inner product is estimated using trajectory data, and \(W\) and \(\Lambda\) are \(k\times k\) matrices. The matrix \(\Lambda\) is diagonal, and the \(a\)-th eigenvalue \(\lambda_{a}\) of \(\mathcal{A}\) is approximated by \(\Lambda_{aa}\); the corresponding eigenfunction \(v_{a}\) is approximated by \(\sum_{b=1}^{k}W_{ab}\,\varphi_{\theta^{s}}^{b}\).

Even when sampling is not required to estimate the matrices in (19) and (20), the numerical rank of \(C^{\tau}\) becomes very small as the eigenfunctions become increasingly aligned with the single dominant eigenfunction. To overcome this issue, we apply an
orthogonalization step between iterations (or every few iterations). Just as the matrices \(C^{0}\) and \(C^{\tau}\) can be estimated using trajectory data, the orthogonalization procedure can also be implemented approximately using data.

Finally, in our experiments we find it advantageous to damp large parameter fluctuations during training by mixing the operator \(\mathcal{A}\) with a multiple of the identity, i.e., we perform our inexact subspace iteration on the operator \((1-\alpha_{s})\mathcal{I}+\alpha_{s}\mathcal{A}\) in place of \(\mathcal{A}\) itself. This new operator has the same eigenfunctions as \(\mathcal{A}\). In our experiments, decreasing the parameter \(\alpha_{s}\) as the number of iterations increases results in better generalization properties of the final solution and helps ensure convergence of the iteration. For our numerical experiments we use

\[\alpha_{s}=\begin{cases}1&s<\sigma\\ 1/\sqrt{s+1-\sigma}&s\geq\sigma\end{cases} \tag{21}\]

where \(\sigma\) is a user chosen hyperparameter that sets the number of iterations performed before damping begins.

The details, including estimators for all required inner products, in the case of the eigenproblem (\(\mathcal{A}=\mathcal{T}^{\tau}\)) are given in Section 7 and Algorithm 1. For the prediction problem with \(\mathcal{A}\) as in (14), they are given in Section 8 and Algorithm 2.

## 5 Alternative Loss Functions

As mentioned above, the inexact application of the operator \(\mathcal{A}\) can be accomplished by minimizing loss functions other than (17). The key requirement for a loss function in the present study is that \(\mathcal{A}\) appears in its gradient only through terms of the form \(\left\langle f,\mathcal{A}g\right\rangle_{\mu}\) for some functions \(f\) and \(g\), so that the gradient can be estimated using trajectory data. As a result, we have flexibility in the choice of loss and, in turn, the representation of \(u\). In particular, we consider the representation \(u_{\theta}=\zeta(z_{\theta})\), where \(\zeta\) is an increasing function, and \(z_{\theta}\) is a function parameterized by a neural network. An advantage of doing so is that the function \(\zeta\) can restrict the output values of \(u_{\theta}\) to some range. For example, when computing a probability such as the committor, a natural choice is the sigmoid function \(\zeta(x)=(1+e^{-x})^{-1}\).

Our goal is to train a sequence of parameter values so that \(u_{\theta^{s}}\) approximately follows a subspace iteration, i.e., so that \(\zeta(z_{\theta^{s+1}})\approx\mathcal{A}u_{\theta^{s}}\). To this end, we minimize with respect to \(\theta\) the loss function

\[\mathbb{E}_{X^{0}\sim\mu}\left[V(z_{\theta})-z_{\theta}\mathcal{A}u_{\theta^{ s}}\right], \tag{22}\]

where \(V\) is an antiderivative of \(\zeta\), and \(\theta^{s}\) is fixed. The subscript \(X^{0}\sim\mu\) in this expression indicates that \(X^{0}\) is drawn from \(\mu\). Note that, as desired, \(\mathcal{A}\) appears in the gradient of (22) only in an inner product of the required form, and an approximate minimizer, \(\theta^{s+1}\), of this loss satisfies \(\zeta(z_{\theta^{s+1}})\approx\mathcal{A}u_{\theta^{s}}\). This general form of loss function is adapted from variational expressions for the divergence of two probability measures [32, 40].

The \(L_{\mu}^{2}\) loss in (17), which we use in several of our numerical experiments, corresponds to the choice \(\zeta(x)=x\) and \(V(x)=x^{2}/2\). The choice of \(\zeta(x)=(1+e^{-x})^{-1}\) mentioned above corresponds to \(V(x)=\log(1+e^{x})\); we refer to the loss in (22) with this choice of \(V\) as the "softplus" loss. We note that in the context of committor function estimation, in the limit of infinite lag time the "softplus" loss corresponds directly to the maximum log-likelihood loss (for independent Bernoulli random variables) that defines the logistic regression estimate of the probability of reaching \(B\) before \(A\). Logistic regression has previously been used in conjunction with data sets of trajectories integrated all the way until reaching \(A\) or \(B\) to estimate the committor function[41, 42, 43, 44, 45, 46].

## 6 Test Problems

We illustrate our methods with two well-characterized systems that exemplify features of molecular transitions. In this section, we provide key details of these systems.

Figure 1: Müller-Brown potential energy surface. The orange and red ovals indicate the locations of states \(A\) and \(B\) respectively when computing predictions. Contour lines are drawn every 1 \(\beta^{-1}\).
### Muller-Brown potential

The first system is defined by the Muller-Brown potential [47] (Figure 1):

\[V_{\rm MB}(y,z) =\frac{1}{20}\sum_{i=1}^{4}C_{i}\exp[a_{i}(y-y_{i})^{2}\] \[+b_{i}(y-y_{i})(z-z_{i})+c_{i}(z-z_{i})^{2}]. \tag{23}\]

The two-dimensional nature of this model facilitates visualization. The presence of multiple minima and saddlepoints that are connected by a path that does not align with the coordinate axes makes the system challenging for both sampling and analysis methods. In Sections 7.7.1 and 8.8.1, we use \(C_{i}=\{-200,-100,-170,15\}\), \(a_{i}=\{-1,-1,-6.5,0.7\}\), \(b_{i}=\{0,0,11,0.6\}\), \(c_{i}=\{-10,-10,-6.5,0.7\}\), \(y_{i}=\{1,-0.27,-0.5,-1\}\), \(z_{i}=\{0,0.5,1.5,1\}\). In Section 8.8.2, we tune the parameters to make transitions between minima rarer; there, the parameters are \(C_{i}=\{-250,-150,-170,15\}\), \(a_{i}=\{-1,-3,-6.5,0.7\}\), \(b_{i}=\{0,0,11,0.6\}\), \(c_{i}=\{-10,-30,-6.5,0.7\}\), \(y_{i}=\{1,-0.29,-0.5,-1\}\), \(z_{i}=\{0,0.5,1.5,1\}\). For prediction, we analyze transitions between the upper left minimum (\(A\)) and the lower right minimum (\(B\)) in Figure 1; these states are defined as

\[A =\{y,z:6.5(y+0.5)^{2}-11(y+0.5)(z-1.5)+6.5(z-1.5)^{2}<0.3\} \tag{24}\] \[B =\{y,z:(y-0.6)^{2}+0.5(z-0.02)^{2}<0.2\}.\]

To generate a data set, we randomly draw 50,000 initial conditions \(X_{j}^{0}\) uniformly from the region

\[\Omega=\{y,z:-1.5<y<1.0,\ -0.5<z<2.5,\ V_{\rm MB}(y,z)<12\} \tag{25}\]

and, from each of these initial conditions, generate one trajectory according to the discretized over-damped Langevin dynamics

\[X_{j}^{t}=X_{j}^{t-1}-\delta_{t}\,\nabla V_{\rm MB}(X_{j}^{t-1})+\sqrt{\delta _{t}\,2\beta^{-1}}\,\xi_{j}^{t} \tag{26}\]

where \(1\leq t\leq\tau\), the \(\xi_{j}^{t}\) are independent standard Gaussian random variables, and the timestep is \(\delta_{t}=0.001\). We use an inverse temperature of \(\beta=2\), and we vary \(\tau\) as indicated below (note that \(\tau\) is an integer number of steps, but we present our results for the Muller-Brown model in terms of the dimensionless scale used for \(\delta_{t}\) immediately above). For each test, we use independent random numbers and redraw the initial conditions because it is faster to generate the trajectories than to read them from disk in this case. All error bars are computed from the empirical standard deviation over 10 replicate data sets.

To validate our results, we compare the neural-network results against grid-based references, computed as described in Section S4 of the Supplementary Material of Ref. [11] and the Appendix of Ref. [48] (our notation here follows the former more closely). The essential idea is that the terms in the infinitesimal generator of the transition operator can be approximated on a grid to second order in the spacing by expanding both the probability for transitions between sites and a test function. To this end, we define the transition matrix for neighboring sites \(x=(y,z)\) and \(x^{\prime}=(y\pm\delta_{x},z)\) or \((y,z\pm\delta_{x})\) on a grid with spacings \(\delta_{x}\):

\[P(x^{\prime}|x)=\frac{1}{1+e^{\beta[V(x^{\prime})-V(x)]}} \tag{27}\]

(this definition differs from that in Ref. [11] by a factor of \(1/4\), and we scale \(P-I\), where \(I\) is the identity matrix, accordingly to set the time units below). The diagonal entry \(P(x|x)\) is fixed to make the transition matrix row-stochastic. We use \(\delta_{x}=0.0125\); the grid is a rectangle with \(y\in[-1.5,1.0]\), and \(z\in[-0.5,2.0]\). The reference invariant subspaces are computed by diagonalizing the matrix \(2(P-I)/\beta\delta_{x}^{2}\) with a sparse eigensolver; we use scipy.sparse.linalg. We obtain the reference committor \(\hat{q}_{+}\) for grid sites in \((A\cup B)^{\epsilon}\) by solving

\[{\rm diag}(\hat{1}_{(A\cup B)^{\epsilon}}) (P-I)\hat{q}_{+}\] \[=-{\rm diag}(\hat{1}_{(A\cup B)^{\epsilon}})(P-I)\hat{1}_{B} \tag{28}\]

and for grid sites in \(A\cup B\) by setting \(\hat{q}_{+}=\hat{1}_{B}\). Here, \(\hat{q}_{+}\) and \(\hat{1}_{B}\) are vectors corresponding to the functions evaluated at the grid sites.
### \(\text{AIB}_{9}\) helix-to-helix transition

The second system is a peptide of nine \(\alpha\)-aminoisobutyric acids (\(\text{AIB}_{9}\); Figure 2). Because AIB is achiral around its \(\alpha\)-carbon atom, \(\text{AIB}_{9}\) can form left- and right-handed helices with equal probabilities, and we study the transitions between these two states. This transition was previously studied using MSMs and long unbiased molecular dynamics simulations [49; 50]. \(\text{AIB}_{9}\) poses a stringent test due to the presence of many metastable intermediate states.

The states are defined in terms of the internal \(\phi\) and \(\psi\) dihedral angles. We classify an amino acid as being in the "left" state if its dihedral angle values are within a circle of radius \(25^{\circ}\) centered at \((41^{\circ},43^{\circ})\), that is

\[(\phi-41^{\circ})^{2}+(\psi-43^{\circ})^{2}\leq(25^{\circ})^{2}.\]

Amino acids are classified as being in the "right" state using the same radius, but centered instead at \((-41^{\circ},-43^{\circ})\). States \(A\) and \(B\) are defined by the amino acids at sequence positions 3-7 being all left or all right, respectively. We do not use the two residues on each end of \(\text{AIB}_{9}\) in defining the states as these are typically more flexible [50]. The states can be resolved by projecting onto dihedral angle principal components (dPCs; Figure 2, right) as described previously [51].

Following a procedure similar to that described in Ref. [50], we generate a data set of short trajectories. From each of the 691 starting configurations in Ref. [50], we simulate 10 trajectories of duration 20 ns with initial velocities drawn randomly from a Maxwell-Boltzmann distribution at a temperature of 300 K. The short trajectory data set thus contains 6,910 trajectories, corresponding to a total sampling time of 138.2 \(\mu\)s. We use a timestep of 4 fs together with a hydrogen mass repartitioning scheme [52], and configurations are saved every 40 ps. We employ the AIB parameters from Forcefield_NCAA [53] and the GBNeck2 implicit-solvent model [54]. Simulations are performed with the Langevin integrator in OpenMM 7.7.0 [55] using a friction coefficient of \(1\,\text{ps}^{-1}\). To generate a reference for comparison to our results, we randomly select 20 configurations from the data set above and, from each, run a simulation of 5 \(\mu\)s with the same simulation parameters. For all following tests on \(\text{AIB}_{9}\), batches consist of pairs of frames separated by \(\tau\) drawn randomly with replacement from the short trajectories (i.e., from all possible such pairs in the data set). From each frame, we use only the atom positions because the momenta decorrelate within a few picoseconds, which is much shorter than the lag times that we consider. However, in principle, the momenta could impact the prediction functions [56] and be used as neural-network inputs as well.

## 7 Spectral estimation

In this section, we provide some further numerical details for the application of our method to spectral estimation and demonstrate the method on the test problems. For our subspace iteration, we require estimators for inner products of the form \(\langle f,\mathcal{T}^{\tau}g\rangle_{\mu}\). For example, the gradient of the loss function (17) involves inner products of the form

\[\left\langle\nabla_{\theta}\varphi_{\theta}^{a},\mathcal{T}^{\tau}\varphi_{ \theta}^{b}\right\rangle_{\mu}. \tag{29}\]

For these, we use the unbiased data-driven estimator

\[\left\langle f,\mathcal{T}^{\tau}g\right\rangle_{\mu}\approx\frac{1}{n}\sum_{ j=1}^{n}f(X_{j}^{0})g(X_{j}^{\tau}). \tag{30}\]

As discussed in Section 4, applying the operator \(\mathcal{T}^{\tau}\) repeatedly causes each basis function to converge to the dominant eigenfunction and leads to

Figure 2: Helix-to-helix transition of \(\text{AIB}_{9}\). (left) Left- and right-handed helices, which we use as state \(A\) and \(B\), respectively, when computing predictions. Carbon, nitrogen, and oxygen atoms are shown in yellow, blue, and red, respectively; hydrogen atoms are omitted. (right) Potential of mean force constructed from the histogram of value pairs of the first two dihedral angle principal components; data are from the 20 trajectories of 5 \(\mu\)s that we use to construct reference statistics (see text). The left-handed helix corresponds to the left-most basin, and the right-handed helix corresponds to the right-most basin. Contour lines are drawn every 2 \(\beta^{-1}\) corresponding to a temperature of 300 K.
numerical instabilities. To avoid this, we orthogonalize the outputs of the networks with a QR decomposition at the end of each subspace iteration by constructing the matrix \(\Phi_{ia}=\varphi_{\theta}^{a}(X_{s}^{s})\) and then computing the factorization \(\Phi=QR\) where \(Q\) is an \(n\times k\) matrix with orthogonal columns and \(R\) is an upper triangular \(k\times k\) matrix. Finally, we set \(\tilde{\varphi}_{s}^{a}=\sum_{b=1}^{k}\varphi_{\theta}^{b}\left(R^{-1}N\right)_{ ba}\), where \(N\) is a diagonal matrix with entries equal to the norms of the columns of \(\Phi\) (before orthogonalization). To ensure that the networks remain well-separated (i.e., the eigenvalues of \(C^{0}\) remain away from zero), we penalize large off-diagonal entries of \(K\) by adding to the loss

\[\gamma_{1}\|K-\text{diag}(K)\|_{\text{F}}^{2}, \tag{31}\]

where \(\gamma_{1}\) allows us to tune the strength of this term relative to others, and \(\|\cdot\|_{\text{F}}\) is the Frobenius norm. We control the scale of each network output using the strategy from Ref. [32]; that is, we add to the loss a term of the form

\[\gamma_{2}\sum_{a=1}^{k}(2\nu_{a}(\langle\varphi_{\theta}^{a},\varphi_{\theta} ^{a}\rangle_{\mu}-1)-\nu_{a}^{2}), \tag{32}\]

where we have introduced the conjugate variables \(\nu_{a}\) which we maximize with gradient ascent (or similar optimization). In general, our numerical experiments suggest that it is best to keep \(\gamma_{1}\) and \(\gamma_{2}\) relatively small. We find that stability of the algorithm over many subspace iterations is improved if the matrix \(K\) is set at its optimal value before each inner loop. To do this, we set

\[K_{1:i,i}=\arg\min_{c}\left\|\sum_{a=1}^{i}\varphi_{\theta}^{a}c_{a}-\mathcal{ T}^{\tau}\tilde{\varphi}_{s}^{a}\right\|+\gamma_{2}\sum_{a=1}^{i-1}c_{a}^{2}. \tag{33}\]

The above minimization can be solved with linear least squares. Finally, we note that in practice any optimizer can be used for the inner iteration steps, though the algorithm below implements stochastic gradient descent. In this work, we use Adam[57] for all numerical tests. We summarize our procedure for spectral estimation in Algorithm 1.

### Muller-Brown model

As a first test of our method, we compute the \(k=3\) dominant eigenpairs for the Muller-Brown model. Since we know that the dominant eigenfunction of the transition operator is the constant function \(v_{1}(y,z)=1\) with eigenvalue \(\lambda_{1}=1\), we directly include this function in the basis as a non-trainable function, i.e. \(\varphi_{\theta}^{1}(y,z)=1.\) To initialize \(\tilde{\varphi}_{1}^{a}\) for each \(a>1\), we choose a standard Gaussian vector \((Y^{a},Z^{a})\in\mathbb{R}^{2}\), and set \(\tilde{\varphi}_{1}^{a}(y,z)=y\,Y^{a}+z\,Z^{a}\). This ensures that the initial basis vectors are well-separated and the first QR step is numerically stable. Here and in all subsequent Muller-Brown tests, batches of trajectories are drawn from the entire data set with replacement. Other hyperparameters are listed in Table 1.

Figure 3 shows that we obtain good agreement between the estimate produced by the inexact subspace iteration in Algorithm 1 and reference eigenfunctions. Figure 4 (upper panels) shows how the corresponding eigenvalues vary with lag time; again there is good agreement with the reference. Furthermore, there is a significant gap between \(\lambda_{3}\) and \(\lambda_{4}\), indicating that a three-dimensional subspace captures the dynamics of interest for this system.

We compare the subspace that we obtain from our method with that from an MSM constructed from the same amount of data by using \(k\)-means to cluster the configurations into 400 states and counting the transitions between clusters. This is a very fine discretization for this system, and the MSM is sufficiently expressive to yield eigenfunctions in good agreement with the reference. The relative error of \(1-\lambda_{2}\) is comparable for the two methods (Figure 4, lower left). To compare two finite dimensional subspaces, \(\mathcal{U}\) and \(\mathcal{V}\), we define the subspace distance as[4]

\[d(\mathcal{U},\mathcal{V})=\left\|(I-P_{\mathcal{U}})P_{\mathcal{V}}\right\|_ {\text{F}}, \tag{34}\]

where \(P_{\mathcal{U}}\) and \(P_{\mathcal{V}}\) denote the orthogonal projections onto \(\mathcal{U}\) and \(\mathcal{V}\), respectively, and \(\|\cdot\|_{\text{F}}\) is the Frobenius norm. Figure 4 (lower right) shows the subspace distances from the reference as functions of lag time. We see that the inexact subspace iteration better approximates the three-dimensional dominant eigenspace for moderate to long lag times, even though the eigenvalues are comparable.

### Aib\({}_{\text{9}}\)

For the molecular test system, we compute the dominant five-dimensional subspace. We compare the inexact subspace iteration in Algorithm 1 with MSMs constructed on dihedral angles ("dihedral MSM") and on Cartesian coordinates ("Cartesian MSM"). We expect the dihedral MSM to be accurate given that the dynamics of AIB\({}_{\text{9}}\) are well-described by the backbone dihedral angles [49, 50], and we thus
use it as a reference. It is constructed by clustering the sine and cosine of each of the backbone dihedral angles (\(\phi\) and \(\psi\)) for the nine residues (for a total of \(2\times 2\times 9=36\) input features) into 1000 clusters using \(k\)-means and counting transitions between clusters. The Cartesian MSM is constructed by similarly counting transitions between 1000 clusters from the \(k\)-means algorithm, but the clusters are based on the Cartesian coordinates of all non-hydrogen atoms after aligning the backbone atoms of the trajectories, for a total of 174 input features. Because of the difficulty of clustering high-dimensional data, we ex
pect the Cartesian MSM basis to give poor results. The neural network for the inexact subspace iteration receives the same 174 Cartesian coordinates as input features. We choose to use Cartesian coordinates rather than dihedral angles as inputs because it requires the network to identify nontrivial representations for describing the dynamics.

As in the Muller-Brown example, we use \(\varphi_{\theta}^{1}=1\) and a random linear combination of coordinate functions to initialize \(\bar{\varphi}_{1}^{a}\) for \(a>1\). Other hyperparameters are listed in Table 1. With these choices, the neural-network training typically requires about 20 minutes on a single NVIDIA A40 GPU; this is much longer than the time required for diagonalization of the \(1000\times 1000\) MSM transition matrix, which is nearly instantaneous. However, the time for neural-network training is negligible compared with the time to generate the data set, which is the same for both approaches.

Taking the dihedral MSM as a reference, the Cartesian MSM systematically underestimates the eigenvalues (Figure 5). The subspace iteration is very accurate for the first four eigenvalues but the estimates for the fifth are low and vary considerably from run to run. A very small gap between \(\lambda_{4}\) and \(\lambda_{5}\) may contribute to the difficulty in estimating \(\lambda_{5}\). In Figure 6, we plot the first two non-trivial eigenfunctions (\(v_{2}\) and \(v_{3}\)), which align with the axes of the dPC projection. The eigenfunction \(v_{2}\) corresponds to the transition between the left- and right-handed helices; the eigenfunction \(v_{3}\) is nearly orthogonal to \(v_{2}\) and corresponds to transitions between intermediate states. It is challenging to visualize the remaining two eigenfunctions by projecting onto the first two dPCs because \(v_{4}\) and \(v_{5}\) are orthogonal to \(v_{2}\) and \(v_{3}\). The estimates for \(v_{2}\) are in qualitative agreement for all lag times tested (Figure 6 shows results for \(\tau\) corresponding to 40 ps), but the subspace iteration results are less noisy for the shortest lag times. Moreover, the estimate for \(v_{3}\) from subspace iteration agrees more closely with that from the dihedral MSM than does the estimate for \(v_{3}\) from the Cartesian MSM. The subspace distance for \(v_{2}\) and \(v_{3}\) between the subspace iteration and the dihedral MSM is 0.947, compared with a value of 0.969 for the subspace distance between the two MSMs. Together, our results indicate that the neural networks are able to learn the leading eigenfunctions and eigenvalues of the transition operator (dynamical modes) of this system despite being presented with coordinates that are not the natural ones for describing the dynamics.

## VIII Prediction

Inexact subspace iteration for \(\mathcal{A}\) in (14) is equivalent to performing the inexact Richardson iteration in (12) on the first basis function \(\varphi_{\theta}^{1}\) and then performing an inexact subspace iteration for the operator \(\mathcal{S}^{\tau}\) on the rest of the basis functions. The iteration requires unbiased estimators of the forms

\[\langle f,\mathcal{S}^{\tau}g\rangle_{\mu}\approx\frac{1}{n}\sum_{j=1}^{n}f(X_ {j}^{0})g(X_{j}^{\tau\wedge T_{j}}) \tag{35}\]

and

\[\langle f,r\rangle_{\mu}\approx\frac{1}{n}\sum_{j=1}^{n}f(X_{j}^{0})\sum_{t=0 }^{(\tau\wedge T_{j})-1}\Gamma(X_{j}^{t}), \tag{36}\]

where \(T_{j}\) is the first time \(X_{j}^{t}\) enters \(D^{\mathsf{c}}\) and \(r\) is the right-hand side of (5), as previously.

The Richardson iterate, \(\varphi_{\theta}^{1}\), must satisfy the boundary condition \(\varphi_{\theta}^{1}(x)=\Psi(x)\) for \(x\notin D\). The other basis functions should satisfy \(\varphi_{\theta}^{a}(x)=0\) for \(x\notin D\). In practice, we enforce the boundary conditions by explicitly setting \(\varphi_{\theta}^{1}(x)=\Psi(x)\) and \(\varphi_{\theta}^{a}(x)=0\) for \(a>1\) when \(x\notin D\).

When the boundary condition is zero, as for the MFPT, we find an approximate solution of the form

\[u_{\theta}=\sum_{a=1}^{k}w_{a}\varphi_{\theta}^{a} \tag{37}\]

Figure 3: First two non-trivial eigenfunctions of the Müller-Brown model. (top) Grid-based reference. (bottom) Neural network subspace after ten subspace iteration steps, computed with \(\tau=300\) steps (i.e., \(0.3\,\delta_{\mathrm{t}}^{-1}\)).
by solving the \(k\)-dimensional linear system

\[\left(C^{0}-C^{\tau}\right)w=p \tag{38}\]

where, for \(a,b\geq 1\),

\[C_{ab}^{t}=\langle\varphi_{\theta}^{a},\mathcal{S}^{t}\varphi_{\theta}^{b} \rangle_{\mu} \tag{39}\]

for \(t=0,\tau\), and

\[p_{a}=\left\langle\varphi_{\theta}^{a},\mathbb{E}_{x}\left[\rho(X)\right] \right\rangle_{\mu}. \tag{40}\]

Figure 4: Spectral estimation as a function of lag time (in units of \(\delta_{t}^{-1}\)) for the Müller-Brown model. (top left) Second eigenvalue. (top right) Third and fourth eigenvalues; only the reference fourth eigenvalue is shown to illustrate the spectral gap. (bottom left) Relative error in the first spectral gap (i.e., \(1-\lambda_{2}\)). (bottom right) Subspace distance between estimated and reference three-dimensional invariant subspaces.

Figure 5: First five eigenvalues of the transition operator for \(\text{AIB}_{9}\) as a function of lag time. (left) Comparison between eigenvalues computed using the dihedral MSM with 1000 clusters (solid lines) and the inexact subspace iteration (dashed lines). The shading indicates standard deviations over five trained networks for the subspace iteration. (right) Comparison between a dihedral MSM (solid lines) and Cartesian MSMs with 1000 clusters (dashed lines). The standard deviations for the Cartesian MSMs over five random seeds for \(k\)-means clustering are too narrow to be seen.
In (40), we introduce the notation

\[\rho(X)=\sum_{t=0}^{(\tau\wedge T)-1}\Gamma(X^{t}) \tag{41}\]

for use in Algorithm 2.

When the boundary condition is non-zero, as for the committor, we restrict (38) to a \((k-1)\)-dimensional linear system by excluding the indices \(a=1\) and \(b=1\) in (39) and (40) and setting

\[\rho(X)=\varphi_{\theta}^{1}(X^{\tau\wedge T})-\varphi_{\theta}^{1}(X^{0})+ \sum_{t=0}^{(\tau\wedge T)-1}\Gamma(X^{t}). \tag{42}\]

In this case the corresponding approximate solution is

\[u_{\theta}=\varphi_{\theta}^{1}+\sum_{a=2}^{k}w_{a}\varphi_{\theta}^{a}. \tag{43}\]

This approximate solution corresponds to the one given by dynamical Galerkin approximation [10; 11] with the basis \(\{\varphi_{\theta}^{a}\}_{a=2}^{k}\) and a "guess" function of \(\varphi_{\theta}^{1}\).

When the boundary conditions are zero, the orthogonalization procedure and the matrix \(K\) are applied to all basis functions as in Section 7. When the boundary conditions are non-zero, the orthogonalization procedure is only applied to the basis functions \(\{\varphi_{\theta}^{a}\}_{a=2}^{k}\), and \(K_{a1}=I_{a1}\). We summarize our procedure for prediction in Algorithm 2.

### Muller-Brown committor

In this section, we demonstrate the use of our method for prediction by estimating the committor for the Muller-Brown model with a shallow intermediate basin at \((-0.25,0.5)\) (Figure 1). Here the sets \(A\) and \(B\) are defined as in Eq. (24) and \(T\) is the time of first entrance to \(D^{\mathsf{c}}=A\cup B\). In this case, a one-dimensional subspace iteration (i.e., \(k=1\) in Algorithm 2) appears sufficient to accurately solve the prediction problem. Figure 7 shows the largest eigenvalue of the stopped transition operator \(\mathcal{S}^{\tau}\) (the second largest of \(\mathcal{A}\) in (14)) computed from our grid-based reference scheme (Section 6.6.1). Richardson iteration should converge geometrically in this eigenvalue[33], and so, for moderate lag times,

Figure 6: First two non-trivial eigenfunctions of \(\mathrm{AIB}_{9}\) projected onto the first two dPCs (i.e., averaged for bins in the two-dimensional space shown). (top) MSM constructed on sine and cosine of dihedral angles with 1000 clusters and lag time corresponding to 40 ps. (middle) Inexact subspace iteration using Cartesian coordinates and the same lag time. (bottom) MSM constructed on Cartesian coordinates with 1000 clusters and the same lag time.

Figure 7: First eigenvalue of \(\mathcal{S}^{\tau}\) (second of \(\mathcal{A}\) in (14)) for the Muller-Brown model as a function of lag time (in units of \(\delta_{t}^{-1}\)). The gap between this eigenvalue and the dominant eigenvalue, which is one, determines the rate of convergence of the Richardson iteration.
we can expect our method to converge in a few dozen iterations. To initialize the algorithm we choose \(\tilde{\varphi}_{1}^{1}=\mathbb{1}_{B}\). All other hyperparameters are listed in Table I.

We compare the estimate of the committor from our approach with that from an MSM constructed from the same amount of data by using \(k\)-means to cluster the configurations outside \(A\) and \(B\) into \(400\) states and counting the transitions between clusters. In addition to the root mean square error (RMSE) for the committor itself, we show the RMSE of

\[\text{logit}_{\varepsilon}(q)=\log\left(\frac{\varepsilon+q}{1+\varepsilon-q }\right) \tag{44}\]

for points outside \(A\) and \(B\). This function amplifies the importance of values close to zero and one. We include \(\varepsilon\) because we want to assign only a finite penalty if the procedure estimates \(q\) to be exactly zero or one; we use \(\varepsilon=e^{-20}\).

Results as a function of lag time are shown in Figure 8. We see that the Richardson iterate is more accurate than the MSM for all but the shortest lag times. When using the \(L_{\mu}^{2}\) loss, the results are comparable, whereas the softplus loss allows the Richardson iterate to improve the RMSE of the logit function in (44) with no decrease in performance with respect to the RMSE of the committor. Results as a function of the size of the data set are shown in Figure 9 for a fixed lag time of \(\tau=0.1\delta_{t}^{-1}\). The Richardson iterate generally does as well or better than the MSM. Again, the differences are more apparent in the RMSE of the logit function in (44). By that measure, the Richardson iterate obtained with both loss functions is significantly more accurate than the MSM for small numbers of trajectories. The softplus loss maintains an advantage even
for large numbers of trajectories.

Figure 8: Committor prediction for the Müller-Brown system as a function of lag time (in units of \(\delta_{t}^{-1}\)). (left) Comparison of the inexact Richardson iteration using the \(L_{\mu}^{2}\) loss and an MSM with 400 states. (right) Same comparison using the softplus loss in place of the \(L_{\mu}^{2}\) loss.

Figure 9: Committor prediction for the Müller-Brown as a function of number of initial conditions for a fixed lag time of \(\tau=0.1\delta_{t}^{-1}\). (left) Comparison of inexact Richardson iteration using the \(L_{\mu}^{2}\) loss and an MSM with 400 states. (right) Same comparison using the softplus loss in place of the \(L_{\mu}^{2}\) loss.
### Accelerating convergence by incorporating eigenfunctions

As discussed in Section 3, we expect Richardson iteration to converge slowly when the largest eigenvalue of \(\mathcal{S}^{\tau}\), \(\lambda_{1}\), is close to \(1\). More precisely, the number of iterations required to reach convergence should scale with \(-1/\log\lambda_{1}=\mathbb{E}\left[T\right]/\tau\), the mean escape time from the quasi-stationary distribution to the boundary of \(D\) divided by the lag time. With this in mind, we can expect inexact Richardson iteration for the Muller-Brown to perform poorly if we deepen the intermediate basin at \((-0.25,0.5)\) as in Figure 10 (top left). Again, the sets \(A\) and \(B\) are defined as in (24), and \(T\) is the time of first entrance to \(D^{\mathsf{c}}=A\cup B\). In this case, \(-1/\log\lambda_{1}\) is on the order of \(100\) for the lag times we consider and, as expected, inexact Richardson iteration converges slowly (Figure 10, bottom left). Estimates of the committor by inexact Richardson iteration do not reach the correct values even after hundreds of iterations (Figure 10, bottom right).

We now show that convergence can be accelerated dramatically by incorporating additional eigenfunctions of \(\mathcal{S}^{\tau}\) (i.e., \(k>1\) in Algorithm 2). For the Muller-Brown model with a deepened intermediate basin, the second eigenvalue of \(\mathcal{S}^{\tau}\) is roughly \(1\times 10^{-4}\) for a lag time of \(\tau=1000\) steps or \(1\,\delta_{t}^{\,-1}\) (while the first is near one as discussed above). We therefore choose \(k=2\) with \(\tilde{\varphi}_{1}^{2}\) initialized as a random linear combination of coordinate functions as in previous examples. We run the subspace iteration for four iterations, compute the committor as a linear combination of the resulting functions, and then refine this result with a further ten Richardson iterations (i.e., \(k=1\) with the starting vector as the output of the \(k=2\) subspace iteration). To combine the functions, we use a linear solve which incorporates memory (Algorithm 3) [58, 59]. We find that the use of memory improves the data-efficiency substantially for poorly conditioned problems. For our tests here, we use two memory kernels, corresponding to \(\tau_{\text{mem}}=\lfloor\tau/4\rfloor\).

The bottom row of Figure 11 illustrates the idea of the subspace iteration. The second eigenfunction (Figure 11, center) is peaked at the intermediate. As a result, the two neural-network functions linearly combined by the Galerkin approach with memory can yield a good result for the committor (Figure 11, bottom right). Figure 12 compares the RMSE for the committor and the RMSE for the logit in (44) for Algorithm 2 with \(k=1\) (pure Richardson iteration) and \(k=2\) (incorporating the first non-trivial eigenfunction), and an MSM with \(400\) states. We see that the Richardson iteration suffers large errors at all lag times; as noted previously, this error is mainly in the vicinity of the intermediate. The MSM cannot accurately compute the small probabilities, but does as well as the subspace iteration in terms of RMSE.

### \(\text{AIB}_{9}\) prediction results

As an example of prediction in a high-dimensional system, we compute the committor for the transition between the left- and right-handed helices of \(\text{AIB}_{9}\) using the inexact Richardson iteration scheme (\(k=1\) in Algorithm 2) with the softplus loss function. Specifically, for this committor calculation \(T\) is the time of first entrance to \(D^{\mathsf{c}}=A\cup B\) with \(A\) and \(B\) defined in Section 6.2. As before, we initialize \(\tilde{\varphi}_{1}^{2}=\mathbbm{1}_{B}\).

To validate our results, we use the \(5\)\(\mu\)s reference trajectories to compute an empirical committor as a function of the neural network outputs, binned into intervals:

\[\bar{q}(s)=\mathbb{P}\left[X^{T}\in B\Big{|}u_{\theta}(X^{0})\in[s,s+\Delta s ]\right] \tag{45}\]

for \(s\in[0,1-\Delta s]\). Here, we use \(\Delta s=0.05\). The overall error in the committor estimate is defined as

\[q\text{ error}=\left(\Delta s\sum_{n=0}^{1/\Delta s-1}\left[\bar{q}(n\Delta s )-n\Delta s\right]^{2}\right)^{1/2}. \tag{46}\]

While this measure of error can only be used when the data set contains trajectories of long enough duration to reach \(D^{\mathsf{c}}\), it has the advantage that it does not depend on the choice of projection that we use to visualize the results.

Results for the full data set with \(\tau\) corresponding to \(400\) ps are shown in Figure 13. The projection on the principal components is consistent with the symmetry of the system, and the predictions show good agreement with the empirical committors. As \(\tau\) decreases, the results become less accurate (Figure 14, top left); at shorter lag times we would expect further increases in the error.

We also examine the dependence of the results on the size of the data set by subsampling the short trajectories and then training neural networks on the reduced set of trajectories (Figure 14, top right). We find that the performance steadily drops as the number of trajectories is reduced and degrades rapidly (Figure 14, bottom) for the data sets subsampled more than \(20\)-fold, corresponding to about \(7\)\(\mu\)s of total sampling.
Finally, we compute the MFPT to reach the right-handed helix using the same data set. For the MFPT calculation \(T\) is the time of first entrance to \(D^{\mathfrak{c}}=B\). Note that the time of first entrance to \(B\) includes

Figure 10: Richardson iteration for the committor converges slowly for a Müller-Brown potential with a deepened intermediate. (top left) Potential energy surface, with states \(A\) and \(B\) indicated. Contour lines are drawn every 1 \(\beta^{-1}\). (top right) Reference committor. (bottom left) Dominant eigenvalue as a function of lag time (in units of \(\delta_{t}^{-1}\)) from an MSM with 400 states, subspace iteration, and the grid-based reference. (bottom right) Example of the Richardson iteration after 400 iterations. Note the overfitting artifacts and lack of convergence near the intermediate state.
long dwell times in \(A\) and is expected to be much larger than the time of first entrance to \(A\cup B\).

We compare against an empirical estimate of the MFPT defined by

\[\bar{m}(s)=\operatorname{\mathds{E}}\Big{[}T\Big{|}u_{\theta}(X^{0})\in[s,s+ \Delta s]\Big{]} \tag{47}\]

for \(s\in[0,m_{\max}-\Delta s]\) where \(\Delta s=3\) and \(m_{\max}=57\) ns. Overall error is defined analogously to Eq. (46).

In Figure 15, we show the MFPT obtained from Algorithm 2 with \(k=1\) (pure Richardson iteration) and \(k=5\). Both calculations use the \(L_{\mu}^{2}\) loss function. Initially we set \(\tilde{\varphi}_{1}^{1}\) equal to an arbitrary positive function (we use \(51_{A}\)) and \(\tilde{\varphi}_{s}^{a}\) for \(a>1\) to a random linear combination of coordinate functions. The horizontal line in Figure 15 indicates a MFPT of about 56 ns estimated from the long reference trajectories. We see that the algorithm with \(k=5\) converges much faster (note the differing scales of the horizontal axes) and yields accurate results at all lag times other than the shortest shown. The need to choose \(k>1\) for this MFPT calculation is again consistent with theoretical convergence behavior of exact subspace iteration. Because the typical time of first entrance to \(B\) from points in \(A\) is very large, we expect the dominant eigenvalue of \(\mathcal{S}^{\tau}\) to be very near to one when \(D=B^{\epsilon}\). In contrast, the committor calculation benefits from the fact that the time of first entrance to \(A\cup B\) is much shorter, im

Figure 11: Illustration of the subspace iteration for the Müller-Brown committor. (top left) Modified Müller-Brown potential. (top center) Reference second eigenfunction. (top right) Reference committor. (bottom left) Neural-network Richardson iterate after four iterations. (bottom center) First non-dominant eigenfunction obtained from the neural network after four iterations. (bottom right) Committor resulting from linear combination of the Richardson iterate and second eigenfunction. Results shown are for \(\tau=1000\) steps (i.e., \(1_{\delta_{t}}^{-1}\)).

Figure 12: Committor for the Müller-Brown potential with deepened intermediate as a function of lag time (in units of \(\delta_{t}^{-1}\)). (left) Comparison of RMSE for subspace iteration as described above, Richardson iteration (as in Section 8.1 but instead with 500 subspace iterations), and an MSM with 400 states. (right) RMSE of the logit function in (44).
plying a smaller dominant eigenvalue of \(\mathcal{S}^{\tau}\) when \(D=\left(A\cup B\right)^{\mathsf{c}}\).

## IX Conclusions

In this work we have presented a method for spectral estimation and rare-event prediction from short-trajectory data. The key idea is that we use the data as the basis for an inexact subspace iteration. For the test systems that we considered, the method not only outperforms high-resolution MSMs, but it can be tuned through the choice of loss function to compute committor probabilities accurately near the reactants, transition states, and products. Other than the Markov assumption, our method requires no knowledge of the underlying model and puts no restrictions on its dynamics.

As discussed in prior neural-network based prediction work[24, 30], our method is sensitive to the quality and distribution of the initial sampling data. However, our work shares with Ref. [24] the major advantage of allowing the use of arbitrary inner products. This enables adaptive sampling of the state space [24, 60] and--together with the features described above--the application to observational and experimental data, for which the stationary distribution is generally unavailable.

In the present work, we focused on statistics of transition operators, but our method can readily be extended to solve problems involving their adjoint operators as well. By this means, we can obtain the stationary distribution as well as the backward committor. The combination of forward and backward

Figure 14: AIB\({}_{9}\) committor for the transition between left- and right-handed helices, as functions of lag time (in ps) and number of initial conditions. (top left) Error in the committor as a function of lag time (in ps). Shading indicates the standard deviation over ten different initializations of the neural-network parameters. (top right) Error in the committor as a function of the number of initial conditions with \(\tau\) corresponding to 160 ps. Shading indicates the standard deviation over ten different random samples of the trajectories. (bottom) Comparison between empirical committors and neural-network committors trained on data sets with (left) 1/2 and (right) 1/20 of the short trajectories. Error bars indicate standard deviations over ten random samples of the trajectories.

Figure 13: AIB\({}_{9}\) committor for the transition between left- and right-handed helices. (left) Averages of \(\mathbb{1}_{B}(X^{T})\) for initial conditions in bins in the first two dPCs computed from 20 long (5 \(\mu\)s) trajectories. (middle) Averages of representative neural-network committors trained on the data set of 6,910 short (20 ns) trajectories; \(\tau\) corresponds to 400 ps. (right) Comparison between empirical committors (as defined in (45)) and the neural-network committors (trained as for the middle panel). Error bars indicate standard deviations over ten different initializations of the neural-network parameters.
predictions allows the analysis of transition paths using transition path theory without needing to generate full transition paths[37, 38, 48] and has been used to understand rare transition events in molecular dynamics[10, 13, 16, 17, 61, 62] and geophysical flows[63, 64, 65, 66, 67]. We leave these extensions to future work.

In cases in which trajectories that reach the reactant and product states are available, it would be interesting to compare our inexact iterative schemes against schemes for committor approximation based on logistic regression and related approaches [35, 41

Figure 16: \(\text{AIB}_{9}\) MFPT to the right-handed helix. (left) Averages of the time to next reach \(B\) for initial conditions in bins in the first two dPCs computed from 20 long (5 \(\mu\)s) trajectories. (middle) Averages of representative neural-network committors trained on the data set of 6,910 short (20 ns) trajectories; \(\tau\) corresponds to 400 ps. (right) Comparison between empirical committors (as defined in (47)) and the neural-network committors (trained as for the middle panel). Error bars indicate standard deviations over ten different initializations of the neural-network parameters.

Figure 15: MFPT between left- and right-handed helices for the \(\text{AIB}_{9}\) system. (top left) Convergence of Richardson iteration. (top right) Convergence of a five-dimensional subspace iteration. (bottom left) MFPT after 100 and 200 subspace iterations as a function of lag time. Shading indicates standard deviations over ten different initializations of the neural-network parameters. (bottom right) Overall error in MFPT.
46, 68]. These schemes are closely related to what is called "Monte-Carlo" approximation in reinforcement learning [31], and also to the inexact Richardson iteration that we propose here with \(\tau\to\infty\).

We have seen that temporal difference (TD) learning, a workhorse for prediction in reinforcement learning, is closely related to an inexact form of Richardson iteration. Variants like TD(\(\lambda\)), have similar relationships to inexact iterative schemes. As we showed, subspace iteration is a natural way of addressing slow convergence. We thus anticipate that our results have implications for reinforcement learning, particularly in scenarios in which the value function depends on the occurrence of some rare-event. Finally, we note that our framework can be extended to the wider range of iterative numerical linear algebra algorithms. In particular, Krylov or block Krylov subspace methods may offer further acceleration. In fact, very recently an approach along these lines was introduced for value function estimation in reinforcement learning[69].

## Author's contributions

J.S., S.C.G., A.R.D., and J.W. conceptualized research. J.S. developed the method. C.L. adapted the linear solve used for prediction from Refs. 58 and 59. J.S. and S.C.G. performed the numerical tests. A.R.D. and J.W. supervised the research. All authors wrote and edited the manuscript.

###### Acknowledgements.

 We are grateful to Arnaud Doucet for pointing our attention to TD methods and the inexact power iteration in Ref. 32, both of which were key motivations for this work. We also thank Joan Bruna for helpful conversations about reinforcement learning. This work was supported by National Institutes of Health award R35 GM136381 and National Science Foundation award DMS-2054306. S.C.G. acknowledges support by the National Science Foundation Graduate Research Fellowship under Grant No. 2140001. J.W. acknowledges support from the Advanced Scientific Computing Research Program within the DOE Office of Science through award DE-SC0020427. This work was completed in part with resources provided by the University of Chicago Research Computing Center, and we are grateful for their assistance with the calculations. "Beagle-3: A Shared GPU Cluster for Biomolecular Sciences" is supported by the National Institute of Health (NIH) under the High-End Instrumentation (HEI) grant program award 1S10OD028655-0.

## Data availability statement

The data that supports the findings of this study are available within the article. Code implementing the algorithms and a Jupyter notebook illustrating use of the method on the Muller-Brown example are available at [https://github.com/dinner-group/inexact-subspace-iteration](https://github.com/dinner-group/inexact-subspace-iteration).
