Learning the Infinitesimal Generator of Stochastic Diffusion Processes

Authors: vladimir r kostic
karim  lounici
hélène  halconruy
timothée  devergne
massimiliano  pontil
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 17:43:59

Abstract: We address data-driven learning of the infinitesimal generator of stochastic diffusion processes, essential for understanding numerical simulations of natural and physical systems. The unbounded nature of the generator poses significant challenges, rendering conventional analysis techniques for Hilbert-Schmidt operators ineffective. To overcome this, we introduce a novel framework based on the energy functional for these stochastic processes. Our approach integrates physical priors through an energy-based risk metric in both full and partial knowledge settings. We evaluate the statistical performance of a reduced-rank estimator in reproducing kernel Hilbert spaces (RKHS) in the partial knowledge setting. Notably, our approach provides learning bounds independent of the state space dimension and ensures non-spurious spectral estimation. Additionally, we elucidate how the distortion between the intrinsic energy-induced metric of the stochastic diffusion and the RKHS metric used for generator estimation impacts the spectral learning bounds.

Main paper: 
 introduction: 
 Continuous-time processes are often modeled using ordinary differential equations (ODEs), assuming deterministic dynamics. However, real systems in science and engineering that are modeled by ODEs are subject to unfeasible-to-model influences, necessitating the extension of deterministic models through stochastic differential equations (SDEs), see [40, 42] and references therein. SDEs are advantageous for modeling inherently random phenomena. For instance, in finance, they specify the stochastic process governing asset behavior, a crucial step in constructing pricing models for financial derivatives [41] . Another compelling application arises in atomistic simulations, where SDEs are used to model the evolution of atomic systems subjected to thermal fluctuations through the Boltzmann distribution [37] .A diverse range of SDEs can be represented as dX t = a(X t )dt + b(X t )dW t , where X 0 = x. Here, W denotes a (possibly multi-dimensional) Brownian motion, and the functions a and b are commonly γ 2 √ n ) O(n -1 4 ) O(n -α 2(α+β) ), α ≥ τ O(n -α 2(β+τ ) ), α < τ Spectral rates ✗ ✗ ✗ ✓ Time complexity O(n 2 + n 3/2 d) O(n 3 d 3 ) O(n 3 d 3 ) O(r n 2 d 2 )Table 1 : Comparison to previous kernel-based works on generator learning. Sample size is n, state-space dimension is d, γ is the regularization parameter of KKR and RRR and r is RRR rank parameter. Our learning bounds are derived in Theorem 2 where the parameters α, β, τ quantify the intrinsic difficulty of the problem and impact of kernel choice on learning IG.known as the drift and diffusion coefficients, respectively. Determining these coefficients from one or more trajectories, whether discretized or continuous, has been a key pursuit in "diffusion statistics" since the 1980s, as seen in works like [23, 30, 31] . However, uncovering the drift and diffusion coefficients alone does not reveal all the intrinsic properties of a complex system, such as the metastable states of Langevin dynamics in atomistic simulations [see e.g. 49, and references therein]. Consequently, there has been a shift and growing interest in the Infinitesimal Generator (IG) of an SDE, as its spectral decomposition offers a more in-depth understanding of system dynamics and behavior, thus providing a comprehensive picture beyond the mere identification of coefficients.Recovering the spectral decomposition of the IG can theoretically be achieved by exploiting the wellstudied Transfer Operators (TO) [see 27, and references therein]. TOs represent the average evolution of state functions (observables) over time and, being linear and amenable to spectral decomposition under certain conditions, they offer a valuable means to interpret and analyze nonlinear systems. However, they require evenly sampled data at a high rate, which may be impractical. Additionally, TO approaches are purely data-driven, complicating the incorporation of partial or full knowledge of an SDE into the learning process. Thus, there is growing interest in learning the IG directly from data, as it can handle uneven sampling and integrate SDE knowledge. The challenge lies in the IG being an unbounded operator, unlike TO which are often well-approximated by Hilbert-Schmidt operators with comprehensive statistical theory [28] . Unfortunately, the existing statistical theory collapses when applied to unbounded operators, prompting the need to completely rethink the problem Related work Extensive research has explored learning dynamical systems from data [see the monographs 7, 32, and references therein]. Analytical models for dynamics are often unavailable, motivating the need for data-driven methods. Two prominent approaches have emerged: deep learning-based methods [5, 12, 36] , effective for learning complex data representations but lacking statistical analysis, and kernel methods [2, 6, 10, 20, 22, 25, 27, 28, 50] , offering solid statistical guarantees for the estimation of the TO but requiring careful selection of the kernel function. A related question, tackling the challenging problem of learning invariant subspaces of the TO, has recently led to the development of several methodologies [22, 34, 38, 46, 52, 51] , some of which are based on deep canonical correlation analysis [3, 29] . In comparison, there have been significantly fewer works on learning the IG and only in very specific settings. In [55] a deep learning approach is developed for Langevin diffusion, while [24] presents an extended dynamical mode decomposition method for learning the generator and clarifies its connection to Galerkin's approximation. However, neither of these two works provides any learning guarantees. In this respect the only previous work we are aware of are [1] , revising the Galerkin method for Laplacian diffusion, [19] , presenting a kernel approach for general diffusion SDE with full knowledge of drift and diffusion coefficients, and [43] addressing Langevin diffusion. As highlighted in Table 1 , the bounds and analysis in these works are either restricted to specific SDEs or are incomplete. Notably, none of these works proposed an adequate framework to handle the unboundedness of the IG, resulting in an incomplete theoretical analysis and suboptimal rates, sometimes explicitly depending on state space dimension. Moreover, the estimators proposed in these works are prone to spurious eigenvalues and do not offer guarantees on the accurate estimation of eigenvalues and eigenfunctions.Contributions In summary, our main contributions are: 1) Proposing a fundamentally new idea to estimate the spectrum of self-adjoint generators of stable Itô SDE from a single trajectory. In contrast to all existing works, we exploit the geometry of the process via a novel energy (risk) functional; 2) In a certain sense we "fight fire (resolvent) with fire (generator)" to derive a new efficient learning method that is able to infer the best approximation of the resolvent of IG on the RKHS, independently of the time-sampling; 3) We prove the first IG spectral estimation finite sample bounds using (imperfect) partial knowledge, which notably, overcome the curse of dimensionality present in classical numerical methods; 4) Each important aspect of our learning method, especially in relation to the most relevant existing works, is empirically demonstrated to complement our theoretical analysis. 
background and the problem: 
 Our drive to estimate the eigenvalues of the infinitesimal generator for an SDE like (1) stems from its crucial role in characterizing dynamics in physical systems. This operator's closed form (3) , relies on the drift a and diffusion coefficient b, where we have partial knowledge: b is assumed to be known, but a is not. To compensate for the lack of prior knowledge about a, we introduce the system's energy as an additional known quantity. Below, we detail the mathematical concepts framing the problem (generator, spectrum, energy) exemplified through the Langevin and Cox-Ingersoll-Ross processes.For detailed list of notation used throughout the paper we refer Table 7 in the appendix.Transfer operator and infinitesimal generator A variety of physical, biological, and financial systems evolve through stochastic processes X = (X t ) t∈R + , where X t ∈ X ⊂ R d denotes the system's state at time t. A commonly employed model for such dynamics is captured by stochastic differential equations (SDEs) of the form dX t = a(X t )dt + b(X t )dW t andEQUATIONwhere x ∈ X , W = (W 1 t , . . . , W p t ) t∈R+ is a R p -dimensional (p ∈ N) standard Brownian motion, the drift a : X → R d , and the diffusion b : X → R d×p are assumed to be globally Lipschitz and sub-linear, so that the SDE (1) admits an unique solution X = (X t ) ⩾0 with values in (X , B(X )). Processes akin to equations like (1) are diverse, spanning models like Langevin and Cox-Ingersoll-Ross processes (see examples below), with broad applications in science and engineering. The process X is a continuous-time Markov process with almost surely continuous sample paths whose dynamics is described by a family of probability densities (p t ) t∈R+ and transfer operators (A t ) t∈R+ such that for all t ∈ R + , E ∈ B(X ), x ∈ X and measurable function f : X → R,P(X t ∈ E|X 0 = x) = E p t (x, y)dy and A t f = X f (y)p t (•, y)dy = E f (X t ) | X 0 = • . (2)Evaluating A t f at x yields the expectation of f starting from x and evolving until time t, making the transfer operator crucial for understanding X dynamics. The family (A t ) t∈R+ satisfies the fundamental semigroup equation A t+s = A t • A s , for s, t ∈ R + . Here, we focus on the transfer operator's effect on the set L 2 π (X ), a choice driven by the existence of an invariant measure π for A t on (X , B(X )) which satisfies A * t π = π for all t ∈ R + . Then, the process X is characterized by the infinitesimal generator L defined for every f ∈ L 2 π (X ) such that the limitLf = lim t→0 + (A t f -f )/t exists in L 2 π (X ).The operator L is closed on its domain dom(L) which is equal to the Sobolev spaceW 1,2 π (X ) = {f ∈ L 2 π (X ) | ∥f ∥ 2 W = ∥f ∥ L 2 π + ∥∇f ∥ L 2 π < ∞}.For SDE dynamics of the form (1), we can prove (see A.2 for details) that L is the second-order differential operator given, for any f ∈ L 2 π (X ), x ∈ X , byEQUATIONwhere∇ 2 f = (∂ 2 ij f ) i∈[d],j∈[p]denotes the Hessian matrix of f . Spectral decomposition Knowing only the drift a and diffusion b is not enough to compute (2) or to understand quantitative aspects of dynamical phase transitions, such as time scales and metastable states. The eigenvalues and eigenvectors of the generator are crucial for capturing these effects. To address the possible unbounded nature of L, one can turn to an auxiliary operator, the resolvent, which, under certain conditions, shares the same eigenfunctions as L and becomes compact. When it exists and is continuous for some µ ∈ C, the operator L µ = (µI -L) -1 is the resolvent of L and the corresponding resolvent set is defined byρ(L) = µ ∈ C | (µI -L) is bijective and L µ is continuousWe assume L has a compact resolvent, meaning ρ(L) ̸ = ∅ and there exists µ 0 ∈ ρ(L) such that (µ 0 I -L) -1 is compact. Under this assumption, and given that (L, dom(L)) is self-adjoint, we can prove the spectral decomposition of the generator (see A.1 for details) as follows:L = i∈N λ i f i ⊗ f i (4)where (λ i )i ∈ N are the eigenvalues of L, and the corresponding eigenfunctions f i ∈ L 2 π (X ), forming an orthonormal basis (f i ) i∈N , are also eigenfunctions of the transfer operator A t .Dirichlet forms and energy To handle the initial lack of knowledge about the drift a, we assume to have access to another quantity, called the energy, defined as E(f ) = lim t→0 X (f (f -A t f ))/tdπ for all functions f ∈ L 2 π (X ) for which this limit exists, defining in the way the domain dom(E). The associated Dirichlet form is the bilinear form defined by polarization for any f, g ∈ dom(E) byEQUATIONFor every f ∈ dom(E), we haveE(f ) = E(f, f ). As for every i ∈ N, 0 ≤ E(f i , f i ) = -X f i (λ i f i )dπ = -λ i ,we check that the eigenvalues of L are negative. To relate L to Dirichlet form, we assume there exists a Dirichlet operatorB = s ⊤ ∇ where s = [s 1 | . . . | s p ] : x ∈ X → s(x) = [s 1 (x) | . . . | s p (x)] ∈ R d×p is a smooth function s.t. Lf = s(s ⊤ ∇f ) = s(Bf ) and so that X (-Lf )gdπ = X (s(x)s(x) ⊤ ∇f (x)) ⊤ ∇g(x)π(dx) = X (Bf (x)) ⊤ (Bg(x))π(dx).We get that for any f ∈ dom(E)EQUATIONwhich is reminiscent of the expected value of the kinetic energy in quantum mechanics [17] .In the following, while we discuss in detail the examples of Overdamped Langevin and the Cox-Ingersoll-Ross processes, we briefly mention other ones that have a Dirichlet form: the Wright-Fisher diffusion (in dimension one), which can be defined in the context of population genetics and can be adapted to model interest rates, see [16] , the geometric Brownian motion which models the price process of a financial asset, the multi-dimensional Brownian motion (a=0) that corresponds to the heat equation, the transport processes associated with advection-diffusion equation, see [33] , and the process associated with Poisson's equation in electrostatics, see [26] .Example 1 (Langevin). Let k b , T ∈ R * + . The overdamped Langevin equation driven by a potential V : R d → R is given by dX t = -∇V (X t )dt + 2(k b T )dW t and X 0 = x, where k b and T respectively represent the coefficient of friction and the temperature of the system. Its infinitesimal generator L is defined by Lf= -∇V ⊤ ∇f + (k b T )∆f , for f ∈ W 1,2 π (X ). Since (-Lf )g dπ = -∇ (k b T )∇f (x) e -(k b T ) -1 V (x) Z g(x)dx = (k b T ) ∇f ⊤ ∇g dπ = f (-Lg) dπ, generator L is self-adjoint and associated to a gradient Dirichlet form with s(x) = (k b T ) 1/2 (δ ij ) i∈[d],j∈[p] . Example 2 (Cox-Ingersoll-Ross process). Let d = 1, a, b ∈ R, σ ∈ R * + . The Cox-Ingersoll-Ross process is solution of the SDE dX t = (a + bX t )dt + σ √ X t dW t and X 0 = x. Its infinitesimal generator L is defined for f ∈ L 2 π (X ) by Lf = (a + bx)∇f + σ 2 x 2 ∆f. By integration by parts, we can check that the generator L satisfies (-Lf )g dπ = σ 2 x 2 f ′ (x)g ′ (x) π(dx) = f (-Lg) dπ, and it is associated to a gradient Dirichlet form with s(x) = σ √ x/ √ 2.Learning in reproducing kernel Hilbert spaces (RKHSs) Throughout the paper we let H be an RKHS and let k : X × X → R be the associated kernel function. We let ϕ : X → H be a feature map [45] such that k(x, x ′ ) = ⟨ϕ(x), ϕ(x ′ )⟩ for all x, x ′ ∈ X . We consider RKHSs satisfying H ⊂ L 2 π (X ) [45, Chapter 4.3] , so that one can approximate L : L 2 π (X ) → L 2 π (X ) with an operator G : H → H. Notice that despite H ⊂ L 2 π (X ), the two spaces have different metric structures, that is for all f, g ∈ H, one in general has ⟨f, g⟩ H ̸ = ⟨f, g⟩ L 2 π . In order to handle this ambiguity, we introduce the injection operator S π : H → L 2 π (X ) such that for all f ∈ H, the object S π f is the element of L 2 π (X ) which is pointwise equal to f ∈ H, but endowed with the appropriate L 2 π norm.Then, the infinitesimal generator restricted to H is simply LS π which can be estimated by S π G for some G ∈ HS (H). This approach is based on the embedding ℓϕ : X → H of the generator in the RKHS that can be defined for kernels k ∈ C 2 (X × X ) whenever one knows drift and diffusion coefficients, see B, so that the reproducing property ⟨ℓϕ(x), h⟩ H = [LS π h](x) holds. Based on this observation, [19] developed empirical estimators of LS π that essentially minimize the risk∥LS π -S π G∥ 2 HS(H,L 2 π ) = E x∼π ∥ℓϕ(x) -G * ϕ(x)∥ 2 H .In scenarios where drift and diffusion coefficients are not known, then ℓϕ becomes non-computable. However if the process has the Dirichlet form (6) , one can still empirically estimate the Galerkin projection (S * π S π ) † S * π LS π onto H, as considered in [1] , which in fact minimizes the same risk. Yet this approach is problematic due to the unbounded nature of the generator and the associated estimators typically suffer from a large number of spurious eigenvalues around zero, making the estimation of physically most relevant eigenfunctions unreliable even for self-adjoint generators. Conversely, classical numerical methods can compute the leading part of a spectrum without spuriousness issues, suggesting that data-driven approaches should achieve similar reliability. In this paper, we address this problem by designing a novel notion of risk, leading to principled estimators designed to surmount these challenges.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes] ,Justification: see the conclusion section.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
