Inexact iterative numerical linear algebra for neural network-based spectral estimation and rare-event prediction

Strahan, John; Guo, Spencer C.; Lorpaiboon, Chatipat; Dinner, Aaron R.; Weare, Jonathan

doi:10.1063/5.0151309

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 dataset 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.

I. 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–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–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.^10,11

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–15 The performance of such a linear approximation depends on the choice of basis,^10,11,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¹⁴) 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.^5,18–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–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 datasets.

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.³¹ 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,³² 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.³³ 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 dataset. This opens the door to treating a wide range of observational, experimental, and computational datasets.

The remainder of the paper is organized as follows: In Sec. II, we describe the quantities that we seek to compute in terms of linear operators. In Secs. III and IV, we introduce an inexact subspace iteration algorithm that we use to solve for these quantities. Section V illustrates how the loss function can be tailored to the known properties of the desired quantity. Section VI 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 Sec. VII, we explain the details of the invariant subspace iteration and then demonstrate its application to our two examples. Finally, Sec. VIII 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.

II. 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

T^{t}

for a Markov process

X^{t} \in R^{d}

⁠, defined as

T^{t} f (x) = E_{x} [f (X^{t})],

(1)

where f is an arbitrary real-valued function and the subscript x indicates the initial condition X⁰ = x. The transition operator encodes how expectations of functions evolve in time. The right eigenfunctions of

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) = E_{x} [Ψ (X^{T}) + \sum_{t = 0}^{T - 1} Γ (X^{t})],

(2)

where T is the first time X^t ∉ D from some domain D, and Ψ and Γ 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) = E_{x} [T]

(3)

and the committor

q (x) = E_{x} [1_{B} (X^{T})] = P_{x} [X^{T} \in B],

(4)

where A and B are disjoint sets (“reactant” and “product” states) and

