1 Introduction↩︎

Many problems in scientific computing and engineering disciplines involve modeling and prediction of dynamical systems, with applications including weather [1]–[3], environmental and ecological science [4]–[7], biology [8]–[10], fluid dynamics [11]–[14], finance [15], [16], and traffic [17]–[19]. However, many physically relevant problems remain challenging due to high-dimensionality, complex or chaotic dynamics, lack of known underlying dynamics, and noisy or low-fidelity observational data.

Purely data-driven methods have emerged as a strong option for learning dynamical systems [20]–[24]. Such methods circumvent the need to form dynamical equations and typically seek to represent unknown dynamics with a high-fidelity reduced-order model. A variety of popular approaches have proven to be competitive in this context. Dynamic mode decomposition (DMD) and variants leverage a Koopman framework that lifts finite-dimensional nonlinear data to an infinite-dimensional linear representation, often leading to a simplified low-dimensional surrogate model [24]–[29]. Sparse identification of nonlinear dynamical systems (SINDy) and variants discover dynamical equations from data by choosing a suitable, often sparse, linear combination of dictionary functions [30]–[33]. Neural ordinary differential equations (NODEs) learning underlying dynamics by training a neural network to represent a continuous-time vector field [34]–[39]. Various machine learning techniques such as Long Short-Term Memory [40]–[43], reservoir computing [44]–[47], and autoencoders [48] have also gained prominence. Kernel-based approaches such as kernel ridge regression (KRR) mitigate the curse of dimensionality with the so-called “kernel trick” [49]–[53]. KRR is especially attractive because it does not require a dictionary of functions, and is straightforward to implement. Recent work has shown that KRR outperforms multiple baseline methods in data-driven dynamical system learning and forecasting problems over a wide range of data sets [52]. While these approaches typically perform well when applied to clean data, their performance often degrades significantly in the presence of measurement or observational noise, especially when the underlying dynamics are chaotic.

Several methods have been proposed to mitigate lack of robustness in the presence of noise, including data assimilation and filtering methods [54]–[56], and Gaussian process approaches [57]–[59]. The approach developed in [55], coined RAFDA, has proven to be competitive for low-dimensional systems, but does not scale well to high-dimensions. A different approach is to combine a learning strategy with a weak formulation, which involves integrating data residuals over a family of test functions. Classical (strong) approaches enforce pointwise consistency between given observational data, and typically perform poorly for noisy data because pointwise errors are magnified by erroneous fluctuations. In contrast, weak approaches relax pointwise loss strategies in favor of orthogonality constraints between residuals and test functions. Loss functions involving weak formulations have been observed to perform more robustly in the presence of noise, and have been incorporated successfully into both the SINDy framework [60]–[63] and the NODE framework [36]. Both [36] and [61] point out that a weak formulation acts as a filter for noisy data, suggesting a mechanism for its observed noise robustness.

In the present manuscript, we will build on these previous insights to show explicitly that the weak formulation acts to filter noisy data. We will consider a family of test functions that arise from uniform translations of a generating function, and will interpret the resulting weak formulation as an orthogonal projection procedure which filters by projecting noisy data onto the subspace spanned by these test functions. Moreover, we provide a brief bias-variance error decomposition of signal filtering via the weak form. In particular, we provide an exact expression for the variance, i.e., the error that arises due to noise corruption. This interpretation connects to well-known results in information theory and signal processing [64]–[68]. Shannon’s sampling theorem, which states that an ideal band-limited function can be perfectly reconstructed with suitably sampled data, may be interpreted as an orthogonal projection onto a subspace of band-limited functions [69]. Subsequent work has extended this orthogonal projection analysis to classes of functions arising from integer shifts of generating functions [64]. While the error analysis that we provide is not novel, we include it to justify the use of a weak formulation and to clarify the mechanism by which it provides noise robustness.

Motivated by the simplicity and recent success of kernel-based learning methods and the noise robustness of weak formulations, we propose Weak-form Kernel Ridge Regression (WKRR) as a noise robust, data-driven learning framework. We remark that WKRR does not require a suite of dictionary functions, in contrast to Weak SINDy approaches [60]–[63] that critically need the underlying functions or vector fields to be spanned by the dictionary functions. Moreover, WKRR does not require extensive hyperparameter tuning as many machine learning frameworks, such as NODE, require [36]. We will show numerically that the forecasting horizon of WKRR is similar to that of a strong KRR formulation, but much less computationally expensive. We will also demonstrate the effectiveness of WKRR with different choices of kernel function, including the standard Gaussian kernel and the Diffusion Maps (DM) kernel, which has recently received attention for its excellent performance across a wide range of chaotic and experimental datasets [52]. The success of WKRR with various kernel functions broadens the method’s applicability and suggests that the practitioner may select kernel functions that favor speed or accuracy as the situation warrants.

The remainder of the manuscript is organized as follows. In §2, we review the classical KRR method and the Gaussian and DM kernel functions. In §3, we review the weak formulation, explain how it acts as a filter, and provide a brief bias-variance error decomposition for the filtering procedure. The main contribution of this work is provided in §4, where we propose WKRR as a noise-robust learning method for dynamical systems. We propose a validation procedure to select appropriate model hyperparameters, and summarize the steps needed to implement WKRR in practice. We apply WKRR to several numerical examples in §5, including chaotic baseline systems in up to 64 dimensions, and real-world turbulent fluid data made available by [14] as part of a Community Challenge. We conclude with a brief discussion in §6.

2 Kernel Ridge Regression Review↩︎

In this section, we review the classical kernel ridge regression (KRR) framework for learning solution operators of dynamical systems. We will refer to this approach as the “strong approach” throughout the manuscript.

2.1 Kernel Ridge Regression↩︎