novel statistical learning framework: 
 In this section, we tackle the challenges in developing suitable generator estimators highlighted earlier. To this end, we introduce a risk metric for resolvent estimation that can be efficiently minimized empirically, leading to good spectral estimation. Since L and (µI -L) -1 share the same eigenfunctions, the main idea is to learn the (compact) resolvent, which can be effectively approximated by finite-rank operators, instead of learning the generator directly. However, this approach is challenging due to the lack of closed analytical forms for the action of the resolvent.First, given µ > 0, in order to approximate the action of the resolvent on the RKHS by some operator G : H → H, we introduce its embedding χ µ : X → H via the reproducing property⟨χ µ (x), h⟩ H = [(µI -L) -1 S π h](x), formally given by χ µ (x) = ∞ 0 E[ϕ(X t )e -µt | X 0 = x]dt, see B for details. Using this notation, we aim to estimate [(µI -L) -1 S π h](x) ≈ [Gh](x), h ∈ H, i.e. the objective is to estimate χ µ (x) ≈ G * ϕ(x) π-a.e.An obvious metric for the risk would be the mean square error (MSE) w.r.t. distribution π of the data. However, this becomes intractable since in general χ µ is not computable in closed form when either full or partial knowledge of the process is at hand.To mitigate this issue, we introduce a different ambient space in which we study the resolventW µ π (X ) = {f ∈ dom(L) | ∥f ∥ 2 W µ π (X ) ≡ E µ [f ] = ⟨f, (µI-L)f ⟩ L 2 π < ∞},where the norm now balances the energy of an observable f : X → R w.r. . IndeedEQUATIONwhere the last equality holds for Dirichlet gradient form (6) , in which case W µ π (X ) is simply a weighted Sobolev space. Importantly, this energy functional can be empirically estimated from data sampled from π, whenever full knowledge, that is drift and diffusion coefficients of the SDE (1), or partial knowledge, i.e. the diffusion coefficient and Dirichlet operator B in (6) , is at hand. With that in mind, instead of the standard MSE, we introduce the energy-based risk functional as minEQUATIONDenoting by Z µ : H → W µ π (X ) the canonical injection, (8) can be equivalently written asEQUATIONπ ) , (9) where we abbreviated W = W µ π (X ) and used Z * µ = S * π (µI-L), recalling that Hilbert-Schmidt and spectral norms for operators A : H → W, are ∥A∥ HS(H,W) = tr(A * (µI-L)A) and ∥A∥ H→W = λ 1 (A * (µI-L)A) ≥ µ -1/2 ∥A∥ H→L 2 π . Therefore, (9) implies that the regularized energy norm, while dominating the classical L 2 π one, exerts a balancing effect on the estimation of the resolvent. This leads us to the first general result regarding the well-posedness of this framework. Proposition 1. Given µ > 0, let H ⊆ W µ π (X ) be the RKHS associated to kernel k ∈ C 2 (X ×X ) such that Z µ ∈ HS (H, W µ π (X )), and let P H be the orthogonal projector onto the closure of Im(Z µ ) ⊆ W µ π (X ). Then for every ε > 0 there exists a finite rank operator G :H → H such that R(G) ≤ ∥(I -P H )(µI -L) -1 Z µ ∥ 2 HS(H,W) + ε. Consequently, when k is universal, R(G) ≤ ε.The previous proposition reveals that whenever the hypothetical domain is dense in the true domain and the injection operator is Hilbert-Schmidt, there is no irreducible risk and one can find arbitrarily good finite rank approximations of the generator's resolvent. Note that Z µ ∈ HS (H, W µ π (X )) is equivalent to Z * µ Z µ = S * π (µI -L)S π being a trace class operator, which is assured for our Examples 1 and 2, see the discussion in E. Now, to address how minimization of the risk impacts the estimation of the spectral decomposition, let us define the operator norm error and the metric distortion functional, respectively, asEQUATION) Proposition 2. Let G = i∈[r] (µ-λ i ) -1 h i ⊗ g i be the spectral decomposition of G : H → H,where λ i ≥ λ i+1 and letEQUATIONwhere gap i is the difference between i-th and (i + 1)-th eigenvalue of (µI -L) -1 .Note that the estimation of the eigenfunctions is first obtained in the norm with respect to the energy space, and then transformed to theL 2 π -norm, as ∥f ∥ W ≥ √ µ∥f ∥ L 2 π , f ∈ W µ,γ. Therefore, by controlling the operator norm error and metric distortion (refer to C), we can guarantee accurate spectral estimation. Consequently, this allows us to approximately solve the SDE (1) starting from an initial conditionEQUATIONQuestion: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes] Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: as for other theoretical papers no particular concerns are present. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to 
empirical risk minimization: 
 In this section we address empirical risk minimization (ERM), deriving two main estimators. The first one minimizes empirical risk with Tikhonov regularization, while the second introduces additional regularization in the form of rank constraints.To present the estimator, we denote the covariance operator w.r.t. L 2 π by C = S * π S π = E x∼π [ϕ(x) ⊗ ϕ(x)], the cross-covariance operator T = S * π LS π , capturing correlations between input and the outputs of the generator, and the covariance operator W µ = Z * µ Z µ = S * π (µI-L)S π w.r.t. energy space W. All operators can be estimated from data, depending on the available prior knowledge. In this work we focus on the case when Dirichlet gradient form is known, i.e. we can define the embedding of the Dirichlet operatorEQUATIONMoreover, defining wϕ : X → H p+1 by wϕ(x) = [ √ µϕ(x), d 1 ϕ(x), . . . d p ϕ(x)] ⊤ ∈ R p+1 , we getEQUATIONRegularized risk Let us first introduce the regularized risk defined for some γ > 0 byEQUATIONwhich, after some algebra, can be written asR γ (G) = tr G * (µC-T +µγI)G-2CG+S * π (µI -L) -1 S π = ∥W -1/2 µ,γ C -W 1/2 µ,γ G∥ 2 HS(H) variance term + E µ [χ µ (•)] -∥W -1/2 µ,γ C∥ 2 HS(H) bias termwhere W µ,γ = W µ + µγI is the regularized covariance w.r.t. W.Hence, assuming the access to the dataset D n = (x i ) i∈[n] made of i.i.d. samples from the invariant distribution π, replacing the regularized energy E µ with its empirical estimate E µ leads to the regularized empirical risk functional expressed asEQUATIONwhere W µ,γ and C are estimated by their empirical counterparts W µ,γ and C, respectively, via (13) .Therefore, our regularized empirical risk minimization approach reduces toEQUATIONand we analyze two different estimators, the first one G µ,γ is obtained by minimizing regularized empirical risk (17) over all G ∈ HS (H), and, hence, the name Kernel Ridge Regression (KRR) of the generators resolvent. The second one G r µ,γ minimizes (17) subject to the (hard) constraint that G is at most of (a priori fixed) rank r ∈ N and, hence, is named Reduced Rank Regression (RRR) of the generator's resolvent. Notice that when r = n, the two estimators coincide. After some algebra, one sees that both minimization problems have closed form solutionsEQUATIONwhere [[•]] r denotes the r-truncated SVD of a compact operator.To conclude this section, we show how to compute the eigenvalue decomposition of (18) . To this end, we define the sampling operators S : H → R n and Z µ :H → R (1+p)n by ( Sh) i = 1 √ n h(x i ), i ∈ [n],and( Z µ h) kn+i = √ µ √ n h(x i ), k = 0, i ∈ [n], 1 √ n s k (x i ) ⊤ Dh(x i ), k ∈ [p], i ∈ [n]. Further, let K = n -1 [k(x i , x j )] i,j∈[n]∈ R n×n be kernel Gram matrix, and introduce the Gram matrices N ∈ R n×pn and M ∈ R pn×pn whose elements, for k ∈EQUATION)We note that although we have introduced the above matrices via inner products in H, they can be readily computed via the kernel and its gradients knowing the Dirichlet form, see D. Theorem 1. Given µ > 0 and γ > 0, let J µ,γ = K-N(M+γµI) -1 N ⊤ +γI. Let ( σ 2 i , v i ) i∈[r] be the leading eigenpairs of the following generalized eigenvalue problemEQUATIONDenoting V r = [v 1 |. . .| v r ] ∈ R n×r and Σ r = diag( σ 1 , . . . , σ r ), if (ν i , w ℓ i , w r i ) i∈[r] are eigentriplets of matrix V ⊤ r V r Σ 2 r ∈ R r×r , then the eigenvalue decomposition the RRR estimator G r µ,γ = Z µ U r V ⊤ t S is given by G r µ,γ = i∈[r] (µ -λ i ) -1 h i ⊗ g i , where λ i = µ -1/ν i , g i = ν -1/2 i S * V r w ℓ i and h i = Z * µ U r w r i for U r = (µγ) -1 [µ -1/2 I | -N(M+γµI) -1 ] ⊤ (KV r -µV r Σ 2 r ) ∈ R (1+p)n×r .The main computational cost of our method, in view of ( 12), to solve SDE (1) lies in the implicit inversion of J µ,γ when solving (20) . When computed with direct solvers this inversion is of the order O(n 3 p 3 ), however leveraging on the fact that µJ µ,γ is Schur's complement of the (1+p)n×(1+p)n symmetric positive definite matrix and using classical iterative solvers, like Lanczos or the generalized Davidson method, when r ≪ n this cost can significantly be reduced to O(r n 2 p 2 ), c.f. [18] . 
