Self-Supervised Evolution Operator Learning for High-Dimensional Dynamical Systems
Authors: Giacomo Turri
Luigi Bonati
Kai Zhu
Massimiliano Pontil
Pietro Novelli
Venue: ICLR 2026 Poster (2026)
Keywords: Operator
Koopman
Transfer
Contrastive
Self-Supervised
DMD
Modes
Dynamics
ENSO
Climate
Molecular Dynamics
Protein
TICA
Slow Modes
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: We introduce an end-to-end approach to learn the evolution operators of large-scale non-linear dynamical systems, such as those describing complex natural phenomena. Evolution operators are particularly well-suited for analyzing systems that exhibit spatio-temporal patterns and have become a key analytical tool across various scientific communities. As terabyte-scale weather datasets and simulation tools capable of running millions of molecular dynamics steps per day are becoming commodities, our approach provides an effective tool to make sense of them from a data-driven perspective. The core of it lies in a remarkable connection between self-supervised representation learning methods and the recently established learning theory of evolution operators. We deploy our approach across multiple scientific domains: explaining the folding dynamics of small proteins, the binding process of drug-like molecules in host sites, and autonomously finding patterns in climate data. Our code is available open-source at: https://github.com/pietronvll/encoderops.

Published as a conference paper at ICLR 2026
SELF-SUPERVISED EVOLUTION OPERATOR LEARNING FOR HIGH-DIMENSIONAL DYNAMICAL SYSTEMS
Giacomo Turri!
Luigi Bonati
Kai Zhu
Massimiliano Pontil 1.4
Pietro Novelli
1 "Computational Statistics and Machine Learning, Italian Institute of Technology, Genoa, Italy
2 AAtomistic Simulations, Italian Institute of Technology, Genoa, Italy
Colleee of Pharmaceutical Sciences, Zhejiang University, Hangzhou, China "AI Centre, University College London, London, United Kingdom
ABSTRACT
We introduce an end-to-end approach to learn the evolution operators of large- scale non-linear dynamical systems, such as those describing complex natural phenomena. Evolution operators are particularly well-suited for analyzing systems that exhibit spatio-temporal patterns and have become a key analytical tool across various scientific communities. As terabyte-scale weather datasets and simulation tools capable running millions of molecular dynamics steps per day are becoming commodities, our approach provides an effective tool to make sense of them from a data-driven perspective. The core of it lies in a remarkable connection between self- supervised representation learning methods and the recently established learning theory of evolution operators. We deploy our approach across multiple scientific domains: explaining the folding dynamics of small proteins, the binding process of drug-like molecules in host sites, and autonomously finding patterns in climate data. Our code is available open-source at: https: https://githuss com/ rompll iencoderops
INTRODUCTION
Dynamical systems are fundamental to understanding phenomena across vast range scientific disciplines. from physics and biology to climate science and engineering. Traditionally, scientists have modeled these systems by formulating differential equations from first principles. However, as systems grow in scale and complexity, this approach quickly becomes computationally burden- some and difficult to interpret (Anderson) 1972). hindering the study of large-scale phenomena. Simultaneously, advancements in data collection techniques and computational power have led to
an explosion of available data from experiments (Hersbach et al. 2020 Chanussot et al. 202 and
high-fidelite simulations (Harvey al. 2009; Abraham et al. 2015 Eastman et al. [201]] Bauer
ct al. 2015). This abundance of data makes data-driven approaches increasingly appealing for
studying complex dynamics, with machine learning (Shalev-Shwartz. & Ben-David (014) becoming a dominant paradigm for learning dynamical systems, largely focusing on predictive tasks such as forecasting. The recent revolution in data-driven weather modeling (Kurth et al., [023] Bi al.)
2022, Lam et al. [2023] Kochkov et al. 2024) stands as a paradigmatic example of ML's power in handling complex spatio-temporal dynamics. Similarly, reinforcement learning (Sution Barto 1998) has reimagined control theory by leveraging data-drive strategies to optimize system behavior. While these data-driven methods excel at prediction and simulation, there remains a significant gap in approaches that offer interpretability. In scientific contexts, merely predicting system behavior is often insufficient; understanding why system evolves a certain way paramount. For instance, comprehending the dynamical shorteuts and bottlenecks happening through atomistic interaction is crucial for understanding why a drug binds to specific target or fails to do so, a level of insight not typically provided by black-box predictive models.
modeling paradigm particularly well-suited for interpretability is that of evolution operators Lasota & Mackey 1994 Applebaum| (009). Under mild assumptions, dynamical systems and stochastic processes can be represented by a linear operator a mathematical entity that maps functions to other functions. This operator-based approach offers multiple advantages. First, linearizes the dynamics,

---

Published as a conference paper at ICLR 2026
greatly simplifying tasks like forecasting and controller design. Second, these operators possess spectral decomposition (Reed & Simon 1972). which expresses the system' complex dynamics as linear combination of fundamental, coherent apatio-temports modes (Molgedey & Schuster, 1994). Each mode represents a distinct, intrinsic pattern associated with a unique spatio-temporal structure defined in terms of growth or decay rates and oscillation frequencies. By identifying and analyzing these principal modes, researchers gain deep insights into the underlying mechanisms driving the system's s macroscopic behavior, offering a structured, physically meaningful understanding.
Building on the understanding that evolution operators provide a powerful framework for inter- pretable analysis, significant effort has been directed towards learning these operators directly from data Kovachki et al. (2023). Data-driven approaches for this task emerged already in the early 2000s, including pioneering work utilizing transfer operators for analyzing stochastic processes in computational biophysics (Schütte et al. [000], as well as the dynamic mode decomposition family of methods (Schmid) 2010) for deterministic systems via the Koopman operator. In the ensuing years, there has been a significant acceleration in machine learning methods for evolution operator learning, encompassing theoretical advances through kernel methods and powerful end-to-end deep learning approaches.
Contributions. In this work, we build upon these recent foundations, showing how evolution operator learning can be scaled to structured and -ighndimensional dynamical systems. We formalize a principled end-to-end protocol that is amenable to GPU training and prove its equivalence to a self-supervised representation learning problem. Leveraging this link, also show the transferability of our trained models in both molecular dynamics and climate settings. Code, data, and weights are made available open-source.
2 EVOLUTION OPERATORS AND HOW TO LEARN THEM
Evolution operator learning is data-driven approach to characterizing dynamical systems, either stochastic, Zt+1 ~ p(-[2z), , or deterministic, Zt+1 ~ 5(- - F(x:)). Throughout, we assume the æ alone and not on the states at
dynamics to be Markovian, times s < t. If this assumption is not satisfied by 3t, a standard trick to re-define the state as a context c f (æt, F2-1, ................................ It-H) with history length H, where f can be simple concatenation, or a learned sequence model (e.g., a recurrent neural network or transformer).
that the evolution of 2; depends
Evolution operators defined as follows: for every function of the state of the system, (x.) is the expected value of f one step ahead in the future, given that at time / the system was found in zt
-f f p(dyla)) = F 1 (x))
(E)) =
=
Notice that E is an operator because it maps any function f to another function, 2t + (z2) and is linear because E(f +g)) = Eff aEg. When the dynamics is deterministic, is known as the Koopman operator (Koopman 1931). while in the stochastic case it is known as the transfer
operator (Applehaum) 2009).
Evolution operators fully characterize the dynamical system because knowing E allows us to recon- struet the dynamical law p(-(az). Indeed, for any subset of the state space C X, applying the indicator function of B, we have
(EI - p(dylæ,) = P [X++1 E B|xr].
An advantage of the operator approach over dealing directly with the conditional probability p(- is that E acts linearly on the objects to which is applied. This means that operators unlock an arsenal of tools from linear algebra and functional analysis, which would be unavailable otherwise. Arguably the most important of them is the spectral decomposition, allowing us to decompose E, and hence the dynamics, into linear superposition of dynamical modes. These ideas lie at the core of the celebrated Time-lagged Independent Component Analysis (Molgedey & Schuster 1994 Pérez-Hernnndez et al. 2013). and Dynamical Mode Decomposition (Schmid) 2010) Kutz al. 2016))
A generalization of the eigenvalue decomposition of matrix.
2

---