Suppose we are given noisy data \(\mathbf{u}(t_i)\equiv\mathbf{u}_i = (u_i^{(1)},\dots,u_i^{(n)})\in \mathbb{R}^n\), \(i=1,\dots,N\), sampled at times \(t_i = (i-1)\cdot \Delta t\). We assume that the data is of the form \[\label{eq:32data} \mathbf{u}_i = \mathbf{x}_i + \sigma \boldsymbol{\xi}_i,\tag{1}\] where \(\mathbf{x}_i\) denotes clean data generated by an autonomous dynamical system of the form \(\mathbf{x}'= \mathbf{f}(\mathbf{x})\), and \(\boldsymbol{\xi}_i\) is i.i.d. Gaussian with zero mean and covariance matrix \(\boldsymbol{\Sigma} = \operatorname{diag}(\eta_1^2,\dots,\eta_n^2)\). We define \(\eta_\ell\) to be the root-mean-square (RMS) of the \(\ell\)th component of the clean signal \[\eta_\ell = \sqrt{\frac{1}{N} \sum_{i=1}^N \left|x_i^{(\ell)}\right|^2}.\] The parameter \(\sigma \geq 0\) represents the signal-to-noise ratio and will be varied in our numerical experiments.

In this section, we assume that \(\sigma = 0\) so that the given data is not corrupted by noise. Let \(\mathcal{M}\subset \mathbb{R}^n\) denote the forward invariant set of the underlying system. Define the flow map \(\boldsymbol{\Phi}_\tau = (\Phi_{\tau}^{(1)},\dots,\Phi_{\tau}^{(n)}): \mathbb{R}^n \to \mathcal{M}\), which advances the system state over a time increment \(\tau>0\), i.e., \(\boldsymbol{\Phi}_\tau(\mathbf{x}(t)) = \mathbf{x}(t+\tau)\).

We will use KRR to learn solution operators with either a direct-connection or a skip-connection scheme. Both frameworks are useful in practice; direct-connection often outperforms skip-connection when the underlying dynamics are stiff. Otherwise, the skip-connection scheme typically yields superior performance [52]. We employ both frameworks in our numerical examples. We now describe each learning framework in turn.

2.1.1 Direct-Connection Scheme↩︎

Let \(k(\cdot,\cdot;\epsilon):\mathbb{R}^n\times \mathbb{R}^n \to \mathbb{R}\) denote a scalar-valued kernel function, where \(\epsilon >0\) denotes a tunable bandwidth parameter. The key idea of the direct-connection scheme is to approximate each component of the flow map directly, \[\label{eq:32KRR32approx} \Phi^{(\ell)}_\tau(\mathbf{x}) \approx \sum_{i=1}^N k(\mathbf{x},\mathbf{u}_i;\epsilon)\alpha_i^{(\ell)},\tag{2}\] for reconstruction coefficients \(\{\alpha_{i}^{(\ell)}\}_{\ell=1,\ldots, n}\) to be determined. To proceed, we collect the available data in snapshot pairs \(S=(\mathbf{X},\mathbf{Y})\). The input matrix \(\mathbf{X}\in \mathbb{R}^{ (N-1)\times n}\) is defined by \[\label{eq:32X32matrix32clean} \mathbf{X}= \begin{bmatrix} \mathbf{u}_1 & \mathbf{u}_2 & \dots & \mathbf{u}_{N-1} \end{bmatrix}^\top,\tag{3}\] while the output matrix \(\mathbf{Y}\in \mathbb{R}^{(N-1)\times n}\) is defined by \[\label{eq:32direct-connection32matrix} \mathbf{Y}= \begin{bmatrix} \mathbf{u}_2 & \mathbf{u}_3 & \dots & \mathbf{u}_{N} \end{bmatrix}^\top.\tag{4}\]

Let \(K(\cdot,\mathbf{X};\epsilon) = [k(\cdot,\mathbf{u}_1;\epsilon),\dots,k(\cdot,\mathbf{u}_{N-1};\epsilon)]:\mathbb{R}^n \to \mathbb{R}^{N-1}\) be the kernel function evaluated over the input data \(\mathbf{X}\). In this notation, our goal is to model the flow map \[\label{eq:32kernel32direct-connection} \mathbf{u}_{i+1} = \boldsymbol{\Phi}_{\Delta t}(\mathbf{u}_i) \approx K(\mathbf{u}_i,\mathbf{X};\epsilon)\boldsymbol{\alpha},\tag{5}\] where \(\boldsymbol{\alpha}\in \mathbb{R}^{(N-1)\times n}\) collects the reconstruction coefficients. With the goal of expressing 5 in matrix notation, define the Gram matrix \(\mathbf{K}(\epsilon)\in \mathbb{R}^{(N-1)\times (N-1)}\) whose \((i,j)\) entry is \(k(\mathbf{u}_i,\mathbf{u}_j;\epsilon)\) for \(i,j=1,\dots,N-1\). We then arrive at a linear least squares problem which enforces pointwise consistency between the snapshot data pairs \[\label{eq:32least32squares32problem} \mathbf{Y}= \mathbf{K}(\epsilon) \boldsymbol{\alpha}.\tag{6}\] We remark that 6 is often ill-conditioned. In practice, a solution \(\hat{\boldsymbol{\alpha}}\) is determined by solving the related KRR problem \[\label{eq:32KRR32problem} \hat{\boldsymbol{\alpha}} = \arg\min_{\boldsymbol{\alpha}} \left(\|\mathbf{K}(\epsilon)\boldsymbol{\alpha} - \mathbf{Y}\|_F^2 + \lambda \boldsymbol{\alpha}^\top \mathbf{K}(\epsilon)\boldsymbol{\alpha} \right),\tag{7}\] where \(\lambda\) is a user-specified regularization parameter. The problem 7 admits a unique solution \[\hat{\boldsymbol{\alpha}} = \left(\mathbf{K}(\epsilon) + \lambda \mathbf{I}\right)^{-1}\mathbf{Y},\] where \(\mathbf{I}\) denotes the \((N-1)\times (N-1)\) identity matrix. Once \(\hat{\boldsymbol{\alpha}}\) has been determined, out-of-sample model forecasting can be achieved via \[\label{eq:32forecasting32direct} \tilde{\mathbf{u}}_{i+1} = K(\tilde{\mathbf{u}}_i,\mathbf{X};\epsilon)\hat{\boldsymbol{\alpha}},\tag{8}\] where \(\tilde{\mathbf{u}}_i\) is the kernel approximation to the given data.

2.1.2 Skip-Connection Scheme↩︎

Instead of directly modeling the flow map, a skip-connection scheme models each component of the residual map \[\Phi^{(\ell)}_\tau(\mathbf{x}) - \mathbf{x}^{(\ell)} \approx \sum_{i=1}^N k(\mathbf{x},\mathbf{u}_i;\epsilon)\alpha_i^{(\ell)}.\] We define \(\mathbf{X}\) according to 3 as before, but now define the output matrix to be \[\label{eq:32skip-connection32matrix} \mathbf{Y}= \begin{bmatrix} \Delta\mathbf{u}_1 & \Delta\mathbf{u}_2 & \dots & \Delta\mathbf{u}_{N-1} \end{bmatrix}^\top, \quad \Delta \mathbf{u}_i = \mathbf{u}_{i+1} - \mathbf{u}_i.\tag{9}\] In this setup, our goal is to model \[\label{eq:32kernel32skip-connection} \mathbf{u}_{i+1} - \mathbf{u}_i = \boldsymbol{\Phi}_{\Delta t}(\mathbf{u}_i) - \mathbf{u}_i \approx K(\mathbf{u}_i,\mathbf{X};\epsilon)\boldsymbol{\alpha}.\tag{10}\] In the same notation as above, we arrive at a least squares problem 6 , but with \(\mathbf{Y}\) now defined according to 9 . Solving for the coefficients \(\hat{\boldsymbol{\alpha}}\) proceeds as above, and forecasting is performed according to \[\label{eq:32forecasting32skip} \tilde{\mathbf{u}}_{i+1} = \tilde{\mathbf{u}}_i + K(\tilde{\mathbf{u}}_i,\mathbf{X};\epsilon)\hat{\boldsymbol{\alpha}}.\tag{11}\]

Notice that both the direct- and skip-connection schemes require the user to specify a kernel function, a bandwidth parameter \(\epsilon\), and a regularization parameter \(\lambda\). We discuss appropriate choices for the kernel function below. We defer the discussion of the tunable parameters \(\epsilon\) and \(\lambda\) to §4.2, where we extend a previously developed validation strategy [52] to handle noisy data using a weak formulation.

2.2 Choice of Kernel Function↩︎

The choice of kernel function in the KRR approach described above is crucial to the success of the method. While any kernel function can be used, a specific choice of kernel may be more, or less, appropriate for a given problem.

Suppose we are given a dataset \(\mathbf{U}= \{\mathbf{u}_1,\dots,\mathbf{u}_N\}\) consisting of observations sampled from the forward invariant set \(\mathcal{M}\subset \mathbb{R}^n\). Unless otherwise stated, in this work we will consider a standard Gaussian kernel \(k_{\text{RBF}}:\mathcal{M}\times \mathcal{M}\to \mathbb{R}\), \[\label{eq:32Gaussian32kernel} k_{\text{RBF}}(\mathbf{x},\mathbf{y};\epsilon) = \exp\left(-\frac{\|\mathbf{x}-\mathbf{y}\|_2^2}{4\epsilon} \right),\tag{12}\] for any \(\mathbf{x},\mathbf{y}\in \mathcal{M}\). The parameter \(\epsilon>0\) is a scalar-valued bandwidth parameter that should be determined by the practitioner. We describe an approach to choose \(\epsilon\) in §4.2 where we extend the KRR approach to noisy data.

We also consider the Diffusion Maps (DM) kernel, a data-driven kernel function based on the Diffusion Maps algorithm [70] that has recently been employed with success over a wide range of datasets [52]. In particular, recent work has shown that the DM kernel can exhibit superior forecasting performance compared to the standard Gaussian kernel when applied to clean data, especially when the dimension of the invariant set \(\mathcal{M}\) is much lower than its ambient dimension, \(n\) [52].

We will compare the Gaussian and DM kernels over two baseline examples in §5 when only noisy data is available. We will show numerically that both kernels have nearly identical performance for a low-dimensional chaotic system, but that the DM kernel exhibits superior performance on a high-dimensional chaotic system with low intrinsic dimension. These results suggest that systems with low intrinsic dimension or special geometry may benefit from the DM kernel, even when observational data is corrupted by noise. Otherwise, the Gaussian kernel may be more appropriate due to lower computational complexity and ease of implementation.

We now describe the numerical implementation of the DM kernel. Following [52], we begin with the standard Gaussian RBF kernel 12 . We then normalize 12 in two stages to arrive at the DM kernel, \(k_{\text{DM},N}:\mathcal{M}\times\mathcal{M} \to \mathbb{R}\), by first writing \[k_{s,N}(\mathbf{x}, \mathbf{y};\epsilon) = \frac{k_{\text{RBF}}(\mathbf{x}, \mathbf{y};\epsilon)}{s_{N}(\mathbf{x};\epsilon) s_{N}(\mathbf{y};\epsilon)}, \qquad s_{N}(\mathbf{x};\epsilon) = \frac{1}{N}\sum_{j = 1}^N k_{\text{RBF}}(\mathbf{x}, \mathbf{u}_j;\epsilon),\] and then writing \[\label{eq:32our32discrete32DM32kernel} k_{\text{DM},N}(\mathbf{x},\mathbf{y};\epsilon) = \frac{{k}_{s,N}(\mathbf{x},\mathbf{y};\epsilon)}{\sqrt{{q}_{N}(\mathbf{x};\epsilon){q}_{N}(\mathbf{y};\epsilon)}},\qquad {q}_{N}(\mathbf{x};\epsilon) = \frac{1}{N}\sum_{j = 1}^N k_{s,N}(\mathbf{x}, \mathbf{u}_j;\epsilon),\tag{13}\] for any \(\mathbf{x},\mathbf{y}\in \mathcal{M}\). Note that the DM kernel is not the same as the DM matrix [70], \[\label{eq:32DM32matrix} P_{N} (\mathbf{x},\mathbf{y};\epsilon) = \frac{ k_{s,N}(\mathbf{x}, \mathbf{y};\epsilon)}{\sum_{j=1}^N k_{s,N}(\mathbf{x}, \mathbf{u}_j;\epsilon)},\tag{14}\] that approximates a Markov transition kernel of a reversible Markov chain on \(\mathcal{M}\). Observe that 14 defines a Markov matrix because the normalizing denominator term here is a summation, not an average as in 13 .

One can verify that \[\label{DMsymmetrickernel} \sqrt{N{q}_{N}(\mathbf{x};\epsilon)} \frac{P_{N} (\mathbf{x},\mathbf{y};\epsilon)}{\sqrt{N{q}_{N}(\mathbf{y};\epsilon)}} = \frac{{k}_{s,N}(\mathbf{x}, \mathbf{y};\epsilon)}{\sqrt{N{q}_{N}(\mathbf{x};\epsilon) N{q}_{N}(\mathbf{y};\epsilon)}}= \frac{1}{N}k_{\text{DM},N}(\mathbf{x},\mathbf{y};\epsilon),\tag{15}\] which suggests that the Gram matrix \(\mathbf{P}(\epsilon)\) corresponding to \(P_{N}\) is diagonally conjugate to the normalized Gram matrix \(\mathbf{K}(\epsilon)/N\) corresponding to \(k_{\text{DM},N}/N\). In particular, the eigenvalues of \(\mathbf{P}(\epsilon)\) and \(\mathbf{K}(\epsilon)/N\) are the same. We refer to [71] for a theoretical discussion of the DM kernel in conjunction with KRR.

3 Filtering Noise Through a Weak Formulation↩︎

The key learning mechanism underlying classical KRR is to enforce pointwise consistency across snapshot pairs of data. While this strong formulation is often appropriate for clean data, its performance often noticeably worsens in the presence of noise because the regression may fit erroneous fluctuations instead of underlying dynamics.

The weak formulation is a popular approach to address this lack of noise robustness. Weak approaches integrate noisy data over a family of test functions, relaxing pointwise consistency in favor of orthogonality constraints. While such approaches have recently been observed to offer increased robustness to noise [36], [60]–[63], [72], the connection to filtering was recently pointed out in [36] and [61]. Building on this intuition, the purpose of this section is to provide a brief overview of the filtering mechanism behind the weak formulation.

3.1 Problem Statement and Notation↩︎

A classical solution to the differential equation \(\mathbf{x}'(t) = \mathbf{f}(\mathbf{x}(t))\) is a solution \(\mathbf{x}(t)\) which satisfies \[\mathbf{R}(\mathbf{x}(t)) \equiv \mathbf{x}'(t) - \mathbf{f}(\mathbf{x}(t)) = \boldsymbol{0}.\] In contrast, a weak formulation involves integrating the residual over smooth test functions \(\varphi(t)\) \[\int_{\mathbb{R}} \text{d}t \; \mathbf{R}(\mathbf{x}(t))\varphi(t) = \boldsymbol{0},\] which may be interpreted as an \(L^2\) inner product enforcing the orthogonality of the residual with the test functions. This integral formulation, with appropriate choice of compactly supported test functions, has been observed to increase robustness to noise [36], [60]–[63], [72].

In the following, we will study weak formulations applied to an arbitrary noisy signal \(\mathbf{u}(t)\), \[\label{eq:32integral} \int_{\mathbb{R}} \text{d}t \; \mathbf{u}(t)\varphi_j(t) = \boldsymbol{0},\tag{16}\] for an appropriate family of test functions \(\{\varphi_j\}\). While we will ultimately employ the weak formulation in conjunction with KRR, the analysis in this section does not depend on a specific learning method.

Let \(\varphi(t)\) be a smooth test function with compact support on \([0,L]\), and define \(\varphi_j(t) = \varphi(t-jh)\) to be a family of uniformly translated test functions with compact support on \([L_{j}, L_{j+1}]\). The parameter \(h\) determines the distance between the adjacent test functions, and together with \(L\), describes the extent to which adjacent test functions overlap. For a dataset \(\mathbf{u}_i\equiv\mathbf{u}(t_i)\) of finite length and fixed \(h\), the indices \(j\) for which the support of \(\varphi_j(t)\) lies in the sampling time is restricted. If the signal consists of \(N\) points sampled uniformly with step \(\Delta t\), we adhere to the convention that \(j\in [1,2,\dots,k^*]\), where \[k^* = \left\lfloor \frac{N\Delta t-L}{h} \right\rfloor.\]

In practice, the integral 16 must be computed with quadrature over the available data \[\label{eq:32quadrature32integral32friend} c_j^{(\ell)} = \left\langle \varphi_j(t), u^{(\ell)}(t) \right\rangle = \int_{L_j}^{L_{j+1}} \text{d}t\; \varphi_j(t) u^{(\ell)}(t) \approx \sum_{i=1}^{N} w_i \varphi_j(t_i) u_i^{(\ell)}, \quad \ell = 1,\dots,n, \quad j=1,\dots,k^*,\tag{17}\] where \(w_i\) are appropriate quadrature weights. We will find it convenient to express these quadrature approximations as a matrix product. To that end, collect the translated test functions into the rows of a matrix \(\boldsymbol{\Psi}\in \mathbb{R}^{k^* \times N}\) whose \((j,i)\) entry is \(\varphi_j(t_i)\). Note that \(\boldsymbol{\Psi}\) is often sparse due to the compact support of the test functions. Further define \(\mathbf{W}=\operatorname{diag}(w_1,\dots,w_N)\in \mathbb{R}^{N\times N}\) to be a diagonal matrix of quadrature weights. We assume that \(\mathbf{W}\) is invertible throughout the manuscript. Note that \(\boldsymbol{\Psi}\mathbf{W}\) is the matrix whose rows are test functions scaled by the quadrature weights.

With this notation, we may discretize 17 as a linear system \[\label{eq:32more32coefficient32nonsense} \mathbf{C}= \boldsymbol{\Psi} \mathbf{W}\mathbf{U},\tag{18}\] where \(\mathbf{C}\in \mathbb{R}^{k^*\times n}\) is the matrix of inner product coefficients \(c_j^{(\ell)}\), and \(\mathbf{U}= [\mathbf{u}_1,\dots,\mathbf{u}_N]^\top\in \mathbb{R}^{N\times n}\) collects the data at available sample times. Notice that in the typical case that \(k^*<N\), the coefficients \(\mathbf{C}\) represent a compressed representation of the original signal \(\mathbf{U}\).

3.2 The Weak Formulation as a Filter↩︎

To see the weak formulation as a filter, we consider the process of reconstructing a signal from a set of coefficients \(\mathbf{C}\). Let \({\hat{\mathbf{u}}}(t)=(\hat{u}^{(1)}(t),\dots,\hat{u}^{(n)}(t))\) denote the reconstructed signal, which we now express as a linear combination of the available test functions \[\label{eq:32expansion32nonsense} \hat{u}^{(\ell)}(t) = \sum_{j=1}^{k^*} a_j^{(\ell)}\varphi_j(t), \quad \ell=1,\dots,n,\tag{19}\] where \(a_j^{(\ell)}\) are unknown coefficients to be determined.

In the event that the \(\varphi_j(t)\) are orthonormal, taking inner products on both sides of 19 is sufficient to isolate the coefficients. Without orthogonality, we may still recover the coefficients by requiring that the inner product of the reconstruction with test functions \(\varphi_q(t)\) agree with the given coefficients, \[\label{eq:32reconstruction32nonsense} \mathbf{C}_{q,\ell}=\langle\hat{u}^{(\ell)}(t),\varphi_q(t) \rangle = \sum_{j=1}^{k^*} a_j^{(\ell)}\langle\varphi_j(t),\varphi_q(t)\rangle.\tag{20}\]

We will find it useful to express 20 as a matrix equation. Let \(\mathbf{A}\in \mathbb{R}^{k^* \times n}\) be the collection of unknown coefficients whose \((j,\ell)\) entry is \(\alpha_j^{(\ell)}\), and let \(\mathbf{G}= \boldsymbol{\Psi} \mathbf{W}\boldsymbol{\Psi}^\top \in \mathbb{R}^{k^*\times k^*}\) be the matrix whose \((q,j)\) entry is a quadrature approximation of the inner product \(\langle \varphi_q(t), \varphi_j(t)\rangle\). Under the assumption that the test functions are linearly independent and overlapping, \(\mathbf{G}\) is invertible. We can then express 20 as the linear system \(\mathbf{C}= \mathbf{G}\mathbf{A}\). It follows that \[\label{eq:32coefficient32nonsense} \mathbf{A}= \mathbf{G}^{-1}\mathbf{C}.\tag{21}\]

Let \(\hat{\mathbf{U}}= [\hat{\mathbf{u}}_1,\dots,\hat{\mathbf{u}}_N]^\top \in \mathbb{R}^{N\times n}\) denote the matrix collecting the reconstructed signal at the available discrete sample times. We can now express 19 in matrix notation to arrive at an expression for \(\hat{\mathbf{U}}\), \[\hat{\mathbf{U}} = \boldsymbol{\Psi}^\top \mathbf{A}= \boldsymbol{\Psi}^\top \mathbf{G}^{-1}\mathbf{C}= \boldsymbol{\Psi}^\top \mathbf{G}^{-1} \boldsymbol{\Psi} \mathbf{W}\mathbf{U}= \boldsymbol{\Psi}^\top ( \boldsymbol{\Psi}\mathbf{W}\boldsymbol{\Psi}^\top)^{-1}\boldsymbol{\Psi}\mathbf{W}\mathbf{U},\] where the second equality follows from 21 and the third follows from 18 . Define \(\mathbf{P}\equiv \boldsymbol{\Psi}^\top ( \boldsymbol{\Psi}\mathbf{W}\boldsymbol{\Psi}^\top)^{-1}\boldsymbol{\Psi}\mathbf{W}\in \mathbb{R}^{N\times N}\). Note that \(\mathbf{P}\) is a projection matrix because \(\mathbf{P}^2=\mathbf{P}\). Moreover, we see that \[\langle \mathbf{P}\mathbf{u},\mathbf{v}\rangle_{\mathbf{W}}:= (\mathbf{P}\mathbf{u})^\top\mathbf{W}\mathbf{v}= \mathbf{u}^\top \mathbf{P}^\top \mathbf{W}\mathbf{v}= \mathbf{u}^\top\mathbf{W}\boldsymbol{\Psi}^\top \mathbf{G}^{-1}\boldsymbol{\Psi} \mathbf{W}\mathbf{v}= \mathbf{u}^\top \mathbf{W}(\mathbf{P}\mathbf{v})= \langle \mathbf{u},\mathbf{P}\mathbf{v}\rangle_{\mathbf{W}},\] for any \(\mathbf{u},\mathbf{v}\in \mathbb{R}^N\), which implies that \(\mathbf{P}\) is self-adjoint with respect to the inner-product weighted by the quadrature weight matrix \(\mathbf{W}\). In particular, we have \[\label{eq:32self-adjoint} \mathbf{P}^\top\mathbf{W}= \mathbf{W}\mathbf{P}.\tag{22}\] These observations imply that \(\mathbf{P}\) is an orthogonal projection under the inner product \(\|\mathbf{P}\|_W^2 = \text{tr}(\mathbf{P}^\top \mathbf{W}\mathbf{P})\).

The formulation above suggests that the weak-form reconstruction procedure can be interpreted as an orthogonal projection onto the span of a family of test functions, \[\label{eq:32weak32reconstruction32formula} \hat{\mathbf{U}} = \mathbf{P}\mathbf{U},\tag{23}\] which allows weak formulations to be interpreted in light of classical results in information theory and signal processing, such as Shannon’s sampling theorem and related results which use the language of orthogonal projection to study filtering [64], [65], [69].

3.2.1 Choice of Test Function↩︎

In this section, we consider a specific choice of test function that was found to be successful in the literature [61], \[\label{eq:32poly32test32function} \varphi(t) = \begin{cases} C_p (t-a)^p(b-t)^p, & a<t<b \\ 0, & \text{otherwise} \end{cases}, \qquad C_p = \left(\frac{2}{L}\right)^{2p},\tag{24}\] where \(L=b-a\) is the support length and \(p\in \mathbb{Z}^+\) is a parameter which specifies the degree of the polynomial. Note that \(\varphi(t)\) is smooth and compactly supported. We will use 24 in all of our numerical examples in §5.

Recall that we construct the matrix \(\boldsymbol{\Psi}\mathbf{W}\) by generating a family of linearly independent test functions, \(\varphi_j(t) = \varphi(t - jh)\), where \(h>0\) dictates the overlap between test functions. Consequently, the polynomial test functions depend on three key parameters: (i) the polynomial degree \(p\), (ii) the support length \(L\), and (iii) the overlap \(h\). Appropriate values for the triple \((p,L,h)\) may vary significantly across datasets and noise levels. In practice, an appropriate choice of these parameter values should balance speed and reconstruction accuracy.

There are many approaches to empirically measure the accuracy of the weak-form reconstruction. If clean or high-quality reference data is available, one could select parameters such that the reconstruction of the noisy data most closely aligns with the clean reference data under a suitable metric, e.g., RMSE. If no clean reference data is available, a common approach is to define a score metric that can be evaluated from only the observed data. Such metrics commonly measure whiteness, variance ratio, or smoothness ratio [36], [73]–[77]. We refer to [36] for a discussion and numerical implementation of this approach. Meanwhile, the number of test functions \(k^*\) is inversely proportional to the overlap parameter \(h\). We will show in §4 that taking \(h\) large has the potential to greatly reduce computational complexity requirements.

In this work, we select polynomial parameters that are observed to filter given noisy data across a range of error metrics. As a qualitative example, we consider here the Lorenz-63 system, discussed in more detail in §5, with 10% noise corruption. Figure 1 visualizes a typical filtering result for the \(x\)-component with \((p,L,h)=(5,\,0.5,\,0.1)\). We remark that filtering with appropriately chosen polynomials can provide competitive or superior results to standard approaches, such as wavelet-based filtering. We refer to Appendix 7 for a quantitative comparison of polynomial and wavelet filtering across a range of error metrics, noise levels, and model systems. The success of the filtering procedure justifies this choice of test function and explains the noise robustness of the weak formulation.

3.3 Bias-Variance Decomposition↩︎

To analyze the error in the general reconstruction 23 , let \(\boldsymbol{\mathcal{E}}=\hat{\mathbf{U}}-\mathbf{U}_{\text{clean}}\) denote the difference between the underlying clean data and the reconstruction of the noisy signal. It follows from the decomposition \[\hat{\mathbf{U}} = \mathbf{P}\mathbf{U}= \mathbf{P}(\mathbf{U}_{\text{clean}} + \sigma \boldsymbol{\xi}),\] where \(\boldsymbol{\xi} = [\boldsymbol{\xi}_1,\dots,\boldsymbol{\xi}_N]^\top \in \mathbb{R}^{N\times n}\) is the collection of noise vectors, that we may write \[\label{eq:32bias32variance32decomposition} \mathbb{E}\left[\|\boldsymbol{\mathcal{E}}\|_W^2\right] = \sum_{\ell=1}^n B_\ell^2 + V_\ell^2, \quad \quad B_\ell^2 = \|(\mathbf{P}- \mathbf{I})\mathbf{U}_{\text{clean}}^{(\ell)}\|_W^2, \quad V_\ell^2 = \sigma^2\eta_\ell^2\text{tr}(\mathbf{W}\mathbf{P}),\tag{25}\] which gives a bias-variance decomposition of the reconstruction error. The terms \(B_\ell^2\) correspond to the “bias” of the reconstruction, which describe the fidelity with which a clean signal can be reconstructed. It measures how well the clean signal can be approximated in the span of the test functions. The terms \(V_\ell^2\) correspond to the “variance” of the signal, which describes additional error arising from noise corruption.

The variance term admits significant simplification in the case of uniform quadrature, \(\mathbf{W}= c\mathbf{I}\) for \(c\in \mathbb{R}\), \[\label{eq:32trace} \text{tr}(\mathbf{W}\mathbf{P}) = c\cdot \text{tr}(\mathbf{P}) = c\cdot\text{rank}(\mathbf{P}) = c k^*,\tag{26}\] which holds because \(\mathbf{P}\) is a projection matrix. Equation 26 implies that filtering error due to noise is governed by the number of test functions, \(k^*\), and is independent of the form of the test function. Notice that the bias term is similarly scaled by \(c\) in this case, so that a specific choice of this constant does not artificially inflate either the bias or variance error relative to the other.

In contrast, the bias term generically depends on the choice of test function and is generally difficult to analyze. Performing computations in Fourier space has the potential to simplify the form of the bias. One can show that in the ideal case of periodic data and uniform quadrature, \(\mathbf{W}=\mathbf{I}\), the bias and variance have the form \[\label{eq:32bias32variance32simplification} B_\ell^2 = \left\|(\mathbf{H}- \mathbf{I})\mathcal{F}\left[\mathbf{U}_{\text{clean}}^{(\ell)}\right]\right\|_F^2, \quad V_\ell^2 = \sigma^2 \eta_\ell^2 k^*.\tag{27}\] where \(k^* = N\Delta t/ h\) (due to the periodicity of the data) and \(\mathcal{F}[\cdot]\) denotes the Fourier transform. The matrix \(\mathbf{H}\in \mathbb{R}^{N\times N}\) is a Fourier representation of \(\mathbf{P}\), \[\label{eq:32H} \mathbf{H}_{mn} = \begin{cases} \frac{\hat{\varphi}_m \overline{\hat{\varphi}_n}}{\sum_{\ell=0}^{h-1} |\hat{\varphi}_{r+\ell k^*}|^2}, & m = n \;(\text{mod } k^*), \\ 0, & \text{otherwise}, \end{cases}\tag{28}\] where \(r=m\text{ mod } k^*\) and \(\hat{\varphi}_m\) is the \(m\)th component of the Fourier transform of the untranslated test function. Note that \(\mathbf{H}\) is sparse with nonzero banded diagonal entries. The entries on the off-diagonal corresponding to aliasing artifacts. We remark that 28 holds for a general test function. However, further simplification for a specific test function is often difficult or impossible. We provide a brief derivation of 28 in Appendix 8. Finally, we remark that approximation theory may be used to bound the bias under certain conditions in terms of the sampling step of the data. We refer to [64], [65], [78], [79] for further theoretical discussion of the bias reconstruction error.

In the remainder of this section, we numerically study the bias and variance in 27 with the specific choice of test function 24 . We assume the ideal case of periodic data and uniform quadrature, \(\mathbf{W}= \mathbf{I}\). Similarly to [61], we observe that in certain parameter regimes the bias term may act as an ideal low-pass filter and is well approximated by the tail of the Fourier modes of the underlying clean signal, \[\label{eq:32tail} B_\ell^2 \approx \sum_{m>k^*} \left|\mathcal{F}[\mathbf{U}_{\text{clean}}^{(\ell)}]\right|^2.\tag{29}\] To perform numerical experiments, we consider a noisy dataset, \(\mathbf{U}= [\mathbf{u}_1,\dots,\mathbf{u}_N]\), and its clean counterpart, \(\mathbf{U}_{\text{clean}}\), that have been generated from the Lorenz-63 system 39 (see §5). We take \(N=1000\), uniform timestep \(\Delta t = 0.01\), and \(\sigma = 0.1\). To proceed, we numerically observe the total error, \[\mathbb{E}\left[\|\boldsymbol{\mathcal{E}}\|_F^2\right] = \mathbb{E}\left[\|\mathbf{P}\mathbf{U}- \mathbf{U}_{\text{clean}}\|_F^2\right],\] and the bias error \[\label{eq:32oberved32bias} \mathbb{E}\left[\|\boldsymbol{\mathcal{E}}_{\text{bias}}\|_F^2\right] = \mathbb{E}\left[\|\mathbf{P}\mathbf{U}_{\text{clean}} - \mathbf{U}_{\text{clean}}\|_F^2\right].\tag{30}\] The observed variance error is then computed according to \[\label{eq:32observed32variance} \mathbb{E}\left[\|\boldsymbol{\mathcal{E}}_{\text{var}}\|_F^2\right] = \mathbb{E}\left[\|\boldsymbol{\mathcal{E}}\|_F^2 - \|\boldsymbol{\mathcal{E}}_{\text{bias}}\|_F^2\right].\tag{31}\] Above, the expectation is empirically approximated over 100 identical trials corrupted with independent noise.

Results are reported in Figure 2 for different test-function parameters. We demonstrate that the form of \(V_\ell^2\) in 27 captures the scaling of the observed variance, and observe that the bias may be well-approximated by 29 for large \(L\) and small \(h\).

4 Proposed Approach: Weak-form Kernel Ridge Regression↩︎

In this section, we propose Weak-form Kernel Ridge Regression (WKRR) as a noise robust extension of the strong KRR approach. We will refer to this framework as the “weak” approach throughout the manuscript.

4.1 Weak-form Kernel Ridge Regression↩︎

To understand WKRR, we follow the formulation of the strong case as in §2.1 up to the formulation of the least squares problem 6 , which enforces pointwise consistency across the given data. In the context of a weak formulation, we will interpret the least squares problem as a residual, \(\mathbf{E}= \mathbf{Y}- \mathbf{K}(\epsilon)\boldsymbol{\alpha}\) with \(\mathbf{Y}\) and \(\mathbf{K}(\epsilon)\) now depending on noisy data, and select \(\boldsymbol{\alpha}\) to enforce orthogonality of the residual \(\mathbf{E}\) against a family of test functions.

As in §3, let \(\varphi_j(t)\) denote a family of smooth test functions, and let \(\boldsymbol{\Psi}\mathbf{W}\in \mathbb{R}^{k^*\times (N-1)}\) be the matrix whose rows are the test functions scaled by quadrature weights. Notice here that \(\boldsymbol{\Psi}\mathbf{W}\) has \(N-1\) rows to be consistent with the snapshot data \(S\). The core idea behind WKRR is to enforce the constraint \(\boldsymbol{\Psi}\mathbf{W}\mathbf{E}= \boldsymbol{0}\). In a literal sense, this procedure approximates the integral of the residual over test functions by a quadrature approximation. In the sense described in §3, this procedure amounts to projecting the noisy data onto the span of the test functions, which we interpret as filtering the noisy data. Accordingly, we consider the modified least squares problem \[\label{eq:32least32squares32problem32weak} \mathbf{Y}_\mathbf{W}= \mathbf{K}_\mathbf{W}(\epsilon)\boldsymbol{\alpha},\tag{32}\] where \(\mathbf{Y}_\mathbf{W}= \boldsymbol{\Psi}\mathbf{W}\mathbf{Y}\in \mathbb{R}^{k^*\times n}\) and \(\mathbf{K}_\mathbf{W}(\epsilon) = \boldsymbol{\Psi}\mathbf{W}\mathbf{K}(\epsilon)\in \mathbb{R}^{k^*\times (N-1)}\).

Naively, classical KRR regularization would involve the term \(\lambda \boldsymbol{\alpha}^\top \mathbf{K}_\mathbf{W}(\epsilon) \boldsymbol{\alpha}\); however, this expression is not defined because the weak-form kernel matrix is no longer square. Instead, we regularize in the weak coefficient space in terms of the matrix \(\widetilde{\mathbf{K}}_\mathbf{W}(\epsilon) = \boldsymbol{\Psi}\mathbf{W}\mathbf{K}(\epsilon)(\boldsymbol{\Psi}\mathbf{W})^\top \in \mathbb{R}^{k^*\times k^*}\), leading to the problem \[\label{eq:32first32weak32KRR32equation} \hat{\boldsymbol{\alpha}}_\mathbf{W}= \arg\min_{\boldsymbol{\alpha}_\mathbf{W}} \left(\| \widetilde{\mathbf{K}}_\mathbf{W}(\epsilon)\boldsymbol{\alpha}_\mathbf{W}- \mathbf{Y}_\mathbf{W}\|_F^2 + \lambda \boldsymbol{\alpha}_\mathbf{W}^\top \widetilde{\mathbf{K}}_\mathbf{W}(\epsilon)\boldsymbol{\alpha}_\mathbf{W}\right).\tag{33}\] Note that \(\boldsymbol{\alpha}_\mathbf{W}\in \mathbb{R}^{k^*\times n}\) is a weak-form representation of the coefficients \(\boldsymbol{\alpha} \in \mathbb{R}^{(N-1)\times n}\) in the original coordinates. To obtain a solution in the original coordinates, we study 33 with the substitution \(\boldsymbol{\alpha} = (\boldsymbol{\Psi}\mathbf{W})^\top \boldsymbol{\alpha}_\mathbf{W}\), which leads to the equivalent problem \[\label{eq:32weak32H32KRR32equation} \hat{\boldsymbol{\alpha}} = \arg\min_{\boldsymbol{\alpha}} \left(\| \mathbf{K}_\mathbf{W}(\epsilon)\boldsymbol{\alpha} - \mathbf{Y}_\mathbf{W}\|_F^2 + \lambda \boldsymbol{\alpha}^\top \mathbf{K}(\epsilon)\boldsymbol{\alpha} \right),\tag{34}\] with unique solution given by \[\label{eq:32weak32H32alpha32hat} \hat{\boldsymbol{\alpha}} = (\boldsymbol{\Psi}\mathbf{W})^\top (\widetilde{\mathbf{K}}_\mathbf{W}(\epsilon) + \lambda \mathbf{I})^{-1}\mathbf{Y}_\mathbf{W}.\tag{35}\] Once \(\hat{\boldsymbol{\alpha}}\) is available, forecasting is performed via 8 or 11 . We emphasize that forecasting is performed in the original coordinates, not in the weak-coefficient space.

Note that the weak least squares formulation arises from the strong formulation left-multiplied by the matrix \(\boldsymbol{\Psi}\mathbf{W}\). Although this procedure appears very simple, it grants a surprising number of benefits. In addition to filtering, we will find that performing KRR in the weak-coefficient space reduces computational complexity and may increase model robustness compared to a strong implementation.

We remark that classically, a weak formulation involves integration-by-parts when integrating both sides of a differential equation over a family of test functions. Here, when learning discrete-time solution operators, we do not integrate by parts. Nevertheless, we continue to use the term “weak” due to the fact that these approaches enforce the differential equations in the integral or weak sense. We note that the approach described here can easily be extended to learn the vector field of a dynamical system, and emphasize that the interpretation of the weak form as a filter remains unaltered in a continuous-time context.

The success of the WKRR approach hinges on the choice of test function used to construct the matrix \(\boldsymbol{\Psi}\mathbf{W}\), and the hyperparameters \(\epsilon\) and \(\lambda\). In this work, we consider only the test function 24 described in §3.2.1. Developing a systematic or optimal framework for choosing test function parameters \((p,L,h)\) is beyond the scope of the present work. Instead, our objective is to demonstrate the robustness of the WKRR approach across a representative range of parameter values. In the numerical experiments (see §5), we select fixed representative values of \((p,L,h)\) and demonstrate their success in the WKRR framework. Additional diagnostics are provided in Appendix 7. The following section describes a validation procedure to choose the parameters \(\epsilon\) and \(\lambda\), which extends the methodology implemented in [52].

4.2 Validation↩︎

The success of WKRR depends on appropriate selection of the bandwidth and regularization parameters \((\epsilon,\lambda)\). In this section, we propose a validation procedure that we observe to be successful across a range of chaotic, high-dimensional, and experimental datasets; see §5.

Typically, one reserves a portion of given observation data for use in validation. This validation data is not seen during training, and is used as a form of ground truth to select appropriate model hyperparameters. In this section we let \(\mathbf{u}_{\text{train},i}\) for \(i=1,\dots,N\) and \(\mathbf{u}_{\text{val},i}\) for \(i=1,\dots,N_V\) denote noisy training and validation data, respectively. While training and validation can be performed on the raw data, in practice we observe that filtering prior to training and validation yields marginally better performance. We refer to §4.3 for a more detailed discussion of our pre-processing procedure. We denote \(\hat{\mathbf{u}}_{\text{train},i}\) and \(\hat{\mathbf{u}}_{\text{val},i}\) as filtered training and validation data, respectively. In the present discussion, we assume training and validation are performed on the filtered data. In our numerical experiments in §5, we consider training on both noisy and filtered data.

We assume that the matrix \(\boldsymbol{\Psi}\mathbf{W}\) is available. For a fixed \((\epsilon,\lambda)\) pair, let \(\mathbf{K}(\epsilon)\) be the kernel Gram matrix formed over the filtered training data \(\hat{\mathbf{u}}_{\text{train},i}\), and let \(\hat{\boldsymbol{\alpha}}\) be computed according to 35 . Denote \(\tilde{\mathbf{x}}_i\) as the corresponding kernel forecast over the filtered validation data \(\hat{\mathbf{u}}_{\text{val},i}\) via either 8 or 11 . We assume that \(\tilde{\mathbf{x}}_1 = \hat{\mathbf{u}}_{\text{val},1}\), i.e., that the forecast is initialized with the first step of known validation data. We now review two error metrics that we will use for validation in the numerical examples considered in this manuscript.

1. Valid Prediction Time (VPT): We consider the VPT metric for chaotic dynamical systems \[\text{VPT} = \Lambda\cdot t_\mathcal{I}, \quad \mathcal{I} = \max \{\mathfrak{I}\mid E_i \leq \gamma \text{ for all } 1\leq i \leq \mathfrak{I}\},\] where \(\Lambda\) is the maximal Lyapunov exponent of the system, \(\gamma\) is a user-defined threshold, and \(E_i\) is a normalized error metric \[\label{eq:32VPT} E_i = \sqrt{\frac{1}{n}\sum_{\ell=1}^n \left(\frac{\tilde{x}_i^{(\ell)}-\hat{u}_{\text{val},i}^{(\ell)}}{\zeta_\ell}\right)^2 },\tag{36}\] where \(\zeta_\ell\) is the standard deviation of the \(\ell\)th component of the filtered validation data.

   In words, VPT measures how long the normalized error between the kernel predication and the validation data stays below a given threshold \(\gamma\). The VPT metric is ideal for chaotic systems because the forecast is expected to diverge after a short time. Normalization by the Lyapunov exponent \(\Lambda\) ensures the metric is unitless and enables fair comparison with other chaotic systems with faster or slower dynamics.

2. Normalized Mean Square Error (NMSE): For all other dynamical systems considered in this work, we use the NMSE metric \[\label{eq:32NMSE} E_i = \frac{\sum_{\ell=1}^n |\tilde{x}_i^{(\ell)}-\hat{u}_{\text{val},i}^{(\ell)}|^2 }{\sum_{\ell=1}^n |\hat{u}_{\text{val},i}^{(\ell)}|^2}.\tag{37}\]

Our goal now is to select appropriate hyperparameters \((\epsilon,\lambda)\) for use in the weak KRR problem 34 . While systematic grid searches are possible, they are often computationally expensive. Instead, to save computation time, we follow a data-driven heuristic method, outlined in Appendix 9, to compute a reference bandwidth \(\epsilon^*\) and a reference regularization parameter \(\lambda^*\). These reference values are used to confine the parameter search to a grid \(H=[\epsilon^*\Delta \epsilon, \epsilon^*/\Delta \epsilon] \times [\lambda^*,\lambda^*/\Delta \lambda]\), where \(\Delta \epsilon = 10^{-4}\) and \(\Delta \lambda = 10^{-5}\).

Apart from the grid, the length of data over which forecasts are performed also plays a crucial role. In the case of clean data, it is often sufficient to validate over a few long trajectories that match or exceed the expected forecast horizon, as in [52]. However, such a procedure may yield unstable results when the available data is noisy; predictions may diverge much more rapidly due to fluctuations, causing the model to fit noise instead of underlying dynamics. To address this concern, here we propose to validate over many short segments, and select parameters that optimize average error over all segments. To reduce computational complexity, we first select parameters from a coarse grid \(H_c\). We then zoom in on the optimal parameter choice from the coarse grid, and form a second refined grid \(H_f\). Validation over all \((\epsilon,\lambda)\in H_f\) yields an optimal parameter pair, which is subsequently used to form the final model.

We now provide implementation details for the proposed validation procedure.

Coarse validation:

1. Form the coarse grid: Select reference parameters \((\epsilon^*,\lambda^*)\), form the coarse grid \(H_c=[\epsilon^*\Delta \epsilon, \epsilon^*/\Delta \epsilon]\times [\lambda^*,\lambda^*/\Delta \lambda]\), and divide \(H_c\) into an \(N_{\text{mesh}}\times N_{\text{mesh}}\) mesh with uniformly spaced gridpoints in log-space.

2. Divide validation data: Divide the given validation data into \(N_c\) (possibly overlapping) segments of length \(\nu_c\), with initial conditions chosen uniformly and randomly from the full validation dataset.

3. Learn and forecast: For each \((\epsilon,\lambda)\) pair in \(H_c\), learn the kernel solution operator via equation 35 , forecast over the \(N_c\) validation segments, and compute the average error over each segment.

4. Select optimal coarse parameters: Select the hyperparameter pair \((\epsilon_c,\lambda_c)\) that results in the best error over the \(N_{\text{mesh}}^2\) measurements.

Fine validation:

1. Form the fine grid: Define a zoom-in factor \(F\) and a fine grid \(H_f = [\epsilon_c/F,\epsilon_c\cdot F]\times [\lambda_c/F,\lambda_c\cdot F]\). Break \(H_f\) into an \(N_{\text{mesh}}\times N_{\text{mesh}}\) mesh with uniformly spaced gridpoints in log-space.

2. Divide validation data: To prevent overfitting, divide the given validation data into \(N_f\) (possibly overlapping) segments of length \(\nu_f\), with initial conditions chosen uniformly and randomly from the full validation dataset. This partition is independent from the previous partition in step 2.

3. Repeat steps 3-4 above using the fine mesh \(H_f\) to extract final hyperparameters \((\epsilon_f,\lambda_f)\).

Unless otherwise stated, we take \(N_{\text{mesh}}=20\), \(N_c=20\), \(N_f=30\), \(F=10^{3/4}\), and \(\nu\equiv\nu_c=\nu_f\). The lengths \(\nu_c\) and \(\nu_f\) depend on the available data and should be chosen on a problem specific basis.

Notice that in the above procedure, the strong formulation requires inverting an \(N\times N\) matrix for each parameter pair, which is computationally expensive when \(N\) is large. In contrast, WKRR involves inverting a \(k^*\times k^*\) matrix, which greatly reduces computational cost in the typical case that \(k^* \ll N\).

4.3 Signal Pre-Processing↩︎

Application of a weak formulation applied to KRR can be interpreted as filtering prior to regression. As we will show numerically in §5, some form of filtering is essential prior to model training in the presence of noisy data, whether the filtering is accomplished via explicit pre-processing or implicitly via a weak formulation. In our numerical experiments, we compare three approaches: (i) a strong formulation trained on pre-processed (filtered) data, (ii) a weak formulation trained on noisy data (but implicitly filtered via projection), and (iii) a weak formulation trained on pre-processed data (and further filtered implicitly via projection). In all cases, we validate on filtered data for the sake of consistency with the training data.

Because the weak formulation acts to filter data, it is natural to compare the proposed WKRR approach to a strong approach applied to filtered data. We now explain that WKRR grants several benefits beyond just filtering.

• Accuracy: We will show numerically in §5 across several examples that WKRR exhibits similar or superior forecasting performance over strong methods when trained on identical filtered data.

• Speed: WKRR leads to a compressed representation of observational data in the typical case of fewer test functions than samples, \(k^*<N\). Consequently, training and validation involve inverting matrices of size \(k^*\times k^*\), instead of \(N\times N\) required for the strong formulation.

• Validation robustness: In our numerical experiments, several strong models performed poorly during testing due to non-optimal validation. In contrast, WKRR models trained on identical data did not produce such outliers. This observation suggests that the weak-coefficient representation of the data may robustify the validation process, making specific hyperparameter choices less critical for success.

These observations lead us to incorporate data pre-processing in the WKRR pipeline. We remark that increased performance due to filtering followed by a weak formulation has been reported in the literature for other learning methodologies and is not unique to WKRR [62]. However, we emphasize that while pre-processing often enhances WKRR forecasting performance, WKRR without pre-processing still performs robustly across a range of examples.

4.4 Summary↩︎

We summarize the implementation of the proposed WKRR approach in Algorithm 3. We assume that noisy observational data \(\mathbf{u}_i\) of the form 1 is given, which may be divided into training data \(\mathbf{u}_{\text{train}}\) and validation data \(\mathbf{u}_{\text{val}}\).

5 Numerical Study↩︎

In this section, we apply the proposed WKRR method to two baseline systems: the 3-dimensional Lorenz-63 (L63) system, and a 64-dimensional discretization of the Kuramoto-Sivashinsky (KS) system. We compare the performance of the standard Gaussian kernel 12 with the DM kernel 13 . We also consider 15,000-dimensional experimental fluid data that corresponds to turbulent cavity flow (see Challenge 2.1 in [14]). When feasible, we compare WKRR with competitive baseline methods, such as strong KRR [52], RAFDA [55], and LSTM [40].

For the two baseline systems, we specify the signal-to-noise ratio of the given data 1 by selecting intensity \[\label{eq:32noise32intensity} \sigma \in \{0.01,\, 0.05,\, 0.1,\, 0.2\},\tag{38}\] which corresponds to 1%, 5%, 10%, or 20% noise corruption.

5.1 Lorenz-63↩︎

We first consider the Lorenz-63 (L63) dynamics, which are governed by the following equations [80] \[\label{eq:32Lorenz326332equations} \begin{align} x' &= \rho_1(y-x), \\ y' &= x(\rho_2-z) - y, \\ z' &= xy - \rho_3z. \end{align}\tag{39}\] Here, \(\rho_1 = 10\), \(\rho_2 = 28\), and \(\rho_3=8/3\), and 39 admits chaotic dynamics with Lyapunov exponent \(\Lambda \approx 0.91\) [36], [52], [81].

For data generation, we randomly generate an initial condition from the standard normal distribution. We integrate using MATLAB’s ode113 syntax with relative error tolerance \(10^{-10}\), absolute error tolerance \(10^{-12}\), and timestep \(\Delta t = 0.01\) to generate a trajectory consisting of \(N=68,000\) steps. Then, we corrupt the signal with a fixed intensity \(\sigma\) from 38 . The first \(10,000\) steps of the trajectory are pruned to eliminate transients and the resulting segment is divided into 2500 steps of noisy training data, 5500 steps of noisy validation data, and 50,000 of clean testing data. We save only the clean testing data as ground truth to gauge predictive performance. We repeat this procedure 100 times independently, which results in 100 segments consisting of training, validation, and testing data.

We then pre-process the training and validation data separately by either: (i) leaving it unfiltered, (ii) filtering with wavelets, or (iii) filtering with the polynomial test functions \(\{\varphi_j(t):=\varphi(t -jh)\}\) as defined in 24 . To filter with wavelets, we use MATLAB’s wdenoise syntax with a sym12 wavelet. To filter with polynomials, we reconstruct the noisy signal via 23 using test functions generated from 24 . The polynomial parameters \((p,L,h)\) are fixed for each noise intensity. Appendix 7 contains diagnostics which empirically show that the selected parameters act as effective filters, justifying their use in this step. After filtering, we prune the first and last 250 steps to avoid boundary artifacts arising from the polynomial reconstruction. For each noise intensity and each filtering procedure, this process results in 100 segments consisting of 2000 steps of filtering training data, 5000 steps of filtered validation data, and 50,000 steps of clean testing data.

We train, validate, and test a distinct WKRR model over each of the 100 trajectories using a skip-connection scheme (see §2). We show results obtained from both the Gaussian kernel 12 and the DM kernel 13 . For each model, validation is performed according to the process described in §4.2 using the VPT error metric with threshold \(\gamma=0.3\), and hyperparameters \(N_c=20\), \(N_f=30\), and \(\nu = 200\). The reference parameters \((\epsilon^*,\lambda^*)\) are computed according to the heuristic described in Appendix 9. Typical validation landscapes are shown in Appendix 10.

For testing, we randomly select 500 integers in the interval \([1, \,47,000]\) to be used across all models. We use these indices as starting points to extract 500 segments of length 3000 from each of the 100 testing trajectories. For each model, we forecast over each of these 500 segments. This process results in a mean VPT for each of the 100 models which are computed using 50,000 total VPTs. We show a typical forecasting result obtained with the Gaussian kernel at 5% noise in Figure 4.

We now compare the following frameworks across noise intensities: (i) strong KRR trained and validated on unfiltered data, (ii) RAFDA trained on unfiltered data, (iii) strong KRR trained and validated on wavelet filtered data, (iv) strong KRR trained and validated on polynomial filtered data, (v) proposed WKRR trained on unfiltered data and validated on polynomial filtered data, (vi) proposed WKRR trained and validated on polynomial filtered data, (vii) proposed WKRR with a DM kernel trained on unfiltered data and validated on polynomial filtered data, and (viii) proposed WKRR with a DM kernel trained and validated on polynomial filtered data. Note that RAFDA does not receive filtered data due to its predictive mechanism. Note that knowledge of the noise level \(\sigma\) is required for RAFDA, but is not needed for either weak or strong KRR. The case of strong KRR applied to unfiltered training and validation data is included as a baseline. Each framework is tested over 100 models with identical data.

Figure 5 reports the empirical mean VPT density for each of the above frameworks in the top four panels, and reports the corresponding statistics and test function parameters in the table below. Across all noise levels, WKRR trained and validated on filtered data (Gaussian kernel - dark purple, DM kernel - dark red) consistently yields the best performance, while WKRR trained on unfiltered data (light purple and light red, respectively) is competitive or better than the strong approaches. Notice that several VPT density tails for strong KRR applied to filtered data (e.g., Figure 5, panel b) are quite long, indicating several models which produced VPTs well below the mean. In contrast, WKRR did not produce such outliers in any of our numerical experiments. This observation indicates that model validation may be more stable using a weak formulation. RAFDA exhibits strong performance for small noise, but its predictive performance worsens as noise increases.

We observe that the Gaussian kernel and DM kernel give nearly identical results for the L63 system, which has an ambient dimension of three and an intrinsic dimension of roughly 2.06 on the buttefly attractor [82]. This result suggests that the RBF kernel may be appropriate for low-dimensional systems, or for systems with ambient and intrinsic dimension roughly equal.

We remark that VPT measures short-term predictive capability, but does not capture long-term behavior. To study long-term forecasting fidelity, we train WKRR and RAFDA on identical data and forecast over a segment of length 40,000. We bin the forecasts and the ground truth in a histogram, which are shown in Figure 6 for each noise intensity. Both WKRR and RAFDA closely track the long-time marginal densities at small noise. However, while WKRR recovers most qualitative features of marginal densities at high noise, the performance of RAFDA deteriorates.

5.2 Kuramoto-Sivashinsky↩︎

We now apply WKRR to the Kuramoto-Sivashinsky (KS) system \[\label{eq:32KS32system} u_t + u u_x + u_{xx} + u_{xxxx} = 0, \quad x\in [0,J], \quad t\geq 0,\tag{40}\] with periodic boundary conditions \(u(t,0)=u(t,J)\). Following [52], we take \(J=22\) and discretize 40 in space over a uniform grid consisting of 64 points. Accordingly, we view the KS system as a 64-dimensional ordinary differential equation. In this setup, the KS system admits chaotic dynamics with Lyapunov exponent \(\Lambda \approx 0.043\) and has Kaplan-Yorke dimension approximately \(5.2\) [83], much smaller than the ambient 64-dimensional space.

To generate data we select \(u(0,x) = \sin(16\pi x/J)\) as an initial condition, and integrate in time using the Exponential Time Differencing with RK4 (ETDRK4) method [52], [84] with stepsize \(\Delta t = 0.01\) for \(N=3.05\cdot 10^7\) steps. The data is then downsampled by a factor of 10, and the first 50,000 points are pruned to avoid transient behavior. We then divide the remainder into 50 consecutive segments of 60,000 points. Each of these 50 segments consists of 4500 steps of training data, 5500 steps of validation data, and 50,000 steps of testing data. This data is corrupted with independent noise as in the L63 case.

The pre-processing step is performed identically to the L63 case, and results in 50 segments consisting of 4000 steps of filtered training data, 5000 steps of filtered validation data, and 50,000 steps of clean testing data.

To validate and test, we employ a direct-connection scheme (see §2). Otherwise, the validation process is similar to the L63 case. We validate according to the procedure described in §4.2 using the VPT error metric 36 with \(\gamma=0.5\), and hyperparameters \(N_c=20\), \(N_f=30\), and \(\nu=200\). The reference parameters \((\epsilon^*,\lambda^*)\) were computed according to the method in Appendix 9. We scale \(\lambda^*\) by a factor of 1000 for 5%, 10%, and 20% noise to promote more stable validation. Typical validation landscapes are shown in Appendix 10. Testing was performed as in the L63 case, resulting in 50 mean VPTs across 50 WKRR models. Figure 7 illustrates a typical forecast under 5% noise with Gaussian kernel.

In our numerical experiments, we consider the same frameworks as for the L63 case, with the exception of RAFDA which was omitted due to prohibitively large computational time.

Figure 8 depicts violin plots of the mean VPT densities across each of the frameworks and noise intensities. The table below provides the corresponding statistics and test function parameters. In all cases, WKRR with a Gaussian kernel applied to filtered data (dark purple) achieved comparable performance to the strong approaches. However, WKRR with a DM kernel applied to filtered data (dark red), achieved superior performance across noise levels. This result suggests that the DM kernel may be more appropriate for this problem. We suspect that the DM kernel exhibits superior performance because the flow map of the KS system is rather smooth, relating to a recent finding [71] that suggests the DM kernel is superior to the Gaussian kernel in identifying smooth labels, especially those spanned by lower-order eigenfunctions of the Laplace-Beltrami operator on a submanifold of \(\mathbb{R}^n\). For completeness, we provide results which show the superior performance of WKRR over KRR using the DM kernel in Appendix 11.

5.3 Validation Length Study↩︎

The validation procedure forms a substantial part of the WKRR framework. In this section, we perform a sensitivity analysis with respect to the validation segment length \(\nu\). Our goal is to demonstrate that the proposed method is robust across several validation strategies and does not require specific parameters to achieve good performance.

Our previous experiments averaged \(N_c=20\) trajectories consisting of \(\nu=200\) points over a coarse parameter grid, and subsequently averaged \(N_f=30\) trajectories consisting of \(\nu=200\) points over a fine parameter grid. This modeling choice provided consistent forecasting performance across noise levels for both the L63 and KS systems.

Here, we compare this baseline validation strategy to four others: (i) \(N_c=80\), \(N_f=120\), and \(\nu=50\); (ii) \(N_c = 40\), \(N_f=60\), and \(\nu=100\); (iii) \(N_c=10\), \(N_f=15\), and \(\nu=400\); and (iv) \(N_c=8\), \(N_F=12\), and \(\nu=500\). Otherwise, the experimental procedure is identical to the baseline cases. These frameworks were considered to provide validation segments whose lengths are both larger and smaller than the expected model forecast. To enact a direct comparison, we define \(\nu_{\text{Lyap}} = \nu \cdot \Delta t \cdot \Lambda\), which has VPT units and may be directly compared to mean VPT. Here, we only show results with the Gaussian kernel since the same conclusion is also valid for the DM kernel.

Figure 9 (a) reports the mean and standard deviation of the mean VPTs at 1% noise (left) and 20% noise (right) over 100 models for the L63 system, while panel (b) reports these statistics over 50 models for the KS system. We remark that the mean VPT remains stable across the validation strategies, highlighting the robustness of the proposed validation procedure. We emphasize that comparable VPTs were achieved from validation segments which are both larger and smaller than the expected testing VPT. While validating on segments whose length exceeds the expected forecast horizon is a typical approach to support appropriate generalization, our results suggest that WKRR still achieves comparable generalization even if it uses shorter validation segments.

5.4 Clean Data Performance↩︎

In this section, we use WKRR and classical strong KRR to study the L63 system 39 and the KS system 40 as above, but with clean \((\sigma = 0)\) training and validation data. The training and validation data are not pre-processed, but are otherwise constructed identically to the above cases.

The validation scheme used here differs from that described in §4.2 only in that we now utilize \(N_c=N_f=3\) (possibly overlapping) validation segments consisting of \(\nu=1500\) points, following the setup in [52]. Validating over longer segments improves performance because the expected forecast horizon under clean data is significantly longer than the noisy cases. For WKRR, we scale the reference regularization parameter \(\lambda^*\) by a factor of \(10^{-6}\) to promote stable validation.

Figure 10 reports the empirical VPT density over 100 models for the L63 system, and 50 models for the KS system, respectively. The tables below provide corresponding quantitative statistics. We remark that WKRR loses only a small amount of accuracy compared to the strong KRR approach, indicating its robustness across both clean and noisy data. We remark that WKRR with a DM kernel exhibits superior performance for both systems. Previous work showed that classic strong KRR with a DM kernel often exhibits superior performance compared to a Gaussian kernel over clean data [52]. The present results suggest the robustness of the DM kernel with a weak formulation.

One might argue that selecting appropriate test functions for WKRR make the method difficult to use in practice. However, the selection is straightforward when clean data is available. One procedure is as follows. For a given family of test functions, reconstruct a portion of the given training data via 23 . Then, compare the available clean data to its reconstruction using a suitable error metric, and select the best-performing test function.

5.5 Experimental Data: Fluid Dynamics↩︎

We now apply WKRR with both the Gaussian and DM kernel to experimental fluid data, made available by [14] as a Community Challenge. The data consists of streamwise velocity \(\mathbf{u}\) and wall-normal velocity \(\mathbf{v}\) measured on cavity flow. The dataset includes measurement noise as well as fluctuations due to the inherent turbulence in the flow. Furthermore, the flow is strongly convective and heavily depends on the time-varying inlet boundary condition on the left, that is time-varying and not precisely known due to the experimental setup. Due to the unknown inlet boundary condition during prediction, it is challenging to achieve long forecasting horizon in this setting. We refer to Problem 2.1 in [14] for details of the experimental procedure and data collection process.

We set \[\mathbf{U}= \begin{bmatrix} \mathbf{u}& \mathbf{v} \end{bmatrix} \in \mathbb{R}^{N \times n},\] where \(n=2\cdot 7854=15708\) is the spatial dimension and \(N\) denotes the number of available data timeslices. We are given \(N_{\text{train}}=12,800\) samples as training data, and 32 segments consisting of \(N_{\text{seg}}=70\) of validation data. The original challenge as described in [14] is to forecasting 30 timesteps past the validation data for each of the segments. However, because this data is not available to us, we instead break apart the 32 segments into \(N_{\text{val}}=40\) for model validation, and \(N_{\text{test}}=30\) for model testing. This procedure allows us to quantitatively analyze the performance of WKRR and compare it to other baseline methods. To fix notation, let \(\mathbf{U}_{\text{train}}=[\mathbf{u}_{\text{train}} \; \mathbf{v}_{\text{train}}]\in \mathbb{R}^{12800 \times 15708}\) denote the given training data, let \(\mathbf{U}^{(q)}_{\text{val}}=[\mathbf{u}_{\text{val}} \; \mathbf{v}_{\text{val}}]\in \mathbb{R}^{40 \times 15708}\) denote the validation data, and let \(\mathbf{U}^{(q)}_{\text{test}}=[\mathbf{u}_{\text{test}} \; \mathbf{v}_{\text{test}}]\in \mathbb{R}^{30 \times 15708}\) denote the testing data for \(q=1,\dots,32\).

We do not directly apply WKRR to the given data, because its size makes computations prohibitively expensive. Instead, we adopt the following pre-processing procedure. Set \(\mu_u, \mu_v \in \mathbb{R}^{1\times 7854}\) to be the time-average of \(\mathbf{u}_{\text{train}}\) and \(\mathbf{v}_{\text{train}}\), respectively. Similarly set \(\sigma_u,\sigma_v\in\mathbb{R}^{1\times 7854}\) to be the standard deviation of the training data. Then for both training and validation data \(\mathbf{u},\mathbf{v}\), we normalize \[\hat{\mathbf{U}}= \begin{bmatrix} \frac{\mathbf{u}-\mu_u}{\sigma_u} & \frac{\mathbf{v}-\mu_v}{\sigma_v} \end{bmatrix},\] where division is interpreted component-wise. We then implement a POD procedure by taking a truncated SVD \[\hat{\mathbf{U}}= \mathbf{U}_r \mathbf{S}_r \mathbf{V}_r^\top,\] where \(r\) modes are retained. The normalized data is then projected onto the POD modes \[\tilde{\mathbf{U}}= \hat{\mathbf{U}}\mathbf{V}_r = \begin{bmatrix} \tilde{\mathbf{u}}& \tilde{\mathbf{v}} \end{bmatrix}.\] Finally, we normalize per mode by writing \[\mathbf{U}^\dagger = \begin{bmatrix} \frac{\tilde{\mathbf{u}}-\mathfrak{u}_u}{\varsigma_u} & \frac{\tilde{\mathbf{v}}-\mathfrak{u}_v}{\varsigma_v} \end{bmatrix},\] where \(\mathfrak{u}_u, \mathfrak{u}_v \in \mathbb{R}^{1\times r}\) and \(\varsigma_u, \varsigma_v \in \mathbb{R}^{1\times r}\) are the mean and standard deviation of the POD modes that arise from the training data. We emphasize that the mean and standard deviations are computed only from the training data, although these values are used to normalize both the training and validation data.

Polynomial parameters are chosen according to the following procedure. We select a segment of training data consisting of 1000 timesteps. The data is normalized according to the above procedure. We then filter and reconstruct the normalized signals over a grid of parameter triples \((p,L,h)\). The reconstruction is the unnormalized per-mode, projected onto original coordinates by the action of \(V_r^\top\), and multiplied by the standard deviation. The result is a mean-subtracted signal. We then compare this signal to the ground truth (also mean-subtracted) using the NSME metric 37 , and select the parameter pair which minimizes the error.

We then apply WKRR to \(\mathbf{U}^\dagger_{\text{train}}\). Both the training and validation data is normalized according to the above procedure. Validation is performed according to the procedure described in §4.2, but with the following modifications. The coarse grid is initialized to \(H_c = [10^2,10^4]\times [10^{-16},10^{-13}]\), and NMSE (computed over the mode-normalized coefficients) is averaged over \(N_c=3\) and \(N_f=5\) segments of length \(\nu = 15\) for the coarse and fine grids, respectively. To alleviate computational time, we validate on only one of the 32 validation segments.¹ The resulting model is tested on all 32 testing segments. This testing data is unnormalized to mean-subtracted physical coordinates and compared to the mean-subtracted ground truth using the NMSE metric. This metric is consistent with that considered in the Community Challenge paper [14].

We compare the performance of WKRR with LSTM, used as a baseline in the Community Challenge paper [14]. The LSTM implementation, provided by [14], was modified to use the same normalized data as WKRR. The baseline results in [14] using LSTM were computed with \(r=25\), noting that increasing \(r\) may destabilize the training of the LSTM networks. In the following, we will compare WKRR and LSTM with both \(r=25\) and \(r=100\).

Figure 11 shows a typical forecast for each method with \(r=25\) over 30 timesteps. We observe that WKRR with the Gaussian and DM kernel gives visually indistinguishable results, and so we report only results using the Gaussian kernel for this case. The top block of four rows correspond to \(\mathbf{u}\) data, while the last block of four rows corresponds to \(\mathbf{v}\) data. In each of the two blocks, we compare the ground truth (first row) with the \(r=25\) POD representation (second row), the WKRR forecast (third row), and the LSTM forecast (fourth row) over one of the 32 testing segments. The error plots in Figure 12 report the NMSE for each of the 32 trials (light gray), the average error over all 32 trials (solid colored line), and the corresponding standard deviation (light shaded band). The dashed colored lines correspond to the trajectory represented in the two blocks.

Both methods produce predictions that reasonably reproduce the POD representation of the data. Our results indicate that both methods are competitive in this regime.

We now repeat the forecasting procedure taking \(r=100\) POD modes. In this case, we report both the Gaussian and DM kernel results. Figures 13 and 14 report the results in the same style as above. The Gaussian kernel results are shown in purple, while the DM kernel results are shown in red. We observe that WKRR with both kernels exhibits noticeably more consistent qualitative agreement and smaller mean error over short-term horizons. This improved performance for higher-dimensional data highlights a key advantage of WKRR.

The performance of WKRR with the Gaussian and DM kernel differs in this regime. Over longer horizons, WKRR with the DM kernel has noticeably smaller quantitative error. However, this metric may be misleading. Qualitatively, forecasts with the DM kernel tend to produce overly smooth predictions which may fail to capture the long-term statistical behavior of the dynamics. In contrast, while predictions with the Gaussian kernel produce larger quantitative error, the forecasts appear to be qualitatively more consistent with the underlying dynamics over long time periods. Over short time periods, the two kernels perform similarly, with a slight edge to the DM kernel.

6 Discussion↩︎

In this paper, we propose Weak-form Kernel Ridge Regression (WKRR) as a data-driven, noise-robust learning framework. The proposed approach is computationally cheaper than classical strong-form KRR, and demonstrates competitive performance in the presence of noise over a range of chaotic, high-dimensional, and experimental data-sets. While selecting appropriate model hyperparameters via validation is often a difficult task that may be sensitive to modeling assumptions, we perform sensitivity studies to demonstrate that WKRR can achieve robust performance with multiple strategies across several baseline systems and noise levels. Furthermore, we show that with appropriately chosen test functions, WKRR applied to clean data greatly reduces computational complexity with only marginal loss in predictive performance. This observation positions WKRR as a flexible framework across a range of datasets when the underlying noise level is unknown. Finally, the success of kernel-based learning frameworks depends strongly on the choice of kernel function. We demonstrate that WKRR achieves competitive performance with the standard Gaussian kernel over a range of baseline and experimental data, highlighting the method’s simple implementation. We also consider the Diffusion Maps (DM) kernel [52], [70] as an alternative choice, and show that the DM kernel can lead to increased predictive performance for systems whose invariant set dimension is much lower than the ambient dimension, even with noisy observational data. The success of WKRR with multiple kernels broadens the applicability of the method and demonstrates that its noise robustness is not tied to a specific choice of kernel function. Moreover, these findings suggest that the practitioner has a modeling choice: a Gaussian kernel reduces runtime and appears to be effective for low-dimensional systems, while a DM kernel may lead to increased forecast horizon for systems with complicated geometry. We emphasize that the proposed WKRR approach has advantages over existing learning methods which have incorporated a weak formulation [36], [60]–[63] because it (i) does not require a choice of dictionary functions that span the target vector field, and (ii) does not require precise parameter tuning common to many machine learning architectures.

Despite the success of WKRR, several open questions remain. We observed that model validation using WKRR often produced fewer outliers with poor forecast horizons than strong KRR over pre-processed data. We hypothesize that validating over weak-form coefficients may lead to improved validation landscapes than using data in the original coordinates. It would be fruitful to pursue this observation, which may lead to more efficient validation strategies.

Selecting appropriate test functions is a challenge common to all weak-form learning approaches. The bias-variance decomposition in §3.3 suggests that using a smaller family of test functions reduces error due to noise. In contrast, one generally expects that using a larger family of test functions reduces error due to signal reconstruction, although this analysis is complicated and typically depends on the specific choice of test function. Developing a systematic approach to balance reconstruction and noise error constitutes a challenging yet worthwhile goal.

We demonstrated numerically that training over noisy observational data is significantly more effective with WKRR than with classical KRR. While we observed that pre-processing training data often leads to marginally improved forecast horizons, the role of pre-processing in training and validation is not yet well-understood. Incorporating multiple bandwidth parameters, as opposed to a single scalar considered in this work, may also influence training and validation. Establishing theoretical convergence results for WKRR in these contexts would strengthen the proposed framework.

The present manuscript assumes that the given data captures all states of an underlying dynamical system. However, in practice, given data may be only partially observed. Additionally, many physical systems may depend on parameters or external non-autonomous forcings. Extending WKRR to handle partially observed data, time-dependent data, or data which depends on parameters would broaden the applicability of WKRR.

Acknowledgment↩︎

This work is partially supported under the NSF grants DMS-2505605 and CMMI-2340266, and the ICDS Penn State seed grant.

Data Availability↩︎

The code to produce the figures in this paper is publicly available at https://github.com/MaxKreider/WKRR.

7 Test Function Parameters↩︎

In this section, we provide additional diagnostics for polynomial test functions as filters. We also compare with wavelet pre-processing as a standard baseline filter.

In our numerical experiments in §5, we apply WKRR to both the L63 system 39 and the KS system 40 . In both cases, the training and validation data were corrupted with various noise intensities. In Figures 5 and 8, we demonstrated that both strong and weak forms applied to polynomial filtered data resulted in higher mean VPT than the unfiltered case for a specific set of polynomial parameters \((p,L,h)\). Here, we numerically demonstrate that the weak formulation with these parameters effectively filters the data, explaining the mechanism underlying the success of WKRR in these scenarios.

For each of the 100 L63 models or 50 KS models, we explicitly reconstruct the noise-corrupted training data with polynomial test functions via 23 . We also filter this training data with wavelets using MATLAB’s wdenoise syntax with a sym12 wavelet. In both cases, we prune the first and last 250 data points to avoid boundary artifacts. In this section, we denote the pruned ground truth data as \(\mathbf{u}_i = [u_i^{(1)},\dots,u_i^{(n)}]^\top\) and the pruned filtered data as \(\hat{\mathbf{u}}_i = [\hat{u}_i^{(1)},\dots,\hat{u}_i^{(n)}]^\top\) for \(i=1,\dots,\mathfrak{N}\). We measure the filter performance using three error metrics: \[\label{eq:32three32error32metrics} \begin{align} \mathcal{E}_{RMSE} &= \frac{1}{n}\sum_{\ell = 1}^n \sqrt{\frac{1}{\mathfrak{N}}\sum_{i=1}^{\mathfrak{N}} \left|u_i^{(\ell)} - \hat{u}_i^{(\ell)}\right|^2}, \\ \mathcal{E}_\theta &= \frac{1}{\mathfrak{N}}\sum_{i = 1}^\mathfrak{N} \arccos\left(\frac{\mathbf{u}_i \cdot \hat{\mathbf{u}}_i}{\|\mathbf{u}_i\|_2\|\hat{\mathbf{u}}_i\|_2}\right), \\ \mathcal{E}_{LSD} &= \frac{1}{n}\sum_{\ell=1}^n \sqrt{\frac{1}{M} \sum_{m=1}^M \left|10\log_{10}\mathcal{P}[\mathbf{u}^{(\ell)}]_m - 10\log_{10}\mathcal{P}[\hat{\mathbf{u}}^{(\ell)}]_m\right|^2}, \end{align}\tag{41}\] where \(\mathcal{P}[\cdot]_m\) is the \(m\)th component of the power spectral density of the signal, assumed to be of length \(M\).

The metric \(\mathcal{E}_{RMSE}\) computes RMSE over time, and averages these errors over spatial dimensions. It provides a standard measure of pointwise error. The metric \(\mathcal{E}_{\theta}\) measures the angle between the truth and reconstruction at a fixed sample time, and averages the result over all available sample times. Loosely speaking, it provides a measure of correctness of direction. The metric \(\mathcal{E}_{LSD}\) computes the log-spectra distance between the components of the true and reconstructed signals, and averages these errors over spatial dimension. It provides a standard measure of spectral difference between two signals.

We report the mean and standard deviation of these error metrics averaged over 100 trials for the L63 system in Figure 15 and averaged over 50 trials for the KS system in Figure 16. In all cases, the data \(\hat{\mathbf{u}}\) is filtered with wavelets (blue), polynomial test functions (orange), or left unfiltered (gray). Across noise levels, both wavelets and polynomials give comparable results in terms of RMSE and angle differences. The polynomial filter produces filtered data with smaller LSD errors. The resulting error metrics are significantly improved relative to the unfiltered baseline, indicating that both wavelets and polynomials effectively filter the data. In particular, the success of the polynomials in filtering the data justifies the \((p,L,h)\) parameter selection used in our numerical experiments.

8 H Matrix↩︎

Suppose that data \(\mathbf{U}= [\mathbf{u}_1,\dots,\mathbf{u}_N]^\top\in \mathbb{R}^{N\times n}\), with corresponding clean data \(\mathbf{U}_{\text{clean}}\), is given. We assume that the underlying clean data is periodic, which implies that \(k^* = N/(h/\Delta t)\). We also assume that \(\mathbf{W}= \mathbf{I}\).

Our main goal in this section is to show that the general expression for the bias given in 25 , under the above ideal assumptions, can be expressed in the form \[\label{eq:32bias32goal} B_\ell^2 = \|(\mathbf{P}- \mathbf{I})\mathbf{U}_{\text{clean}}^{(\ell)}\|_F^2 = \left\|(\mathbf{H}- \mathbf{I})\mathcal{F}\left[\mathbf{U}_{\text{clean}}^{(\ell)}\right]\right\|_F^2,\tag{42}\] in terms of the Fourier transform of the \(\ell\)th component of the true underlying signal, and the matrix \(\mathbf{H}\in \mathbb{R}^{N\times N}\) as defined in 28 . This approach is inspired by the approach taken in [64], which performs computations in continuous time.

We first observe that when \(\mathbf{W}=\mathbf{I}\), we have that \(\mathbf{P}= \boldsymbol{\Psi}^\top ( \boldsymbol{\Psi}\boldsymbol{\Psi}^\top)^{-1} \boldsymbol{\Psi}\). Notice that the entries of the term in parenthesis are quadrature approximations of the integral \[\begin{align} [( \boldsymbol{\Psi}\boldsymbol{\Psi}^\top)]_{ij} &= \int_{\mathbb{R}} \text{d}t\; \varphi_i(t) \varphi_j(t) = \int_{\mathbb{R}} \text{d}\tau\; \varphi(\tau) \varphi(\tau - (j-i)h), \end{align}\] which shows this matrix is Toeplitz, and leveraging periodicity, circulant.

This observation is important because the discrete Fourier Transform (DFT) diagonalizes circulant matrices. Let \(\mathbf{F}\in \mathbb{R}^{N\times N}\) with entries \(\mathbf{F}_{mn} = 1/\sqrt{N}\exp(-2\pi m ni/N)\). Notice that \(\mathbf{F}^*\mathbf{F}= \mathbf{I}_{N\times N}\). Consider the frequency domain representation of the test functions \(\hat{\boldsymbol{\Psi}} = \mathbf{F}\boldsymbol{\Psi}^\top\in \mathbb{C}^{N\times k^*}\). Leveraging periodicity and straightforward manipulations, one can show that the \((m,j)\) entry of this object has the form \[\begin{align} \hat{\boldsymbol{\Psi}}_{mj} &= \sum_{p=0}^{N-1} \frac{1}{\sqrt{N}} \exp(-2\pi m p i /N) \varphi[p-jh] = \exp(-2\pi m j i /k^*) \hat{\varphi}_m, \end{align}\] where \(\hat{\varphi}_m\) is the \(m\)th component of the DFT applied to the mother polynomial function. Let \(\mathbf{E}\in \mathbb{C}^{N\times k^*}\) be the matrix whose \((m,j)\) entries are \(\exp(-2\pi m j i /k^*)\), and let \(\mathbf{D}\in \mathbb{C}^{N\times N}\) be the diagonal matrix whose \(m\)th entries is \(\hat{\varphi}_m\). Then, we have \[\hat{\boldsymbol{\Psi}} = \mathbf{F}\boldsymbol{\Psi}^\top = \mathbf{D}\mathbf{E}.\] Moreover, note that the columns of \(\mathbf{E}\) are \(k^*\)-periodic, so the matrix \(\mathbf{E}\) is composed of \(h\) blocks of size \(k^*\times k^*\) stacked on top of each other. Each of these blocks is the \(k^*\times k^*\) DFT matrix \(\sqrt{k^*}\mathbf{F}_{k^*\times k^*}\). Letting \(\mathbf{J}\) be the block identity matrix of size \(N\times k^*\), we can write \[\hat{\boldsymbol{\Psi}} = \mathbf{F}\boldsymbol{\Psi}^\top = \sqrt{k^*}\mathbf{D}\mathbf{J}\mathbf{F}_{k^*\times k^*},\]

Now that this setup is complete, we can perform several substitutions to recover the \(\mathbf{H}\) matrix above. Consider \[\begin{align} (\mathbf{P}-\mathbf{I}) &= (\boldsymbol{\Psi}^\top (\boldsymbol{\Psi} \boldsymbol{\Psi}^\top)^{-1} \boldsymbol{\Psi}-\mathbf{I}) \\ &= (\boldsymbol{\Psi}^\top (\boldsymbol{\Psi} \mathbf{F}^*\mathbf{F}\boldsymbol{\Psi}^\top)^{-1} \boldsymbol{\Psi}-\mathbf{F}^*\mathbf{F}) \\ &= (\boldsymbol{\Psi}^\top (k^*\mathbf{F}_{k^*\times k^*}^* [\mathbf{J}^* |\mathbf{D}|^2 \mathbf{J}] \mathbf{F}_{k^*\times k^*})^{-1} \boldsymbol{\Psi}-\mathbf{F}^*\mathbf{F}). \end{align}\] Let \(\mathbf{S}= \mathbf{J}^* |\mathbf{D}|^2 \mathbf{J}\in \mathbb{R}^{k^*\times k^*}\). The \(m\)th entry of the diagonal matrix \(|\mathbf{D}|^2\) is \(|\hat{\phi}_m|^2\). The action of \(\mathbf{J}\) on either side is to sum up the “blocks”, so \(\mathbf{S}\) is a diagonal matrix whose \(m\)th entry is \(\sum_{\ell=0}^{h-1} |\hat{\phi}_{m+\ell k^*}|^2\). Continuing, we have \[\begin{align} (\mathbf{P}-\mathbf{I}) &= \frac{1}{k^*}(\boldsymbol{\Psi}^\top\mathbf{F}^*_{k^*\times k^*} \mathbf{S}^{-1}\mathbf{F}_{k^*\times k^*}\boldsymbol{\Psi}-\mathbf{F}^*\mathbf{F}) \\ &= \frac{1}{k^*}(\mathbf{F}^* \hat{\boldsymbol{\Psi}} \mathbf{F}^*_{k^*\times k^*}\mathbf{S}^{-1}\mathbf{F}_{k^*\times k^*}\hat{\boldsymbol{\Psi}}^*\mathbf{F}-\mathbf{F}^*\mathbf{F}) \\ &= \frac{1}{k^*}\mathbf{F}^*\left([\sqrt{k^*}\mathbf{D}\mathbf{J}\mathbf{F}_{k^*\times k^*}\mathbf{F}^*_{k^*\times k^*}] \mathbf{S}^{-1}[\sqrt{k^*}\mathbf{D}\mathbf{J}\mathbf{F}_{k^*\times k^*}\mathbf{F}_{k^*\times k^*}^*]^*-\mathbf{I}\right)\mathbf{F} \\ &=\mathbf{F}^*\left((\mathbf{D}\mathbf{J}) \mathbf{S}^{-1}(\mathbf{J}^*\mathbf{D}^*)-\mathbf{I}\right)\mathbf{F}. \end{align}\] Define \(\mathbf{H}= (\mathbf{D}\mathbf{J}) \mathbf{S}^{-1}(\mathbf{J}^*\mathbf{D}^*)\in \mathbb{R}^{N\times N}\). One can show, using the block structure of \(\mathbf{J}\), that the entries of \(\mathbf{H}\) are of the form \[\mathbf{H}_{mn} = \begin{cases} \frac{\hat{\varphi}_m \overline{\hat{\varphi}_n}}{\sum_{\ell=0}^{h-1} |\hat{\varphi}_{r+\ell k^*}|^2}, & m = n \;(\text{mod } k^*), \\ 0, & \text{otherwise}. \end{cases}\] where \(r = m \text{ mod } k^*\). This computation recovers the form 42 .

9 Validation Heuristic↩︎

In this section, we review a validation heuristic to compute reference bandwidth and regularization parameters \((\epsilon^*,\lambda^*)\) which was utilized in [52] for strong KRR applied to clean data. Here, we modify it slightly to improve robustness in the case of noisy data. We emphasize that the procedure described below is a heuristic which has been observed to provide reasonable results over a range of examples [52], [85]. While it often provides a good starting point, it should not be interpreted as an optimal approach that is appropriate for all problems.

Suppose that training data \(\mathbf{u}_i\) for \(i=1,\dots,N\) of the form 1 is given. This data could be filtered, but is not required to be. Let \(L_*\) be the maximum pairwise \(L_2\) distance between these data points. Let \[\rho(\mathbf{x},\mathbf{y};\eta) = \exp(-\|\mathbf{x}-\mathbf{y}\|_2^2/(L_*^2 \eta)\] be a standard Gaussian RBF kernel with bandwidth \(\eta\). Further define \[S(\eta) = \frac{1}{N^2}\sum_{i,j=1}^N \rho(\mathbf{u}_i,\mathbf{u}_j;\eta).\] Following [52], [86], this fact leads one to consider \[V(\eta) = 2\frac{\eta}{S(\eta)}\frac{\text{d}S(\eta)}{\text{d}\eta},\] and to consider the problem \[\eta^* = \arg\max_\eta V(\eta).\] For clean data, \(V(\eta)\) is often uni-modal, and an appropriate choice of \(\eta^*\) is unambiguous. For noisy data, we numerically observe that \(V(\eta)\) may admit multiple peaks with similar values. In practice, if \(\eta^*_r\) denotes the arguments at which the local maxima occur, we recommend selecting \(\eta^* = \max_r \eta^*_r\). The goal of this procedure is to avoid selecting an extremely small bandwidth, which could lead to poor generalization.

Once \(\eta^*\) has been computed, we use it to define a reference bandwidth \(\epsilon^*\) via the following formula \[\epsilon^* = 250 L_*^2\eta^*.\]

To select the reference regularization parameter \(\lambda^*\), let \(k(\mathbf{x},\mathbf{y};\epsilon) = \exp(-\|\mathbf{x}-\mathbf{y}\|_2^2/\epsilon)\) be the Gaussian RBF kernel, and let \(\mathbf{K}_{\text{RBF}}(\epsilon)\) be the corresponding Gram matrix whose \((i,j)\) entry is \(k(\mathbf{u}_i,\mathbf{u}_j;\epsilon)\). We choose \(\lambda^*\) to be the minimum eigenvalue of \(\mathbf{K}_{\text{RBF}}(\epsilon)\).

10 Validation Landscapes↩︎

In this section, we depict typical WKRR validation landscapes across noise levels for the L63 system 39 in Figure 17, for the KS system 40 in Figure 18, and for the experimental fluid data in Figure 19.

Note that in the case of the experimental fluid data, over-regularizing may improve quantitative performance in terms of the error metric at the cost of qualitative fidelity, i.e., over-regularizing often causes WKRR to rapidly converge to a mean flow that does not capture important qualitative features of the data.

11 KS Diffusion Maps Kernel Data↩︎

Here, we repeat the experimental procedure for the KS system described in §5.2 with the DM kernel 13 instead of the Gaussian kernel 12 . The results are shown in Table 1, and demonstrate the superior performance of WKRR over the strong formulation across all noise levels.

5pt

References↩︎

---

1. Several validation segments were randomly chosen and forecasting results did not change significantly, indicating that the particular choice of validation segment does not greatly alter the outcome.↩︎