spectral learning bounds: 
 Recalling Prop. 2, in order to obtain the bounds on eigenvalues and eigenfunctions of the generator, it suffices to analyze the learning rates for the operator norm error E and metric distortion η. For this purpose, we analyze the operator norm error of empirical estimator G r µ,γ using the decompositionE( G) ≤ ∥(µI -L) -1 Z µ -Z µ G µ,γ ∥ H→W regularization bias +∥Z µ (G µ,γ -G r µ,γ )∥ H→W rank reduction bias +∥Z µ (G r µ,γ -G r µ,γ )∥ H→W estimator's variance, where G µ,γ = W -1 µ,γ C is the minimizer of the full (i.e. without rank constraint) Tikhonov regularized risk andG r µ,γ = W -1/2 µ,γ [[W -1/2 µ,γ C]] ris the population version of the empirical estimator G r µ,γ . Note that, while the last two terms in the above decomposition depend on the estimator, the first term depends only on the choice of H and the regularity of L w.r.t. H. In this work we focus on the classical kernel-based learning where one chooses a universal kernel [45, Chapter 4] and controls the regularization bias with a regularity condition. For details see Rem. 2 of App. E.2. Let µ > 0 be a prescribed parameter, we make following assumptions to quantify the difficulty of the learning problem:(BK) Boundedness. There exists c W >0 such that ess supx∼π ∥wϕ(x)∥ 2 ≤ c W , i.e. wϕ ∈ L ∞ π (X , H 1+p ); (RC) Regularity. For some α ∈ (0, 2] there exists c α > 0 such that C 2 ⪯ c 2 α W 1+α µ ; (SD) Spectral Decay. There exists β ∈ (0, 1] and c β > 0 s.t. λ j (W µ ) ≤ c β j -1/β , for all j ∈ J.The above assumptions, discussed in more details in App E.1, are in the spirit of state-of-the-art analysis of statistical learning theory of classical regression in RKHS spaces [14] , recently extended to regression of transfer operators [35, 28] . In our setting, instead of relying on the injection into L 2 π as in the case of transfer operators, the relevant object is the injection Z µ into the energy space W.The first assumption (BK) is the main limiting factor of our approach. Indeed, since ∥wϕ(x)∥ 2 = µ∥ϕ(x)∥ 2 + k∈[p] ∥d k ϕ(x)∥ 2 , apart from needing the kernel to be bounded, we also need the Dirichlet form embedding to be bounded. Essentially, this means that the Dirichlet coefficients are not growing too fast w.r.t. the kernel's gradients decay. Having this, we assure that Z µ ∈ HS (H, W µ π (X )) which implies the operator norm error can be controlled.Another key difference between generator and transfer operator regression is that the covariance w.r.t. the domain of the operator becomesW µ = Z * µ Z µ instead of C = S * π S π .On the other hand, the "cross-covariance" that captures now RKHS cross-correlations of the resolvent w.r.t domain is simplyZ * µ (µI -L) -1 Z µ = S * π (µI-L)(µI -L) -1 S π = C.With this in mind, (RC) corresponds to the regularity condition in [28] and it quantifies the relationship between the hypothesis class (bounded operators in H) and the object of interest (µI -L) -1 . Indeed, if L has eigenfunctions that belong to α-interpolation space between H and W µ π (X ), (RC) holds true. In particular, if f i ∈ H for all i ≥ 2 (constant eigenfunction excluded), one has that α ≥ 1 (c.f. Proposition 7 in App. E.1.Finally, (SD) quantifies the "size" of the hypothetical domain H within the true domain W µ π (X ) via the effective dimension tr(W -1 µ,γ W µ ) ≤ c β (µγ) -β , which, due to (BK), leads to another notion, known as the embedding property, quantifying the relationship between H and essentially bounded functions in the domain of the operator (KE) Kernel Embedding. There exists τ ∈ [β, 1] and such that c τ = ess supEQUATIONwhere Z µ = j∈J σ j z j ⊗ h j is the SVD of the injection operator Z µ : H → W µ π (X ). Using the above assumptions, we prove generalization bound for the RRR estimator, notably addressing the general case when the prior knowledge of the Dirichlet coefficient is inexact, i.e.(DF) Dirichlet form. For some ϵ∈[0, 1) there existss ϵ : R d →R p so that (1-ϵ)s ⊤ s⪯s ⊤ ϵ s ϵ ⪯(1+ϵ)s ⊤ s holds π-a.e., where s : R d →R p is such that L=B * B for B=s ⊤ ∇ : L 2 π (X ) → [L 2 π (X )] p .Recalling the form of IG in Equation 3, the Dirichlet form in (DF) of a self-adjoint generator exists whenever the positive definite diffusion part satisfies uniform ellipticity conditions and the drift term allows integration by parts, leading to s(x)=b(x)/ √ 2. Thus, to obtain the partial knowledge we need, it is enough to estimate the diffusion function s ϵ with the relative error bound in (DF), see App. E.1. Theorem 2. Let (DF), (RC), (SD), and (KE) hold for some ϵ ∈ [0, 1), α ∈ (0, 2], β ∈ (0, 1] and τ ∈ [β, 1], respectively, and let cl(Im(S π )) = L 2 π (X ). Denoting λ ⋆ k = λ k (S * π (µI -L) -1 S π ), k ∈ N,and given δ ∈ (0, 1) and r ∈ [n], if RRR estimator is built from the Dirichlet coefficients s ϵ : R d → R p and γ≍n -1 α+β and ε ⋆ n =n -α 2(α+β) when α ≥τ, or γ≍n -1 β+τ and ε ⋆ n =n -α 2(β+τ ) when α < τ, (22) then there exists a constant c > 0, depending only on H and gap λ ⋆ r -λ ⋆ r+1 > 0, such that for large enough n ≥ r with probability at least 1 -δ in the i.i.d. draw of D n from π it holds thatEQUATION)Proof sketch. The regularization bias is bounded by c α γ α/2 by Prop. 9 of App. E.2, the rank reduction bias is upper bounded by λ r+1 (S * π (µI -L) -1 S π ), while in the exact knowledge case (ϵ=0) the bounds on the variance terms critically rely on the well-known perturbation result for spectral projectors reported in Prop. 4, App. A. The latter is then chained to Pinelis-Sakhanenko's inequality and Minsker's inequality for self-adjoint HS-operators, Props. 12 and 13 in App. E.3.1, respectively. When the knowledge is not exact, that is ϵ > 0 in (DF), this relative bound implies that ϵ W µ ⪯ W µ -W ϵ µ ⪯ ϵ W µ , where the empirical covariance with the inexact Dirichlet coefficient s ϵ is denoted by W ϵ µ . This allows one to control the additional approximation error in the analysis of variance, paying the price of additive term ϵ. Combining the bias and variance terms, we obtain the balancing equations for the regularization parameter and then the next result follows. □ First, note that the learning rate (23) implies the L 2 π -norm learning rate. Moreover, for α ≥ τ , it matches information theoretic lower bounds for transfer operator learning upon replacing parameters α, β and τ related to the space W with their L 2 π analogues [28] , see App. E.6. This motivates the development of the first mini-max optimality for the IG learning, for which our results are an important first step. Next, remark that Theorem 2 guarantees the reliability of fully data-driven methods when the diffusion coefficients are not known but estimated. Furthermore, when b is constant or linear (e.g. Overdamped Langevin and CIR), the classical estimation bounds coincide with the relative error bound of assumption (DF).To conclude this section, we address the spectral learning bounds steaming from the Prop. 2. The main task to do so is to control the metric distortions, which we show how in App. E.5. In this context, an additional assumption α ≥ 1 is needed, since otherwise the metric distortions can blow-up due to eigenfunctions being out of the RKHS space. Importantly, our analysis reveals thatEQUATIONwhere s i = σ i η i is the empirical spectral bias that informs how good is the estimation of the particular eigenpair is (see Fig. 1 a )), σ i being given in (20) andEQUATIONreveals that our data-driven method for spectral decomposition of differential operator L does not suffer from the curse of dimensionality as present in the classical numerical methods, see App. E.6. 
experiments: 
 In this section, we showcase the key features of our method outlined in Table 1 . We demonstrate that our approach: (1) avoids the spurious effects noted in other IG methods [19, 1] , (2) is more effective than transfer operator methods [27] , and (3) validates our bounds in a prediction task for a model with a non-constant diffusion term. Further details are available in Appendix F.We first investigate the overdamped Langevin dynamics in a potential that presents four different wells, two principal wells and then in each of them two smaller wells, given by V (x) = 4(x 8 +0.8 exp(-80(x 2 ))+0.2 exp(-80(x-0.5) 2 )+0.5 exp(-40(x+0.5) 2 )). This leads to three relevant eigenpairs: the slowest mode corresponds to the transition between the two principal wells, while the others two capture transitions between the smaller wells. In Fig. 1 panel a), we show that the empirical bias s 1 = σ 1 η 1 allows us to choose the hyperparameters of the model, that is higher empirical bias coincides with unreliable estimation of the operator's eigenfunction. In panel b) we observe how it varies w.r.t. land-scale (y-axis) and regularization γ (x-axis) hyperparameters showcasing the robustness of the model, see also Fig. 3 of App. F for hyperparameter µ. Further, in panel c) we show the consistency of our model with the true Boltzmann distribution. Namely, we use our model to forecast the conditional probability density function (pdf) of the system. We perform the same procedure with prefect knowledge and imperfect diffusion coefficient estimated from data. We also report that if the same approach is used with the method described in [19, 43, 1] , no dynamical quantity can be forecast due to the presence of numerous spurious eigenpairs, which prevent the system from relaxing towards the Boltzmann distribution. This issue is further illustrated in panel d) where we show how, contrary to KRR method of [19, 1] , we avoid spuriousness in the estimation of eigenvalues.CIR model Next, with the CIR model we show that our method is not limited to Langevin process with constant diffusion. For this process, the conditional expectation of the state X t is analytically known. We can thus compare the prediction of our model with respect to this expectation using We have trained our method on a real 30-year US mortgage rates dataset and contrasted it with the fitted CIR model using continuous ranked probability scores that are estimated from the forecasts obtained by of each of them, see panel f) of Fig. 1 . Each model has been trained using data from January 2009 to December 2016. The initial condition was the last week of December 2016 and the predictions were made for the years 2017 and 2018. Since the dataset is real, we used the imperfect partial knowledge, that is, for our method, we estimated the diffusion coefficient only via a least squares calibration of a CIR model over the training set. This allows more flexibility on the drift term in our model.We next study Langevin dynamics under more challenging conditions: the Muller-Brown potential. Panels g)-h)-i) of Fig. 1 depict the second eigenfunction obtained by our method compared to the ground truth one, as well as the one found by the transfer operator approach, with the same number of samples. Notably, our physics informed approach outperforms transfer operator learning for this task. Note that with different lag times, we were able to recover this second eigenfunction. 
