Estimating Koopman operators with sketching to provably learn large scale dynamical systems
Authors: Giacomo Meanti
Antoine Chatalic
Vladimir R Kostic
Pietro Novelli
Massimiliano Pontil
Lorenzo Rosasco
Venue: NeurIPS 2023 poster (2023)
Keywords: dynamical systems
kernel methods
koopman operator
sketching
molecular dynamics
efficient machine learning

Abstract: The theory of Koopman operators allows to deploy non-parametric machine learning algorithms to predict and analyze complex dynamical systems.
Estimators such as principal component regression (PCR) or reduced rank regression (RRR) in kernel spaces can be shown to provably learn Koopman operators from finite empirical observations of the system's time evolution. 
Scaling these approaches to very long trajectories is a challenge and requires introducing suitable approximations to make computations feasible. 
In this paper, we boost the efficiency of 
different kernel-based Koopman operator estimators using random projections (sketching).
We derive, implement and test the new ``sketched'' estimators with extensive experiments on synthetic and large-scale molecular dynamics datasets. 
Further, we establish non asymptotic error bounds giving a sharp characterization of the trade-offs between statistical learning rates and computational efficiency.
Our empirical and theoretical analysis shows that the proposed estimators provide a sound and efficient way to learn large scale dynamical systems.
In particular our experiments indicate that the proposed estimators retain the same accuracy of PCR or RRR, while being much faster.

Estimating Koopman operators with sketching to provably learn large scale dynamical systems
Giacomo Meanti Istituto Italiano di Tecnologia giacomo.. meanti@iit.it
Antoine Chatalic Università di Genova
antoine chatalic@dibris unige.it
Vladimir R. Kostic Istituto Italiano di Tecnologia University of Novi Sad vladimir. kosticliit.it
Pietro Novelli
Istituto Italiano di Tecnologia pietro. novellieiit.
Massimiliano Pontil Istituto Italiano di Tecnologia University College London massimil liano. wtw.iooooooniiiitt.... pontil@iit.it
Lorenzo Rosasco Università di Genova Istituto Italiano di Tecnologia
Massachusetts Institute of Technology Irosasco@mit... edu
Abstract
The theory of Koopman operators allows deploy non-parann non-parametric machine learn- ing algorithms to predict and analyze complex dynamical systems. Estimators such as principal component regression (PCR) reduced rank regression (RRR) in ker- nel spaces can be shown to provably learn Koopman operators from finite empirical observations of the system's 's time evolution. Scaling these approaches to very long trajectories is challenge and requires introducing suitable approximations to make computations feasible. In this paper, we boost the efficiency of different kernel- based Koopman operator estimators using random projections (sketching). We de- rive, implement and test the new "sketched" estimators with extensive experiments synthetic and large-scale molecular dynamics datasets. Further, we establish non asymptotic eтror bounds giving a sharp characterization of the tradc-offs between statistical learning rates and computational efficiency. Our empirical and theoretical analysis shows that the proposed estimators provide a sound and efficient way to learn large scale dynamical systems. In particular proposed estimators retain the same accuracy of PCR or RRR, while being much faster. Code is available at https: /tattps://githbbe htw.tthb.ccoooceee htt///ioodooo//// ooiooooooooyysstooooooooaaa
experiments indicate that the
Introduction
In the physical world, temporally varying phenomena everywhere, from biological processes in the cell to fluid dynamics to electrical fields. Correspondingly, they generate large amounts of data both through experiments and simulations. This data often analyzed in the framework of dynamical systems, where the state system æ observed at a certain time t, and the dynamics is described by function if which captures its evolution in time
Bt++ f(x:).
'Equal contribution
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