Published as a conference paper at ICLR 2026
2.1 LEARNING E AND ITS SPECTRAL DECOMPOSITION FROM DATA
We now review the main approaches to learn the evolution operator and its spec- tral decomposition from finite dataset of observations, with emphasis on the least squares approach, which is essential to understand every other method as well. A core idea of operator learning is that operators are defined by
4= ((a)
E
how they act on a suitable lin- ear space of functions, similarly
Latent Space
Linear Evolution
E
E
to how matrices are defined by their action on basis of vectors. Of course, not every function f interesting, and this nicely par- allels with the matrix example,
Non-Linear Representation
most of the variance in the data. Observations
where the most "interesting" di- rections are those that recover
Learning E, therefore, is usually cast the following problem:
Letting (xx) E Rd be learned fixed encoder of the state, find the best approximation of E restricted to the d-dimensional linear space of functions generated by y. given the data.
In practice, the data is usually a collection of transitions D (z,, 3)) where is intended that T. P[X,] are sampled from a distribution of initial states, while yı p(-[ar))
Least squares. In this approach the encoder yp frozen, that is non-learnable, dictionary of functions, and we are interested in approximating the action of E on functions of the form f(z) (w, (xx)) for every w E Rdd To this end, one minimizes the empirical error between the true conditional expectation Eym Ix, [ww, (pv)) |r]] and linear model (Ew, spx)), where the matrix E E Raxd identifies the restriction of the evolution operator the linear span of the dictionary:
N Σ ((w, ((()) - (Ew, ((x,)) N
N ll(())) - E ((x)))1 + A|EE Σ
(2)
On the right-hand side, we assumed |wi|| < 1, used the ridge penalty. The minimizer of (2) can be computed in closed form (Korda & Mezić 2018 Kostic et al. 2022 and references therein) as
inequality, and added
1 E ((x))p())) and Cx Cxx.
Ep (Cx + A Cxx, with CXY
In the limit of infinite data, N + 0o. and infinitely dimensional encoders, d - oo, the least squares estimator converges (Korda & Mezić 2018) in the strong operator topology to the evolution operator E, and similar (but weaker) asymptotic convergence results are proved for its spectrum.
Mode decomposition. The spectral decomposition of E is approximated by expressing the least- squares estimator in its eigenvectors' basis E.p- = QAQ where the columns of = 191. qa] are the eigenvectors of Ep, and A is a diagonal matrix of eigenvalues. In this basis, the expected value in the future for function (xx (w, (ox)) is expressed
d
(xxXmixx,((x)) Re (E,w, p(x)) = (QAQ w. p(x)) = Σ A (g., x(x)) (Q ww).
(4)
The spectral decomposition expresses the transition Iy + F++1 as a sum of modes of the form A (g., (fx)) (Q-'w), each of which can be broken down into three components:
1. The eigenvalues X; determine the time scales of the transition. Indeed, applying the evolution operator times analyze the transition 2; + Fi+ leaves (4) unchanged, except that each A; becomes X₁. Writing X; = pfe in polar coordinates, reveals that the modes decay exponentially over time with rate Ps, while oscillating at frequency wi.
3

---

Published as a conference paper at ICLR 2026
2. The initial state influences the decomposition through the factor (x(x) = (gt, p(x)). This coefficient captures how strongly the state æ aligns with the i-th mode. When 9; corresponds to an cigenvalue with slow decay, i... |r,] = 1, the term (z) serves as natural quantity for clustering states into coherent or metastable sets.
3. The coefficient ww) in turn, indicates how the function represented by the vector w relates the i-th mode. This connection makes possible link the dynamical patterns to specific functions - or observables - thereby deepening our understanding of the system.
Kernel methods. Leveraging the kernel trick, one can learn evolution operators by deriving a closed-form solution of (2) in terms of kernel matrices whose elements are of the form xz) = (xI)) with k(;;)) a suitable kernel function. Thanks to the theory of reproducing kernel Hilbert spaces, this class of methods is backed up by statistical learning guarantees, such as the ones derived in (Kostic et al. 2022 2023 Nüske et al. (022). Similarly to the least-squares approach, one also approximates the spectral decomposition of via kernel methods, and this task captured quite a lot of attention from researchers in this area, see (Williams et al. 2015; Kawahara 2016 Klus et al.
2020
Das & Giannakis 2020: Alexander & Giannakis 2020; Meanti et al. 2023
Deep learning. In contrast to the previous approaches, where the encoder yp is prescribed, a number of methods proposed to approximate from data with end-to-end schemes including & as a learnable neural network. Since learning E ultimately entails learning its action on the linear space spanned
by şo, it is appealing to choose an encoder capturing the most salient features of the dynamics. To
this end, one can train oo via encoder-decoder scheme proposed in (Takeishi et al. 2017 Lusch
eet al. 2018) Otto & Rowley 2019; Azencot et al.) 20200] Wehmeyer & Noē 2018 [Frion et al. 2024)
or with encoder-only approaches in (Li et al.) |
MMrdt et
al. 2018 Yeung al. 2019 Kostic
iet al. 2024b Federici et al. 2024 Jeong et al.) 2025).
In encoder-decoder schemes, şp is trained alongside a decoder network, minimizing combination of prediction and reconstruction error Yet, while minimizing reconstruction loss biases the model towards accurate forecasts of the near future, recent work on model-based reinforcement learning, where one is instead interested in long-term behaviours, suggests that the presence of decoder detrimental (Lyu al., Similarly, (Balestriero & LeCun 2024) showed how features learned by reconstruction are both
2023 Schwarzer al. 2020;; Hansen et al., (002)) to control tasks,
uninformative for perception and hardly transferable.
On the other hand, the competitive advantage of evolution operators over techniques such as (uurth jet al., 2023 Lam et al. 2023, Piaff et al., 2021 (Sanchez-Conzalez et al., 2020; Li et al. 2020) lies in their spectral decomposition, useful for interpretabilityy reduced order modeling, and control tasks. Encoder-only approaches follow this intuition and prioritize approximating the spectral decomposition of E over the raw forecasting performances. Concretely, this is accomplished via loss functions that are minimized when spans the leading singular space of E. Clearly, once an encoder has been trained, it can be transferred to similar dynamical systems, as we demonstrate in 4.2 and 4.3
In this work, we propose an encoder-only method based on a loss function originally designed for self-supervised representation learning. Our approach is numerically stable, scales efficiently, enables transfer across related systems, and can incorporate structural priors, e.g., graph-based encoders, beyond the reach of classical DMD approaches. Though our approach is broadly applicable, we mainly focus on applications involving interpretability and model reduction of scientific dynamical systems, highlighting how ML evolution operators can help in advancing fundamental science. Recent works in RL (Lyu et al., 2023, [Schwarzer et al., 2020; Rozwood et al. 2023 Novelli et al.) 2024 suggest that our approach can be relevant for control tasks, but we leave this for future work.
3 LEARNING EVOLUTION OPERATORS VIA SELF-SUPERVISION
As discussed above, we are interested in the evolution operator
p(yyar)) p(y)
(E)))x) = I [(w))aa.] = Ey
[(v)
Notice that trying to minimize the prediction error 2 alone immediately leads to representation collapse with to mapping every input 0 to obtain prediction error of O.
4

---

Published as a conference paper at ICLR 2026
where in the last equality expressed the expectation in the form of an importance sampling estimator with respect to the probability of the future state P(X(+1). In so doing, we link the evolution operator E to the density ratio
r(r,, x(+1) p(14+1)(2))
(6)
p(z(+1)
The core of our operator-learning scheme, Alg. D is to optimize a model for the density ratio (6) parametrized as the bilinear form (f((x)) ( Pp(x=++))) )), similar to van den Oord et al. (2019). Here, +p
is a d-dimensional encoder, while is a linear predictor layer which, as discussed below, equals the
action of E on the linear subspace of functions spanned by i,, up a known linear transformation.
We minimize the L error between the density ratio and our bilinear model P((e)+))((
E P) = xxxx)))xxxxx.
(r(2, - (1(xx),
-[(x,))x,xxx:
(v())) (xx ((p(()))PP(u)) cst., (7)
where Ex,0X is the expected value between the product of the marginals X₂ and X Estimating the squared term in ( via U-statistics (Hoeffding) 1992) and foregoing the constant term, we finally get to the empirical loss
N
Σ (v(x;), E(pp P)= N(N 1) ンØ1
Σ N (o(a;))
2
(8)
The loss function ( was originally proposed for self-supervised contrastive learning in HaoChen
et al. 2021 2022). Indeed, noticing that (f((r). Pcp(y)) can be interpreted as a measure of similarity between J and y. the first term of (8) minimizes the similarity between randomly chosen i # (i.e.,
negative) pairs, while the second term maximizes the similarity of consecutive (i.c., positive) pairs.
The loss (8) has also been applied in reinforcement learning (Ren al. 2023). causal estimation (Sun)
et al. 2025), and recently Lu et al. (2024) showed that it belongs to wide class of contrastive
learning losses defined by Csiszár d-veeeggesess. Concurrently Wang et al. (2022): ;
Ryu al. (2024) Kostic et al. (2024a); Jeong et al. 2025) applied it to approximate the SVD of linear operators.
3.1 THEORETICAL PROPERTIES OF OUR APPROACH.
Our model (p((x)), 1)( )) is characterized by the presence a linear predictor P. and by shared encoder between 2, and 21+1, in contrast the more agnostic choice ((p(x.). ), th (2++11) adopted by HaoChen et al. (2021): | Wang et al. (2022) : Ryu et al. (2024); Kostic et al. (2024a] Jeong now prove number of theoretical results highlighting
ett al. (025). Our choice deliberate, and
how our parametrization is particularly apt for evolution operator learning, with the predictor layer P having key role. Every lemma in this section is proved in Appendix [ The first observation was already noticed in Wang et al. (2022). and provides a direct link between 7 and the evolution operator regression formalism developed in Kostic al. (2022).
Let sp X + Rd be an encoder whose components are squure--ntterubbe with respect to both μ and v, and let E be Kilbert-Schmile evolution operator. Then, the loss function 7 is equivalent to the following operator learning loss:
Lemma 1.
=||E-
&
We highlight that when the operator E is not Hilbert-Schmidt, the loss function (8) can still be linked to operator learning. In this more general scenario, the loss promotes encoders oo displaying both strong dynamical response Exp, and a good approximation of the true dynamics, see Appendix A..
Plugging our model back into (5). the evolution operator gets parametrized as (E,) x Ey-xx.,, ((p(x)), Pop(w) /(y)). and if the function f is in the linear span generated by the encoder f(v) (p(w), w)) we can simplify the expression above
(Ef) (
P (Eg-xen (e(w) w) = (p(x)) PCrw))
"That is, the product measure P[X,] & P[Xe+1)
5

---

Published as a conference paper at ICLR 2026
where introduced the covariance of the futures Cy = Ey-.x. (x(())())))) Thus, the linear predictor P parametrizes the approximation the evolution operator E over the diite-dimensionll space generated by the state representation (p. We remark that to be sure that prescribed function lies in the span of the encoder, one can add as a non-trainable component of the architecture done in Appendix B.1 Alternatively,
p(()) = [>) [==),
can compute its least-squares
approximation f on the features spanned by şp, and use that in place of f. The following Lemma shows that when the predictor P is optimal, one recovers the least squares estimator (3
Lemma 2. For any fixed the predictor P minimizing 7 can be computed in closed form P. "x'CxyCy and the model for the evolution operator is given by
E = P.Cr = Cx = Eq. 3 with X - 0,
coinciding with the least-squares estimator
The final interesting fact about 7 is relation to the VAMP score (Wu & Noë 2020), originally introduced for representation learning of molecular kinetics. In particular, the VAMP-2 score can be defined in terms covariances as
VAMP = | Cxy 1/2 llis
(10)
Lemma For any fixed up, let P. the optimal predictor of E(p,, P), as Lemma Then, the following holds true:
c(p, P.) - Cxx Cx Cxx
1/2 = VAMP2(g).
Our loss function, therefore, matches the negative VAMP-2 score when P optimal. Compared to methods that directly maximize the VAMP score, such (Mardt et al.) 2018), however, our approach does not require matrix inversions in the computation of the loss (Wu Noč 2020), an operation which is unwieldy and prone to instabilities in large-scale applications. Instead, the loss function (8) written terms of simple matrix multiplications, making it perfect for GPU-based training.
3.2 PRACTICAL IMPLEMENTATION
for & 1 to nuw_sutepp steps do
B + {(I, yi) D forall i do 2, t p(x;)) and qs t Py(())
xc
end
(z, 9) E(v,
rij f- (zi, 9i)
B(B-1) Enjj 2 6, - (5p, P.di dP)
de, dP + V
B Σ, ra)
end
Algorithm 1: A pair of consecutive observations (æ, y) from dynamical system are mapped to representations and q via an embedding function y. The representation q is also processed by predictor P. The algorithm iteratively optimizes and P using the contrastive objective based on the similarity (z.g).
The implementation of our method, summarized in Alg. T follows standard self-supervised learning procedures (Chen al. 2020; Grill et al.) 2020 Zbontar al. [021 Chen & He (202T). There, positive pair of data- Sata-points, in our case pair of consecutive observations of the dynamical system, processed through an encoder network sp, and optionally a predictor network, which in our case is simple linear layer P. We apply simplicial normalization (Lavoie et al.) (002)) to the outputs of the embedding sp. To keep our implementation as close the theoretical insights as possible, we didn
concatenate additional projection heads to the encoder yo, as suggested in (Chen et al. 2020 Grill
al. 2020 Zbontar et al. 2021 Chen & He, 202T). Furthermore, because of the identity we kept P linear, but it is worth mentioning that tiny MLPs might be employed as predictors instead.
"Backpropagation through inversions may lead gradient explosion (Goluh & Pereyra 1973)

---

Published as a conference paper at ICLR 2026
A
B
Implite Limencale (ndl thigher better)
This work DPNETS
LORA
VAMP
1.0
(U] unfolded
1.0
JointLLRA
-- SespLoRA
RMSD
$
Folded
20'
-0.5
0.5
-1.0
(F) Folded
0.0
12 14
10 10
50
75
100
125
150
175
200
Time (us)
Unfolded
Figure I: Trp-Cage foldding A: Implied timescales derived from different baseline methods. Higher implied timescales associated with more accurate approximations of the alow-modes B: Time series of the leading eigenfunction v (red, left axis) alongside RMSD (gray, right axis), capturing transitions between folded and unfolded states. Representative snapshots of each state are shown. In the folded structure, key hydrogen bonds identified as relevant by the LASSO model are highlighted.
Once representation şo is learned, we model the evolution operator E via the least-squares estimator Eg from (3). To compute it, one can make use of the closed form expression (3) by computing the covariances at the end of the training of (p, done in Kostic et al. (2024b); Jeong al. 2025). This two-step procedure, however, requires a full forward pass over the full training dataset, which is impractical for large problems. Another option is use (9). but this again requires the evaluation of the covariance, and might be suboptimal whenever P isn' yet converged to the true minimizer P.. In implementation, instead, kept buffers for Cx and Cxy, which updated online during the training loop via an exponentially moving average of the batch covariances. At the end of the training, we use buffers compute Eip as in (3.. In Appendix [B.5] we show that covariances updated online during training converge to an accurate approximation of the true covariances, and yield identical (or slightly better) results compared re-evaluating the covariances from scratch.
4 EXPERIMENTS
We now put to the test our method on high-dimensional dynamical systems from both molecular dynamics and climate domains. Our focus is assessing the capability of the method to decompose complex dynamics and to evaluate the generalizability of the learned representations. In Appendix[B.. [B.]] we also report an additional experiment the Lorenz '63 system (Lorenz, 1963). illustrating how the method can also used for small forecasting tasks.
4.1 HIGH-RESOL.UTION DYNAMICAL MODELING OF PROTEIN FOLDING
The Trp-Cage miniprotein widely studied benchmark for protein folding due its small size and fast dynamics (Lindorf'-Land Larsen et al. 2011). Previous works, including SRV-based Markov State Models (Sidky et al. (219)) and Graph VAMPNet (Ghorbani et 2022). have modeled Trp-Cage dynamics using coarse-grained representations, where the state of the system is defined by the small subset of 20 C. atoms the backbone of the protein. Our approach allows
scale to a more expressive molecular representation based on all 144 heavy atoms, employing the SchNet (Schütt al. 2017) graph neural network architecture as the encoder s. After training, we calculate the eigenvalue decomposition of the evolution operator as described in Sec. 2 As shown in Fig. [3, the leading eigenfunction V (x)) = (q,, sp(x)) correlates strongly with the system's deviation (RMSD) from the folded structure, confirming that V encodes the Dolding-to-unfoloe transition. Clustering the molecular configurations according to the values of V reveals a clear separation between folded and unfolded ensembles (see snapshots in Fig. BB)
not-mean-squan
To interpret the nature this slow mode, we regress 1 against a library of physically meaningful descriptors specciiccll.....l.l...... hydrogen bond interactions across residue tatrs--using sparse LASSO model (Brunton et al. 2016)] Zhang et al. 20244) Novelli et al., (022). This analysis reveals a network of hydrogen bonds stabilizing the folded state, including contributions from side-chain interactions
7

---