conclusion: 
 We developed a novel energy-based framework for learning the Infinitesimal Generator of stochastic diffusion SDEs using kernel methods. Our approach integrates physical priors, achieves fast error rates, and provides the first spectral learning guarantees for generator learning. A limitation is its computational complexity, scaling as n 2 d 2 . Future work will explore alternative methods to enhance computational efficiency and investigate a broader suite of SDEs beyond stochastic diffusion. 
 
 Appendix: This appendix includes additional background on SDE and RKHS, proofs of the results omitted in the main body and information about the numerical experiments. notation meaning notation meaning ∧ minimum ∨ maximum [[ • ]] r r-A.1 Basics on operator theory for Markov processes.We provide here some basics on operator theory for Markov processes. Let X ⊂ R d (d ∈ N) and (X t ) t∈R+ be a X -valued time-homogeneous Markov process defined on a filtered probability space (Ω, F, (F t ) t∈R+ , P) where F t = σ(X s , s ≤ t) is the natural filtration of (X t ) t∈R+ . The dynamics of X can be described through of a family of probability densities (p t ) t∈R+ such that for all t ∈ R + , E ∈ B(X ),P(X t ∈ E|X 0 = x) = E p t (x, y)dy.Let G be a set of real valued and measurable functions on X . For any t ∈ R + the transfer operator (TO) A t : G → G maps any measurable function f ∈ G toEQUATIONIn theory of Markov processes, the family (A t ) t∈R+ is referred to as the Markov semigroup associated to the process X.Remark 1. A possible choice is G = L ∞ (X ).Here, we are interested in another choice of G related to the existence of an invariant measure for A t i.e., a σ-finite measure π on (X , B(X )) such that P * t π = π for any t ∈ R + . In that case, we can choose G = L 2 π (X ) so that, for all f ∈ L 2 π (X ), P t f converges to f in L 2 π (X ) as t goes to 0. Note that P 0 f = f and P ∞ f = f dX for a suitable integrable function f : X → R.Within the existence of this invariant measure π, the process X is then characterized by the infinitesimal generator (IG) L : L 2 π (X ) → L 2 π (X ) of the family (A t ) t∈R+ defined byEQUATIONIn other words, L characterizes the linear differential equation ∂ t A t f = LA t f satisfied by the transfer operator. The domain of L denoted dom(L) coincides with the the Sobolev space W 1,2 π (X )W 1,2 π (X ) = {f ∈ L 2 π (X ) | ∥f ∥ 2 W = ∥f ∥ L 2 π + ∥∇f ∥ L 2 π < ∞}.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, under certain conditions, shares the same eigenfunctions as L and becomes compact. The following result can be found in Yosida's book ( [53] , Chap. IX) : For every µ > 0, the operator (µI -L) admits an inverse L µ = (µI -L) -1 that is a continuous operator on X and(µI -L) -1 = ∞ 0 e -µt A t dt.The operator L µ is the resolvent of L and the corresponding resolvent set of L is defined byρ(L) = µ ∈ C | (µI -L) is bijective and L µ is continuous .In fact, ρ(L) contains all real positive numbers and (µI -L) -1 is bounded. Here, we assume that L has compact resolvent, i.e. ρ(L) ̸ = ∅ and there exists µ 0 ∈ ρ(L) such that L µ0 is compact. Note that, through the resolvent identity, this implies the compactness of all resolvents L µ for µ ∈ ρ(L). Let then µ ∈ ρ(L). As (L, dom(L)) is assumed to be self-adjoint, so does L µ , so that L µ is both compact and self-adjoint. Then, its spectrum Sp(L µ ) = C \ ρ(L µ ) is purely discrete and consists of isolated eigenvalues (λ µ i ) i∈N such that |λ µ i | → 0 associated with an orthonormal basis (f i ) i∈N (see [53] , chapter XI). In other words, the spectral decomposition of the resolvent writesL µ = i∈N λ µ i f i ⊗ f iwhere the functions f i ∈ L 2 π (X ) are also eigenfunctions of the operator L we getL = i∈N λ i f i ⊗ f i . Example 3 (Langevin). Let (k b T ) ∈ R.The overdamped Langevin equation driven by a potential V : R d → R is given by dX t = -∇V (X t )dt+ 2(k b T )dW t and X 0 = x. Its invariant measure of the solution process X is the Boltzman distribution π(dx) = z -1 e -(k b T ) -1 V (x) dx where z denotes a normalizing constant. Its infinitesimal generator L is defined by Lf = -∇V ⊤ ∇f + (k b T )∆f , for f ∈ W 1,2 π (X ), and if ∇V ∈ L 2 π (X ) is positive and coercive, it has compact resolvent. Finally, since(-Lf )g dπ = -∇ (k b T )∇f (x) e -(k b T ) -1 V (x) Z g(x)dx = (k b T ) ∇f ⊤ ∇g dπ = f (-Lg) dπ, generator L is self-adjoint and associated to a gradient Dirichlet form with s(x) = (k b T ) 1/2 (δ ij ) i∈[d],j∈[p] . Example 4 (Cox-Ingersoll-Ross process). Let d = 1, a, b ∈ R, σ ∈ R * + .The Cox-Ingersoll-Ross process is solution of the SDE dX t = (a + bX t )dt + σ √ X t dW t and X 0 = x. Its invariant measure π is a Gamma distribution with shape parameter a/σ 2 and rate parameter b/σ 2 . Its infinitesimal generator L is defined for f ∈ L 2 π (X ) by Lf = (a + bx)∇f + σ 2 x 2 ∆f . Note that by integration by parts, we can check that the generator L satisfies (-Lf )g dπ = σ 2 x 2 f ′ (x)g ′ (x) π(dx) = f (-Lg) dπ, and it is associated to a gradient Dirichlet form with s(x) = σ √ x/ √ 2.After defining the infinitesimal generator for Markov processes (see A.1), we provide its explicit form for solution processes of equations like (1) . Given a smooth function f ∈ C 2 (X , R), Itô's formula (see for instance [4] , B, p.495) provides fort ∈ R + , f (X t ) -f (X 0 ) = t 0 d i=1 ∂ i f (X s )dX i s + 1 2 t 0 d i,j=1 ∂ 2 ij f (X s )d⟨X i , X j ⟩ s = t 0 ∇f (X s ) ⊤ dX t + 1 2 t 0 Tr X ⊤ s (∇ 2 f )(X s )X s ds.Recalling (1), we getEQUATIONProvided f and b are smooth enough, the expectation of the last stochastic integral vanishes so that we getE[f (X t )|X 0 = x] = f (x) + t 0 E ∇f (X s ) ⊤ a(X s ) + 1 2 Tr b(X s ) ⊤ (∇ 2 f (X s ))b(X s ) X 0 = x dsRecalling that L = lim t→0 + (A t f -f )/t, we get for every x ∈ X ,EQUATIONwhich provides the closed formula for the IG associated with the solution process of (1).Recalling that for a bounded linear operator A on some Hilbert space H the resolvent set of the operator A is defined as ρ(A) = λ ∈ C | A -λI is bijective and (A -λI) -1 is continuous , and its spectrum Sp(A) = C \ {ρ(A)}, let λ ⊆ Sp(A) be isolated part of spectra, i.e. both λ and µ = Sp(A) \ λ are closed in Sp(A). Than, the Riesz spectral projector P λ : H → H is defined byEQUATION)where Γ is any contour in the resolvent set Res(A) with λ in its interior and separating λ from µ. Indeed, we have that P 2 λ = P λ and H = Im(P λ ) ⊕ Ker(P λ ) where Im(P λ ) and Ker(P λ ) are both invariant under A and Sp(A | Im(P λ ) ) = λ, Sp(A | Ker(P λ ) ) = µ. Moreover, P λ + P µ = I and P λ P µ = P µ P λ = 0.Finally if A is compact operator, then the Riesz-Schauder theorem, see e.g. [44] , assures that Sp(T ) is a discrete set having no limit points except possibly λ = 0. Moreover, for any nonzero λ ∈ Sp(T ), then λ is an eigenvalue (i.e. it belongs to the point spectrum) of finite multiplicity, and, hence, we can deduce the spectral decomposition in the formEQUATIONwhere geometric multiplicity of λ, r λ = rank(P λ ), is bounded by the algebraic multiplicity of λ. If additionally A is normal operator, i.e. AA * = A * A, then P λ = P * λ is orthogonal projector for each λ ∈ Sp(A) and P λ = r λ i=1 ψ i ⊗ ψ i , where ψ i are normalized eigenfunctions of A corresponding to λ and r λ is both algebraic and geometric multiplicity of λ.We conclude this section with well-known perturbation bounds for eigenfunctions and spectral projectors of self-adjoint compact operators. Proposition 3 ([11] ). Let A be compact self-adjoint operator on a separable Hilbert space H.Given a pair ( λ, f ) ∈ C × H such that ∥ f ∥ = 1, let λ be the eigenvalue of A that is closest to λ and let f be its normalized eigenfunction.If g = min{| λ -λ| | λ ∈ Sp(A) \ {λ}} > 0, then sin(∢( f , f )) ≤ ∥A f -λ f ∥/ g.). Let A and A be two compact operators on a separable Hilbert space. For nonempty index set J ⊂ N letgap J (A) = min {|λ i (A) -λ j (A)| | i ∈ N \ J, j ∈ J}denote the spectral gap w.r.t J and let P J and P J be the corresponding spectral projectors of A and A, respectively. If A is self-adjoint and for some ∥A -A∥ < gap J (A), then ∥P J -P J ∥ ≤ ∥A -A∥ gap J (A) .In section 2, we have defined the infinitesimal generator of a diffusion process and specified its form when associated with Dirichlet forms. These operators act on L 2 π (X ) or specific subsets of it. To develop our learning procedure, we need to understand these operators' actions when embedding into the RKHS, and define their versions for feature maps. IG and Dirichlet operator in RKHS. As a reminder, we consider H be an RKHS and let k : X ×X → R be the associated kernel function. The canonical feature map is denoted by ϕ(x) = k(x, •) for x ∈ X k(x, x ′ ) = ⟨ϕ(x), ϕ(x ′ )⟩ for all x, x ′ ∈ X . Assuming that k is square-integrable with respect to the measure π, we define the injection operator S π : H → L 2 π (X ) and its adjointS * π : L 2 π (X ) → H by S * π f = X f (x)ϕ(x)π(dx).As a preliminary step, we need to define the reproducing partial derivatives in RKHS, which we introduce via Mercer kernels. Definition 1 (Mercer kernel). A kernel function k : X × X → R is called a Mercer kernel if it is a continuous and symmetric function such that for any finite set of points {x 1 , . . . , x n } ⊂ X , the matrix (k(x i , x j )) n i,j=1 is positive semi-definite.Several standard kernels satisfy the Mercer property with s ≥ 1, including the Gaussian kernel which we will consider subsequently.For s ∈ N and m ∈ N, we define the index set I m s = {α ∈ N m : |α| ≤ s} where |α| = s j=1 α j , for α = (α 1 , . . . , α m ) ∈ N m . For a function f : R m → R and x = (x 1 , . . . , x m ) ∈ R m , we denote its partial derivative d α f at point x (if it exists) asD α f (x) = m j=1 D αj j f (x) = ∂ |α| ∂ α1 x 1 • • • ∂ αm x m f (x).For a function k ∈ C 2s (X × X ) with X ⊂ R d and α ∈ I d s , we defineEQUATIONandD (0,α) k(x, y) = ∂ |α| ∂ α1 y 1 • • • ∂ αm y m k(x, y).Theorem 3 (Theorem 1 in [56] ). Let s ∈ N, X ⊆ R m and k be a Mercer kernel such that k ∈ C 2s (X × X ) with corresponding RKHS H. Then the following hold:EQUATION)ii. A partial derivative reproducing property holds for α ∈ I m sEQUATIONTheorem 3 allows us to introduce the first and second order operators D and D 2 that act on any feature map ϕ(x) asDϕ(x) = ((D ei k) x ) i∈[d] and D 2 ϕ(x) = ((D ei+ej k) x ) i,j∈[d]where the (D ei k) x and (D ei+ej k) x can be defined via (32) . Then, we define the operator d that maps any feature map ϕ(x)to dϕ(x) = s(x) ⊤ Dϕ(x). Denote s ⊤ = [s 1 | . . . |s d ] : x ∈ X → s ⊤ (x) = [s 1 (x)| . . . |s d (x)] ∈ R p×d . We have s(x) ⊤ Dh(x) = d i=1 si (x)∂ i h(x) = d i=1 si (x)⟨D ei ϕ(x), h⟩ H = ⟨dϕ(x), h⟩ H , so that we can define the embedding of the Dirichlet operator B = s ⊤ ∇ : L 2 π (X ) → [L 2 π (X )] p into RKHS dϕ : X → H p via the reproducing property as ⟨dϕ(x), h⟩ H = [BS π h](x) = s(x) ⊤ Dh(x) ∈ R p , h ∈ H. In fact, dϕ is a p-dimensional vector with components d k ϕ : X → H given via ⟨d k ϕ(x), h⟩ H = s k (x) ⊤ Dh(x), k ∈ [p]. Then, we can defineT = -E x∼π [dϕ(x) ⊗ dϕ(x)] = -k∈[p] E x∼π [d k ϕ(x) ⊗ d k ϕ(x)].which captures correlations between input and the outputs of the generator in the RKHS. Defining wϕ : X → H p+1 by wϕ(x) = [ √ µϕ(x), d 1 ϕ(x), . . . d p ϕ(x)] ⊤ ∈ R d+1 , we can considerW µ = Z * µ Z µ = S * π (µI-L)S π = E x∼π [µϕ(x) ⊗ ϕ(x) + dϕ(x) ⊗ dϕ(x)] = E x∼π [wϕ(x) ⊗ wϕ(x)]which corresponds to the RKHS covariance operator w.r.t. energy space W. Examples. One way to ensure that the essential assumption ((KE)) holds is to show that wϕ fulfils the boundedness condition ((BK)), i.e. wϕ ∈ L ∞ π (X , H 1+p ). Let's show that the property holds true for the Damped Langevin (see Example 1) and the CIR process (see Example 2) if we consider the Radial Basis Function (RBF) kernel k(x, y) = k x (y) = exp(-κ∥x -y∥ 2 ) where κ > 0 is a free parameter that sets the "spread" of the kernel. As a reminder, for every x ∈ X , we have∥wϕ(x)∥ 2 = µ∥ϕ(x)∥ 2 + p k=1 ∥d k ϕ(x)∥ 2 with d k = s k (x) ⊤ Dϕ(x).Recalling that for the overdamped Langevin process s(x) = (k b T ) 1/2 (δ ij ) i∈[d],j∈[p], ss ⊤ is diagonal so that ⟨s(x) ⊤ Dk x (•), s(x) ⊤ Dk x (•)⟩ = 0 for every x ∈ X . As ∥ϕ(x)∥ 2 = 1, we get that for every x ∈ X , ∥wϕ(x)∥ 2 ≤ µ =: c H . Consider now the CIR process. We have d = p = 1 and s(x) = σ √ x/ √ 2 for any x ∈ X . For the very same reasons, ⟨s(x) ⊤ Dk x (•), s(x) ⊤ Dk x (•)⟩ = 0 for every x ∈ X , so that ∥wϕ(x)∥ 2 ≤ µ =: c H . In both Langevin and CIR cases, we have wϕ ∈ L ∞ π (X , H 1+p ) when considering an RBF kernel.In this section, we prove key perturbation result and discuss the properties of the metric distortion.We conclude this section with the approximation bound for arbitrary estimator G ∈ B r (H) that is the basis of the statistical bounds that follow. This result is a direct consequence of [27] and Davis-Khan spectral perturbation result for compact self-adjoint operators, [11] .In the framework of Koopman operator learning [28] , spectral bounds are expressed in terms of a distortion metric between the RKHS H and L 2 π (X ), corresponding to the cost incurred from observing the operator's action on the H rather than on its domain L 2 π (X ). Aligned with the risk definition ( 9), here we measure in a certain way the distortion between the H and W µ π (X ) as given in definitions in (10) .EQUATIONwhere gap i is the difference between i-th and (i + 1)-th eigenvalue of (µI -L) -1 .Proof. We first remark that∥((µI -L) -1 -(µ-λ i ) -1 I L 2 π (X ) ) -1 ∥ -1 H→W ≤ ∥((µI -L) -1 Z µ -Z µ G) h i ∥ H→W /∥Z µ h i ∥ ≤ E( G)η( h i ).Then, from the first inequality, using that (µI -L) -1 is self-adjoint as operator W µ π (X ) → W µ π (X ), we obtain the first bound in (11) .(µ-λ) -1 ∈ Sp((µI -L) -1 )\{(µ-λ i ) -1 }, |(µ-λ i ) -1 -(µ-λ) -1 |≥|(µ-λ i ) -1 -(µ-λ) -1 |-|(µ-λ i ) -1 -(µ-λ i ) -1 | we conclude that |(µ-λ i ) -1 -(µ-λ) -1 |≥|(µ-λ i ) -1 -(µ-λ) -1 |-E( G) η( h i ), and min{|(µ-λ i ) -1 -(µ-λ) -1 | | (µ-λ) -1 ∈ Sp((µI -L) -1 ) \ {(µ-λ i ) -1 }} ≥ gap i -E( G) η( h i ).So, applying Proposition 3, we obtainsin(∢( f i , f i )) ≤ ∥(µI -L) -1 f i -(µ-λ i ) -1 f i ∥ [gap i -E( G) η( h i )] + ≤ ∥((µI -L) -1 Z µ -Z µ G) h i ∥/∥Z µ h i ∥ [gap i -E( G) η( h i )] + ≤ E( G) η( h i ) [gap i -E( G) η( h i )] + . Since, clearly ∥ f i -f i ∥ 2 ≤ 2(1 -cos(∢( f i , f i )) ≤ 2 sin(∢( f i , f i )), the proof of the second bound is completed.Next, we adapt the [28, Proposition 1] to our setting as follows.Proposition 5. Let G ∈ B r (H). For all i ∈ [r] the metric distortion of h i w.r.t. energy space W µ π (X ) can be tightly bounded asEQUATIONThe Reduced Rank Regression (RRR) estimator is the exact minimizer of (17) under fixed rank constraint. Specifically, RRR is the minimizer G r µ,γ of R γ (G) within the set of bounded operators HS r (H) on H that have rank at most r. The regularization term γ∥G∥ 2 HS is added to ensure stability. The closed form solution of the empirical RRR estimator isEQUATIONwhile its population counterpart is given by Gr µ,γ = W -1/2 µ,γ [[W 1/2 µ,γ C]]] r .In order to prove Theorem 1, recall kernel matrices in (19) and defineEQUATION)Now we provide the explicit forms of the matrices N and M in the case of Langevin (see Example 1) and CIR (see Example 2) processes, considering an RBF kernel k(x, y) = k x (y) = exp(-κ∥x-y∥ 2 ).As a reminder (see (36) ), N ∈ R n×pn and M ∈ R pn×pn are Gram matrices whose elements, fork ∈ [1 + p], i, j ∈ [n] are given by N i,(k-1)n+j = n -1 ⟨ϕ(x i ), d k ϕ(x j )⟩ H and M (k-1)n+i,(ℓ-1)n+j = n -1 ⟨d k ϕ(x i ), d ℓ ϕ(x j )⟩ H ,where d k ϕ(x j ) = s k (x j ) ⊤ Dϕ(x j ) and Dϕ(x j ) = Dk(x j , •) is defined by (31) .For k ∈ [p], i, j ∈ [n], we have ⟨ϕ(x i ), d k ϕ(x j )⟩ H = ⟨k xi , s k (x j ) ⊤ Dk xj ⟩ H = k xi , s k (x j ) ⊤ lim h→0 k(•, x j + he l ) -k(•, x j ) h l∈[d] = lim h→0 k xi , s k (x j ) ⊤ k(•, x j + he l ) -k(•, x j ) h l∈[d] = s k (x j ) ⊤ lim h→0 k(x i , x j + he l ) -k(x i , x j ) h l∈[d] = s k (x j ) ⊤ D (0,e l ) k(x i , x j ) l∈[d] = 2γs k (x j ) ⊤ (x (l) i -x (l) j )k(x i , x j ) l∈[d]where we have used the continuity of the inner product to get the third line and the reproducing property to obtain the following one. Similarly, for k, ℓ ∈ [p], j ∈ [n], we get⟨d k ϕ(x i ), d ℓ ϕ(x j )⟩ H = s k (x i ) ⊤ D (e l ′ ,0) k(x i , •) l ′ ∈[d] , s ℓ (x j ) ⊤ lim h→0 k(•, x j + he l ) -k(•, x j ) h l∈[d] = lim h→0 s k (x i ) ⊤ D (e l ′ ,0) k(x i , •) l ′ ∈[d] , s ℓ (x j ) ⊤ k(•, x j + he l ) -k(•, x j ) h l∈[d] = lim h→0 D (e l ′ ,0) k(x i , •) l ′ ∈[d] , s k (x j )s ℓ (x j ) ⊤ k(•, x j + he l ) -k(•, x j ) h l∈[d] =s k (x i )s ℓ (x j ) ⊤ lim h→0 D (e l ′ ,0) k(x i , x j + he l ) -D (e l ′ ,0) k(x i , x j ) h l ′ ,l∈[d] =s k (x i )s ℓ (x j ) ⊤ (D l ′ ,l k(x i , x j )) l ′ ,l∈[d], where we have used the partial derivative reproducing (32) property and where we define for l ′ ̸ = l,EQUATION) and for l ∈ [d], D l,l k(x i , x j ) = lim h→0 D (e l ,0) k(x i , x j + he l ) -D (e l ,0) k(x i , x j ) h = -2γ lim h→0 (x (l) i -x (l) j -h)k(x i , x j + he l ) -(x (l) i -x (l) j )k(x i , x j ) h = 2γ -4γ 2 (x (l) i -x (l) j ) 2 k(x i , x j ) = 2γ[1 -2γ(x (l) i -x (l) j ) 2 ]k(x i , x j ). (38) For the overdamped Langevin (see Example 1), for k ≤ d, s k (x i ) = (k b T ) 1/2 e k and s k (x i )s ℓ (x j ) ⊤ = (k b T )I so that N i,(k-1)n+j = n -1 ⟨ϕ(x i ), d k ϕ(x j )⟩ H = 2γ(k b T ) 1/2 n -1 (x (k) i -x (k) j )k(x i , x j ) and M (k-1)n+i,(ℓ-1)n+j = n -1 ⟨d k ϕ(x j ), d ℓ ϕ(x j )⟩ H = (k b T )n -1 (D l ′ ,l k(x j , x j )) l ′ ,l∈[d] ,where the elements of D are given by ( 37) and (38) . For the CIR process (see Example 2) in dimension d = 1, we have s(x) = σ √ x/ √ 2. Then, N i,(k-1)n+j = n -1 ⟨ϕ(x i ), d k ϕ(x j )⟩ H = √ 2σγn -1 √ x j (x i -x j ) exp(-γ|x i -x j | 2 )andM (k-1)n+i,(ℓ-1)n+j = n -1 ⟨d k ϕ(x i ), d ℓ ϕ(x j )⟩ H = σ 2 γn -1 √ x i √ x j exp(-γ|x i -x j | 2 ).Based on the previous formulas, using Theorem 1, which we prove next, one can estimate generator's eigenpairs in practice.Theorem 1. Given µ > 0 and γ > 0, let J µ,γ = K-N(M+γµI) -1 N ⊤ +γI. Let ( σ 2 i , v i ) i∈[r]be the leading eigenpairs of the following generalized eigenvalue problemEQUATION)Denoting V r = [v 1 |. . .| v r ] ∈ R n×r and Σ r = diag( σ 1 , . . . , σ r ), if (ν i , w ℓ i , w r i ) i∈[r] are eigentriplets of matrix V ⊤ r V r Σ 2 r ∈ R r×r , then the eigenvalue decomposition the RRR estimator G r µ,γ = Z µ U r V ⊤ t S is given by G r µ,γ = i∈[r] (µ -λ i ) -1 h i ⊗ g i , where λ i = µ -1/ν i , g i = ν -1/2 i S * V r w ℓ i and h i = Z * µ U r w r i for U r = (µγ) -1 [µ -1/2 I | -N(M+γµI) -1 ] ⊤ (KV r -µV r Σ 2 r ) ∈ R (1+p)n×r .Proof. First, note that µ J µ,γ ≻ 0 is exactly Schurs's complement w.r.t. second diagonal block of F µ,γ ≻ 0, and that, due to block inversion lemma [18] , we have thatEQUATION)where A is some np × np matrix. The first step in computing the RRR estimator lies in computing a truncating SVD ofW -1/2 µ,γ C, that is C W -1 µ,γ Cq i = σ 2 i q i , i ∈ [r]. Now, using the low-rank eigenvalue problem formulation [18] , we have thatq i = S * v i and S W -1 µ,γ C S * v i = σ 2 i v i . Now, recalling that S = [µ -1/2 | 0] Z µ we obtain and that Z µ ( Z * µ Z µ + µγI) -1 =( Z µ Z * µ + µγI) -1 Z µ , we obtain [µ -1/2 I | 0]F -1 µ,γ F µ [µ -1/2 I | 0] ⊤ Kv i = σ 2 i v i, which after some algebra, using (39)µ -1 (I -γJ -1 µ,γ )Kv i = σ 2 i v i , i.e. µ -1 (I -γJ -1 µ,γ )KV r = V r Σ 2 r .By normalizing right singular value functionsq i of W -1/2 µ,γ C, that is by asking that ⟨q i , q i ⟩ H = v ⊤ i Kv i = 1, we obtain that [[ W -1/2 µ,γ C]] r = W -1/2 µ,γ CQ r Q * r , for Q r = [q 1 | . . . | q r ].In other words, we haveG µ,γ = W -1 µ,γ Z * µ [µ -1/2 I | 0] ⊤ S S * V r V ⊤ r S = Z * µ F -1 µ,γ [µ -1/2 I | 0] ⊤ KV r V ⊤ r S = Z * µ U r V ⊤ r S,where we used that J-1 µ,γ KV r = γ -1 [KV r -µV r Σ 2 r ].Once with this form, we apply [27, Theorem 2] to obtain the result.Finally, next result provides the reasoning for using empirical metric distortion of Theorem 23. Proposition 6. Under the assumptions of Theorem 1, for every i ∈ [r]EQUATION)EQUATION)Proof. First, note that (40) follows directly from Theorem 1. Next, since for every i ∈ [r],( η i ) -2 -(η( h i )) -2 = ⟨ h i , ( W µ -W µ ) h i )⟩ ∥ h i ∥ 2 ≤ ∥ W µ -W µ ∥, we obtain η -1 i -(η( h i )) -1 ≤ η -2 i -(η( h i )) -2 (η( h i )) -1 ∨ η -1 i ≤ η( h i ) ∧ η i ∥ W µ -W µ ∥.Recalling the form of IG in Equation 3, the Dirichlet form in (DF) for a self-adjoint L exists whenever the positive definite diffusion part satisfies uniform ellipticity conditions and the drift term allows integration by parts, leading to s(x) = b(x)/ √ 2. Thus, to obtain the partial knowledge we need, it is enough to estimate the diffusion function s ϵ that satisfies the relative error bound in (DF). When the sample paths are observed continuously, such s ϵ can be directly identified from these observations, making it a non-statistical problem. For discrete realizations of the process, the diffusion coefficient can be estimated non-parametrically using various methods, such as pathwise estimation by computing the variance of the increments over small intervals [21] , kernel-based methods [15] , local polynomial regression [13] . Remark, however, that the estimation of the drift (needed for the full knowledge) is a much more demanding task. Different methods are reviewed in [31] and references therein. A more recent approach [9] drawing inspiration from particle systems, consists in constructing estimates from several i.i.d. paths of the solution process.Next, observe that (BK) implies Z µ ∈ HS (H, W µ π (X )), which according to the spectral theorem for positive self-adjoint operators, has an SVD, i.e. there exists at most countable positive sequence (σ j ) j∈J , where J = {1, 2, . . . , } ⊆ N, and ortho-normal systems (z j ) j∈J and (h j ) j∈J of cl(Im(Z µ )) and Ker(Z µ ) ⊥ , respectively, such that Z µ h j = σ j z j and Z * µ z j = σ j h j , j ∈ J. Now, given α ≥ 0, let us define scaled injection operator Z µ,α : H → W µ π (X ) asEQUATIONClearly, we have that Z µ = Z µ,1 , while Im Z µ,0 = cl(Im(Z µ )). Next, we equip Im(Z µ,α ) with a norm ∥•∥ H,α to build an interpolation space[H] α =    f ∈ Im(Z µ,α ) | ∥f ∥ 2 H,α = j∈J σ -2α j ⟨f, z j ⟩ 2 W < ∞    ,noting that the inner product in W µ π (X ) is given by bilinear energy functional ⟨f,g⟩ W = µ⟨f, g⟩ L 2 π -⟨f, Lg⟩ L 2 π .We remark that for α = 1 the space [H] α is just an RKHS H seen as a subspace of W µ π (X ). Moreover, we have the following injections[H] α1 → [H] 1 → [H] α2 → [H] 0 = W µ π (X ),where α 1 ≥ 1 ≥ α 2 ≥ 0.In addition, from (BK) we also have that RKHS H can be embedded intoW µ,∞ π (X ) = {f ∈ W µ π (X ) | ∥f ∥ W µ,∞ π = ess sup x∼π [|f (x)| 2 -f (x)[Lf ](x)] < ∞} that is, for some τ ∈ (0, 1] [H] 1 → [H] τ → W µ,∞ π (X ) → W µ π (X ). Now, according to [14], if Z µ,τ,∞ : [H] τ → L ∞π (X ) denotes the injection operator, its boundedness implies the polynomial decay of the singular values of Z µ , i.e. σ 2 j (Z µ ) ≲ j -1/τ , j ∈ J, and the condition (KE) is assured. Assumption (SD) allows one to quantify the effective dimension of H in ambient space W µ π (X ), while the kernel embedding property (KE) allows one to estimate the norms of whitened feature maps, in our generator setting vector-valued since they define rank(1+p) operators on H,EQUATION)This object plays a key role in deriving the learning rates for regression problems, [29] and the following result is bounding it. Lemma 1. Let (KE) hold for some τ ∈ [β, 1] and c τ ∈ (0, ∞). Then,E x∼π ∥ξ(x)∥ 2 H 1+p ≤ c β β 1-β (µγ) -β , β < 1, c τ (µγ) -1 , β = 1.andEQUATIONProof. W.l.o.g. set µ = 1, observing that the only change in the proof is in scaling γ > 0. We first observe that for every j ∈ J from definition of wϕ and fact that h j (x) = [Z µ h j ](x) π-a.e., it holds thati∈[1+p] ⟨w i ϕ(x), h j ⟩ 2 =µ|h j (x)| 2 -h j (x)[LZ µ h j ](x)=µ|[Z µ h j ](x)| 2 -[Z µ h j ](x)[LZ µ h j ](x), π-a.e., implying that i∈[1+p] ⟨w i ϕ(x), h j ⟩ 2 ≤ σ j µ|z j (x)| 2 -z j (x)[Lz j ](x). So, for every τ > 0 ∥ξ(x)∥ 2 H 1+p = j∈J i∈[1+p] ⟨W -1/2 µ,γ w i ϕ(x), h j ⟩ 2 = j∈J i∈[1+p] 1 σ 2 j + γ ⟨w i ϕ(x), h j ⟩ 2 = j∈J i∈[1+p] σ 2(1-τ ) j σ 2 j + γ ⟨w i ϕ(x), h j ⟩ 2 σ 2 j σ 2τ j =γ -τ j∈J i∈[1+p] (σ 2 j γ -1 ) 1-τ σ 2 j γ -1 + 1 ⟨w i ϕ(x), h j ⟩ 2 σ 2 j σ 2τ j ≤γ -τ j∈J µ|h j (x)| 2 -h j (x)[LZ µ h j ](x) σ 2 j σ 2τ j =γ -τ j∈J (µ|z j (x)| 2 -z j (x)[Lz j ](x))σ 2τ j ,and, due to (21) , we obtain ∥ξ∥ 2 ∞ ≤ γ -τ c τ . On the other hand, we also have that tr(E x∼π [ξ(x) ⊗ ξ(x)]) = tr(W -1/2 µ,γ W µ W -1/2 µ,γ ) = tr(W -1 µ,γ W µ ), which is an effective dimension of the RKHS H in W µ π (X ). Therefore, following the proof of Fischer and Steinwart [14, Lemma 11] for classical covariances in L 2 π , we show that the bound on the effective dimension isEQUATIONFor the case β = 1, it suffices to see that∥z j ∥ W = E µ [z j ] ≤ ∥z j ∥ W µ,∞ π = ess sup x∼π [|z j (x)| 2 -z j (x)[LZ µ z j ](x)],and, hence tr(W -1 µ,γ W µ )≤γ -1 j∈N σ 2 j ∥z j ∥ 2 W =γ -1 j∈N σ 2 j E µ [z j ]≤γ -1 c τ .For β < 1 we can apply the same classical reasoning as in the proof of Proposition 3 of [8] . Proof. Note first that the resolvent (µI -L) -1 admits same eigenfunctions as the generator L, meaning that Im((µI -L) -1 Z µ ) ⊆ cl(Im(Z µ,α )). But according to [54, Theorem 2.2], this last condition is equivalent to Condition (RC).Finally we prove here Proposition 8Proposition 8. Given µ > 0, let H ⊆ W µ π (X ) be the RKHS associated to kernel k ∈ C 2 (X ×X ) such that Z µ ∈ HS (H, W µ π (X )), and let P H be the orthogonal projector onto the closure of Im(Z µ ) ⊆ W µ π (X ). Then for every ε > 0 there exists a finite rank operator G :H → H such that R(G) ≤ ∥(I -P H )(µI -L) -1 Z µ ∥ 2 HS(H,W) + ε. Consequently, when k is universal, R(G) ≤ ε.Proof. Recall first that since Z µ ∈ HS (H, W µ π (X )), according to the spectral theorem for positive self-adjoint operators, Z µ admits an SVD. Its form is provided in (42) taking α = 1.Now, recalling that [[•]] r denotes the r-truncated SVD, i.e. [[Z µ ]] r = j∈[r] σ j z j ⊗ h j , since ∥Z µ - [[Z µ ]] r ∥ 2 HS = j>r σ 2 j , for every δ > 0 there exists r ∈ N such that ∥Z µ -[[Z µ ]] r ∥ HS < µδ/3. Consequently since all the eigenvalues of L are non-positive, ∥(µI -L) -1 (Z µ -[[Z µ ]] r )∥ HS ≤ ∥Z µ -[[Z µ ]] r ∥ HS /µ ≤ δ/3. Next since Im(P H (µI -L) -1 Z µ ) ⊆ cl(Im(Z µ )), for any j ∈ [r], there exists g j ∈ H s.t. ∥P H (µI -L) -1 z j -Z µ g j ∥ ≤ δ 3r , and, denoting B r := j∈[r] σ j g j ⊗ h j we conclude ∥P H (µI -L) -1 [[Z µ ]] r -Z µ B r ∥ HS ≤ δ/3.Finally we recall that the set of non-defective matrices is dense in the space of matrices [47] , 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 ∈ B r (H) such that ∥G -B r ∥ HS < δ/(3σ 1 (Z µ )). So, we conclude∥(µI -L) -1 Z µ -Z µ G∥ HS ≤ ∥(I -P H )(µI -L) -1 Z µ ∥ HS + ∥(µI -L) -1 Z µ -[[(µI -L) -1 Z µ ]] r ∥ HS + ∥[[(µI -L) -1 Z µ ]] r -Z µ B r ∥ HS + ∥Z µ (G -B r )∥ HS ≤ ∥(I -P H )(µI -L) -1 Z µ ∥ HS + δ.Example 5. These three conditions depend on the process X (through its generator L) as well as the chosen RKHS. They are satisfied, for example, by choosing for k as a Gaussian kernel. Indeed, the sub-linearity conditions on A and B required to ensure the existence and uniqueness of the solution of the process of (1), also ensure that A and B are sufficiently 'nice' to fulfil, notably, condition ((BK)).Recalling the error decomposition and passing to the H and L 2 π -norms, we have thatEQUATIONand continue to prove the bound of the first term. Note that, while this proof technique is standard for operator learning [29, 35] , we present it here for the sake of completeness.Proposition 9. Let G µ,γ = W -1 µ,γ C for γ > 0, and P H : W µ π (X ) → W µ π (X ) be the orthogonal projector onto cl(Im(Z µ )). If the assumptions (BK), (SD) and (RC) hold, then ∥G µ,γ ∥ ≤ c α cEQUATIONEQUATIONProof. Recalling that P H := j∈J z j ⊗ z j , start by denoting the orthogonal projectors on the subspace of k leading left singular functions of Z µ asP k := j∈[k] z j ⊗ z j , respectively. Next, observe that C=Z * µ (µI-L) -1 Z µ and Z * µ W -1 µ,γ =Z * µ (Z * µ Z µ + µγI) -1 =(Z µ Z * µ + µγI) -1 Z µ . Therefore, (µI-L) -1 Z µ -Z µ G µ,γ =(I-Z µ W -1 µ,γ Z * µ )(µI-L) -1 Z µ =(I-(Z µ Z * µ +µγI) -1 Z µ Z * µ )(µI-L) -1 Z µ =µγ(Z µ Z * µ +µγI) -1 (µI-L) -1 Z µand, hence(µI-L) -1 Z µ -Z µ G µ,γ =   j∈J µγ σ 2 j + µγ z j ⊗ z j   (µI-L) -1 Z µ =   j∈J µγ (σ 2 j + µγ)σ j z j ⊗ (Z µ h j )   (µI-L) -1 Z µ =   j∈J µγ (σ 2 j + µγ)σ j z j ⊗ h j   C. Therefore, for every k ∈ J, ∥P k ((µI-L) -1 Z µ -Z µ G µ,γ )∥ 2 H→W becomes   j∈[k] µγ (σ 2 j +µγ)σ j z j ⊗ h j   C 2   j∈[k] γ (σ 2 j +µγ)σ j h j ⊗ z j   H→W ,which, due to (RC), implies that∥P k ((µI-L) -1 Z µ -Z µ G µ,γ )∥ H→W ≤c α j∈[k] µγ σ α j σ 2 j + µγ z j ⊗ z j W→W .On the other hand,j∈[k] µγ σ α j σ 2 j + µγ z j ⊗ z j =γ α 2 j∈[k] (σ 2 j (µγ) -1 ) α 2 σ 2 j (µγ) -1 + 1 z j ⊗ z j ⪯ (µγ) α 2 j∈[k] z j ⊗ z j ,where the inequality holds due to x s ≤ x + 1 for all x ≥ 0 and s ∈ [0, 1]. Since the norm of the projector equals one, we get∥P k ((µI-L) -1 Z µ -Z µ G µ,γ )∥ ≤ c α (µγ) α 2 . Next, observe that ∥(P H -P k )((µI-L) -1 Z µ -Z µ G µ,γ )∥ 2 H→W = j∈J\[k] µ 2 γ 2 (σ 2 j +µγ) 2 (Z * µ z j ) ⊗ (Z * µ z j ) H→H which is bounded by j∈J\[k] µ 2 γ 2 (σ 2 j +µγ) 2 ∥Z * µ z j ∥ 2 ≤ j∈J\[k] µ 2 γ 2 σ 2α j (σ 2 j +µγ) 2 ≤ j∈J\[k] σ 2α j .Using triangular inequality, for every k ∈ J we have that ∥P H ((µI-L)-1 Z µ -Z µ G µ,γ )∥ H→W is bounded by ∥P k ((µI-L) -1 Z µ -Z µ G µ,γ )∥ H→W +∥(P H -P k )((µI-L) -1 Z µ -Z µ G µ,γ )∥ H→Wand, therefore,∥P H ((µI-L) -1 Z µ -Z µ G µ,γ )∥ H→W ≤ c α (µγ) α 2 + j∈J\[k] (σ 2β j ) α β , and, hence, letting k → ∞ we obtain ∥P H (µI-L) -1 Z µ -Z µ G µ,γ ∥ ≤ c α (µγ) α 2 .Hence, (47) follows from triangular inequality.To estimate the ∥G µ,γ ∥, note that(RC) implies ∥G µ,γ ∥ ≤ c α ∥W -1 µ,γ W 1+α 2 µ ∥ and considering two cases. First, if (RC) holds for some α ∈ [1, 2], then, clearly ∥G µ,γ ∥ ≤ c α c (α-1)/2 W. On the other hand, if α ∈ (0, 1], thenσ 1+α j σ 2 j + µγ = (µγ) -1 σ 2 j (µγ) -1 1+α 2 σ 2 j (µγ) -1 + 1 ≤ (µγ) α-1 2 ,and, thus, ∥G µ,γ ∥ ≤ c α γ (α-1)/2 .Remark 2. Inequality (47) says that the regularization bias is comprised of a term depending on the choice of γ, and on a term depending on the "alignment" between H and Im((µI -L) -1 Z µ ).The term ∥(I -P H )(µI -L) -1 Z µ ∥ can be set to zero by two different approaches. One is choose a kernel which in some way minimizes ∥(I -P H )(µI -L) -1 Z µ ∥. Another is to choose a universal kernel [45, Chapter 4] , for which Im((µI -L) -1 S π ) ⊆ cl(Im(S π )). While we here develop theory for universal kernels, deep learning approaches that leverage on on our approach can be developed, which is the direction to pursue in future.In order to proceed with bounding the bias due to rank reduction for both considered estimators, we first provide auxiliary result.Proposition 10. Let B := W -1/2µ,γ C, let (RC) hold for some α ∈ (0, 2], and for j ∈ N denoteEQUATIONThen for everyEQUATIONProof. Start by observing thatB * B = [(µI-L) -1 Z µ ] * Z µ W -1 µ,γ Z * µ (µI-L) -1 Z µ = [(µI-L) -1 Z µ ] * (µI-L) -1 Z µ -µγ[(µI-L) -1 Z µ ] * (Z µ Z * µ + µγI) -1 (µI-L) -1 Z µ , implies that, using [(µI-L) -1 Z µ ] * =S π and [(µI-L) -1 Z µ ] * (µI-L) -1 Z µ =S * π (µI -L) -1 S π , S * π (µI -L) -1 S π - j∈J µγ σ 2 j + µγ (S * π z j ) ⊗ (S * π z j ) = B * B ⪯ Z * µ Z µ .Next, similarly to the above, for every k ∈ J, we havej∈[k] µγ σ 2 j + µγ (S * π z j )⊗(S * π z j ) ≤c 2 α j∈[k] σ 2α j σ 2 j (µγ) -1 + 1 z j ⊗z j =c 2 α j∈[k] (σ 2 j (µγ) -1 ) α/2 σ α j (µγ) α/2 σ 2 j (µγ) -1 + 1 z j ⊗z j ≤c 2 α γ α/2 ∥W µ ∥ α/2 ,andj∈J\[k] γ σ 2 j + µγ (S * π z j ) ⊗ (S * π z j ) ≤ c 2 α j∈J\[k] γ σ 2 j + µγ σ 2α j ≤ c 2 α j∈J\[k] (σ 2β j ) α/β .So, as before, letting k → ∞ we get the result.As a consequence, we obtain the bound for the rank reduction bias of the RRR method.Proposition 11 (RRR). Let (RC) hold for some α ∈ (0, 2]. Then the bias of G r µ,γ due to rank reduction, recalling (48) is bounded asEQUATIONProof. Observe that∥Z µ (G µ,γ -G r µ,γ )∥ ≤ ∥W 1/2 µ,γ (G µ,γ -G r µ,γ )∥ = ∥B -[[B]] r ∥ = σ r+1 (B) ≤ σ r+1 ((µI-L) -1 Z µ ) while ∥Z µ (G µ,γ -G r µ,γ )∥ ≥ ∥W 1/2 µ,γ (G µ,γ -G r µ,γ )∥ -(µγ) 1/2 ∥G µ,γ -G r µ,γ ∥ ≥ σ r+1 ((µI-L) -1 Z µ ) -c α ∥W µ ∥ α/4 (µγ) α/4 -2c α (µγ) (1∧α)/2 .All the statistical bounds we present will relay on two versions of Bernstein inequality. The first one is Pinelis and Sakhanenko inequality for random variables in a separable Hilbert space, see [8, Proposition 2] . Proposition 12. Let A i , i ∈ [n] be i.i.d copies of a random variable A in a separable Hilbert space with norm ∥•∥. If there exist constants Λ > 0 and σ > 0 such that for every m ≥ 2 E∥A∥ m ≤ 1 2 m!Λ m-2 σ 2 , then with probability at least 1 -δEQUATIONOn the other hand, we recall that in [39] , a dimension-free version of the non-commutative Bernstein inequality for finite-dimensional symmetric matrices is proposed (see also Theorem 7.3.1 in [48] for an easier to read and slightly improved version) as well as an extension to self-adjoint Hilbert-Schmidt operators on a separable Hilbert spaces. Proposition 13. Let A i , i ∈ [n] be i.i.d copies of a Hilbert-Schmidt operator A on the separable Hilbert space. Let ∥A∥ ≤ c almost surely, EA = 0 and let E[A 2 ] ⪯ V for some trace class operator V . Then with probability at least 1 -δEQUATIONwhereL A (δ) = log 4 δ + log tr(V ) ∥V ∥ .Proposition 14. Given δ > 0, with probability in the i.i.d. draw of (x i ) n i=1 from π, it holds thatP{∥ W µ -W µ ∥ ≤ ε n (δ)} ≥ 1 -δ,whereEQUATION)Proof. Proof follows directly from Proposition 13 applied to rank-(1+p) operators wϕ(x i ) ⊗ wϕ(x i ) using the fact that W µ = E[wϕ(x i ) ⊗ wϕ(x i )].Proposition 15. Let (KE) hold for τ ∈ [β, 1]. Given δ > 0, with probability in the i.i.d. draw of (x i ) n i=1 from π, it holds thatEQUATIONwhereEQUATIONandL 1 (γ, δ) = ln 4 δ + ln tr(W -1 µ,γ W µ ) ∥W -1 µ,γ W µ ∥ .Moreover,EQUATION)Proof. The idea is to apply Proposition 13 for operator ξ(x) ⊗ ξ(x), where ξ(x) is defined in (43) .Due to (1), we have that ∥A∥ ≤ ∥ξ∥ 2 ∞ ≤ (µγ) -τ c τ . On the other hand, we have thatE x∼π [ξ(x) ⊗ ξ(x)] 2 ⪯ ∥ξ∥ 2 ∞ E x∼π [ξ(x) ⊗ ξ(x)] = ∥ξ∥ 2 ∞ W -1/2 µ,γ W µ W -1/2µ,γ , and, hence (54) follows. To complete the proof, observe that∥I H -W -1/2 µ,γ W µ,γ W -1/2 µ,γ ∥ = ∥W -1/2 µ,γ (W µ -W µ )W -1/2 µ,γ ∥ ≤ ε 1 n (γ, δ), and, hence for ε 1 n (γ, δ) smaller than one we obtain∥W 1/2 µ,γ W -1 µ,γ W 1/2 µ,γ ∥ = ∥(W -1/2 µ,γ W µ,γ W -1/2 µ,γ ) -1 ∥ ≤ 1 1 -∥I H -W -1/2 µ,γ W µ,γ W -1/2 µ,γ ∥ .Proposition 16. Let (RC), (SD) and (KE) hold for some α ∈ (0, 2], β ∈ (0, 1] and τ ∈ [β, 1]. Given δ > 0, with probability in the i.i.d. draw of (x i ) n i=1 from π, it holdsP ∥W -1/2 µ,γ ( W µ -W µ )W -1 µ,γ C∥ HS ≤ ε 2 n (γ, δ) ≥ 1 -δ,whereEQUATIONProof. As before, w.l.o.g. set µ = 1. First, recall that HS (H) equipped with ∥•∥ HS is separable Hilbert space. Hence, we will apply Proposition 12 for A = ξ(x) ⊗ ψ(x), where ψ(x) = CW -1 µ,γ wϕ(x) ∈ H 1+p . To that end, observe thatE∥A∥ m HS =E [∥ξ(x)∥ m H 1+p ∥ψ(x)∥ m H 1+p ],and thatE[∥ξ(x)∥ m H 1+p ] ≤ 1 2 m! γ -τ /2 √ c τ m-2 tr(W -1 µ,γ W µ ) 2 .Recalling Lemma 1, the task is to bound ∥ψ(x)∥ 2H 1+p = i∈[1+p] ∥CW -1 µ,γ w i ϕ(x)∥ 2 H . Using (RC), we have that ∥ψ(x)∥ 2 H 1+p ≤c α i∈[1+p] ∥W (1+α)/2 µ W -1 µ,γ w i ϕ(x)∥ 2 H =c α j∈J i∈[1+p] ⟨W (1+α)/2 µ W -1 µ,γ w i ϕ(x), h j ⟩ 2 H .But, since⟨W (1+α)/2 µ W -1 µ,γ w i ϕ(x), h j ⟩ H = σ 1+α j σ 2 j + γ ⟨w i ϕ(x), h j ⟩ H ,expanding as in proof of Lemma 1 and using (KE), we have that∥ψ(x)∥ 2 H 1+p ≤c α j∈J i∈[1+p] σ (2+α-τ ) j σ 2 j + γ 2 ⟨w i ϕ(x), h j ⟩ 2 σ 2 j σ 2τ j ≤ γ -(τ -α) , α ≤ τ, c (α-τ ) W , α ≥ τ.But, since,B * B-B * B = C * W -1 µ,γ C-C * W -1 µ,γ C = ( C-C) * W -1 µ,γ C+C * W -1 µ,γ ( C-C)+C * ( W -1 µ,γ -W -1 µ,γ ) C denoting M = W -1/2 µ,γ ( C -C), N = W -1/2 µ,γ ( W µ -W µ ) and R = W 1/2 µ,γ ( G µ,γ -G µ,γ ), we have B * B -B * B = B * M + M * C 1/2 γ G µ,γ -B * N G µ,γ = B * M + (M * C 1/2 γ -B * N )( G µ,γ ± G µ,γ ) = B * M + M * B -B * N G µ,γ + (M * -B * N W -1/2 µ,γ )R = (G µ,γ ) * ( C -C)+( C-C)G µ,γ -(G µ,γ ) * ( W µ -W µ )G µ,γ +(M * + (G µ,γ ) * N * )R.Therefore, recalling that, due to (59EQUATIONNext, observe that• ∥( C-C)G µ,γ ∥≤∥( C-C)W 1/2 µ,γ ∥∥W 1/2 µ,γ C∥≤µ -1 ∥W 1/2 µ,γ ( C-C)∥ is bounded by Proposi- tion 17 • ∥M -N G µ,γ ∥ ≤ ∥W -1/2 µ,γ ( C-C)∥ + ∥W -1/2 µ,γ ( W µ -W µ )W -1µ,γ C∥ is bounded by Propositions 16 and 17,• ∥G * µ,γ ( W µ -W µ )G µ,γ ∥≤µ -1 ∥W -1/2 µ,γ ( W µ -W µ )W -1 µ,γ C∥ is bounded by Proposition 16.Therefore, using additionally Proposition 15 result follows.Remark that to bound singular values we can rely on the fact|σ i ( B) -σ i (B)| = |σ 2 i ( B) -σ 2 i (B)| σ i ( B) + σ i (B) ≤ |σ 2 i ( B) -σ 2 i (B)| σ i ( B) ∨ σ i (B).Recalling the notation B := W n (γ, δ) < 1, then with probability at least 1 -δ in the i.i.d. draw ofEQUATION) Proof. Start by observing that ∥Z µ (G r µ,γ -G r µ,γ )∥ ≤ ∥W 1/2 µ,γ (G r µ,γ -G r µ,γ )∥ and W 1/2 µ,γ (G r µ,γ -G r µ,γ ) =(W 1/2 µ,γ W -1 µ,γ W 1/2 µ,γ )• W -1/2 µ,γ ( W µ -W µ )G µ,γ P r + W -1/2 µ,γ ( C -C) P r + B( P r -P r ) .Using that the norm of orthogonal projector P is bounded by one and applying Propositions 15 and 16 together with Propositions 4 and 19 completes the proof.and we obtainEQUATIONTherefore, denotingγ n ≍ n -1 τ +β , α ≤ τ, n -1 α+β , α ≥ τ,andEQUATION)as a consequence the operator norm error bound in Theorem 2 holds, and we have the following result on the estimation of singular values of (µI-L) -1 Z µ , that is λ ⋆ i s by singular values of B, that is σ i s from Theorem 1. Proposition 21. Let (RC), (SD) and (KE) hold for some α ∈ (0, 2] and β ∈ (0, 1] and τ ∈ [β, 1], and define (68). Then, there exists a constant c > 0 such that for every given δ ∈ (0, 1), large enough n > r and with probability at least 1 -δ in the i.i.d. draw of (EQUATIONProof. The proof is direct consequence of Propositions 10 and 19 using (64)-(66).Finally, we conclude the proof of Theorem 2 showing the concentration of eigenpairs. Recalling Proposition 2, we need to combine the operator norm bound and metric distortion. Since, as indicated in Proposition 9, the population KRR estimator can grow in the operator norm whenever the regularity condition is violated α < 1, leading to possibly unbounded metric distortions w.r.t. increasing sample size, we restrict to the case α ≥ 1.First, combining (65)-( 67) and Proposition 21, we have thatE( G r µ,γ ) ≤ ( σ r+1 ∧ λ ⋆ r+1 ) + c ε ⋆ n ln δ -1 .On the other hand, from Propositions 40 and 14, we have that with failure probability δ( η i ) 2 -(η( h i )) -2 ≤ ε n (δ),which, if we can prove that η i and η( h i ) are bounded, concludes the proof for empirical spectral biases. To that end, recall that, c.f. Proposition 5, for RRR estimator we haveη( h i ) ≤ ∥ G r µ,γ ∥ σ r (Z µ G r µ,γ ) ≤ ∥ G µ,γ ∥ σ r ((µI -L) -1 Z µ ) -E( G r µ,γ ) ≤ ∥G µ,γ ∥ + ∥ G µ,γ -G µ,γ ∥ λ ⋆ r -E( G r µ,γ ) ≤ ∥G µ,γ ∥ + (µγ) -1/2 ∥W 1/2 µ,γ ( G µ,γ -G µ,γ )∥ λ ⋆ r -E( G r µ,γ ) ≲ 1 + n -1/2 γ -(β+1)/2 λ ⋆ r -n -1/2 γ -β/2 = 1 + n -α-1 α+β λ ⋆ r -n -α α+βwhere in the last inequality we have applied Propositions 9 and 18. Thus, using Proposition 40, the proof of Theorem 2 is concluded.We first discuss the learning rates reported in Table 1 . Notice that although the papers we compare to employ a different risk, our comparison remains meaningful because our energy-based risk measure provides an upper bound on their risk measure. This ensures that any upper bound derived with our risk also applies to theirs.• The learning bound for IG obtained is [1] covers only pure diffusion processes (Laplacian with constant weights). Their learning rate is non-parametric and depends on the state space dimension d in a counter-intuitive way O(n -d 2(d+1) ) in [1, Theorem 3], highlighting a potential limitation of their approach. In comparison, when we specify our RRR method with an RBF kernel (i.e. β = 0), we achieve a much faster parametric learning rate O(n -1/2 ).• The recent work of [43] covers Langevin processes via a kernel approach, but they derive a sub-optimal learning bound for IG of order O(n -1/4 ) in [43, Theorem 4.4] . For Langevin diffusions, our RRR method with an RBF kernel achieves a faster parametric rate O(n -1/2 ). Moreover, the computational complexity of their method is O(n 3 d 3 ), which limits its application in realistic molecular dynamics scenarios.• As for [19, Theorem 7] , although they considered general diffusions, they only derived a suboptimal bound for the variance component of their risk, with an explicit dependence on the dimension of the state space.Additionally, note that the mentioned works lack learning guarantees for eigenfunctions and eigenvalues. Notably, their methods are prone to the spurious eigenvalue phenomenon, requiring expert manual review of each eigenpair to select plausible ones.Next, we contrast our IG learning and well established TO learning. First, note that TO methods apply only to equally spaced data and the sampling frequency 1/∆t must be high enough to distinguish all relevant time-scales. Otherwise, since TO eigenvalues are e λi∆t , small spectral gaps complicate learning (see [28, Thm. 3] ). Conversely, our IG method, which uses gradient information, is timescale independent, handles irregularly spaced measurements, and does not rely on time discretizations. Indeed, recalling the risk functional in eq. ( 8), we see that the "label" of the model χµ (x) ≈ G * ϕ(x)is the action of the resolvent. Since this "label" is not computable, we "fight fire (resolvent) with fire (generator)", i.e. we use the energy norm of eq. ( 9) to rewrite regularized problem (8) . Crucially, this allows us to obtain estimators via energy covariance W µ in ( 14) that completely captures infinitesimal nature of the learning problem without needing time-lagged observations. This contrasts with TO methods, where choosing the time-lag ∆t is the major bottleneck in real applications. Indeed, to obtain data from an invariant distribution π one uses the trajectory data after some burn-in time needed to ensure that ergodic mean approximates well π. Then, the problem is reduced to studying only the dependence as is done e.g. in [27] using standard tools of β-mixing and the method of blocks. This allows one to obtain non-parametric learning bounds for TO methods where the effective sample size suffers from the multiplicative effect of the time-lag to achieve approximate independence. So, TO methods "waste" a lot of data negatively impacting statistical accuracy. Contrary to this, our method can be applied to data with larger time-lags (even irregularly spaced) so that effective sample size is close to the true one. All this results in better generalization, as shown in Fig. 1 g)-h)-i) where our IG estimator captures ground truth significantly better than TO for the same sample size, this generalization across all time scales incurs quadratic computational complexity w.r.t. state dimension and not statistical accuracy. Lastly, with our additional discussion on imperfect knowledge, our IG method can also be applied in a fully data-driven regime as TO methods.To conclude this section, we contrast our method with classical numerical methods for the spectral decomposition of differential operators, such as finite-elements (FEM). For FEM, the approximation error is |λ k -λ k | ≤ ch 2p λ k , where p is the polynomial degree used to construct finite elements, h is the mesh size, and c depends on the eigenfunctions' smoothness. As the number of mesh elements grows exponentially with d (i.e., ∼ h -d ), reducing h is the major bottleneck mitigated by computationally demanding adaptive higher order methods. On the other hand, our IG method that requires less or no knowledge has a quadratic impact of the d only on the computational complexity. Indeed, sample complexity depends on the effective dimension of the equilibrium distribution on the domain that can be much lower than d. Therefore, RRR for IG's resolvent doesn't suffer from the curse of dimensionality as the FEM does.All the experiments were performed on a workstation with 125.6 GiB of memory and AMD® Ryzen threadripper pro 3975wx 32-cores × 64 processor, no graphics card was used. The version of python used is Python 3.9.18. The choice of the kernel in all experiments was Gaussian RBF with specified length-scales, and the hyperparameters were chosen via cross-validation. The code to reproduce the results of the experiments can be found in the following repository: https://github.com/DevergneTimothee/GenLearn_kernel Four well potential For this experiment, we used an in-house code to simulate the system. The equations of motions were discretized using the Euler-Maruyama scheme with a timestep of 10 -4 . RRR was fitted using 1000 points, µ = 5 and γ = 10 -5 . The length scales used were 0.05 and 0.5. This experiment was reproduced 100 times leading to very small change in the estimation of the eigenfunctions. In Figure 2 we report the result of one of them. The reduced rank regression was performed with a rank of 5. Further, in Figure 3 we show the robustness of eigenfunctions w.r.t. choice of shift hyper parameter µ, as well as values of empirical bias for different values of length-scale and regularization hyperparameters. Concerning panel c) of Figure 1 where we show the consistency of our model with the true Boltzmann distribution. Namely, we use our model to forecast the conditional probability density function (pdf) of the system being in one bin, given it started at some point and after 100000 steps (so that it relaxes towards the equilibrium distribution). We bin the space into 50 bins and approximate the pdf via eq. ( 12), where h is characteristic function of a bin. In order to achieve an accurate estimation, we use the knowledge that the leading eigenvalue is zero. Additionally, we perform the same procedure with imperfect diffusion coefficient estimated from data (k b T = 0.45 ± 0.01 instead of true value 0.5a.u.) by looking at the variance of increments over 10 steps. Note that the result (green dotted lines) is unchanged and compares well with the analytical Boltzmann distribution (black lines). Finally, we report that if the same approach is used with the method described in [19, 43] , no dynamical quantity can be forecasted due to many spurious eigenpairs (see Figure 1a ) which prevents the system from relaxing towards the Boltzmann distribution.Muller Brown For this experiment, we used an in-house code to simulate the system. The equations of motions were discretized using the Euler-Maruyama scheme with a timestep of 10 -3 and a temperature of 2 (arbitrary units). RRR was fitted using 2000 points, µ = 1 and γ = 10 -5 . The length scale used was 0.6. The reduced rank regression was performed with a rank of 5CIR model For this experiment, we reproduced 100 times the simulations in order to obtain statistical uncertainties. We used a length scale of 0.5, µ = 1 and γ = 10 -6 . The reduced rank regression was performed with a rank of 2. Justification: we also provided a summary table comparing our results to key related works.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 


