Efficient Parametric SVD of Koopman Operator for Stochastic Dynamical Systems
Authors: Minchan Jeong
Jongha Jon Ryu
Se-Young Yun
Gregory W. Wornell
Venue: NeurIPS 2025 poster ()
Keywords: Koopman operator
spectral decomposition
singular value decomposition
neural networks
low-rank approximation
DMD
VAMP
DPNets
Primary area: machine_learning_for_sciences

Abstract: The Koopman operator provides a principled framework for analyzing nonlinear dynamical systems through linear operator theory. Recent advances in dynamic mode decomposition (DMD) have shown that trajectory data can be used to identify dominant modes of a system in a data-driven manner. Building on this idea, deep learning methods such as VAMPnet and DPNet have been proposed to learn the leading singular subspaces of the Koopman operator. 
However, these methods require backpropagation through potentially numerically unstable operations on empirical second moment matrices, such as singular value decomposition and matrix inversion, during objective computation, which can introduce biased gradient estimates and hinder scalability to large systems.
In this work, we propose a scalable and conceptually simple method for learning the top-$k$ singular functions of the Koopman operator for stochastic dynamical systems based on the idea of low-rank approximation. Our approach eliminates the need for unstable linear-algebraic operations and integrates easily into modern deep learning pipelines. Empirical results demonstrate that the learned singular subspaces are both reliable and effective for downstream tasks such as eigen-analysis and multi-step prediction.

Efficient Parametric SVD of Koopman Operator for Stochastic Dynamical Systems
Minchan Jeong'; Jon Ryu²", Se-Young Yun', Gregory W. Wornell
Kim Jaechul Graduate School of AI, KAIST, Daejeon 34141, South Korea DDeparment of EECS, MIT, Cambridge, MA 02139, United States yunseoonnyy kr, {jonhha gww edu
Abstract
The Koopman operator provides principled framework for analyzing nonlinear dynamical systems through linear operator theory. Recent advances in dynamic mode decomposition (DMD) have shown that trajectory data can be used to identify dominant modes of system in a data-driven manner. Building on this idea, deep learning methods such as VAMPnet and DPNet have been proposed to learn the leading singular subspaces of the Koopman operator. However, these methods require backpropagation through potentially numerically unstable operations on empirical second moment matrices, such as singular value decomposition and matrix inversion, during objective computation, which introduce biased gradi- ent estimates and hinder scalability to large systems. In this work, we propose a scalable and conceptually simple method for learning the top-k singular functions of the Koopman operator for stochastic dynamical systems based on the idea of low- rank approximation. Our approach eliminates the need for unstable linear-algebraic operations and integrates easily into modern deep learning pipelines. Empirical re- sults demonstrate that the learned singular subspaces are both reliable and effective for downstream tasks such as agge-annalssis and multi-step prediction.
Introduction
The Koopman operator theory offers a powerful framework for analyzing nonlinear dynamical systems by lifting them into infinit-diimensional function space, where spectral techniques from linear operator theory can be applied. Recent advances dynamic mode decomposition (DMD) have shown that trajectory data can be effectively used to identify dominant eigen- mogees in
data-driven manner [40, 46 47 18 4]. Inspired by the success of deep learning, recent methods such as VAMPnet [8 24 and DPNet [66] employ neural networks to approximate the leading singular subspaces of the Koopman operator. While shown effective for some benchmark tasks, these approaches often rely numerically unstable operations such singular value decomposition (SVD) or matrix inversion during objective computation. These operations present practical challenges, particularly in computing unbiased gradients and thus scaling high-dimensional systems.
In this work, we propose a conceptually simple and scalable method for learning the top-k singular functions of the Koopman operator for stochastic dynamical systems. Our approach builds on the idea of low-rank approximation, which has recently received increasing attention in the literature due to its favorable optimization structure that aligns well with modern optimization practices; see,
I7 in the machine learning literature. Our method avoids unstable linear-algebraic computations and is easy integrate into modern deep learning pipelines. We demonstrate that it reliably recovers dominant Koopman eigen-subspaces and supports downstream tasks such as prediction and eigen-analysis.
e.g., 21 45 in the numerical optimization literature, and 42 39 32
'The two authors contributed equally to this work.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

2 Problem Setting and Preliminaries
In this section, we introduce the problem setting and establish the foundation for the subsequent discussion. We begin by formulating the problem in the context of discrete-time dynamical systems, which will serve as our primary focus. We then briefly address the continuous-timg timmo case. Next, we review two existing approaches, VAMPnet and DPNet, and highlight their limitations, thereby motivating the need for our proposed method.
Notation. We denote linear operators by stylized script letters, e.g., U, X, L, and I, using the Zapf Chancery font to distinguish them from matrices and scalar functions. Bold lowercase letters such as x and y are reserved for vectors, while bold uppercase letters like X denote vector-valued random variables. Sans-serif uppercase letters, e.g., M and S, are used for matrices. In particular, I denotes the identity matrix.
2.1 Discrete- Dynamical Systems
Consider a stochastic discrete-time dynamical system Xc+1 = E(F((x), Here, F: X' - X' is a possibly nonlinear mapping for a domain X' S Rd. 6 ~ D is an independent random noise X captures independent noise in the process, such as the
variable, and X; + (x)),)) additive white Gaussian noise Assuming that F and the noise distribution D are time-invariant, the process becomes a time-homogeneone
Markov process with transition density p(x' [x), i.e., Pr(X,++ E A = x) = J' p(x') x) dx' for all measurable sets A C X and all 2 0.
In the stochastic setup, the dynamics is fully captured by the transition density p(x' xx), which is induced by (f o (-). and thus the problem becomes analyzing p(x' x) of Markov chain. Here, the Koopman operator becomes the conditional expectation operator, i.e., for an observable g: X + R.
(xg)(x) _ Epxx x) (x)))x))
By the Markov property, repeatedly applying the Koopman operator will correspond to the multi- step prediction in terms of the posterior mean of gX,) given Xo = Xo. i.e., (xx)
Ep(x.| [x())] Throughout the paper, we will assume that K is compact following [16). which is mild assumption that holds for a large class of Markov processes; sec, c.g., [14]. We remark that we cannot directly apply the existing spectral techniques to a deterministic dynamical system, including the technique developed in this paper, since the corresponding Koopman operator is not compact; see [48] Appendix A.5J. The consideration of stochasticity breaks the degeneracy and allows the linear algebraic tech- niques to be applicable. Moreover, assuming stochasticity is not necessarily a restrictive assumption, as physical processes in the real world may be inherently noisy.
Problem Setting, Our goal is to analyze the Markov dynamics of the length-T trajectory (Xo, x Xx) ~ p(xo)p(xi xx)... ..p(TT] xT-1), where p(Xp) is distribution for the ini- tial state. In practice, we assume that we have access to N independent, random trajectories, and collect the consecutive transition pairs {(x0, x1), (x1, x......, (XT-1,XT)} and define the joint distribution over the pair as p(x, x') & + E Pe (x))(((x') x). where Po(x) p(xu) and Pt+1 Px+1(x') e f'p(' (x)xx We let (x() and p (x') denote the marginal distribution over X' of the current and future states. Note that while pi (x') = Epp())]((x' x)] holds, the two distributions, Pu(x) and p1 (x') can differ significantly, particularly when the trajectory length T is relatively short.
Let FX. be a measurable space and let Po and PI be (finite) measures on (X,F) repre- senting the distributions of the current and future states, respectively. For any measure p on F) define 123(2) A {f: X + CI p(xx) < 2o). Equipped with the inner prod- uet (J.g)e A S ((x))(x))((xx), 122(2) is a Hilbert space. The Koopman operator K. then a mapping from L2 (X) to L2 (2)), which is an integral kernel operator with the transition kernel x(x,x') A p(x
mx(x") The adjoint operator X then acts as the backward predictor, i.e., (x* ()()) [x(x|x))[(xx)) where q(x x') A polx)p(x' [x) denotes the conditional distribution pi(x')
induced by Pol (x))((x' xx). Similar to K, repeated application of x yields multi-step backward prediction:
Ep(xo x.) [f (xo)).
2 All machinery developed here can also be adapted for discrete state space, but we focus continuous X.

---

Special Case 1: Stationary Markov Processes, Given p(x'1 x), let π (xx) denote the stationary distribution, i.e., the distribution satisfies (x(x) (p(x') x)] = Txx).. Under mild regularity conditions such as ergodicity and irreducibility, such distribution uniquely exists F In the stationary case, the marginal distributions of current and future states equal, i.e., Pu P1 = n,, and the Koopman operator X becomes a map from L (2() to (2)). Assuming ergodicity, we can collect the time- lagged pairs (x, from long, single trajectory, time averages converge to expectations under the stationary distribution.
Special Case 2: Reversible Markov Processes. A Markov process is time-reversible if and only if it satisfies the detailed balance condition: the joint function P(x,x') x') _ m(x)p(x') x) is symmetric, ie.. P(x, x') P(x',x x) for any x, x' X x X. In this case, the Koopman operator becomes self-adjoint, and thus the eigenfunctions and eigenvalues are real. This is much stronger condition than the normality of Koopman operator and thus the easiest to deal with. For the case of reversible processes, Noé and Nüske (28) proposed a method to approximate eigenfunctions from time-series data, followed by the extended DMD (EDMD) [46].
2.2 Continuous- Time Dynamical Systems
Let (x1)120 be home-hooogeneous Markov process (e.g., Langevin dynamics). The Koopman semigroup {u') |120 acts on observables f: X + as (2t' f)x) & Eptx, xx)) (x)]]( This defines a strongly continuous semigroup satisfying et I (identity operator) and ut' e1 The (infinitesimal) generator L of the Koopman semigroup is given by:
(ut f)(x) f(x)
d
(x))
A lim
dt Ep(x xxx (x))).
where the limit is taken in the strong operator topology.
In general, generators may not be bounded, and thus not compact. Similar to [48] however, one can view that this discrete-time dynamics is discretized version of continuous dynamics with lag time
> O, and the Koopman operator of the discretized system is often compact; see [48]. Moreover, as argued in Kostic al. [16), we can also directly apply the developed technique for special, yet
important continuous-ting dynamical systems such as an overdamped Langevin dynamics. Unless stated explicitly (like in Section 3.1), we will describe our techniques for discrete-time processes.
2.3 Data-Driven Learning of Singular Subspaces: VA.PPnt and DPNet
Learning the dominant singular subspaces of the Koopman operator is key goal for understanding complex dynamical systems, especially in data-driven modeling. First, for non-normal operators, which commonly arise in irreversible or non-equilibrium dynamics, the singular value decomposition (SVD) provides the best possible low-rank approximation measured by Hilbert-Schmidt norm. Second, if process is reversible and stationary (i.e., (x) pi (x) = т (x)), the Koopman operator self-adjoint. In this case, its SVD coincides with its eigenvalue decomposition (EVD). Third, even for irreversible processes, the dominant singular functions still capture essential dynamical features. For example, they can show kinetic distances between states, much like eigenfunctions do in reversible systems [BI, or also find coherent sets in changing Markov processes; these are generalized forms of long- oag-lasting states 13 See also 48 Section 2.3].
These theoretical advantages have motivated recent developments in neural network- beasork- approaches, such as VAMPnet [48] 24 and DPNet [66]. which aim to learn the top-k singular subspace of the Koopman operator directly from trajectory data. Both methods are grounded in variational characterizations of the dominant singular subspaces, yet they differ significantly in their training objectives and optimization strategies. In the following, we briefly contrast these training approaches and highlight their common limitations. Here we focus on the training approaches, and defer the discussion on their inference procedures to Section 3.2) Also, unlike the score convention commonly used in the literature, we will follow the loss convention throughout, whereby the training objective is formulated quantity to be minimized.
Forr example, if the stochasticity in the system is induced by an additive noise, i.c., X = F(x:) + ç for ~ D, with the density of D being positive almost everywhere, then the system is asymptotically stable, i.e., it converges to unique stationary distribution; see [19 Corollary 10.5.1].
3

---

