Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
Authors: Vladimir R Kostic
Karim Lounici
Hélène Halconruy
Timothée Devergne
Pietro Novelli
Massimiliano Pontil
Venue: ICML 2025 poster (2025)
Keywords: continuous Markov processes
statistical learning theory
data-driven dynamical systems
spectral decomposition
kernel methods
Primary area: theory->learning_theory
TL;DR: We develop finite sample learning bounds for infering the spectral decomposition of a strongly continuous semigroup of Markov operators from a single trajectory of non-iid data.

Abstract: Markov processes serve as universal models for many real-world random processes. This paper presents a data-driven approach to learning these models through the spectral decomposition of the infinitesimal generator (IG) of the Markov semigroup. Its unbounded nature  complicates traditional methods such as vector-valued regression and Hilbert-Schmidt operator analysis. Existing techniques, including physics-informed kernel regression, are computationally expensive and limited in scope, with no recovery guarantees for transfer operator methods when the time-lag is small. We propose a novel method leveraging the IG's resolvent, characterized by the Laplace transform of transfer operators. This approach is robust to time-lag variations, ensuring accurate eigenvalue learning even for small time-lags. Our statistical analysis applies to a broader class of Markov processes than current methods while reducing computational complexity from quadratic to linear in the state dimension. Finally, we demonstrate our theoretical findings in several experiments.

Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
Vladimir R. Kostic 12 Karim Lounici 3 Hélène Halconruy *45
Timothée Devergne'
Pietro Novelli' 1 Massimiliano Pontil 16
Abstract
The central mathematical object for describing Markov semi- groups is the Infinitesimal Generator (IG), which governs the evolution of probability distributions over the state space. Its spectra reveal important features such as metastable states, transition statistics, and committor functions, all of which critical for understanding system dynamics. Ac- curately learning the spectral decomposition of the IG is thus pivotal for wide range of applications, including molecular dynamics, time-series clustering, computational neuroscience, and beyond.
Markov processes serve universal models for many real-world random processes. This paper presents a data-driven approach to learning these models through the spectral decomposition of the infinitesimal generator (IG) of the Markov semi- group. Its unbounded nature complicates tradi- tional methods such vector-valued regression and
operator analysis. Existing techniques, including physics-informed kernel re- gression, are computationally expensive and lim- ited in scope, with no recovery guarantees for transfer operator methods when the time-lag is small. We propose novel method leveraging the IG's resolvent, characterized by the Laplace trans- form of transfer operators. This approach is robust to time-lag variations, ensuring accurate eigen- value learning even for small time-lags. Our statis- tical analysis applies to a broader class of Markov processes than current methods while reducing computational complexity from quadratic linear in the state dimension. Finally, we demonstrate our theoretical findings in several experiments.
The field of molecular dynamics has particularly benefited from spectral decomposition methods of Markov semi- groups. Research on Al-augmented molecular dynamics, grounded in statistical mechanics, highlights the importance of accurately identifying spectral gaps (the separation be- tween slow and fast modes of dynamics) in molecular sim- ulations, see (Schütte et al., 2001). Recently, theoretical advancements (Kostic et al., 2024a) were used in (Dev- ergne et al., 2024) to demonstrate the effectiveness of IG- based methods in accelerating simulations and enabling the practical identification of metastable states. The authors phasize that IG methods overcome the limitations of more widely used transfer operator (TO) approaches when ex- tracting dynamical information from biased data, and they underline the scalability to larger proteins as a particularly important advantage of IG methods.
1. Introduction
Markov semigroups play a critical role in modeling dynam- ies of complex systems across various fields, including op- tion pricing in finance (Karatzas & Shreve, 1991), molecular dynamics (Schütte & Huisinga, 2003) and climate modeling (Majda et al., 2009), where understanding long-term behay- ior is essential for accurate forecasting and interpretation.
On the other hand, Klus & Djurdjevac Conrad (2023) intro- duced Transfer Operator (TO)-based spectral clustering method tailored for directed and time-evolving graphs. By leveraging TOs, their approach demonstrates how to iden- tify coherent sets within complex networks, enhancing the analysis of temporal data structures. Furthermore, Cabannes & Bach (2024) emphasize that IG methods open exciting new directions for spectral-based algorithms, be it spectral clustering, spectral embeddings, spectral distances.
'Equal contribution 'CSML, Istituto Italiano di Tecnologia, Genova, Italy FFaculty of Science, University of Novi Sad, Ser- bia 'CMAP, École Polytechnique, Paris, France *SAMOVAR, Télécom Sud-Paris, Institut Polytechnique de Paris, France "MODAL'X, Université Paris Nanterre, France "AI Center, Univer- sity College London, UK. Correspondence to: Vladimir Kostic coladimirkosstic@@ itt.tt
In neuroscience, Marrouch et al. (2020) used TO (also known as the Koopman operator) to analyze brain activ- ity. Their work shows how the operator's spectrum captures the spatiotemporal dynamics of neural signals, providing insights into brain function and possible applications diag- nosing neurological disorders. Ostrow al. (2023) further
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 20255 by the author(s).

---

Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
developed a dynamical similarity analysis based on TOs to distinguish learning rules in unsupervised way, showing that the TO spectrum supports comparative analysis of the temporal structure of computation in neural circuits.
& Bach, 2023; Hou et al., 2023) either have limited scope or offer only partial or suboptimal analysis, as summarized in Table 1. Crucially, none adequately addresses the challenge posed by the unboundedness of the IG, leading to incom- plete frameworks and suboptimal convergence rates, which in some cases depend explicitly on the state space dimension. Moreover, the estimators in these works are susceptible to spurious eigenvalues and do not offer guarantees for ac- curate estimation of eigenvalues and eigenfunctions. The current state-of-the-art (Kostic et al., 2024a) introduces a physics-informed kernel regression method only for Markov processes admitting a Dirichlet form. This approach lever- ages the Dirichlet form to define an energy- -based metric for learning the model and provides a comprehensive statistical analysis with learning guarantees for the spectral decompo- sition of the IG while avoiding spurious eigenvalues. That said, their analysis applies only to self-adjoint IGs, assumes partial knowledge and i.i.d. data, and requires gradients of the feature map, leading to a quadratic scaling with the state space dimension d, thereby hindering broader applicability.
However, the typically unbounded nature of the IG makes designing efficient and reliable estimators challenging. In this paper we address this problem through the lens of Laplace transform. As we shall see, this fully data-driven approach provides a way to bypass both the difficulties of es- timating an unbounded operator and overcomes limitations of the transfer operator based algorithms when trajectory data is sampled at high frequency. In particular, we show theoretically and empirically that our method enables accu- rate and robust estimation of eigenvalues and eigenfunctions even for arbitrarily small time-lags.
Related work. A substantial body of research has focused on using transfer operators to learn dynamical systems from data (see the monographs by Brunton et al., 2022; Kutz et al., 2016, and references therein). This has led to the develop- ment of two primary approaches: deep learning methods (Bevanda et al., 2021; Fan et al., 2021; Lusch et al., 2018), which excel in capturing complex data representations but often lack rigorous statistical foundations, and kernel meth- ods (Das& & Giannakis, 2020; Klus et al., 2019; Kostic et al., 2022; 2023; Williams et al., 2015), which offer strong sta- tistical guarantees for Transfer Operator (TO) estimation but require kernel functtinn selection. A closely related chal- lenge, learning invariant subspaces of the TO, has led to several methodologies (see e.g. Li et al., 2017: Mardt et al., 2018; Tian & Wu, 2021, and references therein), some of which leverage deep canonical correlation analysis (Andrew et al., 2013; Kostic et al., 2024c). Note that TOs share the same eigenfunctions as the IG, which motivates the develop- ment of TO methods aimed at learning the spectral proper- ties of the IG. Yet TOs are sensitive to time-lag choice, with their spectral gap deteriorating as the time-lag decreases, making existing spectral recovery guarantees ineffective for small lags-an an issue observed in practice, see, e.g., (Bonati et al., 2021). This bottleneck is especially problematic in complex tasks like enhanced sampling (Laio & Parrinello, 2002; Shmilovich & Ferguson, 2023). To address this, re- search has focused on learning the IG and its eigenstructure directly. As recently shown (Devergne et al., 2024), IG learning can be combined with enhanced sampling methods to efficiently debias data and reveal true dynamics. Still, research on IG learning remains limited and often ad hoc. For instance, (Zhang et al., 2022) developed deep learning method for Langevin diffusion, while (Klus et al., 2020) ex- tended dynamic mode decomposition to Iearn the generator, connecting it to Galerkin's 's approximation. However, neither of these works provides any formal learning guarantees. To the best of our knowledge. most existing works with theo- retical guarantees (Cabannes & Bach, 2024; Pillaud-Vivier
Contributions. We introduce novel approach for ac- curately estimating the spectral decomposition of the IG for a broad class of Markov semigroups, encompassing all models considered in prior work. Our method leverages useful connection between the resolvent of the IG and the semigroup of TOs via the Laplace transform. Unlike TO methods, it estimates the IG directly, avoiding small time-lag issues and preserving a larger spectral gap for more accurate eigenvalue and eigenfunction learning. We provide sharp statistical guarantees, valid for data sampled from a trajectory in the stationary regime, accounting for slow mix- ing effects. Our results are the first to apply to wide class of Markov semigroups with sectorial IGs. A key technical contribution is our bound on the Bochner integral approx- imation error using Crouzeix's bound, which may be of broader interest. Computationally, our method combines multiple TOs at differrent time-lags through a single matrix product between Toeplitz matrix and the kernel embed- ding, reducing complexity to O(n²d), d)) making it practical for high high-dimensional systems. The complexity can be fur- ther reduced while preserving accuracy by utilizing standard scaling techniques such as random Fourier features. Our ex- periments show a striking performance improvement over bO-based and other IG methods, as predicted by our theory.
2. Background
Various physical, biological, and financial systems evolve through stochastic processes X = (X,) where X, e X' C Rd represents the system's state at time t. We focus on continuous-ting Markov processes with continuous paths, which essential for modeling these systems. This class includes Itô diffusion processes (see Ex. 2.1 and 2.2), re-
2

---

Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
flected time-changed Brownian motions, and processes with local time. Markov processes model phenomena where the future depends only on the present, not the past, and are described by their wwss-meaurres on the path space. This foundational approach to Markov theory defines the process through the infinitesimal generator (IG), a key linear (often unbounded) operator on space of observables (functions defined on the state space).
with eigenvalues (Ax)ieN C C and corresponding left and right eigenfunctions I 9i E T,?,, respectively.
Resolvent operator. Eigenvalues, while informative about long-term behavior, fail to capture transient dynamics of the full time evolution of the process whenever IG is non-normal, that is when LL* + L' (Trefethen & Em- bree, 2020). In contrast, the resolvent of L defined by R2 := (ul / p(L), provides a more comprehensive view of the dynamics, making it the core object of spectral theory of IG. Through its characterization via the Laplace transform, (see for instance (Bakry et al., 2014), equation (A.1.3) as
i p(L) being a shift parameter
Markov theory. The dynamics of continuous-time Markov process X is described by family of probabil- ity densities (pt)CER+
P(X, E B|X =) = SEPeC u)dy,
(1)
R So Aie dtd
(5)
and transfer operators (TO) (A)),ER, such that for all EE
R4, E' E B(X), E X' and measurable function f: X - +R,
it is intrinsically connected to the TO defined for the time- lag / by A, = e¹4, Moreover, the resolvent encodes both the spectrum of L (eigenvalues via its poles and the continuous spectrum) and transient phenomena, such as the system's 's approach to equilibrium. Last, the resolvent is essential for analyzing stability under perturbations and understanding how changes in the IG affect the dynamics.
= =xf=If(()pcc E(())))))))))_E( :.. (2)
= [[f(x)
The transfer operator is essential for understanding the dy- namies of X. We study its action on 2 (X'), the space of functions on that are square-integrable with respect to an in- variant measure x. which satisfies 7 = x for all i E R4. We assume that the Markov process X meets two key proper- ties regarding π: [1] π ensures long-term stability, meaning X converges to π from any initial state x. [2] The process exhibits geometric ergodicity, meaning converges expo- nentially fast to the invariant measure. Finally, the process is characterized by the infinitesimal generator L, defined for I E (A)) by the limit I.f= lim, (A,s-f)/ with I.
Link with SDEs. Itô diffusion processes are a key example of Markov processes, governed by stochastic differential equations (SDEs) of the form:
dX, = a(Xx) +
Xo = 2.
(6)
where x E X, W = (W!.... W? DEER is a standard p-dimensional Brownian motion, the drift a : X + Rd and diffusion b X' - ddxx are globally Lipsehitz and sub-linear. This ensures a unique solution X (Xt)c20 in (X, ((x)). SDEs like (6) include Langevin dynamics processes. The IG I. associated with (6) is second-order differential operator, defined for E L2 (X) and x E X', as:
being closed on its domain.
The class of sectorial generators consists of the operators generating strongly continuous semigroups, analytic in a sector of the complex plane defined by growth conditions in an angular region, i.e., L is a (stable) sectorial operator with angle 0 E 2/2),
and
F(L) C C = (= C R(=) 5
V[())=V[(()T((()) + (V² (7)
[3(1] < R2) tan( (3)
where (2f= f) is the Hessian of f. Its domain is the Sobolev space W22(x)=15 (x) | ||/l + \V/l <oo]. Its spectral de- composition allows one solve SDE (6), that is
where F(L) denotes the numerical range of L. This class covers all time-reversal processes (self-adjoint IG), but also important non-time-reverss processes, such as advection-Diffusion and underdamped Langevin (Kloeden et al., 1992).
E[f(X,) [x=x|= ext (g: fzz())).
(8)
Example 2.1 (Overdamped Langevin). Let o, hn and T E Rt. The overdamped Langevin dynamics of partiele in a potential V :R + satisties (6) with a VVV and b 22(kTT/T). where 7, kb, and T are the friction coefficient, Boltzmann constant, and system temperature, respectively. The invariant measure is the Boltzmann dis- tribution π (dx) x e (Va)/((RTT) dx. In dissipative systems, the IG L is sectorial, with its spectrum usually in the left half-plane
Spectral decomposition. When continuous for some y C, the operator R - -1 is the resolvent of L., and (L))= L is bijective, Ry is continuous] } is called the resolvent set. For sectorial operator, the resol- vent is uniformly bounded outside a sector containing the spectrum. The spectral decomposition of the IG can be
written as
I, = EiEN gi 8
(4)
3

---

Laplace Transform Based Low-Complcxity Learning of Continuous Markoy Semigroups
Aspect
(Cabannes & Bach, 2024)
(Hou al., 2023)
(Kostic et al., 2024a)
Our work
Many Markov Processes
x (only Laplacian)
x (only Ito)
x (only Dirichlet form)
Risk metric
C norm
C norm
Weighted Sobolev norm
C norm
Required prior knowledge
((fll info. needed)
(partial info. needed)
Data
iid
trajectory
iid
trajectory
IG error bound
On 2255))
Var 0 x (spuriousness)
On (04))))
O(n )लोي)
Spectral rates
X (spuriousness)
/ (self-adjoint)
(sectorial)
((rVd) n²)
Computational complexity
O Onn2+n
d())d))
O(r)Pd))
Table 1. Comparison to previous kernel-based works on generator learning. State-space dimension is d, N is the number of features (possibly -o))), is sample size and r < max(n, .) is estimator's rank. Our learning bounds derived in Thm. 6.2 where parameters a E [1, 2] and B E (0, 1] quantify the intrinsic difficulty of the problem and the impact of kernel choice on IG learning.
(OU) processes). The OU process with drift is governed by the SDE (6) with a(x) = Ar and b = B, where A E Radxd is the drift ma- trix and B E Rdxd is the diffusion matrix. This models systems like the Vasicek interest rate and neural dynamics, where fluctuations return to equilibrium. If the real parts of A's eigenvalues are negative, the OU process has in- variant Gaussian distribution with covariance 2x satisfying Lyapunov's 's equation: Alxx + A =BBB.
Example 2.2
(MSE) w.r.t. stationary distribution
R(G) = Exo (lx) - G ((x))
(10)
and we can show, c.f. App. A.2, the universal approxi- mation result for its minimizers over bounded operators in H. Namely, since Rex
R(G) ming R(G)
Sy if H is dense in 2 (X) and the injec- IIS tion operator Hilbert-Schmide then one can find arbitrarily good finite-rank approximations of Sy.S However, learn- ing the IG alone is insufficient for forecasting the process, and estimating the spectral decomposition of L is of greater interest. But, as noted in (Kostic et al., 2023), as the es- timator's 's rank increases, metric distortion between H. and C (X) hinders learning.
3. Problem formulation
Let H be an RKHS with kernel k: X' x X + R, and : X 7 H be a feature map such that k(x,x') (o(x), ((x')) for all x,, E X. We assume that H C C (X)), enabling us to approximate L: (x) (X) with an operator G H + H (Kostic et al., 2022). Al- though H is subset of (2)), they have different metric structures, so for f.g E H. f. + (f.9) c2- To resolve this, introduce the injection operator Sx H - (2). which maps each f E H to its pointwise equivalent in C (X') with the appropriate L norm. For bounded kernels this operator is Hilbert-Schmide allowing one to efficiently learn bounded operators on L via finite rank approxima- tions.
4. Approach and main results
The main bottleneck in the risk functional (10) is the inte- gral computation in (9), which hinders standard operator regression methods. In this section, starting from simple idea to approximate (9) with numerical integration schemes, we present a novel fully data-driven method that addresses this difficulty. Namely, consider
While regressing directly the generator might lead to learn- ing spurious spectra due to its unbounded nature, the resol- vent operator is bounded for μ > 0, and, hence, recalling (5), the standard regression risk is well defined via
em (Xo) = ΣS= o mx)((,,)
(11)
where m = (m²) are real weights given by the famous trapezoid rule with e 2 1 points and dine-isscetizzatinn Al O, that is
R(G)= ExXoN 1 hi (tt (xhal(()))))
KEN
Al c [ 2
4 = jj and my=
Atc if1 sj (12)
where (he)kEN is any orthonormal system of H. Next, if we define the target feature via Bochner integral as
So, we estimate Rp by learning
(X,) = 10xh ((X,) e vf dt.
(9)
the risk can be equivalently written as the mean square error
Rm = Ei
(13)

---

Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
in which case the problem of minimizing the risk (10) over
Analyzing (16), when the sampling frequency is 1/ 1/Atzm!/2 and the hyperparameters are chosen such that 1aann-1/6 114 and μ= / the learning rate |444 /2n y=n
G: H + H. transforms to
min R(G)=EX,ma (Xo) ((Xu)l) (14)
for the operator norm error -1/2 Ins/ (n). This matches the learning rates in (Kostic et al., 2023; 2024a) for TO and Ry, respectively.
G: 7
which can be solved by Reduced Rank Regression (RRR) estimator proposed in (Kostic et al., 2022). Namely, to learn (5), we constrain (14) to rank-r-RKHS operators G E
Further, the error/risk analysis led to spectral learning rates, specifically for estimating the eigenvalues and eigenfunc- tions of L. The key difference in the analysis is that the hypothetical domain H typically has a different geometry (norm) than the true domain L (Kostic et al., 2022; 2023). To control the potential deterioration of spectral learning rates relative to the risk/emror, one must analyze the met- ric distortion of the estimator's eigenfunctions, defined as (h) = || / ||A|llgg for h E H. When the metric distor- tion is uniformly bounded (as can occur in well-specified settings), the eigenvalue bound becomes
(?) = {G: H - H| rank(G) < rr. and obtain the
solution
(15) mzZi, is the ag- where Hm = Sy S, EE
mmyy 1/1/21,,1/211
gregated cross-covariance obtained by combining -,, (xx) @ o(X,,)] being the cross-covariance operators in RKHS H w.r.t. invariant measure ग. and - is the regularized covariance, since C- =xx==== @ (X)]-
While the computational details for deriving the empirical version of (15), denoted by Gm are presented in Sec. 5, the main challenge in the statistical analysis of the risk/error bounds, compared to the standard TO case, lies in address- ing both the bias from approximating the integral and the variance from non-iid data collected along single trajectory sampled at frequency 1/A.. While the general case is dis- cussed in Sec. 6, we focus here on well-specified learning problems using any universal bounded kernel, specifically for eigenvalue estimation of self-adjoint operators.
A - (
SAL In 3/2
(n22 (8) PVn
(17)
I A--Alle A1 (R
where 7, = 1/0; are estimates of the generator's eigen- values Xe for E [r). with - =iein 17s hú 60 Ds being the spectral decomposition. Note that (17) can be trans- formed into eigenfunction bounds using standard arguments, see App. C.5 for technical details.
Due to the unbounded nature of the generator, (13) with the choice of (12) may not always provide good approxima- tion of the integral transform (5). However, for a large class of problems with sectorial IG, such as Examples 2.1 and 2.2, we are able to prove, see App. C.2.4, that [R Rm S ctt, where e a constant when EpAit is sufficiently large, and, consequently, obtain that the differ- ence between the true risk and its approximation is bounded in terms of the time-lag parameter. Concerning the vari- ance, the main challenge is accounting for the unavoidable data dependence by aggregating concentration inequalities at multiples of the initial time-lag, leveraging the mixing time from geometric ergodicity. Our approach reveals the impact of key parameters (shift μ > 0, regularization y > O, L eigenvalues, and time-lag At) on the variance.
Investigating the spectral learning rate in (17), we find that tuning ji in an unbounded manner is prohibitive; if yu is too large, the resolvent's 's eigenvalues collapse to zero. Thus, for spectral estimation, we must fix a small μ > O, which influences the rate and the optimal relationship between and At. Specifically, the spectral learning rate becomes n- 1/3 n-113 In 1/2 (n2//) corresponding to a sampling frequency of 1/ = n 1/3, As expected, the estimation bias depends on the singular value gap of the resolvent operator restricted to the RKHS, given by Ry = Rn Finally, note that metric distortion of eigenfunctions can be estimated from data, allowing for spectral bounds with empirical biases for each eigenpair (see Sec. 6).
In conclusion, while TO methods can learn the IG's eigen- functions by learning As: e I. there are currently no theoretical guarantees for small At. This motivated methods that learn the IG directly. Table contrasts these methods with our contribution, which is applicable to more general settings and guarantees learning the leading eigenvalues with lower computational complexity under mild conditions. Compared to (Kostic et al., 2024a), which relies on the Dirichlet form, our estimator is more general (c.g., applica- ble to underdamped Langevin diffusion) and offers linear complexity in terms of state dimension, albeit with lower bias and higher variance. For details, see Sec. 6.
Putting both analyses together, we bound the excess risk of G = , in the operator norm. That for the the operator with probability at least - 5 over samples drawn at frequency 1/A., we obtain that
norm error E(ā) S. G 1-222
In 3/2 (n/0)
E(Gmm)) S
1 Al
(16)
MV /
holds for large enough n = 28, where the regularization parameter is chosen as 7 xxx 1/(n
5

---

Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
5. Learning algorithms
Algorithm 1 Primal LaRRR
Require: dictionary of functions (z))EE(N;; hyperparame-
In this section we assume the access to the dataset D, (d)-1)=ejn) obtained by sampling the process (X)) at some sampling frequency 1/åt for At being typically small in order observe all the relevant scmels of the process. To simplify analysis, we will assume stationarity, i.c. Xo ~ T,, and, hence. Xمa= ~ т. Clearly, to derive empirical risk we need to replace the expectation in (10) with the empirical mean, which leads to estimating cross- covariance operators by their empirical counterparts
ters μl > 0, Y > 0 and r E [n].
1: Compute Z. [=(xr) ... 1=( E RN**
if self-adjoint then {using Toeplitz matrix (19))
3: Symmetrize M < (M+M )/2
4: end if
5: Solve eigenvalue problem HH U/ 0;= U2 i E [r]. where 2 = and C, 1 ZZ ty
6: Normalize U2 + 0/((,, C, [r]
7: Form V, f0 ... |r3 E RNx
Hm Ef Ei (II)) & ((z+4)),
my
(18)
8: Compute eigentriples (v,, w uf)) of V HV,
HV,,u, !u and h4=z V 9: Construct 0 THV
noting that for the time-lag jAl we can only observe n-j pairs from the joint distribution PaDc. that is @ j-0...,,
10: Compute eigenvalues X; = = 1-11///
Ensure: Estimated eigentriples (7e, Di, kxE of L
TaA:- -
Therefore, we can construct the empirical RRR estimator of Rm as G " f/ In particular, when = n, it coincides with the standard ridge regression estimator,
Alternttvely. to derive dual Alg. 2, applicable to dimensoadi H. we perform computations in "sample" space, relying on the reproducing property h(z)- ((x)) and kernel Gram matrices K (1 E R*X" and K, +
In the reminder of this section we show how to compute the estimator and its eigenvalue decomposition in both settings, when the finite dictionary of N features spans H. (Algorithm 1) and when 7 is infinite dimensional RKHS (Algorithm 2). To derive them we follow general construction of vector- valued RRR estimator developed in (Turri et al., 2023), detailed in App. B. To that end, recall the definition of the sampling operator S: H-+ RD and its adjoint S*: R" -
A Algorithm 2 Dual LaRRR
Require: kernel k; hyperparameters μ 0, y > 0 and
r E [m].
1: Compute K= n ( my =;)l ER
2: if self-adjoint then {using Toeplitz matrix (19)
3: Symmetrize M < (M+M )/2
Sh (h(x+1)))) and s====== -= (iEm
4: end if
5: Solve eigenvalue problem Ky K,,,,ōūKK, K,KKus i E [r].
nー 3. (11) CR" is the standard basis. Then, using (18), we ob- tain that Am=3 MSS where MEIR**" is a Toeplitz matrix (i.e. has constant diagonals) given by
implying that S- 133, where
where K, MKM
6: Normalize uj + u4/(u,7 KK, u,) , E [r]
7: Form U, = [w] I.. ur] E Ruxr and V KU, Compute eigentriples (1,, wt, wx) of V MV,
9: Construct g. = S- M V, and h4 = -
{ (nmns)(n-i)) {
i E |n).0siie, (19) otherwise.
Mra4j
10: Compute eigenvalues 2, = μ 1/V, Estimated eigentriples (A, 9., hx) of L
Ensure:
When the process is time-reversal invariant, meaning that Te, and consequently Rm, are self-adjoint, we can enforce symmetry in the empirical objects by estimating Tj RV (Tjar + Tjjr). which can be done at no cost by replacing M with [M++TT. Consequently, both formulations of the algorithm solve a symmetric eigenvalue problems, resulting in real eigenvalues and avoiding additional numerical errors.
In both algorithms, the most expensive computation is in line 5. Naive computations result in cubic complexity w.r. feature dimension N (primal) or sample size n (dual). How- ever, using classical iterative solvers, like Lanczos or the generalized Davidson method to compute the leading eigen- values of the generalized eigenvalue problem, when r & n this cost can significantly be reduced, c.f. (Hogben, 2006). Namely, assuming that complexity of computing a kernel is linear in d, we obtain the complexity of primal LaRRR to be O(rr V d) nN V r(ne N2)), while the complexity of dual one is O( V d) d))²)
In dines--dimeniional H. we have that (<))=2(x) 2(-), where == ,N]] =N] is a vector of features that span
H. that is H = v. 2 E R. Thus, operator (18) becomes isometrically isomorphic to a matrix computed by
replacing o by 2. Thus, estimator Gin,, can be expressed as a N. x N matrix in basis (2))se(N)-
In light of (8), even when the SDE in (6) is unknown, Al- gorithms and 2 enable the construction of approximate
6