The function f must capture the whole dynamics, and as such it may be non-linear and even stochastic for instance when modeling stochastic differential equations, or simply noisy processes. Applications of this general formulation arise in fields ranging from robotics, atomistic simulations, epidemiology, and many more. Along with a recent increase in the availability of simulated data, data-driven techniques for learning the dynamics underlying physical systems have become commonplace. The typical approach of such techniques is to acquire a dataset of training pairs (æ,, yt. = xt+1) sampled in time, and use them to learn model for which minimizes a forecasting error. Since dynamical systems stem from real physical processes, forecasting is not the only goal and the ability to interpret the dynamics is paramount. One particularly important dimension for interpretation is the separation of dynamics into multiple temporal scales: fast fluctuations can e.g. be due to thermodynamical noise or electrical components in the system, while slow dynamics describe important conformational changes in molecules or mechanical effects.
Koopman operator theory [2 (8 provides an elegant framework in which the potentially non-linear dynamics of the system can be studied via the Koopman operator
(x))) = E [+ (x()))
which has the main advantage of being linear but is defined on typically infiniteddimensional set of observable functions. The expectation in (T) is taken with respect to the potential stochasticity of f. Thanks to its linearity, the operator K. can e.g. be applied twice to get two-steps-ahead forecasts, and one can compute its spectrum (beware however that K. is not self-adjoint, unless the dynamical process is (ime-reversible) Accurately upproximating the Koopman operator and its spectral properties is of high interest for the practical analysis of dynamical systems. However doing so efficiently for long temporal trajectories remains challenging. In this paper are interested in designing estimators which are both theoretically accurate and computationally efficient.
Related works Learning the spectral properties of the Koopman operator directly from data has been considered for at least three decades [39 40 resulting in large body of previous work. Among the different approaches proposed over time (see Mezić [4] for a recent review) it is most common to search for finite dimensional approximations to the operator, from which part of the spectrum and the Koopman modes [8] can be obtained. Dynamic mode decomposition (DMD) [56] 63
time-lagged independent component analysis (ICAA [42 49] and many subsequent extensions
for example can be seen minimizers of the forecasting error when & is restricted to be linear function of the states [52]. Extended DMD (cDMD) 666 25 and the variational approach for
conformation dynamics (VAC)
45
46
instead allow for (potentially learnable, as in recent deep dictionary of non-linear functions th. Kernel DMD 67 26
learning algorithms 32 37
69
62
and kernel tICA (57) are further generalizations which again approximate the Koopman operator but using an infinite dimensional space of features vi, encoded by the feature map of a reproducing kernel. While often slow from a computational point of view, kernel methods are highly expressive and can be analyzed theoretically, to prove convergence and derive learning rates of the resulting estimators [29] 30) Approximate kernel methods which are much faster to run have been recently used for Koopman operator learning by Baddoo et al. [6] where an iterative procedure is used to identify the best approximation the full kernel, but no formal learning rates are demonstrated, and by Ahmad et al. [3] who derive learning rates in Hilbert-Schmidt norm (while we consider operator norm) for the Nyström KRR estimator (one of the three considered in this paper).
Contributions In this paper we adopt the kernel learning approach. Starting from the problem of approximating the Koopman operator in a reproducing kernel Hilbert space, we derive three different estimators based on different inductive biases: kemel ridge regression (KRR) which comes from Tikhonov regularization, principal component regression (PCR) which is equivalent to DMD and its extensions, and reduced rank regression (RRR) which comes from a constraint on the maximum rank of the estimator [24]. We show how to overcome the computational scalability problems inherent in full kernel methods using an approximation based on random projections which is known as the Nyström method [58] 65 The approximate learning algorithms scale very easily to the largest datasets, with computational complexity which goes from O(n²) for the exact algorithm to ((2)) for the approximate one. We can further show that the Nyström KRR, PCR and RRR estimators have the same convergence rates as their exact, slow counterparts which are known to be optimal under our assumptions. We provide learning bounds in operator norm, which are known to translate to bounds for dynamic mode decomposition and are thus of paramount importance for applications. Finally, we thoroughly validate the approximate PCR and RRR estimators on synthetic dynamical systems, comparing efficiency and accuracy against their exact counterparts [29). as well as recently
2

---

proposed fast Koopman estimator streaming KAF 20) To showcase realistic scenario, we train on molecular dynamics simulation of the fast-folding Trp-cage protein [35].
Structure of the paper We introduce the setting in Section 2 and define our three estimators in Section In Section 4 provide bounds the excess risk of our estimators, and extensive experiments on synthetic as well as large-scale molecular dynamics datasets in Section
Background and related work
Notation We consider a measurable space (B where X' corresponds the state space, and denote 2 == L ( B, T) the L space of functions on X' w.r.t to a probability measure :, and L the space of measurable functions bounded almost everywhere. We denote HS(H) the space of Hilbert-Schmidt operators on a space H.
Setting The setting we will consider is that of Markovian, time-homogeneous stochastic process {Xi on X. By definition a Markov process, X2 only depends on X; and not on any previous states. Time-homogeneity ensures that the transition probability TP [ E B| = æ] for any measurable set B does not depend on f, and can be denoted with p(æ, B). This implies in particular that the distribution of (X,, X does not depend on t, and we denote it in the following. We further assume the existence of the invariant density which satisfies π (BB) [ TTx))p(x, B) dæ. This classical assumption allows one to study a large class of stochastic dynamical systems, but also deterministic systems on the attractor, see e.g. The Koopman operator K, (X) (x) is a bounded linear operator, defined by
(x)) / p(x, (((.,)))())) v)g(y) dy E
m], g (2)
(2)
We are in particular interested in the eigenpairs (A. E C x L that satisfy
Kxpp, = Appi.
(3)
Through this decomposition it is possible to interpret the system by separating fast and slow pro- cesses, or projecting the states onto fewer dimensions [15 19 7 In particular, the Koopman mode decomposition (KMD) allows to propagate the system state in time. Given an observable g X - R such that E spantyp; ji E N. the modes allow to reconstruct g(x) with a Koopman eigenfunction basis. The modes n Cd are the coefficients of this basis expansion:
(K> [[o(X.) æ = Σ Aspi(a)m
This decomposition describes the system's 's dynamics in terms of stationary component (the Koopman modes), temporal component (the eigenvalues X₂) and a spatial component (eigenfunctions sP:).
Kernel-based learning In this paper we approximate K₂ with kernel-based algorithms, using operators in reproducing kernel Hilbert spaces (RKHS) H. associated with kernel k : X x X + R and feature map á X + H. We wish to find an operator A H + which minimizes the risk
RHS A) = E, [AA, (x, y))]
where f(A, (x, y)) = ( - A
(5)
The adjoint of A, denoted by A", should thus be understood
estimator of the Koopman operator Kr in H. as will be clarified in (I5). In practice π and p are unknown, and one typically has access to dataset { } sampled from p. where each pair (x,, y; = f(x;)) may equivalently come from a single long trajectory or multiple shorter ones concatenated together. We thus use the empirical risk
1 Σ ((A, (æ;, y:)) ==1
RHS (A)
(6)
as a proxy for Since minimizing eq. (6) may require finding the solution to very badly conditioned linear system, different regularization methods (such as Tikhonov or truncated SVD) can be applied on top of the empirical risk.
3

---

Remark 2.1 (Connections to other learning problems); The problem of minimizing eqs. 5 and (
has strong connections to learning conditional mean embeddings (59 22 44 33) where the predic- tors and targets are embedded in different RKHSs, and to structured prediction fI2 13) which is an even more general framework. On the other hand, the most substantial difference from the usual kernel regression setting [9] is the embedding of both targets and predictors into a RKHS, instead of just targets.
We denote the input and cross covariance C = [o(z()) @ o(z)] and Cyxx Ey(o(y) @ d(z)), and their empirical counterparts as Ć > We also use the abbreviation Cx == C == j (x;) @ o(x;)) and Čyx Σ o(y)) @ of ).
XI. Minimizing the empirical risk (6) with Tikhonov
regularization [ yields the following KRR estimator
Ax arg min (A) + A|AAIis = CCx(C ACHS(H)
Eq. 7 can be computed by transforming its expression with the kernel trick [23] to arrive at a form where one must invert the kernel matrix a X n matrix whose i,j-th entry is k( a.y)-
This operation requires ((z²) time and (n²)) memory, severely limiting the scalability of KRR to
n S 100000 points. Improving the scalability of kernel methods is a well-researched topic, with the most important solutions being random features [50] 51
68 21 and random projections [58 65 21].
In this paper we use the latter approach, whereby the kernel matrix is assumed to be approximately low-rank and is sketched to a lower dimensionality, In particular we will use the Nyström method to approximate the kernel matrix projecting it onto a small set of inducing points, chosen among the training set. The sketched estimators the training trajectories become longer. For example, the state of the art complexity for solving (non vector valued) approximate kernel ridge regression is (n//n) time instead of ( [38 10).
much more efficient than the exact ones, increasingly as
Furthermore, when enough inducing points are used (typically on the order of vT), the learning rates of the exact and approximate estimators are the same, and optimal [5] [33. Hence is possible- and
in this paper we show it for learning the Koopman operator to obtain large efficiency gains, without losing anything in terms of theoretical guarantees of convergence.
3 Nyström estimators for Koopman operator regression
In this section, we introduce three efficient approximations of the KRR, PCR and RRR estimators of the Koopman operator. Our estimators rely on the Nyström approximation, i.e. on random projections onto low-dimensional subspaces of H spanned by the feature-embeddings of subsets of the data.
We thus consider two sets of m < n inducing points { C (0x)]= and {ū, ) C {yr } sampled respectively from the input and output data. The choice of these inducing points (also sometimes called Nyström centers) is important to obtain a good approximation. Common choices include uniform sampling, leverage score sampling [17 55), and iterative procedures such as the one used in 6 to identify the most relevant centers. In this paper we focus on uniform sampling for simplicity, but we stress that our theoretical results in Section
I can easily be extended leverage scores sampling by means of [53, Lemma 7]. To formalize the Nyström estimators, we define operators &x by I - H as xxw = == ())((ā))))) and Ty Σ,= ", and denote
Px and Py the orthogonal projections onto span Φx and span by respectively.
In the following paragraphs we apply the projection operators to three estimators corresponding to different choices of regularization. For each of them a specific proposition (proven in Section C) states efficient way of computing it based on the kernel trick. For this purpose we introduce the kernel matrices K x,x; Ky E Rmxxn between training set and inducing points with entries (K (xxx)); = k( I.), (Ky,r), = k(ū W:), and the kernel matrices of the inducing points K & K E R mm with entries (K xxxx)ax k )x)) and (K xxxxx)a ) k( (x).
Kernel Ridge Regression (KRR) The cost of computing A defined in Eq. 7 is O(n²) 29 which
prohibitive for datasets containing long trajectories. However, applying the projection operators to each side of the empirical covariance operators, we obtain an estimator which additionally depends on the m inducing points:
AKRR = PyCxXPxPPPCPP + +XI- : - H.
(8)
4

---

H. is infinite dimensional, Eq. (8) cannot be computed directly. Proposition 3.1 provides a computable version of the estimator.
proven in Section
Proposition 3.1 (Nyström KRR): The Nyström KRR extimator 8 can be expressed as
A ARRR = bx Ky y.y Ky y x (Kxxx X X x
The computational bottlenecks are the inversion of an X m matrix and a large matrix multiplication, which overall need O(2m² + 2m²n) operations. In particular, in Section 4 we will show that /Ti is sufficient guarantee optimal rates even with minimal assumptions, leading to a final cost of ( Note that similar estimator was derived (3).
Please note that the O(n²) cost is for a straightforward implementation, and can indeed be reduced via iterative linear solvers (possibly preconditioned, to further reduce the practical running time), and randomized linear algebra techniques. In particular, we could leverage results from Rudi et al. (54) to reduce the computational cost to (nv/n)
Principal Component Regression (PCR) Typical settings in which Koopman operator theory is used focus on the decomposition of dynamical system into small set of components, obtained from the eigendecomposition of the operator itself. For this reason, a good prior on the Koopman estimator is for to be low rank. The kernel PCR estimator ÁPCR = Ey x [C] formalizes this concept (29 677, where here [-, denotes the truncation to the first components of the spectrum. Again this is expensive to compute when n. is large, but the estimator can be sketched as follows:
APCR PyCyx|PxxPxst
(10)
The next proposition provides an efficiently implementable version of this estimator.
Proposition 3.2 (Nyström PCR): The sketched PCR estimator 10 satisfies
PCR कy y
K x X,X
(1I)
requiring O(2m² + 2m²n) operations, i.e. optimal rates again be obtained cost of most O(n²) operations.
Note that with analog forecasting (KAF) [.] The sketched estimator of Proposition 3.2 was also recently derived in [6), albeit without providing theoretical guarantees.
APCR is equivalent to the kernel DMD estimator (67. also known as kernel
Reduced Rank Regression (RRR) Another way to promote low-rank estimators is to add an explicit rank constraint when minimizing the empirical risk. Combining such constraint with Tikhonov regularization corresponds to the reduced rank regression [24 29 estimator:
ARRB arg min (cHS)r
Rus( () + A|AA|
(12)
Minimizing Eq. 12 requires solving a n X n generalized eigenvalue problem. The following proposition introduces the sketched version of this estimator, along with a procedure compute which instead requires the solution of a m x m eigenvalue problem. For m x /n, which is enough to guarantee optimal learning rates with minimal assumptions (see Section 4 this represents reduction from ((n²) to (n)/n) time.
Proposition 3.3 (Nyström RRR): The Nyström RRR estimator can be written as
ARRR = [Py (xPPxPx + 1 12xxxx + XI)
(13)
To compute it, solve the m x eigenvalue problem
Kxx + +nAKx x) Kx K Ky K. v. S.rskyy Ky, K xx,x = o?
for the first r eigenvectors W, = [w),... wrl. appropriately normalized. Then denoting D, := y Ki y. Kg,y y K x W, and F, == (K x K x ~>>K xxx)) KKKK x.x xyx D, it holds
ARRR Ty D
(14)
5

---

4 Learning bounds in operator norm for the sketched estimators
In this section, we state the main theoretical results showing that optimal rates for operator learning with KRR. PCR and RRR can be reached with Nyström estimators.
Assumptions We first make two assumptions on the space H used for the approximation, via its reproducing kernel k:
Assumption 4.1 (Bounded kernel): There exists K 0o such that ess
ess sup g s (l 1 K.
Assumption 4.1 ensures that H is compactly embedded in I 61 Lemma 2.3], and we denote T the embedding operator which maps any function in H to its equivalence class r-llmost everywhere in I.
0 H
Assumption 4.2 (Universal kernel); The kernel is universal, i.e. cl(ran()))) L
We refer the reader to 160, Definition 4.52] for a definition of universal kernel. The third assumption on the RKHS is related to the embedding property from Fischer and Steinwart (18). connected the embedding of interpolation spaces. For detailed discussion see Section A..
Assumption 4.3 (Embedding property): There exists t E [0, and C
0 such that
suppa~ ||Cx 1/2 ((xx)) c,A
Next, make an assumption on the decay of the spectrum of the covariance operator that is of paramount importance for derivation of optimal learning bounds. In the following, (x(A) and (v(A) always denote the eigenvalues and singular values of operator A (in decreasing order).
Assumption 4.4 (Spectral decay): There exists 3 E [0, and > 0 such that >c S 1/8
This assumption common in the literature, and we will see that the optimal learning rates depend on B. It implies the bound (Atr(A) == CC( S X on the effective dimension, which is key quantity in the analysis (both statements are actually equivalent, see Section E.2). Note that (A) Ezsr /// <(x) < ess sup ||/= ((x)) and thus it necessarily holds B S For a Gaussian kernel, both 3 and T can be chosen arbitrarily close to zero.
Finally, we make an assumption about the regularity of the problem itself. A common assumption occurring in the literature is that E[](X,) | X, define the Koopman operator directly on the space H. i.e. the learning problem is well-specified. However, this assumption is often too strong. Following 30 D.1] we make a different assumption on the cross-covariance remarking that, irrespectively of the choice of RKHS, it holds true whenever the Koopman operator is self-adjoint (i.e. the dynamics is reversill)....
E H for every / E H. meaning that one can
Assumption 4.5 (Regularity of K.;): There exists a > such that Cxy $ a²c².
Rates The risk can be decomposed as Rus = Ens (A) + Rys,o where R1S,0 is a constant and (As = [K x - x A' llis corresponds to the excess risk (more details in Section B Optimal learning bounds for the KRR estimator the context of CME (i.e. in Hilbert-Schmidt norm) have been developed in [33] under Assumptions 4.1 to 4.4 in well-specified and misspecified settings. On the other hand, in the context of dynamical systems, Kostic et al. [29 Theorem I] report the importance of reduced rank estimators that have a small excess risk in operator norm
((A) = ||, - 1'
(15)
The rationale behind considering the operator norm is that it allows to control the error of the eigenvalues approximation and thus of the KMD B). 4 as discussed below. Optimal learning bounds in operator norm for KRR, PCR and RRR are established in (30]. In this work show that the same optimal rates remain valid for the Nyström KRR, PCR and RRR estimators. According to [29] and (30] these operator norm bounds lead reliable approximation of the Koompan mode decomposition of Eq.
We now provide main result.
Theorem 4.6 (Operator norm error for KRR, i.i.d. data): Let assumptions 4.1 to 4.5 hold. Let
(a,; y1 sisn be i.i.d. samples, and let Py Px be the projection induced by m Nyström landmarks drawn uniformly from (wi)isisn without replacement. Let X cxn -1/(+8) where ex a constant

---

given in the proof, and assume 2 (cx/
Then holds with probability least 1 8
E(A ) - n 2ा
1/2
provided
m 2 max(1, (11(+))))) log(1/))
The proof is provided in Section E.2 but essentially relies decomposition involving the terms II EY1I (C- I], as well as bounding the quantity |P where Px denotes the projection on the orthogonal of ran( Px)- bounded using variants of the Bernstein inequality. Note that results can easily be extended to leverage score sampling of the landmarks by bounding term || C1/21| by means of [53] Lemma 7]; the same rate could then be obtained using a smaller number m of Nyström points.
ran(Px). All these terms are
The rate n -=1/2((1-8) is known to be optimal (up the log factor) in this setting by assuming an additional lower bound the decay of the covariance's ' eigenvalues of the kind (C 2 5-1/8,
30 Theorem 7 in D.4]. One can see that without particular assumptions (B = = 1), we only need the number of inducing points be of the order of (( in order get optimal rates. For r fixed, this number increases when / decreases (faster decay of the covariance 's spectrum), however note that the optimal rate depends on and also improves in this case. The dependence in is particularly interesting, as for instance with a Gaussian kernel it is known that 7 can be chosen arbitrarily closed to zero 33 18 order of (log.)).
In that case, the number of inducing points can be taken on the
Note that bound for the Nyström KRR estimator has been derived in Hilbert-Schhie norm by Ahmad et al. 3] Using the operator norm however allows to derive bounds on the eigenvalues (see discussion below), which is of paramount importance for practical applications. Moreover, now provide a bound on the error of PCR and RRR estimators, which are not covered in (3.
Lemma 4.7 (Operator norm error for PCR and RRR, i.idd data): Under the assumptions of Theorem 4.6 taking X = cxn 1/((+9) with ex as in Theorem 4.6 n 2 (cx/K2)1+8, and provided
m 2 max(1, (1++)))) log(m/8),
it holds with probubility least 1 8
El {(AmR)) 1//2 S CRRR n गग-क,, for s.t. (+x1((xX;;) min((, ((4T))
and E(APOR) S CPCR n भमंग, for r > n मााक
where CRRR = (0, (4x and crcn (o,(4x) 0x+1(4x))- are the problem dependant constants.
Note that when rank of K-z is then there is no restriction for the RRR estimator, while for PCR the choice of depends the spectral decay property of the kernel. In general, if > ,,ीक,,, then 0,+1 ( KC= < o,+1(4x S n (1(((())))) which implies that RRR estimator achieve the same rate of PCR but with smaller rank. Again the rate is sharp (up to the log factor) in this setting (30].
Koopman mode decomposition According to 29 Theorem 1], working in operator norm allows us to bound the error of our estimators for dynamic mode decomposition, as well as to quantify how close the eigenpairs (7,, Φ of an estimator Á are to being eigenpairs of the Koopman operator. Namely, recalling that for function pú. the corresponding candidate for Koopman eigenfunction in ( //)l1/l(( &all < 1/ || &all- While
L space is 4 Pú, one has ||C,, ( eigenvalue and eigenfunction learning rates were studied, under additional assumptions, in [30], where the operator norm error rates were determinant, here, in Section we empirically show that the proposed estimators accurately learn the Koopman spectrum. We refer the reader Section for the details on computation of eigenvalues, eigenfunctions and KMD of an estimator in practice.
Dealing with non-i.i.d. data The previous results hold for i.i.d. data, which is not very realistic assumption when learning from sampled trajectories. Our results can however easily be extended to 9-mixing processes by considering random variables Z: = ード X,+i (thus representing portions of the trajectory) sufficiently separated in time to be nearly independent. We now consider a trajectory
7

---

FullRRR
NyRRRR
FUIPCR NyPCRR
0.75
Elapsed time (s)
1.2
nRMSE
0.50
nRMSE
1.0
Time (s)
10P
0.5
NysRRR
sKAF
0.25
10-1
NysPCR
0.5
10-2
0.0
sKAF NysRRR NysPCR
Number of samples
Number tamples
Forecast time (s)
Figure I: Full and Nyström estimators trained on L63 with increasing Error (left) and running time (right) are plotted to show efficiency gains without accuracy loss with the Nyström approxi- mation. RBF = 3.5) kernel, r = 25 principal components and 250 inducing points.
Figure 2: Nyström and sKAF estimators trained on L63 for increasing forecast horizons; the error (left) and overall running times (right) are shown. We used a RBF kernel with 3.5, 50, m = 250 (for Nyström methods) and v log n random features (for sKAF).
J.8 (re-stated from (29] which allows to translate concentration results on the Z; to concentration the X; by means of the 3-mixing coefficients defined as (x)))= := sup (B( - (7 x (B) where Ps denotes the joint probability of Xx+x). Using this result the concentration results provided in appendix can thus be generalized to the 3-mixing setting, and apart from logarithmic dependencies we essentially obtain similar results to the i.i.d. setting except that the sample size n replaced by p R n/(2)..
$1.... Dn+1 with T1 ~ π and Ø(+1 ~ -(x for E [I. and use Lemma
5 Experimental validation
In this section show how the estimators proposed in section 3 perform in various scenarios, ranging from synthetic low dimensional ODEs to large-scale molecular dynamics simulations. The code for reproducing all experiments available online. Our initial aim is to demonstrate the speed of NysPCR and NysRRR, compared to the recently proposed alternative Streaming KAF (sKAF) [20]. Then we show that their favorable scaling properties make it possible train on large molccular dynamics datasets without any subsampling. In particular we run metastability analysis of the alanine dipeptide and the Trp-cage protein, showcasing the accuracy of our models' eigenvalue and eigenfunction estimates, as well as their efficiency on massive datasets (> 500 points)
Efficiency Benchmarks Lorenz '63 The chaotic Lorenz '63 system 36 consists of ODEs with no measurement noise. With this toy dynamical system we can easily compare the Nyström estimators two alternatives: the corresponding exact estimators and 2. the sKAF algorithm which also uses randomized linear algebra to improve the efficiency of PCR. In this setting we sample long trajectories from the system, keeping the first points for training (the number of training points varies for the first experiment, and is fixed to 000 for the second, see fig. 2). and the subsequent ones for testing. In fig. we compare the run-time and accuracy with of NysPCR and NysRRR versus their full counterparts. To demonstrate the different scaling regimes we fix the number of inducing points 250 and increase the number of data points n. The accuracy of the two solvers (as measured with the normalized RMSE metric (nRMSE) (20 on the first variable) is identical for PCR and close for RRR, but the running time of the approximate solvers increases much slower with n than that of the exact solvers. Each experiment is repeated 20 times to display error bars over the choice of Nyström centers. In the second experiment, shown in fig. 2 we reproduce the setting of (20] by training at increasingly long forecast horizons. Plotting the nRMSE we verify that sKAF and NysPCR converge to very similar accuracy values, although NysPCR is approximately 10 times faster. NysRRR instead offers slightly better accuracy, the expense of higher running time compared to NysPCR. Error bars are the standard deviation of nRMSE over successive test sets with 10 000 points each. In section K we show additional experiments including comparison to the Nyström KRR estimator.
Molecular dynamics datasets An important application of Koopman operator theory is in the analysis of molecular dynamics (MD) datasets, where the evolution of molecule's atomic positions as they evolve over time modelled. Interesting systems are very high dimensional, with hundreds or thousands of atoms. Furthermore, trajectories generated at very short time intervals (< 1 ns) but

---

A 0.924
Aa 0.236
b
metastable state
77
Figure 3: Dynamics of the alanine dipeptide liggtimee 100), Nyström RRR model. On the left the Ramachandran plot of all protein conformations is colored with the value of the first two non-constant eigenfunctions evaluated on the 45-d space. On the right color indicates the discrete state obtained by clustering the dynamics projected onto the first eigenfunctions (using PCCA+ with 3 states). The obtained clusters match the main areas the Ramachandran plot.
interesting events (e.g. protein folding/unfolding) occur timescales on the order of at least 10 so that huge datasets are needed to have a few samples of the rare events. The top eigenfunctions of the Koopman operator learned on such trajectories can be used to project the high-dimensional state space onto low- www-dimensional coordinates which capture the long term, slow dynamics.
We take three 250 ns long simulations sampled ps of the alanine dipeptide [64), which is often taken as a model system for molecular dynamics [47 46]. We use the pairwise distances between 45-dimensional space. We train NysRRR model with 0000 centers on top of the full dataset (449 940 points are used for training, the rest for validation and testing) with lag time 100 ps. and recover 2-dimensional representation which correlates well with the Φ, vp backbone dihedral angles of the molecule, known to capture all relevant long-term dynamics. Figure 3 shows the top two eigenfunctions overlaid onto Ø, 1h the first separates the slowest transition between low and high 6; the second separates low and high v. The implied time-scales from the first two non- coi-trivial eigenvalues are 1262 ps and 69 ps. which are close the values reported by Nüske et al. [47] (1400 ps and 70 ps) who used a more complex post-processing procedure to identify time-scales. We then train PCCA+ [66] model on the first three eigenfunctions to obtain three states, as shown in fig. Pp. PCCA+ acts on top of a fine clustering (in our case obtained with k-means, k 50), to find the set of maximally stable states by analyzing transitions between the fine clusters. The coarse clusters clearly correspond to the two transitions described above.
heavy atoms as features, yielding 45
Finally take 208 μs long simulation of the fast-folding Trp-cage protein (35, sampled every 0.2 ns. Again, the states are the pairwise distances between non-hydrogen atoms belonging to the protein, in 296 dimensions. A NysRRR model is trained 626 370 points, using 5000 centers in approximately 10 minutes. Note that without sketching this would be completely intractable problem. Using lag-time 10 we observe spectral gap between the third and fourth eigenvalues. hence we train PCCA+ model on the first 3 cigenfunctions to obtain the states shown in fig. 4 The first non-trine con-trivial Koopman eigenvector effectively distinguishes between the folded (state and unfolded states as is evident from the first row of fig. 4 The second one instead can be used to identify a partially folded state of the protein (state 0), as can be seen from the insets in fig. 4
Conclusions
We introduced three efficient kernel-based estimators of the Koopman operator relying on random projections, and provided bound on their excess risk in operator norm which of paramount importance to control the accuracy of Koopman mode decomposition. Random projections allow to process efficiently even the longest trajectories, and these gains come free as our estimators still enjoy optimal theoretical learning rates. We leave for future work the refinement our analysis under e.g. an additional source condition assumption in the misspecified setting. Another future research direction shall be to devise ways to further reduce the computational complexity of the estimators.

---

State 0
State 1
State 2
Eigenfunction
Eigenfunction 2
Figure 4: First cigenfunctions for Trp-cage dynamics, colored according to the membership probabil- ity for each state PCCA+ model. The bottom insets show a few overlaid structures from each state. The first eigenfunction exhibits a strong linear separation between state (folded) and the other states, which highlighted by the black lines. The second separates between state (partially folded) ant the rest. NysRRR model trained with 5000, = 10, RBF(a 0.02) kernel, À 10-10
7 Acknowledgements
This paper is part of project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 819789). L. R. acknowledges the financial support of the European Research Council (grant SLING 819789), the AFOSR projects FA9550-18-1 18-1-7009 FA9550-17-1-0390 and BAA-AFRL-AFOSHE 2016-0007 (European Office of Aerospace Research and Development), the EU H2020-MSCA-RISE project NoMADS DIV-777826, and the Center for Brains, Minds and Machines (CBMM), funded by NSF STC award CCF-1231216 M. P., V. K. and P. N. acknowledge financial support from PNRR MUR project PE0000013-FAIR and the European Union (Projects 951847 and 101070617).
References
[1] Amirhesam Abedsoltan, Mikhail Belkin, and Parthe Pandit. Toward large kernel models, 2023.
arXiv: 230200000000260 02605
LG]
[2] Thomas N. E. Greville Adi Ben-Israel. Generalized Inverses: Theory and Applications. CMS Books in Mathematics. Springer, 2nd edition, 2003.
[3] Tamim EI Ahmad, Luc Brogat-Motte, Pierre Laforgue, and d' Alché-Bue Buc Florence. Sketch In, Sketch Out: Accelerating both Learning and Inference for Structured Prediction with Kernels, 2023. arxiv: 2312.10128 10128
[4] Romeo Alexander and Dimitrios Giannakis. Operator-theoretic framework for forecasting nonlinear time series with kernel analog techniques. Physica D: Nonlinear Phenomena, 409, 2020. doi: https://doi.or//I 10. 101600 pwysd.2000.11...... 1020.132520.
[5] Francis Bach. Sharp analysis of low-rand kernel matrix approximations. Journal of Machine Learning Research, 30, 2013.
[6] J.ter J. Baddoo, Benjamin Herrmann Beverley J. McKeon, and Steven I. Brunton. Kernel learn- ing for robust dynamic mode decomposition: linear and nonlinear disambiguation optimization. Proceedings of the Royal Society A, 2022. doi: wttp://orrgg...000//aa....2..2
10

---

[7] Steven L. Brunton, Marko Budišić, Eurika Kaiser, and J. Nathan Kutz. Modern koopman theory for dynamical systems. SIAM Review, 64(2):29--340, 2022. doi: 10.1137/21M1401243.
[8] Marko Budišic, Ryan Mohr, and Igor Meziè. Applied koopmanism. Chaos, 22(4), 2012. doi: 2..1066/1.4721991
[9] Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squaree algorithm. Foundations of Computational Mathematicx, 7(3):31 368, 2007. doi: 10.1007/ s10208-006-0196-8
[10] Luigi Carratino, Stefano Vigogna. Daniele Calandriello, and Lorenzo Rosasco. ParK: Sound and efficient kernel ridge regression by featurespace partitions. In Advances in Neural Information Processing Systems 34, pages 6430-6441, 2021.
[11] Antoine Chatalic, Nicolas Schreuder, Lorenzo Rosasco, and Alessandro Rudi. Nyström Kernll Mean Embeddings. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 3006-3024. PMLR, 2022.
[12] Carlo Ciliberto, Lorenzo Rosasco, and Alessandro Rudi. A Consistent Regularization Approach for Structured Prediction. In Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016.
[13] Carlo Ciliberto, Lorenzo Rosasco, and Alessandro Rudi. A general framework for consistent structured prediction with implicit loss embeddings. Journal of Machine Learning Research, 21 (1), 2022.
[14] G. Da Prato and J. Zabezyk. Ergodicity for Infinite Dimensional Systems. London Math- ematical Society Lecture Note Series. Cambridge University Press, 1996. doi: 10.1017/ TB0978051166288
[15] Michael Dellnitz and Oliver Junge. On the approximation of complicated dynamical behavior. SIAM Journal on Numerical Analysis, 36(2):491-515) 1999. doi: 10.11137/SS3333114299333333
[16] Peter Deuflhard and Marcus Weber. Robust perron cluster analysis in conformation dynamics. Linear Algebra and its Applications, 398: 1114144 2005, doi: https://doi.org/1..11100. la.. 2004.10.026.
[17] Petros Drineas, Malik Magdon-Ismail, Michael W. Mahoney, and David P. Woodruff. Fast approximation of matrix coherence and statistical leverage. Journal of Machine Learning Research, 13(1), 2012.
[18] Simon Fischer and Ingo Steinwart. Sobolev norm learning rates for regularized least-squares algorithms. Journal of Machine Learning Research, 21(1):8464 8501, 2020.
[19] Gary Froyland, Georg A. Goutwal, and Andy Hammerlindl. A computational method to extract macroscopic variables and their dynamics in multiscale systems. SIAM Journal on Applied Dynamical Systems, 13(4):1816-14) 2014. doi: 10. 11711337799
[20] Dimitrios Giannakis, Amelia Henriksen, Joel A. Tropp, and Rachel Ward. Learning to forecast dynamical systems from streaming data. SIAM Journal on Applied Dynamical Systems, 22(2): 527-558, 2023, doi: 10. 1
[21] Alex Gittens and Michael W. Mahoney. Revisiting the nyström method for improved large-scale machine learning. Journal of Machine Learning Research, 17:3977 2016.
[22] Steffen Grünewälder, Guy Lever, Luca Baldassarre, Sam Patterson, Arthur Gretton, and Mas- similano Pontil. Conditional mean embeddings as regressors. In Proceedings of the 29th International Conference on Machine Learning, pages 1803- 1810, 2012.
[23] Thomas Hofmann, Bernhard Schölkopf, and Alexander J. Smola. Kernel methods in machine learning. The Annals of Statistics, 36(3):1171 1220, 2008. doi: 10.1214/ 009053607000000677

---

[24] Alan Julian Izenman. Reduced-rank regression for the multivariate linear model. Journal of Multivariate Analysis, 5(2):248- 264,, 1975. doi: https://doi.org/ wttp://.........0000000022090919 90042-1.
[25] Stefan Klus, Péter Koltai, and Christof Schüitte. On the numerical approximation of the Perron- Frobenius and Koopman operator. Journal of Computational Dynamics, 3(1):51-79, 2016. ISSN 2158-2491. doi: h..3934/jd.2001600
[26] Stefan Klus, Ingmar Schuster, and Krikamol Muandet. Eigendecomppoitiion of transfer opera- tors in reproducing kernel hilbert spaces. Journal of Nonlinear Science, 30(1):283-315, 2020. doi: 1010077/00332--19--9544444
[27] B. O. Koopman. Hamiltonian systems and transformation in hilbert space. Proceedings of the National Academy of Sciences, 17(5):315-318, 1931. doi: 10.1073/pnas. 17.53// 17.5.315.
[28] B. O. Koopman and J. v. Neumann. Dynamical systems of continuous spectra. Proceedings of the National Academy of Sciences, 18(3):255 263, 1932. doi: 10.1073/pnas. 18.3.255.
[29] Vladimir Kostic, Pietro Novelli, Andreas Maurer, Carlo Ciliberto, Lorenzo Rosasco, and Massi- miliano Pontil. Learning dynamical systems via Koopman operator regression in Reproducing Kernel Hilbert Spaces. In Advances in Neural Information Procexsing Systems 35, 2022.
[30] Vladimir Kostic, Karim Lounici, Pietro Novelli, and Massimiliano Pontil. Koopman operator learning: Sharp spectral rates and spurious eigenvalues, 2023. arxiv: 2302. 02004 [cs.LG]
[31] J. Nathan Kutz, Steven L. Brunton, Binghi W. Brunton, and Joshua Proctor. Dynamic Mode Decomposition:
Modeling of Complex Systems. SIAM, 2016.
[32] Qianxiao Li, Felix Dietrich, Erik M. Bollt, and Ioannis G. Kevrekidis. Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the koopman operator. Chaos: An Interdisciplinary Journal of Nonlinear Science, 27(10), 2017. doi: 1..11666149999555
[33] Zhu Li, Dimitri Meunier, Mattes Mollenhauer, and Arthur Gretton. Optimal Rates for Regular- ized Conditional Mean Embedding Learning. In Advances in Neural Information Processing Systems 35, 2022.
[34] Junhong Lin and Volkan Cevher. Optimal convergence for distributed learning with stochastic gradient methods and spectral algorithms. Journal of Machine Learning Research, 21(147): 1-63, 2020.
[35] Kresten Lindorff-Larsen, Stefano Piana, Ron O. Dror, and David E. Shaw. How fast-folding proteins fold. Science, 334(6055):517-520, 2011. doi: 10. 126/science. 1 1208351.
[36] Edward N. Lorenz. Deterministic nonperiodic flow. Journal of Atmospheric Sciences, 20(2): 130 - 141, 1963. doi:
[37] Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton. Deep learning for universal linear embeddings of nonlinear dynamics. Nature Communications, 9(1), 2018 doi: 10.1038/ s41467-018-07210-00
[38] Giacomo Meanti, Luigi Carratino, Lorenzo Rosasco, and Alessandro Rudi. Kernel methods through the roof: handling billions of points efficiently. In Advances in Neural Information Processing Systems 32, 2020.
[39] Igor Meziè. On the geometrical and statistical properties of dynamical systems: Theory and applications. PhD thesis, ProQuest LLC, California Institute of Technology, 1994.
[40] Igor Mezič. Spectral properties of dynamical systems, model reduction and decompositions. Nonlinear Dynamics, 41:309-325, 2005. doi: 10. 10700002107--00--282244
[41] Igor Mezić. Koopman operator, geometry, and learning of dynamical systems. Notices of the American Mathematical Society, 68(7):1087-1105 2021.
12

---

[42] L. Molgedey and H. G. Schuster. Separation of mixture of independent signals using time delayed correlations. Phys. Rev. Lett., 12:3634-3637, 1994. doi: 10.1100/hhysReeLLett2....23 2......6.4.
[43] Mattes Mollenhauer, Ingmar Schuster, Stefan Klus, and Christof Schütte. Singular value decomposition of operators on reproducing kernel hilbert spaces. In Advances in Dynamics, Optimization and Computation, pages 109-131, 2020.
[44] Krikamol Muandet, Kenji Fukumizu, Bharath Sriperumbudur, and Bernhard Schölkopf. Kernel mean embedding of distributions: review and beyond. Foundations and Trends in Machine Learning, 10(1-2):1-141, 2017. doi: 10. 10.1561/2200000066)
[45] Frank Noé and Feliks Nüske. A variational approach to modeling slow processes in stochastic dynamical systems. Multiscale Modeling Simulation, 11(2):635-655) 2013. doi: 10.1 137/ 118588116.
[46] Feliks Nüske, Bettina Keller, Guillermo Herzz-Hernnnde, Antonia S. J. S. Mey, and Frank Noé. Variational approach to molecular kinetics. Journal of Chemical Theory and Computation, 10(4):1739- 1752, 2014. doi: /C:1091/:
[47] Feliks Nüske, Hao Wu, Jan- Hendrik Prinz, Christoph Wehmeyer, Cecilia Clementi, and Frank Noé. Markov state models from short non-equilibrium simulations - analysis and correction of estimation bias. Journal of Chemical Physics, 146(9), 2017. doi: 00.1063/1.497655
[48] I. F. Pinelis and A. I. Sakhanenko. Remarks on inequalities for large deviation probabilities. Theory of Probability & Its Applications, 30(1):143-144, 1986. doi: 10.1137/1 130013.
[49] Guillermo Identification of slow molecular order parameters for markov model construction. Journal of Chemical Physics, 139(1), 07 2013. doi: 0063//1.8811 1489.
Fabian Paul, Toni Giorgino, Gianni De Fabritiis, and Frank Noé.
[50] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances Neural Information Processing Systems 20, 2008
[51] Ali Rahimi and Benjamin Recht. Weighted sums of random kitchen sinks: Replacing minimiza- tion with randomization in learning. In Advances in Neural Information Processing Systems 21, 2009.
[52] Clarence W. Rowley, Igor Meziè, Shervin Bagheri, Philipp Schlatter, and Dan S. Henningson. Spectral analysis of nonlinear flows. Journal of Fluid Mechanics, pages 115-127, 2009, doi: 10.1017/500221 12009992059.
[53] Alessandro Rudi, Raffaello Camoriano, and Lorenzo Rosasco. Less is more: Nyström com- putational regularization. In Advances in Neural Information Processing Systems 15, pages 1657-1665, 2015.
[54] Alessandro Rudi, Luigi Carratino, and Lorenzo Rosasco. FALKON: An optimal large scale kernel method. In Advances in Neural Information Processing Systems 30, 2017.
[55] Alessandro Rudi, Daniele Calandriello, Luigi Carratino, and Lorenzo Rosasco. On fast leverage score sampling and optimal learning. In Advances in Neural Information Processing Systems 31, 2018
[56] Peter Schmid. Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics, 656:5-28, 2010. doi: 10. 10.1077//002111 12010001111111111
[57] Christian R. Schwantes and Vijay S. Pande. Modeling molecular kinetics with tICA and the kernel trick. Journal of Chemical Theory and Computation, 11(2):600-608, 2015. doi: 10. 10.05377730
[58] Alex Smola and Bernhard Schölkopf. Sparse greedy matrix approximation for machine learning. In ICML 17, 2000.
13

---

[59] Le Song, Jonathan Huang, Alex Smola, and Kenji Fukumizu. Hilbert space embeddings of conditional distributions with upplications to dynamical systems. In Proceedings of the 26th Anmual International Conference on Machine Learning, page 961-968, 2009. doi: 10. 145/
1553374.
[60] Ingo Steinwart and Andreas Christmann, Support Vector Machines. Springer Science & Business Media, 2008. URL https://linkk swringerrsprnngrrrreeeee
[61] Ingo Steinwart and Clint Scovel. Mercer's theorem on general domains: On the interaction between measures, kernels, and RKHSs. Constructive Approximation, 35(3):363-417, 2012.
[62] Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi. Learning koopman invariant sub- spaces for dynamic mode decomposition. In Advances Neural Information Processing Systems, page 1130-1140, 2017.
[63] Jonathan H. Tu, Clarence W. Rowley Dirk M. Luchtenburg, Steven L.. Brunton, and J. Nathan Kutz. On dynamic mode decomposition: Theory and applications. Journal of Computational Dynamics, 1(2):391-421, 2014. doi: 0..393/jjcd.21441..3
[64] Christoph Wehmeyer and Frank Noé. Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics. Journal of Chemical Physics, 148(24), 2018, doi: 10.1063/1.501 1399.
[65] Christopher K. I. Williams and Matthias Seeger. Using the Nyström method speed up kernel machines. In Advances in Neural Information Processing Systems 13, 2001.
[66] Mattew O. Williams, Ioannis G. Kevrekidis, and Clarence W. Rowley, A data driven approxima- tion of the Koopman operator: Extending dynamic mode decomposition. Journal of Nonlinear Science, 25(6): 1307 - 1346, 2015. doi: 10. 1007//00032--15--9233
[67] Matthew O. Williams, Clarence W. Rowley, and loannis G. Kevrekidis. A kernel-based method for data-driven Koopman spectral analysis. Journal of Computational Dynamics, 2(2):247-265, 2015. ISSN 1158-24911 doi: 10.3934/jed.20155000
[68] Tianbao Yang, Yu-Feng Li, Mehrdad Mahdavi, Rong Jin, and Zhi-Hua Zhou. Nyström method vs Random Fourier Features: A theoretical and empirical comparison. In Advances in Neural Information Processing Systems 24, 2012.
[69] Enoch Yeung, Soumya Kundu, and Nathan Hodas. Learning deep neural network representations for koopman operators of nonlinear dynamical systems. In 2019 American Control Conference (ACC), 2019. doi: 10.239 202.3999//CC.21999 8.3399999.999
[70] Vadim Yurinsky. Sums and Gaussian Vectors. Lecture Notes in Mathematics 1617. Springer, 1tt edition, 1995.
14

---

Appendices for: "Estimating Koopman operators with sketching to provably learn large scale dynamical systems"
A Setting and notations
A.1 Operators and notations
We define the following operators:
Φ I > H. defined by Φx = fx f())(()) dz(x) for any 1
2x H. + L defined by 2x = (h, o(-)) for any h E H (i.e. the embedding operator mapping function to its r-equivalence class L2).
Y : - + H. defined by Φylx = 1x
Yix H L2z defined by Vixx K2 Dx
C: H + H defined as Ez~=v(()) @ o(x) & x satisfying tr(C) S K2, Note
that under our assumptions, this also corresponds the covariance of Y.
Cxx [(a.))~p(1)) @ Φ(v) Tx "5 x
As well as the following discretized variants:
x I + H. defined by 2xv E (2)(2;) for any [ui....,
.a......,,, Val E R"
x H. + R" defined by i [(0(x1), h) (Ō(I)). n) for any h H.
o _ + H. defined by v -xxx==E=1 (p(())) for any (u1...., vnl IR".
Yix H 1R defined by 4 (((()) (0 ..... (0(yn), H) for any h E H.
= x = = Σ o(x,) 2 o(x;) E L(H) is the empirical covariance.
The Nyström discretized operators are obtained by applying the kernel map to m < n inducing points ( (ar,) and {ū {y
x
x w ΣØ (),( R such that क
. Ty R" + 7 such that T = Σラ= ())(ū;)
Furthermore denote by Px and Py the orthogonal projections onto span कx and span Ty respec- tively.
One important quantity derive the rates is the so-called effective dimension, defined
= = CC).
where C₁ C
A.2
Conditional mean embedding
For any I E X' we denote (1) the conditional mean embedding associated to the transition kernel defined as
Mp (x) = E IXx = x] =]= / d(v)p(z,dg) dy)
The following lemma provides characterization of 4 x in terms of the conditional mean embed- ding.
15

---

Lemma A.1: We have the following relations:
Dylxf / (()))) fe La
(16)
( (x) = (Si1((z)), JEH
(17)
$ T1xx Exwr(p(() @ 1p(x))
(18)
Proof of Lemma A.7
For the first property:
( (1()x))(z) (Xx(f)))
(19)
/ (ox(l())(odd
(20)
= f f(y)p(zdd))
(21)
(f. dl))((z, dy)) = (f. up(x))
(22)
where we used that f and & coincide π alalmsst everywhere. The second property is a direct consequence of the definition of the adjoint. For (18), simply (17) and the definition of Yyix to get
e (3.p(()))4)))))
(f ((o(t))p(())d dr(z)
A.3 Power spaces
We now define the a-power space [H]] in order provide some intuition regarding Assumption 4.3
By Assumption 4.1 tr(C)= / tr(((=)) @e(x)) dı(x)) < K2 and thus C' is trace-class (and compact).
By 18). there exists non-increasing summable sequence (pi)ier for at most countable index
I, family (ex)er E H s.t. (i)ee is an orthonormal basis of span tx S I and (14) e)ier is an orthonormal basis of (ker xx)+ C H such that
C = Σ meki, 1/2 ee)uu
1/3 4z
IE!
For a 2 0, we now define the a-power space as
[H]
Σ
a:pi
(at)ier E E2(1) } ELL
REl
equipped with norm
Σ
ampe
e
[H]
We can now make the following assumption regarding the embedding of the power spaces into L. Assumption A.2 (Embedding): There exists 7 [8, such that Cy || [H]] + 1
We stress that Assumption A.2 implies in particular Assumption 4.3 and is common assumption in the literature, for instance 18]
B Expression of the risk
We have the following risk decomposition.
16

---

Lemma B.1: The risk can alternatively be written
RAs(A)
E(z,))~~loo(0 - Ad(x))1
RHS,o EAs(AA
where
Rus = | 1 llAs
/ [|up(2) ((y)11 dd(x.
and (A) |0xix AΦx llís
/ [1)p(x)) Ab(x)) dr(a).
where inf (E)S( EHs = O, and thus we interpret Eus as the excess risk.
Proof of Lemma B.1 Let (hi)kEN be an orthonormal basis of H. Then
Eus(A) = 1 A liAs
Σ
/ Ah,h,,1122
iEN
Σ
(A",,, ()))) dt(x)
iEN
Σ / (1hr - (hs, Ao(=)) dT(
(by 17
ieN
/ [I)p(x) Ad(x)]1 dd))).
holds
/ /11p()) -(y)))
RHS,0
-
- 2(1p(()) / (115(())1-22((()))_)))) tr =()())(())) dp(x,,)
/ (1)) dr(x)= / (40(0)) (ow)l(e,d)))
dr(x) + / tr (lv))(())'(( dy) dx(x)
H
(IIA(())1 dx(z)) (T(6(v))())) ar())
tr ((uost))pla dr(s)) + tr(()
(x1 + tx(
17

---

where we used the invariance property of π in (2) and Lemma A.1 for the last inequality. Then one can easily check that the sum of both corresponds the full risk defined in (
/ [11)(()) o(v)ll dp(a, v) + / /llpp0o) Ab(x)) (
RHS,0 + (A)
(())(z), d(u)p(,, dy)) + | (( plol)) dw) dr(x)
+ / - ((4)))))))22))) Ad(x)) |A(()))? dx(x)
/ ( /10(w))Pp(s,dy)) 2( d(v))(s, Ab(r)) [|Ad())11²) dz(x)
(16(w))12 - 2(0(v). Ab(x)) + 1 [| dp(x,
/ (100)) Ab(=))1 dp(e. (A)
C Expression of the estimators
In this section we give proofs of propositions 3.1 to 3.3 on how to efficiently compute the Nyström estimators.
For all three - KRR, PCR and RRR estimators, the starting point is their respective full estimator which can be derived by following the first-order optimality criterion for the following minimization problems
Full KRR:
İRR arg min AEN-H
A&x
1iEs + A|
(23)
Full PCR:
APCR
arg min
Prix
4 liis
(24)
AEH
Full RRR:
RRR
arg min
A& 11 + A|AAS
(25)
where II, is the orthogonal projection onto the top-r eigenvectors of Ć.
To derive the Nyström estimators, we project the embedded data Φx. Φvix onto the span of the embedded inducing points Px Ωx. Py P and then express the resulting estimators as Φ W with W Rm This form is particularly useful for later computing forecasts, eigenfunctions and Koopman modes with the estimator. In particular the following equalities for projection (shown here for Px but equivalently exist for Py)
= Ex Ex = 2 = x
and the characterization of Px through the SVD of कx U'V*, that Px = UU*
C.1 Nyström KRR
We begin with the Nyström KRR estimator, providing an alternative but equivalent description in lemma C
Lemma C.1 (Expression of the KRR regularization): Let U be such that Px = -U0UU0 Then holds
9KRR (KRR(C)== (Ō)== Px(PxCPx =
(26)
18

---

Proof of Lemma C.1 Using U 1, holds (U*ČU + XNY = -UU-CCU + A7) and
thus U (UU* (UU-CUU CUU
AT) = (U* + - As consequence,
(C) Px(PxxCxx
=UUM(UU-UU" +
-U(U-ČU +
Then we can provide the computatable formulas for Nyström KRR
Proposition C.2 (Nyström KRR): The Nyström KRR estimator, obtained by projection of eq.
23
AKRR
PvCxxPx +
By K K ky x.X X nAK xx
Proof of Proposition
C.2
Using the definition eq. 26 and lemma C.1
we have
A KRR
PyOy (C)
m,
+yCxxU(VUUU+MIUUU
-UCY -
Now using the fact that E, V. and ČU +1 are full-rank, holds 2 eq. (20)
AKRR = A
= Px ČUEV" " -
Finally, by definition of Py. Cyx and C
ÁKRR
x ( $y ) 10 4
2 xx ( xx + xx)
कy k 4,
K
xx mAK 16xx.
Remark C.1 (Alternative derivation of the Nyström KRR estimator): Note that the Nyström KRR estimator can equivalently be derived as the solution to a variational problem similar to eq. (23). where the operator is restricted to operate between spaces Hx := span 4 and Hy == span Dy.
C.2 Nyström PCR
Define the following filter the spectrum of PxêPx: (CC [PxCPx] which truncates it to the first r components before taking the pseudo-inverse. The Nyström PCR estimator, obtained by projection of eq. 24 is
PCR m
(27)
The next proposition provides an efficiently implementable version of the PCR estimator.
Proposition C.3 (Nyström PCR): The sketched PCR estimator eq.
satisfies
ACCR Ty KY y. Kyyy x [K] R.x
K Kx,
K
1x
(28)
19

---

Proof of Proposition C We begin by computing the decomposition of PxCPx which necessary to obtain 9PCR (C)) The following expressions are equivalent 43 Proposition 3] for determining its eigenvectors and eigenvalues A:
U Ah
Ú'ÙUn= = Nh,
h Uh.
Let the truncated eigenvalues be A, -iaddig((11,---- Ar] and the eigenvectors be H, [M -........ -.. Then H, UH, must be normalized such that H1, HIU*UII, I. The rank-r truncation [PxCPx] is projection onto H,
[PxC = -UU
-UHAH HATHH",HHÜUR H,HUU*) UH,A,
where we used that U-ŪU = HAH"
Now substitute U 4xVE-11 to simplify the eigendecomposition of Č'
Σ- X k VV VE h Ah
VE 2 xx x d= Xd,
h= = xx xx
(29)
where 1/--2** Kt x,x Denote by D, = .d....., dz] the truncated eigenvectors such that H,
E V x x x.8 Dr, normalized such that 11 H = D K xx xx K1 K K x ==1, =1
UH, HH; = = X.X K xxxxx K xx x D.A," D Kx XX
= D
[x[ x x K K x1.
Finally, we can plug the pieces together get
PyCyxIPxCPxl! = T 1 Ky K
[K] 8 xx xx,
Remark C.2 (Variational problem for Nyström PCR): Note that, unlike the NysKRR estimator, the variational problem for NysPCR where the operator restricted A H + Hy is not equivalent to the one obtained in proposition C.3 by projecting the covariance operator. In fact, the former does not take the full covariance into account when computing the low-rank projection, but just the Nyström points.
C.3 Nyström RRR
The Nyström RRR estimator does not correspond to a specific spectral filter. We can nonetheless compute it starting from the expression of the exact empirical estimator (29), projecting the covariance operators, and rearranging the expression to result in fiitteddieeniionll procedure.
Proposition C.4 (Nyström RRR): The sketched RRR estimator can be written
ARRR = PxPPPPP'
X1)
(30)
To compute it, solve the m eigenvalue problem
Kt Xx, + +nAK Kx.x y.y Ky K = o?w;
for the first T eigenvectors Wr = I. Then let D, = KJ K K x W, and
[W wel. normalized such that
W W x x...Ky Ky Kt & Ky
E, = (KxxxXxxx + -nkK m
xxx Ky KyyuU such that the following holds A,
D, D,
(31)
20

---

Proof of Proposition C.4
B n
1 CxPPxxxxxPPPPP X1)-1/2 The computationally
intensive part for this estimator is in evaluating the rank-r truncation [B]].. Its singular values and left singular vectors can be obtained by solving the symmetric eigenvalue problem BB* qi = 0.q. We rewrite DB
BB = P
P x x
= P Px xx n\1) - Pr
Pr
X.N Kx,x -nd) - y
= =P
+ -ndd) K Xx,xxxx 1 x
Py- x
x. + x )) K x.x di
where the second and fourth equalities are applications of the push-through identity, the third by definition of projections and kernel matrices, and the last by collecting K Xx By construction, the non-trivial eigenfunctions of BB are in the range of Py x therefore we can set q. P xxxx W; for some w, E Rm, and solve the following eigenvalue problem instead Py xix K x. (K x xxxXxxx Xx nAK & x) K & V Py K x =o: W (K X nAK x x) x yy K/, y Ky xx 4 ofw;
where we have simplified the left term of both sides of the equation. The eigenfunctions of BB* are therefore qi Py K x. Wi, which must be normalized as
|1 wi K x Ky K Ky x =.
Thanks to this normalization, the projector onto the leading left singular vectors of B is Q,Q;, where Qr .q...... 4r). Then the NysRRR estimator can be written
where
+ + - -
= = x (K x + x
with the same techniques we used for rewriting BB Finally, let D, and E, in the statement. We can apply the projection to obtain
Q,Q;B(Psxxx n Ty D,
D Forecasting & Koopman Modes
The three estimators considered in Section C are all of the form
Å by W
W Rm
We will use this generic form to provide expressions for the following operations:
1. producing forecasts of the dynamical system a future time,
2. computing the approximate eigenvalues and eigenfunctions of the Koopman operator,
3. computing the Koopman modes.
D.1 Forecasting
Given a new data-point E X' and an observable function y E H (note that this can simply be the identity function), approximate the me-step-ahead expectation E x
21

---

(x) by using the obtained estimators A'. Note that by the reproducing property Tig
=: 9m. then
(1) = (कx W ( -
W Im (x)) Σ (WT gm )
D.2 Eigenfunctions and eigenvalues
We wish to compute the eigenfunctions & E H. as well as the cigenvalues A Å. The left eigenfunctions satisfy År;; : = X and the right eigenfunctions satisfy Asi = Agt the following
will use Mollenhauer al. 43 Proposition 3] manipulate the eigendecomposition of operators in H..
Consider the decomposition W = U.V. with U,, V. E Cm which is available for all considered estimators with m. For example, simply take U, D, and V,, = E.. For the Nyström KRR estimator instead, m and we can take the whole W as our U, and V, =.
the Nyström RRR estimator proposition C.4 we can
To compute the right eigenfunctions vhs such that (UyU,V x3x Jobe Aeth, consider the following equivalent eigendecompossitin
2 Uy = Xrg. where vi; by
Note that क = X a miite-dimensional object which can easily computed. The eigenfunctions uh must be normalized such that v/ eh for every i, so we must have
9. UU U.
very similar process can be followed to obtain the left eigenfunctions &, such that x VyU, drg, Here we consider instead
U क V.h, = Xyhy, where &, = V2h4-
U
where once again, कx KP and the eigenfunctions must be normalized such that h; v. x 1 for every Finally, & and & must be orthogonal to each other: must have for i E [r] that (1/s Ej) = 6ij (where 6i is a Dirac delta equals to 1 when i j and O otherwise). We can compute
(ws, Ed)x =
>yhiāh
and note that h0 8ij. but must normalize such that
& xx ViVilu/
D.3
Koopman modes
Given the eigendecomposition of any estimator Å A, -E dets @ Es. for an observable g we have the following
Å X (g, 3i)H
where (g. dinu 7 are the Koopman modes. Expanding the definition of s/ we get
7 = (g. vi)n = " U. U, 9m E C
which we can efficiently compute.
Excess risk of the Nyström KRR estimator
E.1 Almost-sure decomposition of the KRR excess risk
22

---

Lemma E.1 (Excess risk decomposition in operator norm for KRR): Let Assumptions 4.7 to 4.3
and 4.5 hold. Then the Nyström KRR estimator
8 satisfies almost surely
E(AM)) & a/ + adQII(C' - Cx)C 1/2 IBA - (f Oxx) 11800
+ a0,020021 P /2 118//1 + P |18870
where 01
1/2 2/2/ 02 = 11 - 1/ 0g = ||0 Cxll. and a is the constant of
Assumption 4.5
Proof of Lemma E. Let 0, == 1|0 0g == 110 As in Lemma C.1
define
9KRR
(C) ;- CU + A1) -U" We have
E (AAMRR)
|1 - AKRR 1 x (BB(2,
|0p1x A (18(11..
(A CYXOKRR (x)0xxxx111)
|1( x IKRR (c) AKRR) 2x (l
r A x (1813,1
[AA XY -CYx
||(CYxxKKRR (C) AKRR ///31
(32)
where we used the polar decomposition 8 WC for some partial isometry W : H - L2 The first term is
f - xx = =1 - Cyx xxll ||
sax + |[[ Pk) Pyjxl
where we used the definition of Øylx and applied Lemma
The second term of our decomposition 32) can be bounded as follows: It holds
B == 9KRR (1))
(by Lemma H.2
I Sal|C(T 9KRR (2))
Sa |0[0x* : +||C(CT 9KRR (1))
We now bound the terms B₁ and B2 separately.
B1 = =I[C(,,1 - :j 1
-I|C0 - CG)
|100, III('s 1//
1/2 ||
s Of][CC - 1/2
Let P = 1/2
OKRR (1)0 We recall that SKRR (C) = PX9KRR(() so that
p = C22² (OKRR( (()d)p)))KRR(()C)))))
0//2 P (C)))R(())))) (2)
Ps.
23

---

This implies Px = P3 Px. Hence Px is an orthogonal projection, and defining /x P3
holds ||2x [B00 - 1. We can thus bound B2 as follows:
B2 ||0(0= 9KRR (C))((()))//1111
-I- - -I
(by Lemma I.1
-1
- [100 (1000 11scHol1C 1
1|Bolll 1/
|1
(3(70))
01020s1 |
For the third term, due to Lemma
C.1
C x P (C)CA/21
(33)
s |lCyx Pr 1/2 12 1110 1/2
|| PAllscHolI 1/2
|
|B00
(34)
s ff - -y 1/2 |8000
(35)
- o8 ( B|B0 + ll(yrx Oxx) (1(mm)
(36)
where we used Lemma J.7 for the last inequality.
Starting again from 32) and putting everything together, we get
((ARRR) 1/2 S ax1/2
aBfl(C's C.)C, 18H0
+ ffl1(rxx Crx)ox 1/2 (1(0
+ 00,020s||PPt
|1000
+ 0 // (1()))
E.2 Excess risk rates for KRR
In order to control the terms appearing in our decomposition, we recall that Assumption 4.4 implies
B <
daa(A) I C3X- where C3
1-8
B
(37)
K²
.8=1
where e is the constant of Assumption 4.4 see 9 Proposition 3 with + 1/3 and 8 c] and 18
Lemma 11] which shows that the existence of a constant C3 such that the first part of (37) holds i-1/
implies in return Ac(C)
Proof of Theorem 4.6
By Lemma E.1 taking Px = Py. it holds almost surely
E(AKRR) 1/2
1/2 aB? [[(,s - C2 1/2 |()) + @f - Oxx)Cx 1/2 1
a0, 100020011 Px / |8:0 + 0.11 1
|1
and we recall that 01 == 1/2 ///1// and 0g ||0 Call- We bound separately the terms appearing in this expression.
Bound of 0 and 0z. We control these term by bounding |1 1/2 C)C By Lemma J.3 it holds for any 8' Ejo, 1[ and any EjO, [[C]|B(]]]] with probability
c (C--C C)C 1/2
8
+
22, πA=
where
log
8K2 8'A
(38)
3nA
A sufficient condition to bound the right hand side of the previous expression by 1/4 is to have
nx > 32c, (in which cases both terms bounded by 1/8). Assuming this holds, T 1/2 (C-
24

---

c) is invertible and we also have
8 = - c 12
10, - 0 - -
$ =1+11 11+ - (C)),,
1.25
and thus
0g S 1.12
while 89 8 -|=
1/2
1/12 [(O) 1/2 1/ CxCx 112)1111-111111
(1)
1 (1 |0x 1/(((2 -C) 1/2 1)-1
1.34
and thus 01 $ 1.16
where (i) can be obtained by taking the Neumann expansion of 1 - 1/2 (C- 00-0)
Bound for 03. By Lemma J.5 it holds with probability 1- 8'
2K (og((/6)) A
2K ExcttrCC (C)
IICC - ( ol S
(39)
Both terms in the above rhs are bounded by 1/4 provided
\ 8K log(2/6"
n 32K 2 (1++) 1og(?/8)
where we used [0) + the case, we have ||(C CCC 1lpp 1 1/2 < and the operator 1 - -Ox invertible.
Cgt- When this
0g ||(},
(|( (Ōs -x)
(0)
$ (1 |1(0, - -1
2.
where (i) can be obtained by considering the Neumann expansion of I- -C2 -CA- -,)
Bound for | By Lemma
J.6 provided A ejo, [[C]|s[o]] it holds with probability 5'
|
C
L/2
|B00 I V3A
provided m 2 max (7, 5 ess sup zow ||
1/2 x(x))1
ax² log Adr which by Lemma H.3 is ensured if
m 2 max(67, 5 log 4K2 A
Bound for |I(c "
and I| CC, 1/211-
By Lemma J.4 for any 8' Ej0, 1..
each of the following events holds with probability 28':
2K /Cr log(3/8)
2K log((/88)
( - 1/2)
+
Xr/2n
(40)
By Eq. 37 we have (A) S CaA-3.
Choosing d' 8/3, we get via union bound with probability 1 that 0,0203 < 2.6, 03 < 1.34 and
2K Ve, log(2/8)
2K log(2/8')
((AKRR) 1/2
< 2/2 + 1.34(a++))
+ (2.6a + 1.34) 332
XT/2
nd
S 1/21/2 + 62d TT2
+ 1.0001111-00 -B/2
where:
= (5.5a + 2.33)
== 1.34(a + 1)2K /Cr log(?/6)
3 1. + 1) (ok3(glog((//))
25

---

for any À and m satisfying the constraints
8K4 1/7
À n (/2 log 8x
A 2SK Cc (G)()///()))))
A n 11-1/))))))) (32K* log (og((/)))) ((4)))
(41)
A 1
C7 log
4K A8
m 2 max(67,
(uniform sampling)
We pick A cxn cxn 11(+8)))) which is asymptotically the saturating constraint (given that 1/(1 + <1 <)<1 < 2/ 2 1/T), where is a constant choosen to enforce the following equations
(which are sufficient conditions for eq. (41) to hold):
n (1=T1/() > (32c,) ) log
8'ex
(r+1)/2 n -1 1)/((2)))))))) > 8K c log(?/8")
(42)
2 (32K² log(2/8)) 1((+3))
Cxn c2n-1/1+0) K²
As - (T 1)/(2(1++ B)) > 0, a sufficient condition for the second equation is
cx 2 (8K V& log(f/+))3//t+))
Assuming 2 8K2, sufficient condition to satisfy the first constraint is
exn (n11-/(1+ (32c, ( (n (1(((++)))
log (18)-1))
which is in particular ensured (noting that log(n)/n" - 1/(e)) for any n, > whenever
ex > (64e, (aa((aaxx(e((aaaeee 8 - log log((/)))))
Noting that + S get that
cx > (64c, -)) (og
also sufficient. We recall that 1/11 <
< 2/ + 1) - 1/T, so that we can choose
cA log(?/")) max((3233)1((493) (8K //7
while the last constraint (c2/K2)1+) satisfied by assumption.
E( ((KKRR) 1/2 S 1/2 + 02 1/2
C3 - 8/22
$ sac 1/2 12 1/12(1++))
+ C2c, 7/2 7/2(1(()) + C3C2 8/2 2/221222//111
1/2 1/12(1+8) 02C,
7/2 (1-24(000)0)- + c3c, 8/2 -/2 11(((1++)
< c
< (c + 02C T/2
1/2
6/2 11(2(1++))
c3C,
which gives the claimed result. The last constraint (on m) is satisfied by the assumptions of the lemma.
26

---

Excess risk of the Nyström RRR estimator
Recalling (30)
NyströmRRR estimator is of the form ARRR
[B1- where B
m.d
Cxx(Cx) for Črx X and C2 PxCPx AI. While the population ver-
sion is RRRR_ == [ 1/2 where B = Cxx((C))-112
In this section we follow the approach in [30] and decompose the operator norm excess risk in the following way:
A 1/2 AB +1 ( 11 - x |B(L2
m.
Then, recalling that A CYxCr' and AKRR BCJ1/22 also have A РвΛΛ and ARER PB- AKRR, where PB and PB are orthogonal projectors onto leading r left singular vectors of B and B, respectively.
Thus,
E( C(ARRR)1/2 1/2
2 a\ )- + [R
A m, 4 x B(R
= a2 + 1 + |((BAA - P - xx B(H)
sa2 00-1((91xx) + [[(PBB -PB) $
BSCH) + (1 AKRR f ) ([1(C00)
$ 2//2 + 0 + K |IB - BB B(H + |I(AA 2x
o2((B) 0 +
where the last inequality is due | S and 30 Proposition 4]. Recalling Lemma we observe that
|BB 1 BB"| oZ(B) o
E( E(ARRR) 1/2 S ( K
+ 1/2 + 1(Ax - ÁKRR Φx
Sr..s. of the bound in Lemma] E.1
Therefore, to prove Lemma 4.7 for the RRR estimator we just need to bound | BB* (B8H)) To that end, observe that, after some algebra, obtains
B B.-
and, consequently,
[B BB S 2all - -yxllaxn + a21PPxXPx ClIBCH)
+110,1/30,0 B
| /2)
follows using that ||A Sa.
On the other hand,
Crx) ||Cyx Cyx (I I -11B)/Cxx-Cxx)xxxleee +
+ +|Crxxx
which implies that
|[xxx- - Cyx & 1lCyx - B0H + 2 a K E1.
where E := maxx[1|p C2/2 18H) ||
| ( Similarly, we obtain
||Cs - CAlB80 < | Cl| + K E
(43)
But, E can be bounded by Lemma J.6 Indeed, provided À EJ0, |[]]))))))) it holds with probability 8'
E1 S V3A
provided 2 max 07,56 5 ess sup sup || 1/2 /1)0(() log 1 28 which by Lemma H.3 is ensured max(67, 5 5z) log A
1K2
27

---

Additionally,
||0 (0-1122,00 0,,00 1800 1 PACPX 1/2 (1BC) + X P 1/2 (B80) 1870
s PX
+1
P. 212 +1
s 03 1011 Px(xx(sx)) Px 1lB00 + 1,
implies that
| "0,,""0,00''1 1/2
B00 S 02 + +15115 2.25.
(44)
provided, as above, that n2A 32c,
Therefore, setting Eo max coxx/allCxx Cyx(lam) a²I1C Cl for all i E [m] we have
[B) - - s ||BB BB*|| 3eg + 6.93 K a² 2 + 2.25
(45)
where E2 = ||A - A C /2 the variance of Nyström KRR estimator, and conclude that
E(A E(A 1/2 S 50+41(911x + K 35D + 6.93 K a² X + 2.25 + a 1/2 + +
dB( o
Therefore, the proof of Lemma 4.7 for RRR estimator directly follows from the bound on a X 1/2 + given the proof Theorem 4.6 and the fact that, see e.g. [29]. Eo S n 1/2 $ X1/2.
G Excess risk of the Nyström PCR estimator
Recalling Eq. 27 NyströmPCR estimator is of the form
APCR = [ = AKRR Pex
for x 0 and with Pús being the orthogonal projector onto leading eigenspace of Ö2. So, prove Lemma 4.7 for PCR estimator, denote APCR := AKER Pcj for À 2 0, and let us define the population version APTR_ A where Pc, being the orthogonal projector onto leading eigenspace of C)-
As in the previous section we start with decomposition
E( E(APCR) = 0 A2 (B + 11(AA
APCR J (IB((BLL
[[(APER
A
APCR 2x (BBLX + [[(APCR APCR x (xlB((2).))
The first and the second term easily bounded by || A
S s/ and
II(Ax ) o 1x115122117 = ||Ax - Pc, } (1(12,00
S a [ Pos Cc 118000 S (xx))
For the third term, start by observing that
[[ PCR
APCR ) x
= II(A,Po AKRR Pú; /C1/211B(8)
< ||x, -Po, (23))C/2211860 + |[AA
)
C/1/2 (1000
(180)
< K PCccs - Po, (1ECH) + (1(AAA AKRR)
Pc,
1 ([H00)
-1/2
|[O' - CA (BH
< a K
+ 11(AAA AARR
Pex
C1/211
-2 oZ (xx)
Eo + 2a K E1
<
+ 11(Ax
AKRR) (2/2((BBH)
( o (xx)
28

---

where the second last inequality is due to 30 Proposition 4] and the last one uses Eq. 43 Moreover, we have that
IIAA A 01/2111000
11Ax A C 1/2 1|8(H0 |1 0/ 110000
|1x 12pp C3/22
S E2 V 1 + 03 C (121211/111000
82 03 |1 |BB0,
where we have used Eq. 44) and the fact that Pa is the spectral projector of 0/1/2 Therefore, due
1/2
| 12
1/2
BBCH
1 + PxI0, 1 & || Px
xxx11/1
BO) +
PX
||00
0211C Px
1/2 1 + |1 P (anoll(; 1/2 (8(0)
02|ICx Px 6//2
1/2
||C
Px
B0H) -1/2
B0H
02 + =1
we obtain
[I(A) - AKRR P Pa. C1/211/rH)) I E2 (1.68 + 1-5x-1122 11)
provided that nX > 32,,B.
Finally for the last term, observe that [[C]]] and [Co]l share the same eigenvectors, and hence [Col] [(SI X[OxCo]f. Hence, holds that
[[A APCR) x (15(21112,1)) =
[C,l)) 11a0 = X1|Ōxx0 [[ol " B00
S O l |a 1000 1
|18000
" 10, //2 |15000
1[166]t1l)))) [ 110,
1/2
c 1/2
||HH))
Now, recalling Eq. (45 we can bound
| Bib S Mx 1/2 |BB0 + |BB" BB (|
s a²²k2 + 350 6.93 a² 1/2 2.25
and,
11/2
11/2
11/2
X4/2 [ = (xx(Pxx Ar CCas - CXllB0H)
02((xx) - 1 - Cl B(H)) K E1
Thus, consequently, obtain
[[P m,
x
x (B()) 1/2 (02 + (12-1111 (a2K² + 360 6.93 K d²x1/2 2.25
11/2
x3(4x) |1 C((B(K) 2 K
To conclude, observe that r >> n काोम due to Assumption 4.4 implies that 0 (1x) S n 2ーーM..
Therefore, collecting all the terms, under the assumptions of Lemma 4.7 we obtain
E ((APCR) S CPCR n 20+3)
where CPCR (0,(0x) of the problem dependant constant.
H Auxiliary results
29

---

Lemma H.1 (330, Proposition 2] with a = 1); Under Assumption 4.2 it holds
x Cyx xx1 - 12
Lemma H.2: Let A be a bounded operator: Under Assumptions and 4.5 it holds
||Cxx S a||CA||
(46)
Proof of Lemma H.2 Note that under Assumption 4.5 as Oxx & a²C² it also holds A x AC2, and thus:
All
= =11AAOrxCxxxI]]
5 C²
= =a|CAA|].
The next lemma is a consequence of Assumption 4.3 and will be used in our concentration inequalities.
Lemma H.3: Under Assumption 4.3 it holds n-llmost surely for any
| (1-1)/2 x(x) || S c-x-1r-v14 22
The two following corollaries can be obtained picking and -:
c 1/2
1| S x and
2 | 6(x)) ||
CT
S
X-+1
Proof of Lemma H.3
By 18
Theorem 9], it holds o, |A: 12 sup Exer pf [ex(x))2
(where (e;) is defined in Section A.3 and recall that (Vite;)iEN is orthonormal basis of H.
Denoting μ4 == X (C), it holds
Σ iEl
(i A) (1-(11 (Vpres) (Vites) d(xx)
Σ GET
(1)( (14 1)
Σ ut T + +A) 144 11 - 1Te((z)²
(Mi (mi+))" me((())
A
iel
< Σ
(14 + A)
(T-
ie/
S -- k (1(
where we used sup 1ud S R2
30

---

I
Deterministic sketching results
Lemma 11: Denoting := 1 0/2 IKRR (C)O2 it holds
ROX² RCX PX.
Proof of Lemmal I.1
This is direct consequence of the fact that (xRR(C)(OP) Px:
9KRR R Px -222p Px - -2 9KRR (C)CxPx = 0.
J Concentration results
J.1 Generic concentration lemmas
All our concentration results derive from two versions of the Bernstein inequality. We first state an inequality for sums of random variables in Hilbert space based on 770 Theorem 3.3.4], which itself derives from a result of [8
Lemma J.1: Let (At) be i..dd copies of a random variable A in a separable Hilbert space (H, I1-1)) Assume EA p and FF 0, 3. 0, Vp 2 E|| nll I ploa TP-2, Then for any
E00,11 [ holds:
2L log(2/6) <
20 log(2/6) +
A
2 21-6
(47)
n
n
n
The assumption on the moments holds in particular when A s L/2 and ||[||| < o².
Proof of Lemma J.1 See proof of Lemma E.3] for precise derivation based on 70 Theorem 3.3.4).
We now state version of the Bersstein concentration inequality for self-ajoint operators in operator norm, which a restatement of [34] Lemma 24]. In the following, we denote refr( (A) := [||]|||||] the effective rank of nonnegative definite operator A.
acting on Hilbert): Let H be separable Hilbert space and A; be i.i.d. copies of a random variable A taking values the space of self-adjoint operators on H. Assume EA = 0, ess ess Alp < for some >> 0 (where Il-llop denotes the operator norm) and that there exists a positive semi-defiitee trace class operator such that E|
Lemma J.2 (Bernstein for self-ajoint operators
V. Then for any ó e)0, 1 and 2 it holds
2cB
2
B = log (4rg(V))
P
A
+
< where
(48)
3n
Proof of Lemma
J.2
See
34
Appendix B.7, Lemma 24].
31

---

J.2 Applied concentration lemmas
Lemma J.3: Let Assumption 4. hold. Let 6 EJ0, 1. Then for i..d samples (2i, Mi) and any
Ej0, B(/))) it holds
4cT,
20, 3
P
where
8 == log
SK-
(49)
3nX7
Proof of Lemma J.3
We apply Lemma J.2 on the random variables A, = E(Xi) @ E(Xi)
where {(X,)
CT
1/2 (X,)). It holds
ess A, BBH I 2 ess sup supplc (X,
2c,
(by Lemma H.3
d
B[A
<
B(IIE(X)))
((X()))((X)))
x-
CC
x-
Thus applying Lemma J.2
with c =
2c. x
and V = =#00,
we get
P
A
4c, 3
+
88
where
B = =loog]
SK3
(50)
3nXT
X'n
where we used the fact that ||CCT| BO) S and controlled the effective rank using tr tr(CC K2/> and ||CO BoB0 = ||C|| (ICC///BB + 1) 2 1/2 because X I ||C|) by assumption.
Lemma J.4: Let Assumptions we get
4. and 4.3
hold. Let 5 EJO, 1[ Then for i.i.d. samples (Ii, vi)gggn
P[I[C- - 1/2 1/2// lop S sc(,,
2 2-8
(51)
and P[l(Crx- M -xxx) 1/2 1mp S ((1, 8)
2 21-8
(52)
2K /cr log(2/)) XT/2n
2K den (A) log (3/6)
where c(A,8)
n
(53)
Proof of Lemma J.4 We first write the proof for the eq. 51 For this result, use the fact that [I - o Boa & [[ 1/ and bound the Hilbert-Schmidd norm.
As (1-())),lllllsss)) is a Hilbert space, apply Lemma J.I on the random variables A, ((re)) s8 E(x;) where {(x) 1/2 1/
ess
= csssupllo@))lll(())))
< KVG
(by assumption 4. and lemma H3
\1/2
E A||As| }[
(dn
32

---

Thus applying Lemma
J.1
with
Kyer
and o² = don()) gives Xt/2
2K /c, log(2/8)
2K 2K log(3/)) 6.
P
A
n
17/2n
n
HS
This yields the desired result via the inequality between operator and Hilhert-Schmidt norms. For the bound eq. (52) on the cross-covariance, we take A, = x(y:) @ E(x,) but the rest of the proof inchanged.
Lemma J.5: Let Assumptions 4.1 and 4.3 hold. Let EJ0, 1s Then for i.i.d. samples (I;, Ik yilisin get
2K /c log(2/6) 1
2K 2K² (r(CT() log(2/))
P [I(C - C)
2 21-1 8.
(54)
Proof of Lemma J.S
For this result, we use the fact that II(C C I1800 S [(C ((C-O)0x" and bound the Hilbert-Schmidt norm. As (HS(H), (l-)ss)())) is a Hilbert space, apply Lemma on the random variables A; = (x(x)) @ wxx;) where w(x) Cx o(x))
HS
J.1
ess supl Allss ess supllo(=)]]lla///
S K K (by assumption 4.1 and lemma H.
X
E B||sA| (x)l (|)()))
S K (
= K tr2
Thus applying Lemma J.1 with L
K K and 0² = K2 Cr(C gives xc+1i
2K 2K V log(3/)
2K Ir(C log(2/6)
p
A,
2 1 - d.
HH
This yields the desired result via the inequality between operator and Hilbert-Schmidr norms.
J.3 Probabilistic inequalities
The following lemma is restatement from 153 Lemma 6].
Lemma J.6 (Uniform Nyström approximation): Let Assumption 4.1 hold. Let P: -> H denote the orthogonal projection span Φ(ā,) $ m }. where the landmarks (#;)1sism drawn uniformly without replacement from the empirical data. Then for any E]O, ||C') ((())) we have
|P 12 - 3A
with probability at least 1 provided
2 max(67, 5 ess sesss 1/2 ((x))(1))))))))))))))))))))))( (1) log 4K
A8
33

---

J.4
Concentration lemmas for the sketched operators
Lemma J.7: It holds almost surely
Il 121 + -
Proof of Lemma 1.7
It holds
Cyx - CCxx- -Cxx-PrCxxx
=
Thus
[I -xCxx) CyX Py (xxCxx-Crx))
1 11 -xx)
$1 "1 |l(yyx - -xx)
Eventually it holds |00 1/2 CyxC xxxxx1122 1/2 & 1. Indeed, as π is invariant, it holds that
/ (())(t, dx(x) 1 fErg:l/f
||
sup
and denoting 4x C1/22 the polar decomposition of 4x for some partial isometry I + H. and using +6 xx we get
|| 1/2 112111xxx, Cxx
||=
x 122
< ||1,U,1,,,,,
< 1.
J.5 Concentration for mixing processes
Lemma J.8 (Kustic et al. 229 Lemma 1)): Let X be strictly stationary with values in a normed
space (X, II-II) and assume 2pk with p, E N. Let Zı,..., Zp be p independent copies of
Z1 = Σk= X,, Then for s > 0:
X 8
2P
Σ 2n
s/2 + +(p ()3)
j=1
K Additional experiments
In this section extend the experiment of Figure (on the Lorenz' 63 dataset) to include the Nyström KRR setting. Since the KRR estimator does not use rank truncation as a regularizer, it will always output the full rank estimator. In Figure S we use the same settings as in the main text (m 250 and r = 50), which lead to KRR being very close to RRR. No particular instability is encountered with these settings, and KRR is about 2x faster than RRR. In Figure 6 we inerease the number of centers, without changing the amount of regularization. Now KRR is very unstable, since it needs to estimate larger number of eigenvalues which have a very small magnitude. RRR on the other hand maintains the rank regularization and has no instability. obtaining a small performance boost thanks to the increased number of centers m. To improve the regularization from X = 10-1 to X 5 10-3 However, notice that a) the instability is not completely fixed and b) accuracy noticeably reduced to the point where it's nearly the same that of PCR with fewer centers.
stability KRR, Figure 7 we increase
34

---

1.2
1.0
sKAF
nRMSE
NysPCR
NysRRR(M-2250 M 250)
RRMSE
NyRRRR
0.6
Ny 1000)
NysKRR
NysKRR(M=1000
0.4
0.4
0.2
0.2
0.0
Forecast time
Forecast time (s)
Figure 5: Forecasting error (nRMSE) versus fore- cast time. KRR has similar performance RRR. We used rank = 50 for scalable KAF, NysPCR, NysRRR and no rank for NysKRR. We used m = 250 for the Nyström estimators. For Nys-
Figure 6: NysKRR is very unstable. Here
r = 50 for both NysRRR experiments, 1000 for NysKRR and NysRRR(m=1000) Nys- RRR(m=250) is shown for comparison with Fig- ure 0.0001.
RRR and NysKRR we set 0.0001.
1.2
1.0
sKAF
ARMSE
MysPCR(M
NysRRR(M=1100)
NysKRR(MM-000))
0.2
0.0
time
Figure 7: We increase the regularizer for NysKRR À 0.005. Instability reduced, but aceu- racy is also noticeably worse; in fact it is nearly equal that of NysPCR (with fewer centers) and sKAF, Parameters for NysRRR were 50,
1000, À 0.0001; for NysPCR: = 50, m 250.
Connected to the instability problems, KRR is not suited for Koopman operator learning since without rank regularization, a large number of spurious eigenpairs pop up during estimation.
35


Reviews:
[{"review_id": "anLnhBWZrl", "reviewer": "Reviewer_FXhE", "summary": "The paper talks about estimating Koopman operators. Sketching (random projections), specifically the Nyström method is used as a computational tool to make things tractable. Furthermore, it is shown that these approximates have the same convergence rate as their slower counter parts. Further error bounds are provided in terms of the operator norm.", "questions": "1) What is $A^\\ast$ is line 104? The equation 5 just above has only $A$, are they the same?\n2) In experimental evaluation, only PCR and RRR are considered, what about KRR ?", "limitations": "I don't see a direct potential negative social impact.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The authors do a good job in positing the current work within the literature which is very wide otherwise.\n\nThe flow of the paper is very good and reads well. \n\nAssumptions in Section 4 are clearly mentioned before stating the main result.\n\nI appreciate authors commenting on the non iid case, apart from detailed analysis of the iid scenario.\n", "weaknesses": "I view the papers contributions as two fold. \n\nThe first part is the use of sketching for approximation of the three scenarios of KRR, PCR, RRR considered. In this case (as authors have acknowledged, references [3] and [6] therein), there are  other works who have pointed this out. \n\nAgain authors mention the special case of their sketched PCR as kernel DMD. These make the contributions of the current work from this point of view a bit restricted. Having said that, I greatly appreciate authors mentioning them in the paper!\n\nThe second is in terms of learning bounds via operator norm, where I think the contribution is solid. \n\nHowever, it is hard to have empirical evidence to support this and is understandable. It would be good to have some discussion around this, as Koopman operators are heavily used in practical settings and having some evidence would greatly improve the impact of the paper over it's current form.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their helpful comments, which we address below.\n\n**Novelty and comparison with refs. [3], [6] and kernel DMD.**\n\nReferences [3] and [6] deal with the Nystrom KRR model, with [3] being more concerned with theoretical considerations (in a slightly different setting from us: [3] bounds convergence rates in Hilbert-Schmidt norm, while our work is in operator norm) and [6] with the physical interpretation of the resulting model, without providing theoretical guarantees. \nOur contributions are a generalization of these works to two different estimators: PCR and RRR. We argue that these estimators are particularly useful for Koopman operator learning where the desired rank $r$ can be much smaller than the number of centers $m$ (see for example our last experiment), and it is only necessary to estimate the first $r$ components of the Koopman operator's spectrum. Using rank-reduced estimators reduces the risk of spectrum pollution and (at least in the case of $r \\ll m$) reduces the computational complexity).\n\nIn fact, kernel PCR is equivalent to the kernel DMD algorithm. Compared to kernel DMD, our contribution is the introduction of sketching, which greatly improves the computational complexity making long trajectories (or more generally a large total number of time-steps) easy to handle.\n\n**Empirical evidence to support the learning bounds.**\n\nWith regards to a qualitative demonstration of the learning bounds, we can use the l63 system as an example with relatively small amounts of noise. In Fig. 4 of the rebuttal pdf (and correspondingly in the paper) we show  that -- as the number of training samples increases the test error decreases as $\\sim n^{-1/2}$. While we increase $n$, we must also increase $m$ to satisfy the assumptions (in particular, we set $m = 5\\sqrt{n}$). \nWe perform a qualitative fit of the test error to set two parameters $a$ and $b$ which contain the unknown constants in our theorems ($\\mathcal{E}(\\hat{A}) \\sim \\frac{a}{\\sqrt{n}} + b$), to show that the error approximately follows the expected rate.\nUnfortunately this experiment is too coarse to distinguish between the rates of PCR and RRR, which are also influenced by hyperparameters such as $\\lambda$, and constants depending on the particular dataset. Note that KRR becomes unstable as $n$ (and correspondingly $m$) increases, since it needs to estimate more and more eigenvalues.\n\n**What is A\\* line 104?**\n\nYes, they are the same operator. $A^*$ denotes the adjoint of $A$. Due to the way we formulate the problem in eq. 5 $A$ is the estimator of the adjoint of the Koopman operator, or equivalently $A^*$ is the estimator of the Koopman operator itself. If the passage is a bit unclear we could change the text to say:\n    \n    The operator $A^*$ should thus be understood as an estimator of the adjoint of the Koopman operator ...\n\n\n**In experimental evaluation, only PCR and RRR are considered, what about KRR ?**\n\nWe will add KRR to the Lorenz '63 experiment. Please see the \"Author rebuttal\" (and the accompanying figures in the PDF) for how this would look like, and some motivation for why this was not done in the first place.\n"}, {"review_id": "2qBrueQuG3", "reviewer": "Reviewer_Hh1S", "summary": "This paper improves the computational efficiency of existing kernel method for Koopman operator approximation, including  KRR, PCR, RRR based on the Nystrom estimator. In this way, the authors reduce the original computational complexity from $\\mathcal{O}(n^3)$ to $\\mathcal{O}(n\\sqrt{n})$. They provide theoretical learning bounds and several numerical studies to demonstrate the efficacy of the proposed method.", "questions": "[1] In the experiment of Molecular dynamics datasets, do the $\\phi,\\psi$ correspond to the first two eigenfunctions with largest eigenvalues of the approximated Koopman matrix?  And why not demonstrate the predictive performance of the proposed methods in the original 45-D state space in this task?", "limitations": "The authors have not provided experiments to demonstrate the efficacy of the proposed methods with the non-i.i.d. data, but the current work is complete anyway.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "[1] The theoretical estimation of the learning bounds guarantee the performance of the proposed methods in practical scenarios.\n\n[2] The experiments validate the efficacy of the proposed methods can scale to large dynamical systems.", "weaknesses": "[1] The idea of introducing random projections to the existing kernel methods is not novel enough, since the similar idea has been used in the  deep learning community, such as generative models [R1] and neural ODEs [R2].\n\n[2] The last two experiments seem to be interesting but the authors do not describe them clearly.\n\n[3] In line 62, the authors have typos about 'decomposition'.\n\n__References__\n\n[R1] Song, Y., Garg, S., Shi, J., & Ermon, S. (2020, August). Sliced score matching: A scalable approach to density and score estimation. In Uncertainty in Artificial Intelligence (pp. 574-584). PMLR.\n\n[R2] Grathwohl, W., Chen, R. T., Bettencourt, J., Sutskever, I., & Duvenaud, D. (2018). Ffjord: Free-form continuous dynamics for scalable reversible generative models. arXiv preprint arXiv:1810.01367.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would first of all like to thank the reviewer for taking the time to read this paper and helping us improve its current state. We will address the main points raised below.\n\n**Novelty of random projections & comparison to random projections for generative models and neural ODEs.**\n\nRandom projections are indeed an old idea, since they provide a computationally efficient way of performing dimensionality reduction. They have been used in many different areas of computer science, including more recently deep learning. In kernel methods, random projections often go under the name of the Nyström method when they can be interpreted as subsampling operations on the kernel matrix.\nHowever, there are far too many instances of random projections in machine learning to mention them all, or even to give a partial overview. For this reason we only cite the ones which are more relevant to the models being considered, similarly to also for example [R1] which does not cite the Nyström method in kernel ridge regression.\n\n**Clarity of last two experiments**\n\nFor a description of the experiment on the alanine dipeptide please see the next point. Regarding the last experiment, we will try to explain it better here. Please let us know if some aspects remain unclear.\nAfter learning the Koopman operator with the Nystrom RRR model on a training set, we compute the first three eigenfunctions on the test dataset. \nWe then train a clustering model on top of the eigenfunctions (PCCA+ model) with three states. \nIn the plot, for each eigenfunction-state pair, we have time along the trajectory on the x-axis, the eigenfunction value on the y-axis, and color indicating the probability of that particular state. So eigenfunction 1 clearly distinguishes between states 1 and 2, while eigenfunction 2 clearly distinguishes between states 0 and 2. \nFinally we plot some representative structures for each state showing that indeed we have identified three different conformations that the protein can take throughout the trajectory, distinguishing them in a very high-dimensional input state space.\n\n**Meaning of $\\phi, \\psi$ in the alanine dipeptide experiment. Why not predictive performance on original state space?**\n\nThe alanine dipeptide can be described in two ways: either through a 45d state (pairwise distances between heavy atoms), or through a 2d state (two angles between backbone atoms, called $\\phi$ and $\\psi$).\nIt is well known that the important dynamics of the protein (changes in the conformation which are stable over time) can be characterized fully from the 2d state. In particular there are 2 transitions between the left and right parts of the plot and top and bottom parts. \nThrough our experiment we wish to show that even though we learn the Koopman operator on the larger state, the first two eigenfunctions can capture the stable transitions of the protein well. \nTo show this we can overlay the eigenfunction values on top of the 2d-state plot, and verify that the eigenfunction changes are highly correlated with changes in the angles $\\phi$ and $\\psi$.\nOn this task, predictive performance (i.e. forcasting $x_{t+1}$ from $x_t$) is not a common metric to use. In fact the main reason for learning a dynamical system through the Koopman operator is to obtain the spectrum. If one were not interested in that, and just in forecasting, there are lots of other ways to tackle the problem. This is the reason behind our focus on the eigenfunctions here.\n"}, {"review_id": "Gm4Kz7EiBp", "reviewer": "Reviewer_jbb3", "summary": "This paper proposes sketching algorithms for approximating three different kernel-based estimators of Koopman operators from large scale temporal data. New non-asymptotic error bounds are derived. ", "questions": "- What is the complexity with respect to the dimension of the state space? Is it linear? \n- Does the Markovian assumption in line 84 present a limitation for real-world applications? Is there a way to adapt Koopman operator learning (not necessarily sketched) to the non-Markovian case? ", "limitations": "Limitations are discussed. ", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- The paper leverages the concept of sketching to speed up the estimation of Koopman operators without sacrificing accuracy. The newly proposed algorithms can potentially have a high impact in several scientific disciplines where the analysis of dynamical systems is of interest.  \n- The authors explain and contextualize very well their contributions and the related work is comprehensive. \n- The paper is very well written and organized. \n- New theoretical results are rigorously derived to provide robust guarantees for the proposed methodology.\n- The empirical analysis is sound and well-conducted. To showcase various aspects of the method, two distinct applications are considered. The first experiment effectively demonstrates the tradeoff between time and accuracy. The second experiment is particularly interesting as it provides a practical example where it is necessary to consider full, non-subsampled long trajectories. Furthermore, it focuses on the interpretation of dynamics and scientific discovery. It also highlights the scalability with respect to the state space dimension, which serves as motivation for the study of Koopman operator methods.", "weaknesses": "- While the idea of using sketching is not original, this contribution has the potential to significantly advance dynamical systems analysis by providing efficient learning algorithms for Koopman operators. Sketches for operators mapping between Hilbert spaces have recently been studied in [3]. However, the authors outline the unique aspects of their contribution and provide a clear motivation for deriving new results.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the anonymous reviewer for their evaluation of this paper, and for their kind comments. We will address the main questions below.\n\n**Complexity with respect to the dimension of the state space?**\n\nYes, essentially the complexity with respect to state space dimensionality is contained in the kernel function: computing an entry of the kernel typically depends linearly on the dimension of the states $d$, e.g. think about the Gaussian kernel $k(x, x') = exp(-\\gamma \\lVert x - x' \\rVert^2)$.\n\n**Limitations of the Markovian assumption, and adaptation to the non-Markovian case**\n\nThe Markovian assumption is essentially always used in molecular dynamics, where often even more restrictive assumptions are considered (such as time-homogeneity and time-reversibility). The definition of the Koopman operator could be generalized beyond this setting, for instance by considering the dynamics of the features of an history of k>1 consecutive states, and one could relatively easily extend the proposed algorithm. However the analysis in this setting would be radically different and is out of the scope of this paper; for instance the spectrum of the Koopman operator has been well studied in the Markovian setting and is known to be of key interest, while it is unclear how to generalize and leverage such a spectral analysis in a non-Markovian setting.\n"}, {"review_id": "isXNrSn9eh", "reviewer": "Reviewer_8Gvs", "summary": "This paper utilizes theory developed around random projections to develop fast/efficient numerical methods for approximating the Koopman mode decomposition of data sets with long-time trajectories, that additionally come with learning rates that are optimal. ", "questions": "To summarize, the Weakness section above, \n\n1. Why would one want to use very long-time trajectory data when computing the Koopman mode decomposition? (You explain this some, but making it more explicit would be helpful)\n\n2. What are the differences between KRR, PCR, and RRR in their performance on the different data sets studied (at least the Lorenz '63)?\n\n3. How are the random $\\tilde{\\Phi}$ are generated? ", "limitations": "The authors do a sufficient job explaining their limitations (the assumptions they make) and what kinds of data sets their proposed methods work best on. The only exception to this is that it should be clarified that the proposed methods improve Koopman mode decomposition estimates of long-time trajectory data and not large dimension data. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "1. This paper developed new fast/efficient algorithms for computing the Koopman spectra, the effectiveness of which was demonstrated on very large data sets with long-time trajectories. These include a traditional benchmark (Lorenz '63), as well as molecular data sets, which are a growing area of interest to the machine learning and Koopman communities. Additionally, the molecular data sets were of a size that could not have been studied previously. \n\n2. This paper provided learning rates for their new algorithms, which they show to be optimal. \n\n3. This paper discusses how different methods for computing the Koopman mode decomposition (e.g., Nystrom RRR vs. PCR) differ in their learning rate dependencies. \n\n4. The paper was well-written. \n", "weaknesses": "1. Motivation for why very long-time trajectories of data would be needed for some systems was missing from the Introduction. This was partially discussed at the end in Sec. 5 \"Molecular dynamics datasets\", but developing this more is important. Additionally, the title and some places in the text describe the sketching approach as being useful for \"large scale dynamical systems\". As I understand it, the methods are beneficial for \"long-time trajectory dynamical systems\". Of course, the molecular examples shown have both large dimension and long-time trajectories, but the title and text should be corrected to accurately emphasize that it is long-time trajectory dynamical systems that this method is useful for. \n\n2. The discussion on how Nystrom PCR and RRR differ in their learning rates was very interesting. Including more on this (maybe in table form) would be nice. Additionally, it would be helpful to see all 3 methods compared, so KRR should be included in Figure 1. Additionally, the reason for using just RRR in Figs. 3 and 4 should be noted.\n\n3. I found it a little confusing that going from Eq. 8 to Eq. 9, the P_Y \\hat{C}_{YX} P_X  became K^\\dagger_{\\tilde{Y}, \\tilde{Y}} K_{\\tilde{Y}, Y} K_{X, \\tilde{X}}, but that going from Eq 10 to Eq 11 the P_Y \\hat{C}_{YX}  also led to K^\\dagger_{\\tilde{Y}, \\tilde{Y}} K_{\\tilde{Y}, Y} K_{X, \\tilde{X}}. I assume this has to do with the way the \\dagger was absorbed into the [ \\cdot ]_r?\n\n4. How the random $\\tilde{\\Phi}$ are generated should be explained in the text. \n\nMinor comments:\n\n1. Mezic, 2005 should be cited when discussing the Koopman mode decomposition and Budisic et al. 2012 should be cited when discussing how Koopman mode decomposition has been leveraged in the past (line 95). \n\n2. Abbreviations should be explained before using them (DMD, tICA, VAC, etc.). \n\n3. The box plot denoting the sKAF data in Fig. 2 has a strange shape in its lower box (an inverted U). Why is this? \n\n4. It would be helpful and appropriate to have citations when noting that \"An important application of Koopman operator theory is in the analysis of molecular dynamics\". \n\n5. The colors of Fig. 3 should be explained in more detail. \n\n6. It might be helpful to have a line across the 3 states in Fig. 4 to show how Eigenfunction 1 has linear separation between State 0 and the others, and Eigenfunction 2 has linear separation between State 1 and the others. \n\n7. There are a few typos: \ni. \"allows to deploy\" (line 1)\nii. \"as theirs exact\" (line 68)\niii. \"to study large class\" (line 90) \niv. \"]0, 1]\" (Assumption 4.3)\nv. \"]0, \\tau]\" (Assumption 4.4)\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for the helpful comments. We address the main ones below, while minor comments and typos will directly be fixed in the revised version.\n\n**Motivation for long trajectories**\n\nWe stress that our contribution allows to tackle large datasets, which can in practice be built using multiple trajectories in order to simply improve the overall accuracy. However, there are some applications (see molecular dynamics paragraph) where it is indeed physically necessary to record very long trajectories to capture events occurring at large scales, but that could be missed if simply subsampling the dataset. We will add a sentence to clarify this in the introduction.\n\n**Large scale vs. Long time**\n\nThere is indeed a difference between the dynamical systems nomenclature, where large-scale is understood in terms of dimension of the data, and the regression setting where it usually refers to the number of data points. \nAlthough our contribution only focus on improving the complexity wrt the number of points, kernel methods have a linear complexity with respect to the dimension, which makes of our method an good candidate to deal with both large $n$ and $d$. We will however clarify the fact that sketching only reduces the complexity w.r.t. $n$ in order to avoid confusion.\n\n**Comparison between KRR, PCR and RRR**\n\nWe added as requested KRR to the experimental results for the Lorenz' 63 dataset (see PDF). \nExplaining the behavior of Nyström KRR can be complex since it does not have a notion of rank, and will always output the full-rank estimator.\nWe repeat the same experiment as in the paper, with M=250, r=50 adding the KRR estimator. In this case KRR is on par with RRR: no instability is encountered at this point, and KRR is about 2x faster than RRR. (figure 1)\nIf we increase the number of centers, without increasing regularization, KRR becomes highly unstable while RRR maintains the rank regularization and hence has no problem, obtaining a small performance boost thanks to the increased M (figure 2). \nTo make KRR more stable, we can increase regularization (from λ=1e-4 to λ=5e-3) and the instability is not completely fixed, while the accuracy is noticeably reduced, to the point where it's nearly the same as PCR with fewer centers (figure 3).\n\nAnother reason why KRR is not a great estimator for Koopman operator learning, is that having the all $m$ components in the spectrum leads to many spurious eigenpairs popping up.\nWe will condense this rather long discussion for the camera ready paper.\n\nRegarding theoretical rates of PCR and RRR, we indeed provide explicit constants in each case and discuss the role of the rank, however both estimators overall enjoy the same (optimal) rate. For this reason, it seems that adding a comparative table to the paper would not bring much. We will however stress this point in the introduction.\n\n**Estimator derivation**\n\nYes, unfortunately the process of going from the infinite-dimensional formulation to the kernel formulation is a bit tricky. Part of the question can be answered considering $P_Y = \\tilde{Φ_Y} K_{\\tilde{Y}, \\tilde{Y}}^† (\\tilde{Φ_Y})^*$ (and similarly for X). The computations to get from the first to the second formulation involve moving part of what's in the inverse outside of it for KRR. For PCR we first simplify the $[ · ]_r^†$ expression where some of the inverses cancel out, and a similar expression outside of $[·]_r$ appears when plugging the two pieces together. The full derivation is in Appendix C.\n\n**Generation of $\\tilde{Φ}_X$**\n\nThe Nyström feature map comes from evaluating the normal feature map $\\phi: \\mathcal{X} \\to \\mathcal{H}$ on the so-called \"Nystrom centers\", which are randomly selected. We mention line 139 that we sample them uniformly, while many other strategies have been studied (e.g. leverage scores sampling, k-means).\nSomething that is not mentioned (and will be added) is that the centers for X and Y are sampled in the same way (i.e. if $x_t$ is chosen as Nyström center in $\\tilde{Φ_X}$ then $y_{t+1}$ is chosen as Nyström center in $\\tilde{Φ_Y}$). \nWe also tried empirically to take $P_X=P_Y$, which yielded a similar empirical performance.\n\n**Citations and other related work**\n\nThanks for pointing this out, we realize that the \"related work\" section is missing a lot of references and we will add the two suggested references for the camera ready version.\nSimilarly for the usage of Koopman operators in molecular dynamics we will reference as examples [1], [2], [3] and as a review [4]\n\n[1] \"A Small Molecule Stabilizes the Disordered Native State of the Alzheimer’s Aβ Peptide\", Löhr et al. (2022).\n\n[2] \"Markov State Models and tICA Reveal a Nonnative Folding Nucleus in Simulations of NuG2\", Schwantes et al. (2016).\n\n[3] \"Slow dynamics of a protein backbone in molecular dynamics simulation revealed by time-structure based independent component analysis\", Naritomi et al. (2011).\n\n[4] \"Machine learning for protein folding and dynamics\", Noé et al. (2020).\n\n**The colors of Fig. 3 should be explained in more detail.**\n\nWe will expand the figure caption to explain that the color of the two left panels corresponds to the value of the first two eigenfunctions, evaluated at different protein conformations. \nThe eigenfunctions are evaluated on the 45-d space, but the plot is overlayed on top of a 2-d subspace which is known to represent the protein dynamics well. \nThe third panel shows discrete colors after clustering of the dynamics projected onto the first eigenfunctions (using PCCA+ with 3 states), showing that the obtained clusters correspond to the main areas of the Ramachandran plot.\n\n**Helpful to have a line across the 3 states in Fig. 4 to show linear separation**\n\nThanks, we have updated the figure and the caption correspondingly (see attached PDF).\n\n**Artifact in Fig. 2**\n\nThis seems to be an artifact from the matplotlib boxplot, there just is not much variance there."}]

Author rebuttal:
We thank all the anonymous reviewers for their thoughtful comments.
We provide here a pdf file with additional experimental results, and address directly the comments of the reviewers in the individual rebuttals.

In particular, the experiments concern the addition of Nystrom KRR to a comparison with the other two (PCR and RRR) estimators, which was aked for by two reviewers. We repeat here the comments made on the rebuttal for Reviewer 8Gvs since they may be of interest to the other Reviewers.

**Comparison between KRR, PCR and RRR**

We added as requested KRR to the experimental results for the Lorenz' 63 dataset (see PDF). Explaining the behavior of Nyström KRR can be complex since it does not have a notion of rank, and will always output the full-rank estimator. We repeat the same experiment as in the paper, with M=250, r=50 adding the KRR estimator. In this case KRR is on par with RRR: no instability is encountered at this point, and KRR is about 2x faster than RRR. (figure 1) If we increase the number of centers, without increasing regularization, KRR becomes highly unstable while RRR maintains the rank regularization and hence has no problem, obtaining a small performance boost thanks to the increased M (figure 2). To make KRR more stable, we can increase regularization (from λ=1e-4 to λ=5e-3) and the instability is not completely fixed, while the accuracy is noticeably reduced, to the point where it's nearly the same as PCR with fewer centers (figure 3).

Another reason why KRR is not a great estimator for Koopman operator learning, is that having the all components in the spectrum leads to many spurious eigenpairs popping up. We will condense this rather long discussion for the camera ready paper.

Meta-review: The paper presents novel "sketched" estimators for learning Koopman operators.. The authors place their work within the wider literature effectively, which was well received by the reviewers. Kernel-based Koopman operator estimators are commonly used to predict and analyze complex dynamical systems, but their application to lengthy trajectories leads to computational challenges. This paper proposes to use random projections (sketching) to address this issue.

Beyond its theoretical contribution, the paper demonstrates practical utility through a series of well executed simulations. The empirical evaluation involved diverse settings, highlighting the trade-off between time and accuracy, and included a case where non-subsampled long trajectories were necessary.

The proposed methodology is convincing - a combination that could be impactful across fields where the analysis of dynamical systems is central. Given its potential impact in the field of large-scale dynamical systems, this paper is unanimously recommended for acceptance.
Decision: Accept (poster)