Published as a conference paper at ICLR 2026
1/
(B) Bound
2/2
Learning from scratch Transfer Learning
[5] Somi bound
(U) Unbound
24
26
Time [ns]
Time [ns]
Figure 2: Calixarene binding. Eigenfunctions V 1 (left) and 222 (right) capture ligand transitions from unbound (U) to semi-bound (S) and bound (B) states. The model using a representation transferred from other ligands (solid line) closely matches one trained from scratch (dashed).
that would be invisible to coarse-grained dynamical models such as (Ghorbani et al.. 2022 Finally, we note that the implied timescale derived from the leading eigenvalues of the learned operators are influenced by both the choice of representation and learning loss. The LoRA baselines from [Jonn iet al. (2025) use similar loss function, but do not share the encoder between I and 2 VAMP, from Mardt et al. (2088) and DPNETS from Kostic et al. (2024b), on the other hand, minimize different loss functions. According to the variational principle for Markov processes (Wu & Noë 2020, Noé & Nüske 2013), higher implied timescales indicate a better approximation of the system's s true slow dynamics (see Fig. [A).
4.2 LEARNING TRANSFERABLE REPRESENTATIONS FOR THE BINDING OF SMALL MOLECULES
Our second case study focuses on the binding of small molecules to a calixarene- based system Yin ett al. 2017), which is often used as simplified model to study the dynamical processes relevant, for instance, in drug design. Our baseline is obtained by using Alg. Ioo train an encoder yp on molecular dynamics data describing the binding dynamics of single molecule (G2) to the host system. As in the previous example, we employ a SchNet architecture for y. As shown in Fig. 2 the slowest dynamical mode, captured by the dominant eigenfunction 1/, is associated with transition between a semi-bound configuration and the fully bound state. Structural inspection reveals that this intermediate state corresponds to misaligned pose of the guest, caused by the presence of water molecule occupying the binding pocket. The second eigenfunction V2 instead resolves the unbound-to-ben -o-bunn transition. Our findings align with previous works (Rizzi et al., 2021), where water occupancy was identified as a key kinetic bottleneck in host-guest interactions.
We now turn key question: can representation trained one set of molecular systems generalize to others? This capability is essential for scalable modeling in applications like drug discovery, where retraining a model for every new compound is prohibitive. To test the transferability of the representations & trained with our method, we trained the encoder on molecular dynamics simulations for two molecules (Gl and G3), and used it to analyze the binding dynamics of a different ligand (G2). Using the frozen encoder, we compute the evolution operator of (G2) via (3). and examine its dominant eigenfunctions. Remarkably, the transferred representation successfully recovers the key dynamical modes of the binding process of (G2) without having seen it during the representation learning phase. In particular, it recovers both the entry of the guest molecule into the host cavity and its final locked configuration Fig. 2 This result illustrates that our self-supervised model learns features that are not only informative but also transferable across molecular systems.
4.3 PATTERNS IN GLOBAL CLIMATE
Finally, we test our method on climate data. Specifically, we aim to retrieve El Niño- -oouhSru Oscillation (ENSO), one of the most influential sources interannual climate variability (Diaz & Markgraf 2000 Callahan & Mankin 2023), arising from coupled ocean-atmosphere dynamics in the tropical Pacific (Bjerknes) 1969 Philander 1983). Characterizing ENSO remains a central goal in climate science, particularly in the context of its potential changes under global warming
'The implied timescale can be computed from the eigenvalues as 7, AAt/log(A)), where At is the time lag between two consccutive observations.
8

---

Published as a conference paper at ICLR 2026
B
ONI
WENSO
90
80
6
30
0
0
30'S
60'S
90%5
60
120'E 180 1
60
1980 1985 1990 1995 2000 2005 2010 2015 2020 Time
Figure 3: ENSO mode retrieved with our method. A: Mode associated with the second leading eigenfunction, highlighting dominant activation in the tropical Pacific. Boxes indicate standard ENSO monitoring zones. B: Right eigenfunction corresponding to the second leading eigenvalue, compared to the ONI index. The vertical line marks the split between training and validation sets.
(ETT).NSSY ENSO conventionally characterized by monthly-averan sea surface temperature (SST) anomalies, denoted as SST*, computed following the procedure described (NCP Center). The SST fields are obtained from the ORAS5 reanalysis (Zuo et al. 2019) and provided through the ChaosBench dataset (Nathaniel et al. 2024). However, the dataset comprises only 540 snapshots, which may hinder effective model training. To overcome this, analogous to the drug design experiment described in Sec. 4.2 we adopt transfer learning strategy using longer synthetic trajectory generated by the Community Earth System Model (CESM) (Hurrell et al. 2013), consisting of 12,598 samples. The objectives of this experiment twofold: (i) to determine whether our method can retrieve ENSO dynamics, and (ii) to assess whether representations & learned from simulated data can be effectively transferred real-world climate observations.
(McPhaden et al., 2006, Cai et al. 202T).
A convolutional neural network-based encoder is trained using the simulated SST* fields, after which the learned representation @o is applied to real data. The transfer operator is then estimated following (3) using the period 1979-2016 for training and 2017-2023 for validation, and subsequently subjected to spectral decomposition to extract the dominant modes. As expected, this procedure recovers modes corresponding to known climate periodicities, such as annual oscillations (see Tab, 2) Remarkably, one of the leading nontrivial modes (second in magnitude) clearly reflects ENSO dynamics. The associated right eigenfunction exhibits strong Pearson correlation (r = 0.82, p .001)
with the Oceanic Niño Index (ONI) (Fig. 3B), a widely used metric for ENSO monitoring (Glantz &
Ramirez 2020), while the associated spatial mode shows dominant activation over the tropical Pacific (Fig. AA) Importantly, our method generalizes effectively to unseen data, successfully detecting the 2023 El Niño event within the validation set. It is worth noting that training the same model directly on observational data also recovers the ENSO mode; however, the correlation between the associated right eigenfunction and ONT is weaker against VAMPNets (Mardt et al. 20188) and DPNets (Kostic et al. 22024b) with results indicating that our approach achieves stronger correlation in capturing the ENSO mode (see Appendix B.4).
0.71, p P .01).. Additionally,
compared our method
This experiment underscores the model's ability to autonomously identify complex climate phenom- ena in an unsupervised manner without prior localization (unlike previous approaches (Froyland et al. 2021 Lapo et al. (025)). Importantly. the transfer learning approach enables the model to leverage knowledge from large, high-quality simulations to mitigate for the scareity of observational data, thereby enabling a more robust extraction of complex patterns such as ENSO. These findings highlight the ability of our approach to learn a robust and generalizable representation, effectively transferring knowledge from synthetic simulations real-world observations
5 CONCLUSION
In this work, we proposed an end-to-end framework for learning evolution operators and their spectral decomposition. Our method scales effectively to large and complex systems, making it practical tool for uncovering physically meaningful patterns in their dynamics. By leveraging connection between contrastive learning objectives and the spectral properties of evolution operators, we break new ground on the transfer of dynamical representations. Our experiments atomistic and climate

---

Published as a conference paper at ICLR 2026
systems demonstrate the versatility of our approach and its generalization capabilities. Looking ahead, this connection opens the door to more expressive learning architectures, robust training strategies, and broader applications in scientific discovery and control.
Limitations. Due to the nature of our experiments, evaluation was more qualitative than typical in ML: benchmarks specifically targeting the accuracy of the spectral decomposition are, to the best of our knowledge, not yet available.
REPRODUCIBILITY STATEMENT
Theory. Our theoretical claims are supported by complete proofs provided Appendix [A] Code. All code used in this study is available at btttps://gith/i com/ httpl//coopppeetteeeeeeeeeeeeeee pietronvi1/ encoderops] Detailed experimental procedures and implementation are described in Appendix B Datasets, Instructions for generating, downloading, requesting the datasets are included in Appendix B
ACKNOWLEDGMENTS
This work was partially funded by the European Union NextGenerationEU initiative and the Italian National Recovery and Resilience Plan (PNRR) from the Ministry of University and Research (MUR), under Project PE0000013 CUP J53C22003010006 "Future Artificial Intelligence Research (FAIR)), and European Project ELIAS N. 101 120237. We acknowledge ISCRA for awarding this project access to the LEONARDO supercomputer, owned by the EuroHPC Joint Undertaking, hosted by CINECA (Italy).
REFERENCES
Mark James Abraham, Teemu Murtola, Roland Schulz, Szilárd Páll, Jeremy C. Smith, Berk Hess, and Erik Lindahl. Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX, 1-2:19-25, 2015. ISSN 2352-7110. doi: https: //doi.org/10. wttp//...1116...sseeeeoo..ooo oftpsoootx.005.00666
Romeo Alexander and Dimitrios Giannakis. tperator- framework for forecasting nonlinear time series with kernel analog techniques. Physica Nonlinear Phenomena, 409:132520, 2020. ISSN 0167-2789. doi: https://doi.org// 10t.://..oo/1111111....../.. p011
P. W. Anderson. More is different. Science, 177(4047):393-300) 1972. doi: 10.1126/science.11 177. 4047.393.
David Applebaum. Lévy Processes and Stochastic Calculus. Cambridge Studies Advanced Mathematics. Cambridge University Press, edition, 2009.
Omri Azencot, N. Benjamin Erichson, Vanessa Lin, and Michael Mahoney. Forecasting sequential data using consistent koopman autoencoders. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 475-485. PMLR, 13-18 Jul 2020.
Randall Balestriero and Yann LeCun. How learning by reconstruction produces uninformative features for perception. Proceedings of the 41st International Conference on Machine Learning, ICML'24 24. JMLR.org, 2024.
Caitlin C. Bannan, Kalistyn H. Burley, Michael Chiu, Michael R. Shirts, Michael K. Gilson, and David Mobley, Blind prediction of cyclohexane-ve water distribution coefficients from the sampl5 challenge. Journal of Computer-Aided Molecular Design, doi: 10.10077410822--01---
2016. ISSN 1573-4951.
Peter Bauer, Alan Thorpe, and Gilbert Brunet. The quict revolution of numerical weather prediction. Nature, 525(7567):47-55, 2015. ISSN 1476-4687. doi: 10. 10000////tree1995555
Christopher I. Bayly, Piotr Cieplak, Wendy Cornell, and Peter A. Kollman. A well-behaved electro- static potential based method using charge restraints for deriving atomic charges: the resp model. The Journal of Physical Chemistry, 97(40): 10269-10229000 1993. doi: 10. 0.1001/10014220004
10

---