Decision: The reviewers agree that this is an interesting paper that addresses a different and important problem from a typical NeurIPS submissions. Counterbalancing its strengths as a theoretically grounded estimation framework, the reviewers list some concerns about motivation, hyperparameter tuning, experimental evaluation and presentation. I think the authors have addressed the first two well enough (including additional results). I think there is some truth to the latter two points, and while I personally though the  presentation in the main body was fine, the over submission crams a lot of theory material into the supplementary material, and I do not think the reviewers have or can be expected to go carefully over it. As such, I think this paper sits just at the borderline.

Reviews: review: paper_summary: paper_summary: This paper discusses learning the generator of stochastic diffusion processes in reproducing kernel Hilbert space. In particular,
section 2: background on the generator, Dirichlet form, energy, learning in RKHS, empirical risk in the Hilbert-Schmidt norm
section 3: introduce an energy-based risk functional for resolvent estimation which leads to spectral estimation
section 4: empirical risk minimization with Tikhonov regularization and rank constraints.
strength_weakness: strength_weakness: The weaknesses are mainly in the motivation (please see question 1 below) and the experiment (see questions 3 and 4 below). The authors did a great job in the theoretical comparison with previous work on TO and IG learning, however, the empirical comparison focuses on the spuriousness of eigenvalues and the recovery capability of the metastable states. To sufficiently demonstrate the practical advantage of the proposed method, other comparisons of the learned dynamics, phase transition, and invariant measure should also be included.
questions: questions: 1. In the paper, 'full knowledge' refers to given drift and diffusion coefficients, and 'partial knowledge' refers to the given diffusion coefficient and energy/Dirichlet operator. In practice, is it necessarily easier to obtain diffusion coefficient and energy/Dirichlet operator than to obtain drift and diffusion coefficients? 
In a fully data-driven scenario, one can first estimate drift and diffusion coefficients, and then use physics-informed method with 'full knowledge' obtained. Does this make the 'full knowledge' methods cited in the paper more flexible than the 'partial knowledge' method proposed?