---

Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
solutions from a single (long) simulated trajectory by esti- mating the dominant spectrum of L. Specifically, they allow approximating conditional expectations
(BK) Boundedness. There exists cH>0 such that ess sup w-wa
[(((())) ScH, i.e. L (2 (2,
(RC) Regularity. For some E (0.2) there exists
Ca > 0 such that H H I (Ca/p) ( with H₂ Hy=StRySs= Jo Ttt-
E(((X.) | [o==]= x Eietr) ext (Tr () (x)), (20)
where (ū,, h) can be computed on the training set via the kernel trick, see (Kostic al., 2022). This approach is particularly interesting for high-dimensional state spaces, where classical numerical methods become unfeasible due to the curse of dimensionality, making data-driven methods a key tool in fields like molccular dynamics, (Schütte et al., 2023). We prove in Sec. 6 that the precision of our method does not depend on the state dimension, but on intrinsic effective dimension of the process, which, together with its linear complexity w.r.t. d, makes it an attractive approach in such problems. Finally, we remark that (20) enables forecasting of full state distributions and not just the mean (e.g. I can be an indicator function), noting that this formula extends to all 2 functions, at the price of an additional projection error, which leads to predicting the evolution of distributions, c.f. (Kostic al., 2024b).
(SD) Spectral Decay. There exists 3 (0, 1] and cg > 0 s.t. Xy(C)) S j for all j E J.
These assumptions, which originate from the ofate-of-the-art statistical learning theory for regression in RKHS (Fischer & Steinwart, 2020), have been extended to TO regression (Li al., 2022) and to learning self oeff-adjoint IG of diffusions (Kostic et al., 2024a). Condition (BK) ensures that H C C2 while (SD) quantifies the regularity of H. Similar to the regularity condition in (Kostic et al., 2023), (RC) quanti- fies the relationship between the hypothesis class (bounded operators in H) and the object of interest, Ry. Specifically, (RC) holds if L has eigenfunctions in the a-interpolation space between H and (x). fi E H for all E N, then a - 1 (see App. C.I). Since the worst-case bound on metric distortion in Prop. 6.1 depends
the estimator's 's norm, must be restricted to [1, 2] to avoid vacuous bounds, though this restriction isn't 't needed for empirical metric distortions. We present a E [1, 2] here and analyze < 1 in App. C.6.
6. Statistical learning guarantees
Theorem 6.2. Let L be sectorial operator such that Let (BK), (RC) and (SD)
We derive now statistical bounds for estimating IG's re- solvent using the Laplace transform-basse Reduced Rank Regression algorithm (LaRRR). We then derive learning rates for IG's eigenvalues and eigenfunctions, assuming the RKHS is generated by universal kernel k, hence using Alg. 2 to compute Ci
W, (-=->a(L+L")/253 hold for some [12] and 3E(0, 1), respectively, and ((Im({x))=(((z)) Given 8E(0, 1) and rE[n], let
In (n/5)
a+o
In (n/8)) 21++8)
(0)=
(21)
n Atw.
W.
We start with the following auxiliary result, essentially proven in (Kostic et al., 2022). It shows that estimated cigenvalues in the operator regression are guaranteed to lie in the c-pseudospectrum Sp, of the true operator (union of all spectra of perturbed operators), where depends on the operator norm error E(G)=||R,, S.,S--S and the metric distortion r((h)= := ||/]] |||||//Il //Allce. h E H. of es- timated eigenfunctions. To obtain the result, the latter is either uniformly bounded or empirically estimated.
11
n2 W.W.)
Al= and 1/E=o(=ã); then there exists a constant e> 0, depending only on H and o - (Ry,S,) O. such that for large enough n²r with probability at least 1 - 8 in the draw of Dr it holds that
8(Gmms)) S mmxx (ōrt1, (R + (8). (22)
Let G Apn- Gp--0,1/1CC Halr and
Proposition 6.1. Let G = Eie(r) is h4 & qi be the spectral decomposition of G: H. + H. and denote the empirical metric distortions as f/ = ||h.ll/||Shaal| i E [r]. Then for
Proof sketch. Gi - 1/2 [/2 1/2 " and set for brevity |I-I1 =
I1 1111x-24 Then we have
every p > 0, At > 0, l 2 and i ie [r].
E(Gmm) S R 1| -G Guaill
(B.I-R.) 1I S f max (E(G))(h))
E S.)
(1 regularization bias
(I) rank reduction bias
+ S=(( Juy Gm (III) integration bins
| |||ss Sa (G - G Gmma
implying that 0ú belongs to p:-sseuospeccrum of Ryy
(IV) estimator variance
To bound the operator norm version of the excess risk, E(G)= sumptions that quantify the complexity of learning problem and suitability of the chosen RKHS:
x we make the following as-
Regularity assumption (RC) guarantees that (I) S 0/20/2 From definition of RRR, we immediately get (II) S 0rt1 (R) Applying proposition in App. C.2.4 yields
#
7

---

Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
( S Al. Results of App. C.3 gives the control on (IV). Hence get w.p.a.l. 1 1- 6,
assumptions above for clarity. In their equation (24), they obtained:
E E5 2 (R) "V (n--)) At w. 78²
(n()
D X
@rtz(Bp) Sn" 2iiee In(5-)).
| 1 -
(x((
noting that the same result holds when 0 Sx is re- placed by Ø,+1 computed in line 5 of Algorithms 1-2. By balancing with respect to first and then with respect to At, we derive the final bound. D
Note that their bias term (V Or (By)) ( (Ry) is larger than ours (a,+42(R)))/ (R₂)), )) while their variance term (x n -=/2((++))) is smaller than ours (xx n=-0(3a-20)) This is expected, as their method is restricted to self-adjoint IGs exploiting their structure, while our results are structure agnostic and apply to the broader class of sectorial IGs.
Since Prop. 6.1 transforms the estimation error via metric distortion into the pseudospectral perturbation level, it is starting point to derive estimation of cigenvalues and eigen- vectors of Ry, and consequently of L. The quality of such bounds depends on the properties of L, in particular on the conditioning of eigenvalues, that is on the angles between its eigenfunctions, or more generally its spectral projectors. The nicest case is for normal generators, where we derive the following.
7. Experiments
In this section, we demonstrate that our LaRRR Algorithms and 2 successfully recover the eigenvalues and eigenfune- tions of the process's ' IG. We present both 1D and 2D experiments, where we measure error against the ground truth, as well as a higher dimensional molecular dynamics experiment, where we obtain results consistent with the independent study of Mardt et al. (2018). Finally, to em- pirically support our claims in Tab. 1, in Fig. 5 in App. D we provide an additional experiment, showcasing that our method, despite not being pyssicc-inffrmme as (Kostic et al., 2024a), does not suffer from spurious estimation of IG's eigenvalues as baselines (Cabannes & Bach, 2024; Hou et al., 2023) do. The implementation of the methods, as well as all experiments are available in the GitHub repository.
Theorem 6.3. Under the assumptions of Thm. 6.2, lel (At, gi, hy) ke(r be the output of Algorithm 2. If L=L==L" then for large enough nzr with probability at least 1 6 in the draw of Dn,
1 Se n, and |
|
< 2
1 ---l
24 gapp ,,+++
where =6ō,-iã A (R,I(
Sx)
F. = Szh; [|S,halle and gap; is the difference between i-th and (i + 1)-th eigenvalue of R, E [r].
Without further details on the rich theory of spectral pertur- bations, note that in non-normal setting, the conditioning of eigenvalues typically comes as multiplicative factor, e.g. Bauer-Fike theorem (App. A.3).
Non-normal process in 2D. We evaluate the performance of our primal LaRRR algorithm on a one-dimensional Ornstein- Uhlenbeck with (non-normal) 2 x 2 matrix. As features we use 1000 random Fourier features, and test if we can recover eigenvalues of the non-normal linear drift. In Fig. I we compare to TO RRR estimator over ten trials and for two different dime-discrttizations. While high variability in the estimation for both methods is expected due to sensitiveness of the eigenvalues perturbations due non-normality, we can observe that LaRRR achieves better estimation.
Comparison to other approaches. Assume for sim- plicity that the RKHS aligns perfectly with 22 (X) and that I=="", meaning that o1(Aa1S))=0yAAa) for all j5r+1. TO methods
and ( learn Ast = eArl, which shares the same eigen- functions as the generator L. Analysis from (Kostic et al., 2023, Thm. 3 and Eq. 8) for TO ensures that 22 A, denoting the empirical estimate of this product from data. However this bound becomes vacuous as At-+ since (Ast)) but let.AáeeArát | > 0 for finite sample size. In contrast, our generator-based approach is not sensitive to time-lag, as gap, (R,) is independent of At, guaranteeing recovery of eigenfunctions even At- which is crucial for capturing fast dynamics. We note that this analysis of the bounds indeed captures the reality of the both estimators in practice see Fig. 2 bellow and Fig. 4 in App. D.
Overdamped Langevin dynamics in 1D. We evaluate the performance of our dual LaRRR algorithm on a one-dimensional Langevin process with a triple-well poten- tial (Schwantes & Pande, 2015). Specifically, we simulate the dynamics with inverse temperature 3=1 and potential U (x))=4(x++088 Ree--02ee -004 -00e-00000--000000 50(x=0.5)) (.tte--..5++0.... +0.5e- -500=+000)))) that consists of three Gaussian-like wells located at x E {-0.5, O, 0.5). with a bounding term proportional to x to confine the equilibrium distribution primarily within the interval [1. We generate 10 independent trajectories from this process and apply LaRRR to each trajectory in order to estimate the leading eigenvalues of the generator. By fitting 10 separate models, we assess the distribution of the relative errors in the eigenvalue
[[eap,
We now compare our method to the physics innfore- ap- proach of (Kostic et al., 2024a), using the same simplifying
8

---

Laplace Transform Based Low-Complcxity Learning of Continuous Markov Semigroups
Estimated Aup for d:=0.01
Estimated Aup for cb=0.001
et al., 2017) to train and validate the model, respectively, The recovery of metastable states corresponding to the lead- ing two (non-trivial) eigenfunctions of IG is shown in Fig. 3: the eigenfunctions' values (color) are represented in the 2D plane of dihedral angles. Note that approximately con- stant values (red and blue) reveal two metastable states that align with the taateoofthh-aa expert knowledge in molecular dynamics, (Wehmeyer & Noé, 2018).
0.30
0.50
imaginary hart
0.25
0.25-
-3.5
2%
-0.50
-0.50
-2.0
-1.3
-1.0
-0.5
Reel Part
heal Part
Figencshies
LaRRR
PAR
Figure/ The eigenvalues of the (non-normal) 2 2 drift matrix in Ceatein-Uhlenbeek process are estimated via primal LaRRR and TO RRR methods with 1000 random Fourier features in ten trials for two different siiee-disrretiztions. Number of samples for Al le- 2 is 1e4, while for At le- - is le5. This large number of samples is anticipated due to additional sensitivity of the true eigenvalues perturbations.
First Eigenfunction
Second Eigenfunction
2-
-2
-3
estimates, defined as AAs-AA1 for each eigenvalue As Dual algorithm with RBF kernel is assessed in Fig. 2. Similar results obtained by the primal version using random Fourier features is shown in Fig. 4 in App. D.
Figure 3. The first two trivial eigenfunctions of the alanine dipeptide in water trained on a 250 ns simulation with At = 50 ps and displayed on an independent 250 ns test simulation. The color of the points corresponds to the values of the cigenfunctions.
Comparison Dual T aRRR to for ectimating 12
8. Conclusion
2.00
1.75
Relaiize srmt 11, [
We presented a first-ofiits-kind method to learn continu- ous Markov semigroups, offering both theoretical guaran- tees at any time-lag and linear computational complexity in the state dimension, enabling efficient exploration of high-dimensional complex systems. Notably, our method applies to a broad range of Markov processes previously unaddressed, and overcomes the problem of -O-bedd meth- ods failing to capture slow dynamics when trained on data with high sampling frequency. The main limitation of the current results lies in the assumption of uniform sampling of the (full) state of the system. While the theory can be seam- lessly adapted to multiple observations sharing the same non-uniform sampling, an important challenge is to extend it to non-uniformly sampled single trajectory data, as well as to the case of partially observed systems.
0.
1.00-
0.0
A
A
A
0.75
0.0
025-
0.00
-2.25
-2.00
-1.75
-1.30
-1.25
-1.00
-0.73
-0.39
Figure 2. Comparison of the LaRRR dual algorithm the TO RRR one on the task of estimating the slowest timescales of the process. As predicted by our theory, while LaRRR (blue) not affected by decreasing the time-lag At, the error of TO method (red) explodes as At ++0 0. The main figure shows three quartiles relative error for A2 across 10 independent trajectories of a ID Langevin process on trille-well potential, while the inset figure shows distributions for the lcading threc eigenvalucs for At- 0.005.
Acknowledgments. This work was supported by the EU Project ELIAS (grant No. 101 101 120237), and by the Euro- pean Union NextGenerationEU and the Italian National Recovery and Resilience Plan through the Ministry of Uni- versity and Research (MUR), under Project PF0000013 CUP J53C22003010066, KL acknowledges support from the French National Research Agency for the DECATTLON project LANR-24-CE40-333 We thank the anonymous reviewers for their insightful and valuable feedback.
Alanine Dipeptide. Alanine dipeptide in water serves as a benchmark for studying dynamics due to its simplicity and metastability, where the system predominantly occupies several metastable states, that can be identified using the dihedral angles @ and db. Under reasonable assumptions, the long-term dynamics can be treated as Markovian by integrating out the water degrees of freedom. In this study, we apply our method to alanine dipeptide, using the inter- atomic distances between heavy atoms as input. Notably, the state space dimension is d = 45, which is intractable for the kernel based generator learning method of (Kostic et al., 2024a). We use two independent trajectories from (Nüske
Impact Statement
This paper aims to advance the field of Machine Learning. While the work may have societal implications, none require specific mention here.

---

Laplace Transform Based Low-Complexity Learning of Continuous Markoy Semigroups
References
Hou, B., Sanjari, S., Dahlin, N., Bose, S., and Vaidya, U. Sparse learning of dynamical systems in RKHS: An operator-theoretic approach. In International Conference on Machine Learning. 2023.
Andrew, G., Arora, R., Bilmes, J., and Livescu, K. Deep canonical correlation analysis. In International Confer- ence on Machine Learning, pp. 1247-1255. PMLR, 2013.
Karatzas, I. and Shreve. S. Brownian Motion and Stochastic Calculus. Graduate Texts in Mathematics. Springer, 1991.
Bakry, D., Gentil, I., and Ledoux, M. Analysis and Geome- try of Markov Diffusion Operators. Springer, 2014.
Kato, T. Perturbation Theory for Linear Operators. Classics in Mathematics. Springer, 2012.
Bevanda, P., Beier, M., Kerz, S., Lederer, A., Sosnowski, S., and Hirche, S. KoopmanizingFlowws diffeomorphi- cally learning stable Koopman operators, arXiv preprint arXiv:2112.04085 2021.
Kloeden, P. E.,, Platen, E., Kloeden, E., and Platen, Stochastic Differential Equations. Springer, 1992.
Bonati, L., Piccini, G., and Parrinello, M. Deep learning the slow modes for rare events sampling. Proc. National Academy of Sciences, 118(44):e21 13533118, 2021.
Klus, S. and Djurdjevac Conrad, N. Koopman- based spec- tral clustering of directed and time-evolving graphs. J. Nonlinear Sci., :3(1):8 2023.
Bradley, R. Introduction to Strong Mixing Conditions. Num- ber vol. 1,2,3. Kendrick Press, 2007.
Klus, S., Schuster, I., and Muandet, K. Figendecomposi- tions of transfer operators in reproducing kernel Hilbert spaces. Nonlinear Sci., 30(1):283-33)) 2019.
Brunton, L., Budišié, M.. Kaiser, E., and Kutz, J. N. Modern Koopman theory for dynamical systems. SIAM Review, 64(2):229-340, 2022.
Klus, S., Nüske, F., Peitz, S., Niemann, J.-H., Clementi, C., and Schüttee C. Data-driven approximation of the Koop- man generator: model reduction, system identification. and control. Physica D, 406: 4406:132416, 2020.
Cabannes, V. and Bach, F. The Galerkin method beats graph- based approaches for spectral algorithms. In Proceedings of The 27th International Conference on Artificial Intelli- gence and Statistics, volume 238, pp. 451-459 2024.
Kostic, V., Novelli, P., Maurer, A., Ciliberto, C.,, Rosasco, L., and Pontil, M. Learning dynamical systems via Koopman operator regression in reproducing kernel Hilbert spaces. In Advances in Neural Information Processing Systems, 2022.
Crouzeix, M. and Palencia, C. The numerical range is a (1+2)-spectral set. SIAM Journal on Matrix Analysis and Applications, 38(2):649 655, 2017.
Kostic, V. R., Lounici, K., Novelli, P., and Pontil, M. Sharp spectral rates for Koopman operator learning. Ad- vances in Neural Information Processing Systems, 2023.
Das, S. and Giannakis, D. Koopman spectra reproduc- ing kernel Hilbert spaces. Applied and Computational Harmonic Analysis, 49(2):573 607,5777--77550 2020.
Kostic, V. R., Lounici, K., Halconruy, H., Devergne, T., and Pontil, M. Learning the infinitesimal generator of stochastic diffusion processes. In Advances in Neural Information Processing Systems, 2024a.
Davis, C. and Kahan, W. M. The rotation eigenvectors by perturbation. iii. SIAM Journal on Numerical Analysis, 7(1):1-46, 1970.
Devergne, T., Kostic, V. R., Parrinello, M., and Pontil, M. From biased unbiased dynamics: An infinitesimal gen- erator approach. In Advances in Neural Information Pro- cessing Systems, 2024.
Kostic, V. R., Lounici, K., Inzerilli, P., Novelli, P., and Pontil, M. Consistent long-term forecasting of ergodic dynam- ical systems. In International Conference on Machine Learning, pp. 25370-25395. PMLR, 2024b.
Fan, F,, Yi, B., Rye, D., Shi, G., and Manchester, I. R. Learning stable Koopman embeddings. arXiv preprint a0X/v.2110.06500, 2021.
Kostic, V. R., Novelli, P., Grazzi, R., Lounici, K., and Pontil, M. Learning invariant representations of time- homogeneous stochastic dynamical systems. In Interna- tional Conference on Learning Representations, 2024c.
Fischer, S. and Steinwart, I. Sobolev norm learning rates for regularized least-squares algorithms. Journal of Machine Learning Research, 112205):1-38, 2020.
Kutz, J. N., Brunton, S. L., Brunton, B. W., and Proctor, J. L. Dynamic Mode Decomposition. SIAM, 2016.
Hogben, L. (ed.). Handbook of Linear Algebra. CRC Press, 2006.
Laio, A. and Parrinello, M. Escaping free-energy minima.
Proc. Natl. Acad. Sci. USA, 9922 2002.
10

---

Laplace Transform Based Low-Compicxity Learning of Continuous Markov Semigroups
Li, Q. Dietrich, F., Bollt, E. M., and Kevrekidis, I. G. 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.
Schüttte, C., Huisinga, W., and Deufthard, P. Transfer opera- tor approach to conformational dynamics in biomolecular systems. In Ergodic Theory, Analysis, and Efficient Simu- lation of Dynamical Systems. Springer, 2001.
Schüitte, C., Klus, S., and Hartmann, C. Overcoming the timescale barrier in molecular dynamies: Transfer oper- ators, variational principles and machine learning. Acta Numerica, 32:517-673, 2023.
Li, Z., Meunier, D., Mollenhauer, M., and Gretton, A. Op- timal rates for regularized conditional mean embedding learning. In Advances in Neural Information Processing Systems, 2022.
Schwantes, C. R. and Pande, S. Modeling Molecular Ki- neties with tICA and the Kernel Trick. Journal of Chemi- cal Theory and Computation, 1(2))600-608, 2015.
Lusch, B., Kutz, J. N., and Brunton, S. L.. Deep learning for universal linear embeddings of nonlinear dynamics. Nat. Commun., 91),, 2018.
Shmilovich, K. and Ferguson, A. L. Girsanov reweighting enhanced sampling technique (grest); On-the-fly data- driven discovery of and enhanced sampling in slow col- lective variables. J. Phys. Chem. A, 127((5))3477-311, 2023.
Majda, J., Franzke, C., and Crommelin, D. Normal forms for reduced stochastic climate models. Proc. Natl, Acad. Sci. USA, 106( (10)::00)393666
2009.
Mardr, A., Pasquali, I.., Wu, H., and Noé, F. VAMPnets for deep learning of molecular kinetics. Nat. Commun., 9(1,, 2018.
Tian, W. and Wu, H. Kernel embedding based variational approach for low-dimensional approximation of dynami- cal systems. Comp. Meth. Appl. Math., 21(3):635-65) 2021.
Marrouch, N., Slawinska, J., Giannakis, D., and Read, L. Data-driven Koopman operator approach for computa- tional neuroscience. Annals of Mathematics and Artificial Intelligence, 88(11): 1555 1173, 2020.
Trefethen, L. N. and Embree, M. Spectra and Pseudospec- tra: The Behavior of Nonnormal Matrices and Operators. Princeton University Press, 2020.
Nüske, F., Wu, H., Prinz, J.-H., Wehmeyer, C., Clementi, C., and Noé, F. Markov state models from short non- equilibrium simulations-anally analysis and correction of esti- mation bias. J. Chem. Phys., 146(9):094104, 2017.
Turri, G., Kostic, V., Novelli, P., and Pontil, M. A random- ized algorithm to solve reduced rank operator regression. arXiv preprint arXiv:2312.
2023.
Ostrow, M., Eisen, A., Kozachkov, L., and Ficte, I. Beyond geometry: comparing the temporal structure of computa- tion in neural circuits with dynamical similarity analysis. In Advances in Neural Information Processing Systems, 2023.
Wehmeyer, C. and Noé, F. Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinet- ics. Chem. Phys., 148(24):2417033 2018
Williams, M. O., Rowley, C. W., and Kevrekidis, I. G. A kernel-based method for data-driven Koopman spectral analysis. J. Comput. Dyn., 2(2):247-2655 2015.
Pillaud- Vivien, L. and Bach, F. Kernelized diffusion maps. In The Thirty Sixth Annual Conference on Learning The- ory. pp. 5236-5259. PMLR, 2023.
Yosida, K. Functional Analysis. Springer, 2012.
Yu, B. Rates of convergence for empirical processes of stationary mixing sequences. Ann. Probab., 22(1):94- 116, 1994.
Quarteroni, A., Sacco, R., and Saleri, F. Numerical Mathe- matics, volume 37. Springer, 2010.
Rahimi, A. and Recht, B. Random features for large-scale kernel machines. Advances in Neural Information Pro- cessing Systems, 20, 2007.
Zabezyk, J. Mathematical Control Theory: An Introduction. Springer, 2020.
Zhang, W., Li, T., and Schütte, C. Solving eigenvalue pdes of metastable diffusion processes using artificial neural networks. J. Comput. Phys., 465:111377, 2022.
Reed, M. and Simon, B. Functional Analysis. Academic Press, 12 1980.
Schütte, C. and Huisinga, W. Biomolecular conformations can be identified as metastable sets of molecular dynam- ics. In Special Volume, Computational Chemistry, vol- ume 10 of Handbook of Numerical Analysis, pp. 699-744. 2003.
Zwald, L. and Blanchard, G. On the convergence of eigenspaces in kernel principal component analysis. In Advances in Neural Information Processing Systems, 2005.

---

Laplace Transform Based Low-Compicxity Learning of Continuous Markov Semigroups
Supplementary Material
notation
meaning
||
notation
meaning
A
minimum
V
maximum
R
real part of complex number
S
imaginary part of complex number
[
r-truncated SVD of an operator
I
identity operator
HS
(H. G)
space of Hilbert-Schmidt operators H. +
Br
set of finite rank-r operators on H
[A|
operator norm of an operator A
All
Hilbert-Schmidt norm of an operator A
F(A)
numerical range an operator A i-th singular value of an operator
V
Hessian of function
o.(.)
Ac(-)
i-th eigenvalue of an operator
X
state space of the Markov process
(Xe)120
time-homogeneous Markov process
p
transition kernel of the Markov process
π
invariant measure of the Markov process
drift of the Itô process
b
diffusion of the Itô process
(X)
1.2 space of functions on X' w.r.t measure π
A
transfer operator on 2 (2)) for time-step t
(X
Sobolev space w.r.t. measure on X'
L
generator of the semigroup on (X)
1
shift parameter
Al
time discretization
R
resolvent of
Rm
approximated resolvent
Ry
resolvent operator restricted to
Rn(L))
resolvent of L
k(x, y)
kernel
Φ
canonical feature map canonical injection H. + (
H
reproducing kernel Hilbert space
S
1
function in (X) constantly equal to 1 true risk
?
regularization parameter
R
E
operator norm error
Rex
excess risk, i.c. HS norm error
Ro
irreducible risk
R
approximated risk Bochner integral
R
empirical risk
approximated Bochner integral
sampling operator w.r.t. (X))
P
joint invariant measure the Markov process
C
covariance operator on L (X)
C
empirical covariance operator on C (X')
C
regularized covariance operator on (x))
C
regularized emp. covariance operator on L (X)
Tr
cross-covariance operator at time step /
T.
emp. cross-covariance operator at time step
Hμ
integrated cross-covariance
Hm
empirical aggregated cross-covariance
Hm
aggregated cross-covariance
M
Toeplitz matrix associated to Hm regularized kernel Gram matrix
K
kernel Gram matrix
K.
G
population estimator of Ry on H
G
empirical estimator of Ry on H
Gay
population KRR estimator
Gmn
empirical KRR estimator empirical RRR estimator
G
population RRR estimator
G m.
Gmn
approximated population RRR estimator
P
approximated projector
P
spectral projector
P
empirical spectral projector
1)
metric distortion
empirical metric distortion
A
i-th generator eigenvalue
X
i-th eigenvalue of the empirical estimator
fi
i-th generator right eigenfunction in (x) fr
i-th empirical right cigenfunction in C (x')
gi
i-th generator left eigenfunction in (x)))
i-th empirical left eigenfunction in (X)
h,
i-th generator right eigenfunction in (x)
h;
i-th right empirical eigenfunction
i-th empirical eigenvalue
boundedness constant
regularity parameter
Ca
regularity constant
3
spectral decay parameter
Cg9
spectral decay constant
T
embedding parameter
6-
embedding constant
12

---

Laplace Transform Based Low-Compicxity Learning of Continuous Markov Semigroups
The supplementary material is organized follows. App. A reviews key results on Markov semigroups, RKHS operator regression, and spectral perturbation theory. App. B covers prior reduced rank regression methods and extends LaRRR to non-uniformly sampled data from multiple trajectories. App. C presents proofs of Theorems 6.2 and 6.3, including guarantees for LaRRR with non-unnnonnn uniform sampling. The table above summarizes the notation.
A. Background
A.1. Markov semigroups
Markov processes, where the future depends only on the present, are key to modeling physical, biological, and financial systems. This class includes Itô diffusions, time-dependent singular drift processes with continuous paths, and jump processes like Poisson, Lévy, and Hawkes. Here, we focus continuous Markov processes, primarily Itô diffusions. All Markov processes can be defined both as time-dependent random functions and by their laws, via measures on path space. A key tool is the infinitesimal generator (IG), a linear operator on observables.
We provide here some basics on operator theory for Markov processes. Let X' C Rd Ed and (Xt)ER, be a X'-valued time-homogeneous Markov process defined on a filtered probability space (2, (Fi)EE,, where Jr = a(x. S t) the natural filtration of (Xt)rER+. The dynamics of a continuous-time Markov process X described by a family of probability densities (Pt)tER+
P(X, E B|X,- =x) 1 pc v)dy,
(23)
and transfer operators (TO) (Ar)Ee,, such that for all t E R4, E B(x). x E X and measurable function f : - R,
= A.f = f (((( X.
| =
(24)
The transfer operator essential for understanding the dynamics of X. We examine its action on (X)) noting the presence of invariant measure v, which satisfies Ar n for all i. E R4. In theory of Markov processes, the family (At)eER, is referred to as the Markov semigroup associated to the process X. The process X is then characterized by the infinitesimal generator (IG) L: (²(2 + L (X') of the family (Ar)ER, defined by
A t
L= lim
(25)
In other words, L characterizes the linear differential equation A,A,J LAi satisfied by the transfer operator. The spectrum of the IG can be difficult to capture due to the potential unboundedness of L.. To circumvent this problem, one can focus on an auxiliary operator, the resolvent, which shares the same eigenfunctions as L and becomes compact under certain conditions. The following result can be found in Yosida's book ((Yosida, 2012), Chap. IX) : For every μ >> 0, the operator (uT - T.) admits an inverse R,, = (μT T.)-1 that is continuous operator on X' and
(μI - (uI-L)-1 J -ut
e-' Acdt.
The operator Lp is the resolvent of L and the corresponding resolvent set of L is defined by
P(L) = {u+ E C (μI L) is bijective and L is continuous].
In fact, P(L) contains all real positive numbers and (μl - L)-1 operator (see (38) below) is uniformly bounded outside a sector containing the spectrum. Its spectral decomposition writes
1 is bounded. In particular, the resolvent of a sectorial
= ΣiEN &
(26)
with eigenvalues (At)ien C C and corresponding left and right eigenfunctions fú, gi L2, respectively. We detail the two examples of processes with sectorial IG discussed in the paper: the overdamped Langevin process and the Maastein-Uhlenbed process.
Example 2.1 (Overdamped Lagenvin detailed) Let kb and T Ry. The overdamped Langevin dynamics of particle potential V IRd -> R is given by the SDE: dx,- - -VV +
+ /2(kTT/77)dW) and X, = x,
13

---

Laplace Transform Based Low-Complcxity Learning of Continuous Markov Semigroups
where 7, bb. and T are the friction coefficient, Boltzmann constant, and system temperature, respectively. The in- variant measure is the Boltzmann distribution π m(dx) Z--e- (-te-V(e)/(ka)ddd da, where Z is normalizing constant. Its infinitesimal generator is Lf -7 -VV' Vf + (k.TH))&f, for f E (2) Since [(-Lf)g du -s [ ((ITT e (//T//-/via g(()dz (kaT/>) f V dm (gg dr, the generator is self- adjoint. In dissipative systems, the IG L is sectorial, with its spectrum usually in the left half-plane For confining potentials, the spectrum is discrete, featuring O as the largest eigenvalue, which corresponds to the distribution π.
Example 2.2 (Ornstenn-Uhlenbect detailed) The Ornstein-
(OU) process with drift is governed by the SDE dx,- AX, + BdW, with Xo = x, where E Rdxd is the drift matrix and B E Rdxd is the diffusion matrix. This models systems like the Vasicek interest rate, damped harmonic oscillators, and neural dynamics, where fluctuations return to equilibrium. If the real parts of A's eigenvalues
negative, the OU process has invariant Gaussian distribution with covariance Ex satisfying Lyapunov's ' equation: A& A - BB". Its infinitesimal generator is defined, for E T.2, by I.f (x) V fxx) Ax + 2 [Tr(BT (x)) The IG has a discrete spectrum with eigenvalues related to the drift matrix A, where O corresponds to the distribution and the rest are negative, reflecting relaxation rates.
A.2. Operator Regression in RKHS spaces
Recalling the operator regression problem for learning the resolvent of the generator, for the reader's 's convenience we state that the learning problem is well posed, the proof of which is essentially the same as in (Kostic et al., 2022). Proposition A.1. Given μ 0, ler H C (X) be the RKHS associated to kernel XxX + R such that S E HS (H, C2), and let PM be the orthogonal projector onto the closure of (m{))) ) _ (X). Then for every ε >> 0 there exists a finite rank operator G: H - H such that
R(G)S S | - |R + [I ([T-Px) R, 1H +
Ra
Rex(G) (G)
Consequently, when ls is universal, the excess risk can be made arbitrarily small.
Proof. First, note that decompose
R(G) S ||St|ins
||
|R + +I - S.G.H
R.
(C)
as in (Kostic et al., 2022, Proposition 4) but now applied with additional integration. Next, since S₂ E HS (H. W (X)), according to the spectral theorem for positive self-adjoint operators, S₂ admits an SVD S, = EjEN @ hy. Recalling that [+) denotes the r-truncated SVD, i.e. [[ Ejelr ozlj @ hy. since ||S 11 = Ejsr o7 for every > O there exists r E such that ||St - [S+]rI|n < <10/3 Consequently since all the eigenvalues of L have non- positive real part, (S, [- S |S [s*IrI /p S 8/3. Next since Im (m(P)) S cl(Im(S,)), for any j E [r]. there exists g₁ E H s.t. [[ - Zug;il S $. and, denoting B, = Σje(n Jggj @ h2 we conclude [[PRRR,SssI Sx B, Ils S 8/3. Finally we recall that the set of non-defective matrices is dense in the space of matrices (Trefethen & Embree, 2020), implying that the set of non-defective rank-r linear operators is dense in the space of rank-r linear operators on a Hilbert space. Therefore, there exists a non-defective G E (H) such that ||G -G-AA -.llss < 66((0,, (8+)). So, we conclude
+IS Sr + RR - S ss I |I[IIPw))R, R, Sx + +IRR... [RpSaIillss + |IR
B,)ll
< [[ [1- Pxx Rr Sx IHS +
A.3. Spectral perturbation theory
Recalling that for a bounded linear operator A on some Hilbert space H. the resolvent set of the operator A is defined as
Res(A)= := {} E C: A - A/ is bijective), and its spectrum Sp(A) == C (R)s(A)), let X C Sp(A) be isolated part of spectra, i.e. both and == Sp(A) \ > are closed in Sp( Than, the Riesz spectral projector PA: H + H is defined by
1 2T J
/
PA
(zl A)z
(27)
14

---

Laplace Transform Based Low-Compicxity Learning of Continuous Markov Semigroups
where P is any contour in the resolvent set Res(A) with À in its interior and separating À from μ. Indeed, we have that P3 = PA and H Im PP) @ Ker( (P) where Im Px) and Ker(Px) are both invariant under A and Sp(A))er)) X, Sp(AI = μ.. Moreover, Px + P.. = and Px =
Finally if is compact operator, then the Shesud-S theorem, see e.g. (Reed & Simon, 1980), assures that Sp(T) is a discrete set having no limit points except possibly À = O. Moreover, for any nonzero À E Sp(T), then À is an eigenvalue (i.c. it belongs to the point spectrum) of finite multiplicity, and, hence, we can deduce the spectral decomposition in the form
=
Σ Px
(28)
AES((A) A)
where geometric multiplicity of X, ra- == rank(P.), is bounded by the algebraic multiplicity of X. If additionally A is normal operator, i.c. AA* A A, then PA Px is orthogonal projector for each À E Sp(A) and Px Σ 1i @ 1h, where 14 are normalized eigenfunctions of A corresponding to À and rx is both algebraic and geometric multiplicity of X.
Next we review well-known perturbation bounds for eigenfunctions and spectral projectors of normal compact operators, that is when AA" = A.
Proposition A.2 (Davis & Kahan, 1970)). Let A be compact self-adjoint operator on a separable Hilbert space H. Given L, let be the eigenvalue of A that closest X and let f be its normalized eigenfunction. Ifğ == minf[ā AI>E E Sp(A) (2}} > 0, then sin( (F) - -AA--A-
a pair f) C x H. such that |ل||
Proposition A.3 (Sinclair's theorem, see (Zwald & Blanchard, 2005)). Let A and À be two compact operators on separable Hilbert space. For nonemply index set J C N lef
gap,(A) == min ((>((A) M(AA)llEENN
denote the spectral gap w.r.t J and let Ps and Ps be the corresponding spectral projectors of A and À, respectively, If A is normal and for some ||A - Al < (A), then
AA - A
[P, Pyll
gap,(A)
For non-normal operators, spectral perturbation becomes much more involved. Two core objects in it are the pseudospectrum and the numerical range, which we review next.
Definition A.4 (Pseudospectrum), Given e > O, the t-psendospectrumm of a bounded operator A on H (w.rt. the spectral norm) is set in the complex plane that consists of all eigenvalues of all e-perturbations A (w.r.t spectral norm). Formally, is defined
Sp, (A) {2 E C Sp(A+ + E), for some operator E s.t.
|E| S <
Since, it provides insights on the sensitivity of eigenvalues, pseudospectrum is the tool of choice in the study of non-normal operators. Equivalently, it can be expressed via resolvent
Sp. (4) {= CI [())) se).
implying that the lower bounds on the resolvent's norm lead to bounds on the eigenvalue sensitivity to the perturbations. Whenever the operator is normal, that is commutes with its adjoint, indeed we have that dist( Sp(A)) = minxesp(A)
A| for all z E 0 Sp. (A). On the other hand, when operators are not normal, distance of the eigenvalues of the perturbed operator to the spectrum is typically amplified. A general result showing this is the following.
Proposition A.5 (Bauer-Fike theorem, see (Trefethen & Embree, 2020)), Let A be diagonalizable bounded operator on separable Hilbert space, that is there exists bounded operator X with bounded inverse X such that X-1AX AX is diagonal operator. If B. := {z E C || & {) denotes the e-hall in C, then
Sp(A) + B. C (A) C Sp(A) + Bex(x),
where x(X)= =
||||||--||| is the condition number of X. Consequently, (((()) all =Ed Sp. (p.()).
15

---

Laplace Transform Based Low-Compicxity Learning of Continuous Markov Semigroups
We recommend book by (Trefethen & Embree, 2020) for in depth study of nonormality, pseudospectra and related quantities.
Definition A.6 (Numerical Range). field of values, known also as numerical range, of a bounded operator A on H is set in the complex plane that is closely related to the spectrum of A. Formally, it is defined as
F(A) = (Av, u) (0, -UEH,
Due to Hausdorff's theorem, numerical range is a convex subset of C and that its closure contains the spectrum of A. It is one of the key objects in the study of operator semigroups, since it is known that
max ()))(A)) < ||eta| S e¹ =ax(B((=]=EE
Finally, we conclude this review with a key result on bounding analytic functions of an operator using its numerical range. Theorem A.7 (Crouzeix's s Theorem, see (Crouzeix & Palencia, 2017)). Let A be bounded linear operator on Hilbert space IH, and let f be any analytic function. Then
22) [If(A))l S (1+ VL+
max [f(z)),
2(4)
where cl ((A) denotes the closure of numerical range of A.
B. Methods and their extension to non-uniform sampling
In this section we discuss the Algorithms I and 2, together with their simple extensions to the case of learning from multiple independent copies of the process over a finite time-horizon.
For computationally efficient primal and dual form of the solution of general vector-valued reduced rank regression problem we refer to (Turri et al., 2023, Proposition 2.2). Coupling this result with classical result on solving low rank eigenvalue problems, see e.g. (Kostic et al., 2022, Theorem 2) for application the context TO, we obtain the final form two algorithms.
Now, let discuss the setting where dataset D, (xet))-)))ee))jje(e-- consisting of trajectories sampled at the same (possibly non-uniform) times 0 to < 4 <... < 4.. This is typical when we observe particle systems where n copies of the same process are tracked.
In this setting, learning transfer operators becomes difficult task since, due to diverse time-lags one ends up with the non-convex optimization problems. Instead, learning IG of the process by LaRRR is direct. Namely, we simply have that Gm C, 1/2 12, 1/2 1/2 Hmlr where in this setting
C₂ 1((+1)
Σ SEIn| JE[(+1]
Σ
((IItt-)) 2 (( and Am
Σ mj-1 - H (a 8 (( jEl(+1]
xE[m]
where now mg = _ e uéo my - and my In this case, the sampling operator becomes
44--4-1 e- uty otherwise for 1 Sjse-1.
2
Sh (h( (M vn(e+1)
implying that the matrix M E pn((+1)xn((+11) is defined as
Mijots = { ( + (e+1) m)
E In].0siie
(29)
otherwise.
Therefore, we can readily implement Algorithms 1 and 2 also in this setting of non-uniformly sampled data. Moreover, analyzing such an estimator, as we discuss in App. C.7, from the statistical perspective, becomes an easier task than a single trajectory analysis. This is due to the fact that we observe n iid copies allows the use of classical concentration inequalities without the need to asses the impact of mixing.
16

---

Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
C. Statistical learning theory
In this section, we provide the proof' the main results from Section 6, along with additional discussions for the reader's convenience. We begin by stating the main assumptions, then present results on controlling the approximation error and the estimator's variance, and finally prove the main results on operator norm error and spectral learning bounds. To this end, we first decompose the operator norm error E(G"m.)) S| -2 S, -m,ylltt- of the Laplace transform- based Reduced Rank Regression algorithm (LaRRR) Ginyy as
E(G"mA) S (I ( S, S.in
( -mmm
(0) representation bias
(1) regularization bias
(lI) rank reduction bias
(III) integration bias
St ( G (ma))
(30)
(IV) estimator variance
1C,1/2 where PH is the orthogonal projection in (2) onto the cl(mm(Sx)) GH. C H₂ and G -C,,-/////100 10 1/2 100T///1110001 Hp are
the population KRR and RRR models for the risk (10), respectively, while Gm,y C- and G C |C,1/2 Hm are the population KRR and RRR models for the approximated risk (14), and, for simplicity, we abbreviate ||-|| == || 2-
C.I. Assumptions
In this section we discuss the assumptions that are necessary study the learning problem. They concern the interplay between the resolvent of IG and the chosen RKHS space, which is encoded in the injection operator S and the restriction of the resolvent to the RKHS Zy = S..S..
We start by observing that Sz E HS (H, (X)), according to the spectral theorem for positive self-adjoint operators, has an SVD, i.c. there exists at most countable positive sequence (α,))EJ, where J := {1,2..... } S. N, and ortho-norma systems (Ej)jes and (hz)jes of cl (),)) and Ker(S,) ) respectively, such that S,h, o, and S = ,,hj. j E J. Now, given a 2 0, let us define scaled injection operator S H C (X) as Sr,a := Σee, oft, @ hy. Notice that S 1.,,, while Im ST el cl(Im(S))) Next, we equip Im(Sz,a) with a norm ll-lla to build an interpolation space.
{
[H|a = f (m((s,a)))))) 11
20 00 {
JEJ
We remark that for a = 1 the space [H] is just an RKHS H seen as a subspace of C (x). Moreover, we have the following injections
[Mm.
[?]] + [M] + 1H] (x))
where on 2 1 2 022 2 0.
In addition, from (BK) we also have that RKHS can be embedded into L (X)), i.e. for some E (0,
[1| + [H], - 1 (x) + (x))
According to (Fischer & Steinwart, 2020), if Sa [7(,, + (X)) denotes the injection operator, its boundedness implies the polynomial decay of the singular values of 5n, i.e. o7 (2 j-11/1 E J, and the following condition holds
(KE) Kernel embedding property: there exists T [1, such that
G = || 22 sup
ΣofT)( < +00.
(31)
2~0
JEJ
Assumption (SD) allows to quantify the effective dimension of 7 in ambient space C (2)), while the kernel embedding property (KE) allows one to estimate the norms of whitened feature maps
{(x) C, 1/2 (xx) E H,
(32)
as the following result states.
17

---

Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
Lemma C.1 (Fiischrr& & Steinwart, 2020)). Let (KE) hold for some 7 E [3,1] and c7 E (0, Then,
1-8Y , < 1, and ||E = css sup 202
E 15(x)]1 -
[<(x))] S C07
(33)
B =1
Next, we discuss assumption (RC) that quantifies the regularity of the problem w.r.t. the chosen RKHS. Typical form in which well-specifiedness of the operator regression problem is expressed, see e.g. (Li et al.,, 2022), is to ask that H is invariant under the action of the operator. In our setting, this reads as Im(Zp) C Im(S)). This can be relaxed (or straightened) by using interpolation spaces [H]a = Im an(Sx.)). which leads to (RC). Indeed, according to (Zabczyk, 2020][ 2.2], the condition (RC) is equivalent to Im Zn) = p SR) C Im(S i.e. G := fuSt n.a S R = μuSt n,a H₂ is bounded operator on H and SR,.S,+ = Sa.a Gz- In this case we have that μH,,, = S (HR,, )S, = S G% and, thus, H 3 (I|Ggll lI Cp)+C We note that re-scaling with j is motivated by the fact that | S 1.
We conclude this section discussing the assumptions for the case when RKHS is build as span of some finite dictionary of functions (zy)JE(N), 2z X + R. j [N].
H = hhw = EEENN u= (u=(u,...... uN) E R
(34)
The choice of the dictionary, instrumental in designing successful learning algorithms, may be based on prior knowledge on the process or learned from data (Kostic et al., 2024c; Mardt et al., 2018). The space H is naturally equipped with the geometry induced by the norm || == Σソ= uz. Moreover, every operator A: H can be identified with matrix A E RN*N by Ah, = 2 Au. In the following, we will refer to A and A as the same object, explicitly stating the difference when necessary.
N
In this setting all spaces [H),, are finite dimensional, and, hence, Im C Im(S²) implies also Im( (Zp))) C (m((S,,)) for every a 0. This choice makes (RC) for different a E (0, 2] all equivalent asking for H to contain all eigenfunctions of the generator. Moreover, in this setting we can set
and observe that in the limit y + 0 we have
tr( tr(C,+ + N and | - ess
ess sup (Φ(x)), C (x))) < 00.
C.2. Approximation error analysis
In this section we study the approximation errors, i.c. bias terms in (30).
C.2.1. REPRESENTATION BIAS
The first term is the representation error that one incurs only when the hypothesis space H. is not dense in the true domain C (2)). That is, if one uses universal kernels, such RBF Gaussian kernel k(z, v) uia with some length-scale >> 0, then [I Pxt RI, || = 0. On the other hand, when one uses the finite-dimensional RKHS (34), this term quantifies the loos of information due to the restriction of the model to H. In recent work by (Kostic et al., 2024c) it has been shown how to learn dictionary of functions with neural networks so that the minimal representation error is incurred from perspective of TOs. While, one can use the same method for learning an appropriate kernel for the estimation of the resolvent, an interesting future directing would be to develop representation learning based on the Laplace transform.
C.2.2. REGULARIZATION BIAS
We control this term using the regularity assumption. The result is stated in the following proposition whose proof we omit since it follows exactly the same lines as (Kostic et al., 2023, Proposition 5) with the only difference that (RC) assumption holds for Ry instead of TO with fixed lag-time.
Proposition C.2. Let Gp.-- C- H for > 0, and PH: 2x - (X) be the orthogonal projector onto cl(Im(S.)) (Sx)). If the assumptions (BK), (SD) and (RC) hold, then
(ca/1)c(a-11)) (2-1) a E [1,2],
||G <
and |Pre - SrGurll S (ca/u))
(35)
(ca/H)) aa E (0, 1],
18

---

Laplace Transform Based Low-Complcxity Learning of Continuous Markoy Semigroups
While the previous result gives us the means to study the learning bounds when H. is dense in (2),, whenever the RKHS finte-ddimensional and not learned, the regularization bias is usually coupled with the representation bias and studies the decay of [ St Gy.ll as function of > 0 w.r.t N + oo for the choice of dictionary that forms basis of C (Z') in the limit.
C.2.3. RANK REDUCTION BIAS
To control of this term, observe that
S (G - G,,))|= - 10 H
- 11 P.) - 1r+1(B),
s
(36)
where P. is the orthogonal projector onto the leading right singular space of operator B = =C C7.11211p
This bound on the bias can be further analyzed using the (RC) assumption, as indicated by the following proposition
essentially proven in (Kostic et al., 2023, Propositions and 7).
Proposition C.3. Let B == C.y H let (RC) hold for some a E (0, Then every E J,
(z, - (ca/N)² a/2 a s (b) S oz)(2
CH
(37)
C.2.4. INTEGRATION BIAS
Diversely the previous terms analyzed up to now, the control of the integration bias term truly novel technical contribution in the study of TO and IG data-driven methods. This necessary step allows one connect empirical objects to the population ones by approximating the integral with the finite sum.
We start with the key result showing that, for the class of sectorial IGs, we can control the difference between the resolvent and its sumite-sum approximation via TOs of different times-lags using Crouzeix's 's theorem A.7.
Proposition C.4. Let be a (stable) sectorial operator with angle 0 E π 1/2) that is
F(L) S Cp = {= E R(=) SO A [3(() < R(=) tan())),
(38)
where F(L) (LS.) I1122 f dom(L) \ {0}} denotes the numerical range of L, and we let (m;))= be given by the trapezoid rule
f1 pdo
. =,
my
b e
2
(39)
- - pll
. =j
where 0 = to <1 <... < be is (possibly non-uniform) discretization in time. Then, for every 0,
R [.m-Rm1 s 21,
Rmm|
(1+v2)=²* x2 18c cos (0)
A
ptc F
(40)
where At = maxje(411(t; - 1)-1) is the maximal time-step. te is time-horizon and k = max t/-1)/ http//minneeentttt/ - 1 E [1, +00) the conditioning of discretization that measures the amount of non non-uniformity.
Consequently, for uniform discretization t, jAt, O....,L, it holds
R2, Rm < (2 Ise (OT + pat
(1+/2)=²
At = (t.e(Dt).
(41)
Proof. Given USiSE 1, remark that the functions
[/(=) = ((hm (4tt e c((-=>)a di 1-1-t (e -(m-2)) - (-(-=)))=-1) and gu(2)) = / ((-=)))
are analytic whenever R((H)= 2 R(). Thus, since L is stable, we have that
e-
R2, Riml
Σ (L( + (L)
< $1>52 +
Σ |1 + (L)I
==1
19

---

Laplace Transform Based Low-Complcxity Learning of Continuous Markov Semigroups
Next, using that ||A < for all 2 0. we bound the first and last term as
| - ("' e + 4/2 (1+ +e =
ft
1 -- E pul
1 e
=
t S 21
2
and
11 (L)l - ( dt
Eule
e- dt =
f
Finally, we bound (T) using Crouzeix's theorem A.7. To that end, we need to replace the unbounded operator T. with its bounded Yoshida approximation L. := s²(s -) > O, for which we know that for every E dom(L) L.f + Lf. as + oc, (Kato, 2012). This implies that for every > 0, there exists large enough s > 0 so that F(L,) S Co + E. Now, since fu+s is analytic over Cp + e, applying the Crouzeix's theorem A.7 we conclude that
S
1+y2
sup =ECIF(I..)
(f==(())
S
L+v2 2
sup [f(()))
=EC,
So, since S/. is the error of approximating the integral of complex valued function 1. + (x(=22 on the real interval [i (j + 1)Dd] by the trapezoid rule with two points, we have that (see for instance (Quarteroni et al., 2010), equation (9.12))
[1(=) = (1i+1 - 15) ogle- (M- (1, - + E(4+1)) = (1it1 - 14,3 (n -02--10-=) -10=))00(--0)5555)
6
6
for some E (0, 1). Hence, for z E Cju
[f(x) < (1i+1- - 1,)3 11 - 21²e (M-)))(())))
(ts+1 - -t,)3
(1 tan 0) - R(2)) ( (M-
<
6
(1 + tan 2 A) 4e-
2
(t/+1 - -t,)3
2k At
6
e
cos 0
e cos j2
Therefore,
Σ [[A S At
(1+ 22)
1
(1+ At.
I
32a cos2
j2
18c2
cos2 0
j=1
and the proof follows by letting s + oo and noting that sequence of bounded operators (L,)), converges to a bounded operator (L).
Now, appying the above result on the integration bias in (30) for the case without the rank reduction (reasonable when ni,, we obtain
|| (s, S S. S VCN (2+
(1+v2)m 18c2 cos+
At + VCN
(42)
P
which becomes arbitrarily small for At 0 and tt + 20.
On the other hand for the rank reduction case.
(T - Gr Gmar " S.B. RR, Sq P.) & (R,. - S + (P,. (P.))
(43)
that is, |=(( juy - -ā - S VOH 2
(1+/2)=' 1862 (o) At+ VOH e- + -r.- F.
which can be controlled by bounding P, - P. via B=B - D using Theorem A.3. Since in the analysis of the variance term similar construction should
be performed for the orthogonal projector onto the leading singular subspace of the empirical operator å = O, the following we jointly bound last to terms in (30).
1/2 Am in
20

---

Laplace Transform Based Low-Complcxity Learning of Continuous Markov Semigroups
C.3. Estimator's variance
In this section, with the exception of App. C.7, we assume the single trajectory setting with uniform time-discretization, and analyze the last term in (30) developing concentration inequalities in Hilbert spaces for non-iid variables based the notion of beta- beta-mixing and the method of blocks introduced in (Yu, 1994). Before we begin, let us summarize the reminder of the terms, assuming the case of universal kernel together with (BK), (SD) and (RC)
E(Gm,)) S (ca/H) + (B) + S. p.77 Taann-Gii (
(44)
C.3.1. REMINDER ON ERGODIC AND EXPONENTIALLY MIXING MARKOV PROCESSES
We recall some fundamental results ergodic, exponentially mixing Markov Processes.
Definition C.5 (Strict stationarity). A Markov process X = (Xi))ex with values in X' is strictly stationary if for every m, I E N the marginal distribution of (X........ Xm++) is the same as (Xi,..., Xm).
For set I " N and a strictly stationary process X = (Xi)sen we let E, for the o-algebra generated by (Xi)eer and μrr for the joint distribution of {X. Notice that pi+s = paz. In this notation = M(1) and Pr = ((1,++))- We can now introduce the 3-mixing coefficients
3x
(T)
sup ((LL4T) (B) - P(1) P(L+)) (B) BEΣςY
which by the Markov property is equivalent to
8x (T) =
sup BEE &
| (B) p1 x pJ (B)),
where J C IN with j > i for all i E I and j E J. The latter is the definition of the mixing coefficients for general strictly stationary processes.
is well-known that a stationary, geometrically ergodic Markov process is 3-mixing with 8x (r) - 0 exponentially fast as 7 co, that is 3x (F) S me 77, for some 7), 7 E (0, 00) for all E N. See e.g. (Bradley, 2007, vol. 2 Theorem 21.19 pp 325).
The following result exploits a block-process argument to extend concentration bounds from the iid setting to strictly stationary 3-mixing Markov processes.
Lemma C.6 (cf. (Kostic et al., 2022), Lemma 1). Let X be strictly stationary with values in a normed space II-I) and let F, m E N such that m is the largest integer satisfyong n 2 2m7. Moreover, let Zi..... Zen+1 be m. + independent copies of Z- = Ei X4. Thenfor s > O.
R'
P
Σ
X
8) S P
+
23
Σ x
51 2
+ 2m3x (F
(F-1). 1).
==1
where we define &s E |0,27 27,27-13 1] such that k n 22m +F with l {0, 13 and k' E [0, 10,5-11. and X = X(2m+1)++1
for all i E [1. KT
Proof. Recall first the definition of the blocked variables
(2/-1)7
2jF
Y, =
Σ
X₁ and Y
Σ
X
(45)
i= 2 1 +)
i- 12
Then, we have,
k
Fr
X
Σ
X X2mp+i +i
1 1=1
Y2 +
X X2mr+i
+
Y
X2m4+11+4
s=1
t=1
1=1
21

---

Laplace Transform Based Low-Complcxity Learning of Continuous Markov Semigroups
where we let k n 2mF lF + k'' with / {0,1] and k' [0, 1]. and use the convention Σ 0. Thus, for 0,
n Σ x
m+1
m
P
>
Σ
23
> $1 21
(m - 1+ 1+1)) -
Σ
Yj
+
x
2
8
+ +m3x(F 1)
i=1
j=1
j=1
1=1
m+
m
$
Σ 23
>
2 +
$1
z3 +
Σ x > 2 + 2m,3x (F - 1),
2
j=1
j=1
i=1
where we have defined X = X (2m++1)++ for all E [1, and we have used (Kostic al., 2022), Lemma to get the first inequality.
C.3.2. CONCENTRATION FOR TRANSFER OPERATORS
We consider the transfer operators Tja St Sy for any j to...., O and their empirical versions
Tz
o ((Xis)) @ (X((n)))
We prove several concentration results on operators C' and Tj j E [1] and Hm
Proposition C.7. Assume that n - j 2 2mnF Let d > (m 1)8x (F - 1). Let Assumption (KE) be satisfied. With probability at least - 8 in the draw Xo ~ n,, Xisi p(X ((X -), i [n].
1,, Tsae) $8 /22 In 8 2(m - 1)3x (7 1) (1
G
(46)
B)
m²nr
1/2 Proof of Proposition (C.7). Set E(x) = C 1/ ((x).. Assume first for simplicity that n -j j 22m.. We have
-T):: =
{(Xaar) @ 6(X ((((+))))) - E (<(Xise) @ 6 ((((+)))))))
Set A, == ((X(i-1) @ o( ((((1-+)) at). Let Zi...., Um m iid copies of Z = E AG. For any 2 2, by
triangular inequality and convexity of a + xh, we have that || < (+ ||Ailln)) next that
k < E | Observe
E|I|A = E (||((X
E [[
199
(x-12 Vcr cx) k-2 (Vex tr(
2
= || k-2 |
tr(c, C) S 2 k!
In view of (Fischer & Steinwart, 2020, Lemma 11), we have under Condition (SD) when 3 < that
(C) <
8 Vy > 0.
3
Consequently
B
1 k! S
(7-1/2
k-2
(Vox
E||| 7 His
E I1
VC+CH)
--
3
(47)
1 B
22

---

Laplace Transform Based Low-Complcxity Learning of Continuous Markoy Semigroups
We apply (Kostic al., 2023, Proposition 9) get with probability at least - 8,
2 8
Cg + GT V (
Z
< 24v2cH In
3) m²yy
1=1
Replacing 8 by - (m 11 ae (7 1) and union bound combining the last display with Lemma C.6 gives the result.
We assume now that the Markov process X. (Xisr)iEN is ergodic, exponentially mixing, that is, there exists Cuix E (0, oo), such that for all T E N,
x (7) Cmix Arwar VF 2 1,
(48)
where we recall that w. 1o,=-)((4+))/2200
For any n 2 1 and 8 E (0, 1), define the rate
In2 (g) At W. n
In (o)
(6) = c
+
(49)
v W.w.
where c c(Cmixx is large enough numerical constant.
Proposition C.8. Let 8 E (0, 1). Assume that (n j) At u+ 2 4ln 22 (m
Let Condition 48 be satisfied. With
probability at least 5 in the draw Xo ~ n, Xist ~ p( ((X((-1)Ar))) -), [n].
Tyst Tase| S En-)(o).
(50)
Proof. We start from the following result (Kostic et al., 2022, Proposition 3) assuming n - 2 2mn
P [ - Tsae S m
4m
||C||
4m
In (ā (m - 1x (7 17) 12
2 1 d.
m δ (m 1) 3x ai (7 1 1)
(51)
For any j E [1], we take integers my. Fj 2 such that n - j 2 2y and 8 2 2 2m, - 1) 3x. (7) 1). Hence we can pick
(n - () At w.
1
2e Ceeemim
m, : 2
In and ;= = Al A.
(52)
2e Cmiy (1-)) (1-
Assuming that is large enough such that (n j) At 2 4ln
2e (m)
we get
[e (ee (in--)
mj ス + (n - j)) Attw./ In
We also have in view of (48) that (my - 1)3x ae (F - 1) s 5/2. Replacing these quantities in (51), get the result.
Define for any n 2 1
Cr In(no-1)
Cr In(nd-1) At 7r
(7,
C (1,8) +
C (7, (1,
(53)
ntt u4 70
and
(1:,6)==11 ( In
)
tr(C tr(c,,0)
||
where the numerical constant e only depends on Cunxx-
23

---

Laplace Transform Based Low-Compicxity Learning of Continuous Markov Semigroups
Proposition C.9. Let Assumption (KE) and Condition (48) be satisfied. Then, with probability at least 8 in the draw Xa n, Xist (( [n].
(0, (7,6))
In addition, for (0,1) and n large enough such that e/h (7: 8) E (0,
c, 1 - 1 (1,
1
Finally, for any j E [1] such that n - 2 2mF, we have with probability at least in the draw Xo ~ n, Xise ~ ((X((-1))t:)) ) [n].
|, T2 1/2
- (1,0)
Proof of Proposition C.9. Define for any m 2 1 and (0,1)
e (7, 8) == 4e, L' (7, +
2cr (7,5),
(54)
3my²
mm
where
tr(C++)) c
C (7, = In ( + In
Combining Lemma C.6 with (Kostic al., 2023, Proposition 13), obtain with probability at least 1 6,
| (C (0), // < 1 (7. 8/. (m 8x. (F 1)).
(55)
Exploiting again the 3-mixing assumption in (48), we take F as the smallest integer such that
4cmmxx n1
121+
Atw. In
(56)
With this choice of and picking m as the largest integer such that n 2 2m7, then we get
(0(2)) (7, 8/2 (m 1) 3x. (7-1)) s 2=!! (7, 8/4) s (7,8),
provided that the numerical constant e c(Cmix) > large enough.
Finally, for d E (0, 1) and large enough such that E(l' (7, 6) E (0, and observing that
1 - - 11,,,,,,,,,1 1/2 = c --
1 -
get
02/2 0 14/2
1,,
= 1/220,1000- 1/12
1
| -,
1/2 C,CC 1/2
(1) En (7,5)
(,,1/2 TT11) follows from the argument.
The proof for the last result on
Proposition C.10. Let the assumptions of Proposition C.7 and Condition (48) be satisfied. Assume in addition that At E (0, 1) and w, At < 1. Then with probability at least 1 8 in the draw Xo n,, Xist ~ -(X(((x -), i i [n],
2(2) (7,
c. 1/2 Am)) S
(57)
1
where
In2 n((++))
n((++) 3/
(2) (7, =
(58)
Al A
n
At w
8 n - 1)
where e c(Cmixx CHH, Cr, 03, 8) > O is a large enough numerical constant.
24

---

Laplace Transform Based Low-Compiexity Learning of Continuous Markoy Semigroups
Proof of Proposition C.10. For any j []] we take integers my: Ty 2 such that n. j 2 2m and 6 2
4 (m-- 1) 9x..s (7g 1). We the same choices as in (52).
Replacing these quantities in (46), we get with probability at least 1 5,
| |/ Tsan)
S 88 8222 In
8
(V
26
In
22 -j)
3)7 (n j) AlW. W4
v 22,
In
2e Cmix (n2 - ( ))
(2) [n)
(59)
(n j) At w. 77/2
8
Next, by definition of Hm, Hm and elementary union bound, get with probability at least 1 - 8,
| 1/2
Hm)
|
S
At
Σ
1,,1/((TAar
- Tiaii/
W/
1=0
At
W/
(12) (5/(1 1
3/2 m((+1)
e PAJA Aa
In²
1((+1)
u Ar
e
C
Vtt
v /w- //
vn
w.yT/2
==0
j=0
for some large enough numerical constant C = C(CH, Cry Cmix, Cnixxca,,?
that can depend only CHt Cr Crnix; CB. B.
We apply the well-known Abel transform formula
j fi) Σ gk:
Σ fagh fi Ei 9k
(fs+)
5=0
to series E/zo
cuser (n-f)a+ with E {1/2, 13. Since we assumed that μ E (0, 1), elementary computations give
Σ j=0
1
1
1
n j)a
(n 1)a HI
(n - j) (n jn-j-11 1)a
)=
1
2
(n 1)a 1 - e Håt
HAt (n - 1)a
where we have used the inequality -= < - x + 222 true for any x E (0, in the last line.
Proposition C.11. Let Assumption (KE) and Condition (48) be satisfied. Fix E (0, 1) and assume in addition that
n is large enough such that n At w. 2 4ln 2x ggint). Then, with probability at least 1 - 5 in the draw Xo n,
Xist p(X((-1)A0))) i [n].
- 1,,12(0-
(60)
Proof of Proposition C.I.. Similarly the proof of Proposition C.7, we obtain, assuming that n > 2mr, with probability at least 1 - 8 in the draw Xo n. Xisi (i-(i=)) +), [n].
1, - & 8 /
Cr
(61)
8 - 2(m - 1)3xx ai (7 1)) (1 8 m²ag
We pick integers m. F 2 1 such that > 2m and 8 2 2 (m 1) 3x.sr (F - 1). Hence we pick
n At W.
1
22 22e"Cmix (m
m
and =:=
At In
(62)
2n
20
25

---

Laplace Transform Based Low-Compicxity Learning of Continuous Markov Semigroups
Assuming that is large enough such that n At w, 2 4n (==)) we get m 2 f At w./ In ( qux=))))) We also have
in view of (48) that (m 1))8x (7 - 1) - 8/2. Replacing these quantities in (61), we get the result.
C.3.3. VARIANCE OF SINGULAR VALUES
We introduce the approximated and empirical KRR models as Gm,y C. Hm and Gm. - - Hm. respectively, 1/2
We recall that Ciny C. 1,,,"1,,111 HmIr mml and Gin, C 12" HmIr. We also introduce B == C. C, H₂ and 1/2 B 0/1/2 Hm. let denote P. and P. denote the orthogonal projector onto the subspace of leading r right singular vectors of B and B, respectively. Then we have [B], = BP, and [B], BP,, and, hence Gmn,,y Gm,,P, and G Gm..Pr...
Let us recall that (,, Hp
1] and E [3, 1]. Let B = C, Proposition C.12. Let (RC), (SD) and (KE) hold for some a E [1,2], 3 E (0,
1/2 H
and B- == C-1/2 Hm- Given > if 2 En
(7,3 < 1/2, then with probability at least in the draw Xo ~ T.
Xiat ((X(i-1) ), E [n].
| BB-B-BBBI
2 (Vox + 3 2 + Von (o (= (7, (7,8/3) + ven μ s (At)
(63)
122
where (,, = 1,, 8/3) (2) (7, 8/3)ca 11/2 ( (7, 8) is defined in (58), and sn (At) in (41). Conse-
( (7, 6/3)/ | + VCH Sp S vox it holds that quently, for every i E [n]. when en
[ (B - 2/CH (1: (/3) + von 1 Sp
(64)
Proof. We start from the Weyl's inequalities for the square of singular values
lo²( - - B B i i [x)
But, since,
B- - -BG--BB=1l H C Hm - (im - H C. + H HH,) Hz[O - cr
denoting M = -- ( Hµ), N 1 - C) and R= -C-/² (G Gurd) we have
M + +M' +M' -2/221111"0 B N) = Gwn) B = M C2 C1 Gm,, + B" M- B" NG = B
= B + B 13" NG (M* -M--"" NO, 1
= =GGnn) (Hm - H₂) + (fm - Hn
M.y - (Gr,a) ) -CG + +M* + (Gnn)
Note next by definition of Gj,yy and G Gm,y, we have
Gmn GH.TT C 1/2 1/2 H₂) - 1C "C f/
(65)
Therefore, due to (65), R = 0 12 - C 1/2 (M - NGpa), we conclude
B B²B = (Gpa)" (H - -Hp) + (A - Hp)'Ga (Gpn)"
+ (M "O/=0"0000000 -NWm))
(66)
Hence,
1B. BBBBBB BBB
+ | + || [| |NT
$ 2 B|| + C / ([|MII + |I + || |[
Now, noting that
|M < ( Hm)) + // - Hnي < || - -Hm)| VON Rm - Ryy
and applying Propositions C.4, C.9, C.10 and 1.11 we obtain (63), and, therefore, (64) follows.
26

---

Laplace Transform Based Low-Compicxity Learning of Continuous Markov Semigroups
We remark that to bound singular values we can rely on the fact
lo,(B)- [3 - (o:0))-a000--00000000000
اا (B + oB V (
(67)
C.3.4. VARIANCE OF RRR ESTIMATOR
Proposition C.13. Let (RC), (SD) and (KE) hold some a E [1, 3 E (0, 1] and E [8, 1]. Given > 0 and y > if (1)) (7,5) < 1/2, then with probability at least 1 - 8 in the draw Xo ~ n, Xiar p( ) i [n].
2 (7,3))
cH/H (voul (1) (1, 5/3) + sn At) (ca)) /12
Sa Sa(Gm, -(u))
-
- 1
E.
+ VCH us (Ste(
(68)
CH 7al2
a/2
St | - Gnn) $ (C/2(C/nn - -Gp and
Proof. Start by observing that
C2 GGpn) G -CC
(=-1 + ( - Hp)] P, + B(P, P.).
(69)
Taking the norm, using that the norm of orthogonal projector P, is bounded by one and applying Proposition A.3, we obtain
B-
(mm | 'c |0 1| +||05 (8m Hい Hm)|]| + |B||
(
To complete the proof it suffices to apply Propositions C.2, C.3, C.9, 1. C.11 A.3 and C.12.
We remark that the previous proof simplifies when we do not have the rank reduction, that is we can estimate
1x (,m, Gnri ( s c c
||| ||| $ (C (C-c) || +1.11++/// |||c5 (A7m H Hm.))
27 m (r. 8/3)
(3)
+ 2// sp.e (At).
(70)
and further obtain that
1a-1)/2
5/3) 22,
1 (6mal s ||Guzill + Vi ( -Guri
Ca CH 1
+ 22/0 sne(At))
(71)
VT
1
C.4. LaRRR operator norm error bounds
Theorem 6.2. Let L be sectorial operator such that w, XL Let (BK), (RC) and (SD) hold for some a and (0, 1). respectively, and cl(mm()))== Given 6' (0, 1) and re[n], let
In (n/o))
atj
In (n/0)
(0)
(21)
AL w.
W.
At=s; and 1/E=o(E")) then there exists constant c> O, depending only H and or that for large enough ni2r with probability at least 1 - 6 in the draw of D it holds that
0, such
E(Gm.) S max (ār+1, or+ I (6).
(22)
Proof of Theorem 6.2. Note first that under the gap condition o, (R )a) - (Ry > 0 and by choice of in (73), we have for n large enough that dZ( - (Zp) - (ca/1)² cH
a/2 20/2 > -1 0.
27

---

Laplace Transform Based Low-Compicxity Learning of Continuous Markoy Semigroups
Next we set n/22. By definition of the rates en (6), En (7, 6), ( (7/ 6) and / (7, 8), and since 7 E (8, 1}, we have for n large enough, =/"h 6/3) E (0, 1/2), npAi >> and
ln
33) m. (7,3 VOHU s. (At)
At 709
We note indeed that / = n1/2 satisfies the condition 1/l o(E for large enough. Hence, recalling the decomposition of the operator norm error in (44) and Proposition C.13, we get with probability at least 1 - 6,
E(Gmm, ( Or+1
- (22/22 1
+1I).
(72)
At W. 7
By balancing (A) with respect to y first and then (B) with respect to At, we derive the upper bound on (B) S (6) for the choices of y and At given in (73). Finally, if 0 (B) = then we proved the result. Otherwise (B) > O, then an union bound combining (72) with (64) and (67) gives the final bound.
C.5. LaRRR spectral operator bounds
First, recalling Proposition 6.1, essentially proven in (Kostic et al., 2022; 2023), with the only difference that instead of applying it to TO, we apply it to the resolvent of IG. This exposes the perturbation level such that the pseudospectrum of R, contains eigenvalue of its estimator. The bound the perturbation level comes in two forms: one specific the targeted eigenvalue and one uniform over all leading eigenvalues. Since we have resolved the bound on the operator norm error, the only term that remains is the metric distortion for which we have the following
|
(h)) - TF GS
and
T n & (g((h) A 7.) n(hp)) ni 1 dcl
Now, applying this to the LaRRR estimator, since we have that lo, lo, Gr or (S Gpm)) S 1| "S=(,mn, Gb,))]l and ž or (B) v ||
ar CT, or or(Ck/ -5
using Propositions C.3, C.13 together with (71) we obtain that for every [r] metric distortion (in the worst case) bounded by
(ha) S (ca/H) ()) (a-1)/ 1)/2 + -/-1/2
CH
(6)
Sw)
07 (R
which further ensures that r/hh))= - n S (8). Hence, checking that
(1+ ata
7- -11= = (In (n/(()
In In²(n (n/5)
In 3 (n/6)
In (n/6)
At w.
/
2(a+3)
n2
W.
/ 71
2ia
that is
2n 2ug 21a 8[[3+422]
2u
In IIn (n/0)
1/2 (2)) x
11
1
W.
remains bounded w.r.t. n + oo for 2 1, we can apply Theorem 6.2 and Proposition A.2 to obtain the statement of Theorem 6.3. Alternatively, we could apply Proposition A.5 to obtain that for general (non-normal) sectorial operator the estimated eigenvalues satisfy
A
Ad Ad |u Ao
(141(1 +5 (8), or (RPSS)
S ||FI||FF-1|
(ā-1
|
where the bounded operator with bounded inverse such that F-1FF diagonal operator.
28

---

Laplace Transform Based Low-Complcxity Learning of Continuous Markov Semigroups
C.6. LaRRR bounds in the misspecified setting
In this section, we briefly discuss the case where the eigenfunctions of the generator do not belong to RKHS, i.e. the misspecified learning setting when a < 1. First, we note that the misspecified setting for the learning of transfer operators was considered in (Li et al., 2022) in terms of the excess risk Wiilbert-Schmine norm) of kernel ridge regression (KRR), while the error (operator norm) of KRR, principal component regression (PCR), and reduced-rank regression (RRR) was addressed in (Kostic et al., 2023). In both works, risk/error could be made arbitrarily small for every E (0, 1), however at arbitrarily slow rate as a + 0. By adapting our proofs, notably in (68) by avoiding the splitting of Gp,7 from the concentration result in Proposition C. CO, one can also derives of the operator norm error for learning of Ry.
On the other hand, recalling Proposition 6.1, the main challenge of the misspecified setting arises when developing spectral learning rates. Namely, since given any eigenfunction E 2 (X') [H], for every sequence (Ne) S HH of functions in H that approximates it, i.e., such that || S, 2 O, metric distortions explode, that is 17(hs) + oc. Consequently, with current approaches, it is impossible to derive learning bounds for an eigenvalue whose corresponding right eigenfunction lies outside of the RKHS.
C.7. LaRRR bounds for multiple unevenly sampled trajectories
Here we only derive a bound on the operator norm excess risk when we have access to multiple iid trajectories sampled unevenly, noting that spectral bounds then readily follow in the same fashon as above. Here n stands for the number of trajectories and E is the number of samples on each trajectory.
Theorem C.14. Let L be sectorial operator such that w. -02 (2)( Let (BK), (RC) and (SD) hold some |E[1, 2] and (0, 1]. respectively, and (== Given 6E 1) and re[n], les
In2 (E8
a+3
(nn (nn
I(S+0]
(0)=
(73)
W.
n2 W.
then there exists a constant e >> 0, depending only on H. and gap o (Ru() (R > 0, such for large enough n2r and 2 1 with probability at least 1 - in the draw of Dr it holds that
E(Gm.)) S T A ( + (
(E,)1)) + +12
(74)
where At = max t)-1) is the maximal time-step, te is time-horizon and conditioning number k maxje(++i[(tz - 1// minje(++1](t - 1)) E (o,+ quantifies the unevenness of the sampling.
Proof. We proceed with the same decomposition (30) of the excess as in the proof of Theorem 6.2. The steps (0), (I) and (II) remain unchanged. The integration bias in (III) is again tackled with Proposition C.4. As for the estimator variance in step (IV), we now use concentration result for multiple iid trajectories instead of the results we derived in App. C.3 for a single 3-mixing trajectory. Specifically, Proposition C.7 is replaced by Proposition 14 in (Kostic et al., 2023). Proposition 11 is replaced by Proposition 13 in (Kostic et al., 2023) with n replaced by nıl 1). Proposition 14 is replaced by Proposition 17 in (Kostic et al., 2023). Hence instead of Proposition 15, we obtain the following control on the variance term. For n and
é large enough, in view of (40), we have with probability at least 1 1-6 6,
In (1) W. 708
Sx (a - Gr Gpai
+ + A -
(75)
Thus we get w.p.a.l. 1 8,
E( E(Gmm)) or+1
S 1( /0/2
1
In (1)
+ x²2 At -
(76)
1
v /w. 70
(A)
We can balance (A) as in the proof of Theorem 6.2.
Next set At- (in 1)-1). Note we have the following conservative bound for I large enough
PAU At
S At.
PAt
The end of the proof follows from the same argument as in the proof of Theorem 6.2.
29

---

Laplace Transform Based Low-Complexity Learning of Continuous Markov Semigroups
D. Experiments
Comparison
estimating
0.0
1.75
1.5
relaive -2///h)
1.50
1,0
0.5
As
A
075
-2.25
-2.00
-2.75
-3.50 -1.23 -1.00
-0.15
+0.50
Figure 4. Comparison of the primal LaRRR algorithm the TO RRR method for estimating the slowest timescales of the process in the basis of random Fourier features (Rahimi Recht, 2007). predicted by our theory, LaRRR (blue) remains stable as - 0, whereas the error of TO method (red) diverges. The main plot shows the three quartiles of relative error for λ2 across 10 independent trajectories of a ID Langevin process on a triple-well potential. The inset displays the distributions of the top three eigenvalues for Åt 0.005.
25
10
5-
45
-30
-55
20
-15
1.0 -5
45 40 35
30
-25
-20
-35
100
Figure setting experiment on Overdamped ID Langevin dynamics, we show how LaRRR (presented red) compares to two physics-informed IG baselines: energy based regression of Dirichlet form by Kostic al. (2004) (magenta), and Galerkin projection by Hou et al. (2023) (blue). The eigenvalues are shown on the horizontal axis, and the vertical axis indicates the histogram hight of those estimated by Galerkin projection. While both the method of Kostic al. (2004) (left) and LaRRR (right) recover the eigenvalues of IG (black lines), Galerkin projection (due to the unbounded nature of IG in 2 produces many spurious eigenvalues near the origin.
30


Reviews:
[{"review_id": "USH5WA6ByP", "reviewer": "Reviewer_PJzG", "summary": "This paper studies a new class of non-parametric learning algorithms for continuous-time Markov processes, specifically for learning the eigenfunctions and eigenvalues of their infinitesimal generator (IG) of the semigroup of transfer operators (TO). While existing methods tend to focus on learning the TO which share the same eigenfunctions to IG, the authors criticize their deteriorating spectral gap as the sampling frequency of the data increases (i.e. as the time-lag decreases). Then for more recent methods that directly learn the IG and its eigenstructure, the authors point out that the challenges coming from the unboundedness of the IG is not well addressed. Thus, the goal of the paper is to propose a learning algorithm for IG that can properly handle its unboundedness. To this end, the authors propose to leverage an auxiliary operator, called the resolvent operator, which has the same eigenfunctions as the IG, and can be obtained through the Laplace transform of the TO. For a class of Markov processes with sectorial IG (this includes all time-reversal-invariant Markov processes and some important non-time-reversal processes), the resolvent is uniformly bounded outside a sector containing the spectrum, which addresses the unboundedness problem. The authors then extend an established prior method (reduced rank regression; RRR) for operator learning in RKHS to learn the resolvent operator, and provide statistical learning guarantees. A particularly new contribution is bounding the integration bias, i.e., the error of Laplace transform under possibly irregular sampling intervals, using spectral perturbation theory. The authors demonstrate the proposed algorithm in two time-reversal-invariant processes. The first experiment demonstrates that while TO learning is sensitive to small time-lag, the proposed method learning the resolvent is robust. The second demonstrates that the proposed method successfully recovers the leading two eigenfunctions of IGs from a molecular dynamics data, in a setup where a previous IG learning method is intractable due to the high dimensionality.\n\n## update after rebuttal\n\nThe authors provided detailed clarifications in response to my concerns. I think this is a solid work in the domain of nonparametric learning of dynamical systems and would like to retain my supportive rating.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer’s insightful evaluation and valuable comments. Below, due to space limitations, we briefly address the highlighted weaknesses and respond to the reviewer’s questions, committing to incorporating all feedback in our revision. If needed, we can elaborate further on each point.\n\n## Additional empirical evidence: \nFollowing the reviewer's suggestion, we have expanded our empirical evaluation to include IG baselines and support the claims summarized in Table 1. For a brief discussion on the [results](https://green-chantal-92.tiiny.site) of our additional experiments, please see our reply to reviewer __nvbC__ . \n\n## Algorithms\nWe will expand Appendix B with detailed derivations, including a discussion of the sampling operator $\\hat{S}$. Briefly, by the definition of $\\hat{S}$ for $h \\in \\mathcal{H}$, we derive $\\hat{S}^* \\hat{S} h = n^{-1/2} \\hat{S}^* [h(x_1) \\dots h(x_n)]^{\\top} = n^{-1} \\sum_{i \\in [n]} \\phi(x_i) h(x_i) = [n^{-1} \\sum_{i \\in [n]} \\phi(x_i) \\otimes \\phi(x_i)] h$, concluding that the empirical covariance can be written as $\\hat{S}^* \\hat{S}$.\n\n## Transient dynamics\nAs discussed in Appendix A.3, when $L^* L \\neq LL^*$, eigenvalues alone may not fully explain the evolution of linear dynamics. The norm of the resolvent relates to transient growth in stable systems via the pseudo-spectrum (Trefethen & Embree, 2020). We will expand this discussion in the revision.\n\n## TO as $\\Delta t \\to 0$\nThe core idea is that $A_{\\Delta t} = e^{\\Delta t L }\\to I$, which implies that the spectral gap vanishes, thereby rendering the eigenfunction bounds vacuous. \n\n## Minor issues\nWe will correct all typos. Additionally, we remark:\nThe presence of $G^*$ in the risk can be understood from two perspectives. One follows from TO’s risk formulation in the space of observables, which, due to the kernel trick, translates into vector-valued regression (Kostic et al., 2022). The other, as the reviewer notes, considers the Perron-Frobenius operator, (Klus et al., 2019), on probability distributions, linking the estimator to the MMD metric induced by a characteristic kernel. We will expand Section A.2 to introduce regression-based operator learning for both TO and IG methods.\n$\\rho_t$ denotes the distribution of $(X_s, X_{s+t})$ and $\\mathrm{gap}_i$ is introduced in the last line of Theorem 6.3. A table of major notations will be added at the beginning of the appendix.\n\n## Questions\n[Q1] When $X_0$’s distribution is not invariant, variance analysis requires adapting the method of blocks for mixing with Bernstein inequalities in Hilbert spaces for independent (but non-identically distributed) variables. This complicates the effective dimension, as the covariance operator w.r.t. the invariant distribution is replaced by one w.r.t. the ergodic mean. While technically feasible, this adds complexity, so we opted for a simpler approach to highlight key contributions.\n\n[Q2]  Knowing IG implies knowledge of the SDE, requiring a numerical solution for process realizations. In contrast, given IG’s spectral decomposition, Equation (8) provides directly the solution, since it enables forecasting of full state distributions beyond just the mean (e.g., $f$ as an indicator function). See reply to __CrN8__.\n\n[Q3] To our knowledge, our technique for controlling Bochner integral approximation errors for linear operators is novel, possibly extending beyond TO/IG methods. We will emphasize this in the contributions paragraph.\n\n[Q4] As the reviewer rightly points out, the Laplace transform is a well-known analytical tool in the study of continuous operator semigroups. In the context of deterministic dynamical systems, it was used by Mohr & Mezic, (2014) to investigate the spectral decomposition of the Koopman operator (TO for deterministic dynamical systems). Their work underscores its potential but leaves efficient numerical methods as an open problem. In contrast, we consider SDEs and not only design numerical methods based on the Laplace transform, but also develop statistical learning theory, providing sharp bounds on spectral estimation. Additionally, we remark that the results of (Bevanda et al. 2023) are also limited to deterministic dynamical systems. There, the Laplace transform served more as inspiration than an explicit component of the method, as the authors built an RKHS via finite-horizon integration of kernel features over trajectories to formulate Koopman operator regression and study learning bounds.", "claims_and_evidence": "Most of the claims are supported by theoretical or empirical evidence. Some claims listed below could benefit from additional clarification.\n\n- In pages 2, 3, and Table 1, it is claimed that current methods that directly learn the IG are susceptible to spurious eigenvalues due to the unboundedness of the IG. This claim is not explicitly supported by theoretical results or experiments as far as I can confirm.\n- As a related note to the above, while a proper empirical comparison against TO learning method is made in the first experiment, comparison against IG learning methods are only done theoretically, not empirically.\n- In page 3, 8, and Table 1, the authors claim that the proposed method is structure-agnostic and apply to the broad class of sectorial IGs which covers not only self-adjoint IGs but also important non-self-adjoint ones. However the two experiments both concern self-adjoint IGs, as far as I understand.\n- In page 5, the authors claim that the definition of the sampling operator and its adjoint implies that the empirical estimation of the covariance can be expressed using them. It would be better if the derivation is made more explicit.\n- In page 4, the authors claim that while eigenvalues are informative about long-term behavior, they are static properties and fail to capture transient dynamics or the full time evolution of the process. Can the authors elaborate a bit on this? One might argue that eigenvalues decaying rapidly over time can provide information about the transient dynamics.\n- In page 7, the authors claim that the eigenfunction learning bound for the considered TO learning method becomes vacuous as the time-lag goes to 0. Can the authors elaborate a bit more on why? I was not able to fully understand the description.", "methods_and_evaluation_criteria": "The proposed method is suitable for the problem setup considered in the paper. For the benchmark dataset, the scope of the experiment could be improved to include non-time-reversal-invariant Markov processes, which could add a strong supporting evidence to the paper; please see the Claims and Evidence section.", "theoretical_claims": "I carefully checked the soundness of the problem setup, including Markov semigroups, transfer operator and their infinitesimal generator, the resolvent operator, and their spectral decomposition including the compatibility of eigenfunctions and they are sound. I did go over the derivation of the bound of the integration bias but could not carefully verify the correctness due to my limited expertise in spectral perturbation theory. I did not carefully check the proofs for the other bounds.", "experimental_designs_or_analyses": "Please see the Methods and Evaluation Criteria section.", "relation_to_broader_scientific_literature": "Resolvent operators of Markov semigroups and their expression as Laplace transform has been studied theoretically (Engel & Nagel, 1999) but has not been leveraged in machine learning context as far as I know. The closest ideas I am aware of are identifying Koopman eigenfunctions directly using Laplace transform or similar integral operators (Mohr & Mezic, 2014; Bevanda et al., 2023), which I think are related but not exactly the same to the approach in this work.\n\nEngel & Nagel, One-Parameter Semigroups for Linear Evolution Equations (1999)\n\nMohr & Mezic, Construction of Eigenfunctions for Scalar-Type Operators via Laplace Averages with Connections to the Koopman Operator (2014)\n\nBevanda et al. Koopman Kernel Regression (2023)", "essential_references_not_discussed": "Important related works are properly discussed in the paper.", "other_strengths_and_weaknesses": "Strengths\n- Originality: The paper is original in its use of the resolvent operator and its expression as Laplace transform of the Markov semigroup in the non-parametric operator learning context.\n- Significance: The paper is significant in its generality (sectorial IGs, and possibly handling irregular sampling intervals) and addressing of the issues of the current TO and IG learning methods (sensitivity to small time-lag and spuriousness due to the unboundedness of IGs, respectively).\n- Clarity: The paper thoroughly introduces the operator theory backgrounds necessary to understand the theoretical results.\n\nWeaknesses:\n- Clarity: The readability of the paper could be improved, especially in its description of the kernel learning algorithm in Sections 4 and 5. Currently it requires that the reader is already familiar with Kostic et al. (2022), Kostic et al. (2024), and Turri et al. (2023).\n- For other weaknesses, please see the previous sections.\n\nKostic et al. Learning Dynamical Systems via Koopman Operator Regression in Reproducing Kernel Hilbert Spaces (2022)\n\nKostic et al. Learning the Infinitesimal Generator of Stochastic Diffusion Processes (2024)\n\nTurri et al. A Randomized Algorithm to Solve Reduced Rank Operator Regression (2023)", "other_comments_or_suggestions": "Minor typos or ambiguities:\n- Page 1: trajecotry -> trajectory\n- Page 1: In the introduction it is stated that current kernel methods for TO learning require kernel function selection, but this is a shared limitation with the proposed method.\n- Page 2: \"the space of functions on [?] that are square-integrable\"\n- Page 3: Below equation (4), which space is $L^2$? The same goes for the first line of page 13.\n- Page 3: Below equation (5), should \"the spectrum of $G$\" be \"the spectrum of $L$\"?\n- Page 3: In equation (8), the symbol $f$ refers to both the observable and right eigenfunctions of $L$\n- Some parts in the appendix denote the resolvent operator by $L_\\mu$, while the main text uses $R_\\mu$\n- Page 3: In equation (9), the reason why the adjoint $G^*$ is used is too implicit. As far as I understand, it is because technically we are regressing the embedded Perron-Frobenius operator (under the terminology of Klus et al. 2019), not the transfer operator (Koopman operator).\n- Page 3: \"an universal approximation result holds true the minimizers of for the minimizer (9)\"\n- Page 5: Below equation (18), which joint distribution does $\\rho_{j\\Delta t}$ represent?\n- Page 5: \"Ridge Regression solution of the regularized risk , that is without\"\n- Page 7: What does $\\mathrm{gap}_j(R_\\mu)$ precisely represent? I might have missed it but it is not defined in the main text.\n- Page 13: \"the RKHS associated to kernel $k\\mathcal{X}\\times\\mathcal{X}\\to\\mathbb{R}$\"\n- Page 15: \"(?)Theorem 2]Kostic2022\"\n\nKlus et al. Eigendecompositions of Transfer Operators in Reproducing Kernel Hilbert Spaces (2019)", "questions_for_authors": "- In page 5, the authors assume stationarity $X_0\\sim\\pi$ to simplify the analysis. How restrictive is this, and what can be done if we were to remove this assumption?\n- In page 1, the authors mention that understanding long-term behavior is essential for accurate forecasting and interpretation, and in page 4 that learning the IG alone is insufficient for forecasting the process and estimating the spectral decomposition is of greater interest. Indeed, the experiments focus on learning the eigenstructure. Can the authors comment on what is additionally needed to use the proposed method for forecasting?\n- In Appendix C.2.4., I am not sure if I understood the first paragraph correctly. Am I right in understanding that the control of the integration term is a novel technical contribution of this paper, which was not established or considered in the previous TO and IG learning methods?\n- Out of curiosity, can the authors comment on how the proposed approach is related to the use of Laplace transforms or similar integral operators in Mohr & Mezic, (2014) and Bevanda et al., (2023)?\n\nMohr & Mezic, Construction of Eigenfunctions for Scalar-Type Operators via Laplace Averages with Connections to the Koopman Operator (2014)\n\nBevanda et al. Koopman Kernel Regression (2023)", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "1IrKNWD4nw", "reviewer": "Reviewer_nvbC", "summary": "The paper deals with learning continuous-time Markovian dynamics. While existing methods focus on learning transfer operators, here the authors suggest to learn a spectral decomposition of the semigroup's generator, under some assumptions. This is done by finding a (finite-rank) approximation of the resolvent in an internal RKHS, thus exploiting its appealing properties (e.g. boundness).\nThis approach is applicable to a relatively broad class of processes, and gives raise to accurate and efficient data-driven algorithms.\n\n## update after rebuttal\n\nThe paper and the author's answers are convincing, and my recommendation to accept the paper remains.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer’s insightful evaluation and valuable comments. Below, due to space limitations, we briefly address the highlighted weaknesses and respond to the reviewer’s questions, committing to incorporating all feedback in our revision. If needed, we can elaborate further on each point.\n\n## Additional empirical evidence (results can be found [here](https://green-chantal-92.tiiny.site) ): \n\nWhile we still feel that our work is mainly theoretical, following the reviewer's suggestion, we have expanded our empirical evaluation to include IG baselines and support the claims summarized in Table 1. In particular: \n\n### 1D Langevin dynamics experiment: \nWe extend the original experiment by comparing it to __two IG baselines__ that use prior knowledge of IG:\n-  **Galerkin projection estimate (Hou et al., 2023)**, and \n- **Energy based Dirichlet form regression (Kostic et al., 2024)**. \n\nIn Figure 4 (see the above link) we plot true eigenvalues as black vertical lines, while the estimated eigenvalues are plotted as magenta (our dual method) and red (Dirichlet form regression) dashed lines. Further, since the Galerkin projection estimate has numerous spurious eigenvalues, we plot their empirical distribution in the form of histogram, stressing out that expert knowledge is needed for this estimator to extract good estimates from the spurious ones.  Finally, we report that **our results are comparable to the SOTA** (Kostic et al., 2024) estimator,  **despite not using the explicit IG knowledge**,  **while being one order of magnitude faster to train**.\n\n### Non-normal sectorial IG:\nWe further conduct an experiment using the 2D Ornstein-Uhlenbeck process with non-symmetric drift, estimating the leading nontrivial complex-conjugate eigenvalue pair. We compare our primal method based on random Fourier features with the corresponding TO, noting that **Kostic et al. (2024) is inapplicable** in this setting and that **Galerkin projection estimators (Hou et al., 2023) are inefficient**, since, as observed above, they typically result in over 50 eigenvalues in the zone of interest. In Figure 5 (see the link above) we show comparison for 10 random trials using 1000 random features, two different time discretizations $\\Delta t=0.01$ and $\\Delta t = 0.001$, and corresponding sample sizes $n=10^4$ and $n=10^5$. Note that, to obtain small errors, sample sizes are much higher than in the self-adjoint case (Langevin). This is due to non-normality of the generator, as predicted by our theoretical analysis, see Appendix C.5, where one can see that our method consistently has estimates in the tighter $\\varepsilon$-pseudospectrum of the operator.\n## Presentation\nAs suggested, we will include a table of major notations at the beginning of the appendix.", "claims_and_evidence": "All claims made in the submission seem to be supported by a sound theoretical analysis and exact proofs (albeit naturally it was impossible to really get into details within the short review period). \nThe experimental section, however, does not show a clear case of superiority over existing methods. A comparison is made only on one example and only to one baseline, where even there the results are quite arguable.", "methods_and_evaluation_criteria": "Basically yes, but more numerical experiments are needed. The inset frame in Fig.1 is unclear to me.", "theoretical_claims": "I briefly went over proof, though the details are too deep to cover during the short review period.", "experimental_designs_or_analyses": "See my answer above.", "relation_to_broader_scientific_literature": "Compared to existing literature, this work extend the analysis to a broad class of dynamics, and find solutions to priorly known drawbacks (e.g. short sampling intervals). The key idea of studying the resolvent seems to be new and mind opening.", "essential_references_not_discussed": "The review of related work seems adequate.", "other_strengths_and_weaknesses": "The submission seems original and novel, with strong theoretical contributions. The empirical contribution, however is not established in my opinion.", "other_comments_or_suggestions": "Although overall clarity of the paper is not bad, it is still very easy to get lost with all symbols and notations. It would be recommended to encapsulate notations and definitions into a single table  (in the Appendix, maybe).\n\nIt is crucial to add extra experiments and comparisons with more baselines.", "questions_for_authors": "None", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Q1HHiTvXIz", "reviewer": "Reviewer_CrN8", "summary": "The authors present an approach for learning continuous Markov semigroups. Notably their approach comes with theoretical guarantees at any time-lag. In addition, their approach scales linearly in the state dimension opening the door to apply their methods on high-dimensional problems. Finally they demonstrate their approach on two test problems.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the reviewer’s insightful evaluation and valuable comments. Below, due to space limitations, we briefly address the highlighted weaknesses and respond to the reviewer’s questions, committing to incorporating all feedback in our revision. If needed, we can elaborate further on each point.\n\n## Additional empirical evidence: \nFollowing the reviewer's suggestion, we have expanded our empirical evaluation to include IG baselines and support the claims summarized in Table 1. For a brief discussion on the [results](https://green-chantal-92.tiiny.site) of our additional experiments, please see our reply to reviewer __nvbC__ . \n\n## Presentation\nTo help the reader, we will include a table of major notations at the beginning of the appendix, and expand our appendix with more context on the kernel based operator learning and derivation of the algorithms. Further, if accepted, we will use the extra page to include smoother transitions between challenging concepts. Finally, concerning the discussion of the Alanine-Dipeptide experiment, since there is no ground truth for this system, the only thing we could show is that the result aligns with the state-of-the-art expert knowledge in molecular dynamics, (Wehmeyer & Noe, 2018). In the revised manuscript we will discuss more on the interpretation of this figure.\n\n## Questions\nWe thank the reviewer for their questions, which suggest potentially interesting extensions of our framework. \n\n- Partially observed systems: The major problem in partially observed systems, in the context of TO/IG kernel-based learning, lies in the fact that one might lose the universality of representation, and hence encounter the difficulty to properly predict distributions. One possible way to overcome this is to rely on Takens-type theorems (deterministic and/or stochastic) (Sauer et al. ”Embedology”, 1991; ) Koltai and Kunde, “A Koopman–Takens theorem”, 2024), that essentially guarantee, under certain assumptions, that by augmenting the current measurements with the past ones we decrease the loss of information and improve estimation. \n\n- Fractional Brownian motion: Consider the SDE $dX_t = b(X_t) dt + \\sigma(X_t) dB_{t}^H$ where $B^H =(B\\_{t}^H )\\_{t}$ is a fractional Brownian motion. When $H \\neq 1/2$, the fractional Brownian motion (fBm) is no longer a semimartingale, which means Itô calculus does not apply directly. Although the equation can still be given a meaning—via the Stieltjes integral for $H > 1/2$ or the Skorokhod integral for $H \\in (1/4,1/2)$—the resulting process is no longer Markovian. Our generator-based approach no longer applies in this setting.  However, as you suggest, we could adapt our method to estimate the generator for a Markovian approximation  associated with the autonomous version of the fractional SDE introduced in [1]:  $dX_t = b_\\theta (X_t) dt + \\sigma_\\theta (X_t) dB_{t}^H.$  This process can be approximated using a finite linear combination of Ornstein-Uhlenbeck processes, $\\hat{B}^H(t)$, leading to the so-called Markov-Approximate fractional SDE (MA-fBMSDE):  $dX_t = b_{\\theta} (X_t) dt + \\sigma_{\\theta} (X_t) d\\hat{B}\\_{t}^H,$  where  $d\\hat{B}\\_{t}^H = \\sum_{k=1}^{K} \\omega_k dY_t^k, \\quad \\text{with} \\quad dY_t^k = -\\gamma_k Y_t^k dt + dW_t.$  By applying Proposition 2, the process $X_t$ can be augmented with a finite number of Markov processes $Y_t^k$ (which approximate $B^H$), forming a higher-dimensional state variable:  $Z_t = (X_t, Y_{t}^1, \\dots, Y_{t}^K) \\in \\mathbb{R}^{D(K+1)}.$ The process $Z = (Z\\_{t})\\_{t}$ is Markovian and can be described by an ordinary SDE:  $dZ_t = h_{\\theta} (Z_t) dt + \\Sigma_{\\theta} (Z_t) dW_t.$  The infinitesimal generator is given by:  $$\\mathcal{L} f(z) =\\left( b (x) - \\sigma(x) \\sum_{k=1}^{K} \\omega_k \\gamma_k y^k \\right) \\frac{\\partial f}{\\partial x} + \\sum_{k=1}^{K} (-\\gamma_k y^k) \\frac{\\partial f}{\\partial y^k}+ \\frac{1}{2} \\sum_{i,j} (\\Sigma_{\\theta} \\Sigma_{\\theta}^\\top)_{i,j} \\frac{\\partial^2 f}{\\partial z_i \\partial z_j}.$$ This satisfies our assumptions, allowing us to apply our procedure.\n\n## Forecasting \nGiven IG’s spectral decomposition, Eq (8) provides directly the solution since it enables forecasting of full state distributions beyond just the mean (e.g., $f$ as an indicator function). Hence, we can use our method to forecast $\\mathbb{E}[h(X_t)\\vert X_0 = x] \\approx \\sum_{i\\in[r]} e^{\\hat \\lambda_i} \\langle \\hat{g}\\_{i}, h \\rangle\\_{\\mathcal H} \\hat{f}\\_{i}(x)$, for $h\\in\\mathcal{H}$, noting that $\\hat{g}\\_{i}, h \\rangle\\_{\\mathcal H}$ can be computed on the training set via kernel trick, see (Kostic et al 2022). Further, note that this formula extends to all $\\mathcal{L}^2_\\pi$ functions, at the price of an additional projection error, and that we can, hence, predict evolution of distributions, see e.g. (Klus et al., 2019) and (Kostic at al. \"Consistent long-term forecasting of ergodic dynamical systems\", 2024). If useful, we can include an empirical example in the main body or appendix.", "claims_and_evidence": "- *Novel approach for learning the spectral decomposition of the IG for Markov semigroups*: The proposed approach is novel as far as I'm aware. This is supported by a solid literature review and a detailed description of the approach. I personally found the algorithm blocks to be particularly helpful. \n- *Statistical guarantees for their approach*: The authors prove tight error bounds for their approach which they validate with numerical studies.\n- *Numerical stability for $\\Delta t \\to 0$*: One of the major claims the proposed approach is that it remains stable as $\\Delta t \\to 0$ (something which is not true for TO RRR). They show that this is indeed the case in the numerical study in Figure 1.", "methods_and_evaluation_criteria": "The authors demonstrate their approach on two problems:\n- Estimating the eigen values of the generator for a 1D triple well problem. \n- Estimating the eigen functions of the alanine dipetide in water with a dimension of 45. It wasn't clear to me what was expected for the ground truth of the first and second eigenfunctions. \n\nWhile I think these problems do a decent job of sketching out some of the main ideas and motivations for learning the IG, some more impressive numerical studies would have definitely strengthened the paper.", "theoretical_claims": "- Theorem 6.2 & 6.3: The sketch of the proof seems correct but I haven't gone through all the details carefully. \n- The discussion on the comparison between bounds from [Kostic'23a] and [Kostic'24] was really helpful for understanding the motivation of the present work.", "experimental_designs_or_analyses": "- To reiterate, it wasn't clear to me what to read from Figure 2. The authors mention that the eigen functions have different constant values in the expected metastable states but some estimate for the accuracy of the approach would have been useful.", "relation_to_broader_scientific_literature": "This work relates to the broad problem of learning models of Markov processes from data.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "Strengths:\n- Your introduction and review of the relevant literature was really well-written and easy to follow. I appreciate the effort and care that went into this section and I think it will make your work more broadly accessible to a wider ML audience. \n\nWeaknesses:\n- While I think your approach is a strong contribution some sections of your write-up were challenging to parse. While some of this comes with the territory and the nature of the problem you are tackling, I think adding an appendix on notation / simplifying some of your notation would make your work more accessible. For example, I had to dig quite a bit to understand the meaning of $[\\cdot]_r$ and $\\otimes$.", "other_comments_or_suggestions": "- under Section 3: \"result holds true the minimizers of for the...\"", "questions_for_authors": "- I'm wondering how much of challenge it might be to apply your approach for estimating latent dynamics (i.e. where not all relevant states are directly observed)? \n- As a follow up to this question, might it be possible to apply your approach to estimate the generator for Markov-approximate fractional Brownian motion? [1]\n- Can you discuss how you might apply your approach for forecasting? \n\n[1] Daems, Rembert, et al. \"Variational inference for SDEs driven by fractional noise.\" arXiv preprint arXiv:2310.12975 (2023).", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "COgi4jVZYk", "reviewer": "Reviewer_iVn3", "summary": "The paper presents theoretical developments to learning the infinitesimal generators of markov semigroups using Laplace transform.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful feedback. Due to space constraints, we address the key issues and questions briefly, committing to incorporate all suggestions in our revision. Further details are available upon request.\n\n# Empirical evidence:\nFor our additional [results](https://green-chantal-92.tiiny.site), please see reply to __nvbC__ .\n\n# Presentation\nFollowing your suggestions, we'll emphasize motivations in Sec 3, address unboundedness effects, expand the app. on kernel-based operator learning for TOs and IGs, and add a notation table. If accepted, we'll use the extra page to enhance transitions between challenging concepts.\n\n# Questions\n- In stochastic process theory, TOs form the Markov semigroup associated with the process, acting as direct image functors in measurable spaces. When a measurable func pushes a measure forward, TO describes density transformations. Neural-network learning methods for TO relate to FNO but are distinct concepts, (Nakao and Mezic, 2020).\n\n- TOs are linear (evolution) operators defined on func spaces. We can say that $A_t$ acts as the pushforward of a function $f$ under the transition probability defined by $p_t$, making them crucial for understanding $(X_t)\\_t$ dynamics. On the other hand, IGs, defined via Eq 25, are linear differential operators related to the eqs of motion via backward Kolmogorov eq. $\\partial_t \\mathbb{E}[f(X_t)|X_0=\\cdot ]=L f$ for observable $f$ and forward Kolmogorov eq. (Fokker-Planck) $\\partial_t q_t =L^* q\\_t$ for probability density $q_t$ of $X_t$ (* being adjoint). Indeed, Eq 7 is obtained by applying Itô's lemma to $f(X_t)$ where $X$ is the solution of Eq 6. \n\n- Learning the spectral decomposition of TO/IG is equivalent to learning the solution of an associated SDE, so we can predict the evolution of distributions, see reply to __CrN8__. Working with IG allows us to do this reliably in continuous time.\n\n- Aș explained in Sec 5, the data is the observed trajectory of states. \n\n- Galerkin projection methods (empirical risk minimization) suffer from spectral pollution when applied to unbounded operators, (Kato 2012, Kostic et al. 2024). When the process evolves in a bounded domain $\\mathcal X$, the IG are still typically unbounded, i.e. the IG doesn’t map all funcs $\\mathcal{X}\\to\\mathbb{R}$ (observables of the process) in the $\\mathcal{L}^2_\\pi$ space to other funcs in that space - it’s only defined on a proper subset (domain of IG). \n# Q. for Authors\n\n[1] $(p_t)_t$ is a family of transition probability densities $p_t(x,y)dy=P(X_t \\in dy|X_0=x)$ for all $x,y\\in\\mathcal{X}$. In this sense, $p_t(x,y)dy$ quantifies the probability that the process $(X_t)_t$ is in the infinitesimal set $dy$ at time $t$, given that it started at $x$.\n\n[2-3] See above.\n\n[4-5] While the action of TO can be defined for any measurable func $\\mathcal{X}\\to\\mathbb{R}$, to characterize dynamics via semigroup, one needs the space of such funcs to be invariant under the action of TOs, and a typical choice is $\\mathcal{L}_\\pi^2$. Note that invariance of measure $\\pi$ (both marginals of $P[X_t,X_0]$ are $\\pi$), doesn't mean that TOs (conditional laws $P[X_t|X_0]$ ) are trivial.\n\n[6] In fact, for all $t$, $A_t$ is an operator that turns a func $f$ into a func $A_t f$ s.t. for any $x\\in\\mathcal{X}$, $A_t f(x)=\\mathbb{E}[ f(X_t )|X_0 =x]$. While $L$ is an operator that maps a func $f\\in dom(L)=\\mathcal{W}\\_\\pi^{1,2}\\subset \\mathcal{L}\\_\\pi^2$ to a func $Lf \\in\\mathcal{L}\\_\\pi^2$, meaning that $\\int_{\\mathcal{X}} |L f(x)|^2\\pi(dx)<\\infty$.\n\n[7-8] In the RKHS theory it is standard to assume $\\mathcal{H}\\subset\\mathcal{L}^2_{\\pi}(X)$, $\\pi$ being the probability of data samples (e.g. Gaussian kernel), and the main issue in the learning bounds is the difference between $\\mathcal{H}$ (chosen) and $\\mathcal{L}^2_{\\pi}(X)$ (unknown) norms, e.g. (Steinwart&Christmann,2008).\n\n[9] For conciseness, the defs of $A$, $L$, and $\\pi$ in the Ornstein-Uhlenbeck case are in App A.1. If space permits we can include them in the main text.\n\n[10-11] The funcs $\\phi$ are not unknown but are provided when choosing the RKHS, e.g. $\\phi(x)=e^{-\\|x-\\cdot\\|^2/(2\\sigma^2)}$ for Gaussian kernel.\n\n[12-15] In Eqs 9-10, $X_0$ is not the SDE's initial cond. but any random variable with the invariant measure, see reply to Q1 of __PJzG__ for extra clarity. Note that $\\phi(X_t),\\psi(X_0), G^*\\phi(X_0) \\in \\mathcal{H}$ are funcs, not  scalars.\n\n[14 & 16-17] $\\mathbb{E} [\\psi(X_0)|X_0=\\cdot]:\\mathcal{X}\\to\\mathcal{H}$ is the Riesz representation of lin. operator $R_\\mu$ in RKHS $\\mathcal{H}$, and our learning objective. It is the regression func related to the excess risk, see Prop A.1, where the risk is given via the target feature $\\psi(X_0)$. Since we cannot observe it precisely (due to the integral form) in Eqs 11-12 we show how to approximate it via quadrature. \n\n[18-19] $A_{t_j}$ is TO for time-delay $t_j$ (see Eq 2), $h$ was used as a generic func in RKHS $\\mathcal{H}$ and $v$ as a vector in $\\mathbb{R}^N$.", "claims_and_evidence": "The paper's claims all theoretically substantiated. Empirical evidence backs the claims, but the experiments are of very small scale.", "methods_and_evaluation_criteria": "Ok", "theoretical_claims": "I read the entire paper carefully, and understood some. I see no issues with correctness.", "experimental_designs_or_analyses": "The experiments are small scale, but seem reasonable.", "relation_to_broader_scientific_literature": "The work is positioned sufficiently to its specific domain, but the broader picture is lacking.", "essential_references_not_discussed": "No issues.", "other_strengths_and_weaknesses": "S: The paper seems to present a strong theoretical advancement to operator learning\n\nS: The theoretical presentation is convincing.\n\nS: The empirical results support the contributions.\n\nW: The presentation lacks exposition and the higher-level story gets drowned in all the technicalities. It was difficult to follow what was happening, and why. The technical presentation is dense and convoluted.\n\nW: The results are consistent with the contributions, but they have small scale and are thus anecdotal. The method is only compared to one other method. The results show that one can learn the first two eigenfunctions from SDE dynamics. It is not demonstrated how is this beneficial, or what problem this solves. Why do we need these eigenfunctions? There are no benchmarks or result tables. I feel like the paper is not really realising its own potential, or had no room left for showing how the contributions lead to improved methods and results.", "other_comments_or_suggestions": "I want to start by stating that I’m not a theoretician or expert in the paper’s domain, but have knowledge on many subtopics of this paper. Hopefully my review can still be useful in showing how the manuscript opens up to an interested outsider.\n\nThe paper is very theoretical and technical with abundance of math. The paper does little to make the material accessible to wider audience. The presentation is technically excellent, although some parts show some looseness. I was able to follow maybe third of the material, and rest was too deep or cursory to be understood (which is ok). \n\nHowever, one should be able to follow the higher-level story regardless: what problem are we solving, why math is introduced, and what does the math achieve. The paper is not particularly clear on these aspects, and there is insufficient motivation and introduction to the theoretical presentation to be even following the higher-level story. There is a lot of math’itis where math happens but it’s not clear why. The material is also not clearly contextualised into wider literature. The paper is presenting probably too much material (experiments start at last page), and the manuscript would have improved by prioritising some of the material (and some to appendix) to allow for room to breath.\n\nI'm rating the paper weak accept due to its strong theoretical contributions, although the presentation has lots of room for improvement. I'm looking forward to hearing author comments and revise my review accordingly.\n\nQuestions:\n- Is TO a standard wording for this concept in the literature? How does this relate to pushforwards, or to FNO’s or to operator learning?\n- What’s the connection of IG to Fokker-Planck or to Ito’s lemma? Eq 7 looks like an expected Ito. Can you contextualise this a bit more?\n- The scope of the work is abstract: what task are you solving when you \"learn markovian semigroups\"? I think this paper is about learning operators, especially in SDE context, but what problems can we then solve by these operators or by learning semigroups? Are we learning the solution operators of PDEs, or something more general? \n- What kind of observations do we need?\n- The open problem is quite vague. The introduction states that unboundedness of IGs makes designing “estimators” “challenging”. What specific problems does the unboundedness pose? Why does it make things challenging? It seems that earlier methods can learn PDE operators just fine: how are they deficient?\n- In stochastic processes we usually operate within some boundary constraint, or within some bounded domain. Is the unboundedness of IG then a problem in those contexts? Can you give a bit more explanation what the unboundedness of IG means?", "questions_for_authors": "- Why the $p_t$ has two arguments? What is $y$? Should this be perhaps interpreted as $p_{t,s}(X_t=x, X_s=y)$ for $t > s$?\n- Is A_t a pushforward? Can you elaborate a bit more on what this means conceptually, and how it connects with IG\n- What does A^* mean?\n- The paper claims to study the transfer operators on a space of invariant measures. But in that space the TO seemingly does nothing: what is here to study?\n- $f$ is first just a functional, and then later it’s from the space of invariant measures. Can you explain if these are two different $f$’s? \n- In sec 2 L was an operator that turns functional f(x0) into a derivative of expected functional E[f(xt) | x0]. In sec 3 L is a mapping between two invariant measures. Are these two descriptions of L complementary, or did we just define a new, different L in sec 3?\n- I don’t understand why the H and \\calL have different metric structures. Surely if f,g are in both H and L, their inner products should be the same. Can you elaborate?\n- Why can you assume that the RKHS is a subset of calL? The connection between the RKHS and stationary measures is not clear here. Does the kernel span only invariant measure functions? If so, how can you know such a kernel exists or how is it defined?\n- It would be useful to expand example 2.2. to also describe it wrt A/f/L/etc \n- I’m a bit confused of the setting of sec 3 first paragraph. We introduce a bunch fo stuff, but there is little motivation or introduction. For instance, it’s not clear what the \\calL is, why we want the L, and why H is a subset of calL, and why we introduced kernels. All of this happens, but it's purpose is not described: what is the task/problem/context?\n- Eq 9 has phi, but usually that is unknown. How do we handle this?\n- In eq 9 we draw X0 from pi. But pi is a invariant distribution, which is where the process converges at some (late) time. Surely it’s not the initial time: so why do we then write X0?\n- The eq 9 has a H-norm. But the square error seems to be over just scalar evaluations of measures: why do we need the Hilbert norm for scalars?\n- Eq 10 shows that some function of X0 is an integral of Xt’s. What does this mean conceptually? Is this perhaps the solution operator in the hilbert space?\n- I struggle to follow eqs 9 and 10. I don’t see much connection to earlier material since notation changes from A/f/L/S to phi/psi/pi. Can you elaborate a bit more here?\n- The eq 10 has the exp form for the solution of the process. Is this a limitation? Does this apply to all kinds of SDEs or semigroups?\n- Is eq 11 a quadrature?\n- What are A_t_j? Where do these come from, and are they known?\n- What is h or v?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: **Summary.**

The topic of this work is estimation of the Infinitesimal Generator (IG) of a continuous-time Markov processes in $\mathbb{R}^d$ from a trajectory of observations.

Specifically, the authors introduce a new learning algorithm for the problem (with theoretical guarantees), which estimates a spectral decomposition of the (typically unbounded) IG under additional assumptions (e.g. geometrically ergodicity, sectoriality).
The contribution contrasts with existing methods which focused on estimating transfer operators instead.
The authors demonstrate the applicability and robustness of their approach through experiments.

**Strengths.**

* The reviewers asked for numerous technical clarifications in their reviews, but the response was very effective and no soundness concerns were reported.
* Strong theoretical analysis (nvbC, iVn3, PJzG).
* The method overcomes a key drawback of previous methods: it remains stable even for short time-lags.
* The central idea of focusing on the resolvent of the IG was considered to be novel and insightful.
* The proposed approach is applicable to a broad class of Markov processes (including but not limited to self-adjoint IG) and to high-dimensional spaces (an experiment was conducted with d=45).

**Weaknesses.**

* There were mixed opinions about clarity: the paper is neatly written (CrN8) but is too dense (iVn3); its readability could be improved (PJzG, nvbC).
* Experiments fail to demonstrate superiority of the method (nvbC, CrN8, iVn3).
* Impact on the general audience of ICML is questionable (PJzG, iVn3).

**Discussion and reviewer consensus.**

The reviewers unanimously agree that this paper should be accepted. No major concern remains after the author response. The author response to the only remaining concern about the broader impact of the results (PJzG, iVn3) was sufficient to convince me.

**Overall evaluation.**

This is serious work; I recommend acceptance. In my opinion, the theoretical contribution warrants publication, and the authors have done an honest effort to back-up their claims experimentally; remaining experimental shortcomings should not disqualify the paper.

**Recommendation to the authors.**

If the paper is accepted, I recommend the authors to make good use of the additional space to increase the clarity of their writing; additional experiments may be added in the supplementary material (refer also to the suggestions of nvbC regarding adding a table of notations).
Decision: Accept (poster)