Learning Kolmogorov Backward Operators for Spectral Analysis and Forecasting: A Data-Driven Variational Galerkin Framework
Authors: Yuanchao Xu, Claire Valva, Ben Gao, Isao Ishikawa
Organizations: Department of Mathematics, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan · Department of Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA 91125, USA · Laboratoire Hubert Curien, UMR 5516, CNRS, IOGS, INRIA, Université Jean Monnet, 18 Rue du Professeur Benoît Lauras, 42000 Saint-Étienne, France
For a time-homogeneous Markov process, the fixed-lag solution operator (Tτf)(x)=E[f(Xt+τ)∣Xt=x] of the Kolmogorov backward equation evolves observables linearly and deterministically, even when the dynamics are stochastic. We approximate this operator using a two-stage method where we first approximate this operator directly from data, without knowledge of the governing equation (e.g., Fourier Neural Operator), and then compress the Galerkin projection onto a finite-dimensional observable space that is fixed in advance or selected variationally. The error analysis traces the learning error of the estimator, the sampling error of the Galerkin inner products, and the invariance and reconstruction properties of the observable space to the computed spectra and forecasts. The sampling error enters only through a multiplicative factor. Near an isolated spectral cluster, the invariance defect of a subspace controls its distance from the target invariant subspace. The forecast error is controlled by the one-step error of the reduced model and the projection error of the observable; when the measure is invariant and the reduced model preserves constants and is stable on mean-zero observables, the bound is uniform in the horizon. For deterministic dynamics, a learning rate for Fourier neural operators (FNOs) gives one explicit bound for the learning error. Numerical experiments illustrate spectral approximation and forecasting, including long-rollout stability and a comparison with extended dynamic mode decomposition (EDMD) when fewer transition pairs than basis functions are observed.
Figures & tables
Figure 1: OU spectral approximation at τ=0.15 . Top: analytic generator rates and the BSO (learned dictionary) estimates, labeled “trained basis” in the legend, each matched to the nearest analytic rate. Bottom: postprocessed BSO (learned dictionary) eigenfunctions and the analytic Hermite modes on ∣x∣≤3 .
Figure 2: Triple-well spectral approximation at the lag 3τ=0.9 . Top: leading BSO (learned dictionary) multipliers and finer-grid finite-difference references. Bottom: the stationary, left–right, and upper–lower modes, paired by one-to-one matching; the panel labels give the continuous-time rates λ=(3τ)−1logμ corresponding to the multipliers μ of the top panel ( section 2.2 ).
Figure 3: Pooled relative L2 forecast error of BSO (fixed dictionary) and the autoregressive FNO over 64 held-out origins for Lorenz–63 (top) and Lorenz–96 (bottom). The short-term Lorenz–63 panel uses a logarithmic vertical axis; the remaining panels use linear axes with different vertical scales.
Figure 4: Global SST evaluation. Top: forecast skill of BSO (fixed dictionary), FNO, SFNO, and the three baselines over 96 held-out 2020–2022 origins, showing area-weighted RMSE and MAE through day 30 and RMSE and centered anomaly ACC through day 90. The SFNO curve ends at the first lead time for which scores are not available for all 96 origins. Bottom: area-weighted SST RMSE in the 540-day numerical-stability test over 48 origins, shown over the full and a zoomed logarithmic range.
Method
Error E40EEG
Reduction with BSO
BSO
0.749
–
Direct FNO
0.874
14.2%
EDMD-DL
0.941
20.4%
Kernel EDMD
0.945
20.7%
PCA–VAR
0.963
22.2%
Table 1: EEG experiment: relative L2 error accumulated through step 40 of the four-second rollout, E40EEG , and the relative reduction achieved by BSO with respect to each comparator, computed from unrounded errors.
Figure 5: Cumulative relative L2 error of BSO (learned dictionary) and the comparators for the four-second autonomous rollout on the S029R02 resting EEG record. Each step advances 0.1 seconds.
Figure 6: One-step prediction error on held-out transition pairs (relative L2 ) versus the number nop of observed transition pairs for Lorenz–96 with 40 sites and m=40 (vertical dashed line). BSO (fixed dictionary) is available for every nop shown. The dashed and dotted curves are the ridge-regularized and minimum-norm m=40 EDMD fits, which coincide to within 10−4 ; the solid EDMD curve uses the first min{nop,m} dictionary coordinates. In the legend, n stands for nop .
Appendix figures & tables1 asset
Supplementary material from the paper’s appendix.
Appendix
Field
Source (variable)
Level / depth
Unit
SST (target)
WB2 (ERA5 SST)
surface
K
Zonal wind stress
WB2 (from 10 m winds)
surface
N m -2
Meridional wind stress
WB2 (from 10 m winds)
surface
N m -2
Sea-surface height
GODAS ( sshg )
surface
m
Potential temperature
GODAS ( pottmp )
55 m
∘ C
Potential temperature
GODAS ( pottmp )
105 m
∘ C
Appendix
Table 2: Input fields of the SST experiment on the 32×64 grid. WB2: WeatherBench 2 archive of ERA5. Units are those of the fields as used, before standardization. Only SST is verified.
We establish approximation and learning guarantees for Fourier neural operators (FNOs) applied to time-T solution operators of dissipative evolution equations. The analysis builds on the premise that FNOs can efficiently approximate and learn solution operators whenever these operators admit stable and accurate spectral discretizations. To formalize this idea, we introduce classes of evolution operators defined through spectral methods and derive FNO approximation bounds and polynomial sample complexity guarantees for these classes. For equations with polynomial nonlinearities, the learning rates depend primarily on the smoothness of the input space and the dimension of the physical domain. Our results hold uniformly over broad families of dissipative equations, rather than for a single fixed PDE, and apply in particular to the Navier--Stokes, Allen--Cahn, and Cahn--Hilliard equations. For equations with non-polynomial smooth nonlinearities, we prove that polynomial sample complexity still holds with rates that now additionally depend on the smoothness of the nonlinear terms and the dissipation strength. Overall, we connect classical spectral approximation theory with modern operator learning and explain when FNOs can learn nonlinear evolution operators efficiently.
Nisha Chandramoorthy, Daniel Sanz-Alonso, Nathan Waniorek
Koopman operator theory provides a global linear representation of nonlinear dynamics and underpins many data-driven methods. In practice, however, finite-dimensional feature spaces induced by a user-chosen dictionary are rarely invariant, so closure failures and projection errors lead to spurious eigenvalues, misleading Koopman modes, and overconfident forecasts. This paper addresses a central validation problem in data-driven Koopman methods: how to quantify invariance and projection errors for an arbitrary feature space using only snapshot data, and how to use these diagnostics to produce actionable guarantees and guide dictionary refinement? A unified a posteriori methodology is developed for certifying when a Koopman approximation is trustworthy and improving it when it is not. Koopman invariance is quantified using principal angles between a subspace and its Koopman image, yielding principal observables and a principal angle decomposition (PAD), a dynamics-informed alternative to SVD truncation with significantly improved performance. Multi-step error bounds are derived for Koopman and Perron--Frobenius mode decompositions, including RKHS-based pointwise guarantees, and are complemented by Gaussian process expected error surrogates. The resulting toolbox enables validated spectral analysis, certified forecasting, and principled dictionary and kernel learning, demonstrated on chaotic and high-dimensional benchmarks and real-world datasets, including cavity flow and the Pluto--Charon system.
Gustav Conradie, Nicolas Boullé, Jean-Christophe Loiseau +2
Centre for Mathematical Sciences, University of Cambridge, UK · Department of Mathematics, Imperial College London, UK · Laboratoire DynFluid, Arts et Métiers Institute of Technology, France +1
Machine learning has opened new frontiers in purely data-driven forecasting of dynamical systems; the resulting methods are showing considerable promise in some application domains. However, in contrast to model-driven algorithms, analysis of these data-driven methods is poorly developed. In this paper we address this issue, developing a theory to underpin data-driven methods to solve forecasting problems. The theoretical framework comprises two key components: (i) establishing the existence of a mapping defining the forecasting problem; (ii) showing that it is possible to approximate this mapping by supervised learning techniques, to arbitrary accuracy. We work in the continuous-time and noise-free setting as this defines a simple canonical setting in which to establish an underpinning theory; the continuous-time setting leads to an operator formulation of the problem in which the underlying forecasting map is defined on pathspace. As a byproduct of our framework we show that analogous existence and approximation results may be obtained for the smoothing problem of recovering unobserved components of dynamical systems, from observed components, using supervised learning. The theoretical results, for both forecasting and smoothing, are illustrated with experiments studying the Lorenz 63, Lorenz 96 and Kuramoto-Sivashinsky dynamical systems.
Edoardo Calvello, Elizabeth Carlson, Nikola Kovachki +2
California Institute of Technology, Pasadena, CA · Oregon State University, Corvallis, OR · NVIDIA Corporation, Santa Clara, CA +1