Thank you for your insightful feedback. Due to space constraints, we address the key issues and questions briefly, committing to incorporate all suggestions in our revision. Further details are available upon request.

Empirical evidence:

For our additional results, please see reply to nvbC .

Presentation

Following your suggestions, we'll emphasize motivations in Sec 3, address unboundedness effects, expand the app. on kernel-based operator learning for TOs and IGs, and add a notation table. If accepted, we'll use the extra page to enhance transitions between challenging concepts.

Questions

• In stochastic process theory, TOs form the Markov semigroup associated with the process, acting as direct image functors in measurable spaces. When a measurable func pushes a measure forward, TO describes density transformations. Neural-network learning methods for TO relate to FNO but are distinct concepts, (Nakao and Mezic, 2020).

• TOs are linear (evolution) operators defined on func spaces. We can say that $A_t$ acts as the pushforward of a function $f$ under the transition probability defined by $p_t$, making them crucial for understanding $(X_t)\_t$ dynamics. On the other hand, IGs, defined via Eq 25, are linear differential operators related to the eqs of motion via backward Kolmogorov eq. $\partial_t \mathbb{E}[f(X_t)|X_0=\cdot ]=L f$ for observable $f$ and forward Kolmogorov eq. (Fokker-Planck) $\partial_t q_t =L^* q\_t$ for probability density $q_t$ of $X_t$ (* being adjoint). Indeed, Eq 7 is obtained by applying Itô's lemma to $f(X_t)$ where $X$ is the solution of Eq 6.

• Learning the spectral decomposition of TO/IG is equivalent to learning the solution of an associated SDE, so we can predict the evolution of distributions, see reply to CrN8. Working with IG allows us to do this reliably in continuous time.

• Aș explained in Sec 5, the data is the observed trajectory of states.

• Galerkin projection methods (empirical risk minimization) suffer from spectral pollution when applied to unbounded operators, (Kato 2012, Kostic et al. 2024). When the process evolves in a bounded domain $\mathcal X$, the IG are still typically unbounded, i.e. the IG doesn’t map all funcs $\mathcal{X}\to\mathbb{R}$ (observables of the process) in the $\mathcal{L}^2_\pi$ space to other funcs in that space - it’s only defined on a proper subset (domain of IG).

Q. for Authors

[1] $(p_t)_t$ is a family of transition probability densities $p_t(x,y)dy=P(X_t \in dy|X_0=x)$ for all $x,y\in\mathcal{X}$. In this sense, $p_t(x,y)dy$ quantifies the probability that the process $(X_t)_t$ is in the infinitesimal set $dy$ at time $t$, given that it started at $x$.

[2-3] See above.

[4-5] While the action of TO can be defined for any measurable func $\mathcal{X}\to\mathbb{R}$, to characterize dynamics via semigroup, one needs the space of such funcs to be invariant under the action of TOs, and a typical choice is $\mathcal{L}_\pi^2$. Note that invariance of measure $\pi$ (both marginals of $P[X_t,X_0]$ are $\pi$), doesn't mean that TOs (conditional laws $P[X_t|X_0]$ ) are trivial.

[6] In fact, for all $t$, $A_t$ is an operator that turns a func $f$ into a func $A_t f$ s.t. for any $x\in\mathcal{X}$, $A_t f(x)=\mathbb{E}[ f(X_t )|X_0 =x]$. While $L$ is an operator that maps a func $f\in dom(L)=\mathcal{W}\_\pi^{1,2}\subset \mathcal{L}\_\pi^2$ to a func $Lf \in\mathcal{L}\_\pi^2$, meaning that $\int_{\mathcal{X}} |L f(x)|^2\pi(dx)<\infty$.

[7-8] In the RKHS theory it is standard to assume $\mathcal{H}\subset\mathcal{L}^2_{\pi}(X)$, $\pi$ being the probability of data samples (e.g. Gaussian kernel), and the main issue in the learning bounds is the difference between $\mathcal{H}$ (chosen) and $\mathcal{L}^2_{\pi}(X)$ (unknown) norms, e.g. (Steinwart&Christmann,2008).

[9] For conciseness, the defs of $A$, $L$, and $\pi$ in the Ornstein-Uhlenbeck case are in App A.1. If space permits we can include them in the main text.