2. As transfer operator A_t and the infinitesimal generator L satisfy the relation A_t = e^{Lt}, is there a consistency between the learned TO from cited TO methods with the learned IG from the proposed method?

3. What are the kernels used for the 3 examples in experiments? How sensitive is the experimental result to the kernel selection, regularization, and rank parameter?

4. Can the authors provide the learned dynamics of the 3 examples compared to the true dynamics, the invariant measure, and time scales?
limitations: limitations: Limitations have been discussed.
paper_summary: paper_summary: In this paper, the authors consider the problem of learning the infinitesimal generator a Stochastic Diffusion Process (SDP). Compared to existing approaches such as [1] they tackle the unbounded nature of the generator by introducing a novel statistical framework which is based on the Dirichlet form associated with the SDP. In this framework, the authors estimate the resolvent of the generator which can be approximated with finite-rank operators. They consider a regularized and rank-truncated version of the regression loss. Theorem 1 provides a way to compute the eigenvalues and eigenvectors of the estimated operator. Spectral learning bounds are derived in Theorem 2. The theoretical results of the paper are completed with three examples: learning the overdamped Langevin generator for a one-dimensional four well potential, the overdamped Langevin generator for a Muller Brown potential and finally for a Cox-Ingersoll-Ross process (CIR process).