Published as a conference paper at ICLR 2026
Soumendranath Bhakat and Pär Söderhjelm. Resolving the problem of trapped water in bind-
ing cavities: prediction of host-guest binding free energies the sampl5 challenge by fun- nel metadynamics. Journal of Computer-Aided Molecular Design, 31(1):119-132, 2017. doi: 10.
Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian. Pangu-weather: A 3d high-resolution model for fast and accurate global weather forecast. arXiv preprint arXiv:2211.02556 2022.
Bjerknes. Atmospheric teleconnections from the equatorial pacific. Monthly Weather Review, 97 (3):163 172, 1969. doi: 10.11577300449991667777776666111
Luigi Bonati, Enrico Trizio, Andrea Rizzi, and Michele Parrinello. A unified framework for machine learning collective variables for enhanced sampling simulations: micolvar. The Journal of Chemical Physics, (59(1))01480) 07 2023. ISSN 0021-9606.
doi: N0.1063/5.0156330
Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15):3932-333) 2016. doi: 10. 10.1073/pnas. 151738411
Giovanni Bussi, Davide Donadio, and Michele Parrinello. Canonical sampling through velocity rescaling. The Journal of Chemical Physics, 126(1):01410,, 01 2007. ISSN 0021-96066 doi: 1..1063/..2408450
Wenju Cai, Agus Santoso, Matthew Collins, Boris Dewitte, Christina Karamperidou, Jong-Seong Kug, Matthicu Lengaigne. Michael McPhaden, Malte F. Stuecker, Andréa S. Taschetto, Axel Timmermann, Lixin Wu, Sang- Wook Ych, Guojian Wang, Benjamin Ng, Fan Jia, Yun Yang. Jun Ying, Viao-Tong/ Zheng. Tobias Bayr, Josephine R. Brown, Antonietta Capotondi, Kim M. Cobb, Bolan Gan, Tao Geng, Yoo-Geun Ham, Fei-Fei Jin, Hyun-Su Jo, Xichen Li, Xiaopei Lin, Shayne McGregor, Jae-Heung Park, Karl Stein, Kai Yang. Li Zhang, and Wenxiu Zhong. Changing el nifo-southern oscillation in warming climate. Nature Reviews Earth & Environment, 2(9): doi: 00.1038/43011--22-0009999
628-644, 2021. ISSN 2662-138X.
Christopher W. Callahan and Justin S. Mankin. Persistent effect of niño on global economic growth. Science, 380(6649); 10691069,00066000000066006911161 2023. doi: 10. 10.11266siiencceadd22953.
Lowik Chanussot, Abhishek Das, Siddharth Goyal, Thibaut Lavril, Muhammed Shuaibi, Morgane Riviere, Kevin Tran, Javier (ieras-Domingo. Caleb Ho. Weihua Hu, Aini Palizhati, Anuroop Sriram, Brandon Wood, Junwoong Yoon, Devi Parikh, C. Lawrence Zitnick, and Zachary Ulissi, Open catalyst 2020 (oc20) dataset and community challenges. ACS Catalysis, 11(10):6559-6072 2021. doi: 10. /aca//cccaaal.00052525
Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton. A simple framework for contrastive learning of visual representations. Proceedings of the 37th International Conference on Machine Learning, ICML'20. JMLR.org, 2020.
Xinlei Chen and Kaiming He. Exploring simple siamese representation learning. In 2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 15745-15753, 2021. doi: 10.1 10099/CVPPP444..2222.111.1
Matthew J. Colbrook, Lorna J. Ayton, and Máté Szöke. Residual dynamic mode decomposition: robust and verified koopmanism. Journal of Fluid Mechanics, 955:A21, 2023, doi: 10. 1017fmm 1052.1052.
Suddhasattwa Das and Dimitrios Giannakis. Koopman spectra in reproducing kernel hilbert spaces. Applied and Computational Harmonic Analysis, 49(2):573-607, 2020. ISSN 1063-5203, doi: https://oii.org org/ ac....0................oo
Henry F Diaz and Vera Markgraf. El Niño and the Southern Oscillation: multiscale variability and global and regional impacts. Cambridge University Press, 2000.
11

---

Published as a conference paper at ICLR 2026
Peter Eastman, Jason Swails, John D. Chodera, Robert T. McGibbon, Yutong Zhao, Kyle A. Beauchamp, Lee-Ping Wang. Andrew C. Simmonett, Matthew P. Harrigan, Chaya D. Stern, Rafal P. Wiewiora, Bernard R. Brooks, and Vijay S. Pande. Openmm 7: Rapid development of high performance algorithms for molecular dynamics. PLOS Computational Biology, 13(7):1-17, 07 2017. doi: 10.1371 1..771/jounll..chh. 1005659.
Marco Federici, Patrick Forré, Ryota Tomioka, and Bastiaan S. Veeling, Latent representation and simulation of markov processes via time-lagged information bottleneck. In The Twelfth International Conference on Learning Representations, 2024.
Anthony Frion, Lucas Drumetz, Mauro Dalla Mura, Guillaume Tochon, and Abdeldjalil Aïssa El Bey. Neural Koopman prior for data assimilation. JEEE Transactions on Signal Processing, 72: 4191-4206, 2024. doi: 10.1 10.1779922222244411655
Gary Froyland, Dimitrios Giannakis, Benjamin R. Lintner, Maxwell Pike, and Joanna Slawinska. Spectral analysis of climate dynamics with operator -perator-theoretic approaches. Nature Communications, 12(1))6570, 2021. ISSN 2041-1723, doi: 10.10388/441447-021226333
Damien Garreau, Wittawat Jitkrittum, and Motonobu Kanagawa. Large sample analysis of the median heuristic. arXiv preprint 1rXiv:1707.07260 2017.
Mahdi Ghorbani, Samarjeet Prasad, Jeffery B. Klauda, and Bernard R. Brooks. Graphvampnet, using graph neural networks and variational approach to markov processes for dynamical modeling of biomolecules. The Journal of Chemical Physics, 156(18):18410,, 05 2022. ISSN 0021-9606. doi: 10.1063/5.0085007.
Michael H. Glantz. and Ivan J. Ramirez. Reviewing the oceanic niño index (oni) to enhance societal readiness for niño's impacts. International Journal of Disaster Risk Science, 11(3):394 403, 2020. ISSN 2192-6395, doi: 10. 10.10077/11373-020000077 222-00075-w.
G. H. Golub and V. Percyra. The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate. SIAM Journal on Numerical Analysis, 10(2):11-433 1973. doi: 10. 10.113700 137707100066
Jean-Bastien Grill, Florian Strub, Florent Altché, Corentin Tallec, Pierre H. Richemond, Elena Buchatskaya, Carl Doersch, Bernardo Avila Pires, Zhaohan Daniel Guo, Mohammad Gheshlaghi Azar, Bilal Piot, Koray Kavukcuoglu, Rémi Munos, and Michal Valko. Bootstrap your own latent new approach to self-supervised learning. In Advances in Neural Information Processing Systems, Red Hook, NY, USA, 2020. Curran Associates Inc. ISBN 9781713829545
Nicklas Hansen, Hao Su, and Xiaolong Wang. Td-mpe2: Scalable, robust world models for continuous control, 2024.
Jeff Z. HaoChen, Colin Wei, Adrien Gaidon, and Tengyu Ma. Provable guarantees for self-supervised deep learning with spectral contrastive loss. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, volume 34, pp. 5000-5011. Curran Associates, Inc., 2021.
Jeff Z. HaoChen, Colin Wei, Ananya Kumar, and Tengyu Ma. Beyond separability: Analyzing the linear transferability of contrastive representations to related subpopulations. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35. pp. 26889-26902. Curran Associates, Inc., 2022.
M. J. Harvey, G. Giupponi, and G. De Fabritiis. Acemd: Accelerating biomolecular dynamics in the microsecond time scale. Journal of Chemical Theory and Computation, 5(6):1632-166) 1639, 2009. doi: 10. /c.10//8890 PMID: 26609855.
Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, Adrian Simmons, Comel Soci, Saleh Abdalla, Xavier Abellan, Gianpaolo Balsamo, Peter Bechtold, Gionata Biavati, Jean Bidlot, Massimo Bonavita, Giovanna De Chiara, Per Dahlgren, Dick Dee. Michail Diamantakis, Rossana Dragani, Johannes Flemming, Richard Forbes, Manuel Fuentes, Alan Geer, Leo Haimberger, Sean Healy, Robin J. Hogan, Elías Hölm, Marta Janisková, Sarah Keeley, Patrick Laloyaux, Philippe
12

---

Published as a conference paper at ICLR 2026
Lopez, Cristina Lupu, Gabor Radnoti, Patricia de Rosnay, Iryna Rozum, Freja Vamborg, Sebastien Villaume, and Jean-Noël Thépaut. The era5 global reanalysis. Quarterly Journal of the Royal Meteorological Society, (96733))199922044 June 2020. ISSN 1477-870X. doi: 10.1002/qj.3803.
Wassily Hoeffding. A Class of Statistics with Asymptotically Normal Distribution, pp. 308- 334. Springer New York, New York, NY, 1992. ISBN 978-1-4612-0919-3. doi: 10.1007/
978-1-4612-0919-5-
James W. Hurrell, M. M. Holland, P. R. Gent, S. Ghan, Jennifer E. Kay, P. J. Kushner, J.-F. Lamarque, W. G. Large, D. Lawrence, K. Lindsay, W. H. Lipscomb, M. C. Long, N. Mahowald. D. R. Marsh, R. B. Neale, P. Rasch, S. Vavrus, M. Vertenstein, D. Bader, W. D. Collins, J. J. Hack, Kiehl, and Marshall. The community earth system model: framework for collaborative research. Bulletin of the American Meteorological Society, 94(9):1339 1360, 2013. doi: 10.1175/BAMS- 12-00121.1
Minchan Jeong, J Jon Ryu, Se- Young Yun, and Gregory W Wornell. Efficient parametric svd of koopman operator for stochastic dynamical systems. arXiv preprint arXiv:2507.0722) 2025.
William L. Jorgensen, Jayaraman Chandrasekhar, Jeffry D. Madura, Roger W. Impey, and Michael Klein. Comparison of simple potential functions for simulating liquid water. The Journal of Chemical Physics, 79(2):926-95) 07 1983. ISSN 0021-9606.
doi: 10. 1063/1 445869..
Yoshinobu Kawahara. Dynamic mode decomposition with reproducing kernels for koopman spectral analysis. In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 29, Curran Associates, Inc., 2016.
Stefan Klus, Ingmar Schuster, and Krikamol Muandet. Eigendecompositions of transfer operators in reproducing kernel hilbert spaces. Journal of Nonlinear
W Science, 30(1):283-315, 2020. ISSN
1462-1467. doi: 10. 00.1007/s00332--1--099774
Dmitrii Kochkov, Janni Yuval, Ian Langmore. Peter Norgaard, Jamie Smith, Griffin Mooers, Milan Klöwer, James Lottes, Stephan Rasp, Peter Düben, Sam Hatfield, Peter Battaglia, Alvaro Sanchez- Gonzalez, Matthew Willson, Michael P. Brenner, and Stephan Hoyer. Neural general circulation models for weather and climate. Nature, 63208007):1000--0066 2024. ISSN 1476-4687. doi:
B. O. Koopman. Hamiltonian systems and transformation in Hilbert space. Proceedings of the National Academy of Sciences, 315):315-318, 1931. doi: 2..1073/paas11.55.11
Milan Korda and Igor Mezić, On convergence of extended dynamic mode decomposition to the koopman operator. Journal of Nonlinear Science, 28(2)687-710, 2018. ISSN 1432-1467. doi:
10.10077600332-017--4225
Vladimir R Kostic, Pictro Novelli, Andreas Maurer, Carlo Ciliberto, Lorenzo Rosasco, and Mas- similiano Pontil. Learning dynamical systems via koopman operator regression in reproducing kernel hilbert spaces. In Koyejo. S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 4017-4031. Curran Associates, Inc., 2022.
Vladimir R Kostie, Karim Lounici, Pietro Novelli, and Massimiliano Pontil. Sharp spectral rates for Koopman operator learning. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (eds.), Advances in Neural Information Processing Systems, volume 36, pp. 32328-32339, Curran Associates, Inc., 2023.
Vladimir R. Kostic, Karim Lounici, Grégoire Pacreau, Giacomo Turri, Pietro Novelli, and Massi- miliano Pontil. Neural conditional probability for uncertainty quantification. In A. Globerson, I.. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomezak, and C. Zhang (eds.), Advances in Neural Information Processing Systems, volume 37, pp. 60999-61038 Curran Associates, Inc., 2024a. doi: 10.52202/079011--955.
Vladimir Kostic, Pietro Novelli, Riccardo Grazzi, Karim Lounici, and massimiliano pontil. Learning invariant representations of time-homogeneout stochastic dynamical systems. In B. Kim, Y. Yue, S. Chaudhuri, K. Fragkiadaki, M. Khan, and Y. Sun (eds.), International Conference on Learning Representations, volume 2024, pp. 14329-553411 2024b.
13

---

Published as a conference paper at ICLR 2026
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadeneshelii, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: learning maps between function spaces with applications to pdes. J. Mach. Learn. Res., 24(1), January 2023. ISSN 1532-4435.
Thorsten Kurth, Shashank Subramanian, Peter Harrington, Jaideep Pathak, Morteza Mardani, David Hall, Andrea Miele, Karthik Kashinath, and Anima Anandkumar, Fourcastnet: Accelerating global high-resolution weather forecasting using adaptive fourier neural operators. Proceedings of the Platform for Advanced Scientific Computing Conference, PASC '23, New York, NY, USA, 2023. Association for Computing Machinery. ISBN 9798400701900. doi: 141144559999993555444
J. Nathan Kutz, Steven L.. Brunton, Bingni W. Brunton, and Joshua L. Proctor. Dynamic mode decomposition: data-driven modeling complex systems. SIAM, 2016.
Remi Lam, Alvaro Sanchez-Gonzaler Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, Alexander Merose, Stephan Hoyer, George Holland, Oriol Vinyals, Jacklynn Stott, Alexander Pritzel, Shakir Mohamed, and Peter Battaglia. Learning skillful medium-rann raagee global weather forecasting. Science, 382(6677): 1416-1421, 2023. doi: 10.1 httpsscciinn.a.ii333....
Karl Lapo, Sara M. Ichinaga, and J. Nathan Kutz. A method for unsupervised learning of coherent spatiotemporal patterns multiscale data. Proceedings of the National Academy of Sciences, 122 (7x:22115786122, 2025. doi: 10. 20.1073/pnas.241578111
Andrzej Lasota and Michael C. Mackey, Chaos, Fractals, and Noise: Stochastic Aspects of Dynamics. Springer New York, New York, NY, 1994. ISBN 978-1-4612-4286-4. doi: 10.1007/978- 1-4612-4286-4.
Samuel Lavoie, Christos Tsirigotis, Max Schwarzer, Ankit Vani, Michael Noukhovitch, Kenji Kawaguchi, and Aaron Courville. Simplicial embeddings in supervipe learning and down- stream classification. arXiv preprint arXiv:2204.c 00616, 2022.
Qianxiao Li, Felix Dietrich, Erik M. Bollt, and loannis 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, 271))11031 10 2017. ISSN 1054-1500. doi: 10. 10.31/ 49993544.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint urXi:22010.08895, 2020.
Vittorio Limongelli, Massimiliano Bonomi, and Michele Parrinello. Funnel metadynamics as accurate binding free-energy method. Proceedings of the National Academy of Sciences, 110(16):6358- 6363, 3013. doi: 10.1073/pnas. 130333113 10.
Kresten Lindorff-Larseen Stefano Piana, Ron O. Dror, and David E. Shaw. How fast-folding proteins fold. Science, 334(6055):517-520, 201 doi: 10. 1166science. 1208555
1208351.
Edward N. Lorenz. Deterministic nonperiodic flow. Journal of Atmospheric Sciences, 20(2):130- 141, 1963. doi: 10.
Yiwei Lu, Guojun Zhang, Sun Sun, Hongyu Guo, and Yaoliang Yu. f-micl: Understanding and generalizing infonce-based contrastive learning. arXiv preprint arXiv:2402. 2001:00000.000500 2024.
Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton. Deep learning for universal linear embeddings of nonlinear dynamics. Nature Communications, 91))4955,, 2018. ISSN 2041-1723. doi: 10. 1038/ 111467-018-0721000
Xubo Lyu, Hanyang Hu, Seth Siriya, Ye Pu, and Mo Chen. Task-oriented Koopman- based control with contrastive encoder. In 7th Annual Conference Robot Learning, 2023.
Andreas Mardt, Luca Pasquali, Hao Wu, and Frank Noé. Vampnets for deep learning of molec- ular kinetics. Nature Communications, 9(1):5, 2018. ISSN 2041-1723. doi: 10.1038/
541467-017-02388-1
14

---