[10-11] The funcs $\phi$ are not unknown but are provided when choosing the RKHS, e.g. $\phi(x)=e^{-\|x-\cdot\|^2/(2\sigma^2)}$ for Gaussian kernel.

[12-15] In Eqs 9-10, $X_0$ is not the SDE's initial cond. but any random variable with the invariant measure, see reply to Q1 of PJzG for extra clarity. Note that $\phi(X_t),\psi(X_0), G^*\phi(X_0) \in \mathcal{H}$ are funcs, not scalars.

[14 & 16-17] $\mathbb{E} [\psi(X_0)|X_0=\cdot]:\mathcal{X}\to\mathcal{H}$ is the Riesz representation of lin. operator $R_\mu$ in RKHS $\mathcal{H}$, and our learning objective. It is the regression func related to the excess risk, see Prop A.1, where the risk is given via the target feature $\psi(X_0)$. Since we cannot observe it precisely (due to the integral form) in Eqs 11-12 we show how to approximate it via quadrature.

[18-19] $A_{t_j}$ is TO for time-delay $t_j$ (see Eq 2), $h$ was used as a generic func in RKHS $\mathcal{H}$ and $v$ as a vector in $\mathbb{R}^N$.

We appreciate the reviewer’s insightful evaluation and valuable comments. Below, due to space limitations, we briefly address the highlighted weaknesses and respond to the reviewer’s questions, committing to incorporating all feedback in our revision. If needed, we can elaborate further on each point.

Additional empirical evidence:

Following the reviewer's suggestion, we have expanded our empirical evaluation to include IG baselines and support the claims summarized in Table 1. For a brief discussion on the results of our additional experiments, please see our reply to reviewer nvbC .

Presentation

To help the reader, we will include a table of major notations at the beginning of the appendix, and expand our appendix with more context on the kernel based operator learning and derivation of the algorithms. Further, if accepted, we will use the extra page to include smoother transitions between challenging concepts. Finally, concerning the discussion of the Alanine-Dipeptide experiment, since there is no ground truth for this system, the only thing we could show is that the result aligns with the state-of-the-art expert knowledge in molecular dynamics, (Wehmeyer & Noe, 2018). In the revised manuscript we will discuss more on the interpretation of this figure.

Questions

We thank the reviewer for their questions, which suggest potentially interesting extensions of our framework.

• Partially observed systems: The major problem in partially observed systems, in the context of TO/IG kernel-based learning, lies in the fact that one might lose the universality of representation, and hence encounter the difficulty to properly predict distributions. One possible way to overcome this is to rely on Takens-type theorems (deterministic and/or stochastic) (Sauer et al. ”Embedology”, 1991; ) Koltai and Kunde, “A Koopman–Takens theorem”, 2024), that essentially guarantee, under certain assumptions, that by augmenting the current measurements with the past ones we decrease the loss of information and improve estimation.

• Fractional Brownian motion: Consider the SDE $dX_t = b(X_t) dt + \sigma(X_t) dB_{t}^H$ where $B^H =(B\_{t}^H )\_{t}$ is a fractional Brownian motion. When $H \neq 1/2$, the fractional Brownian motion (fBm) is no longer a semimartingale, which means Itô calculus does not apply directly. Although the equation can still be given a meaning—via the Stieltjes integral for $H > 1/2$ or the Skorokhod integral for $H \in (1/4,1/2)$—the resulting process is no longer Markovian. Our generator-based approach no longer applies in this setting. However, as you suggest, we could adapt our method to estimate the generator for a Markovian approximation associated with the autonomous version of the fractional SDE introduced in [1]: $dX_t = b_\theta (X_t) dt + \sigma_\theta (X_t) dB_{t}^H.$ This process can be approximated using a finite linear combination of Ornstein-Uhlenbeck processes, $\hat{B}^H(t)$, leading to the so-called Markov-Approximate fractional SDE (MA-fBMSDE): $dX_t = b_{\theta} (X_t) dt + \sigma_{\theta} (X_t) d\hat{B}\_{t}^H,$ where $d\hat{B}\_{t}^H = \sum_{k=1}^{K} \omega_k dY_t^k, \quad \text{with} \quad dY_t^k = -\gamma_k Y_t^k dt + dW_t.$ By applying Proposition 2, the process $X_t$ can be augmented with a finite number of Markov processes $Y_t^k$ (which approximate $B^H$), forming a higher-dimensional state variable: $Z_t = (X_t, Y_{t}^1, \dots, Y_{t}^K) \in \mathbb{R}^{D(K+1)}.$ The process $Z = (Z\_{t})\_{t}$ is Markovian and can be described by an ordinary SDE: $dZ_t = h_{\theta} (Z_t) dt + \Sigma_{\theta} (Z_t) dW_t.$ The infinitesimal generator is given by: $\mathcal{L} f(z) =\left( b (x) - \sigma(x) \sum_{k=1}^{K} \omega_k \gamma_k y^k \right) \frac{\partial f}{\partial x} + \sum_{k=1}^{K} (-\gamma_k y^k) \frac{\partial f}{\partial y^k}+ \frac{1}{2} \sum_{i,j} (\Sigma_{\theta} \Sigma_{\theta}^\top)_{i,j} \frac{\partial^2 f}{\partial z_i \partial z_j}.$ This satisfies our assumptions, allowing us to apply our procedure.