[1] Cabannes et al. (2023) -- The Galerkin method beats Graph-Based Approaches for Spectral Algorithms
strength_weakness: strength_weakness: Before expressing my main concerns regarding this paper, I want to emphasize that I'm not an expert in this domain and therefore my understanding of the main competitors and methods used in the paper is lacking. 

* My main concern is regarding the experiments. I know this is a theoretical paper but a new statistical framework and a new learning procedure should be clearly validated against competing methods like [1,2,3]. Even though some qualitative conclusions are drawn I would have liked to see a more extensive study (for different choices of hyperparameters for both the target Langevin diffusion and the method introduced by the authors). 

* Indeed, there are little details in the paper about how to choose the crucial parameters of the algorithm like $\mu$, $r$ and $\gamma$. Establishing the robustness of the procedure regarding these hyperparameters seem crucial to validate the methodology.

* As mentioned before, even though I appreciate the rigour of the paper it is notation heavy. Having to constantly rely on the (massive) table of page 13 to remember the notation hindered my understanding of the paper. Is there a way to either reduce the notational load or to incorporate a lightweight version of Table 2 in the main paper? This would greatly help the reader.

* All assumptions are commented except (KE). What are the main limitations of this assumption? I am not familiar with the domain so a light explanation in the main text would also be appreciated.

