Optimization of randomized neural networks for transfer operator approximation Mohammad Tabish∗1 and Stefan Klus2
arXiv:2605.23689v1 [cs.LG] 22 May 2026
1
Maxwell Institute for Mathematical Sciences, University of Edinburgh and Heriot–Watt University, Edinburgh, UK 2 School of Mathematical & Computer Sciences, Heriot–Watt University, Edinburgh, UK
Abstract RaNNDy is a randomized neural network architecture for the data-driven approximation of transfer operators associated with complex dynamical systems. The weights and biases of the hidden layers of the network are randomly initialized and kept fixed, only the output layer is trained. This has several advantages over fully optimized neural networks, notably a closed-form solution for the output layer and significantly lower training costs. Despite these advantages, RaNNDy is restricted to the initial selection of weights and biases that parametrize the basis functions required for the operator approximation. Since the basis functions are determined by the activation function, choosing an appropriate activation function for the hidden layers is crucial. In this work, we propose an algorithm that optimizes the activation function itself, while keeping the weights and biases in the randomized neural network fixed, providing a more suitable dictionary. We illustrate the efficacy of the approach using various benchmark problems, including stochastic differential equations and random walks on graphons.
1 Introduction Data-driven methods for the approximation of transfer operators have become an area of intense research, as they enable the analysis of the global behavior of complex dynamical systems from data. Spectral decompositions of these operators can be used to analyze the long-term behavior of dynamical systems, e.g., metastable states of molecules or coherent structures in fluid flows, see, e.g., [1] for an overview of different applications. Existing data-driven methods such as extended dynamic mode decomposition (EDMD) [2, 3] and the variational approach of conformation dynamics (VAC) [4, 5] learn a finite-dimensional representation of transfer operators using a fixed set of basis functions, also known as a dictionary. The dictionary could, for instance, comprise monomials, trigonometric functions, Gaussian functions, or indicator functions. The optimal choice of basis functions, however, is highly system-dependent and in general an open problem. Neural networks (NNs) have been used to represent the dictionary required to approximate the operators, and several architectures have been proposed in the literature, see, e.g., [6, 7, 8]. The training of these networks is performed using backpropagation, which iteratively minimizes a loss ∗
Corresponding author: [email protected]
1
function to update the networks’ parameters. This approach is computationally expensive for complex networks and can also lead to problems such as slow convergence, numerous local minima, exploding/vanishing gradients, and sensitivity to hyperparameters such as the learning rate [9, 10]. Randomized neural networks avoid using backpropagation by fixing the weights and biases of the hidden layers and training only the output layer for a particular task [11, 12], making the training computationally cheaper. Based on this idea, we recently proposed RaNNDy [13], a randomized neural network for the data-driven approximation of transfer operators. RaNNDy computes eigenfunctions of the operator directly via a closed-form solution for the output layer. Fixing the hidden layers’ weights and biases, however, leaves no scope for the optimization of the randomized basis. That is, if the randomly generated basis is not suitable, i.e., it is not able to represent the eigenfunctions of the operator accurately, it could lead to poor approximations. Hence, we propose an algorithm to optimize hyperparameters of parametrized activation functions while keeping the parameters of the model fixed. In other words, we optimize the network while maintaining the computational advantages of fixed weights and biases. We use the variational principle to optimize the activation parameters. The main contributions of this work are: 1. We propose a novel algorithm to optimize the randomized basis of RaNNDy by tuning the hyperparameters of the activation function, while keeping the network weights and biases fixed. 2. We present various applications, including the Bickley jet, high-dimensional protein folding processes, and random walks on graphons. The proposed algorithm can be regarded as a compromise between fully randomized and fully trained neural networks for transfer operator approximation. The paper is structured as follows: In Section 2, we will introduce transfer operators and randomized neural networks. In Section 3, we will discuss the proposed algorithm to optimize the randomized basis used by RaNNDy. Numerical results will be outlined in Section 4 and a conclusion will be presented in Section 5.
2 Background We start by briefly introducing transfer operators for reversible and non-reversible dynamical systems and randomized neural networks.
2.1 Transfer operators Transfer operators offer an alternative viewpoint for analyzing stochastic processes in terms of observables and probability densities instead of states of the system. The Koopman operator, for example, describes the evolution of observables and the Perron–Frobenius operator the evolution of the probability densities of the system. These operators are infinite-dimensional but linear representations of the underlying finite-dimensional but nonlinear dynamics. Since we cannot work with infinite-dimensional operators in practice, the goal is to find suitable finite-dimensional approximations of these operators. Let {Xt } be a stochastic process defined on the state space X ⊂ Rd that is governed by a stochastic differential equation (SDE) of the form dXt = b(Xt ) dt + σ(Xt ) dWt ,
2
(2.1)
where b : Rd → Rd is the drift term, σ : Rd → Rd×d is the diffusion term, and Wt is a d-dimensional Wiener process. In what follows, let pτ : X × X → R be the conditional probability of Xt+τ = y given Xt = x, where τ is a fixed lag time. Additionally, let Lr denote the space of (equivalence classes of) r-Lebesgue integrable functions and Lrµ the corresponding µ-weighted space, with 1 ≤ r ≤ ∞. Definition 2.1 (Transfer operators). Let τ > 0 be a fixed lag time. 1. The Perron–Frobenius operator P τ is defined by Z τ P ρ(x) = pτ (y, x) ρ(y) dy. X
2. The Koopman operator Kτ is defined by Z pτ (x, y)f (y) dy. Kτ f (x) = X
3. The reweighted Perron–Frobenius operator T τ that propagates densities w.r.t. a reference density µ is defined by Z 1 T τ u(x) = pτ (y, x)u(y)µ(y) dy, ν(x) X where ν = P τ µ. It is well known that the eigenvalues λi and the corresponding eigenfunctions φi of these transfer operators—or compositions thereof—can help us identify, for example, metastable states in molecular dynamics or coherent sets in fluid flows [14, 15, 16, 3, 17]. Definition 2.2 (Invariant density). The eigenfunction of P τ corresponding to the eigenvalue λ = 1 is called stationary density or invariant density, usually denoted by π, i.e., it satisfies P τ π = π. The operators introduced above are well-defined as P τ : L21/π → L21/π , T τ : L2π → L2π and Kτ : L2π → L2π , if π exists. We refer the reader to [18] for more details on the domains. Example 2.3. For the overdamped Langevin equation given by p dXt = −∇V (x) dt + 2β −1 dWt , where V (x) is a confining potential and β > 0 is the inverse temperature, there exists a unique invariant density given by π(x) = e−βV (x) . Definition 2.4 (Reversibility). The process Xt is called reversible if the so-called detailed balance condition is satisfied, i.e., π(x) pτ (x, y) = π(y) pτ (y, x) for all x, y ∈ X. For reversible processes, the Koopman operator Kτ and the Perron–Frobenius operator are self-adjoint with respect to suitably reweighted inner products [5]. However, for non-reversible dynamical systems, they are in general not self-adjoint. That is, their eigenvalues and eigenfunctions may be complex-valued. In this case, we approximate the singular values and singular functions of these operators to analyze the system. The dominant singular functions of these operators allow us to detect coherent sets in fluid flows [17].
Pτ
3
2.2 Randomized neural networks for transfer operators In randomized neural networks, the weights and biases of the hidden layers are randomly initialized and kept fixed, only the output layer is trained. The advantage is that the network can be trained without using backpropagation [12]. The hidden layers act as a random feature map (RFM) through which the input data are transformed into a random feature space, and the output is represented using the output layer. For a single hidden layer and an output layer, the model can be understood as a function f : Rd → Rn that produces n outputs for a d-dimensional input data point and is defined by f(x) = Wo · σ(W x + b), where σ is the activation function applied component-wise to the vector, W ∈ RN ×d are the weights of the hidden layer for N neurons, b ∈ RN is the bias term, and Wo ∈ Rn×N is the matrix containing the output weights for the n outputs. Training the network now involves finding the optimal output layer weights Wo [12, 11].
3 Optimizing the activation function in RaNNDy Consider a set of N basis functions ψi : X → R, with i = 1, . . . , N , written as a vectorvalued function ψ(x) = [ψ1 (x), ψ2 (x), . . . , ψN (x)]⊤ . Provided that these functions are linearly independent, the basis spans an N -dimensional subspace for the finite-dimensional approximation of the Koopman operator Kτ projected onto this subspace, denoted by Kψτ . RaNNDy utilizes the random feature map (RFM) of randomized neural networks as a basis and approximates the n dominant eigenfunctions of the operator, which are then represented by the output layer. RaNNDy is trained using the variational principle to find the weights of the output layer. For a given set of hidden layer weights and biases, the randomized basis can be written as ψ(x) = σ(W x + b). However, as mentioned above, fixing the hidden layers does not leave any scope for the optimization of the basis, i.e., for systems where the eigenfunctions cannot be represented by the fixed basis due to the poor initialization of the hidden layers’ weights and biases. We thus want to optimize the basis to obtain more accurate approximations of the eigenfunctions. We can provide flexibility in the selection of randomized basis functions in RaNNDy by introducing a set of tunable parameters in the selected activation function and the distribution of the weights and biases. These parameters could, for instance, be the scale of the distribution from which the hidden weights are sampled. That is, the randomized basis now depends on a parametric activation function σ, i.e., ψ(x, ω) = σ(W (ωW )x + b(ωb ), ωa ), where ω = [ωa , ωW , ωb ]⊤ contains the parameters to be optimized, W (ωW ) : RpW → RN ×d and b(ωb ) : Rpb → RN are the functions that depend on the parameters ωW and ωb , respectively. For example, ωW could be the scale of the standard normal distribution, etc., and ωa ∈ Rpa depends on the selected activation. For instance, a parametric version of the tanh activation function would be ψ(x, ω) = tanh(W (ωW )x + b(ωb )).
4
Self-adjoint Koopman operators To find the best hyperparameters ω for the randomized basis ψ(x, ω), we can define a loss function using the Rayleigh variational principle, proposed in [4], that provides a principled way to approximate spectral decompositions of transfer operators. Theorem 3.1. Let φ bi be an approximation of the ith true eigenfunction φi of the Koopman operator Kτ . Assuming that φ bi is normalized and orthogonal to the previous i − 1 eigenfunctions, we have ⟨Kτ φ bi , φ bi ⟩ ≤ λi . Hence, at each training epoch, for a given basis parametrized by ω, we approximate the eigenvalues of the projected operator and, using the above theorem, maximize Pn their sum to find the best approximation of the eigenfunctions. For a function f (x) = i=1 wi ψi (x, ω), we can optimize the Rayleigh quotient max
⟨f, Kψτ f ⟩ ⟨f, f ⟩
f
= max w
w⊤ C01 (ω)w , w⊤ C00 (ω)w
where [C01 (ω)]ij = ⟨ψi ( · , ω), Kψτ ψj ( · , ω)⟩,
[C00 (ω)]ij = ⟨ψi ( · , ω), ψj ( · , ω)⟩.
The optimal solutions w of the above problem are given by the dominant eigenvectors of −1 (ω) C01 (ω) of the projected operator Kψτ , see [4, 13] the matrix representation A(ω) = C00 for details. These eigenvectors are orthonormal with respect to the inner product weighted by the matrix C00 (ω) [4]. Hence, the sum of the eigenvalues, i.e., n X
λi = max
ψ1 ,...,ψn
i=1
s.t.
n X ⟨ψi (x, ω), Kψτ ψi (x, ω)⟩, i=1
⟨ψi (x, ω), ψj (x, ω)⟩ = δij ,
can be optimized to tune the activation parameters ω using, for example, gradient descent or Bayesian optimization [19]. After each training step, we get an optimized set of basis function parameters ω. In general, we cannot compute the required integrals in the above equations and thus have to approximate them from data. For a fixed lag time τ , let xi , yi ∈ X denote the states of the system, where xi = Xt and yi = Xt+τ . The data can be obtained by direct measurements or simulations, for example, by numerically solving (2.1). The obtained data can then be stored in the matrices X, Y ∈ Rd×m , defined by X = [x1 , x2 , . . . , xm ]
and
Y = [y1 , y2 , . . . , ym ].
The matrices C00 (ω) and C01 (ω) can be approximated from data using b00 (ω) = 1 C m b01 (ω) = C
m X
1 Ψ0 (ω) Ψ0 (ω)⊤ , m
ψ(xi , ω)ψ(yi , ω)⊤ =
1 Ψ0 (ω) Ψ1 (ω)⊤ , m
i=1
m 1 X
m
ψ(xi , ω) ψ(xi , ω)⊤ =
i=1
5
(3.1)
b where Ψ0 (ω) = ψ(X, ω) ∈ RN ×m and Ψ1 (ω) = ψ(Y, ω) ∈ RN ×m , see [13]. Let A(ω) = + b b C00 (ω) C01 (ω), finally the loss function to train the network can be written as b max tr A(ω) . ω
Let ω ∗ be the optimal hyperparameters of the basis ψ(x, ω) after training the hidden layers to construct the covariance and cross-covariance matrices, and let Wo ∈ RN ×n be the matrix of the output layer weights for n outputs. To obtain the top n eigenvectors of b ∗ ) for the approximation of n dominant eigenfunctions of the operator, we can solve A(ω the following optimization problem and find a closed-form solution for the output layer b01 (ω ∗ ) Wo , max tr Wo⊤ C Wo ∈RN ×n
b00 (ω ∗ )Wo = I. Wo⊤ C
s.t.
Using the Lagrange multiplier method and the fact that the Koopman operator is selfadjoint, we get the output layer as the solution of the following eigenvalue problem b + (ω ∗ ) C b01 (ω ∗ ) Wo = Wo Λ, C 00 see [13] for more details on the solution of the above optimization problem.
Non-self-adjoint Koopman operators For non-reversible systems, i.e., those for which the Koopman operator is in general not self-adjoint, RaNNDy approximates the singular values and singular functions of the operator using the variational principle for the singular values [20]. Following the above procedure and [13, 20], the projected forward–backward operator (a composition of the b b00 (ω)+ C b01 (ω)C b11 (ω)+ C b10 (ω) and Koopman operator and its adjoint) is given by A(ω) =C we obtain the following loss function to tune the activation function of the network and to find the best hyperparameters ω b max tr A(ω) . ω
Let the columns of Wo′ ∈ RN ×n and Wo ∈ RN ×n denote the top n left and right singular vectors, respectively, that approximate the left and right singular functions of the operator. Then, with the optimal parameters ω ∗ of the activation function, RaNNDy solves the following optimization problem to approximate the singular functions corresponding to the largest n singular values max
Wo ,Wo′ ∈RN ×n
s.t.
tr(Wo⊤ C01 (ω ∗ )Wo′ ), Wo⊤ C00 (ω ∗ )Wo = I, Wo′⊤ C11 (ω ∗ )Wo′ = I,
see [13] for a detailed derivation. Again, using the Lagrange multiplier method, we get the following eigenvalue problem to obtain Wo b00 (ω ∗ )+ C b01 (ω ∗ )C b11 (ω ∗ )+ C b10 (ω ∗ ) Wo = Wo Λ2 . C Using Wo and Lagrange optimality conditions, we can find the weights Wo′ [13]. The procedure of tuning the activation functions of the randomized neural networks is summarized in Algorithm 1.
6
Algorithm 1 RaNNDy optimization Initialization: • Given the training data X, Y ∈ Rd×m .
• For k = 0, parameters ω (0) , learning rate η, for tuning σ (activation function) of the network. • Sample the initial weights and biases of hidden layers from a given distribution with parameters ωW and ωb to get an initial random basis ψ(x, ω). while not converged do Construct Ψ0 (ω (k) ), Ψ1 (ω (k) ). b00 (ω (k) ), C b01 (ω (k) ), C b10 (ω (k) ), C b11 (ω (k) ). Compute C if K is self-adjoint then b00 (ω (k) )+ C b01 (ω (k) ) A(ω (k) ) = C else b00 (ω (k) )+ C b01 (ω (k) )C b11 (ω (k) )+ C b10 (ω (k) ) A(ω (k) ) = C end if Define L(ω (k) ) = tr(A(ω (k) )). Compute gradient ∇ω L(ω (k) ). Update ω (k+1) = ω (k) + η∇ω L(ω (k) ). Set k ← k + 1. end while ▶ With optimized ω ∗ , solve eigenproblem: if K is self-adjoint then b00 (ω ∗ )+ C b01 (ω ∗ )Wo = Wo Λ C else b00 (ω ∗ )+ C b01 (ω ∗ )C b11 (ω ∗ )+ C b10 (ω ∗ )Wo = Wo Λ2 C end if ▶ Select the top n eigenvectors from the sorted eigenvectors in Wo (according to the eigenvalues) to get the approximation of the n dominant eigenfunctions/singular functions.
4 Numerical Experiments We will now show how the proposed approach can be applied to different types of timeseries data.
4.1 Graphons Although we only introduced transfer operators for stochastic differential equations, it has recently been shown in [21] that these definitions can be easily extended to random walks on graphons. Graphons, which can be interpreted as limits of convergent sequences of graphs, are described by measurable functions g : [0, 1]2 → [0, 1], where g(x, y) represents the weight of an edge between the vertices x, y ∈ [0, 1] (or its probability). A graphon is said to be symmetric or undirected if g(x, y) = g(y, x). Otherwise, it is called unsymmetric or directed. Spectral decompositions of transfer operators associated with symmetric graphons are given by X X Pρ = λi ⟨φ bi , ρ⟩ 1 φ bi and Kf = λi ⟨φi , f ⟩π φi , i
π
i
7
(a)
(b)
(c) 1.0
0.20 0.15 0.10 0.05 0.00 1.0 0.0
0.5 0.5
x
y
1.0 0.0
2
x(i)
4
0 1.0 0.0
0.5 0.5
x
0.5
0.0 0
y
1000
2000
i
1.0 0.0
Figure 1: (a) Symmetric graphon g. (b) The associated transition density function p. (c) A random walk on the graphon. where φi is an eigenfunction of the Koopman operator and φ bi = πφi the corresponding eigenfunction of the Perron–Frobenius operator. This also allows us to reconstruct the transition probability density p(x, y) and the graphon g(x, y), i.e., X X λi φi (x)φ bi (y) and g(x, y) = Z λi φ bi (x)φ bi (y), p(x, y) = i
i
R1 where π(x) = d(x) Z , with Z = 0 d(x)dx, is an invariant density and d(x) is the degree function, see [21] for details. Hence, provided that there is a spectral gap after the ith eigenvalue such that λi+1 ≈ 0 (this is the case if there are i clearly separated clusters), we can approximate the graphon g and the corresponding transition density function p using the dominant eigenfunctions. As a basic example, we consider the symmetric graphon defined in [21], given by g(x, y) = 0.2e−
(x−0.2)2 +(y−0.2)2 0.02
+ 0.1e−
(x−0.5)2 +(y−0.5)2 0.02
+ 0.2e−
(x−0.8)4 +(y−0.8)4 0.0005
,
and visualized in Figure 1 along with the corresponding transition density function and a random walk sampled from the graphon that provides the data for the numerical approximations. To demonstrate the application of the proposed algorithm, we initialize RaNNDy with a parametric tanh activation function. We choose the randomized basis ψ(x, ω) = tanh(W (ωW )x + b(ωb )), where ωW and ωb are the scales of normal distribution and the parameters to be optimized with the proposed algorithm. Using [256, 512, 256] neurons in the hidden layers and initializing ωW = ωb = 0.1, we use Algorithm 1 to numerically approximate the top five eigenvalues and eigenfunctions of the associated transfer operators. In Figure 2a, we can see that the initial basis with ω = 0.001 was not suitable and further optimization was required. The algorithm converges within few (less than 20) epochs. As expected, we see a spectral gap after the third eigenvalue in Figure 2b, indicating the presence of three “metastable sets” in the system, which can also be seen from the eigenfunctions of the Koopman operator in Figure 2c. Using these three dominant eigenfunctions of the Koopman and the Perron–Frobenius operators, we reconstruct the graphon and the corresponding transition density, shown in Figures 2e and 2f. Comparison: In Figure 3a, we show a grid-search method to select the most suitable distributions for the weights and biases for the initialization of RaNNDy, illustrating that
8
(a)
(b)
(c)
1.0
0.5
2.25
λi
2
ωW
tr(A)
2.50
0.5
1
2.00
0.0
0 2 4 6 8 10 12 14 16 18
0.0
epochs
1
2
3
4
5
0.0
i
(d)
0.5
1.0
x
(e)
(f )
0.5
−0.5 0.0
0.5
1.0
x
0.4
0.6 0.4 0.2 0.0 1.0
0.0
0.0
0.5 0.5
x
0.2 0.0 1.0 0.0
y
0.5 0.5
x
1.0 0.0
y
1.0 0.0
Figure 2: (a) Training loss for the optimization of RaNNDy with the tanh activation parameter ω. (b) Five dominant eigenvalues of the Koopman operator after training. (c) Three dominant associated Koopman operator eigenfunctions, where denotes the first, the second, and the third eigenfunction. (d) The associated Perron–Frobenius operator eigenfunctions. (e) Reconstructed graphon g with rank 3. (f) Corresponding transition density p.
(b)
zeros
1.5
ones
2.0 1.5
rm
1.0 0.7 0.3 0.1 0.08 0.001
0. 0 0. 01 0 0. 2 0. 1 0. 3 0. 5 1. 8 2. 0 3. 0 5. 0 6. 0 8. 0 0
ni fo tu
ro glo
he
lec
un
no
no
rm
rm
un al ifo r un m ifo rm
1.0
Weights Distribution
2.5
2.5
ωW
1.0 0.5
VAMP Score
2.0
VAMP Score
uniform
ωb
2.5
(c)
VAMP Score
normal
al
Bias Distribution
(a)
2.0 1.5 RaNNDy VAMPnet
1.0 0
1
2
3
4
5
6
7
8
Time (s)
Figure 3: Selecting the distributions for the hidden layers’ weights and biases in RaNNDy. (a) Grid search for different distributions. (b) Loss surface for different values of scales of the normal distributions. (c) Comparison of training the whole network (VAMPnets with learning rate lr = 10−2 ) vs. only the scales of the distributions for the initializations of RaNNDy. We can see that the optimizer is stuck in a local minimum, while RaNNDy converges even for a larger learning rate of lr = 0.5.
9
the effect of the bias distribution is negligible in this case. In Figure 3b, we show the loss surface of the RaNNDy approximation for different values of the scales of the normal distributions for sampling the hidden layers’ weights and biases. Figure 3c shows the comparison of training the full neural network (VAMPnet) vs. tuning the distribution for RaNNDy, showing that RaNNDy converges to a good solution, while the VAMPnet gets stuck in a local minimum. We observed similar behavior for different learning rates and schedulers. For learning rates larger than 10−2 , the VAMPnet approach diverges.
4.2 Bickley Jet We will now show how we can tune RaNNDy to efficiently detect coherent sets [16, 17] and apply the proposed algorithm to the Bickley jet, a simple benchmark problem representing an idealized stratospheric flow [22]. The system is given by ∂Φ − ∂y ẋ = ∂Φ , ẏ ∂x with the stream function y y Φ(x, y, t) = c3 y − U0 L tanh + A3 U0 L sech2 cos(k1 x) L y L + A2 U0 L sech2 cos(k2 x − σ2 t) L y + A1 U0 L sech2 cos(k1 x − σ1 t). L We generate m = 5 000 uniformly sampled initial conditions within the domain [0, 20] × [−4, 4]. Then, using the methods and parameters from the Python library deeptime [23], we simulate the flow from t0 = 0 to t1 = 40. A few snapshots of the flow at different times t are visualized in Figure 4. The particles in red stay close to each other, forming a coherent set, whereas the yellow particles are dispersed by the flow. We again select the tanh activation function with parameter ω. The randomized basis is defined by ψ(x, ω) = tanh(W (ωW )x + b(ωb )). Initializing ωW = 0.001 and ωb = 0.001 as the scales of the normal distribution from which the weights are sampled, we apply the proposed algorithm to optimize the basis and approximate the dominant nine right singular functions of the operator T τ . The results are presented in Figure 5. The algorithm converges in fewer than 10 iterations. The clustering of optimized singular functions into 9 clusters reveals the coherent sets in the Bickley jet.
4.3 Protein NuG2 As a last example, we consider high-dimensional molecular dynamics data, where we aim to understand the folding–unfolding dynamics of a protein, see [24] for more details. We use simulation data of the protein NuG2 acquired from D.E. Shaw Research [25]. Protein G (NuG2) is a 56 residue molecule. The trajectory data has a length of 3.680002 · 108 ps. We subsample the trajectory of 929251 frames to create the training data {xi , yi }m i=1 . With a lag time of 50 frames between each xi and yi , we compute the contact maps (distances between the residue pairs of the molecule) and store them in the data matrices X, Y ∈ R1431×18 585 , here we ignore the near-neighbor residues and self-contact distances.
10
(a)
(b)
(c)
1
1
1 −1
−1
−1 −3
y
3
y
3
y
3
0
3
6
9
−3
12 15 18
−3 0
3
6
9
x
12 15 18
0
3
6
x
9
12 15 18
x
Figure 4: Bickley jet flow at time (a) t = 0, (b) t = 25, and (c) t = 50.
(b)
(a)
1.0
5
λi
2
8
0.2
2 2
4
6
0.6 0.4
1 0
Initial eigenvalues Optimized eigenvalues
0.8
ωW
tr(A)
11
0.0
8
2
epochs
4
6
8
i (d)
(c)
1
1
y
3
y
3
−1 −3
−1 0
3
6
9
−3
12 15 18
x
0
3
6
9
12 15 18
x
Figure 5: Optimization of the activation function for the Bickley jet. (a) The loss function and the activation parameter values. (b) Eigenvalues before and after optimization. (c) Clustering of the initial nine dominant singular functions into nine clusters. (d) Optimized clustering.
11
To approximate the eigenvalues and eigenfunctions of the Koopman operator, we optimize the randomized basis of RaNNDy using the proposed algorithm. With the optimized basis, we approximate the dominant eigenvalues of the Koopman operator, see Figure 6a. The second dominant eigenfunction in Figure 6b represents the two states of the molecule. We can see that the initial randomized basis leads to poor approximation results. In Figures 6c and 6d, we show a few folded and unfolded states along with the contact frequencies in Figures 6e and 6f of the NuG2 protein molecule identified using optimized RaNNDy. The contact map shows the stability of the identified folded and unfolded states.
5 Conclusion and future work We proposed an algorithm to optimize RaNNDy, a randomized neural network framework for transfer operator approximation. RaNNDy allows for a closed-form solution for the output layer of the network representing the eigenfunctions of the transfer operator, making it computationally cheap while mitigating common issues associated with deep learning, such as slow convergence. We further improved the benefits of RaNNDy by making it optimizable, i.e., having parametric activation functions that can be tuned. Instead of training all the weights and biases of the hidden layers of the network, we introduce a few hyperparameters in the activation function and tune them to find suitable weights. This can be regarded as a compromise between fully trained NNs such as VAMPnets and RaNNDy, which is fully randomized. We illustrated the efficiency and accuracy of the proposed approach with the aid of several numerical examples, showing that the optimization of the randomized basis functions is important. Furthermore, tuning the distributions of the hidden weights and biases is easier due to a simpler loss landscape, which could allow us to use algorithms such as Bayesian optimization. It would now be interesting to consider different activation functions with more hyperparameters and analyze the effects of tuning them. Furthermore, a systematic run time comparison for more complex dynamical systems would be helpful to assess the advantages and disadvantages of randomized neural network approaches.
Acknowledgments M.T. was supported by the EPSRC Centre for Doctoral Training in Mathematical Modeling, Analysis and Computation (MAC-MIGS) funded by the UK Engineering and Physical Sciences Research Council (grant EP/S023291/1), Heriot–Watt University and the University of Edinburgh. We would like to thank D.E. Shaw Research for providing the protein molecule data.
References [1] S. Klus and N. D. Conrad. Dynamical systems and complex networks: A Koopman operator perspective. Journal of Physics: Complexity, 5(4):041001, 2024. doi:10. 1088/2632-072X/ad9e60. [2] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley. A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition. Journal of Nonlinear Science, 25:1307–1346, 2015. doi:10.1007/s00332-015-9258-5.
12
(a)
1.0
Initial eigenvalues Optimized eigenvalues
0.8
λi
(b)
0.8
0.6
0.4
0.4
0.0
0.2
−0.4
1 2 3 4 5 6 7 8 9 10
i
0
5000
t
(c)
(d)
(e)
0
10
20
30
(f )
40
Residue
10
50
0
20
0.6
30
0.4
40
0
10
20
30
40
50 0.8
10
0.8
Residue
0
10000 15000
0.6
20
0.4
30 40
0.2
0.2 50
50
Residue
0.0
Residue
0.0
Figure 6: (a) Ten dominant eigenvalues of the approximated Koopman operator. (b) The second dominant eigenfunction to distinguish the folded and unfolded states, where denotes the initial and denotes the optimized eigenfunction evaluated at the data points. (c) & (d) Some folded and unfolded states of the NuG2 protein extracted using the second eigenfunction. (e) & (f) Contact map frequencies between different residue pairs over all the identified folded and unfolded states, respectively.
13
[3] S. Klus, P. Koltai, and C. Schütte. On the numerical approximation of the Perron– Frobenius and Koopman operator. Journal of Computational Dynamics, 3(1):51–79, 2016. doi:10.3934/jcd.2016003. [4] F. Noé and F. Nüske. A variational approach to modeling slow processes in stochastic dynamical systems. Multiscale Modeling & Simulation, 11(2):635–655, 2013. doi: 10.1137/110858616. [5] F. Nüske, B. G. Keller, G. Pérez-Hernández, A. S. J. S. Mey, and F. Noé. Variational approach to molecular kinetics. Journal of chemical theory and computation, 10(4):1739–1752, 2014. doi:10.1021/ct4009156. [6] Q. Li, F. Dietrich, E. M. Bollt, and I. G. Kevrekidis. Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the Koopman operator. Chaos: An Interdisciplinary Journal of Nonlinear Science, 27(10), 2017. doi:10.1063/1.4993854. [7] E. Yeung, S. Kundu, and N. Hodas. Learning deep neural network representations for Koopman operators of nonlinear dynamical systems. In 2019 American Control Conference (ACC), pages 4832–4839, 2019. doi:10.23919/ACC.2019.8815339. [8] A. Mardt, L. Pasquali, H. Wu, and F. Noé. VAMPnets for deep learning of molecular kinetics. Nature communications, 9(1):5, 2018. doi:10.1038/s41467-017-02388-1. [9] M. Gori, A. Tesi, et al. On the problem of local minima in backpropagation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(1):76–86, 1992. doi: 10.1109/34.107014. [10] H. A. B. Te Braake and G. Van Straten. Random activation weight neural net (RAWN) for fast non-iterative training. Engineering Applications of Artificial Intelligence, 8(1):71–80, 1995. doi:10.1016/0952-1976(94)00056-S. [11] L. Zhang and P. N. Suganthan. A survey of randomized algorithms for training neural networks. Information Sciences, 364:146–155, 2016. doi:10.1016/j.ins.2016.01. 039. [12] A.K. Malik, R. Gao, M. A. Ganaie, M. Tanveer, and P.N. Suganthan. Random vector functional link network: recent developments, applications, and future directions. Applied Soft Computing, 143:110377, 2023. doi:10.1016/j.asoc.2023.110377. [13] M. Tabish, B. Leimkuhler, and S. Klus. How deep is your network? Deep vs. shallow learning of transfer operators. arXiv preprint arXiv:2509.19930, 2025. [14] I. Mezić. Spectral properties of dynamical systems, model reduction and decompositions. Nonlinear Dynamics, 41:309–325, 2005. doi:10.1007/s11071-005-2824-x. [15] C. Schütte and M. Sarich. Metastability and Markov state models in molecular dynamics, volume 24. American Mathematical Soc., 2013. URL: https://bookstore. ams.org/cln-24. [16] G. Froyland. An analytic framework for identifying finite-time coherent sets in timedependent dynamical systems. Physica D: Nonlinear Phenomena, 250:1–19, 2013. doi:10.1016/j.physd.2013.01.013.
14
[17] R. Banisch and P. Koltai. Understanding the geometry of transport: Diffusion maps for Lagrangian trajectory data unravel coherent sets. Chaos: An Interdisciplinary Journal of Nonlinear Science, 27(3), 2017. doi:10.1063/1.4971788. [18] P. Koltai, H. Wu, F. Noé, and C. Schütte. Optimal data-driven estimation of generalized Markov state models for non-equilibrium dynamics. Computation, 6(1):22, 2018. doi:10.3390/computation6010022. [19] P. I. Frazier. A tutorial on bayesian optimization, 2018. URL: https://arxiv.org/ abs/1807.02811, arXiv:1807.02811. [20] H. Wu and F. Noé. Variational approach for learning Markov processes from time series data. Journal of Nonlinear Science, 30(1):23–66, 2020. doi:10.1007/ s00332-019-09567-y. [21] S. Klus and J. J. Bramburger. Learning graphons from data: Random walks, transfer operators, and spectral clustering. arXiv preprint arXiv:2507.18147, 2025. doi: 10.48550/arXiv.2507.18147. [22] I. I. Rypina, M. G. Brown, F. J. Beron-Vera, H. Koçak, M. J. Olascoaga, and I. A. Udovydchenkov. On the Lagrangian dynamics of atmospheric zonal jets and the permeability of the stratospheric polar vortex. Journal of the Atmospheric Sciences, 64(10):3595–3610, 2007. doi:10.1175/JAS4036.1. [23] M. Hoffmann, M. Scherer, T. Hempel, A. Mardt, B. de Silva, B. E. Husic, S. Klus, H. Wu, N. Kutz, S.L. Brunton, et al. Deeptime: a Python library for machine learning dynamical models from time series data. Machine Learning: Science and Technology, 3(1):015009, 2021. doi:10.1088/2632-2153/ac3de0. [24] C. Schütte, S. Klus, and C. Hartmann. Overcoming the timescale barrier in molecular dynamics: Transfer operators, variational principles and machine learning. Acta Numerica, 32:517–673, 2023. doi:10.1017/S0962492923000016. [25] K. Lindorff-Larsen, S. Piana, R. O. Dror, and D. E. Shaw. How fast-folding proteins fold. Science, 334(6055):517–520, 2011. doi:10.1126/science.1208351.
15