Forecasting

Given IG’s spectral decomposition, Eq (8) provides directly the solution since it enables forecasting of full state distributions beyond just the mean (e.g., $f$ as an indicator function). Hence, we can use our method to forecast $\mathbb{E}[h(X_t)\vert X_0 = x] \approx \sum_{i\in[r]} e^{\hat \lambda_i} \langle \hat{g}\_{i}, h \rangle\_{\mathcal H} \hat{f}\_{i}(x)$, for $h\in\mathcal{H}$, noting that $\hat{g}\_{i}, h \rangle\_{\mathcal H}$ can be computed on the training set via kernel trick, see (Kostic et al 2022). Further, note that this formula extends to all $\mathcal{L}^2_\pi$ functions, at the price of an additional projection error, and that we can, hence, predict evolution of distributions, see e.g. (Klus et al., 2019) and (Kostic at al. "Consistent long-term forecasting of ergodic dynamical systems", 2024). If useful, we can include an empirical example in the main body or appendix.

We appreciate the reviewer’s insightful evaluation and valuable comments. Below, due to space limitations, we briefly address the highlighted weaknesses and respond to the reviewer’s questions, committing to incorporating all feedback in our revision. If needed, we can elaborate further on each point.

Additional empirical evidence (results can be found here ):

While we still feel that our work is mainly theoretical, following the reviewer's suggestion, we have expanded our empirical evaluation to include IG baselines and support the claims summarized in Table 1. In particular:

1D Langevin dynamics experiment:

We extend the original experiment by comparing it to two IG baselines that use prior knowledge of IG:

• Galerkin projection estimate (Hou et al., 2023), and
• Energy based Dirichlet form regression (Kostic et al., 2024).

In Figure 4 (see the above link) we plot true eigenvalues as black vertical lines, while the estimated eigenvalues are plotted as magenta (our dual method) and red (Dirichlet form regression) dashed lines. Further, since the Galerkin projection estimate has numerous spurious eigenvalues, we plot their empirical distribution in the form of histogram, stressing out that expert knowledge is needed for this estimator to extract good estimates from the spurious ones. Finally, we report that our results are comparable to the SOTA (Kostic et al., 2024) estimator, despite not using the explicit IG knowledge, while being one order of magnitude faster to train.

Non-normal sectorial IG:

We further conduct an experiment using the 2D Ornstein-Uhlenbeck process with non-symmetric drift, estimating the leading nontrivial complex-conjugate eigenvalue pair. We compare our primal method based on random Fourier features with the corresponding TO, noting that Kostic et al. (2024) is inapplicable in this setting and that Galerkin projection estimators (Hou et al., 2023) are inefficient, since, as observed above, they typically result in over 50 eigenvalues in the zone of interest. In Figure 5 (see the link above) we show comparison for 10 random trials using 1000 random features, two different time discretizations $\Delta t=0.01$ and $\Delta t = 0.001$, and corresponding sample sizes $n=10^4$ and $n=10^5$. Note that, to obtain small errors, sample sizes are much higher than in the self-adjoint case (Langevin). This is due to non-normality of the generator, as predicted by our theoretical analysis, see Appendix C.5, where one can see that our method consistently has estimates in the tighter $\varepsilon$-pseudospectrum of the operator.