Published as a conference paper at ICLR 2026
Michael J. McPhaden, Stephen E. Zebiak, and Michael H. Glantz. Enso integrating concept in earth science. Science, 314(5806):1740-17728 2006. doi: 10.1126/science. 1132588.
Giacomo Meanti, Antoine Chatalic, Vladimir R Kostic, Pietro Novelli, Massimiliano Pontil, and Lorenzo Rosasco. Estimating koopman operators with sketching to provably learn large scale dynamical systems. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (eds.), Advances Neural Information Processing Systems, volume 36, pp. 77242-77276. Curran Associates, Inc., 2023.
Igor Mezić. Spectral properties of dynamical systems, model reduction and decompositions, Nonlin- ear Dynamics, (1(1):309-323) 325, 2005. ISSN 1573-269X. doi: 10. 1007/s00ss 1071-005-22---
Thomas Minka. Old and new matrix algebra useful for statistics. 4, 2000.
Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida. Spectral normalization for 2018.
generative adversarial networks. arXiv preprint
L. Molgedey and H. G. Schuster. Separation of mixture of independent signals using time delayed correlations. Physical Review Letters, 72(23):363-3337, June 1994. ISSN 0031-9007. doi: 10.1 10.10000/hyseevlett22...55.
Juan Nathaniel, Yongquan Qu, Tung Nguyen, Sungduk Yu, Julius Busecke, Aditya Grover, and Pierre Gentine. Chaosbench: multi-channel, physics-basee benchmark for mubseasonal-to-seanee 2024.
climate prediction. arXiv preprint
NOAA's climate prediction center NCP Center. Climate Prediction Center ONI
origin..pe.ncee.nnaagggv https://origin: cpc. nwwp.ccppp.......c noaa gov/products//
php
[Accessed 09-05-2025).
Frank Noé and Feliks 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.
Pietro Novelli, Luigi Bonati, Massimiliano Pontil, and Michele Parrinello. Characterizing metastable states with the help of machine learning. Journal of Chemical Theory and Computation, 18(9): 5195-5202, 2022. doi: 10. 1021/acc..jcccc2200333. PMID: 35920063.
Pietro Novelli, Marco Pratticò, Massimiliano Pontil, and Carlo Ciliberto. Operator world models for reinforcement learning. In A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomezak, and C. Zhang (eds.), Advances in Neural Information Processing Systems, volume 37, pp. 111432- 11 1463.. Curran Associates, Inc., 2024. doi: 10. 10.522022079011--3339
Feliks Nüske, Sebastian Peitz, Friedrich Philipp, Manuel Schaller, and Karl Worthmann. Finite-dar data error bounds for koopman-based prediction and control. Journal of Nonlinear Science, 33(1):14, 2022. ISSN 1432-1467. doi: 10. 10000/000032-02--9996666
Samuel E. Otto and Clarence W. Rowley. Linearly recurrent autoencoder networks for learning dynamics, SIAM Journal on Applied Dynamical Systems, 18(1):558-59% 2019, doi: 10.1137/ 18M1177846.
Bette Otto-Bliesner, Esther C. Brady, John Fasullo, Alexandra Jahn, Laura Landrum, Samantha Stevenson, Nan Rosenbloom, Andrew Mai, and Gary Strand. Climate variability and change since 850 ce: An ensemble approach with the community earth system model. Bulletin of the American Meteorological Society, 97(5):735 754, 2016. doi: 10.1 1117//B/M-DD 10-0-44-023311.
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzaler and Peter Battaglia. Learning mesh-based simulation with graph networks. International Conference on Learning Representations, 2021.
S. G. H. Philander. El niño southern oscillation phenomena. Nature, 302(5906):2953300 1983. ISSN 1476-4687. doi: 10. 00.103/302295a0.
Guillermo Pérez-Hemándee Fabian Paul, Toni Giorgino, Gianni De Fabritis, and Frank Noé. Identification of slow molecular order parameters for markov model construction. The Journal of Chemical Physics, 139(1 (9(11:055122 07 2013. ISSN 0021-9606. doi: 10.1063/1.44 1.81 4811489
15

---

Published as a conference paper at ICLR 2026
Michael Reed and Barry Simon. Methods of Modern Mathematical Physics. Academic Press, 1972. ISBN 9780125850018.
Tongzheng Ren, Tianjun Zhang, Lisa Lee, Joseph E. Gonzalez, Dale Schuurmans, and Bo Dai. Spectral decomposition representation for reinforcement learning. In The Eleventh International Conference on Learning Representations, 2023.
Valerio Rizzi, Luigi Bonati, Narjes Ansari, and Michele Parrinello. The role of water in host- guest interaction. Nature Communications, 12(1):93, 2021. ISSN 2041-1723. doi: 10. 10.1038/ s411677-20-2031044
Preston Rozwood, Edward Mehrez, Ludger Pachler, Wen Sun, and Steven Brunton. Koopman- reinforcement learning. In NeurlPS 2023 Al for Science Workshop, 2023.
J. Jon Ryu, Xiangxiang Xu, H. S. Melihcan Erol, Yuheng Bu, Lizhong Zheng, and Gregory W. Wornell Operator svd with neural networks via nested low-rank approximation, 2024.
Alvaro Sanzhez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning simulate complex physics with graph networks. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference on Machine Learning. volume 119 of Proceedings of Machine Learning Research, pp. 8459-8468. PMLR, 13-18 Jul 2020.
Peter J. Schmid. Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics, 656:5-28, 2010. doi: 10.1017/500221 12010001217.
Kristof Schütt, Pieter-Jan Kindermans, Huziel Enoc Sauceda Felix, Stefan Chmiela, Alexandre Tkatchenko, and Klaus- Müller. Schnet: continuous-fllte convolutional neural network for modeling quantum interactions. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30, Curran Associates, Inc., 2017.
Ch. Schütte. W. Huisinga, and P. Deuflhard. Transfer operator approach to conformational dynamics in biomolecular systems. In Bernold Fiedler (ed.), Ergodic Theory, Analysis, and Efficient Simulation of Dynamical Systems, pp. 191-223, Berlin, Heidelberg, 2001. Springer Berlin Heidelberg.
Max Schwarzer, Ankesh Anand, Rishab Goel, R Devon Hjelm, Aaron Courville, and Philip Bach- man. Data-efficient reinforcement learning with self preficredittive representations. arXiv preprint arXiv:2007.00599, 2020.
Shai Shaley-Shwariz and Shai Ben-David. Understanding Machine Learning: From Theory to Algorithms. Cambridge University Press, 2014.
Hythem Sidky, Wei Chen, and Andrew L. Ferguson. High-resolution markov state models for the dynamics of trp-cage miniprotein constructed over slow folding modes identified by state-free reversible vampnets. The Journal of Physical Chemistry B, 123(38):7999 8009, 2019. ISSN 1520-6106. doi: 10.102 htt.bessppp.bbb0557...
Haotian Sun, Antoine Moulin, Tongzheng Ren, Arthur Gretton, and Bo Dai. Spectral representation for causal estimation with hidden confounders. Yingzhen Li, Stephan Mandt, Shipra Agrawal, and Emtiyaz Khan (eds.), Proceedings of The 28th International Conference Artificial Intelli- gence and Statistics, volume 258 of Proceedings of Machine Learning Research, pp. 2719-2727. PMLR, 03-05 May 2025.
R.S. Sutton and A.G. Barto. Reinforcement Learning: An Introduction, volume 9, MIT Press, 1998. doi: 19.1109/N...999..12211 1998.7
Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi. Learning koopman invariant subspaces for dynamic mode decomposition. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
Gareth A. Tribello, Massimiliano Bonomi, Davide Branduardi, Carlo Camilloni, and Giovanni Bussi. Plumed 2: New feathers for an old bird. Computer Physics Communications, 185(2):6046633, 2014. ISSN 0010-4655. doi: http://oooo./.1................0
16

---

Published as a conference paper at ICLR 2026
Warwick Tucker. The lorenz attractor exists. Comptes Rendus de l'Académie des Sciences Series Mathemattics, 128(1))1177-120 1999. ISSN 0764-4442. doi: https://doi.or/10..00116 50764-4442(99)8803333
Aaron van den Oord, Yazhe Li, and Oriol Vinyals. Representation learning with contrastive predictive coding, 2019.
Junmei Wang, Romain M. Wolf, James W. Caldwell, Peter A. Kollman, and David A. Case. Devel- opment and testing of a general amber force field. Journal of Computational Chemistry, 25(9):
1157-1174, 2004. doi: https://doi.org/II 100.000000..0003.......
Ziyu Wang, Yucen Luo, Yueru Li, Jun Zhu, and Bernhard Schölkopf. Spectral representation learning for conditional moment models. arXiv preprint arXiv:2210..6552. 2022.
Christoph Wehmeyer and Frank Noé. Time-lagged autoencoders: Deep learning of slow collective variables for molecular kineties. The Journal of Chemical Physics, 148(24):241703, 03 2018. ISSN 0021-9606. doi: 10.1063/1.5011399
Matthew O. Williams, Clarence W. Rowley, and Toannis G. Kevrekidis. A kernel-based method for data-driven koopman spectral analysis. Journal of Computational Dynamics, 2(2):247-265, 2015. ISSN
doi: 0..3334/jcd.201000
Hao Wu and Frank Noé. Variational approach for learning markov processes from time series data. Journal of Nonlinear Science, 30(1):23-66, 2020. ISSN 1432-1467. doi: 10.1007/
s00332-019-09566--9
Enoch Yeung. Soumya Kundu, and Nathan Hodas. Learning deep neural network representations for koopman operators of nonlinear dynamical systems. In 2019 American Control Conference (ACC), pp. 4832 4839, 2019. doi: N..2319/ACC.222 188199911333.
Jian Yin, Niel M. Henriksen, David R. Slochower, Michael R. Shirts, Michael W. Chiu, David Mobley, and Michael K. Gilson. Overview of the sampl5 host guest challenge: Are we doing better? Journal of Computer-Aided Molecular Design, 31(1):1-19, 2017. ISSN 1573-4951. doi: 1000007/11822-016-9919
Jure Zbontar, Li Jing, Ishan Misra, Yann LeCun, and Stephane Deny. Barlow twins: Self-supervised learning via redundancy reduction. In Marina Meila and Tong Zhang (eds.), Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 12310-12320. PMIR, 18-24 Jul 2021.
Jintu Zhang, Luigi Bonati, Enrico Trizio, Odin Zhang, Yu Kang, TingJun Hou, and Michele Parrinello. Descriptor-free collective variables from geometric graph neural networks. Journal of Chemical Theory and Computation, 20024); 2024. ISSN 1549-9618. doi: 10.1021/acs /acs//c..ec 4c01 197.
Jiawei Zhuang, raphael dussin, David Huard, Pascal Bourgault, Anderson Banihirwe, Stephane Raynaud, Brewster Malevich, Martin Schupfner, Filipe, Charles Gauthier, Sam Levang, André Jüling, Mattia Almansi, RichardScottOZ, RondeauG, Stephan Rasp, Trevor James Smith, Ben Mares, Jemma Stachelek, Matthew Plough, Pierre, Ray Bell, Romain Caneill, and Xianxiang Li. pangeo-data/xesmf: v0.8.10. 10, April 2025.
H. Zuo, M. A. Balmaseda, S. Tietsche, K. Mogensen, and M. Mayer. The ecmwf operational ensemble seeualyis--nnllye system for ocean and sea ice: description of the system and assessment. Ocean Science, 15(3):779-800 2019. doi: 10.5 15-779-2019.
17

---

Published as a conference paper at ICLR 2026
we made aware of a concurrent preprint (Jeong et al. 2025 proposing a similar methodology to learn the evolution operators of stochastic dynamical systems. In particular, Jeong et al. (2025) proposes a variation of the loss 7 with the more agnostic choice (ap ), (((++1)) with both yp, v trainable. Our work, however, differs from Jeong et al. (2025) two key aspects. First, as proved in Sec. [31] our parametrization
Concurrent work. During the preparation of this submission (September 2025).
(p(æz), P(((T++1) recovers the least-squares evolution operator learning framework (Lemma D).
for which both approximation (Korda & Mezić, 2018) and statistical learning Kostic et al. (2022)
Nüske et al. (2022) results have been proved. Furthermore, as discussed in Experiment Sec. 4.1 our
parametrization results in higher-implied timescales which, according to the variational principle (Noé)
& Nüske 2013 Wu & Noé, 2020), are associated with an improved accuracy in the estimation of the evolution operator. Our parametrization is also directly linked the VAMP score through Lemma
As second point of departure from Jeong et al. (2025), our experiments in Sec. 4 focus high-
dimensional dynamical systems, and to the best of our knowledge demonstrate evolution operator learning at scales never reached so far. Further, drawing on the connection of (7) to self-supervised
learning- which not acknowledged in Jeong et al. (2025) despite the prior work Wang et al. (2022)
experiments Sec. 4.2 and Sec. 4.3] demonstrate that representations learned via (7) can be successfully transferred to new and unseen dynamical systems.
A PROOFS OF THE THEORETICAL CLAIMS.
Define v P[X,), the distribution of the initial states in our dataset, and μ Eg P[X++1] | the distribution of the evolved states. In practice, can be the following:
(p-|=))
If a simulator available, v can be any distribution of initial states, and f is obtained by single step of the simulator on data from v.
If one samples trajectories of length from an initial distribution P[Xi]. then v =
T D P[X.).
If - as in molecular dynamics, or the Lorenz 63 example below - one samples from an = P[X++1] n. one has v / π.
invariant distribution ग such that P[X.] π
The evolution operator E maps functions from (H) into (v),,, that is E (²(x)) + (v))). Notice that since we allow for general initial and evolved distributions, respectively and μ., our method does not require E to be associated to stationary, nor ergodic dynamical system, as often the case in the theoretical literature, see e.g. (Mezić 2005 Section 2.3) Kostic et al. (2022). We will now show this simple equivalence:
Lemma 1. Let y + Rd be an encoder whose components are square-integrable with respect to both μ and and let E be equivalent to the following operator learning loss:
evolution operator. Then, the loss function 7 is
elp, P) = I|E
Σ 1
se & Pijssslis
Proof. The Lemma was already proved in Wang et al. (2022 Lemma 4.1). or Kostic et al. 2202aa
Theorem 1). Here we provide self-contained proof. Since the encoder is square-integrable in both μa and si E (1) and (v) we can define the linear operators
4 (p) + Rdd
f fH = (1f 9i) (P()))
8 Rd 13(v))
ZH (()))z
By direct substitution of the definition above, follows that
E Φつ PT
11 = ||E Σ ss 6 Pliss
VVersion I of Jul 2025.
18