* Is the assumption of l.188 that there exists a Dirichlet operator realistic? It is true for Langevin operator and CIR (which already covers quite a lot of ground) but apart from these models? Does this impose anything on the diffusion?

Some minor remarks:

* l.178 "Importantly, this energy functional can be empirically estimated from data sampled from π, whenever full knowledge, that is drift and diffusion coefficients of the SDE (1), or partial knowledge," -- This sentence is not really clear to me. How can one leverage the drift and diffusion coefficient here?

* l.215 typo (spurious new line) 

[1] Cabannes et al. (2023) -- The Galerkin method beats Graph-Based Approaches for Spectral Algorithms

[2] Hou et al. (2023) -- Sparse learning of dynamical systems in RKHS: An operator-theoretic approach

[3] Pillaud-Vivien et al. (2023) -- Kernelized Diffusion Maps
questions: questions: See weaknesses
limitations: limitations: Limitations are addressed in Section 7 ("Conclusion")
paper_summary: paper_summary: The paper considers a time-homogeneous Stochastic Differential Equation (SDE) with known diffusion part and known or unknown drift. The problem is to find properties of this equation from known data, in particular, to find a (low-rank) representation of its infinitesimal generator (IG). For this purpose, the resolvent of this operator is considered and the reduced-rank estimator in reproducing kernel Hilbert spaces (RKHS) together with the energy-based risk metric are used. Accuracy estimates for the found approximation to IG and time-complexity are given. Several model numerical experiments are conducted for proof-of-concept.
strength_weakness: strength_weakness: - The paper uses quite complex concepts and may be difficult to understand in a 9 page format.

- There are not enough practical examples of the application of these results to machine learning


Minor 

The x-axis in Figure 1 (a) is not signed

L120 comma is not needed

L124 space after dot is missing

L206: colon missing at the end of the line

L215 A new line  is redundant (`this work we focus on the case when`)

L275 space after dot is missing

L915 "anan" -> "and"
questions: questions: - Can your method be extended to the case when neither drift nor diffusion coefficients are known?
limitations: limitations: Limitation are described in the paper. The paper is mostly theoretical.
paper_summary: paper_summary: This paper proposes a relevant and sound approach for learning self-adjoint SDE generators via operator learning techniques.

The paper includes a compactification, a novel prior knowledge inclusion, and first-of-its-kind statistical learning guarantees that extend the known ones from discrete Markov processes.
strength_weakness: strength_weakness: 1. Unclear extensions to entirely data-driven regimes (no partial knowledge) and requiring sampling from an invariant distribution.

2. Experiments only consider toy examples that validate theory but do not showcase the practicality of the approach for more general cases.

3. I am not sure that contribution **4)** is well supported by the current presentation, primarily due to hard-to-parse figures. A revised comparison and legends with sharper plots could make the results much easier to understand.
questions: questions: 1. *Are Dirichlet forms commonly deduced from the diffusion part of the SDE?* It would help assess the ease of knowing these forms beforehand for readers unfamiliar with the notion.

2. *Is it reasonable to think you could deal with a bias of the Dirichlet form using imperfect diffusion knowledge?*

3. *Could you highlight what practical gains/losses are introduced w.r.t. discrete-time setting?* I would guess you need less data due to diffusion priors and not requiring as small of a discretization. On the other hand, IG regression complexity seems higher for high-dimensional systems.

4. *Can a sample complexity gain be recognized using prior knowledge in this form for your theoretical results?* Could the performance gain w.r.t. TO learning be recognized using the derived bounds?

5. *How arbitrary is the choice of $\mu$?*

6. *How does the time-sampled data come into play, and does the sampling rate of the data influence any aspect of the approach?* Does sampling the invariant distribution + Dirichlet form allow you to avoid time-derivative or "1-step" dynamics observations (e.g., compared to [A])? Based on the current writing, it is not immediately clear to me (but I might have missed something).

7. *How do you compare (using full knowledge) to classical numerical methods using known drift and diffusion? Would your sample efficiency be better than classical FEM for a given precision?* You mention how, for known operators (via drift and diffusion knowledge), there would be no spuriousness in estimating eigenpairs. This is perhaps of separate interest.


[A] - Meng, Y., Zhou, R., Ornik, M., & Liu, J. (2024). Koopman-Based Learning of Infinitesimal Generators without Operator Logarithm. http://arxiv.org/abs/2403.15688
limitations: limitations: 1. *Possibly strong prior knowledge requirements.*
2. *Limited to self-adjoint generators.* 
3. *Experiments are on 1D toy examples.*

Comments: title: Acknowledgement to the reviewer
Acknowledgement to the reviewer
Acknowledgement to the reviewer
Rebuttal Reply for Submission18874 by Reviewer D2P3
comment: We would like once more to thank the reviewer for their comments, which inspired us to make additional steps and improve our work.We are happy that our rebuttal was helpful and we commit to incorporate it in the revised manuscript.
We would like once more to thank the reviewer for their comments, which inspired us to make additional steps and improve our work.We are happy that our rebuttal was helpful and we commit to incorporate it in the revised manuscript.
I thank the authors for their explanations. After reading the global response and the discussion with other reviewers, I believe that in general my concerns were addressed. I raise my score to 7 "Accept".
Thank you for your thorough  answer and rebuttal. I have increased my score.
We would like once more  to sincerely thank the reviewer for their deep questions, which inspired us to make these additional steps and improve our work. 

We are happy that our rebuttal is appreciated and we commit to incorporate it in the revised manuscript.
I thank the authors for the extensive rebuttal in the general reply and comprehensively addressing my review comments. Also, it is great to see the authors put in the effort to deliver adjusted theory as well as run new experiments.

Many important aspects of the paper got clearer: limits of FEM and TO approaches, fully-data driven regime, energy functional rationale and shift parameter.
Accordingly, I have increased my overall score.

I expect this rebuttal content to be incorporated as part of a final submission.