2.3.1 Common Setup
Suppose that we wish capture the top-k singular subspaces using neural networks x + f(x) - [fr (x)..... ................................ ft (x)]T R and x' + g(x') - [ (x')..... 9 (x')T I² These are sometimes referred to the encoder and the lagged encoder, respectively, which are intended to capture the left and right singular subspaces. Since a Koopman operator, as a special example of canonical dependence kernels, always has the top singular functions Φ1 (x) = 1 and 1/1 (x') = 1 with singular value I, we can simply set (x) - and gi (x') < [32] We refer to this operation as centering, as it ensures that the remaining modes are centered, i.e., orthogonal to the constant modes. We adopt this parameterization by default, although many existing methods do not follow this practice. Properly handling the constant modes is crucial, as mishandling them can lead to misleading results; see Appendix A.2 for further discussion.
In variational optimization framework for SVD, the key quantities are second-momen matrices. We thus introduce the following shorthand notation. The second-moment matrix of vector-valued functions h1 and (-)-) with respect to distribution is denoted as
M, [h,, h2] A Bp(x) (x))2(xx))
and we write M, hi & M, h| for shorthand. The joint second-moment moment matrix of h h,()) and h2 over the joint distribution po (x))(x)((())))))(((()(() x) is defined and denoted as
T[h bo] (x()))p(x) x<) [b (x)
which satisfies the identities T( h2] = p, h1, Kh2] Mp, KKhh,hy h2)-
In what follows, will define population objective function using the population second moment matrices Mp [f], Mp (g). and g]. In practice, these quantities are to be estimated with trajectories,
and their empirical estimates are denoted with hats ("), i..., Mpo [f]. M [g). and T(f.g).
2.3.2 VAMPnet
For an integer r 2 1, Wu and Noé 48 introduce the VAMP-r objective
Cuamp (f.g) A (M ff)- 3 T (M, (g)
(I)
Here, || ||~ denotes the Schatten r-norm. The variational principle behind the VAMP-r objective is
explained in Appendix B.1 for completeness. In practice with finite samples, to avoid the numerical
instability in computing the inverse matrices, CCamp- (f. g) - (M [f] + XI) 1 T g (Mp, [g]
Al) 1 fr. Tuning in practice can be done by a cross validation. In the literature, the use of 24
48 has been advocated.
2.3.3 DPNet
Kostic et al. [6] proposed an alternative objective, called the deep projection (DP) objective,
<l (f.g)
-| ( (g) ! [| 7(R(M po (f) + RMM ig))),
(2)
where they further introduced the metric distortion loss R: 84 + R defined as R(M) A tr(M²
M In M) for M 0. Note that with O. the objective becomes equivalent to the VAMP. objective. To detour the potential numerical instability of the first term, which the VAMP-2 objective, the authors further proposed relaxed objective called the DP-relaxed objective
gll
der (f.g)
+ (R p. [f) + RMM (g))).
(3)
||M pu Ilo Mp. oglloop
where || |lop denotes the operator norm. In both cases, Kostic et al. [16 argued that using > is crucial for improving the quality of the learned subspaces. In contrast, the scheme we propose below does not require such regularization, and we empirically found that it offers no benefit.
*We note that this expression corresponds to the maximal VAMP-r score in the original paper 48 but we call this VAMP-r objective, which is slight abuse, as this is the training objective train neural network basis.

---

2.4 Practical Limitations of VAMPnet and DPNet
Although these objectives are well-founded in the population limit for characterizing the desired singular subspaces, they do not permit efficient optimization via modern mini-batch-base training.
To compute the VAMP-1 objective [24] one must evaluate the matrix square root inverse Mo [f]
and M [g), as well the nuclear norm of the matrix (Mpu [/2)-_//T[f,,
(g))- These
computations involve numerical linear algebra operations such as eigenvalue decomposition and singular value decomposition of empirical second-momen matrices. Such inverses may be ill-defined or numerically unstable when span(f) is nearly rank-deficient during optimization, and further, introduce instability during training
backpropagating through these numerical operations Moreover, since empirical second-moment matrices are estimated from mini-batch samples, the resulting gradients can be highly biased, potentially slowing convergence during optimization.
The VAMP-2 objective [48] suffers from similar issues, as it also requires computation of matrix inside the Frobenius norm. The DPNet objective in Eq. (2) introduces the additional metric distortion loss (M), which inherits both of the aforementioned issues. Similarly, the operator norms in the denominator of the DeNet-relaxed objective in Eq. 3 are subject to the same challenges.
3 Proposed Methods
In this section, we introduce an optimization framework based on the idea of low-rank approximation (LoRA), which circumvents the issues in the existing proposals. Although LoRA has been recently explored in various contexts (1 singular functions of the Koopman operator remains largely unexplored. A detailed discussion of related formulations and prior work on LoRA provided in Appendix [ After describing the learning procedure, we investigate two inference methods that build upon the learned representations.
177, its application to learning parametric
3.1 Learning
To bypass the numerical and optimization challenges, we propose directly minimize the low-rank approximation error || - fa @ gi fillis to find the Koopman singular functions. Suceinetly, the learning objective, whose derivation is provided in Appendix C.T. can be expressed as
Liora (f, g) - 2 tr.(TTf..f....
g))
tr Mf [g))
(4)
The celebrated Eckart-Young-Mrsky theorem 6 27 (or, more precisely, Schmidt's theorem 3d
establishes that this objective precisely characterizes the singular subspaces of K as a global optimum.
Proposition 3.1 (Optimality of LoRA loss; see, e.g., [32 Theorem 3.1]). Let X: L (X) - L (X) be a compact operator having SVD E oe @ ei with 2 global minimizer of (f, g). If o > okti, then E f gi Σ o. @ vh.
02 2 2 0. Let (f", (f",g") denote
Moreover, owing to the simple form of the objective, its lcarning-theood theoretic properties are amenable to analysis. In particular, under a mild boundedness assumption, the empirical objective Liora (f. converges to the population objective (f, g) at the rate Op(N 1/2), when N denotes the number of pairs from the trajectory data. We defer the statement of Theorem C.T to Appendix C.2
notable property of the LoRA objective, compared to the VAMPnet and DPNet objectives, that it is expressed entirely as a polynomial in the second-moment matrices. As result, its gradient can be naturally estimated an unbiased manner, in contrast to VAMPnet and DPNet, which require backpropagation through numerical linear-algebra operations. This property is particularly advantageous when optimizing large-scale models with moderately sized minibatches.
Special Case: Reversible Continuous- Time Dynamics. As argued in [16], we can also directly analyze reversible continuous- time dynamics, which includes an important example of (overdamped)
'In the PyTorch implementation, functions such lstsq, eigh, and matrix norm from the torch. linalg package are used; see, e.g., the PyTorch implementations of VAMPnet in indeettime and kkooplearn libraries. For example, cigendecomposition (torch. linalg, aaalg..iigh eigh) may be numerically unstable when the matrices are ill- conditioned. since the gradients of eigenvectors scale inversely with eigenvalue gaps, which can lead to gradient explosions in nearly degenerate subspaces.
5

---

Langevin dynamics. In this special case, we can apply the spectral techniques under a weaker assumption than compactness; for example, we only require the largest eigenvalue to be separated from its essential spectrum; see, e.g., [12 Section III.4). Our objective in Eq. (4) simplifies to
Csa lora (f) 2 tr pr [f, Lf) + || [f]]
(5)
where Ma [f, Lf] A (C(( is plugged in in place of T[r, gj M [f, Ag). Ryu al. [2] Theorem C.S] show that the LoRA objective applied on a (possibly non-compact) self-adjoint operator can find the positive eigenvalues that are above its essential spectrum. We describe a special
example of stochastic differential equations in Appendix G.3
Nesting Technique for Learning Ordered Singular Functions. As introduced in [32 we can apply the nesting technique to directly learn the ordered singular functions. We note that learning the ordered singular functions is not an essential procedure, given that we only require well-learned singular subspaces during inference, as we explain below. We empirically found, however, that LoRA with nesting consistently improves overall downstream task performance. We conjecture that the nesting technique helps the parametric models to focus on the most important signals, and thus improves the overall convergence.
The key idea of nesting is to solve the LoRA problem for all dimensions i E [k] simultaneouslyy. There are two versions proposed in 32inttont and sequential nesting. On one hand, in joint nesting, we simply aim minimize single objective Cloin g) _ Ek= ai (f g1:;) for any choice of positive weights ..,,..., 0r > 0, while the uniform weighting has been proven to be the most effective [32]. The joint objective characterizes the ordered singular functions as its unique global optima. On the other hand, sequential nesting iteratively update the i-th function pair (fi, gi), using their gradient from Ciora (f....gis)). as if the previous modes (f13-1.g1(-11) were perfectly fitted to the top-(i - 1) singular-subspacee Given that different modes are independently parameterized, one can use an inductive argument to show the convergence of sequential nesting. Both joint and sequential nesting can be implemented efficiently, with almost no additional computational cost compared to the LoRA objective without nesting. We defer the details to Appendix [F]
Ryu et al. [32 advocate using sequential nesting for separately parameterized networks, while recommending joint nesting otherwise. In all our experiments, however, we employed two neural networks lg and go (i.e., assuming joint parameterization). and empirically observed that sequential nesting achieved convergence comparable to joint nesting. Thus, a practitioner may by default adopt joint nesting as a principled choice for jointly parameterized networks, while regarding sequential nesting as another option as a working heuristic. In the molecular dynamics experiment below, we found that the nesting technique can provide a stabilization effect during training; see Section 4 and Appendix G.4
Explicit Parameterization of Singular Values. An alternative, yet natural parameterization is to explicitly parameterize the singular values by learnable parameters Y E R& and plug in f + of and g & v to the LoRA objective in Eq. under the unit-norm constraints Il |1 = gexplicit (7, f,g) 2
1 for all E [k:. Then, we get the explicit LoRA objective Lora
-2 v
gl)++r(M 70 ffMM (gg). In a similar spirit to [16], Kostic al. 17 proposed a regularized
[ 11 + [f]l and objective pexplicit (y, f, g) + A(R [f] R2,, [g]) where R2 [f] A |
R/ [g] is similarly defined. An alternative implementation of the explicit parameterization involves La-batch normalization, as suggested by Deng et al. [3]. We experimented both formulations and empirically observed that the original parameterization in Eq. (4) performs well in practice, and no performance improvement was observed with these regularization approaches.
3.2 Inference
After fitting f(-) and g(-) to the top-k singular subspaces, we can perform the downstream tasks such as (1) finding ordered singular functions, (2) performing anallsss, and (3) multi-step prediction. We describe two approaches, each of which is a slightly extended version from VAMPnet and DPNet, respectively. The two approaches are based on rather different principles: Approach (VAMPnet type) is based on estimating the Koopman operator by LoRA using both left and right singular functions, while Approach 2 (DPNet type) uses one set of basis and performs downstream tasks by projecting the Koopman operator onto the span of the basis. In the experiments below, we will compare and

---

Table Summary of inference procedures given learned singular functions (or basis) f. g.
Approach 1. CCA + LoRA
Approach EDMD (b E {f,
f(xx _ (Mpa (x)
CCA step Whitening
N/A
g(x') _ (M,r (g) (g)))>((((xx
CCA step SVD
T[f, g) USVT
N/A
((x) - 21(2UT((X)
Ordered singular functions
N/A
(xx) = E'(2)
K & M₂s [yb, <b] (right)
Kg & (M,, [b]) + b (right)
Approximate Koopman matrix
K
& Mou [ab, Φ) (left)
Kg (M₂1 [b]) [b))+TT T[b (left)
0. (xxx)"
Forward Ep(x, [xx))
(Kright
1 (le, 1)
b(xo)T(KK")) (MMgg)((MMM)))) [b]) + (h, b)
2.
Backward Epixo! xx [b(x))]
(xx,) T (K
left.
) (h, Φ)
[b]) (h, b)mr
a.p
discuss their pros and cons. For simplicity, we describe the procedures using population quantities; in practice, they replaced with empirical estimates (i.e., M t- M and T T).
3.2.1 Approach 1: CCA + LoRA
In similar spirit to the inference procedure VAMPnet 488] we can perform canonical correlation analysis (CCA) [10] to retrieve the (ordered) singular values and singular functions follows, given that f(-) and g(-) capture left and right singular subspaces. First, we define the whitened basis
functions f(x) _ (M,o f//ff(x and gxx') & (Mp, [g)) 1/2 gx'). Then, we perform the SVD of the joint second moment matrix T[f. g] = UEVT, where U E Raxr E E Rexr, V Rexr, We define aligned singular functions as
x(x) X1/2UT((X) E R and (x)) / (x)) ER"
(6)
and approximate the transition kernel as k(x, x') plx -(*| ~ =(x)
pi(x')
Hence, compared the direct approximation f(x)' the CCA procedure whitens (by Mp [f] and Mo (g)) and corrects (by SVD of T(f,g) the given basis. In a real-world scenario, the CCA alignment can always help f. g better aligned, and improve the quality of singular value estimation.
Given this finite-rank approximation, the eigenfunctions approximate Koopman matrices, based on the theory developed for finite-rank Koopman operators in Appendix D Specifically, the eigenvalues of K and the right eigenfunctions using the matrix Kright Φ.
be reconstructed using the following
be approximated A Mp,l4b, ø]; see Theorem D.2 Similarly, the left eigenfunctions and
corresponding eigenvalues can be approximated using 3. A Mp [ Φ]; see Theorem D.3 We
note that this eigen-analyse using the finite-rank approximation is new compared to [48]
Lastly, given the LoRA p(x' (x)
(xx)) Tx (xx),, the conditional expectation be approximated as
pi(x))
p(x' xx)
Ep(x (x) (xx'|1x)/i [h(x')]
E Ex))
h(x))
Σ (b,, x(x) (h. 3)p
(x
where we let (h, }),, [(h, pi........, (h, Fx) p. E IR. Extending the reasoning, can obtain an approximate sulti-step prediction as
Epix, xx) [xx)] )] Rv ((xx) ( (Krriett e-1 (h, B)ps 3.
f(x0) TT & (M,, ig. (TT(f.g) g)..
Mp, lg. [f, is diagonalizable, we can perform its EVD to make the matrix power computation
more efficient. Similarly, we can perform the multi-step backward prediction using Kleft as follows:
Ep(xo x.) [x)o)())) (xx)) TT f] T -1)
(h, F)
(8)
7

---

3.2.2 Approach 2: Extended DMD
Once a good subspace span(b) has been learned using some basis function b: X - R Kostic et al. [66 proposed to perform the operator regression [14] (also known as principal component regression (S)). which is essentially the EDMD [46). The EDMD approximates the Koopman operator by kg - (M [bb)+ & which can be understood as the best finiteddimensional approximation of the Koopman operator K. restricted on span(b), in the sense that (*b))(x) = Ep(x (x]] se Kb(x). Now, given a function h span(b),
can again choose the least- (M,,(b) [b)) (h, b) po to find the best z such that h(x) R ==b(x) for (x(x). Given this, we can finally approximate the multi-step prediction as
squares solution z (h) _
Ep(x, xx) [x+]] b(xa) (Kg (h) b(xo) (Tb)
(M [b]) (h, b) (9)
Applying the same logic to the adjoint operator K", we can perform the backward prediction as
Exxx x. [x(x))] b ((Mp) [[))T [b)+ (h, b)m-
(10)
We defer a more detailed derivation of the EDMD predictions to Appendix L We note that Kostic al. [66 proposed to use the left singular basis f for the basis b, whereas we empirically found that using the right singular basis yields comparable results.
4 Experiments
We demonstrate the efficacy of the proposed techniques using the experimental suite of [16]. Unless stated otherwise, all experimental settings are mostly identical those in [16] except the molecular dynamics experiment. All technical details and configurations are provided in Appendix G The appendix also includes an additional experiment on an instance of a ID noisy logistic map, whose Koopman operator has finite rank. We defer this result to the appendix. as most methods perform comparably in this simple setting. Our PyTorch 30 implementation is publicly available at ht //github. www.://aaacaaaannog//o/o/////// Kooppaansaooommn
4.1 Ordered MNIST
We considered the ordered MNIST example, a synthetic experiment which was first considered in [14] given an MNIST image x, with digit yt E (O....,4), x is drawn at random from the MNIST images of digit yt-1 = V/ + (mod 5). While this process is not time-reversible, the process clearly normal. We tested VAMPnet-1, DPNet, DPNet-relaxed, and sequentially and jointly nested version of LoRA. For each method, we trained singular basis parameterized by convolutional neural networks with 100 epochs using 10 random seeds.
We evaluated the multi-step prediction performance of each method, by computing (1) the accuracy using an oracle classifier similar to [16], and (2) root mean squared error (RMSE) of the prediction evaluated on the test data. The multi-step prediction performed with EDMD(g) reported in the first row of Figure T. We highlight that LoRA and its variants consistently outperform the other methods, exhibiting robust RMSE performance over range of prediction steps E {-15...., 15). Notably, the nesting techniques helped improve the R.SE,,, especially the joint nesting worked best. We also remark in passing that, unlike the significant failure reported [66, the VAMPnet-1 prediction quality is comparable to DPNet in both metrics.
In the second row, we showed the performance of different prediction methods with the basis learned with LoRAjnt- We note that the quality is very close to each other, and CCA+LoRA prediction seemss follow the trend EDMD(f) for forward prediction and EDMD(g) for backward prediction. We also remark that the prediction quality with CCA+LoRA when = 1 is particularly bad, which we conjecture to caused by the absence of the action of Koopman matrices; see Eq. (7) and Eq.
4.2 Langevin Dynamies
We also tested the performance of LoRAses with the objective in Eq. 5 to learn the eigenfunctions of ID Langevin dynamics, which a continuous-time, time-reversible process; see Appendix G.3 for the stochastic differential equation used in the experiment.

---

Multistep Accuracy w/ EDMD(g)
Multistep RMSE w/ EDMD)
WAMPnet-1
DPNet
0.975
40
DPNet
RASS
LORA
0.950
LoRAus
dgt
0.925
LoRA,n
Multistap Accuracy (LoRApr)
Multistep RMSE (LoRAm)
1,000
accurcyy
0.975
40
EDMDIA
20
EDMD(g)
OGA+LORA
0.925
20
0.900
-15
-10
10
15
-15
-10
-5
10
time steps
time steps
Figure| 1: Summary of the MNIST experiment. The shaded area indicates #1 standard deviation.
\ = 0.00
A 1.20
7.89
A
15.36
0.1
0.0
elovals
-50
-01
As
82.08
Ag
94.93
A7
132.62
Ao 144.39
100
0.1
0.0
150
0.1
10
107
10"
1
1
training steps
Figure 2: Visualization of the cigenfunctions of the ID Langevin dynamics, learned by LoRAseq In the first panel, the learned eigenfunctions across training iterations overlaid, with later iterations displayed with higher opacity (red). The dashed lines indicate the true eigenfunctions.
We drew a sample trajectory of length 7 x 10s by the Marruamma scheme and used for training. We parameterized the top- 10 eigenfunetions using single MLP with hidden units -28-128-122 and CeLU activation. We trained the network using Adam optimizer with learning rate 10 for 50,000 iterations with batch size 128. The differentiation operation in the generator was computed by the autodiff feature of PoTorch. Exponential moving average with decay 0.995 was applied to result in a smoother result. In Figure 2 we report the first eight eigenfunctions and the convergence behavior of the estimated eigenvalues, without any postprocessing other than sign alignment. As demonstrated, the learned eigenfunctions and eigenvalues converge to the ground truths reliably with LoRAseq-
po [f, <f) + RR(M po (f). following the recommended configuration in [16], but we were unable to obtain successful results.
We also implemented the DPNet objective ((7) (f) A tr(M [f] [f
This example demonstrates the effectiveness of the LoRA framework for continuous-time dynamics.
4.3 Chignolin Molecular Dynamics
To assess scalability on high-dimensional data, we apply our method to the analysis of chignolin molecular dynamics. Chignolin is an artificial mini-protein that is considered standard model for
studying rapid folding dynamics due to its complex and fast transitions 20 B] [6]. The underlying physical system is time-reversible, and thus the Koopman operator is self-adjoint. The slowest decay mode chignolin is associated with the folding-unfolding transition, which occurs on a microsecond timescale 20). Unlike [16], we use the public dataset of [25]. To isolate algorithmic stability from variance due to data sampling, we fix a data split chosen to balance the root mean square distance
statistics across splits. All experimental details are provided in Appendix G.4
We first remark on the overall quality of the learned models. As reported in 16 we observe that VAMPnet-2 and DPNet diverge. However, we empirically find that models trained with VAMPnet-1 can converge, albeit with several caveats. First, we observe that VAMPnet-1 is sensitive low- level numerical implementations in torch. lihhllinnggg..mmmmm and diverges under specific PyTorch version. Moreover, we observe that VAMPnet-1 diverges for certain alternative train-test splits that are not reported here. Consequently, the only methods that consistently exhibit convergence in our experiments are our LoRA variants. We nonetheless report VAMPnet-1 results to demonstrate that method outperforms it even in regimes where VAMPnet-1 converges.

---

Table 2: Test VAMP-E scores on the 300 K chignolin dataset [25 reported as the mean 95% confidence interval (CI). The low-rank case reports results over 5 random seeds. In the high-rank case, CI is calculated with the final 5 checkpoints of single run. (Note: We omit the numerically unstable baselines VAMPhet-2 and DPNet.)
Low-rank (k 16)
High-rank (k 64)
Algorithm
I1 64, B 384
H1 = 128, B 384
H = 64, B = 96
DPNet-relaxed
7.36
7.85
6.97
VAMPnet-1
9.54
34.76-0.0
19.71
LoRA
10.27
38.89
37. 37.74
LoRA
10.74
39.37
38.5010.80
LoRAug
+0.07
42.44
37.33
Operator Approximation Performance. 2 reports the VAMP-E scores which measure the operator approximation eπrоr the Hilbert-Schmidt norm. In the low-rank case with 16, LoRA variants consistently outperform baselines, with LoRAsea achieving the highest score and lowest variance. We further assessed stability under a high-rank setting with k 64. As shown in the right columns of Table VAMPnet-1 exhibits significant performance degradation in the restricted batch setting. We attribute this behavior numerical instability arising from ill-conditioned moment matrices when using small batches. In contrast, our LoRA variants maintain strong performance across all regimes, demonstrating superior robustness.
Test-Time Orthogonality Analysis. As another performance measure, we report the test-time orthogonality of the learned basis. To this end, we compute the Gram matrix of the whitened basis with respect to the test data, formally calculated as M Magrin [f] pgges Mf] graan [f] 3 and similarly
for g. Ideally, this quantity should be identity matrix.
As illustrated in Figure 3 LoRA variants maintain a relatively cleaner diagonal structure than the baselines, which exhibit sensitivity to distribution shifts. This confirms that the bases learned by LoRA more effectively capture the invariant geometry the system.
90
25
DeNeetrelaxed
VAMPmett11
LoRA
LoRAyn
LoRA
Figure 3: Test-set orthogonality evaluation of the learned bases whitened by training data. LoRA variants preserve clear diagonal structure, demonstrating superior generalization.
5
Concluding Remarks
The evolution of linear-algebraic tools for Koopman analysis, such as DMD, extended DMD, kernel DMD, VAMPnet, and DPNet, mirrors that of CCA and its kernel and The LoRA framework for Koopman operator can be similarly understood as an adoption of the recent advances in correlation analysis and representation learning, which can naturally avoid unstable decompositions and train efficiently with standard mini-batch gradients. Across diverse benchmarks, LoRA variants deliver superior singular subspaces, eigenfunctions, long-horizon forecasts, and scalability, establishing a practical, reliable path for data-driven modeling of complex dynamics.
variants.
While the LoRA framework represents important step toward scalable deep learning methods for analyzing stochastic dynamical systems, key limitations persist. For example, as shown in the chignolin experiment, the quality of learned dynamics is inherently constrained by the data' 's temporal resolution, which can prevent the recovery the slowest physical processes. Furthermore, robustly identifying coherent structures in highly non-normal or chaotic systems remains open challenge.
Future research should therefore focus on advancing the deep learning methodologies tailored for Koopman analysis. This includes developing specialized neural architectures and establishing stronger theoretical foundations for learning from trajectories of dynamical systems. A deeper understanding of the landscape of the LoRA objective will also be crucial for advancing this line of work.
10

---

Acknowledgments
JJR and GWW were supported in part by the MIT-IBM Watson AT Lab under Agreement No. W1771646. MJ and SYY were supported by the Institute for Information & Communications Technology Planning & Evaluation (IITP) grant funded by the Korean government (MSIT) [No. RS-2022-11220111 Development of Goal-Oriented Reinforcement Learning Techniques for Contact- Rich Robotic Manipulation of Everyday Objects, 80%], [No. RS RS-2024-00457882, Al Research Hub Project, 10%], and [No. RS- RS-2019- Artificial Intelligence Graduate School Program (KAIST), 10%].
References
[1] Galen Andrew, Raman Arora, Jeff Bilmes, and Karen Livescu. Deep canonical correlation analysis. In Proc. Mach. Learn. Res., pages 1247-1255. PMLR, 2013. p. I7
[2] Omri Azencot, N. Benjamin Erichson, Vanessa Lin, and Michael W. Mahoney. Forecasting sequential data using consistent Koopman autoencoders. In Proc. Mach. Learn. Res., 2020.
15
[3] Luigi Bonati, GiovanniMaria Piccini, and Michele Parrinello. Deep learning the slow modes for rare events sampling. Proc. Natl. Acad. Sci. U.S.A., 118(44))c21133333111 2021. doi: 10. 1073/pnas. 21 2113533118, URL http////tt////////.ppssoo///////
pp. and 28
[4] Matthew J Colbrook, Lorna Ayton, and Máté Szōke. Residual dynamic mode decomposition: robust and verified Koopmanism. J. Fluid Mech., 955:A21, 20233, p.
[5] Zhijie Deng, Jiaxin Shi, and Jun Zhu. NeuralEF: Deconstructing kernels by deep neural networks. In Proc. Mach. Learn. Res., volume 162, pages 4976 4992.. PMLR, 17-23 Jul 2022. URL htt.hhttss/prrccedinggg mlr. hwwp./eeeee///e///eg/.2bb...
p. 6
[6] Carl Eckart and Gale Young. The approximation of one matrix by another of lower rank. Psychometrika, 1(3):211-216 September 1936. ISSN 0033-3123, 1860-0980. doi: 10.1007/ BF02288367.
[7] Weiguo Gao, Yingzhou Li, and Bichen Lu. Triangularized orthogonalizzaionn method for solving extreme eigenvalue problems. SIAM J. Sci. Comput., 93(63):1 October 2022. doi:
10.10077:109155-22--00005.
p. [6
[8] Jeff Z HaoChen, Colin Wei, Adrien Gaidon, and Tengyu Ma. Provable guarantees for self- supervised deep learning with spectral contrastive loss. Adv. Neural Inf. Proc. Syst. 34: 5000-5011, 2021. pp. and [6
[9] Shinya Honda, Toshihiko Akiba, Yusuke S. Kato, Yoshito Sawada, Masakazu Sekijima, Miyuki Ishimura, Ayako Ooishi, Hideki Watanabe, Takayuki Odahara, and Kazuaki Harata. Crystal structure ten-amino acid protein. Journal of the American Chemical Society, 130(46): 15327-15331, 2008. doi: 10. 10.1021/j88000533. htt.htts:://doi http://.....rrro.oo.o///////011.
PMID: 18950166. p. 28
[10] Harold Hotelling. Relations between two sets of variates. Biometrika, 28(3/4) 1936.
ISSN 00063444. URL BRLL/http://ww/.. jstor. wtt.://taaabb//3335....3 pp. and 18
[11] Takafumi Kanamori, Shohei Hido, and Masashi Sugiyama. A least-squares approach to direct importance estimation. J. Mach. Learn. Res., 10: 10:1391-1445, 2009. p. T6
[12] Tosio Kato. Perturbation theory for linear operators, volume 132. Springer Science & Business Media, 1980. p.
[13] Péter Koltai, Hao Wu, Frank Noé, and Christof Schütte. Optimal data-driven estimation of generalized Markov state models for non-equilibrium dynamies. Computation, 6(1):22, 2018.
3
11

---

[14] Vladimir Kostic, Pietro Novelli, Andreas Maurer, Carlo Ciliberto, Lorenzo Rosasco, and Massi- miliano Pontil. Learning dynamical systems via Koopman operator regression reproducing kemel Hilbert spaces. In Adv. Neural Inf. Proc. Syst., volume 35, pages 4017-4031, 2022. pp.
and X
[15] Vladimir R Kostic, Karim Lounici, Pietro Novelli, and Massimiliano Pontil. Sharp spectral rates for Koopman operator learning. In Adv. Neural Inf. Proc. Syst., 2023. pp. and
[16] Vladimir R Kostie, Pietro Novelli, Riccardo Grazzi, Karim Lounici, and Massimiliano Pontil. Learning invariant representations of time-homogeneous stochastic dynamical systems. In Int. Conf. Learn. Repr. 2024. pp.
18 25
26 and 28
[17] Vladimir R Kostic, Gregoire Pacreau, Giacomo Turri, Pietro Novelli, Karim Lounici, and Massimiliano Pontil. Neural conditional probability for uncertainty quantification. In Adv. Neural Inf. Proc. Syst., 2024. pp. 3 6 and 16
[18] J. Nathan Kutz, Steven L. Brunton, Bingni W. Brunton, and Joshua L. Proctor. Dynamic Mode Decomposition. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2016. doi: 00.113/11.9716111 19745008. https:/ts:///pubs. htte//pbs.siia.....m..... http://o//bbb/o..111/1 9781611974508 p.
[19] Andrzej Lasota and Michael C Mackey. Chaos, fractals, and noise: stochastic aspects of dynamics, volume 97. Springer Science & Business Media, 2013. p.
[20] Kresten Lindorff-Larsen, Stefano Piana, Ron O. Dror, and David E. Shaw. How fast-folding proteins fold. Science, 334(6055)
October 201 I. doi: 10.1126/science. wawiiiinnee.228.0. 1208351. URL
122.2////sccienncc..81111111 pp. and 28
[21] Xin Liu, Zaiwen Wen, and Yin Zhang. An efficient Gauss Newton algorithm for symmetric low-rank product matrix approximations, SIAM J. Optim., 25()):1511-1088, 2015, pp.
and 16
[22] Yong Liu, Chenyu Li, Jianmin Wang, and Mingsheng Long. Koopa: Learning non-stationary time series dynamics with Koopman predictors. In Adv. Neural Inf. Proc. Syst., volume 36, pages 1221111111111 12290, 2023. URL https wwwriinnnn.nuuus....e..... wtw.eeurrps.ccppppepeeeee
15
[23] Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton. Deep learning for universal linear em- beddings of nonlinear dynamics. Nat. Commun. 9(1):9950, November 2018. ISSN 2041-1723. doi: 10. 10.1088//41447..188-721200
URL http:://oo..or///.....////211200
p. P5
[24] Andreas Mardt, Luca Pasquali, Hao Wu, and Frank Noé. VAMPnets for deep learning of
molecular kinetics. Nat. Commun., 9(1):5, January 2018. ISSN 2041-1723, doi: 10. 1038/
s4 44146-0017-02388
UtL/https://doi..rrr htt.:///w......////......071.00. 1 pp. 5
15
17
18
and 29
[25]
Tim Marshall, Robert M. Raddi, and Vincent Voelz. An evaluation of force field accuracy for the mini-protein chignolin using Markov state models, 2024. UtL/httpss//doi.orrg
26434/chemrxix 2024- xmztm, Preprint. pp. and 28
[26] Stanislav Minsker. On some extensions of Bernstein's 's inequality for self-adjoint operators.
Statist, Probab. Lett., 127: 127:111-119, 2017. p. 20
[27] Mirsky. Symmetric gauge functions and unitarily invariant norms. Q. J. Math., 11(1):50-59,
January 1960. ISSN 0033-5606. doi: 10. 10.1093/qmath/11.1 1.550 5
[28) Frank Noé and Feliks Nüske. A variational approach modeling slow processes in stochastic dynamical systems. Multiscale Model. Simul., 11(2):635-655, 2013. doi: 10.1137/110888116 URL http:/ttsss//oi..g///oooo.. 1ttp://...113//1111111
_eprint: doi.org www.oor///01...... 137/1 10858616.. 3
[29] Andrzej Ostruszka, Prot Pakoński, Wojciech Slomczyáski, and Karol Zyczkowski. Dynamical entropy for systems with stochastic perturbation. Phys. Rev. E, 62(2):2018, 2000. p. 25
12

---

[30] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang. Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An imperative style, high-performance deep learning library. In Adv. Neural Inf. Proc. Syst., Red Hook, NY, USA,
2019. Curran Associates Inc. pp. 23 and 24
[31] Fabian Paul, Hao Wu, Maximilian Vossel, Bert I. de Groot, and Frank Noé. Identification of kinetic order parameters for non-equilibrium dynamics. J. Chem. Phys., 150(16), 2019. p p
[32] Jon Ryu, Xiangxiang Xu, Hasan Sabri Melihcan Erol, Yuheng Bu, Lizhong Zheng, and Gregory W. Wornell. Operator SVD with neural networks via nested low-rank approximation. In Proc. Mach. Learn. Res., 2024. URL http:://oppnreeeeeeete//toorrrrd
pp-
16 and 25
[33] J. Jon Ryu, Samuel Zhou, and Gregory W. Wornell. Revisiting orbital minimization method for neural operator decomposition. Adv. Neural Inf. Proc. Syst., 2025. URL https://
p. T6
[34] Erhard Schmidt. Zur Theorie der linearen und nichtlinearen Integralgleichungen. Math. Ann., 63(4):433 476, December 1907. ISSN 5831, 1432-1807. doi: 10. 10.1007/BF014499700 p.
[35] Kristof T. Schütt, Stefaan S. P. Hessmann, Niklas W. A. Gebauer, Jonas Lederer, and Michael Gastegger. SchNetPack 2.0: A neural network toolbox for atomistic machine learning. J.
Chem. Phys. (1558114)444401, 04 2023. ISSN 0021-9606. doi: 10.1063/5.0138377. URL
https:://oi..oo org/ 0163
p.
30
[36] G w Stewart. On the early history of the singular value decomposition. SIAM Rev. Soc. Ind. Appl. Math., 35(4)5551-566 December 1993. ISSN 0036-1445, 1095-7200. doi: 10.1137/1035132 34.
16
[37] Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi. Learning Koopman invari- ant subspaces for dynamic mode decomposition. In Adv. Neural Inf. Proc. Syst., vol- ume 30, 2017. URL htRp://cpp:: neuriiss wwwrips.cc/pppp........... hwtpc/ppee.lllee/eppeeeeeeeeeeee wwt.ff///pppee/o////iille 3885582557554555ffffff5555858333 Paper.
p. 15
[38] Joel A Tropp. An introduction to matrix concentration inequalities. Found. Trends Mach. Learn., 8(1-2):1-230, 2015. ISSN 1935-8237,1935-8224 doi: 10.1561/2200000048. p.
[39] Yao-Hung Hubert Tsai, Han Zhao, Makoto Yamada, Louis-Phillipe Morency, and Ruslan Salakhutdinov. Neural methods for point-wise dependency estimation. In Adv. Neural Inf. Proc. Syst. 2020. pP- and [6
[40] Jonathan H. Tu, Clarence W. Rowley, Dirk M. Luchtenburg, Steven L. Brunton, and J. Nathan Kutz. On dynamic mode decomposition: Theory and applications. J. Comput. Dyn., 1(2):991- 421, 2014. ISSN 2158-2491. doi: 10.3934//jc..20114..99 URL hig.
I
[41] Giacomo Turri, Luigi Bonati, Kai Zhu, Massimiliano Pontil, and Pietro Novelli. Self-supervised cvolution operator learning for high-dimensiona. dynamical systems. arXiv, 2025. URL
https://arxiv, org/ 2000s///5555.00 18671 p. 16
[42] Lichen Wang, Jiaxiang Wu, Shao-Lun Huang, Lizhong Zheng, Xiangxiang Xu, Lin Zhang and Junzhou Huang. An efficient approach to informative feature extraction from multimodal data. In Proc. AAAI Conf. Artif. Intell., volume 33, pages 5281-5288, July 2019. doi: 10. 1609/aaa
v33101.33015281. pp.
and 16
[43] Rui Wang, Yihe Dong. Sercan O Arik, and Rose Yu. Koopman neural forecaster for time- series with temporal distributional shifts. In Int. Conf. Learn. Repr., 2023. URL https!
p. 15
[44] Ziyu Wang, Yucen Luo, Yueru Li, Jun Zhu, and Bernhard Schölkopf. Spectral representation learning for conditional moment models. arXiv preprint arXiv:2210.1655 2022. p. I6
13

---

[45] Zaiwen Wen, Chao Yang, Xin Liu, and Yin Zhang. Trace-pentt penalty minimization for large-scale eigenspace computation. Sci. Comput. 66: 175-1203, 2016. pp. 5
[46] Matthew O. Williams, Ioannis G. Kevrekidis, and Clarence Rowley. A data-driven approxi- mation of the Koopman operator: Extending dynamic mode decomposition. J. Nonlinear Sci., 25, 08 2014. doi: 10.1007/s00532--155 9258-5, pp.
and 22
[47] Matthew O. Williams, Clarence W. Rowley, and loannis G. Kevrekidis. A kemel-based method for data-driven Koopman spectral analysis.
J. Comput. Dyn., 42):247-265, 2015. ISSN
2158-2491 doi: 10.39 10.3934/jed.2015005. URL
ce535396 f-fe-4fff-0 b6de- 36de-baaf996f6193
pp. and
[48] Hao Wu and Frank Noé. Variational approach for learning Markov processes from time
series data. J. Nonlinear Sci., 30(1):333333 February 2020. ISSN 1432-1467. doi: 10. 1007/ s00332-015 019 01952-0-9-09567-1 URL
pp.
7
5
18
and 29
[49]
Xiangxiang Xu and Lizhong Zheng, Neural feature learning in function space. J. Mach. Learn Res., 25(142):1-76, 2024. pp. 5 and 16
[50] Bernt Øksendal. Stochastic Differential Diifarentiallguuuainss.. Springer Berlin Heidelberg, 1995. doi: 10.
1007/978-3-662--031888 /tLl//tpp//dxx.dd..r////0.... 1087998-- 662 03185-8 p 27
14

---

Appendix
Related Work
15
A.1 Deep Learning Approaches for Koopman Operator Estimation
15
A.2 On Low-Rank Approximation
16
A.3 On Nesting Techniques
16
Printic for VAMPnet and DPNef
17
B. VAMPnet
17
B.2 DPNet
18
Deferred Technical Statements and Proofs
18
C.1
Derivation Low-Rank Approximation Objective
18
Statistical-Learingg Theoretic Guarantee
19
Special Case of Finite-Rank Koopman Operator D.I Singular Value Decomposition
21
D.2 Eigen-analysis
22
Extended DMD and Multi-Step Prediction
22
Implementation
23
Implementation of Nesting Techniques
23
F2
Pseudocode for LoRA Loss with Nesting
24
Deferred Details on Experiments
25
G.1
Noisy Logistic Map
25
G.2
Ordered MNIST
26
G.3
Langevin Dynamics
27
G.4
Chignolin Molecular Dynamics
28
G.5
Computing Resources
30
A
Related Work
In this section, we discuss some related work better contextualize our contribution.
A.1 Deep Learning Approaches for Koopman Operator Estimation
Deep learning models are increasingly being integrated into the Koopman theory framework to enable the learning of complex system dynamics by leveraging the expressive power of neural networks. One of the earliest examples is Koopman Autoencoders (KAEs) for forecasting sequential The KAE approach in general aims to learn latent embedding and global, time-invariant linear operator, guided by reconstruction and latent-space constraints. While these methods laid important groundwork, more recent developments 224 superior performance. Consequently, our experiments focus on these more advanced techniques.
data 47 37 23 2
15 [6] have demonstrated
Recently, various deep learning architectures for learning Koopman operator have been proposed to handle complex real-world environments, such as non-stationarity in time series data. For example, Liu et al. [22 introduce Fourier filters for disentanglement and distinct Koopman predictors for time- variant and invariant components. Another example Koopman Neural Forecaster (43). which utilizes predefined measurement functions, transformer-base local operator, and feedback loop. We remark, however, that while these methods demonstrate impressive empirical results on some datasets, they mainly prioritize forecasting accuracy rather than eigen-analysis essential for understanding the underlying system and achieving robust prediction, which lies the core of the Koopman framework. In this work, focus this core aspect, leaving forecasting extensions for future study.
15

---

A.2 On Low-Rank Approximation
The idea of LoRA as a variational characterization for SVD of a matrix or operator has a century- long history; see [36] for overview. Despite its simple and elegant formulation, dominated by the common Rayleigh quotient maximization framework and its variants, its wide adoption as a variational framework for computing singular subspaces based on optimization started only recently. For example, see
45 45 from the numerical optimization literature, [ 44) in representation learning,
and [22 49 17 for correlation analysis. We refer the reader to 32 Appendix B] for more detailed
overview of the literature.
Our learning technique closely follows the wnw-rann-approximatio NeuralSVD in [32], including the nesting techniques. Ryu et al. [32] proposed the LoRA framework with nesting generie tool to perform SVD of a general compact operator using neural networks including scientific simulation and representation learning. While they also considered decomposing the density ratio p which they call the canonical dependence kernel (CDK), as a special case, ply)
framework dubbed as
they did not explicitly consider its application to dynamical systems. We note passing that there exists a similar squared loss studied in density ratio estimation [1],9]] but it becomes only related p(z.g)
to spectral decomposition when the underlying density ratio is in the form of CDK p(z)p(a))
More recently, in concurrent work, Turri et al. 41 also proposed to use the LoRA loss, emphasizing self-supervised learning perspective on learning singular functions. Specifically., they parameterize the srroder-laggedeeeee
pair by defining the latter as g(x') - Pf(x'), where P kg**k is a learnable matrix and f(.) is a neural network encoder. For a fixed encoder f(-), they show that the optimal under the LoRA objective coincides with a sqaares estimator for the CDK. A systematic study of the benefits of this parameterization is an interesting direction for future work.
During the course of independent study, we also identified an unpublished objective named EYMLoss in the GitH repository of Kostic et al. (16), which effectively corresponds to the LoRA loss investigated in this work; see hhttps://githute https////ggthuuu.cm//aahhahhiihe www.eaai...aae......... Learning Sasaammical--Sseeemm
L79C7-L79C14 The associated commit message indicates that the authors were unable to obtain satisfactory results with this formulation. We attribute this failure to a subtle but critical oversight in handling the constant modes. Recall that a Koopman operator, as a special case of CDK operators, admits constant functions as its leading singular functions. In our implementation, we explicitly account for this structure by prepending constant components (i.e., appending 1's) the outputs of both the encoder and the lagged encoder. This ensures that the constant modes are appropriately represented and preserved throughout training and inference. contrast, EYMLoss incorporates the constant modes by projecting out their contribution from the observables. Specifically, it replaces f(x) with f(x) - (x) [f()) in the LoRA loss, where the expectation is approximated empirically using minibatch samples. This approach implicitly enforces orthogonality to constant modes by requiring the learned singular functions to be zero-mean with respect to the data distribution. While this treatment is valid in principle, it necessitates that the encoded features remain centered not only during training but also at inference time. Moreover, the constant singular functions must be explicitly included in the final model representation. These additional steps do not appear to have been fully implemented in the kooplearn codebase, which may have contributed to the method's lack of empirical success.
Finally, we remark in passing that a recent work studies another unconstrained variational optimization framework, called the orbital minimization method, for decomposing positive definite (PD) operators using neural networks [33]. As the framework is restricted to PD operators, it is not directly applicable to the setup considered in the present manuscript. Developing an extension of this method for analyzing stochastic dynamical systems would be of theoretical interest.
A.3 On Nesting Techniques
The sequential nesting technique proposed in [32 can understood as a function-space version of the triangularized orlhogonalizationn method studied in 77. 'he joint nesting technique was first proposed specifically for CDK by Xu and Zheng [49], and later extended as generic tool in Ryu et al. (32). For more discussion on the history of nesting, refer to (32..
16

---

B Variational Principles for VAMPnet and DPNet
In this section, we review the principles behind VAMPnet and DPNet in more details.
B.1
VAMPnet
VAMPnet (48 24 is based on the following trace-maxiiiizaaiio--yye characterization the top-k: singular subspaces.
Theorem B.1. Ler т: 12 ()) + (2) be a compact operator with singular triplets {(o, o.(.), (.())))))) i.e. (T)))(z) oz),(z), (T*))) = o,4b(y), and 2 2 ... 2 0. Then, for any fixed r 1,2,... the optimizers of the following optimization problem
Ei. T9
(fraa)) CT (2)) (P)
1=1
subject to
M p. [ Mp g
(11)
characterize the top-k singular functions T up an orthogonal transformation.
We refer an interested reader to the proof in 48
For the particular case of 'T being a Koopman
operator x: I2 (X) I.2 (X)), Wu and Noé 48 named the objective
R, gl &
E (f. Kor)
(12)
as the VAMP-r score of f and g. is noted that VAMP-1 score is identical to the objective of deep
CCA M. Wu and Noé [48.. Mardt al. [4 advocated the use VAMP-2 score, based on its
relation to the Tz-approximation error.
In the original VAMP framework 48 24 neural networks not directly trained by the VAMP-r
score. Instead, they considered two-stage procedure as follows:
Basis learning: They parameterize the singular functions
fx) (Tf.
gv,0 VTga(x';
where fx) A (x). fo (x))T= Ro and gø(x') A (ge, (x'),.. 99.k. (x))) & denote some basis functions, which in the trainable basis by neural networks, and U bo and V E R& denote the projection matrix to further express the singular functions using the given basis. Since the jointly optimizing over all 0, U, V is not practical, to train fa and ge. they plug-in fu.a(x) (x) and (x.) into Eq. (II). and consider the best U and by solving the partial maximization problem, i.e.,
is typically parameterized
max R. (b. +) ||| (M [f)) (f))=+ T[f. (g)) U.V
Here, || ll. denotes the Schatten r-norm. This objective what we call the (maximal) VAMP-r score in our paper, which is the objective in the VAMP framework that trains the trainable basis functions.
Inference: Once fa and go are trained, Eq. (11) become an optimization problem over U and V:
max U.V
RU,
subject to UTM no FJU =
VTM igi I.
where
(U 2 Σ (u.
17

---

Wu and Noé 48 proposed to solve this using the (linear) CCA algorithm [00,, and calling the final algorithm feature TCCA. To avoid the degeneracy, basis functions need to be whitened. Strictly speaking, however, the feature TCCA algorithm has nothing to do with the optimization problem. We note that our CCA+LoRA approach is built upon this inference method in the VAMP framework.
Let 01....,0K 0% denote the top-k singular values of the approximate Koopman matrix K & (M,) [M][f)))= [f]) g (M, [g)) After CCA, we can consider the rank-k approximation of the under- lying Koopman operator X. as we considered our Approach I in Section 3.2 If we call the K. denote the corresponding approximate operator. Wu and Noé [88] called the (shifted and negated) approximation error in the Hilbert-Schmide norm
RE
|| || X
which called the VAMP-E score. While this looks essentially identical to low-rank approximation enror consider in this paper, Wu and Noé [8]] and Mardt et al. (24) used this metric only for evaluation, but not considered for training. In Appendix [.44 we provide a precise definition of the VAMP-E score for completeness.
B.2 DPNet
As introduced in the main text, the DPNet objective is given as
c cl) (f,g) - (M (f)) T(f, gl (g)) || ((R (f) + RMM (gl))-
Here, the metric distortion loss R. s + R2 is defined as RMM) - tr(M2 M - In M) for M > 0.
If y and Mp [f] and [o are nonsingular, the objective becomes equivalent the VAMP-2 score.
To detour the potential numerical instability in the first term of the DPNet objective, which is essentially the VAMP-2 objective, Kostic et al. [16] further proposed relaxed objective
[T[fg]]k:
dar (f.
+ pR PM (f)) + RM,
||
f op Mp,
glllpp
Kostic al. [16 proved the following statement:
Theorem B.2 (Consistency of DPNet objectives). Lery 2 0. If is compact,
cla) (f. g) 2 1/p (f.g) 2 - c.(T)
The equalities attained when f(x) and g(y) the top-k singular functions of T. If is Hilbert- Schmidt and ( > 0k+1 (T) and 7 >> 0, the equalities are allained only if f(x) and g(y) are orthogonal rotations of top-k singular functions.
For the continuous- dynamics with a self-adjoint Koopman generator L, the authors proposed to use
E(1) dp (f) é - (r)(M, Lf)) + (R,M,,
C Deferred Technical Statements and Proofs
In this section, provide short derivation of the low-rank approximation (LoRA) loss for com- pleteness, and present a learning-theowe theoretic guarantee for the LoRA objective.
C.1 Derivation of Low-Rank Approximation Objective Proposition C.1.
2
min f.g
x
Σ
fz X
gi
min Ciora (f, g).
HS
f.g
18

---

Proof. Recall that Ciora (f.g)
-2tr(([ff, gl) + tr [f (g)). To prove, first note that
x Σ
f 8
| iis = - Σ (Fs, Kg:) Σ
(fe. Po 9))
HS
Here, the first term can be rewritten as
Σ 1=1
(fi Koh)
=
E px 1x))))xxxx)((((()x(())))))))))))
==
tr (
-r((T[f.g))
Likewise, the second term is
(f. (a Iila tr( (x)))))))))))))(((((()))))))))))
= prr (F (g).
This concludes the proof.
C.2 Statistical-Leariingg Theoreetic Guarantee
We consider empirical estimates of the second moment matrices T[f, g). M [f], Mp. [g].
f(f.g
g)
A
f(x))g(x()
N 1=1
1
M no
[r]
((x))(())),
No
Mp
[g]
Σ B(x))g(x))
N
Here, xx), px(x)p(xx) x) are i.i.d. samples and pa(x) and X pi(x') are i.i.d. samples.
Note that, if we follow the common data collection procedure, X; px(x) and x pi (x) are drawn from single trajectory, and thus the samples (x, x') cannot be independent. The independence assumption between {x:} and (x) } are for the sake of simpler analysis. In practice, however, this can be enforced by splitting the set of given independent trajectories into two subsets and compute M,, [f] and M,, ig] with different subsets of samples.
Theorem C.1. Let [[f(x)]| < and [g(x')ll I R almost surely for some R > 0. Let Vp A
Epotso (()) - {ng(f))²
[ E Rexk and Væ A ((xxte/(ls((')xxxx" M,, { E k
Then, with probability at least
we have
16
6
6
j2 8--1 C R2
log
log
3N
Vellop
12r(Vr)
2R
12r(Vr)
+
log
log
No
3.N
2
Vg
12r(Vg)
2R²
12r
12r(Vg) g)
+
g
op
log
log
N1
5
3.N
8
To prove this, we will invoke Bernstein inequalities for bounded random variables and "bounded" self-adjoint matrices:
19

---

Lemma C.1 (Bernstein's inequality [88 Theorem 1.6.1]). Let X Xx be i.i.d. copies of random variable X such that [X] I almost surely, [[X] 5 E (0, 1), with probability at least 5,
0, and B[x²] I o². Then, for any
n
4u
2 +
40 2
X
log
log
3n
n 8"
Lemma C.2 (Matrix Bernstein inequality with intrinsic dimension [26). 38 Theorem 7.7.1J). Let
......, An be ii.d. copies of a random self-adjoint matrix A, which satisfies || almost
surely, [[A] = 0, and E[A?] = V = 0. Let r(V) & sx(V) denote the intrinsic dimension V. Then, [Vllop
for any 8 E (0, 1), with probability at least 1 5.
2c 3n
4r(V)
4r(V)
log
log
72
8
n
We are now ready to prove Theorem C.1
Proof of Theorem
C.1
Consider
22 4 = [r(- ( 2f(.. + [M tr(- g]] + M,, [.] [g))
tr] g] g]f, + [Mr FJM [g] Mp [f] [g]]]
For the first term [r[[[[]g]] T[f, g))l note that
tr (fr.g) - T(f.g] = N Σ fgxx) Bpo(xinxe (x) [f(x)
1
Here, note that
[f(x) [[x)) - [|f(x)]l (f(x))l-Ilg(x') R ((())) R
Hence, we can apply Bernstein's inequality in Lemma C.1 for g((x)) with +- R.3
o² Re, and + 8', where 8' be decided the end of proof.
Now, we wish bound (tr [f] [g]
[fo[f]MM [g))]-
Let Ap A Moa [f] Mf) and Ag A M [g] M, [g) Then,
tr (r, []Mp, [g] Mp. FM [g))
tr
M p. igl + MA, {]
tr
Ay M
[ + tr
tr(M M [f]
3
A tr
[M [g]) + ||Agllop tr(M, po [f)
1 N.
||
Ex llg
Agllo N;
Σ
|g(()
(
< R² (IArllop |A
where (a) follows from the inequality tt((A)] < |Allop tr(B) for square matrices A, such that B ž 0, and (b) from the boundedness assumption, i.e., [|f(x)] < R and [[g(x')] < R almost surely. Here, we can apply the matrix Bernstein inequality in Lemma C.2 to bound | with <- fxx)f(xxT Mp, [f] with 8
8', and similarly for ||Agloop-
Finally, by applying the union bound on the Bernstein inequality and two matrix Bernstein inequalities with 8' 8/3, we can conclude the desired bound.
20

---

D Special Case of Finite-Rank Koopman Operator
If a given operator is of finite rank, we can analyze the operator using finite-dimensional matrices of the same dimension. We can develop a computational procedure for SVD and eigen-analysis under the finite-rank assumption, and apply the tools once we approximate the target operator using low rank expansion. Some of the results established in this section are used later when (numerically) computing the ground truth characteristies of the noisy logistic map studied in Appendix [
Suppose that the Koopman operator X: 12 (2¹) - 12 (2²) of our interest is of finite rank. Specifi- cally, the corresponding kernel can be expressed in the following factorized form:
p(x' x) (x)) = (x)3((xx))
4xx,X' x')
Pr (x()
Due to the separable form, the kernel has a rank at most r. We aim to find its SVD as follows:
1(x,x') Σ (x)((())
(13)
where the functions satisfy the orthogonality conditions: (6), ky) (wi, Vyyer = 6.g The singular
values 01 =1 1 2 02 2 ... 2 or 2 0. The first singular functions 01 and v1. corresponding to the trivial mode 01 = 1, are constant functions.
D.1 Singular Value Decomposition
We can compute the singular value decomposition (SVD) using Mno [a] and Mp, [3] as follows. Theorem D.1. The SVD of the matrix
Ssart [a, 3] 2 [M)a /3)
shares the same spectrum, i.e, there exist orthonormal matrices U E Rkxr and E Rkxr and Σ diag(o1....,.,,, such that
[a. = UEVT Σ ou..vI.
S Ssqrt
The ordered, normalized singular functions of the kernel k(x,, x') are given as
((x) UT (M)
(14)
(xx') = VT(M,,, 11)
(15)
Proof. First of all, is easy to check that Mpr Mp. I. Next, we can show Eq. [33 To show this,
(x)x()))ssxx)) +f(x
(UT(MMplaa) (1 la, (M,, (0
(M (MM (M) (M)) =xxx)T((x')
=k(x,x').
To show (X)))(xx) = (x)),,)), consider
(x') A / (0x)))))
= =((x'' (al))
21

---

B(x') / (x())(() 1
(M)
-3(x'T (12)))1/2(((00))))) Polpi
[ax, T
Poips
o, [3))-1/2v
We show (x* (x)) 0;0; (x) by the same reasoning.
D.2 Eigen-analysis
can compute the eigenvalues and eigenfunctions of the kernel by eigen-analyzing the matrices Mo [3, ad] and M [3, ox.
Theorem D.2. Let
C be right eigenvector of Mp. 13, ai] with eigenvalue E C, i.e.
Mo 1,, [/ Aw. If we define n(xx) & wTx(xx), then n right eigenfunction of with eigenvalue X, i.e. Kn An).
Proof. Since x(x,x') ((x)),(xx), have (xCn)
/ /
=x(x) / (xx))a(x(TT
Anj(x).
Theorem D.3. Let Ck be a left eigenvector of Mp 13, with eigenvalue C, iee. 'p Mo [3, a]] Az" we define ((x') 2 3(x',, then Ç is a left eigenfunction of K. with eigenvalue X, i.e. X-C = XÇ.
Proof. We have
(xx))(x / 1(x, x(xx)))(())))()xxx
= =B(x'1 / (x(x))((x)xxxl
=3(x')T az)
= =T
= =(((x'.
Interestingly, this implies that spectrum of M,, [3, o)] and M ai] must be identical.
Moreover, if yo denotes the Perron Frobenius eigenvector of Mp, [3, al (i.e., the left eigenvector with eigenvalue 1), then (x) yI3(x) is the stationary distribution.
E Extended DMD and Multi-Step Prediction
In this section, we illustrate the extended DMD (EDMD) procedure [46) using basic ideas as elementary linear regression. Suppose that we are given top-k left and right singular functions f(-) (encoder) and g(-) (lagged encoder), respectively. If h E span(f) & span ff then there exists some & such that l(x) =f(().. Then, denote by Kr the restriction of the Koopman operator onto the span, which acts on /2 as
(Xrh)(x) f(x)
22

---

for some KF E R With finite data, there are two issues in this picture: (1) How can we compute Kp from data? (2) Given span(f), how can we compute the corresponding z R
For the first question, the EDMD aims to find the best finite-dimensional approxima- tion of the Koopman operator K restricted on span(f), in the sense that (Xff)xx Ep(x' x)) [f(')] = Kf(x). As natural choice, we can choose the ordinary least square solution that solves f(x') = (x(()) for (x, x') ~ Po(x)p(x' x), where
k A F+F = (FTF) FFTF' = (Mpa(r)) (T())
where we denote the data matrices by
E IR Nok
and F' [fx)))... E NR
Note that M.. [f] and [[f] are the empirical estimates of the second moment matrices Mn, [f] and T(.).
For the second point, we can again view this as a linear regression problem, since hxx:) 2 (x;) for i 1,....: N.
Hence, can define the best again the OLS solution 2%, [f] & F+H (Mpo fl) FFTH, where
2 [h(x)), h(xN)] E RR
Now, suppose that already computed T(f, g] and 2n from data. Then, can approximate the multi-step prediction as
Ep(x, xx)) (x))) [xxx)] P f(xo) TT[ g)
(16)
The same logic applies to the subspace spanned by the right singular functions vb-).
F
On Implementation
In this section, we explain how we can implement nesting with automatic differentiation and provide pseudocode in PyTorch 30
F1
Implementation of Nesting Techniques
We can implement the sequential nesting by computing the derivative of the following objective using automatic differentiation:
S (f.g) _ tr Agl) + tr [F] [g]).
where we define a partially stop-gradient second moment matrix
fi) (f2, sgl fil) (fs. sgl fill
(sg sg(fil,fa) f2,
sg[ft]. fa)
(sg[f1), fx)
(sg]f2) falp
(sg fo]. fk)
5 g. &
(fs, sgifall
(fs, fa)
Asglfsl. fk)p
(ffx, sg(fi))
(fe, (gu,e))
(fre fu)
This can be implemented efficiently with almost no additional computation. We note that the imple- mentation of the sequentially nested LoRA provide here is more direct and simpler version than the implementation of NeuralSVD in which implemented the sequential nesting by custom gradient.
The joint nesting can also be implemented in an efficient manner via the matrix mask, as ex- plained in (32]. Define the matrix mask P E R as P.j = mmax(i.j) with m, A Σ/ aj. Then, cjoint (f, g: ax) 2 E aiLiora(f1. (fi.. gi:)) Then, we can write
k
cjoint (f, g; a.) Σ
an (fin, = ur ur(Por 2M [f, Kg] + 1g) (s))
i=1
23

---

F.2 Pseudocode for LoRA Loss with Nesting
Based on the nesting techniques described above, here we provide a simple and efficient Py PyTorch 30 implementation of the Nes (i.e., LoRA with nesting) objective.
class NestedLoRALoss:
def i
3
self, use le nesting-None, n_mades=None, ):
self.use_ Less__eernee_seee _aceddssssssssse use_ @ss_leearned_stee svals
assert nesting in [None, 'seq', ']nt']
self. nesting nesting
self.nesting 'jnt':
12
assert n_modes is not None self.vec_mast self.mat_mas. = jos.oont_nestinnnass..
weightss-np.oness(nndes)]] / n __modes,
else:
self.vec_mask, self mat .ma__mask None, None
self. www.koic_eeeulaizatiiooo. kostic_regulariaatee
__call___see f: forch.Tensor. g: torch.Tensor). # [b, k]
21
22
# g:: [b, k]
if self nesting 'jnt':
\sum_ [i=1)"k <f_i, g_i
26
corr_ted term = -2 (self.vec_me) mask. to(f.device) + f g).mean(o). # \sum_[i=a]"K (sum_(j=a)** <f_i, f_j> <g_i, 8_}>
27
tr cov (f_(a:k)) cov(g_(1:k)))
M_f compute compute_second second (o
M_g compute_sede seruutee_econd___m seeosnd_mmmmnt((g)
31
metric term (self .att laa._mmaasss___ wttpeaaasss.fff...e.oeeeooooeeee to(f.device) M_f (_g).sum()
else:
\sum {i=1)"k i, T ggi>
35
corr_ term = -2 (f g) me)an()).sni)))
36
# [jum_(i=a)* \sum_ (j=a) <f_i, f_j> <g_i, g_j>
tr ( cov (f_(ack)) co:(g_(1:k))
M_f compute_ ce seq_nesting=self feettinnneeenngg.....ee 'seq')
39
M_g compute_ second moment(g, seq_nesting mesting .. 'seq" metric term (M_f M_g).sum()
40
corr_term metric_term
44
def compute
compute_second second (
45
f: torch. Tensor, torch. Tensor None, seq_ ne__nesting: bool False
48
-> Norch.Tensor.
29
compute (optionally sequentially nested) second- moment matrix
M_ij <f_i, g_j>
with partial stop-gradient handling when seq_nesting True.
53
54
args
(n, tensor
g: (n, k) tensor or None
58
seq_nesting bool
if is None:
24

---

s.shape[o]
63
if not seq_neeesee nesting: return (f.T a g)
64
65
else:
# partial stop gradient
# tower-triangular: <f_i, sg[g_j]> for i
68
lower torch. corch.tril(f.T @ g.detach(), diagonal=-1)
69
# upper-triangular: ssg[f_i], 8_j> for < upper torch. (ranchhttii(..deach((.....). g, diagonal=+1)
# diagonal:
<f_i, g_i> (no stop-grad)
72
diag forch.diag((1 s)).su((im=oo) return (lower diag + upper) /
73
74
def get get_joint_nng wtt___nstinn_mssssweeeeee
nddrrayyy
vector_mase mask
77
vector_mast torch.
78
matrix_mases mask torch.min minimum
79
pector_mask.unsquee unsqueeze(1), vector_mas eetto__ask.unsquoozeeeee f.float()
return vector mask, matrix_mast
G Deferred Details on Experiments
In this section, we provide details for each experiment. See 3able for an overview of the benchmark problems from [16]. Our implementation builds upon the codebase of [16 the kooplearn package"
and the codebase 32
Table 3: Overview of the experiments.
Examples
Time
Spectral complexity
Stationarity
Noisy logistic map Ordered MNIST
discrete
non-normal
(nearly) stationary stationary
discrete
normal
ID SDE
continuous
self-adjoint
(ncarly) stationary
Molecular dynamics
continuous (discretized)
normal
non-stationary
G.1
Noisy Logistic Map
We consider a noisy logistic map defined as
X.+1 (rXi(1- X.) + &.) mod 1
p(E). where p(5) A CN (OTE is the order N trigonometric noise over & [-0.5, 0.5] 29 where CN 2 // 7/B( M+1 ()) is the (reciprocal) normalization constant.
for E
Although the dynamics non-normal, structure of the polynomial trigonometric noise ensures that the associated kernel is of finite rank N 1. Concretely, the transition density can be written as
p(x'(z) Σ (x())((('',
2=
where
Bx(x)) /ow (M) cos't (Te) sin N ~i(Trr)
25

---

(x) a ((((())).
Let 3,(x') (x')) (1())
such that we can write p(21a) = (x))
Computation of Ground- -ruuh Properties. Since the operator is of finite rank, the underlying singular functions and eigenfunctions can be computed numerically, using the theory developed in Appendix D Let
[ af A E. (x)x)
Ma [1 E (x)))
Note that
aM (x)) (x))
((xx)((x)T((x)((x
(xo())p(x
((x) (x
E (x) = T[Ø].
Then we can compute the approximate Koopman matrix as the ordinary least square regression, which is equivalent the EDMD:
Kols [+)) (T) (Mx [Ø [M 3|T.
Experimental Setup. Setting N 20, we generated random trajectory of length 16384 to construct the pair data. We used multi-layer perceptron (MLP) with hidden units 64- 128-14 and leaky ReLU activation for left and right parametric singular functions. We evaluated VAMPnet-1 (bluc), DPNet (green), DPNet-relaxed (green with marker x), and LoRA (red) across varying batch sizes B E {256, 1024, $192} and learning rates E {10-3 10-*} using the Adam optimizer; see Figure 4 Each configuration with different random seeds and report the average values with standard deviation. We set the number of modes to k 20.
trained for 500 epochs. We conducted five independent runs per configuration
Results. To evaluate the quality the learned singular basis, we first estimated the singular values à; via CCA. Since the ground-truth spectrum is known in this setting, we computed the relative error in
estimating the squared singular values, defined as lo2 |=2-621 for each i E [k.] As shown in the first row of Figure 4 LoRA consistently outperforms the other methods in capturing the singular subspaces according to this metric, except in the configuration with large batch size = 8192 and a small learning rate 10-4. We also note that the performance of VAMPnet- I is highly sensitive to the batch size, whereas DPNets exhibit relatively robust performance across settings.
Next, we assessed the quality of the estimated eigenvalues. Given the CCA-aligned basis (Φ. }b we computed the Koopman matrix via EDMD using the first i basis vectors b and extracted the corresponding cigenvalues (2 Since the true system has three dominant eigenvalues λ1 λg As we evaluated the estimation quality using the directed Hausdorff distance max minlees] P Ayl. following Kostic et al. [6). The results are presented in the second row of Figure
As expected, increasing the number of singular functions improves the estimation quality across all methods. Kostic et al. [16] reported bascline value of 0.06 achieved by DPNet-relaxed (indicated by the gray, dashed horizontal line), and nearly all our configurations outperform this base- line. We observe that DeNet-relaxed and LoRA yield comparable performance, with DPNet-relaxed occasionally achieving the best results. We attribute the discrepancy between singular subspace quality and eigenvalue estimation performance to the non-normal nature of the underlying dynamics.
G.2 Ordered MNIST
Following the setup of [6], we generated two independent trajectories of length 1000. We used one trajectory for training and the other for evaluation. We used Adam optimizer with learning rate 10-3 batch size 64, for 100 epochs. The convolutional neural network we used is same as (16). namely, 16] set the metric deformation loss coefficient be as suggested for DPNet and reNeeerelaxed.
[2]-o[[[dd]][[[]]]]]]]]] [Re [cpppooo[[]]]]]]]]]ppp]]]]]]]p]] [z]=>Linear] [o]]. We
We note that in the original implementation of Kostic et al. [6 the singular values computed as if
p(alla) axxx) ((x)). which led to wrong singular values.
26

---

B-256,
T- B-1024,
3PMee-mbus LP1A
style TONDY: SUMDoy
10-1
Landtiess
Fandtiens
sisguaar funetors
fandtions
Surellors
Figure Summary of the noisy logistic map experiment. The first row reports relative error in singular values, and the second row the directed Hausdorff distance of the estimated eigenvalues to the three most dominant underlying eigenvalues. The shaded area indicates =0.5 standard deviation.
G.3 Langevin Dynamies
Consider the following stochastic differential equation (SDE):
dx, A(X;)d + B(X)),W
(17)
where A: Rd + IR is a drift term, B: R R is the diffusion term, and W, is m-dimensional Wiener process. Since the diffusion process can be modeled as reversible stochastic differential equation (SDE) with respect to its stationary distribution adjoint within the Hilbert space endowed with the inner product (;,)) [50]. For the Itô SDE in Eq. (17), we can write the action of the generator L (also called the Kolmogorov backward operator, or the derivative) of the diffusion on a function f: Rd I can be written as, by
the associated generator becomes self-
formula [50].
1 + tr ((X
(x)) & Vxf(x)TA(x) +
Special ID Example. Consider simple ID case where
A(x) = _ -('(s)) and B(x)= = 2
In this case, we can write
(Lf)(x) - (=0 7
For simplicity, we set KgT = 1 and 7 0.1.
We the Schwantes potential (xx) defined as
U(x) 4 (x 0.8c 80z? +.2e -..2e-80(2-0.35) 80(2)0.)) 0.5c 40(x+0.5))
whose derivative is given
(xx = -3.r 51 2re 80.² 128(x 0.5)e 80(==0.5)
160(x + 0.5)e 40(2+0.5)
27

---

G.4 Chignolin Molecular Dynamies
This section provides the detail of our experimental setup including the data, as well as sup- plementary analysis, which demonstrates the stabilization effect of the nesting techniques during training.
Experimental Setup. Our work utilizes a publicly available trajectory dataset of the classic chignolin peptide (sequence: GYDPETGTWG, PDB: IUAO) simulated 300 K, provided by [25] In contrast, contemporary studies R 16 analyze simulation data at 340K of the stabilized CLN025 variant (sequence: YDPETGTWY, PDB: 5AWL; (9) from [20), which not publicly available. To maintain maximal consistency with the prior works, we selected trajectories from this dataset generated using the CHARM22* force field and TIP3P solvent, but with halved sampling interval of 100 ps. Moreover, unlike the non-public data created from a single trajectory, the full dataset contains approximately 300,000 snapshots from 34 trajectories, initiated from either folded or unfolded states.
Due to the short trajectory lengths, transition events are relatively rare in this dataset [25] Therefore, naive random splitting may lead to disjoint ensembles, for example, a training set lacking transition events while the test set contains them. To mitigate this distribution shift and decouple data sampling variance from optimization variance, we employed
split strategy. Specifically, we generated 200 candidate stratified splits and selected the partition that minimizes the 1- I-Wasserstein distance between the training and test set distributions of the radius of gyration, computed using the 10 Ca atoms. To ensure balanced initial conditions, we further constrained the split such that the 17 folded- inditt-ittttted and 17 unfolded-inition initiatee trajectories were each divided into 13 training and testing trajectories.
Eestribution-matchiin
After fixing the train-test split, all experiments performed with five independent model training runs using this fixed dataset. This allows us to evaluate numerical stability and convergence robustness of the learning algorithms.
We adopted the SchNet-based architecture and hyperparameters from [16. Therefore, we trained the models for 100 epochs with the Adam optimizer whose learning rate is 10-3 but doubled the batch size to 384 to maximize GPU memory utilization. We also used y = 0.01 for the regularization coefficient of DPNet-relaxed. The SchNet architecture used in our experiments is summarized in Table
Table 4: Architecture and hyperparameter configuration of the MchNetModel used in our experiments. The model processes atomic structures to produce fixed-dimensional feature vectors. Key hyperpa- rameters were adopted from [6] for consistency with prior work.
Component// Layer
Key Parameter
Value in Experiments
Input
Dictionary of atomic properties (positions, atomic numhers, etc.)
Pai PairiiseDistances
Cutoff radius (cutoff)
6.0 Á
GaussianRBBF
Number of radial basis functions (n_rbf)
20
SchNet Blocks
Feature dimension (n_aton_base basis)
Number of interaction blocks interactions)
3
Cutoff function
CosineCutoff
Output Layer (nn. Sequential)
Linear Layer
Output dimension (n_ www.fffettteeoeeee (oodee
15 (for k=16 modes with centering)
2. BatchNormid
Use batch normalization (use_batee batchnorn)
True
Output
Per-atom fcature tensor of shupe (Total Atoms, n_ feature _ooeess
Since we utilized batch normalization in the output layer, the learned feature vectors are constrained to have zero mean over the batch. This effectively enforces orthogonality to the constant cigenfunction. To explicitly account for this structure, we applied the centering strategy (Section 2.3) as default for all algorithms, modeling the constant mode separately. Consequently. the output dimension of the SchNet module set to k as we incorporate the fixed constant mode as an additional, separate feature.
15. This differs from the 16 dimensions used by Kostic et al. [16),
Implementation of VAMP-E Score. To evaluate the performance of the learned bases, we compute the VAMP-E score proposed by (48 evaluate the generalization capability of the learned Koopman operator. Unlike the VAMP-r scores which measure the canonical correlations within specific
28

---

dataset, the VAMP-E score estimates the approximation error of the learned model with respect to the true Koopman operator in the Hilbert-Schmidt norm 48
Specifically, maximizing the VAMP-E score is equivalent to minimizing the decomposition error IK-Kill on unseen data. The squared error can be expanded as |K- ||C]] - H- [|ilis Since |[К]] is unknown constant dependent only on the system, minimizing the error equivalent to maximizing the score Rg å 2(C,, K) - |[ Here, the first term measures the consistency between the learned model and the actual test dynamics, while the second term serves as a regularization term penalizing the model complexity on the test distribution.
In our implementation, we calculate the VAMP-E score using a cross- validation strategy where the model parameters are fixed based on the training data, and the score is evaluated on the test moments. Let Mg [f]. Mog" g [g] and Test [f g] denote the empirical second-moment moment matrices calculated from the test set. We first project these moments onto the singular subspace learned from the training data. Let Wu & pgrain [f])-1/²
1/2 and W1 - ppainin (g) -1/2 be the whitening matrices. Furthermore, let U, V and 2 be the singular vectors and the diagonal matrix of singular values, respectively, obtained from CCA with respect to the training data. The projected test moments are defined as:
Coo - U WTW.M P [r] = Σ 1/2 pi (+) - 1/2
Co1 - STW T test [f. g W/V V = 1.-1/2T,a test ( 2//2
C₁1 - V
W V = -E- -E-1/2M pios [ 1/22
Using the singular value matrix Σ computed from CCA with training data, the VAMP-E score is then defined as:
Rg - 2tr(2CCo1)
tr(2ICoo
2 tr (Treeet
Teeestt 10, x5) tr (M, [0] (b)).
(18)
=||Ai|
Intuitively, this formulation penalizes the model if the singular values Σ learned from training are inconsistent with the correlations observed in the test data (Co1), or if the learned basis functions lose orthonormality on the test set (Cao. C11 * I). Therefore, VAMP-E serves as a rigorous metric for validating both the learned singular subspaces and the estimated timescales [24].
On the Stabilization Effect of Nesting Techniques.
Test VAMP-2 Score (without Last-layer BatchNorm)
LoRA
6
LoRApa
Lo RoAAg
VAMP-2 Score
7
5
25
50
75
100
125
150
175
200
Training Epoch
Figure 5: Test VAMP-2 score on the chignolin dataset for LoRA variants trained without the final batch normalization layer. While the standard LoRA objective without nesting fails to converge, both nesting techniques LoRAjnt and LoRAgeg achieve high scores, suggesting improved optimization behavior in this challenging setting.
As an additional study, we report further benefit of nesting techniques during training. To create a more challenging learning setup, we remove the final batch normalization layer from the SchNet architecture and train for 200 epochs. Figure shows that the standard LoRA model fails to learn meaningful dynamics, as indicated by its low and stagnant VAMP-2 score. In contrast, both joint and sequential nesting exhibit stable training behavior and converge to high-scoring solutions. This
29

---

suggests that nesting techniques are not merely post hoc alignment mechanisms, but also provide stronger optimization signals this challenging setting. We leave a more in-depth analysis of this empirical observation for future work.
G.5
Computing Resources
All experiments conducted a server equipped with two Intel(R) Xeon(R) Gold 5220R CPUs, about 500 GiB of total RAM, four NVIDIA A5000 GPUs, and a 1TB NVMe SSD for storage. The chignolin molecular dynamics experiment was the most computationally intensive; each training run required 2-3 hours on two NVIDIA A5000 GPUs. Considering multiple runs across 5 data splits (Appendix G.4). the total computational effort for the reported chignolin experiments was equivalent to approximately 2-3 days. The other experiments (noisy logistic map, ordered MNIST, Langevin dynamics) were less demanding, completing in approximately hour, typically using one A,000 GPU, All experiments utilized PyTorch; the chignolin simulations additionally employed SchNetPack 2.1.1 [35). The total compute for the entire research project, including preliminary experiments, estimated at approximately 3 days of operational time on this hardware.
30

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We clarified the scope in the abstract and introduction. We do not overstate more than the paper.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Yes, discuss the limitation in the concluding remarks.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model sell-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't It acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
31

---

Justification:
Yes, we have the theoretical results and provide full proof in the appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: Yes)
Justification:
We described full details of experiments.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reprodueibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
32

---

Answer: [Yes]
Justification:
We release our code implementation.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.cc/ for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
//nps..c//ublii/giideee
http:/Sbbbmsssooooouooooo sionPollcy) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper speciffy all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer:[Yor Yes]
Justification: Yes, we are opening full training and test details.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: Yes, we ran the experiments with various training seeds.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, at least for the experiments that support main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
33

---

It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of enrors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Answer: [Yes]
Justification:
We included these details in the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g.. preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer:
[Yes]
Justification: Yes, definitely follow the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is a special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to
34

---

generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
there are negative socictal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators original owners of assets (e-g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: We did this for the chignolin experiment.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
assets are released, the license, copyright information, and terms of use the package should be provided. For popular datasets, pwpprswi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
this information is not available online, the authors are encouraged to reach out to the asset's is creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
35

---

Answer: [Yes]
Justification:
We release our code implementation.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
We are not contatining crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
We do not involve crowdsourcing research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Answer: [Yes]
36

---

Justification:
We have not used LLM for our core research components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy htte::///eerris..cccoooeeeeerrrr for what should or should not be described.
37


Reviews:
[{"review_id": "Jodkaa7wYC", "reviewer": "Reviewer_jeGy", "summary": "This paper introduces a scalable and stable method for learning the top-k singular subspaces of the Koopman operator in stochastic dynamical systems using LoRA. Unlike prior methods such as VAMPnet and DPNet, which rely on backpropagation through numerically unstable operations like SVD or matrix inversion, LoRA formulates a variational objective entirely in terms of polynomial forms of moment matrices, enabling unbiased gradient estimation and efficient training using minibatches. The proposed method is validated on multiple benchmark systems, including synthetic chaotic maps, Langevin dynamics, molecular simulations (Chignolin), and Ordered MNIST, demonstrating superior performance in general.", "questions": "N/A", "limitations": "N/A", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\n- The paper is clearly structured, with a coherent flow from motivation, problem formulation, method, theory, and experiments. \n- I’m not much familiar with the field, but the proposed approach’s motivation is convincing. e.g., challenges on VAMPnet and DPnet (e.g., unstable backpropagation, scalability) are clearly discussed.\n- Notation is mostly consistent, and figures are helpful for understanding convergence and performance trends.\n\nWeaknesses and Questions\n- While LoRA seems not been applied to dynamic mode decomposition, it looks like a simple integration of existing LoRA method to existing problem to a certain extent. \n- While the theoretical background and motivation are sound, some of the proof sketches and intuition for why LoRA performs better (especially under non-normal dynamics) could be clarified in the main text. Please correct me if it was presented in section 3.\n- I agree that VAMPnet and DPnet would be more difficult to optimize. However, I'm not sure the experimental results effectively exhibit such limitations. In terms of rel. error in Figure 1, LoRA achieves the lowest error comparing others, but it’s not always true when evaluating Hausdorff (the lowest distance achieved by DPNet).\n- Nesting (seq or jnt) results do not seem consistent. For example, LoRA_jnt is superior in Figure 2, but LoRA_seq is better in Table 1. \n- Also, I expected results showing LoRA is significantly faster and scalable than VAMPnet and DPnet, but such experiments seem not considered.\n\nAs I’m only little familiar with the field, I’ll have to see other reviewers score to finalize mine.", "was_revised": true, "final_justification": "Due to the abovementioned reasons, my final rating is 4: borderline accept.", "rebuttal": "We appreciate the reviewer for their time and effort in reviewing our manuscript, providing thoughtful comments. We will address the reviewer's questions individually below. We will update our manuscript to improve the clarity, accompanying the feedback.\n\n---\n\n1. **On the Novelty of Our Work**\n\nWe emphasize that our key contribution is to raise a key practical bottleneck in the dominant optimization frameworks for Koopman operator learning, and address it by **reconceptualizing the underlying variational objective**.\n\nAs we detailed in the Item 1 of our response to Reviewer bFvN, prior methods like VAMPnet/DPNet, while theoretically sound, are built on objectives that are fundamentally incompatible with modern deep learning optimization due to their reliance on numerically unstable and statistically unfavorable operations such as SVD and inversion of empirical moment matrices. Our work is the first to identify this practical failure mode as a central problem and propose a principled solution via a direct low-rank approximation.\n\n2. **On the Intuition Behind LoRA’s Superiority**\n\nOur method’s core advantage stems from its numerical stability, a crucial aspect that has been largely overlooked in the field of Koopman operator learning. As discussed in the first item, prior methods, including the recent VAMPnet and DPNet, suffer from operations like matrix inversion on empirical moment matrices. They may result in not only numerical instability for ill-conditioned empirical moments, but also biased gradient estimates, which may ultimately hinder their application to complex, large-scale dynamical systems.\n\nOur LoRA approach directly addresses this fundamental problem by reformulating the underlying variational principle of Koopman learning. Instead of relying on objectives that require unstable linear algebraic operations, we propose to directly minimize the operator norm of **low-rank approximation error** of the Koopman operator. \n\nThe resulting objective function consists solely of polynomials of second-moment matrices, thereby eliminating the need for any unstable linear algebraic operations during training. Furthermore, we have theoretically validated the practicality of this objective function by proving its statistical consistency under finite-minibatch settings, as detailed in Appendix E.2. This combination of a stable formulation and theoretical guarantees provides the intuition for LoRA's superior performance and scalability.\n\n3. **On Experimental Results: Stability vs. Performance Metrics**\n\nThe primary goal of the experiments is to demonstrate that LoRA offers superior performance and scalability as promised. As the reviewer noted, in the Logistic map experiment, DPNet-relaxed shows superior performance in terms of the Hausdorff distance for estimated eigenvalues. This suggests that DPNet-type objectives may work as expected for small-scale systems. \n\nOur subsequent experiments reveal, however, that the DPNet(-relaxed) technique suffers to scale to high-dimensional setups. In the ordered MNIST experiment, where the environment is higher-dimensional, all three LoRA methods exhibit far better prediction compared to the baselines, including variants of DPNet. More critically, the continuous-time version of the DPNet objective also fails to accurately capture eigenfunctions of 1D Langevin dynamics. Ultimately, in the Chignolin experiments, VAMPnet and DPNet completely fail to converge, and DPNet-relaxed only achieves a far worse data fit compared to LoRA. In the item 1 of our response above, we further strengthen this scalability argument by showing that DPNet-relaxed becomes to break down (i.e., diverge or significantly underperform) on settings with larger models and smaller batch sizes.\n\n4. **On Nesting Techniques**\n\nWe acknowledge that the reviewer rightfully asked about a winning strategy among the two nesting techniques. As alluded to above (see item 1), we used “joint parameterization” throughout the paper. \n\n* On the one hand, the joint nesting is a principled technique that may guarantee to recover the desired solutions by finding the global optima regardless of the form of parameterization, as alluded to in line 200\\.  \n* On the other hand, sequential nesting may not guarantee convergence for the case of joint parameterization, as it assumes independent parameter updates for different modes.\n\nIn this paper, however, we tested and presented the performance of sequential nesting in this scenario as “a working heuristic” that can occasionally outperform joint nesting (as shown in Table 1). Therefore, as suggested in \\[29\\], we recommend using joint nesting in general for jointly parameterized eigenfunctions, but practitioners may consider applying sequential nesting as a heuristic with some caution of no guaranteed convergence. For the case of separate parameterization where each mode is parameterized by a separate network, sequential nesting may be preferred, as argued in \\[29\\]. We will carefully clarify these points in the manuscript. \n\nOverall, we remark that our main claim about nesting is that nesting can help learning more than without it, leading to better performance in general (as in Figure 2). To provide direct evidence for this claim, we conducted a new, more challenging experiment on the Chignolin dataset. We deliberately made the optimization more difficult by removing the last batch normalization layer from the SchNet architecture and trained for 200 epochs. This setup tests the intrinsic stability of the gradient updates. The plots below, visualizing VAMP-2 score on test data, present a clear benefit of nesting. While the vanilla LoRA (without nesting) struggles to learn meaningful dynamics, in sharp contrast, both nesting (both joint and sequential) enable the model to converge to a high-quality solution, demonstrating that nesting acts as a critical regularizer and stabilizer.\n\n```\n  10 +----------------------------------------+\n   9 |                                        |\n   8 |                           LoRA ******* |\n   7 |                                        |\n   6 |                                        |\n   4 |                                        |\n   3 |****                                    |\n   2 |    ************************************|\n   1 +----------------------------------------+\n     0   20  40  60  80  100 120 140 160 180 200\n  10 +----------------------------------------+\n   9 |                    ********************|\n   8 |             *******    JntLoRA ******* |\n   7 |           **                           |\n   6 |          *                             |\n   4 |         *                              |\n   3 |   *  ***                               |\n   2 |*** **                                  |\n   1 +----------------------------------------+\n     0   20  40  60  80  100 120 140 160 180 200\n  10 +----------------------------------------+\n   9 |                                        |\n   8 |                        SeqLoRA ******* |\n   7 |                 ***********************|\n   6 |               **                       |\n   4 |              *                         |\n   3 |            **                          |\n   2 |************                            |\n   1 +----------------------------------------+\n     0   20  40  60  80  100 120 140 160 180 200\n\t\t\tEpochs\n```\n\n5\\. **On Stability and Scalability Compared to VAMPnet and DPNet**  \nWe agree that the advantages (stability and scalability) of LoRA can be better highlighted in our experiments. \n\n* The latter two experiments (Langevin dynamics and Chignoilin) more directly demonstrate the stability and scalability of LoRA. In the Langevin dynamics example, we were not able to train a model using the DPNet objective. In the Chignolin experiment, VAMPnet and DPNet objectives fail to converge, as shown in Table 2\\. The only working example is DPNet-relaxed, and we demonstrated how LoRA can lead to a better fit of the neural networks to the trajectory data.  \n* To further highlight LoRA's superior scalability and stability against DPNet-relaxed,  we conducted additional experiments on Chignolin dynamics (data split for seed $0$). We compared the test VAMP-E score on our LoRA variants against DPNet-relaxed by increasing the number of modes (**N**) and the SchNet node feature dimension (**H**), while also reducing the batch size (**B**). \n\n    The DPNet-relaxed baseline, previously considered a stable alternative, **completely failed to train under these conditions**, as indicated by the negative VAMP-E scores. In stark contrast, LoRA and its variants learned reliably and remained stable across all these challenging settings. This new experiment corroborates that LoRA solves the fundamental issues of numerical instability and resulting poor scalability in the prior work.\n\n    |(B,N,H)|SeqLoRA|JntLoRA|LoRA|DPNet-relaxed|DPNet-relaxed-Ctr|  \n    |:-:|:-:|:-:|:-:|:-:|:-:|  \n    |(384, 64, 128)|31.61|29.89|29.51|-9.76${}^2$|5.20${}^1$|  \n    |(96, 64, 64)|25.83|23.82|25.51|4.75${}^1$|5.93${}^1$|  \n    |(384, 64, 64)|26.66|27.25|24.11|-1.97${}^2$|4.72${}^1$|  \n    |(384, 32, 128)|17.68|15.66|16.26|1.75|3.08|  \n    |(96, 32, 64)|17.19|15.64|17.36|4.03${}^1$|6.44|  \n    |(384, 32, 64)|16.69|17.02|15.79|4.36|4.91|  \n    |(384, 16, 64)|9.73|8.79|8.80|3.05|4.19|  \n    |(96, 16, 64)|10.22|8.80|8.52|2.74|4.40|\n\n    (VAMP-E scores are calculated on overall test data, not distinguishing folded/unfolded initialization.)  \n    ${}^1$ During training, the VAMP-2 score exceeded the theoretical maximum (number of modes $N$). This is a strong indicator of numerical instability, suggesting a severe breakdown of orthogonality or a failure to learn diverse features.\\\\  \n    ${}^2$ The negative VAMP-E score implies that the approximation error $\\\\|\\\\mathcal{K} \\- \\\\hat{\\\\mathcal{K}}\\\\|$ is larger than the (Hilbert-Schmidt) norm of the true operator $\\\\|\\\\mathcal{K}\\\\|$ itself, rendering the approximation meaningless."}, {"review_id": "UQVAbYfyCL", "reviewer": "Reviewer_bFvN", "summary": "This work studies a deep-learning–based framework for modeling dynamical systems by approximating the Koopman operator through a low-rank singular-value decomposition. Its principal contributions are evaluation of the Low-Rank (LoRa) contrastive loss in the context of the (stochastic) Koopman operator approximation, systematic comparisons with existing related methods, careful elaboration  of implementation details, and  and extensive experiments.\n\nThe chief advantage of LoRa is that its computation avoids the inversion of latent-space covariance matrices—objects that are often ill-conditioned during training. By avoiding these unstable inverses, LoRa produces unbiased estimates of the loss and its gradients, delivering more reliable optimization and superior performance compared with previous state-of-the-art losses.", "questions": "- **Figure 1 (top row):** Why does the singular-value error *increase* as the number of singular functions grows? Intuitively, adding more modes should *reduce* the approximation error. Please clarify the underlying reason.\n- **Left vs. right singular spaces:** The paper compares the span of the left-right-singular-functions on downstream inference tasks, yet theory predicts that both should converge to the same Koopman-invariant subspace. What practical insight is gained from this comparison?.\n- **Scope of the empirical study:** The work largely integrates results from prior literature and promises a thorough ablation of LoRa-based Koopman approximation. In practice, the only clear conclusion is LoRa’s superiority over DPNet and VampNet losses. To strengthen the empirical contribution, consider adding experiments that yield actionable insights on:\n    - **Latent dimensionality:** How do all models scale as the latent space grows? We know that large latent spaces will lead to unstable DPNet and Vamp scores. Is LoRa based optimization unaffected by this?\n    - **Orthonormal regularization:** Is orthonormal regularization influencial for numerical stability and/or convergence speed of the training process?\n    - **Batch-size sensitivity:** How does batch size affect optimization stability of these methods?", "limitations": "Section limitation is present, but the discussion is quite superficial and not discussing limitations per se, but rather plans for future work.", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 1, "strengths_and_weaknesses": "### __Strengths__\n\n- The paper provides a solid review of related deep-learning approaches for stochastic Koopman operators and provides their comparisons.\n- The paper is self contained, well written and accessible.\n- The variants of the proposed algorithm appear to consistently beat state-of-the-art baselines.\n\n### __Weakenesess__\n\n- __Novelty:__ Although the paper cites [29] and [15], Section 3.1 should explicitly acknowledge that the LoRa loss was *first* introduced for low-rank conditional-expectation operator approximation in those works, while paper [14] proposed truncated SVD as the approach for representation learning for stochastic processes estimation. Additionally, I strongly suggest to place LoRa in the wider context of self-supervised contrastive learning and density ratio fitting by citing \n    \n    • Wang *et al.* (2022), “Spectral Representation Learning for Conditional Moment Models”\n    \n    • Sugiyama *et al.* (2012), “Density Ratio Estimation in Machine Learning”\n    \n- __Inconclusive empirical results:__ While the paper promises an experimental evaluation of design choices LoRa based low-rank approximation of the Koopman operator, many results appear inconclusive:\n    - **Left vs. right singular functions**\n    Figure 1 shows no clear performance gap, and Figure 2 reports no statistically significant difference. These results could be summarized in prose, omitting the corresponding curves to reduce visual clutter.\n        \n        The claim\n        \n        > “...we empirically found that using the right singular basis can sometimes yield better results.”\n        > \n        \n        is not convincingly supported by the current data.\n        \n        Since theory predicts both subspaces converge to the same Koopman-invariant space once training stabilizes, this design choice may not warrant prominent discussion; consider moving it to the appendix.\n        \n    - **Joint vs. sequential nesting**\n    The implementation of sequential nesting is unclear, given that all left/right functions are parameterized by the same encoder / delayed-encoder. Please clarify the procedure. Furthermore in the joint nesting, it is unclear how the weights are chosen given that the during learning the left-and-right data representations are not ordered by relevance of the singular space. It seems that if the nesting technique is going to be one of the main contributions of the paper, it warrants both better introduction and more conclusive results.\n\n    - **Hyper-parameter study**\n        In the synthetic experiment, Figure 1 shows no discernible trend for learning rate or batch size. It is hard for the reader to draw any conclusion from these plots, beyond the better performance of LoRa against baselines. Either try to pursue a deeper ablation study with more conclusive results and trends, or simplify this plot to avoid the clutter.  \n  \n- __Minor issues:__  Font sizes in most plots are too small, compromising readability. Use the available whitespace to enlarge labels and legends. The x-axis of Figure 3 is unlabeled; specify the quantity and units. The paper does not report the scalar weighting the DPNet regularization term, nor whether this parameter was tuned.\n\n\n__In summary, I believe that the paper is solid contribution with interesting empirical results, but my concern lies in clear and sufficient novelty w.r.t. the existing literature.__", "was_revised": true, "final_justification": "While I do agree with the authors that this work is meaningful contribution to the community, and I believe that the paper is well written, I remain doubtful about clear acceptance of the presented contribution.  \n\nAs I noted in the discussion with the authors. Previous works already demonstrated the role of the SVD of the transfer operator, as a specific case of the conditional expectation operator (CEO), in reliably learning dynamical systems, and demonstrated how one learns spectral features of CEO via MSE loss on its kernel. So, from my perspective the novelty just remains in **empirically testing existing approaches** on **already proposed  datasets** from dynamical systems. I don't find this strong enough contribution reflecting the score 5.", "rebuttal": "We appreciate the reviewer’s effort in reviewing our manuscript and providing constructive feedback. \nBelow, we address each concern separately.   \n(Note: Just to clarify our convention, we intended it to be LoR**A** for **Lo**w-**R**ank **A**pproximation.)\n\n---\n\n1. **On Novelty/Originality**  \n\nWe thank the reviewer’s comments on the literature (on the truncated SVD view in \\[14\\] and on self-supervised contrastive learning) and will incorporate these in our revision. \n\nHowever, we wish to remark that the idea of LoRA dates back to as early as Schmidt’s paper in 1907, which proposed it to characterize the top-$K$ eigenfunctions of integral operators. Given this, we emphasize that our key contribution is, in the field of data-driven learning of Koopman operators: (1) we first **identify** common, crucial flaws in the dominant optimization frameworks, namely VAMPnet and DPNet, and resolve them by **applying** the LoRA framework and concretely **demonstrating** its benefit by experiments.\n\nIn this work, identifying the challenges in practical optimization, we show that LoRA, albeit a classical idea, can address the issues in the prior work in a principled yet simple way.  This stability is what allows us to reliably learn superior representations in high-dimensional settings, as demonstrated by our results. This is why, even on smaller benchmarks, our method consistently finds better solutions (e.g., higher VAMP-2/E scores and more accurate timescales in Chignolin), as it is hindered by neither the numerical instability nor gradient biased estimates during optimization that plague its predecessors. We will revise Sec 3.1 to explicitly state this distinction and place our work in the context of the suggested literature in a broader context.\n\nWe would like to ask the reviewer to kindly reassess the “originality” of our work given this context.\n\n2. **On the Scope of the Empirical Study: Scalability and Orthogonality**\n\nTo provide more direct empirical evidence, we conducted new scalability and stress-test experiments on the Chignolin experiments, with data split for seed $0$. We compared the test VAMP-E score on our LoRA variants against DPNet-relaxed by increasing the number of modes (**N**) and the SchNet node feature dimension (**H**), while also reducing the batch size (**B**):\n\n|(B,N,H)|SeqLoRA|JntLoRA|LoRA|DPNet-relaxed|DPNet-relaxed-Ctr|  \n|:-:|:-:|:-:|:-:|:-:|:-:|  \n|(384, 64, 128)|31.61|29.89|29.51|-9.76${}^2$|5.20${}^1$|  \n|(96, 64, 64)|25.83|23.82|25.51|4.75${}^1$|5.93${}^1$|  \n|(384, 64, 64)|26.66|27.25|24.11|-1.97${}^2$|4.72${}^1$|  \n|(384, 32, 128)|17.68|15.66|16.26|1.75|3.08|  \n|(96, 32, 64)|17.19|15.64|17.36|4.03${}^1$|6.44|  \n|(384, 32, 64)|16.69|17.02|15.79|4.36|4.91|  \n|(384, 16, 64)|9.73|8.79|8.80|3.05|4.19|  \n|(96, 16, 64)|10.22|8.80|8.52|2.74|4.40|\n\n(VAMP-E scores are calculated on overall test data, not distinguishing folded/unfolded initialization.)  \n${}^1$ During training, the VAMP-2 score exceeded the theoretical maximum (number of modes $N$). This strongly suggests numerical instability, suggesting a severe breakdown of orthogonality or a failure to learn diverse features.  \n${}^2$ The negative VAMP-E score implies that the approximation error $\\\\|\\\\mathcal{K} \\- \\\\hat{\\\\mathcal{K}}\\\\|$ is larger than the (Hilbert-Schmidt) norm of the true operator $\\\\|\\\\mathcal{K}\\\\|$ itself, implying a poor approximation of the operator.\n\nOur experiments reveal that DPNet-relaxed completely fails when scaled to larger models or trained with smaller batch sizes. In direct contrast, our LoRA-based methods demonstrate stability across all settings, even at 4 times more modes. This experiment offers strong empirical proof that LoRA exhibits much more favorable scalability than the prior works.\n\n**Orthonormal Regularization**\n\nWe answer the reviewer's question about the role of orthonormal regularization. Our LoRA method does not require explicit regularization to achieve stable convergence, which is a key advantage. The necessity of such regularization in prior methods stems from two interconnected issues:\n\n* **Theoretical Optimality**: As established by the Eckart-Young theorem, the best rank-K approximation of the Koopman operator is achieved by its top-K orthonormal singular functions. Therefore, learning an orthonormal basis is not just a heuristic but a core requirement for finding the optimal low-rank model.  \n* **Practical Stability**: As the reviewer insightfully pointed out, maintaining this orthonormality is crucial for the stability of prior methods like DPNet. Their objectives rely on numerically unstable operations (e.g., matrix inversion). If the learned basis functions lose orthogonality and become collinear during training, these operations fail, leading to the training collapse we demonstrated in the Chignolin experiment.\n\nOur LoRA objective, by its design, sidesteps this issue. To provide a more direct evidence, we analyzed the orthonormality of the raw learned basis functions (before any post-hoc CCA) on the Chignolin test data (seed $0$). We measured the Frobenius distance between the normalized Gram matrix and the identity matrix, where the maximum is about $15.5$:\n\n|(B, N, H)|SeqLoRA|JntLoRA|LoRA|DPNet-relaxed|DPNet-relaxed-Ctr|  \n|:---:|:---:|:---:|:---:|:---:|:---:|  \n|(384, 16, 64)|0.38 / 0.53|0.30 / 0.26|0.21 / 0.25|6.13 / 6.17|5.37 / 5.32|  \n|(96, 16, 64)|0.50 / 0.48|0.34 / 0.35|0.24 / 0.25|4.31 / 4.22|5.00 / 5.16|\n\nThe results reveal an order-of-magnitude difference, showing that our LoRA variants naturally learn a nearly orthonormal basis directly from the optimization, while DPNet-relaxed fails to do so. This demonstrates that LoRA's stability is not just about preventing training failure but also about guiding the model toward a well-structured, theoretically optimal solution without explicit regularization.\n\n3. **On the Clarity and Implementation of Nesting Techniques**\n\nWe agree that some clarifications are warranted for nesting techniques, as we elaborate below. We will carefully revise our manuscript to strengthen our contribution.\n\n**Clarification on Sequential Nesting with Shared Encoders**   \nAs pointed out by the reviewer, we used “joint parameterization” for estimating singular functions throughout the paper. \n\n- The joint nesting is a principled technique that may guarantee to recover the desired solutions by finding the global optima regardless of the form of parameterization, as alluded to in line 200.  \n- On the other hand, sequential nesting may not guarantee convergence for the case of joint parameterization, as it assumes independent parameter updates for different modes.\n\nIn this paper, however, we tested and presented the performance of sequential nesting in this scenario as a working heuristic that can occasionally outperform joint nesting (Table 1). Therefore, similar to suggested in \\[29\\], we recommend using joint nesting in general for jointly parameterized eigenfunctions, but practitioners may consider applying sequential nesting as a heuristic with some caution of no guaranteed convergence. For the case of “separate parameterization” where each mode is parameterized by a separate network, sequential nesting may be preferred, as argued in \\[29\\]. We will carefully clarify these points in the manuscript. \n\nThe detailed implementations for both versions of nesting are currently described in Appendix F. For sequential nesting, we achieve a sequential learning process by controlling the flow of the gradients. When training the $i$-th function pair ($f\\_i,g\\_i$), we use `torch.no_grad()` (referred to as `sg` in the pseudocode in Appendix F) to stop gradients from flowing back to the parameters that define the previously learned functions ($\\\\mathbf{f}\\_{1:i-1}, \\\\mathbf{g}\\_{1:i-1}$). This procedure effectively \"freezes\" the previously optimized subspaces, ensuring that the network focuses its updates only on learning the new, orthogonal mode. In our revision, we will make this implementation detail explicit in Section 3.1.\n\n**Weight Selection for Joint Nesting**  \nIn our experiments, we used uniform weights ($\\\\alpha\\_i \\= 1/k$ for all $i$), as we did not assume any prior knowledge about the relative importance of different singular modes. We will add this detail to the experimental setup section for full reproducibility.\n\n4. **Responses to Specific Questions**\n\n* **Increase in Singular-Value Error**: As explained in line 284, the error in Figure 1 is the *relative* error, $(\\\\hat{\\\\sigma}\\_i^2 \\- \\\\sigma\\_i^2)/\\\\sigma\\_i^2$ for each index $i$. As the index $i$ increases, the true singular value $\\\\sigma\\_i$ becomes smaller and less dominant. This causes the relative error to increase, even if the absolute error decreases. We will clarify this in the caption.  \n* **Left vs. Right Singular Spaces**: We appreciate the reviewer's suggestion. We intended to provide a complete analysis for practitioners and show that either basis can be used for EDMD with similar performance, given that the literature only considered using the left eigenbasis. However, we agree this point can be de-emphasized to improve clarity. We will move the detailed comparison to the appendix in our revision.  \n* **On Minor Issues and DPNet Regularization**: We agree that the plots can be revised to improve readability and will correct all visual elements (font sizes, axis labels) in the revised manuscript. Regarding the DPNet regularization parameter $\\\\gamma$, we followed the experimental setup of the DPNet paper \\[14\\] to ensure a fair and direct comparison with the reported results in the literature. We will clarify this in our experimental setup section.\n---\nWe hope our responses and planned revisions address the reviewer's concerns. If these clarifications adequately address the reviewer’s concerns, we would greatly appreciate it if the reviewer could consider updating the score to reflect what we believe are meaningful contributions to the community."}, {"review_id": "7oPNIDPNwt", "reviewer": "Reviewer_ZVY9", "summary": "This manuscript presents a novel method, termed LoRA, for learning the top-k singular functions of the Koopman operator for stochastic dynamical systems. The authors identify a critical issue in existing deep learning-based approaches like VAMPnet and DPNet: their reliance on numerically unstable operations such as matrix inversion and singular value decomposition within the training loop, which can lead to biased gradients and optimization challenges. To address this, the proposed LoRA method is based on directly minimizing the Hilbert-Schmidt norm of the low-rank approximation error. This formulation results in a loss function that is a simple polynomial of empirical second-moment matrices, thereby avoiding unstable computations during optimization. The method's efficacy and stability are demonstrated through a series of experiments on synthetic and molecular dynamics datasets, where it is shown to outperform baseline methods.", "questions": "Based on the analysis, we pose the following questions to the authors. First, in Theorem C.3 of Appendix C.2, there is a clear discrepancy. The theorem states that for a left eigenfunction $\\zeta$, the relationship is $\\mathcal{K}^{*}\\zeta=\\lambda\\zeta$. However, the final step of the proof on line 832 derives the result as $\\overline{\\lambda}\\zeta(x')$, where $\\overline{\\lambda}$ is the complex conjugate of $\\lambda$. This conclusion from the proof, $\\mathcal{K}^{*}\\zeta=\\overline{\\lambda}\\zeta$, is mathematically correct for adjoint operators but contradicts the theorem's statement. We ask the authors to clarify and correct this imprecision.\n\nSecond, the authors acknowledge in Appendix G.4.4 that the short lag time used in the Chignolin analysis \"inherently biases the model towards resolving faster dynamic events.\" This raises a critical question about the method's general applicability: for a new system, how should a practitioner select an appropriate lag time to ensure the discovery of the most critical, slow dynamics, rather than being misled by parameter choices into focusing on faster, less relevant modes? Does this imply that the successful application of LoRA requires extensive prior knowledge or tedious sensitivity analyses?\n\nThird, while the stability of LoRA is well-demonstrated, the experimental validation relies significantly on systems that could be considered \"toy\" or methodologically \"tricky.\" The Noisy Logistic Map, for instance, is a 1D system chosen specifically because its finite-rank kernel is analytically tractable. Similarly, the Ordered MNIST experiment uses a synthetic process with simple, deterministic dynamics layered on top of image data. While valuable for controlled comparisons, how confident are the authors that the strong performance on these simplified or contrived benchmarks will translate to more complex, high-dimensional, and chaotic real-world systems where the underlying dynamics are not low-rank or neatly structured?", "limitations": "The most fundamental limitation of this research is that, as a data-driven method, the quality of the learned dynamical model is entirely beholden to the quality and characteristics of the input data. The Chignolin experiment serves as a primary example: if the available data, due to factors like trajectory length or sampling frequency, does not sufficiently contain the signal of the slowest physical processes (e.g., microsecond-scale global folding), the method cannot recover them, regardless of its mathematical optimality or algorithmic stability. Therefore, the method cannot guarantee the discovery of the absolute slowest physical timescale in all cases; rather, it finds the most prominent observable slow process within the given data. This data-dependency defines the principal boundary for the method's applicability when applied to new systems.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The primary strength of this work lies in its clear identification of a key technical bottleneck in existing Koopman analysis methods—namely, numerical instability during optimization—and its proposal of a mathematically rigorous and algorithmically elegant solution. The LoRA objective function is well-motivated, and its stable polynomial form represents a significant advantage, facilitating straightforward implementation and optimization within modern deep learning frameworks. The derivation provided in the appendix confirms its theoretical soundness. The experimental results, particularly the successful application to the complex Chignolin molecular dynamics system where baseline methods reportedly failed to converge, compellingly demonstrate the superior stability and practical utility of the proposed algorithm.\n\nHowever, the rigor and persuasiveness of the work are undermined by several significant flaws in its reporting and execution. Firstly, the reporting on the crucial Chignolin experiment contains multiple errors in the main manuscript, including incorrect units, data splits, and table labels, which, while noted in an appendix errata, reflect a lack of rigor in the report itself. Secondly, the work’s central conclusion relies heavily on a post-hoc physical justification in Appendix G.4.3 to bridge the gap between their experimental data and a reference value from a different temperature; this justification is qualitative and lacks quantitative support, making it appear speculative. Most critically, the authors admit in Appendix G.4.4 that their choice of a short lag time \"inherently biases the model towards resolving faster dynamic events,\" which is a concession of a systemic flaw in the experimental design that makes it unsuitable for validating the theoretical goal of finding the system's slowest physical processes.", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thorough review, which has helped us strengthen our manuscript. We address the reviewer's concerns individually below. \n\n---\n\n1. **On the Chignolin Experiment: Lag Time, Reference Timescale, and LoRA's Robustness**\n\n    1. `On the Choice of Lag Time`: The reviewer points out that `the authors acknowledge in Appendix G.4.4 that the short lag time used in the Chignolin analysis \"inherently biases the model towards resolving faster dynamic events.\"`, and raises a valid concern regarding how lag time should be chosen and whether this necessitates sensitivity analysis or extensive prior knowledge. We apologize for any confusion this statement may have caused, and we would like to clarify its intended meaning below.\n\n        In an MD analysis setup, there exists two important design parameters:   \n         - **(1) Lag time**, which is the resolution to emulate/discretize the continuous time dynamics. If the lag time gets longer, the discretized trajectories will only contain very low-frequency dynamics. As the lag time gets shorter, the trajectories will capture more refined, high-frequency information of the underlying dynamics, but it may capture more stochastic noise which thus could make the analysis more difficult.  \n         - **(2) Trajectory length**, which is the length of the simulated trajectory. As a trajectory gets longer, the experiment setup aims to capture longer dynamics, and vice versa.\n\n        Our experiment analyzes rather short-length trajectories in a short lag time regime (100 ps).   \n        Some prior works in the literature analyze the trajectories of sweeping the lag time for a complete understanding of the data, see, e.g., \\[Fig. 6, 22\\]. In the current work, we aim to demonstrate whether the considered methods can reliably fit to the high-dimensional trajectories in the challenging short lag time regime.\n\n        In this setup, we show that **LoRA can reliably extract the dominant, slow dynamics from high-dimensional trajectories, even in scenarios where DPNet-Relaxed only achieves strictly subpar performance in terms of the both training/test VAMP-2 score, while all other methods completely fail to converge**.\n\n    2. `Regarding the Timescale Analysis`  \n        The reviewer raises a concern on the validity of the timescale analysis of the Chignolin experiment. We agree with the reviewer’s sentiment on the issue, and here we clarify the purpose and validity of such analysis. First of all, we emphasize that the concerned molecular dynamics is normal, and thus accurate estimation of singular functions with a high test VAMP-2 score directly implies accurate learning of the dynamics. Therefore, the subsequently high VAMP-2 score of LoRA variants against the baseline methods should be the main experimental takeaway of the experiment. Given this, the timescale analysis, a standard post-hoc analysis in this field \\[19,22\\], is to provide a physical relevance of the answer from the computational method. Our close estimate to the known reference value mildly suggests that our computational result is on the track, but this does not indicate any scientific conclusion. Moreover, as the reviewer rightfully noted, we have some additional discussions in the appendix to further compensate for the different configurations from the reference.We acknowledge that this physical discussion may distract our main message. We thus will tone down the current emphasis on the timescale analysis, but rather emphasize the good fit of LoRA to the underlying data. \n\n2. **On Reporting Rigor and Theorem C.3**\n\n    We sincerely thank the reviewer for their meticulous reading.\n\n    * **Reporting Errors**: We agree that the errors in the main text (e.g., units, data splits) were oversights. We have corrected all of them to ensure the manuscript is accurate.  \n    * **Regarding Theorem C.3**: After careful reexamination, we have found no inconsistency between the statement and proof in Theorem C.3. Namely, we claim that $\\\\mathcal{K}^{\\*}\\\\zeta=\\\\overline{\\\\lambda}\\\\zeta$ in the statement, and prove it subsequently. If there is anything we have missed, please let us know. \n\n3. **On Performance in More Complex, Real-World Systems**\n\n    - While our benchmarks were chosen for their controlled properties, they were designed to test the very failure modes we aimed to solve. Our experiments validate the core advantage of our method that are crucial for real-world application with complex dynamics. The **Noisy Logistic Map** and **Ordered MNIST** were selected because they exhibit non-normal or structured dynamics where optimization stability is critical. The Chignolin system, while smaller than real-world proteins, represents a real and challenging molecular dynamics problem.\n\n    - Regarding its applicability in more challenging real-world setups, we believe that the normal dynamics can be relatively easily handled by the LoRA framework even if it is not low-rank, as we can increase the number of modes to fully capture the dynamics. We conjecture that dealing with non-normal dynamics such as chaotic systems will remain as great challenges, as accurate estimation of non-normal dynamics in the sense of accurate top-K SVD may fail to capture some dominant modes in its eigen analysis. We note that this is an inherent limitation of available computational tools, not limited to the LoRA framework. That said, since we do not know how to directly compute eigenfunctions of non-normal dynamics, the top-K SVD framework will remain as an alternative solution that provides a systematic learning of dynamics, which could still provide some insights.\n\n4. **Additional Experiments**\n\n    As we explained in our response to **Reviewer bFvN**, we conducted additional ablation study in the Chignolin experiment, by varying model dimensionalities and batch sizes to test the scalability and robustness of LoRA. The results  corroborate that LoRA's advantage translates to more complex, higher-dimensional settings. While further validation on larger systems is an important next step, we are confident that the fundamental stability we have demonstrated may allow the data-driven, deep-learning-based approach to be more widely applicable."}, {"review_id": "jNxEv2QJw4", "reviewer": "Reviewer_acob", "summary": "The authors propose a row rank approximation method to approximate Koopman operators. The simplicity of the algorithm makes it easy for the method to integrate with modern deep learning pipelInes by eliminating the need for unstable linear\nalgebraic operations.", "questions": "The singular vectors should be orthogonal each other. How do you guarantee the orthogonality of the vectors $f_1,\\ldots, f_k$?", "limitations": "Investigating what type of neural networks can be used for constructing $\\mathbf{f}$ and $\\mathbf{g}$ would be interesting.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Finding good low-dimensional space for representing Koopman operators is a nonobvious and important problem. This paper try to tackle this problem, and the topic is related to the community.\n\nThere are existing methods that aim to learn low-dimensional spaces for representing Koopman operators. The advantage of learning a singular space is not clear for me. Could you clarify the advantage of learning a singular space instead of learning general linear spaces for low-dimensional approximation of Koopman operators?\n\nMinor comments:  \n- Line 49: what do you mean by the identity $x_0=x_0$?\n- Line 102: it seems that the limit is taken in the topology of the function space where the Koopman operator is defined, not the strong operator topology for operators.\n- Line 115: \"for even for \" should be \"even for\"?", "was_revised": true, "final_justification": "The authors addressed my comments properly. Thus, I updated my score.", "rebuttal": "We appreciate the effort and time in reviewing our manuscript. We address the reviewer's questions individually as follows. \n\n---\n\n1. **On the Advantage of Learning Singular Spaces**  \nThe reviewer asks why it is preferable to learn a singular subspace instead of a general linear space. While we are not sure if we understand the reviewer’s definition of a \"*general linear space,*\" the choice of learning singular subspaces for linear operator regression can be justified as follows. \n\n    For a general stochastic dynamical system, the corresponding Koopman operator may be non-normal, for which there is no spectral theorem. In such cases, there is no systematic way to recover the eigenfunctions of the Koopman operator that dominate the dynamics. Due to the inherent difficulty, prior works (i.e., VAMPnet and DPNet) instead formulate the problem as first learning the singular subspaces, which can be done systematically even for non-normal operators as long as they are compact, and then performing an operator regression step on top of them. For the special case of normal operators, EVD can be done by SVD, and thus SVD is fully justified. All in all, the choice of learning the singular subspaces is justified by the well-definedness of SVD for compact operators and the available computational techniques for computing SVD.\n\n2. **On Guaranteeing Orthogonality of Singular Vectors**  \n    Our method guarantees orthogonality in a two-step process:\n\n   1. **Implicit Orthogonality via Loss Minimization**: We implicitly enforce orthogonality by minimizing the objective function. Our LoRA loss function $-2\\\\mathrm{tr}(\\\\mathsf{T}\\[\\\\mathbf{f}, \\\\mathbf{g}\\]) \\+ \\\\mathrm{tr}(\\\\mathsf{M}\\[\\\\mathbf{f}\\] \\\\mathsf{M}\\[\\\\mathbf{g}\\])$ is derived from minimizing $\\\\|\\\\mathcal{K} \\- \\\\sum\\_i f\\_i \\\\otimes g\\_i\\\\|\\_{\\\\text{HS}}^2$. The global minimum of this loss is achieved when the learned functions $\\\\{f\\_i\\\\}$ and $\\\\{g\\_i\\\\}$ form bases for the true top-K singular subspaces. While the neural network outputs $\\\\mathbf{f}(x)$ and $\\\\mathbf{g}(x)$ may not be exactly orthogonal during training, the optimization process **drives them towards forming an orthogonal basis** for the target singular subspaces, as this is the configuration that minimizes the approximation error.  \n   2. **Explicit Orthonormalization for Inference**: For downstream tasks (like eigen-analysis or prediction), we do not assume the raw neural network outputs are perfectly orthogonal. Instead, as described in **Approach 1 (CCA \\+ LoRA, Sec 3.2)**, we perform a post-hoc canonical correlation analysis (CCA). This step explicitly computes the whitened bases $\\\\tilde{\\\\mathbf{f}}(x)$ and $\\\\tilde{\\\\mathbf{g}}(x')$ and the rotation matrices (from the SVD of $\\\\mathsf{T}\\[\\\\tilde{\\\\mathbf{f}},\\\\tilde{\\\\mathbf{g}}\\]$) to yield the final, strictly orthonormal singular functions $\\\\phi(x)$ and $\\\\psi(x')$.\n\n    In summary, our framework first learns a useful, near-orthogonal basis via stable loss minimization and then refines it into a strictly orthonormal basis for reliable downstream analysis.\n\n3\\. **Responses to Minor Comments**\n\nWe appreciate the reviewer's comments on typos and clarifications. We will carefully revise our manuscript, incorporating these points.\n\n* **Line 49**: We apologize for the confusion. We meant the posterior mean of $g(\\\\mathbf{x}\\_t)$ given a state $\\\\mathbf{x}\\_0$ at time $t=0$.   \n* **Line 102**: We meant the strong operator topology defined on a function space like Hilbert space $\\\\mathcal{H}$. That is, locally convex topology induced by seminorms $\\\\mathcal{T} \\\\mapsto \\\\|\\\\mathcal{T}(g)\\\\|\\_{\\\\mathcal{H}}$ ($g\\\\in\\\\mathcal{H}$). In this topology, $\\\\mathcal{T}\\_n \\\\rightarrow \\\\mathcal{T}$ iff $\\\\| \\\\mathcal{T}\\_n(g) \\- \\\\mathcal{T}(g) \\\\|\\_{\\\\mathcal{H}} \\\\rightarrow 0$ for all $g \\\\in \\\\mathcal{H}$.  \n* **Line 115**: Thank you for catching the typo. It should be \"even for\".\n\n---\n\nWe hope these clarifications address the reviewer's concerns, and if they do, we would be grateful if the reviewer would consider updating the score to reflect the contribution to the community."}]

Author rebuttal:

Meta-review: The paper addresses a key numerical instability in existing Koopman learning objectives and proposes LoRA to avoid ill-conditioned matrix inversions and for stable training. Reviewers raised concerns about the degree of novelty, clarity of nesting techniques, and reporting rigor, but after the extensive discussion phase were unanimously positive that this submission is a meaningful and practically important contribution. I agree with the eventually overall positive reception and recommend the authors follow through on their updates and promises from the discussion phase.
Decision: Accept (poster)