---

Published as a conference paper at ICLR 2026
Now, let's notice that by direct calculation one obtains
t E, [+(())((x)]]
8, Ep [+(x))<v)]]
where p(dz, dy) = p(dy)=)))))))))))))))))p)) ())(rr) is the joint distribution of (Xt,X++)) By the definition of the Hilbert-Schmidt norm, have
IE P iiss
|E H5 - 2Tr [E "E"0 P + Tr o P P
||E| Its 2Tr ( Tr o P PI
||E|| AAs - 2F (x.w)~ Tr (()))(()) P]] + (Z.)))uV [Tr PT
||E|iss 2E (E. [((()) Py(())) + ((..))uuv (v(())
where we repeatedly used the cyclic property of the trace.
The following Lemma shows that when P is optimal with respect to then it recovers the least squares estimator 3
Lemma 2. For any fixed i. the predictor P minimizing 7 can be computed in closed form P. Cx CxyCr' and the model for the evolution operator is given by
Egp = P.Cr Cx Cxy
Eq.
with X 0,
coinciding with the least-squares estimator
Proof. The proof follows by noticing that c({p, P) is convex in P. Taking the gradient (see, for example (Minka 2000)) one has:
rel(p, P) -2F( [x())<(x) 2F (aw))oov (e(v)<(U) pT PT(())(()T]]
= = PCy
+
As the problem is convex, the global minimum P. is attained when Vpe(p, 0. This condition is equivalent to solve the equation
2Cyx + 2Cy Px =
By multiplying the expression above by Cx' on the right and Cj on the left, re-arranging it, and taking the transpose everything, we finally get
P. Cx Cxy
The following Lemma shows the equivalence of and the VAMP-2 loss of Wu Noč (2000); Mardt et al. (2018)
Lemma 3. For any fixed let P. the optimal predictor of E(sp, P), in Lemma Then, the following holds true:
Ptp, = - 1/2 12 1/2 2 -VMMP2(0).
Proof. By noticing that the loss function can be equivalently rewritten as
c(p, P) Tr(PTT CPPC 2PCyx].
and substituting the optimal predictor P. Cx' from Lemma inside this expression, we immediately obtain the identity.
19

---

Published as a conference paper at ICLR 2026
A.1 ON THE HILBERT-SCHMIDT SCHMIDT ASSUMPTION, AND BEYOND
In the main text, we assumed the evolution operator E to be Hilbert-Schmidt, which immediately guarantees the well-posedness of the proposed loss function 77. Lemma indeed, implies that for any Hilbert-Schmidt E it holds
s(p, -) -- @ Pippy Its < 00. 1.j
The Hilbert-Schmide assumption is valid in broad class of stochastic systems, particularly when the transition kernel exhibits smoothing properties. atomistic simulations, for instance, the presence of finite temperature results in a Gaussian smoothing that makes E Hilbert-Schmidt, As illustrative example, consider the overdamped Langevin dynamics
X X, VV(X,) N(,,
where V is potential function, and ( or) denotes isotropic Gaussian with mean 0 and variance a proportional to the system's 's temperature. Assuming that the data is sampled from the equilibrium distribution π (x))dx x e (-9V(x) dx, we compute
||E| ||E||As hs
Σ 1E Σ
/
/
ply ()) ev)) (dy)
T (dx)
= =f ply x²² (z) dy dz x ʃ /|oxp( ly - ((())(()))))))))))))))))))))(() dy oc
20²
where e; are elements of an orthonormal basis of (TT) and the third equality follows from Parseval's identity.
The Sihtert--chhidt assumption, however, is violated in important deterministic dynamical systems, such those governed by fluid dynamics equations Mezić (2005). Remarkably, the empirical loss (8) interpretation when is merely bounded operator. Indeed, let P. (n) (1) denote the orthogonal projector onto the subspace spanned by the encoder 4 d where D, is in Lemma we define the abstract loss
still admits a precise operator-
p. By definition, P D
With slight abuse of the notation in
c P) = 2r(E
I p Pase,liss
whose empirical estimator exactly coincides with (8). the loss which we actually optimized in our experiments. Now, since y spans finte-diieensional subspace, both P and EP,, e are finite rank, hence Hilbert-Schmidt, In particular, EPg, so is the restriction of to the subspace generated by y. An immediate calculation shows that
( P)
||E e
||EP IIAs:
(11)
where besides basic algebraic manipulations the result is obtained using the ciclicity of the trace, and the relation P. - ( (,,,,,,) +,,) 4,,-
The first term on the right-hand side of (11) is familiar, and represents the error incurred by the model Eig PP: Pujyi in approximating the restriction of E the subspace spanned by 40. This error can be linked to the least-squares approach discussed in Sec. [ see, for example Korda & Mezic (2018 Theorem 1). Minimizing it with respect to leads to representation collapse, since (ix) = for all trivially minimizes it. Our learning objective (T)), instead, avoids collapse through the second term, - -iIEP,, plis. which can be interpreted as follows. Without loss generality, write
P = >i e, e,, where e; form orthonormal basis of span((1........ Pa). By definition of Hilbert-Schmidt norm one has
|I lits
Σ
|E Σ (E.))))
Now notice that interpreting e; as probe we have at our disposal to observe the system the term (z) quantifies the "dynamical response" read by our probe, given that the system was prepared
Functions of the state of the systems commonly referred to as observables,
20

---

Published as a conference paper at ICLR 2026
Table 1: Forecasting errors and training times for the Lorenz '63 example (20 independent runs). Note that for LinLS and KRR is reported the total fitting time while for the other methods the epoch time is reported. Best results are highlighted bold.
Ours
LinLS
KRR
VAMPNets
DPNets
DAE
CAE
RMSE (x 10 3)
0.49.10.24
1.29
2.10
0.78
0.58
0.77.1.11
2.58
MAE (x10
0.32:10.4
0.84=0.00
1.27=0000
0.46:20.8
0.3620.08
0,55+0008
1.95=0.14
Time (ms)
181.1440.1
4t.1
(25.3+0.2)10
165,5+108
190.7#41.
166.8+9.10
408.5=41.9
to be in state x. The quantity EEPolliS: therefore, measures the average strength of such responses, implying that the second term in the loss promotes encoders yp whose span captures observables with the highest possible dynamical variability. To close the discussion, we highlight that for Hilbert- Schmidt operators, the observables with the highest dynamical response are precisely the leading singular functions, and the loss function (7) is indeed minimized when ip spans the leading singular space of the evolution operator E, see (Kostic et al. 2024a, Theorem 1).
B EXPERIMENTAL DETAILS
The experiments have been performed on the following hardware:
Node with 32 cores Ice Lake at 2.60 GHz, x NVIDIA Ampere A100 GPUs, 64 GB and 512 GB RAM.
Node with 20 cores Xeon Silver 4210 2.20 GHz, NVIDIA Tesla V100 GPUs, 16 GB and 384 GB RAM.
A workstation equipped with i7-5930K CPU at 3.50 GHz, 2 x NVIDIA GeForce GTX TITAN X GPUs, 12 GB and 32 GB of RAM.
B.1 ADDITIONAL EXPERIMENT: LORENZ '63
We evaluated our method on the Lorenz '63 system (Lorenz 1963), classical example of chaotic dynamical system governed by three coupled ordinary differential equations. To validate the per- formance of our approach, we tested on a one-step-ahead forecasting task, and we analyzed the learned dynamical modes. Because of the fow-dimensionality of the state zt, appended it non-learnable feature the encoder (xx,) ) = [xL æt] to ensure that the forecasting target-the state itself-lies in the linear space of functions spanned by yp by design. The learnable part of the encoder consisted of small multi-layer perceptron (MLP).
In Tab. [ we compare the performance of the estimator E.p from with encoder oo trained according to Alg. against several baseline models. These include Linear Least Squares (LinLS)) Kernel Ridge Regression (Kostic et al.) 2022) (KRR) with Gaussian kernel, VAMPNets (Mardt
et al. 2018). DPNets (Kostic et al. 2024bb). Dynamic Autoencoder @Lusch et al. 2018) (DAE), and Consistent Autoencoder (Azencot et al. 2020) (CAE). To ensure a fair comparison, we matched the encoder architecture for VAMPNets, DPNets, DAE and CAE, while decoders of DAE and CAE were defined as MLPs symmetric their respective encoders. For KRR, the rank was set equal to the latent dimensionality used in the deep learning models.
The results on the forecasting task demonstrate that, although our model is not specifically designed for prediction, it achieves the best performance among all considered methods. Finally, we verified that the leading eigenfunctions obtained by our approach correctly identify coherent sets on the stable attractor (see Fig. 4
We generated single long trajectory of 15,000 time steps using the kooplearn 1.1.3 implementation of Lorenz '63 dynamical system, with default parameters. To ensure conver- gence to a system's attractor (Tucker 1999). we discarded the first 1,000 time steps. Also, to obtain approximately -ndeeendeenddnt segments for training, validation and testing, further discarded 1,000 time steps between each split. In total, 10,000 time steps were used for training, and 1,000 time steps each for validation and testing.
Training details.
21

---

Published as a conference paper at ICLR 2026
As 087800133
As 1--099881771
0.844
0.005
Our
At 0 9961--0056
A 0.50
LinLS
A 1.000
As O 994--9 093)
As 0.996 DAG
0.9162-0..01
Au 6.973 0218
As 6.900-6.303
KAR
1s 1.000
A 9.08-000931
As 0.801:-0.122
Pa 89944-00.0
VANPNets
As 1.000
As 0.905-30011
As 8092--0000)
Au n 0.001-01000
As (501
As 0.00
0.581
DPNets
Ju 1.000
A, 099311:771)
As 0.002.0.000
As 0.00
0.7755.1771
Au 0.97780.143
- 0.052
DNE
1 0.9932/0011
0000201120
As 0es5e0 173
J.- 0.01500003
As 000
As 00320.001)
CAE
Figure 4: Leading eigenfunctions computed by and baseline approaches. Each corresponds to different method, and each column shows an eigenfunction ordered by decreasing eigenvalue magnitude.
Our encoder consisted of an MLP with an input layer of size 3, two hidden layers with units each, and an 8-dimensional latent space, using ReLU activation functions. The model was trained for 100 epochs using the AdamW optimizer, with an initial learning rate of 10-3 decayed 10-4 via cosine schedule, a batch size of 512, and lag time of 10 time steps.
Baseline methods. We compared our approach against the following baseline models:
Linear Least Squares (LinLS), linear regression model trained directly on the raw input features without any nonlinear transformation or latent representation.
Kernel Ridge Regression (KRR) (Kostic et al., 2022). We trained KRR model with Gaussian kernel, using the bandwidth estimated via the median heuristic (Garreau et al. 2017). The model was trained with rank of 8, Tikhonov regularization parameter of 10 and using Arnoldi iterations.
VAMPNets (Mardt et al., 2018). Trained using the same MLP encoder as ours, with the VAMP-2 loss loss and centered covariances.
DPNets Kostic et al. 20245). Trained using the same MLP encoder as ours, with the relaxed DP loss and centered covariances.
22

---