D = {(A \cup B)}^{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–36 and as a key ingredient in transition path theory statistics.^24,37,38 For any τ > 0, the prediction function u(x) satisfies the linear equation

(I - S^{τ}) u (x) = E_{x} [\sum_{t = 0}^{(τ \land T) - 1} Γ (X^{t})]

(5)

for x ∈ D, with boundary condition

u (x) = Ψ (x)

(6)

for x ∉ D. In ,

I

is the identity operator and

S^{τ} f (x) = E_{x} [f (X^{τ \land T})]

(7)

is the “stopped” transition operator.¹⁰ We use the notation τ ∧ T = min{τ, T}. The parameter τ is known as the lag time. While it is an integer corresponding to a number of integration steps, in our numerical examples, we often express it in terms of equivalent time units.

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 τ. A natural and generally applicable approach to finding an approximate solution to the prediction problem is to attempt to minimize the residual of over parameters θ of a neural network u_θ(x) representing u(x). For example, we recently suggested an algorithm that minimizes the residual norm

\frac{1}{2} {‖(I - S^{τ}) u_{θ} - r‖}_{μ}^{2},

(8)

where r(x) is the right-hand side of and μ is an arbitrary distribution of initial conditions X⁰ (boundary conditions were enforced by an additional term).²⁴ The norm ‖⋅‖_μ is defined by ‖ f‖_μ²=⟨ f,f⟩_μ, where ⟨ f,g⟩_μ= ∫f(x)g(x)μ(dx) for arbitrary functions f and g. The gradient of the residual norm in with respect to neural-network parameters θ can be written as

{⟨(I - S^{τ}) u_{θ} - r, \nabla_{θ} u_{θ}⟩}_{μ} - {⟨(I - S^{τ}) u_{θ} - r, S^{τ} \nabla_{θ} u_{θ}⟩}_{μ} .

(9)

Given a dataset of initial conditions

{X_{j}^{0}}_{j = 1}^{n}

drawn from μ and a single sample trajectory

{X_{j}^{t}}_{t = 0}^{τ}

of X^t for each

X_{j}^{0}

⁠, the first term in the gradient can be approximated without bias as

\begin{align} {⟨(I - S^{τ}) u_{θ} - r, \nabla_{θ} u_{θ}⟩}_{μ} \\ \approx \frac{1}{n} \sum_{j = 1}^{n} (u_{θ} (X_{j}^{0}) - u_{θ} (X_{j}^{τ \land T_{j}}) - \sum_{t = 0}^{(τ \land T_{j}) - 1} Γ (X_{j}^{t})) \nabla_{θ} u_{θ} (X_{j}^{0}), \end{align}

(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 $S^{τ}$ ⁠.^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).³¹ As in the semigradient approximation, the algorithms proposed in this paper only require access to inner products of the form ${⟨f,, A, g⟩}_{μ}$ for an operator $A$ related to $T^{τ}$ or $S^{τ}$ ⁠, and they avoid terms that are non-linear in $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.

III. AN INEXACT POWER ITERATION

To motivate the iterative numerical linear algebra point of view, observe that the first term in is the exact gradient with respect to θ′ of the loss

\frac{1}{2} {‖u_{θ^{'}} - S^{τ} u_{θ} - r‖}_{μ}^{2},

(11)

evaluated at θ′ = θ. The result of many steps of gradient descent (later, “inner iterations”) on this loss with θ held fixed can then be viewed as an inexact Richardson iteration for , resulting in a sequence

u_{θ^{s + 1}} \approx S^{τ} u_{θ^{s}} + r,

(12)

where, for each s,

u_{θ^{s}}

is a parameterized neural-network approximation of the solution to . To unify our discussion of the prediction and spectral estimation problems, it is helpful to observe that can be recast as an equivalent inexact power iteration

{\bar{u}}_{θ^{s + 1}} \approx A {\bar{u}}_{θ^{s}},

(13)

where

{\bar{u}}_{θ} = (\begin{matrix} 1 \\ u_{θ} \end{matrix}) and A = [\begin{matrix} 1 & 0 \\ r & S^{τ} \end{matrix}] .

(14)

Note that (1, u)^⊤ is the dominant eigenfunction of

A

and has eigenvalue equal to 1.

Reference 32 introduced an inexact power iteration to compute the stationary probability measure of $T^{τ}$ ⁠, 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}$ norm in (11) can be replaced by any other measure of the difference between u_θ′ and $S^{τ} u_{θ} + r$ ⁠, as long as an unbiased estimator of its gradient with respect to θ′ is available. This flexibility is discussed in more detail in Sec. V, 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 $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 $A$ in (14) is 1 and its next largest eigenvalue is the dominant eigenvalue of $S^{τ}$ ⁠. For rare-event problems, the difference between these two eigenvalues is expected to be very small. Indeed, when X⁰ is drawn from the quasi-stationary distribution of X^t in D, the logarithm of the largest eigenvalue of $S^{τ}$ is exactly $- τ / E [T]$ ⁠.³⁹ Thus, when the mean exit time from D is large, we can expect the spectral gap of $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 $S^{τ}$ with largest magnitude eigenvalues.³³ Iterative schemes that approximate multiple dominant eigenfunctions simultaneously, such as subspace iteration and Krylov methods, can converge much more rapidly.³³ In Sec. IV, 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 $A$ in (14), it can also be applied with $A = T^{τ}$ to compute multiple dominant eigenfunctions and values of $T^{τ}$ itself.

IV. AN INEXACT SUBSPACE ITERATION

Our inexact subspace iteration for the k dominant eigenfunctions/values of

A

proceeds as follows. Let

{\{φ_{θ^{s}}^{a}\}}_{a = 1}^{k}

be a sequence of k basis functions parameterized by θ^s (these can be scalar or vector valued functions depending on the form of

A

⁠). We can represent this basis as the vector valued function

U_{θ^{s}} = (φ_{θ^{s}}^{1}, φ_{θ^{s}}^{2}, \dots, φ_{θ^{s}}^{k}) .

(15)

Then, we can obtain a new set of k basis functions by approximately applying

A

to each of the components of

U_{θ^{s}}

⁠,

U_{θ^{s + 1}} K^{s + 1} \approx A U_{θ^{s}},

(16)

where K^s+1 is an invertible, upper-triangular k × k matrix that does not change the span of the components of

U_{θ^{s + 1}}

but is included to facilitate training. One way the approximate application of

A

can be accomplished is by minimizing

\frac{1}{2} \sum_{a = 1}^{k} {‖\sum_{b = 1}^{k} φ_{θ}^{b} K_{b a} - A φ_{θ^{s}}^{a}‖}_{μ}^{2}

(17)

over θ and K with θ^s held fixed.

The eigenvalues and eigenfunctions of

A

are then approximated by solving the finite-dimensional generalized eigenproblem

C^{τ} W = C^{0} W Λ,

(18)

where

C_{a b}^{τ} = {〈 φ_{θ^{s}}^{a}, A φ_{θ^{s}}^{b} 〉}_{μ},

(19)

C_{a b}^{0} = {〈 φ_{θ^{s}}^{a}, φ_{θ^{s}}^{b} 〉}_{μ},

(20)

each inner product is estimated using trajectory data, and W and Λ are k × k matrices. The matrix Λ is diagonal, and the ath eigenvalue λ_a of

A

is approximated by Λ_aa; the corresponding eigenfunction v_a is approximated by

\sum_{b = 1}^{k} W_{a b} φ_{θ^{s}}^{b}

⁠.

Even when sampling is not required to estimate the matrices in (19) and (20), the numerical rank of C^τ 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⁰ and C^τ 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

A

with a multiple of the identity, i.e., we perform our inexact subspace iteration on the operator

(1 - α_{s}) I + α_{s} A

in place of

A

itself. This new operator has the same eigenfunctions as

A

⁠. In our experiments, decreasing the parameter α_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

α_{s} = \{\begin{cases} 1 & s < σ \\ 1 / \sqrt{s + 1 - σ} & s \geq σ, \end{cases}

(21)

where σ 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 $(A = T^{τ})$ are given in Sec. VII and Algorithm 1. For the prediction problem with $A$ as in (14), they are given in Sec. VIII and Algorithm 2.

ALGORITHM 1.

Inexact subspace iteration (with $L_{μ}^{2}$ loss) for spectral estimation.

Require: Subspace dimension k, transition data

{X_{j}^{0}, X_{j}^{τ}}_{j = 1}^{n}

⁠, batch size B, learning rate η, number of subspace iterations S, number of inner iterations M, regularization parameters γ₁ and γ₂

1: Initialize

{φ_{θ}^{a}}_{a = 1}^{k}

and

{\{{\tilde{φ}}_{1}^{a}\}}_{a = 1}^{k}

2: for s = 1 … S do

3: for m = 1 … M do

4: Sample a batch of data

{X_{j}^{0}, X_{j}^{τ}}_{j = 1}^{B}

5:

{\hat{L}}_{1} \leftarrow \frac{1}{B} \sum_{j = 1}^{B} \sum_{a = 1}^{k} [\frac{1}{2} {(\sum_{b = 1}^{k} φ_{θ}^{b} (X_{j}^{0}) K_{b a})}^{2} - \sum_{b = 1}^{k} φ_{θ}^{b} (X_{j}^{0}) K_{b a} \{α_{s} {\tilde{φ}}_{s}^{a} (X_{j}^{τ}) + (1 - α_{s}) {\tilde{φ}}_{s}^{a} (X_{j}^{0})\}]

6:

{\hat{L}}_{K} \leftarrow γ_{1} {‖ K - diag (K) ‖}_{F}^{2}

7:

{\hat{L}}_{norm} \leftarrow γ_{2} \sum_{a = 1}^{k} (2 ν_{a} (\frac{1}{B} \sum_{j = 1}^{B} φ_{θ}^{a} {(X_{j}^{0})}^{2} - 1) - ν_{a}^{2})

8:

\hat{L} \leftarrow {\hat{L}}_{1} + {\hat{L}}_{K} + {\hat{L}}_{norm}

9:

θ \leftarrow θ - η \nabla_{θ} \hat{L}

10:

K \leftarrow K - η t r i u (\nabla_{K} \hat{L})

11:

ν_{a} \leftarrow ν_{a} + η \nabla_{ν_{a}} \hat{L}

12: end for

13: Compute the matrix

Φ_{i a} = φ_{θ}^{a} (X_{i}^{0})

⊳

Φ \in R^{n \times k}

14: Compute diagonal matrix

N_{a a}^{2} = \sum_{i} φ_{θ}^{a} {(X_{i}^{0})}^{2}

15: Compute QR-decomposition Φ = QR ⊳

Q \in R^{n \times k}

and

R \in R^{k \times k}

16:

{\tilde{φ}}_{s}^{a} \leftarrow \sum_{b = 1}^{k} φ_{θ}^{b} {(R^{- 1} N)}_{b a}

17:

K_{1 : i, i} \leftarrow \arg \min_{c} \frac{1}{n} \sum_{j = 1}^{n} (\sum_{a = 1}^{i} φ_{θ}^{a} (X_{j}^{0}) c_{a} - {\tilde{φ}}_{s}^{a} {(X_{j}^{τ}))}^{2} + γ_{2} \sum_{a = 1}^{i - 1} c_{a}^{2}

18: end for

19: Compute the matrices

C_{a b}^{t} = \frac{1}{n} \sum_{j = 1}^{n} {\tilde{φ}}_{s}^{a} (X_{j}^{0}) {\tilde{φ}}_{s}^{b} (X_{j}^{t})

for t = 0, τ ⊳

C^{t} \in R^{k \times k}

20: Solve the generalized eigenproblem C^τW = C⁰WΛ for W and Λ

21: return eigenvalues Λ, eigenfunctions

v_{a} = \sum_{b = 1}^{k} W_{a b} {\tilde{φ}}_{s}^{b}

ALGORITHM 2.

Inexact subspace iteration (with $L_{μ}^{2}$ loss) for prediction functions.

Require: Subspace dimension k, stopped transition data

{\{X_{j}^{0}, X_{j}^{τ \land T_{j}}\}}_{j = 1}^{n}

⁠, reward data

{\{ρ (X_{j})\}}_{j = 1}^{n}

⁠, batch size B, learning rate η, number of subspace iterations S, number of inner iterations M, regularization parameters γ₁ and γ₂

1: Initialize

{φ_{θ}^{a}}_{a = 1}^{k}

and

{\{{\tilde{φ}}_{1}^{a}\}}_{a = 1}^{k}

2: for s = 1 … S do

3: for m = 1 … M do

4: Sample a batch of data

{X_{j}^{0}, X_{j}^{τ \land T_{j}}}_{j = 1}^{B}

⁠,

{\{ρ (X_{j})\}}_{j = 1}^{B}

5:

{\hat{L}}_{1} \leftarrow \frac{1}{B} \sum_{j = 1}^{B} \sum_{a = 1}^{k} [\frac{1}{2} {(\sum_{b = 1}^{k} φ_{θ}^{b} (X_{j}^{0}) K_{b a})}^{2} - α_{s} (\sum_{b = 1}^{k} φ_{θ}^{b} K_{b a} ({\tilde{φ}}_{s}^{a} (X_{j}^{τ \land T_{j}}) + ρ (X_{j}) I_{a 1}))]

6:

{\hat{L}}_{2} \leftarrow - \frac{1}{B} \sum_{j = 1}^{B} \sum_{a = 1}^{k} (1 - α_{s}) \sum_{b = 1}^{k} φ_{θ}^{b} K_{b a} {\tilde{φ}}_{s}^{a} (X_{j}^{0})

7:

{\hat{L}}_{K} \leftarrow γ_{1} {‖ K - diag (K) ‖}_{F}^{2}

8:

{\hat{L}}_{norm} \leftarrow γ_{2} \sum_{a = 2}^{k} (2 ν_{a} (\frac{1}{B} \sum_{j = 1}^{B} φ_{θ}^{a} {(X_{j}^{0})}^{2} - 1) - ν_{a}^{2})

9:

\hat{L} \leftarrow {\hat{L}}_{1} + {\hat{L}}_{2} + {\hat{L}}_{K} + {\hat{L}}_{norm}

10:

θ \leftarrow θ - η \nabla_{θ} \hat{L}

11:

K \leftarrow K - η (t r i u (\nabla_{K} \hat{L}))

12:

ν_{a} \leftarrow ν_{a} + η \nabla_{ν_{a}} \hat{L}

13: end for

14: if Ψ(x) = 0 then

15: Compute the matrix

Φ_{i a} = φ_{θ}^{a} (X_{i}^{0})

⊳

Φ \in R^{n \times k}

16: else

17: Compute the matrix

Φ_{i a} = φ_{θ}^{a} (X_{i}^{0})

for a > 1 ⊳

Φ \in R^{n \times (k - 1)}

18: end if

19: Compute QR-decomposition Φ = QR

20: Compute diagonal matrix

N_{a a}^{2} = \sum_{i} φ_{θ}^{a} {(X_{i}^{0})}^{2}

21:

{\tilde{φ}}_{s}^{a} \leftarrow \sum_{b = 1}^{k} φ_{θ}^{b} {(R^{- 1} N)}_{b a}

⊳ if Ψ(x) = 0 exclude a = 1

22:

K_{1 : i, i} \leftarrow \arg \min_{c} \frac{1}{n} \sum_{j = 1}^{n} (\sum_{a = 1}^{i} φ_{θ}^{a} (X_{j}^{0}) c_{a} - {\tilde{φ}}_{s}^{a} {(X_{j}^{τ \land T_{j}}))}^{2} + γ_{2} \sum_{a = 1}^{i - 1} c_{a}^{2}

23: end for

24: Compute the matrix

C_{a b}^{t} = \frac{1}{n} \sum_{j = 1}^{n} φ_{θ}^{a} (X_{j}^{0}) φ_{θ}^{b} (X_{j}^{t})

for t = 0, τ ∧ T_j ⊳

C^{t} \in R^{k \times k}

25: Compute the vector

p_{a} = \frac{1}{n} \sum_{j = 1}^{n} φ_{θ}^{a} (X_{j}^{0}) ρ (X_{j})

26: Solve linear system (C⁰ − C^τ)w = p ⊳ if Ψ(x) = 0 enforce w₁ = 1

27: return

u = \sum_{a = 1}^{k} w_{a} φ_{θ}^{a}

V. ALTERNATIVE LOSS FUNCTIONS

As mentioned above, the inexact application of the operator $A$ can be accomplished by minimizing loss functions other than (17). The key requirement for a loss function in the present study is that $A$ appears in its gradient only through terms of the form ${⟨ f, A g ⟩}_{μ}$ 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_θ = ζ(z_θ), where ζ is an increasing function, and z_θ is a function parameterized by a neural network. An advantage of doing so is that the function ζ can restrict the output values of u_θ to some range. For example, when computing a probability such as the committor, a natural choice is the sigmoid function $ζ (x) = {(1 + e^{- x})}^{- 1}$ ⁠.

Our goal is to train a sequence of parameter values so that

u_{θ^{s}}

approximately follows a subspace iteration, i.e., so that

ζ (z_{θ^{s + 1}}) \approx A u_{θ^{s}}

⁠. To this end, we minimize with respect to θ the loss function

E_{X^{0} \sim μ} [V (z_{θ}) - z_{θ} A u_{θ^{s}}],

(22)

where V is an antiderivative of ζ, and θ^s is fixed. The subscript X⁰ ∼ μ in this expression indicates that X⁰ is drawn from μ. Note that, as desired,

A

appears in the gradient of only in an inner product of the required form, and an approximate minimizer, θ^s+1, of this loss satisfies

ζ (z_{θ^{s + 1}}) \approx A u_{θ^{s}}

⁠. This general form of loss function is adapted from variational expressions for the divergence of two probability measures.^32,40

The $L_{μ}^{2}$ loss in (17), which we use in several of our numerical experiments, corresponds to the choice ζ(x) = x and V(x) = x²/2. The choice of $ζ (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 τ 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 datasets of trajectories integrated all the way until reaching A or B to estimate the committor function.^41–46

VI. 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.

A. Müller-Brown potential

The first system is defined by the Müller–Brown potential⁴⁷ (Fig. 1)

\begin{align} V_{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}] . \end{align}

(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 Secs. VII A and VIII A, 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}, and z_i = {0, 0.5, 1.5, 1}. In Sec. VIII B, 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}, and 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 Fig. 1; these states are defined as

\begin{align} A = \{y, z : 6.5 {(y + 0.5)}^{2} - 11 (y + 0.5) (z - 1.5) + 6.5 {(z - 1.5)}^{2} < 0.3\} \\ B = \{y, z : {(y - 0.6)}^{2} + 0.5 {(z - 0.02)}^{2} < 0.2\} . \end{align}

(24)

FIG. 1.

View large Download slide

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 β⁻¹.

To generate a dataset, we randomly draw 50000 initial conditions

X_{j}^{0}

uniformly from the region

Ω = \{y, z : - 1.5 < y < 1.0, - 0.5 < z < 2.5, V_{MB} (y, z) < 12\}

(25)

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

X_{j}^{t} = X_{j}^{t - 1} - δ_{t} \nabla V_{MB} (X_{j}^{t - 1}) + \sqrt{δ_{t} 2 β^{- 1}} ξ_{j}^{t},

(26)

where 1 ≤ t ≤ τ, the

ξ_{j}^{t}

are independent standard Gaussian random variables, and the timestep is δ_t = 0.001. We use an inverse temperature of β = 2, and we vary τ as indicated below (note that τ is an integer number of steps, but we present our results for the Müller–Brown model in terms of the dimensionless scale used for δ_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 ten replicate datasets.

To validate our results, we compare the neural-network results against grid-based references, computed as described in Sec. 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′ = (y ± δ_x, z) or (y, z ± δ_x) on a grid with spacings δ_x,

P (x^{'} | x) = \frac{1}{1 + e^{β [V (x^{'}) - V (x)]}}

(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 δ_x = 0.0125; the grid is a rectangle with y ∈ [−1.5, 1.0], and z ∈ [−0.5, 2.0]. The reference invariant subspaces are computed by diagonalizing the matrix

2 (P - I) / β δ_{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)}^{c}

by solving

diag ({\hat{1}}_{{(A \cup B)}^{c}}) (P - I) {\hat{q}}_{+} = - diag ({\hat{1}}_{{(A \cup B)}^{c}}) (P - I) {\hat{1}}_{B}

(28)

and for grid sites in A ∪ 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.

B. AIB₉ helix-to-helix transition

The second system is a peptide of nine α-aminoisobutyric acids (AIB₉; Fig. 2). Because AIB is achiral around its α-carbon atom, AIB₉ 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 AIB₉ poses a stringent test due to the presence of many metastable intermediate states.

FIG. 2.

Helix-to-helix transition of AIB9. (Left) Left- and right-handed helices, which we use as states 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 µ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 β−1 corresponding to a temperature of 300 K.

View large Download slide

Helix-to-helix transition of AIB₉. (Left) Left- and right-handed helices, which we use as states 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 µ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 β⁻¹ corresponding to a temperature of 300 K.

The states are defined in terms of the internal ϕ and ψ 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° centered at (41°, 43°), that is,

{(ϕ - 41^{°})}^{2} + {(ψ - 43^{°})}^{2} \leq {(25^{°})}^{2} .

Amino acids are classified as being in the “right” state using the same radius, but centered instead at (−41°, −43°). 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 AIB₉ in defining the states as these are typically more flexible.⁵⁰ The states can be resolved by projecting onto dihedral angle principal components (dPCs; Fig. 2, right) as described previously.⁵¹

Following a procedure similar to that described in Ref. 50, we generate a dataset of short trajectories. From each of the 691 starting configurations in Ref. 50, we simulate ten trajectories of duration 20 ns with initial velocities drawn randomly from a Maxwell–Boltzmann distribution at a temperature of 300 K. The short trajectory dataset thus contains 6910 trajectories, corresponding to a total sampling time of 138.2 µs. We use a timestep of 4 fs together with a hydrogen mass repartitioning scheme,⁵² and configurations are saved every 40 ps. We employ the AIB parameters from Forcefield_NCAA⁵³ and the GBNeck2 implicit-solvent model.⁵⁴ Simulations are performed with the Langevin integrator in OpenMM 7.7.0⁵⁵ using a friction coefficient of 1 ps⁻¹. To generate a reference for comparison to our results, we randomly select 20 configurations from the dataset above and, from each, run a simulation of 5 µs with the same simulation parameters. For all following tests on AIB₉, batches consist of pairs of frames separated by τ drawn randomly with replacement from the short trajectories (i.e., from all possible such pairs in the dataset). 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⁵⁶ and be used as neural-network inputs as well.

VII. 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 with

A = T^{τ}

⁠, we require estimators for inner products of the form

{⟨ f, T^{τ} g ⟩}_{μ}

⁠. For example, the gradient of the loss function involves inner products of the form

{⟨\nabla_{θ} φ_{θ}^{a}, T^{τ} φ_{θ}^{b}⟩}_{μ} .

(29)

For these, we use the unbiased data-driven estimator

{⟨ f, T^{τ} g ⟩}_{μ} \approx \frac{1}{n} \sum_{j = 1}^{n} f (X_{j}^{0}) g (X_{j}^{τ}) .

(30)

As discussed in Sec. IV, applying the operator

T^{τ}

repeatedly causes each basis function to converge to the dominant eigenfunction and leads to 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

Φ_{i a} = φ_{θ}^{a} (X_{i}^{0})

and then computing the factorization Φ = QR, where Q is an n × k matrix with orthogonal columns and R is an upper triangular k × k matrix. Finally, we set

{\tilde{φ}}_{s}^{a} = \sum_{b = 1}^{k} φ_{θ}^{b} {(R^{- 1} N)}_{b a}

⁠, where N is a diagonal matrix with entries equal to the norms of the columns of Φ (before orthogonalization). To ensure that the networks remain well-separated (i.e., the eigenvalues of C⁰ remain away from zero), we penalize large off-diagonal entries of K by adding to the loss

γ_{1} {‖ K - diag (K) ‖}_{F}^{2},

(31)

where γ₁ allows us to tune the strength of this term relative to others, and ‖·‖_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

γ_{2} \sum_{a = 1}^{k} [2 ν_{a} (\frac{1}{n} \sum_{j = 1}^{n} φ_{θ}^{a} {(X_{j}^{0})}^{2} - 1) - ν_{a}^{2}],

(32)

where we have introduced the conjugate variables ν_a, which we maximize with gradient ascent (or similar optimization). In general, our numerical experiments suggest that it is best to keep γ₁ and γ₂ 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} \frac{1}{n} \sum_{j = 1}^{n} (\sum_{a = 1}^{i} φ_{θ}^{a} (X_{j}^{0}) c_{a} - {\tilde{φ}}_{s}^{a} {(X_{j}^{τ}))}^{2} + γ_{2} \sum_{a = 1}^{i - 1} c_{a}^{2} .