Presentation

As suggested, we will include a table of major notations at the beginning of the appendix.

We appreciate the reviewer’s insightful evaluation and valuable comments. Below, due to space limitations, we briefly address the highlighted weaknesses and respond to the reviewer’s questions, committing to incorporating all feedback in our revision. If needed, we can elaborate further on each point.

Additional empirical evidence:

Following the reviewer's suggestion, we have expanded our empirical evaluation to include IG baselines and support the claims summarized in Table 1. For a brief discussion on the results of our additional experiments, please see our reply to reviewer nvbC .

Algorithms

We will expand Appendix B with detailed derivations, including a discussion of the sampling operator $\hat{S}$. Briefly, by the definition of $\hat{S}$ for $h \in \mathcal{H}$, we derive $\hat{S}^* \hat{S} h = n^{-1/2} \hat{S}^* [h(x_1) \dots h(x_n)]^{\top} = n^{-1} \sum_{i \in [n]} \phi(x_i) h(x_i) = [n^{-1} \sum_{i \in [n]} \phi(x_i) \otimes \phi(x_i)] h$, concluding that the empirical covariance can be written as $\hat{S}^* \hat{S}$.

Transient dynamics

As discussed in Appendix A.3, when $L^* L \neq LL^*$, eigenvalues alone may not fully explain the evolution of linear dynamics. The norm of the resolvent relates to transient growth in stable systems via the pseudo-spectrum (Trefethen & Embree, 2020). We will expand this discussion in the revision.

TO as $\Delta t \to 0$

The core idea is that $A_{\Delta t} = e^{\Delta t L }\to I$, which implies that the spectral gap vanishes, thereby rendering the eigenfunction bounds vacuous.

Minor issues

We will correct all typos. Additionally, we remark: The presence of $G^*$ in the risk can be understood from two perspectives. One follows from TO’s risk formulation in the space of observables, which, due to the kernel trick, translates into vector-valued regression (Kostic et al., 2022). The other, as the reviewer notes, considers the Perron-Frobenius operator, (Klus et al., 2019), on probability distributions, linking the estimator to the MMD metric induced by a characteristic kernel. We will expand Section A.2 to introduce regression-based operator learning for both TO and IG methods. $\rho_t$ denotes the distribution of $(X_s, X_{s+t})$ and $\mathrm{gap}_i$ is introduced in the last line of Theorem 6.3. A table of major notations will be added at the beginning of the appendix.

Questions

[Q1] When $X_0$’s distribution is not invariant, variance analysis requires adapting the method of blocks for mixing with Bernstein inequalities in Hilbert spaces for independent (but non-identically distributed) variables. This complicates the effective dimension, as the covariance operator w.r.t. the invariant distribution is replaced by one w.r.t. the ergodic mean. While technically feasible, this adds complexity, so we opted for a simpler approach to highlight key contributions.

[Q2] Knowing IG implies knowledge of the SDE, requiring a numerical solution for process realizations. In contrast, given IG’s spectral decomposition, Equation (8) provides directly the solution, since it enables forecasting of full state distributions beyond just the mean (e.g., $f$ as an indicator function). See reply to CrN8.

[Q3] To our knowledge, our technique for controlling Bochner integral approximation errors for linear operators is novel, possibly extending beyond TO/IG methods. We will emphasize this in the contributions paragraph.

[Q4] As the reviewer rightly points out, the Laplace transform is a well-known analytical tool in the study of continuous operator semigroups. In the context of deterministic dynamical systems, it was used by Mohr & Mezic, (2014) to investigate the spectral decomposition of the Koopman operator (TO for deterministic dynamical systems). Their work underscores its potential but leaves efficient numerical methods as an open problem. In contrast, we consider SDEs and not only design numerical methods based on the Laplace transform, but also develop statistical learning theory, providing sharp bounds on spectral estimation. Additionally, we remark that the results of (Bevanda et al. 2023) are also limited to deterministic dynamical systems. There, the Laplace transform served more as inspiration than an explicit component of the method, as the authors built an RKHS via finite-horizon integration of kernel features over trajectories to formulate Koopman operator regression and study learning bounds.