Published as a conference paper at ICLR 2026
Dynamic Autoencoder (DAE) (Lusch et al., 2018). Trained with the same MLP encoder architecture as in our approach; the decoder was defined symmetrically. The loss components for reconstruction, prediction, and linear evolution were equally weighted (all set to 1).
Consistent Autoencoder (Azencot et al., 2020). Trained with the same MLP encoder architecture as in our approach; the decoder was defined symmetrically. The CAE loss weights for reconstruction, prediction, backward prediction, linear evolution, and consistency were all set to 1.
For all deep learning-based baselines (VAMPNets, DPNets, DAE, and CAE), models were trained for 100 epochs using a batch size of 512, and lag time of 10 time steps. VAMPNets and DPNets used the Adam W optimizer with learning rate of 10 4 and 10 2, respectively; DAE and CAE used the Adam optimizer with learning rate of 10 3 All baselines were implemented using koop 1 ea rn 1.1.3.
Additional analysis. In Fig. we show the leading eigenfunctions of the transfer operators computed using our method and the baseline approaches. These visualizations highlight qualitative differences in the learned spectral structures, offering insight into the dynamics captured by each method. The leading eigenfunction of KRR, VAMPNets, DPNets, and DAE is constant and associated with the stable attractor. Our method, LinLS, and KRR, find an cigenfunction with eigenvalue se .996 which clearly separates the two lobes of the attractor.
B.2 PROTHIN FOLDING
Training details. We used data from (Lindorfi-Larsen et al. 2011). which can be requested directly to De Shaw Research and are available without charge for academic usage. Our encoder consisted of SchNet (Schütt et al., 2017) graph neural network with interaction blocks, 16 RBF functions and an hidden dimension of 64. The model was trained with an Adam optimizer with starting learning var rate of 10-2 decaying to 10-1 with a cosine schedule, using the ml colvar (Bonati et al. 2023 library.
Additional analysis. To understand what mode is associated the leading eigenfunction 17,, in Fig. [ we correlated it with two physical quantities associated with the folding, which are the Root- Beea-SSuare--a -eviiatin (RMSD) and the Radius of Gyration, see Fig. 5 Furthermore, to obtain a finer understanding, we used sparse linear models to approximate the CVs via LASSO regression. This yields a surrogate model which is a linear combination of a few physical descriptors, hence interpretable. To choose the regularization strength, computed the Mean Square Error of the surrogate model versus the number of features, see Fig. 6
We performed LASSO regression
set of contact functions determining the presence of hydrogen bonds. The features selected by this procedure, as well as a snapshot of the protein where these features are highlighted, are reported in Fig. 7 Interestingly, some of the selected features pertain to side-chain interactions, a piece of information that would have been impossible to get using only C. atoms to train the encoder.
B.3 LIGAND BINDING
Simulations details. We selected subset of host-guest systems for the SAMPL.5 challenge (Bannan)
cct al. 2016 Yin et al.) 2017)) to evaluate our method's performance, including three ligands (GI,
G2, G3) and the octa-acid calixarane host (OAMe). Simulations were run in GROMACS 2024.5
(Abraham et al., 2015) patched with PLUMED 2.9.3 (Tribello et al. 2014). Systems were built using
the GAFF (Wang al. 2004) force field with RESP (Bayly al.) 1993) charges, solvated in a cubic 1983) water box 40.27 À length, containing 2100 water molecules. System charge balanced with Na ions. Our timestep is f and the temperature is set to 300 via a velocity rescale thermostat (Bussi al. 2007) with coupling time of 0.1 ps. All simulations aligned the host's ' vertical axis / with the box axis and centered coordinates on virtual atom V1. All production simulations were initiated from the dissociated state of each ligand. Trajectories were terminated when the ligand fully rebounded into the binding pocket (defined as host-guest distance h < Å). For each ligand, we performed 10 independent production trajectories, with coordinates saved every 500 steps.
TIP3P (Jorgensen et al.
23

---

Published as a conference paper at ICLR 2026
Interpretation via LASSO regression
1.8-
0.20
Maporau
Σ dED
Waz
Raliia of Gyralio [M]
1.4
04
Meat, Squar Eree
0.15
C
19, Upparaa + Iholla
0.2
0.0
0.10
-0.2
1.0
0.06
10 Number of features
25
0.25
0.50
0.75
1.00
1.25
1.50
Sqgare Destation
Figure 6: MSE of approximating Vi by LASSO regression meaningful physical de- scriptors. For Trp-cage we constructed a library of hydrogen-bond contact functions. The selected descriptors are reported in Fig.
Figure 5: The value of the leading eigenfunction VI the evolution operator is highly correlated with the RMSD and Radius of Gyration of the Trp-cage protein.
Physical deseriptors (H-bonds)
Normalized Coefficient
GLY10-O SERI3-N
0.307
GLY11-0 ARGI6-N
0.294
TRP6-0 GLY11-N
0.170
TRP6-NEIs (sidechain) ARGI6-0
0.109
GLN5-0 ASP9-N
0.073
TRP6-NEIs (sidechain) PRO17-0
0.044
TRP6-NEIs (sidechain) PRO18-N
0.002
Figure 7: Normalized hydrogen-bow bond coefficients selected by the LASSO model (left) and represen- tative structural snapshot (right) with the features highlighted.
In our simulations, we applied funnel restraint (Limongelli et al. 2013) to limit the conformational space explored by the ligand in the unbound state, turn accelerating the binding process. The parameters are identical to those used in previous studies (Rizzi et al. 2021). We define h as the projection of each ligand along the binding axis, treated as its radial component. For h 2 10 Å, the funnel surface is a cylinder with radius Reyt 2 Á along the vertical axis. For h < 10 Å, the funnel opens into a conical shape with 45 angle, defined by r 12 h. The force acting on displacement v from the funnel surface is harmonic:
Frunnet -k with kF 20 kJ mol À A-
An additional harmonic restraint prevents the ligand from escaping too far from the host, enforcing an upper boundary:
Fupper ky 18) for h 18Å, with by 40 mol-! A
The data will released to ensure the reproducibility of the experiment.
Training details. Our encoder consisted of a SchNet (Shütt et al. (017)) graph neural network with 3 interaction blocks, 16 RBF functions, and hidden dimension of 64 with an AdamW optimizer with starting learning rate of 10-2 decaying to 10-4 with cosine schedule.
Additional analysis. In Fig. g we inspect the two leading eigenfunctions of the evolution operator by correlating them with two physical descriptors connected to the binding: the distance along the
24

---

Published as a conference paper at ICLR 2026
Vi semi- houii-hound-e -ooodd--boni +bound
Va unbound (semm)))sshhonn hound
(Lairing
Caatiing genach)
1.5
Z distance
1.0
0.5
Cleprescutatio transfe)
Repossmratiio transle)
1.5
distance [A]
1.0
0.5
Normalized W value
50
100
150
50
100
150
Angle [deg]
Angle [deg]
Figure 8: Analysis of the leading eigenfunctions in the space of the host-guest distance 2 and the ligand orientation for 1 (left) and v22 (right). The first row contains the results obtained from training from scratch the representation, while the second row contains the case which it transferred from other systems.
direction between the center of of the host and the guest and the angle of the ligand with respect to the z axis (see figure in the inset). These results allow us to correlate the y eigenfunction the transition between the semi-bound pose to the native one, which is due to the presence of trapped water molecules inside the pocket (Rizzi et al., 2021 Bhakat & Söderhjelm 2017). The second eigenfunction V instead associated with the binding process. Furthermore, we compared the eigenfunctions obtained by training the representation from scratch on the G2 ligand with the case in which this transferred from other ligands (Gl and G3), obtaining a remarkable agreement. The ligands GI, G2, and G3 are represented in Fig.
Representation learning
Transfer
Gl
63
G2
Figure The three different molecules studied in the ligand-binding experiment.
B.4 CLIMATE MODELING
Datasets.
Following the methodology outlined in (NCP Center), compute SST* from
surface temperature (SST) data provided by the ORAS5 reanalysis (Zuo et al.) 2019), as made available through the ChaosBench dataset (Nathaniel et al. 2024). The dataset spans 45-year period
25

---