(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⁵⁷ for all numerical tests. We summarize our procedure for spectral estimation in Algorithm 1.

A. Müller–Brown model

As a first test of our method, we compute the k = 3 dominant eigenpairs for the Müller–Brown model. Since we know that the dominant eigenfunction of the transition operator is the constant function v₁(y, z) = 1 with eigenvalue λ₁ = 1, we directly include this function in the basis as a non-trainable function, i.e., $φ_{θ}^{1} (y, z) = 1$ ⁠. To initialize ${\tilde{φ}}_{1}^{a}$ for each a > 1, we choose a standard Gaussian vector $(Y^{a}, Z^{a}) \in R^{2}$ ⁠, and set ${\tilde{φ}}_{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 Müller–Brown tests, batches of trajectories are drawn from the entire dataset with replacement. Other hyperparameters are listed in Table I.

TABLE I.

Parameter choices used in this work.

	Spectral estimation		Committor			MFPT
Hyperparameter	Müller-Brown	AIB₉	Müller-Brown	Modified Müller-Brown	AIB₉	AIB₉
Subspace dimension k	3	5	1	2, 1a	1	5
Input dimensionality	2	174	2	2	174	174
Hidden layers	6	6	6	6	6	6
Hidden layer width	64	128	64	64	150	150
Hidden layer nonlinearity	CeLU	CeLU	ReLU	ReLU	ReLU	ReLU
Output layer nonlinearity	None	tanh	Sigmoid/none	None	None	None
Outer iterations S	10	100	100	4 + 10a	100	300
Inner iterations M	5000	2000	200	5000	2000	1000
σ	2	50	50	0	50	0
Batch size B	2000	1024	5000	2000	1024	2000
Learning rate η	0.001	0.0001	0.001	0.001	0.001	0.001
γ₁	0.15	0.001	⋯	⋯	⋯	0.1
γ₂	0.01	0.01	⋯	⋯	⋯	10
Loss for $φ_{θ}^{1}$	$L_{μ}^{2}$	$L_{μ}^{2}$	$L_{μ}^{2}$ /softplus	Softplus	Softplus	$L_{μ}^{2}$
Loss for $φ_{θ}^{a}$ for a > 1	$L_{μ}^{2}$	$L_{μ}^{2}$	⋯	$L_{μ}^{2}$	⋯	$L_{μ}^{2}$

	Spectral estimation		Committor			MFPT
Hyperparameter	Müller-Brown	AIB₉	Müller-Brown	Modified Müller-Brown	AIB₉	AIB₉
Subspace dimension k	3	5	1	2, 1a	1	5
Input dimensionality	2	174	2	2	174	174
Hidden layers	6	6	6	6	6	6
Hidden layer width	64	128	64	64	150	150
Hidden layer nonlinearity	CeLU	CeLU	ReLU	ReLU	ReLU	ReLU
Output layer nonlinearity	None	tanh	Sigmoid/none	None	None	None
Outer iterations S	10	100	100	4 + 10a	100	300
Inner iterations M	5000	2000	200	5000	2000	1000
σ	2	50	50	0	50	0
Batch size B	2000	1024	5000	2000	1024	2000
Learning rate η	0.001	0.0001	0.001	0.001	0.001	0.001
γ₁	0.15	0.001	⋯	⋯	⋯	0.1
γ₂	0.01	0.01	⋯	⋯	⋯	10
Loss for $φ_{θ}^{1}$	$L_{μ}^{2}$	$L_{μ}^{2}$	$L_{μ}^{2}$ /softplus	Softplus	Softplus	$L_{μ}^{2}$
Loss for $φ_{θ}^{a}$ for a > 1	$L_{μ}^{2}$	$L_{μ}^{2}$	⋯	$L_{μ}^{2}$	⋯	$L_{μ}^{2}$

^a

Four subspace iterations with k = 2 followed by ten iterations with k = 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 λ₃ and λ₄, indicating that a three-dimensional subspace captures the dynamics of interest for this system.

FIG. 3.

View large Download slide

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 τ = 300 steps (i.e., $0.3 δ_{t}^{- 1}$ ⁠).

FIG. 4.

View large Download slide

Spectral estimation as a function of lag time (in units of $δ_{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 − λ₂). (Bottom right) Subspace distance between estimated and reference three-dimensional invariant subspaces.

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 − λ₂ is comparable for the two methods (Fig. 4, lower left). To compare two finite dimensional subspaces,

U

and

V

⁠, we define the subspace distance as⁴

d (U, V) = {‖(I - P_{U}) P_{V}‖}_{F},

(34)

where

P_{U}

and

P_{V}

denote the orthogonal projections onto

U

and

V

⁠, respectively, and ‖·‖_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.

B. AIB₉

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₉ 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 (ϕ and ψ) for the nine residues (for a total of 2 × 2 × 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 expect 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 Müller–Brown example, we use $φ_{θ}^{1} = 1$ and a random linear combination of coordinate functions to initialize ${\tilde{φ}}_{1}^{a}$ for a > 1. Other hyperparameters are listed in Table I. 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 × 1000 MSM transition matrix, which is nearly instantaneous. However, the time for neural-network training is negligible compared with the time to generate the dataset, which is the same for both approaches.

Taking the dihedral MSM as a reference, the Cartesian MSM systematically underestimates the eigenvalues (Fig. 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 λ₄ and λ₅ may contribute to the difficulty in estimating λ₅. In Fig. 6, we plot the first two non-trivial eigenfunctions (v₂ and v₃), which align with the axes of the dPC projection. The eigenfunction v₂ corresponds to the transition between the left- and right-handed helices; the eigenfunction v₃ is nearly orthogonal to v₂ 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₄ and v₅ are orthogonal to v₂ and v₃. The estimates for v₂ are in qualitative agreement for all lag times tested (Fig. 6 shows results for τ corresponding to 40 ps), but the subspace iteration results are less noisy for the shortest lag times. Moreover, the estimate for v₃ from subspace iteration agrees more closely with that from the dihedral MSM than does the estimate for v₃ from the Cartesian MSM. The subspace distance for v₂ and v₃ 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.

FIG. 5.

First five eigenvalues of the transition operator for AIB9 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.

View large Download slide

First five eigenvalues of the transition operator for AIB₉ 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.

FIG. 6.

First two non-trivial eigenfunctions of AIB9 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.

View large Download slide

First two non-trivial eigenfunctions of AIB₉ 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.

VIII. PREDICTION

Inexact subspace iteration for

A

in is equivalent to performing the inexact Richardson iteration in on the first basis function

φ_{θ}^{1}

and then performing an inexact subspace iteration for the operator

S^{τ}

on the rest of the basis functions. The iteration requires unbiased estimators of the forms

{⟨f,, S^{τ}, g⟩}_{μ} \approx \frac{1}{n} \sum_{j = 1}^{n} f (X_{j}^{0}) g (X_{j}^{τ \land T_{j}})

(35)

and

{⟨f,, r⟩}_{μ} \approx \frac{1}{n} \sum_{j = 1}^{n} f (X_{j}^{0}) \sum_{t = 0}^{(τ \land T_{j}) - 1} Γ (X_{j}^{t}),

(36)

where T_j is the first time

X_{j}^{t}

enters

D^{c}

and r is the right-hand side of , as previously.

The Richardson iterate, $φ_{θ}^{1}$ ⁠, must satisfy the boundary condition $φ_{θ}^{1} (x) = Ψ (x)$ for x ∉ D. The other basis functions should satisfy $φ_{θ}^{a} (x) = 0$ for x ∉ D. In practice, we enforce the boundary conditions by explicitly setting $φ_{θ}^{1} (x) = Ψ (x)$ and $φ_{θ}^{a} (x) = 0$ for a > 1 when x ∉ D.

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

u_{θ} = \sum_{a = 1}^{k} w_{a} φ_{θ}^{a}

(37)

by solving the k-dimensional linear system

(C^{0} - C^{τ}) w = p,

(38)

where, for a, b ≥ 1,

C_{a b}^{t} = ⟨φ_{θ}^{a}, S^{t} φ_{θ}^{b}⟩_{μ}

(39)

for t = 0, τ, and

p_{a} = ⟨φ_{θ}^{a}, E_{x} [ρ (X)]⟩_{μ} .

(40)

In , we introduce the notation

ρ (X) = \sum_{t = 0}^{(τ \land T) - 1} Γ (X^{t})

(41)

for use in Algorithm 2.

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

ρ (X) = φ_{θ}^{1} (X^{τ \land T}) - φ_{θ}^{1} (X^{0}) + \sum_{t = 0}^{(τ \land T) - 1} Γ (X^{t}) .

(42)

In this case, the corresponding approximate solution is

u_{θ} = φ_{θ}^{1} + \sum_{a = 2}^{k} w_{a} φ_{θ}^{a} .

(43)

This approximate solution corresponds to the one given by dynamical Galerkin approximation^10,11 with the basis

{φ_{θ}^{a}}_{a = 2}^{k}

and a “guess” function of

φ_{θ}^{1}

⁠.

When the boundary conditions are zero, the orthogonalization procedure and the matrix K are applied to all basis functions as in Sec. VII. When the boundary conditions are non-zero, the orthogonalization procedure is only applied to the basis functions ${φ_{θ}^{a}}_{a = 2}^{k}$ ⁠, and K_a1 = I_a1 the a1 element of the identity matrix. We summarize our procedure for prediction in Algorithm 2.

A. Müller–Brown committor

In this section, we demonstrate the use of our method for prediction by estimating the committor for the Müller–Brown model with a shallow intermediate basin at (−0.25, 0.5) (Fig. 1). Here, the sets A and B are defined as in Eq. (24) and T is the time of first entrance to $D^{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 $S^{τ}$ [the second largest of $A$ in (14)] computed from our grid-based reference scheme (Sec. VI A). Richardson iteration should converge geometrically in this eigenvalue,³³ and so, for moderate lag times, we can expect our method to converge in a few dozen iterations. To initialize the algorithm we choose ${\tilde{φ}}_{1}^{1} = 1_{B}$ ⁠. All other hyperparameters are listed in Table I.

FIG. 7.

View large Download slide

First eigenvalue of $S^{τ}$ [second of $A$ in (14)] for the Müller–Brown model as a function of lag time (in units of $δ_{t}^{- 1}$ ⁠). The gap between this eigenvalue and the dominant eigenvalue, which is one, determines the rate of convergence of the Richardson iteration.

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

{l o g i t}_{ε} (q) = \log (\frac{ε + q}{1 + ε - q})

(44)

for points outside A and B. This function amplifies the importance of values close to zero and one. We include ɛ because we want to assign only a finite penalty if the procedure estimates q to be exactly zero or one; we use ɛ = e⁻²⁰.

Results as a function of lag time are shown in Fig. 8. We see that the Richardson iterate is more accurate than the MSM for all but the shortest lag times. When using the $L_{μ}^{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 dataset are shown in Fig. 9 for a fixed lag time of $τ = 0.1 δ_{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.

FIG. 8.

View large Download slide

Committor prediction for the Müller–Brown system as a function of lag time (in units of $δ_{t}^{- 1}$ ⁠). (Left) Comparison of the inexact Richardson iteration using the $L_{μ}^{2}$ loss and an MSM with 400 states. (Right) Same comparison using the softplus loss in place of the $L_{μ}^{2}$ loss.

FIG. 9.

View large Download slide

Committor prediction for the Müller–Brown as a function of number of initial conditions for a fixed lag time of $τ = 0.1 δ_{t}^{- 1}$ ⁠. (Left) Comparison of inexact Richardson iteration using the $L_{μ}^{2}$ loss and an MSM with 400 states. (Right) Same comparison using the softplus loss in place of the $L_{μ}^{2}$ loss.

B. Accelerating convergence by incorporating eigenfunctions

As discussed in Sec. III, we expect Richardson iteration to converge slowly when the largest eigenvalue of $S^{τ}$ ⁠, λ₁, is close to 1. More precisely, the number of iterations required to reach convergence should scale with $- 1 / \log λ_{1} = E [T] / τ$ ⁠, 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 Müller–Brown to perform poorly if we deepen the intermediate basin at (−0.25, 0.5) as in Fig. 10 (top left). Again, the sets A and B are defined as in (24), and T is the time of first entrance to $D^{c} = A \cup B$ ⁠. In this case, −1/log λ₁ is on the order of 100 for the lag times we consider and, as expected, inexact Richardson iteration converges slowly (Fig. 10, bottom left). Estimates of the committor by inexact Richardson iteration do not reach the correct values even after hundreds of iterations (Fig. 10, bottom right).

FIG. 10.

View large Download slide

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 β⁻¹. (Top right) Reference committor. (Bottom left) Dominant eigenvalue as a function of lag time (in units of $δ_{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.

We now show that convergence can be accelerated dramatically by incorporating additional eigenfunctions of $S^{τ}$ (i.e., k > 1 in Algorithm 2). For the Müller–Brown model with a deepened intermediate basin, the second eigenvalue of $S^{τ}$ is of order 10⁻⁴ for a lag time of τ = 1000 steps or $1 {δ_{t}}^{- 1}$ (while the first is near one as discussed above). We therefore choose k = 2 with ${\tilde{φ}}_{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 three memory kernels, corresponding to τ_mem = ⌊τ/4⌋.

ALGORITHM 3.

Memory-corrected linear solve for predictions.

Require: Stopped transition data

{\{X_{j}^{0}, X_{j}^{1 \land T_{j}}, \dots, X_{j}^{τ \land T_{j}}\}}_{j = 1}^{n}

⁠, guess function h, reward data

{\{ρ (X_{j})\}}_{j = 1}^{n}

⁠, basis set

{f^{a}}_{a = 1}^{k}

⁠, lag between memory kernels τ_mem.

1: for s = 0 … (τ/τ_mem) do

2: Initialize the matrix C^s with zeros ⊳

C^{s} \in R^{(k + 1) \times (k + 1)}

3:

C_{11}^{s} \leftarrow 1

4: for a = 2 … k do

5:

C_{a 1}^{s} \leftarrow \frac{1}{n} \sum_{j = 1}^{n} f^{a} (X_{j}^{0}) ρ (X_{j}^{s τ_{mem} \land T_{j}})

6: for b = 2 … k do

7:

C_{a b}^{s} \leftarrow \frac{1}{n} \sum_{j = 1}^{n} f^{a} (X_{j}^{0}) f^{b} (X_{j}^{s τ_{mem} \land T_{j}})

8: end for

9: end for

10: end for

11: A ← C¹ – C⁰

12: for s = 0 … (τ/τ_mem) − 2 do

13:

M^{s} \leftarrow C^{s + 2} - C^{s + 1} - A {(C^{0})}^{- 1} C^{s + 1} - \sum_{j = 0}^{s} M^{j} {(C^{0})}^{- 1} C^{s - j}

14: end for

15:

A_{mem} \leftarrow A + \sum_{s = 0}^{(τ / τ_{mem}) - 2} M^{s}

16: Solve A_memw = w

17: return

u = h + \sum_{a = 2}^{k} w_{a} f^{a}

The bottom row of Fig. 11 illustrates the idea of the subspace iteration. The second eigenfunction (Fig. 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 (Fig. 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.

FIG. 11.

View large Download slide

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 τ = 1000 steps (i.e., $1 {δ_{t}}^{- 1}$ ⁠).

FIG. 12.

View large Download slide

Committor for the Müller–Brown potential with deepened intermediate as a function of lag time (in units of ${δ_{t}}^{- 1}$ ⁠). (Left) Comparison of RMSE for subspace iteration as described above, Richardson iteration (as in Sec. VIII A but instead with 500 subspace iterations), and an MSM with 400 states. (Right) RMSE of the logit function in (44).

C. AIB₉ 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 AIB₉ 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^{c} = A \cup B$ with A and B defined in Sec. VI B. As before, we initialize ${\tilde{φ}}_{1}^{1} = 1_{B}$ ⁠.

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

\bar{q} (s) = P [X^{T} \in B | u_{θ} (X^{0}) \in [s, s + Δ s]]

(45)

for s ∈ [0, 1 − Δs]. Here, we use Δs = 0.05. The overall error in the committor estimate is defined as

q e r r o r = {(Δ s \sum_{n = 0}^{1 / Δ s - 1} {[\bar{q} (n Δ s) - n Δ s]}^{2})}^{1 / 2} .

(46)

While this measure of error can only be used when the dataset contains trajectories of long enough duration to reach

D^{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 dataset with τ corresponding to 400 ps are shown in Fig. 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 τ decreases, the results become less accurate (Fig. 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 dataset by subsampling the short trajectories and then training neural networks on the reduced set of trajectories (Fig. 14, top right). We find that the performance steadily drops as the number of trajectories is reduced and degrades rapidly for the datasets subsampled more than 20-fold (Fig. 14, bottom right), corresponding to about 7 µs of total sampling.

FIG. 13.

AIB9 committor for the transition between left- and right-handed helices. (Left) Averages of 1B(XT) for initial conditions in bins in the first two dPCs computed from 20 long (5 µs) trajectories. (Middle) Averages of representative neural-network committors trained on the dataset of 6910 short (20 ns) trajectories; τ 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.

View large Download slide

AIB₉ committor for the transition between left- and right-handed helices. (Left) Averages of $1_{B} (X^{T})$ for initial conditions in bins in the first two dPCs computed from 20 long (5 µs) trajectories. (Middle) Averages of representative neural-network committors trained on the dataset of 6910 short (20 ns) trajectories; τ 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.

FIG. 14.

AIB9 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 τ 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 datasets with (left) 1/2 and (right) 1/20 of the short trajectories. Error bars indicate standard deviations over ten random samples of the trajectories.

View large Download slide

AIB₉ 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 τ 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 datasets with (left) 1/2 and (right) 1/20 of the short trajectories. Error bars indicate standard deviations over ten random samples of the trajectories.

Finally, we compute the MFPT to reach the right-handed helix using the same dataset. For the MFPT calculation T is the time of first entrance to $D^{c} = B$ ⁠. Note that the time of first entrance to B includes long dwell times in A and is expected to be much larger than the time of first entrance to A ∪ B.

We compare against an empirical estimate of the MFPT defined by

\bar{m} (s) = E [T | u_{θ} (X^{0}) \in [s, s + Δ s]]

(47)

for s ∈ [0, m_max − Δs], where Δs = 3 and m_max = 57 ns. Overall error is defined analogously to Eq. .

In Fig. 15, we show the MFPT obtained from Algorithm 2 with k = 5 and the $L_{μ}^{2}$ loss function. Initially, we set ${\tilde{φ}}_{1}^{1}$ equal to an arbitrary positive function (we use $5 1_{A}$ ⁠) and ${\tilde{φ}}_{s}^{a}$ for a > 1 to a random linear combination of coordinate functions. In Fig. 16, we examine the convergence of the MFPT from the left-handed helix to the right-handed helix for the MFPT computed with k = 1 (pure Richardson iteration) and k = 5. To obtain the results shown, we train neural networks on the short-trajectory dataset and then average their predictions for structures in the left-handed helix state in the long reference trajectories. The horizontal line in Fig. 16 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 $S^{τ}$ to be very near to one when $D = B^{c}$ ⁠. In contrast, the committor calculation benefits from the fact that the time of first entrance to A ∪ B is much shorter, implying a smaller dominant eigenvalue of $S^{τ}$ when $D = {(A \cup B)}^{c}$ ⁠.

FIG. 15.

AIB9 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 µs) trajectories. (Middle) Averages of representative neural-network committors trained on the dataset of 6910 short (20 ns) trajectories; τ 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.

View large Download slide

AIB₉ 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 µs) trajectories. (Middle) Averages of representative neural-network committors trained on the dataset of 6910 short (20 ns) trajectories; τ 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.

FIG. 16.

MFPT between left- and right-handed helices for the AIB9 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. To obtain the results shown in this figure, we first use the short-trajectory dataset to train neural networks to predict the MFPT; we then use these networks with fixed parameters to evaluate the MFPT for all structures in the long reference trajectories and average the results for those structures in the left-handed helix state. The horizontal lines in the top panels are obtained from averaging the time to the right-handed helix for the same structures.

View large Download slide

MFPT between left- and right-handed helices for the AIB₉ 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. To obtain the results shown in this figure, we first use the short-trajectory dataset to train neural networks to predict the MFPT; we then use these networks with fixed parameters to evaluate the MFPT for all structures in the long reference trajectories and average the results for those structures in the left-handed helix state. The horizontal lines in the top panels are obtained from averaging the time to the right-handed helix for the same structures.

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 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–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–46,68} These schemes are closely related to what is called “Monte Carlo” approximation in reinforcement learning,³¹ and also to the inexact Richardson iteration that we propose here with τ → ∞.

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(λ), 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.⁶⁹

ACKNOWLEDGMENTS

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 the National Institutes of Health (Award No. R35 GM136381) and the National Science Foundation (Award No. 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 No. 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” was supported by the National Institutes of Health (NIH) under the High-End Instrumentation (HEI) Grant Program Award No. 1S10OD028655-0.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author 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.

John Strahan: Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Spencer C. Guo: Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Chatipat Lorpaiboon: Methodology (equal); Writing – review & editing (equal). Aaron R. Dinner: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Jonathan Weare: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support 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 Müller–Brown example are available at https://github.com/dinner-group/inexact-subspace-iteration.

REFERENCES

1.

F.

Noé

and

F.

Nüske

, “

A variational approach to modeling slow processes in stochastic dynamical systems

,”

Multiscale Model. Simul.

11

(

2

),

635

–

655

(

2013

).

https://doi.org/10.1137/110858616

Google Scholar

Crossref

2.

F.

Nüske

,

B. G.

Keller

,

G.

Pérez-Hernández

,

A. S. J. S.

Mey

, and

F.

Noé

, “

Variational approach to molecular kinetics

,”

J. Chem. Theory Comput.

10

(

4

),

1739

–

1752

(

2014

).

https://doi.org/10.1021/ct4009156

Google Scholar

Crossref

PubMed

3.

S.

Klus

,

F.

Nüske

,

P.

Koltai

,

H.

Wu

,

I.

Kevrekidis

,

C.

Schütte

, and

F.

Noé

, “

Data-driven model reduction and transfer operator approximation

,”

J. Nonlinear Sci.

28

,

985

–

1010

(

2018

).

https://doi.org/10.1007/s00332-017-9437-7

Google Scholar

Crossref

4.

R. J.

Webber

,

E. H.

Thiede

,

D.

Dow

,

A. R.

Dinner

, and

J.

Weare

, “

Error bounds for dynamical spectral estimation

,”

SIAM J. Math. Data Sci.

3

(

1

),

225

–

252

(

2021

).

https://doi.org/10.1137/20m1335984

Google Scholar

Crossref

PubMed

5.

C.

Lorpaiboon

,

E. H.

Thiede

,

R. J.

Webber

,

J.

Weare

, and

A. R.

Dinner

, “

Integrated variational approach to conformational dynamics: A robust strategy for identifying eigenfunctions of dynamical operators

,”

J. Phys. Chem. B

124

(

42

),

9354

–

9364

(

2020

).

https://doi.org/10.1021/acs.jpcb.0c06477

Google Scholar

Crossref

PubMed

6.

R. T.

McGibbon

,

B. E.

Husic

, and

V. S.

Pande

, “

Identification of simple reaction coordinates from complex dynamics

,”

J. Chem. Phys.

146

(

4

),

044109

(

2017

).

https://doi.org/10.1063/1.4974306

Google Scholar

Crossref

PubMed

7.

L.

Busto-Moner

,

C.-J.

Feng

,

A.

Antoszewski

,

A.

Tokmakoff

, and

A. R.

Dinner

, “

Structural ensemble of the insulin monomer

,”

Biochemistry

60

(

42

),

3125

–

3136

(

2021

).

https://doi.org/10.1021/acs.biochem.1c00583

Google Scholar

Crossref

PubMed

8.

G.

Pérez-Hernández

,

F.

Paul

,

T.

Giorgino

,

G.

De Fabritiis

, and

F.

Noé

, “

Identification of slow molecular order parameters for Markov model construction

,”

J. Chem. Phys.

139

(

1

),

015102

(

2013

).

https://doi.org/10.1063/1.4811489

Google Scholar

Crossref

PubMed

9.

C. R.

Schwantes

and

V. S.

Pande

, “

Improvements in Markov State Model construction reveal many non-native interactions in the folding of NTL9

,”

J. Chem. Theory Comput.

9

(

4

),

2000

–

2009

(

2013

).

https://doi.org/10.1021/ct300878a

Google Scholar

Crossref

PubMed

10.

J.

Strahan

,

A.

Antoszewski

,

C.

Lorpaiboon

,

B. P.

Vani

,

J.

Weare

, and

A. R.

Dinner

, “

Long-time-scale predictions from short-trajectory data: A benchmark analysis of the trp-cage miniprotein

,”

J. Chem. Theory Comput.

17

(

5

),

2948

–

2963

(

2021

).

https://doi.org/10.1021/acs.jctc.0c00933

Google Scholar

Crossref

PubMed

11.

E. H.

Thiede

,

D.

Giannakis

,

A. R.

Dinner

, and

J.

Weare

, “

Galerkin approximation of dynamical quantities using trajectory data

,”

J. Chem. Phys.

150

(

24

),

244111

(

2019

).

https://doi.org/10.1063/1.5063730

Google Scholar

Crossref

PubMed

12.

W. C.

Swope

,

J. W.

Pitera

, and

F.

Suits

, “

Describing protein folding kinetics by molecular dynamics simulations. 1. Theory

,”

J. Phys. Chem. B

108

(

21

),

6571

–

6581

(

2004

).

https://doi.org/10.1021/jp037421y

Google Scholar

Crossref

13.

N.

Frank

,

C.

Schütte

,

E.

Vanden-Eijnden

,

L.

Reich

, and

T. R.

Weik

, “

Constructing the equilibrium ensemble of folding pathways from short off-equilibrium simulations

,”

Proc. Natl. Acad. Sci. U. S. A.

106

(

45

),

19011

–

19016

(

2009

).

https://doi.org/10.1073/pnas.0905466106

Google Scholar

Crossref

PubMed

14.

An Introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation

, edited by

G. R.

Bowman

,

V. S.

Pande

, and

N.

Frank

,

Volume 797 of Advances in Experimental Medicine and Biology

(

Springer Netherlands

,

Dordrecht

,

2014

).

Google Scholar

Crossref

15.

J.

Finkel

,

R. J.

Webber

,

D. S.

Abbot

,

E. P.

Gerber

, and

J.

Weare

, “

Learning forecasts of rare stratospheric transitions from short simulations

,”

Mon. Weather Rev.

149

(

11

),

3647

–

3669

(

2021

).

https://doi.org/10.1175/mwr-d-21-0024.1

Google Scholar

Crossref

16.

A.

Antoszewski

,

C.

Lorpaiboon

,

J.

Strahan

, and

A. R.

Dinner

, “

Kinetics of phenol escape from the insulin R₆ hexamer

,”

J. Phys. Chem. B

125

(

42

),

11637

–

11649

(

2021

).

https://doi.org/10.1021/acs.jpcb.1c06544

Google Scholar

Crossref

PubMed

17.

S. C.

Guo

,

R.

Shen

,

B.

Roux

, and

A. R.

Dinner

, “

Dynamics of activation in the voltage-sensing domain of Ci-VSP

,”

bioRxiv:10.1101/2022.12.19.521128

(

2022

).

https://doi.org/10.1101/2022.12.19.521128

Google Scholar

18.

G.

Andrew

,

R.

Arora

,

J.

Bilmes

, and

K.

Livescu

, “

Deep canonical correlation analysis

,” in

Proceedings of the 30th International Conference on Machine Learning

(

PMLR

,

2013

), pp.

1247

–

1255

.

Google Scholar

19.

A.

Mardt

,

L.

Pasquali

,

H.

Wu

, and

N.

Frank

, “

VAMPnets for deep learning of molecular kinetics

,”

Nat. Commun.

9

(

1

),

5

(

2018

).

https://doi.org/10.1038/s41467-018-06999-0

Google Scholar

Crossref

PubMed

20.

C.

Wehmeyer

and

F.

Noé

, “

Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics

,”

J. Chem. Phys.

148

(

24

),

241703

(

2018

).

https://doi.org/10.1063/1.5011399

Google Scholar

Crossref

PubMed

21.

B.

Lusch

,

J. N.

Kutz

, and

S. L.

Brunton

, “

Deep learning for universal linear embeddings of nonlinear dynamics

,”

Nat. Commun.

9

(

1

),

4950

(

2018

).

https://doi.org/10.1038/s41467-018-07210-0

Google Scholar

Crossref

PubMed

22.

W.

Chen

,

H.

Sidky

, and

A. L.

Ferguson

, “

Nonlinear discovery of slow molecular modes using state-free reversible VAMPnets

,”

J. Chem. Phys.

150

(

21

),

214114

(

2019

).

https://doi.org/10.1063/1.5092521

Google Scholar

Crossref

PubMed

23.

A.

Glielmo

,

B. E.

Husic

,

A.

Rodriguez

,

C.

Clementi

,

N.

Frank

, and

A.

Laio

, “

Unsupervised learning methods for molecular simulation data

,”

Chem. Rev.

121

,

9722

(

2021

).

https://doi.org/10.1021/acs.chemrev.0c01195

Google Scholar

Crossref

PubMed

24.

J.

Strahan

,

J.

Finkel

,

A. R.

Dinner

, and

J.

Weare

, “

Predicting rare events using neural networks and short-trajectory data

,”

J. Comput. Phys.

488

,

112152

(

2023

).

https://doi.org/10.1016/j.jcp.2023.112152

Google Scholar

Crossref

PubMed

25.

H.

Li

,

Y.

Khoo

,

Y.

Ren

, and

L.

Ying

, “

A semigroup method for high dimensional committor functions based on neural network

,” in

Proceedings of the 2nd Mathematical and Scientific Machine Learning Conference

(

PMLR

,

2022

), pp.

598

–

618

.

Google Scholar

26.

Y.

Khoo

,

J.

Lu

, and

L.

Ying

, “

Solving for high-dimensional committor functions using artificial neural networks

,”

Res. Math. Sci.

6

(

1

),

1

(

2018

).

https://doi.org/10.1007/s40687-018-0160-2

Google Scholar

Crossref

27.

Q.

Li

,

B.

Lin

, and

W.

Ren

, “

Computing committor functions for the study of rare events using deep learning

,”

J. Chem. Phys.

151

(

5

),

054112

(

2019

).

https://doi.org/10.1063/1.5110439

Google Scholar

Crossref

28.

B.

Roux

, “

String method with swarms-of-trajectories, mean drifts, lag time, and committor

,”

J. Phys. Chem. A

125

(

34

),

7558

–

7571

(

2021

).

https://doi.org/10.1021/acs.jpca.1c04110

Google Scholar

Crossref

PubMed

29.

B.

Roux

, “

Transition rate theory, spectral analysis, and reactive paths

,”

J. Chem. Phys.

156

(

13

),

134111

(

2022

).

https://doi.org/10.1063/5.0084209

Google Scholar

Crossref

PubMed

30.

G. M.

Rotskoff

,

A. R.

Mitchell

, and

E.

Vanden-Eijnden

, “

Active importance sampling for variational objectives dominated by rare events: Consequences for optimization and generalization

,” in

Proceedings of the 2nd Mathematical and Scientific Machine Learning Conference

(

PMLR

,

2022

), pp.

757

–

780

.

Google Scholar

31.

R. S.

Sutton

and

A. G.

Barto

,

Reinforcement Learning: An Introduction

,

2nd ed.

(

The MIT Press

,

Cambridge, MA

,

2018

).

Google Scholar

32.

J.

Wen

,

B.

Dai

,

L.

Li

, and

S.

Dale

, “

Batch stationary distribution estimation

,” in

Proceedings of the 37th International Conference on Machine Learning, ICML’20

(

JMLR.org

,

2020

), pp.

10203

–

10213

.

Google Scholar

33.

G. H.

Golub

and

C. F.

Van Loan

,

Matrix Computations

,

3rd ed.

(

The Johns Hopkins University Press

,

1996

).

Google Scholar

34.

R.

Du

,

V. S.

Pande

,

A. Y.

Grosberg

,

T.

Tanaka

, and

E. S.

Shakhnovich

, “

On the transition coordinate for protein folding

,”

J. Chem. Phys.

108

(

1

),

334

–

350

(

1998

).

https://doi.org/10.1063/1.475393

Google Scholar

Crossref

35.

A.

Ma

and

A. R.

Dinner

, “

Automatic method for identifying reaction coordinates in complex systems

,”

J. Phys. Chem. B

109

(

14

),

6769

–

6779

(

2005

).

https://doi.org/10.1021/jp045546c

Google Scholar

Crossref

PubMed

36.

S. V.

Krivov

, “

On reaction coordinate optimality

,”

J. Chem. Theory Comput.

9

(

1

),

135

–

146

(

2013

).

https://doi.org/10.1021/ct3008292

Google Scholar

Crossref

PubMed

37.

W.

E

and

E.

Vanden-Eijnden

, “

Transition-path theory and path-finding algorithms for the study of rare events

,”

Annu. Rev. Phys. Chem.

61

,

391

–

420

(

2010

).

https://doi.org/10.1146/annurev.physchem.040808.090412

Google Scholar

Crossref

PubMed

38.

E.

Vanden-Eijnden

, “

Transition path theory

,” in

Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology

(

Springer

,

2006

), Vol.

1

, pp.

453

–

493

.

Google Scholar

39.

P.

Collett

,

S.

Martinez

, and

J.

San Martin

, “

Quasi-stationary distributions

,” in

Probability and its Applications

(

Springer Berlin Heidelberg

,

2012

).

Google Scholar

Crossref

40.

X.

Nguyen

,

M. J.

Wainwright

, and

M. I.

Jordan

, “

Estimating divergence functionals and the likelihood ratio by convex risk minimization

,”

IEEE Trans. Inf. Theory

56

(

11

),

5847

–

5861

(

2010

).

https://doi.org/10.1109/tit.2010.2068870

Google Scholar

Crossref

41.

B.

Peters

and

B. L.

Trout

, “

Obtaining reaction coordinates by likelihood maximization

,”

J. Chem. Phys.

125

(

5

),

054108

(

2006

).

https://doi.org/10.1063/1.2234477

Google Scholar

Crossref

PubMed

42.

B.

Peters

,

G. T.

Beckham

, and

B. L.

Trout

, “

Extensions to the likelihood maximization approach for finding reaction coordinates

,”

J. Chem. Phys.

127

(

3

),

034109

(

2007

).

https://doi.org/10.1063/1.2748396

Google Scholar

Crossref

PubMed

43.

H.

Jung

,

R.

Covino

, and

G.

Hummer

, “

Artificial intelligence assists discovery of reaction coordinates and mechanisms from molecular dynamics simulations

,” arXiv:1901.04595 (

2019

).

Google Scholar

44.

A.

Chattopadhyay

,

E.

Nabizadeh

, and

P.

Hassanzadeh

, “

Analog forecasting of extreme-causing weather patterns using deep learning

,”

J. Adv. Model. Earth Syst.

12

(

2

),

e2019MS001958

(

2020

).

https://doi.org/10.1029/2019MS001958

Google Scholar

Crossref

PubMed

45.

H.

Jung

,

R.

Covino

,

A.

Arjun

,

C.

Leitold

,

C.

Dellago

,

P. G.

Bolhuis

, and

G.

Hummer

, “

Machine-guided path sampling to discover mechanisms of molecular self-organization

,”

Nat. Comput. Sci.

3

(

4

),

334

–

345

(

2023

).

https://doi.org/10.1038/s43588-023-00428-z

Google Scholar

Crossref

46.

M.

George

,

B.

Cozian

,

P.

Abry

,

P.

Borgnat

, and

F.

Bouchet

, “

Probabilistic forecasts of extreme heatwaves using convolutional neural networks in a regime of lack of data

,”

Phys. Rev. Fluids

8

(

4

),

040501

(

2023

).

https://doi.org/10.1103/PhysRevFluids.8.040501

Google Scholar

Crossref

47.

K.

Müller

and

L. D.

Brown

, “

Location of saddle points and minimum energy paths by a constrained simplex optimization procedure

,”

Theor. Chim. Acta

53

(

1

),

75

–

93

(

1979

).

https://doi.org/10.1007/BF00547608

Google Scholar

Crossref

48.

C.

Lorpaiboon

,

J.

Weare

, and

A. R.

Dinner

, “

Augmented transition path theory for sequences of events

,”

J. Chem. Phys.

157

(

9

),

094115

(

2022

).

https://doi.org/10.1063/5.0098587

Google Scholar

Crossref

PubMed

49.

S.

Buchenberg

,

N.

Schaudinnus

, and

G.

Gerhard Stock

, “

Hierarchical biomolecular dynamics: Picosecond hydrogen bonding regulates microsecond conformational transitions

,”

J. Chem. Theory Comput.

11

(

3

),

1330

–

1336

(

2015

).

https://doi.org/10.1021/ct501156t

Google Scholar

Crossref

PubMed

50.

A.

Perez

,

F.

Sittel

,

G.

Stock

, and

K.

Dill

, “

MELD-path efficiently computes conformational transitions, including multiple and diverse paths

,”

J. Chem. Theory Comput.

14

(

4

),

2109

–

2116

(

2018

).

https://doi.org/10.1021/acs.jctc.7b01294

Google Scholar

Crossref

PubMed

51.

F.

Sittel

,

T.

Filk

, and

G.

Stock

, “

Principal component analysis on a torus: Theory and application to protein dynamics

,”

J. Chem. Phys.

147

(

24

),

244101

(

2017

).

https://doi.org/10.1063/1.4998259

Google Scholar

Crossref

PubMed

52.

C. W.

Hopkins

,

R. C.

Walker

,

S.

Le Grand

, and

A. E.

Roitberg

, “

Long-Time-step molecular dynamics through hydrogen mass repartitioning

,”

J. Chem. Theory Comput.

11

(

4

),

1864

–

1874

(

2015

).

https://doi.org/10.1021/ct5010406

Google Scholar

Crossref

PubMed

53.

G. A.

Khoury

,

J.

Smadbeck

,

P.

Tamamis

,

A. C.

Vandris

,

C. A.

Kieslich

, and

C. A.

Floudas

, “

Ab Initio charge parameters to aid in the discovery and design of therapeutic proteins and peptides with unnatural amino acids and their application to complement inhibitors of the compstatin family

,”

ACS Synth. Biol.

3

(

12

),

855

–

869

(

2014

).

https://doi.org/10.1021/sb400168u

Google Scholar

Crossref

PubMed

54.

H.

Nguyen

,

D. R.

Roe

, and

C.

Simmerling

, “

Improved generalized Born solvent model parameters for protein simulations

,”

J. Chem. Theory Comput.

9

(

4

),

2020

–

2034

(

2013

).

https://doi.org/10.1021/ct3010485

Google Scholar

Crossref

PubMed

55.

P.

Eastman

,

J.

Swails

,

J. D.

Chodera

,

R. T.

McGibbon

,

Y.

Zhao

,

K. A.

Beauchamp

,

L.-P.

Wang

,

A. C.

Simmonett

,

M. P.

Harrigan

,

C. D.

Stern

,

R. P.

Wiewiora

,

B. R.

Brooks

, and

V. S.

Pande

, “

OpenMM 7: Rapid development of high performance algorithms for molecular dynamics

,”

PLoS Comput. Biol.

13

(

7

),

e1005659

(

2017

).

https://doi.org/10.1371/journal.pcbi.1005659

Google Scholar

Crossref

PubMed

56.

P.

Metzner

,

C.

Schütte

, and

E.

Vanden-Eijnden

, “

Illustration of transition path theory on a collection of simple examples

,”

J. Chem. Phys.

125

(

8

),

084110

(

2006

).

https://doi.org/10.1063/1.2335447

Google Scholar

Crossref

PubMed

57.

D. P.

Kingma

and

B.

Jimmy

, “

Adam: A method for stochastic optimization

,” arXiv:1412.6980 (

2014

).

Google Scholar

58.

E.

Darve

,

J.

Solomon

, and

A.

Kia

, “

Computing generalized Langevin equations and generalized Fokker-Planck equations

,”

Proc. Natl. Acad. Sci. U. S. A.

106

(

27

),

10884

–

10889

(

2009

).

https://doi.org/10.1073/pnas.0902633106

Google Scholar

Crossref

PubMed

59.

S.

Cao

,

A.

Montoya-Castillo

,

W.

Wang

,

T. E.

Markland

, and

X.

Huang

, “

On the advantages of exploiting memory in Markov state models for biomolecular dynamics

,”

J. Chem. Phys.

153

(

1

),

014105

(

2020

).

https://doi.org/10.1063/5.0010787

Google Scholar

Crossref

PubMed

60.

D.

Lucente

,

J.

Rolland

,

C.

Herbert

, and

F.

Bouchet

, “

Coupling rare event algorithms with data-based learned committor functions using the analogue Markov chain

,”

J. Stat. Mech.: Theory Exp.

2022

(

8

),

083201

(2022).

https://doi.org/10.1088/1742-5468/ac7aa7

Crossref

61.

Y.

Meng

,

D.

Shukla

,

V. S.

Pande

, and

B.

Roux

, “

Transition path theory analysis of c-Src kinase activation

,”

Proc. Natl. Acad. Sci. U. S. A.

113

(

33

),

9193

–

9198

(

2016

).

https://doi.org/10.1073/pnas.1602790113

Google Scholar

Crossref

PubMed

62.

B. P.

Vani

,

J.

Weare

, and

A. R.

Dinner

, “

Computing transition path theory quantities with trajectory stratification

,”

J. Chem. Phys.

157

(

3

),

034106

(

2022

).

https://doi.org/10.1063/5.0087058

Google Scholar

Crossref

PubMed

63.

J.

Finkel

,

D. S.

Abbot

, and

J.

Weare

, “

Path properties of atmospheric transitions: Illustration with a low-order sudden stratospheric warming model

,”

J. Atmos. Sci.

77

(

7

),

2327

–

2347

(

2020

).

https://doi.org/10.1175/jas-d-19-0278.1

Google Scholar

Crossref

64.

P.

Miron

,

F. J.

Beron-Vera

,

L.

Helfmann

, and

P.

Koltai

, “

Transition paths of marine debris and the stability of the garbage patches

,”

Chaos

31

(

3

),

033101

(

2021

).

https://doi.org/10.1063/5.0030535

Google Scholar

Crossref

PubMed

65.

D.

Lucente

,

C.

Herbert

, and

F.

Bouchet

, “

Committor functions for climate phenomena at the predictability margin: The example of El Niño-Southern Oscillation in the Jin and Timmermann model

,”

J. Atmos. Sci.

79

(

9

),

2387

–

2400

(

2022

).

https://doi.org/10.1175/jas-d-22-0038.1

Google Scholar

Crossref

66.

J.

Finkel

,

E. P.

Gerber

,

S. A.

Dorian

, and

J.

Weare

, “

Revealing the statistics of extreme events hidden in short weather forecast data

,”

AGU Adv.

4

(

2

),

e2023AV000881

(

2023

).

https://doi.org/10.1029/2023av000881

Google Scholar

Crossref

67.

J.

Finkel

,

R. J.

Webber

,

E. P.

Gerber

,

D. S.

Abbot

, and

J.

Weare

, “

Data-driven transition path analysis yields a statistical understanding of sudden stratospheric warming events in an idealized model

,”

J. Atmos. Sci.

80

(

2

),

519

–

534

(

2023

).

https://doi.org/10.1175/jas-d-21-0213.1

Google Scholar

Crossref

68.

J.

Hu

,

A.

Ma

, and

A. R.

Dinner

, “

A two-step nucleotide-flipping mechanism enables kinetic discrimination of DNA lesions by AGT

,”

Proc. Natl. Acad. Sci. U. S. A.

105

(

12

),

4615

–

4620

(

2008

).

https://doi.org/10.1073/pnas.0708058105

Google Scholar

Crossref

PubMed

69.

E.

Xia

and

M.

Wainwright

, “

Krylov-Bellman boosting: Super-linear policy evaluation in general state spaces

,” in

Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, Volume 206 of Proceedings of Machine Learning Research

, edited by

F.

Ruiz

,

D.

Jennifer

, and

J.-W.

van de Meent

(

PMLR

,

2023

), pp.

9137

–

9166

.

Google Scholar

2023

Author(s)

Inexact iterative numerical linear algebra for neural network-based spectral estimation and rare-event prediction

I. INTRODUCTION

II. SPECTRAL ESTIMATION AND THE PREDICTION PROBLEM

III. AN INEXACT POWER ITERATION

IV. AN INEXACT SUBSPACE ITERATION

V. ALTERNATIVE LOSS FUNCTIONS

VI. TEST PROBLEMS

A. Müller-Brown potential

B. AIB₉ helix-to-helix transition

VII. SPECTRAL ESTIMATION

A. Müller–Brown model

B. AIB₉

VIII. PREDICTION

A. Müller–Brown committor

B. Accelerating convergence by incorporating eigenfunctions

C. AIB₉ prediction results

IX. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Inexact iterative numerical linear algebra for neural network-based spectral estimation and rare-event prediction

I. INTRODUCTION

II. SPECTRAL ESTIMATION AND THE PREDICTION PROBLEM

III. AN INEXACT POWER ITERATION

IV. AN INEXACT SUBSPACE ITERATION

V. ALTERNATIVE LOSS FUNCTIONS

VI. TEST PROBLEMS

A. Müller-Brown potential

B. AIB9 helix-to-helix transition

VII. SPECTRAL ESTIMATION

A. Müller–Brown model

B. AIB9

VIII. PREDICTION

A. Müller–Brown committor

B. Accelerating convergence by incorporating eigenfunctions

C. AIB9 prediction results

IX. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

B. AIB₉ helix-to-helix transition

B. AIB₉

C. AIB₉ prediction results