Published as a conference paper at ICLR 2026
(1979-200) (023) spatial resolution of 1.5", resulting in time series of 540 monthly snapshots, each with dimensions 121 x 240. Data from 1979 to 2016 was used for training, while the 2017-2023 period was reserved for validation.
For the transfer learning task. we employed simulations from the CESM Last Mil- lennium Ensemble project (Otto-Bliesner al.) 2016), spanning the years 850-2006 (files b.e11. BBL.11.LLLMRCCCNN..OMC 001 ottpr:aas....................... bege__a.aaacccccc..........o..o.
www..S....0000-089111110 00 htw.::p...........0000000100000.
to
available here: To ensure spatial compatibility between synthetic and observational data, the CESM SST fields were regridded onto the same 1.5 x 1.5 regular latituee grid of ORASS using the XESMF (Zhuang et al. 2025 python package.
b.e11..LWTRCCCN.
https:// ucar http://w.....a...eeeaaaaatttette
Training details. Both models trained with CESM and ORASS data use lightweight CNN encoder with four convolutional layers, batch normalization, and max pooling. A masked global average pooling layer, leveraging binary land-ocean mask, ensures only ocean data contribute the output representation. The pooled features are mapped through a final linear embedding layer.
For the CESM model, the linear layer P maps to 128-dimensional latent space. Training included simplicial normalization (Lavoie et al., 2022) (dimension 2), spectral normalization (Miyato et al.) 2018) the linear layer P, gradient clipping (max norm 0.2), a lag time of one month, 100 epochs, AdamW optimizer, and a cosine-dccayed learning rate from 10 to 10-5 with a batch size of 64. Leading eigenvalues of the transfer operator reported in Tab.
Table 2: Leading eigenvalues of the transfer operator learned ORAS5 data with ş trained on CESM data. Each eigenvalue is expressed in terms of its real (Re), imaginary (Im), and absolute (Abs) components. The associated decorrelation times and oscillation frequencies (in years) are also reported. Eigenvalues are listed in descending order with respect to their absolute value, and those with decorrelation time shorter than 1/12 years, i.e., the sampling frequency, were discarded.
Idx
Re
Im
Abs
Decor (yr)
Freq (yr)
Idx
Re
Im
Abs
Decor (yr)
Freq (yr)
6
0.92
0.00
0.92
1.01
0,00
13
0.60
0.11
0.61
0.17
3.01
4
0.88
0.09
0.89
0.70
5.23
15
0.58
0.13
0.59
0.16
2.37
0.88
-0.09
0.89
0.70
-5.23
16
0.58
-0.13
0.59
-2.37
2
0.76
0.40
0.86
0.54
1,09
17
0.58
0.03
0.58
0.15
9.78
3
0.76
-0.40
0.86
0.54
-1.09
18
0.58
-0.03
0.58
0.15
-9.78
7
0.85
0.00
0.85
0.50
0.00
19
0.52
0.12
0.53
0.13
2.34
8
0.79
0.00
0.79
0.35
0.00
20
0.52
-0.12
0.53
0.13
-2.34
0.41
0.67
0.78
0.34
21
0.47
0.15
0.49
0.12
1.68
1
0.41
-0.67
0.78
0.34
-0,52
22
0.47
-0.15
0.49
0.12
-1.68
9
0.74
0.12
0.75
0.29
3.26
23
0.47
0.03
0.47
0.11
9.27
10
0.74
-0.12
0.75
0.29
-3.26
24
0.47
-0.03
0.47
0.11
-9.27
0.71
0.00
0.71
0.24
0.00
31
0.39
0.00
0.39
0.00
12
0.64
0.00
0.64
0.19
0.00
27
0.35
0.16
0.38
0.09
1.19
14
0.60
-0.11
0.61
0.17
-3.01
28
0.35
-0.16
0.38
0.09
-1.19
For the ORAS5 model, the linear layer P maps to a 256-dimensional latent space. Training details were otherwise identical, except 12-month input history was used.
The hyperparameters reported above were selected via grid search; Tab. summarizes the ranges explored.
Comparisons. We further compared our method VAMPNets (Mardt et al. 2018). DPNets (Kostic et al. 2024b), Linear Least Squares (LinL.S), and Kernel Ridge Regression (KRR) with a Gaussian kernel. For the deep-learning methods, we used identical training parameters across models. For the classical approaches applied to the raw inputs, we selected the best model via grid search over regularization strengths
[10-5, (10-5_g) 10 10-2] We also varied the estimator rank in the set (5, 8, 10, 16, 32, 50, 64, 128) assess if low-rank approximations in the raw space could recover the dynamics. As shown Tab. 6 and Fig. our method outperforms both baselines the evaluated tasks.
[10-? 10 -0- and, for KRR, kernel coefficients
26

---

Published as a conference paper at ICLR 2026
Table 3: Hyperparameter ranges explored during grid search for the climate modeling task.
Hyperparameter
Search Range
Latent dimensions
[32, 64, 128, [None, 0.1, 0.2, 0.5]
1024]
Max gradient clipping norm
Normalization of linear layer
[False, True]
Regularization
[0, 1e-5, .... 1e-2]
Simplicial normalization dimensions
[0, 2,
16]
History length
[0, 1, 2, 3, 12]
B.5 ABLATIONS
In first set of ablations, investigated the dependence of self-supervised scheme the encoder's architecture. Specifically, we studied the scaling of the loss function with respect to (i) the latent dimension and (ii) the overall parameter count of the encoder.
Scaling laws: Graph-NN encoder. We retrained the SchNet architecture Schlitt et al. (011)) 7) on the data from the protein folding experiment Sec. [.T for three different sizes of the encoder, summarized in Tab. 4 and values of the latent dimension from 4 to 256. The results of this comprehensive ablation study are reported in Fig. I0 We observed monotonically improving losses with respect both an increasing number of training dimensions (panel A) and an increasing model size (panel B). This result provides robust confirmation of the good scalability properties of the loss function (8) studied in this work. As a test-time metric, we evaluated the eigenvalue residuals, defined in (Colbrook et al., 2023, Algorithm 1), see panel of Fig. I0 This metric assesses the extent to which the eigenvalues obtained from our model satisfy the eigenvalue equation Eg Àg. The leading eigenvalue λ is the one enjoying the overall best approximation. Larger architectures are associated with smaller residuals across all the leading eigenvalues.
Table 4: Architectural configuration of the three SchNet model sizes used in the ablation study.
Model
Layers
Filters
Hidden Channels
Params
SchNet-S
2
16
32
6,480 33,088
SchNet-M
3
32
64
SchNet-L
3
64
128
125,536
B
C
Training Loss (SchNet-M)
Loss Sca ing Law
Eigenvalue Residuals (d=128)
SchNet-S
d=4
SchNat-M M
dma
SchNet
-6
d6
dm32
10
SchNet 5 SchNet-L
-10
-12
d 256
d= 128
10k
20k
16 Latent Dimension
128
A2
Training Steps
Figure 10: Scaling laws for the protein folding experiment in 4.1 A Training loss dynamics function the number latent dimensions B Final training loss for three different model sizes, function of the number latent dimensions. C Eigenvalue residuals (lower is better), defined in (Colbrook et al. 2023 Algorithm I) for three different model sizes.
Scaling laws: CNN encoder. The same set of ablations for the climate experiment Sec. 4 with convolutional NN encoder, are reported in Fig. II The overall qualitative behavior exactly matches what already observed for the Graph NN encoder: increasing latent dimensions and/or model size
27

---

Published as a conference paper at ICLR 2026
(see Tab. 5) is associated with higher performance. To rule out the possibility that these improvements are linked to overfitting, in Fig. [1 we report the validation loss, instead of the training loss of Fig. TO
Obtaining the same qualitative results across such distinct physical domains provides strong empirical evidence for the generality of the self-supervised method we propose.
Table 5: Architectural configuration of the three CNN model sizes used in the ablation study.
Model
Layers
Hidden Channels
Params
CNN-S
[8, 16, 24, 32]
12,888
CNN-M
[16, 32, 64, 128]
101,760 397,024
CNN-L
4
[32, 64, 128, 256]
A
B
Va idation loss (CNN medium)
Lass Scaling Law (Validation)
Eigenvalue Residuals (d=128)
0.35
CNN-small
-5
16
Residuels
CNN-meddum
dmm32 32
0.25
CNN-
w0p.
-10
0.20
-15
d 128
0.15
-20
-20
CNN-small
0.10
CNN medium
-25
256
CNN large
0.00
0
10k
20k
8
16
32
128
Ag
7/2
Training Steps
Latent Dimersion
Figure 11: Scaling laws for the climate experiment in of the number of latent dimensions Final validation loss for three different model sizes, as a function of the number of latent dimensions. C Eigenvalue residuals (lower is better), defined in (Colbrook et al. 2023 Algorithm for three different model sizes.
4.3 A Validation loss dynamics as a function
Online versus offline covariances. We conducted an ablation study assess the effect of using covariances Cx. Cxy computed either online during training via EMA or offline from the full training set when estimating the evolution operator E.p.
In the Lorenz '63 experiment, we trained the models as in the main Lorenz-63 experiment (see Ap- pendix B.1) except for lag time set to enable direct comparison between covariance estimation methods. The results show that online covariances yielded better performance, with RMSE and MAE of 0.51 + 0.11 and 0.30 0.06, respectively, compared to 0.63 + 0.21 and 0.45 + 0.19 for offline covariances.
In the climate experiment, the ENSO mode is easily recovered with both approaches. Specifically, for the model trained on ORAS5, the Pearson correlation between the right eigenfunction of E,p and the ONI was 0.72 with online covariances and 0.71 with offline covariances, indicating comparable performance. The associated eigenvalues were also very similar: λENSO 0.9531 + 0. 0.1206i (online) and λENSO = 0.9927 0. 277i (offline).
Stability of EMA covariance. To assess how AMAA-ased covariances converge toward their offline counterparts, computed via a full-pass over the entire training set, we measured their discrepancy in terms of Frobenius norm, i.e., ||CEMa Ciull-pass F. during training the Lorenz '63 data. As shown Fig. 12, this difference peaks in correspondence with the step-like drop in the validation loss, which interpret the encoder discovering number of epochs, as the network converges and settles into a stable representation, the discrepancy steadily decreases and approaches zero. These observations demonstrate how EMA offers robust and practical online approximation the offline, full-pass covariance, offering a clear advantage when dealing with large-scale datasets where computing full-pass covariances may be computationally infeasible.
representational directions. For sufficiently large
28

---

Published as a conference paper at ICLR 2026
C[ -Xrut- lF
10a
P
10-1
w
-4
10-2
-6
1k
2k
Steps
Figure 12: Stability of EMA-based covariance during training on Lorenz '63. The blue curve (left y-axis) shows the discrepancy between EMA and full-pass covariances, while the red curve (right y-axis) shows the validation loss.
Table 6: Performance comparison in terms of Pearson correlation between the right eigenfunction associated with the ENSO mode and ONI, alongside the time per training epoch. Best results are highlighted in bold.
Transfer learning task (model trained CESM, evaluated on ORAS5).
Ours
VAMPNets
DPNets
LinLS
KRR
Pearson correlation (r) Time per epoch (s)
0.81
0.56
0.77
N/A
N/A
25.27 0.74
28.44 + 0.79 29.17 - 0.79
N/A
N/A
Model trained and evaluated on ORASS data.
Ours
VAMPNets
DPNets
LinLS
KRR
Pearson correlation (r)
0.72
0.56
0.62 2.05 = 0.15
0.60
0.63
Time per epoch (s)
1.91 + 0.17
2.03 + 0.15
N/A
N/A
29

---

2010 2015 2020 Published as a conference paper at ICLR 2026
VENSO
00
VAMPNets
20
2
DPNets
500
2
1980
1985
1990 1995 2000
2005 Time
Figure 13: Comparison of ENSO modes retrieved using transfer learning by our method, VAMPNets. and DPNets.
30


Reviews:
[{"review_id": "pK0t8TheE1", "reviewer": "Reviewer_NFuM", "summary": "This paper presents a self-supervised approach to learn the evolution operators, which characterize the temporal dynamics of complex stochastic/deterministic systems. Different from the forecasting or reconstruction type models, the method uses a contrastive self-supervised loss from the spectral contrastive learning to estimate the transfer operators and spectral decomposition directly from data.  The main contribution is to theoretically show the equivalence between the self-supervised loss and minimizing the operator regression objective, connecting the least-square estimator, VAMP-2 score, and HS norms. Experiments on three domains, including protein folding, molecular binding, and climate dynamics, show the effectiveness and generalization of the proposed method.", "questions": "1. It would be a benefit to elaborate more on Equation 8. How can this be viewed as a contrastive learning paradigm? Should state more if positive and negative pairs are not explicitly defined.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. Theoretically explain the operator learning theory and self-supervised contrastive learning, which provides justification using contrastive objectives in scientific dynamic systems. \n2. The method avoids computationally expensive matrix inversions and is implemented with simple matrix multiplications and covariance updates, which are very suitable for large-scale high-dimensional data. \n3. The paper demonstrates the effectiveness through diverse and convincing experimental validation.", "weaknesses": "1. Although the qualitative results, such as the eigenfunction visualizations, look good. The paper does not report any standardized quantitative metrics and comparisons with the baseline on some experiments. \n2. This paper is missing an ablation study that compares different encoder architectures or embedding dimensionalities. \n3. Although admitted by authors that the Hilbert-Schmidt assumption is often not held in deterministic systems. The paper lacks a discussion of how the performance will degrade when the assumption fails. ( Can be one of your ablation studies.)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "rRtI8Sqokc", "reviewer": "Reviewer_FVgy", "summary": "This paper introduces an encoder-only method for learning evolution operators of nonlinear dynamical systems by connecting self-supervised learning with operator theory. The approach is demonstrated across multiple domains, including protein folding, molecular binding, and climate pattern discovery.", "questions": "1. How does your method handle cases where the target function *f* is not in the encoder's linear span?\n\n2. What is the key advantage of your method over existing Koopman or DMD-based approaches, and what new problem does it enable you to solve?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The theoretical analysis provides a thorough and rigorous discussion of the paper's core claims.\n\n2. The experimental design is strengthened by the use of multi-disciplinary datasets, ensuring the scenarios are both diverse and representative.\n\n3. The integration of self-supervised representation learning with evolution operator theory represents a novel and promising research direction.", "weaknesses": "1. When the target function *f* does not lie in the linear span of the encoder, it is unclear how the method should be adjusted. The paper lacks discussion on this point—for instance, whether simply increasing the embedding dimension would suffice to satisfy this assumption.\n\n2. The paper lacks practical guidance on how to choose the embedding dimension in real-world systems.\n\n3. The authors should more clearly distinguish their approach from existing methods that combine deep learning with Koopman theory or DMD, and better highlight what additional problems their method can solve.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "6F0WrYmkt5", "reviewer": "Reviewer_3ogW", "summary": "This paper proposes an end-to-end, self-supervised approach to learning evolution operators for high-dimensional nonlinear dynamical systems, with a particular focus on scientific domains such as molecular dynamics and climate modeling. The authors connect contrastive, self-supervised representation learning techniques with the operator-theoretic framework, providing theoretical grounds for their approach and demonstrating its equivalence to established least-squares operator learning methods. Comprehensive experiments are reported on protein folding, molecular binding, and climate data, with claims of scalability, interpretability via spectral decomposition, and transferability of learned representations.", "questions": "- Can the authors provide (or reference) quantitative metrics for spectral decomposition accuracy on real-world, high-dimensional datasets, or propose benchmarks for such evaluation? For example, can clustering purity or precision/recall be reported for folded/unfolded distinctions in protein folding, or event detection in climate data?\n- Can the authors report on representation transfer performance when the target system is more dissimilar to the source (e.g., transfer across molecular families, not just ligands, or different climate regimes)? What are the limits of transferability in practice?\n- In extremely large-scale systems, does online or EMA-based covariance computation pose bottlenecks or stability issues? Can sparse or approximate methods be used safely?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- The authors make a deliberate connection between self-supervised contrastive learning and classical operator-theoretic approaches (notably, the least-squares estimation of evolution operators). This is not only highlighted at an intuitive level but also underpinned with theoretical results.\n- The paper supports claims with diverse, high-dimensional benchmarks, spanning molecular simulations (protein folding, ligand binding) and challenging climate datasets.\n- The spectral decomposition of the learned operator yields interpretable, physically relevant modes, e.g., hydrogen bonding patterns in protein folding.", "weaknesses": "- While the core mathematical exposition is sound, there are places where notation could be clarified for easier accessibility:\n  - In Lemma 2 and Appendix A, the formula for the predictor $P_*$ involves both $C_X^{-1}$ and $C_Y^{-1}$, but the mapping from the loss gradient to the closed-form $P_*$ is sketched rather than fully elucidated.\n  - Equation (9) relates the action of E to the covariance of futures, but implementation choices for non-stationary or non-ergodic systems are not fully discussed.\n- The limitations section (Conclusion) briefly notes the Hilbert-Schmidt assumption and qualitative evaluation. However, broader issues are not meaningfully discussed. For example, the treatment of deterministic versus stochastic systems, failure cases when the operator is not Hilbert-Schmidt, and scalability bottlenecks in training or memory for extremely large state spaces are not addressed.\n- While the online/offline covariance ablation adds value, ablations for the choice of network architectures (GNN vs. CNN vs. MLP), size/sparsity of $P$, or impact of history/context window on model performance would strengthen the empirical narrative.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "lpPBm2tfYW", "reviewer": "Reviewer_eZUJ", "summary": "This paper presents an end-to-end deep learning framework for learning evolution operators of high-dimensional dynamical systems. It establishes a novel connection between self-supervised contrastive learning and operator learning theory.\n\nThrough experiments, the paper demonstrates that the proposed method can extract meaningful dynamical modes that reflect physically relevant processes such as folding-to-unfolding transition, molecular binding, and climate oscillation. It shows that these learned representations can be reused effectively across related systems.", "questions": "(1) Your Lemmas assume the evolution operator $E$ is Hilbert-Schmidt, yet you acknowledge this is often violated by deterministic dynamical systems. The Lorenz '63 system in your experiments is deterministic (or near-deterministic). So, what happens to the approximation quality when $E$ is not Hilbert-Schmidt?\n(2) In the deterministic case, the Koopman operator acts on delta functions. How does your bilinear model $\\langle \\phi(x_t), P\\phi(x_{t+1}) \\rangle$ approximate this structure?\n(3) How does the approximation error related to the latent dimension?\n(4) In the climate experiment, how does your method perform if one applies a kernel-based approach directly to the raw 29,040-dimensional input features, instead of using the learned 128-dimensional embedding? Such a comparison would clarify if standard dimensionality reduction techniques could achieve similar results.\n(5) Could you discuss any scenarios or system characteristics under which the learned representations may fail to transfer effectively?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper established a novel and rigorous theoretical connection between evolution operator learning and self-supervised contrastive learning. A major contribution is the scalability of the proposed framework to high-dimensional system. Empirically, the model demonstrates stable training and meaningful spectral decomposition in different scientific domians, even in terabyte-scale dataset.", "weaknesses": "(1) The finite-sample analysis is not provided;\n(2) The paper lacks a clear comparison against other (Koopman) operator-learning approaches (e.g., kernel EDMD, Neural Operator frameworks).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)