PIKS: Universal Physics-Informed Kernel Methods

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

PIKS: Universal Physics-Informed Kernel Methods

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

PIKS: Universal Physics-Informed Kernel Methods Joachim Bona-Pellissier1 , Giacomo Meanti1 , Matteo Santacesaria2 , Lorenzo Rosasco1,3 1

MaLGa Center, DIBRIS, Università degli Studi di Genova, Genoa, Italy MaLGa Center, DIMA, Università degli Studi di Genova, Genoa, Italy 3 Istituto Italiano di Tecnologia, Genoa, Italy

arXiv:2607.27062v1 [stat.ML] 29 Jul 2026

[email protected], [email protected] [email protected], [email protected]

Abstract Physics-informed machine learning incorporates physical principles—often expressed via differential operators—into data-driven models. While physics-informed neural networks (PINNs) dominate empirical applications, the complexity of neural network architectures and optimization landscapes hinders the development of a corresponding learning theory. In turn, kernel methods offer an appealing alternative with closed-form solutions and analytical tractability, yet existing guarantees primarily cover the well-specified setting where the target belongs to the native Reproducing Kernel Hilbert Space (RKHS). This imposes unrealistic regularity assumptions that physical targets often fail to satisfy. In this paper, we introduce and analyze Physics-Informed Kernel methodS (PIKS). We establish the universal consistency of PIKS for linear differential constraints, proving that for universal kernels (such as Gaussian or Matérn), the estimator asymptotically learns the target while satisfying physical constraints. We further derive finite-sample bounds under suitable source conditions. Our analysis is based on extending classical operator-theoretic analysis of kernel methods to physics-informed machine learning. Numerical experiments demonstrate that PIKS can be competitive with PINNs and traditional finite element methods.

Keywords: physics-informed machine learning, kernel methods, universal learning, statistical learning theory, scientific machine learning

Introduction

Physics-informed machine learning (PIML) combines data-driven statistical learning with additional knowledge about a problem’s physical properties. Indeed, in many applications governed by physical principles, the target function is constrained by relations involving its derivatives or, more generally, the action of differential operators (Cuomo et al., 2022). These relations provide further information that can be incorporated into the learning problem as functional constraints, such as prescribed gradients, partial differential equation (PDE) residuals, or conservation laws (Raissi et al., 2019; Karniadakis et al., 2021; Rackauckas et al., 2020). There can be several benefits in leveraging structural knowledge for a learning task. When the amount of labeled data is limited, for instance because measurements are challenging to collect, physical knowledge about the target function can help circumvent data scarcity. Furthermore, simply approximating the target (as in classical regression) may not be enough in certain settings 1

where the physical consistency of a solution can be as important as its accuracy. For example even a small violation of the conservation of energy in the force-field function may crash a molecular dynamics simulation (Fu et al., 2023). Augmenting the learning problem with structural constraints that the unknown target function must satisfy has been shown to solve this challenge effectively, dramatically improving sample efficiency, stability, and out-of-distribution robustness (Cuomo et al., 2022; Quarteroni et al., 2025). Starting with the introduction of physics-informed neural networks (PINNs) (Raissi et al., 2019), the use of machine learning — in particular of neural networks — for tackling forward and inverse physics-informed problems has surged in popularity. Engineering and applied aspects of the problem have taken center stage: improved neural network architectures and new ways of framing the problems have resulted in numerous success stories across domains (see e.g. Toscano et al., 2025; Zhao et al., 2024, among many examples). As is common with deep learning, however, theoretical understanding has lagged behind empirical progress. In this paper we propose to tackle a basic yet fundamental question concerning the universal consistency of learning methods for linear physicsinformed problems. Given a model which is fixed a-priori, can we guarantee that for any target function, the model will learn the target as the data grows to infinity? This question has been studied in standard supervised learning for many families of models including neural networks and kernel methods, and relies on the universality of the corresponding function class, i.e. when a hypothesis class is rich enough to approximate any function, typically in a L2 or L∞ sense. However, as we will show, standard notions of universality fall short in the physics-informed setting, as the additional structure requires a more careful analysis. In this paper, we study the universal consistency of kernel methods for physics-informed machine learning. Kernel methods provide a tractable framework for obtaining learning guarantees (see e.g. Smale and Zhou, 2007; Caponnetto and De Vito, 2007; Steinwart and Christmann, 2008; Blanchard and Mücke, 2018). They are also naturally suited to incorporating linear functional constraints in the learning objective, a property which has been used at least since the 1970s (Kimeldorf and Wahba, 1971) and has been applied to Hermite-Birkhoff interpolation (Kimeldorf and Wahba, 1971), meshless methods for PDEs (Fasshauer, 1996; Wendland, 2004), or self-supervised learning with manifold constraints (Belkin et al., 2006) among others. We review these works extensively in Section 2.2, but to the best of our knowledge, certain fundamental properties such as universal consistency have not yet been proven in the physics-informed case. Indeed, the error analyses in the literature usually hold in the well-specified setting, in which the target function u∗ belongs to the same reproducing kernel Hilbert space (RKHS) H from which the estimator is chosen. Instead, universal consistency is a property that is relevant in the misspecified setting, i.e. the case u∗ ∈ / H. This setting is important in practice: for example a Sobolev space H s (Ω) is an RKHS only when s > d/2 but the target function for real-world problems in high dimensions may not be smooth enough (i.e., it may belong to a Sobolev space with exponent s ≤ d/2). Analyses in the misspecified setting allow us to obtain learning guarantees even when the target is not smooth enough to belong to any RKHS, or when perfect knowledge about the physical process is not available, leading to an inaccurate choice of kernel (and thus of space H). We focus on a specific kernel-based algorithm for physics-informed machine learning which we refer to as PIKS (for Physics Informed Kernel methodS ). In Section 2 we introduce the physicsinformed learning problem, which amounts to regression with an additional structural constraint. Like in PINNs, the constraint is enforced by adding a physics-informed term to the training loss. Throughout the paper we emphasize that our general formalism can apply to diverse physicsinformed settings. In particular, we show in Section 2.1 that it covers the setting of pure PDE solving, where measurements of the target values are only available on the boundary of the domain (while the PDE provides information inside the domain). But in general, measurements of the 2

target can be available anywhere in the domain, as PIKS covers any problem of learning with linear constraints. In Section 3, we precisely formalize the PIKS algorithm and in Section 4 proceed with its theoretical analysis. In the latter, after listing the working assumptions, we state our main result, Theorem 1, which shows consistency of PIKS. While Theorem 1 is formulated for general linear operators, in the following subsection, we apply this result in the context of differential operators on Sobolev spaces, which is very common in applications. We then focus more specifically on the setting of PDE solving (already described in Section 2.1), where measurement data is only available on the boundary of the domain. Finally, we show that under additional assumptions, we can go further than the consistency, and obtain finite-sample bounds, as stated in Theorem 2. In Section 5 we study the empirical performance of PIKS in different settings, and show that it compares well to other methods in the literature. In summary, the main contributions of this work are: 1. We establish a set of general assumptions under which the PIKS estimator is a universal learner : it is able to both learn the target and approximate the physical constraint asymptotically, even in the misspecified setting in which the target u∗ does not belong to the native RKHS of the estimator. We analyze in more detail the typical case of differential operators on Sobolev spaces. 2. We make this general result concrete in two specific settings: one in which target measurements are available on the domain boundary (typical of PDE problems) and one in which they are available inside the domain. In the case of elliptic differential operators we show how to use the PDE’s regularity with our main theorem to obtain stronger convergence results. 3. We provide convergence rates under stronger assumptions on the target, which take the form of a source condition defined as powers of an integral operator. We illustrate the rates with an example of the Laplacian on periodic functions. 4. We demonstrate empirically that PIKS is competitive with recent kernel-based methods as well as physics-informed neural approaches on representative PDE solving tasks. We further study in practice the impact of model misspecification on the performance of the algorithm. When the target function is rougher than the base kernel, the convergence is slower, but remains competitive with classical FEM solvers.

Physics-informed statistical learning

We consider learning problems with an additional constraint in the form of a linear operator which encodes the problem’s physics. More precisely, let u∗ ∈ F where F is a space of maps from X ⊆ Rd to R. Let G be another space of maps from X to R and let D : F → G be a linear operator. Consider input data random variables X ∼ ρX and Z ∼ ρZ in X as well as noise random variables ϵ and η (centered and independent from X and Z respectively), and let Y = u∗ (X) + ϵ,

W = Du∗ (Z) + η.

(1)

Two datasets (xi , yi )ni=1 and (zj , wj )m j=1 can be obtained by sampling i.i.d. copies of (X, Y ) and of (Z, W ) respectively. The objective is to minimize the expected risk min R(u), R(u) = E (u(X) − Y )2 + E (Du(Z) − W )2 (2) u∈F

using the n + m training samples. Note that the target function u∗ is a minimizer of R, and the problem can equivalently be viewed as learning u∗ satisfying (1). We will see later how this point of view is useful when considering PDEs. 3

The PIML approach approximates the expectations in (2) with empirical estimates and restricts the hypothesis space to a smaller H ⊂ F resulting in the following empirical risk minimization problem n m 1X 1 X b b R(u) = min R(u), (u(xi ) − yi )2 + (Du(zj ) − wj )2 . (3) u∈H n m i=1

j=1

In the following, we focus in particular on the setting where D is a known linear differential operator (e.g., divergence, Laplacian). The function space F is assumed to be a Sobolev space H s (X ) with s > 0 such that D is well defined from H s (X ) to L2 (X ). Note that the probability distributions of X and Z may not have the same support and, more generally, we need not observe u∗ and Du∗ at the same points. This allows us to connect the framework directly to PDE solving.

2.1

Solving Partial Differential Equations with Machine Learning

PDEs are among the most important problem classes addressed by physics-informed learning, with applications ranging from climate modeling (Kashinath et al., 2021) to the cardiovascular system (Kissas et al., 2020) to permanent magnets (Kovacs et al., 2022). In particular, boundary value problems are ubiquitous in applied mathematics: they combine differential equations governing the interior of a domain with conditions imposed on its boundary. This setting fits in our framework by allowing the laws of X and Z to be concentrated on different subsets: X is sampled from the boundary and Z is sampled from the interior. More precisely, consider the goal of estimating the solution u∗ of a boundary value problem on a bounded Lipschitz domain Ω ⊂ Rd ( Du∗ (x) = q(x) x ∈ Ω (4) u∗ (x) = h(x) x ∈ ∂Ω, where we have access to D, q and h. Then, let X = Ω and consider ρX , ρZ such that ρX (∂Ω) = 1 and ρZ (Ω) = 1. Thus, samples of X are boundary points whereas samples of Z are interior collocation points. We can then proceed as in the previous section, possibly considering the noiseless case (ϵ, η = 0) if we have perfect knowledge of q and h. In this view, minimizing (3) can be interpreted as numerically solving the PDE (4) using a random discretization. This choice of sampling is standard in machine learning theory and can be contrasted to classic deterministic discretizations in PDE, such as meshes or collocation points (see e.g. Wendland, 2004, and references therein). A simple example of a boundary-value problem is the Poisson equation on a domain Ω: ( ∆u(x) = q(x) x ∈ Ω (5) u(x) = h(x) x ∈ ∂Ω. If Ω, q and h are regular enough, we can consider F = H 2 (Ω) and take D = ∆, where ∆ : H 2 (Ω) → L2 (Ω) is the Laplace operator (Evans, 2010). Before describing our approach, we provide an in-depth overview of related problems and results in the literature.

2.2

Related Settings

Here we give an overview of the settings which are related to the one introduced above and have appeared in different contexts.

Scientific Machine Learning and PINNs. Machine learning methods for solving PDEs have attracted a lot of interest under the umbrella terms physics-informed machine learning or scientific machine learning (Rackauckas et al., 2020), especially using methods such as physics-informed neural networks (PINNs) (Raissi et al., 2017b) and Sobolev training (Czarnecki et al., 2017). PINNs have become a common approach for both forward and inverse PDE problems from fluid dynamics (Cai et al., 2021), to geophysics (Rasht-Behesht et al., 2022), and medical sciences (Sahli Costabal et al., 2020). In our setting, PINNs correspond to considering Eq. (3) while choosing neural networks as the hypothesis space H. Recent work has made significant progress in establishing learningtheoretic guarantees and error estimates for PINNs (Shin et al., 2020; Mishra and Molinaro, 2023; De Ryck et al., 2024; Zeinhofer et al., 2025; Doumèche et al., 2025c). However, fully characterizing the behavior of neural network based estimators remains challenging due to the highly non-convex nature of the underlying optimization problem. Such non-convexity is non trivial and can lead to failures in practice, as has been well observed in the literature (Krishnapriyan et al., 2021; Wang et al., 2021; Rathore et al., 2024). Operator Learning. A growing body of work focuses on operator learning, where the objective is to approximate infinite-dimensional mappings between function spaces (e.g., mapping parametric PDE coefficients). Examples include Deep Operator Networks (DeepONets) (Lu et al., 2021) and Fourier Neural Operators (FNOs) (Li et al., 2021). However, these operator-based methods require massive offline datasets of pre-computed, high-fidelity PDE solutions to learn the underlying physical mapping. In contrast, PIKS addresses the single-instance problem (analogously to standard PINNs) where the goal is to infer the solution of a specific PDE using the governing equations, boundary conditions, and sparse empirical measurements, and thus requires only small datasets. Gaussian processes and kernel methods. In parallel with the development of PINNs, Gaussian Process (GP) based methods for scientific machine learning have also been proposed (Owhadi, 2015; Raissi et al., 2017a, 2018). These approaches estimate solutions, and potentially their uncertainty, by minimizing the marginal log-likelihood (instead of the empirical risk Eq. (3)). Some extensions to non-linear PDEs have also been proposed using the Gauss-Newton algorithm (Chen et al., 2021) although convergence can only be locally guaranteed (Batlle et al., 2025). Recently Baptista et al. (2025) used the same algorithm to solve PDEs in the weak form with a non-smooth forcing term. Doumèche et al. (2024) showed that empirical risk minimization with linear differential constraints is equivalent to defining a new, physics-informed kernel and then performing kernel ridge regression, yielding theoretical evidence that physical constraints can improve convergence rates, characterized through the effective dimension of the new kernel. Such a kernel is defined by a continuous constraint (instead of pointwise evaluations as is the case in the present work) and is a priori not available in closed form, which is why subsequent works focus on practical approximations (Doumèche et al., 2025b) and fast implementations (Doumèche et al., 2025a). Hermite-Birkhoff interpolation. Despite the recent revival under the hat of physics-informed learning, the problem of learning a function from its values and its derivatives can be traced back to Hermite and Borchardt (1878) and Birkhoff (1906). It was formulated as a splines Pmproblem in the2 1D case by Kimeldorf and Wahba (1971), who considered the problem minu∈H j=1 (Lj u − wj ) with Lj a linear functional defined on an RKHS (thus including the setting of differential operators Lj u = Du(zj )). It was also considered in a multivariate setting using Radial Basis Functions (RBFs) (Zongmin, 1992). More recently Shi et al. (2010) derived learning rates for the case of Lj the gradient operator, showing that a weaker source condition is sufficient compared to standard

regression. Meshless methods for PDEs. Such methods have been widely used to approximate PDE solutions in a meshless way (i.e. without the strict mesh requirements of finite elements solvers). Kansa (1990); Fasshauer (1996) among others studied how to combine pointwise evaluations of differential constraints and boundary conditions using kernel bases. This forms part of the more general literature on scattered data approximation (Wendland, 2004), which is mainly concerned with the interpolation problem: the constraints on the data are to be enforced exactly instead of weakly as in Eq. (3). The interpolant is decomposed into a basis generated by a type of kernel known as Radial Basis Function (RBF). On the theoretical side, the analysis typically considers settings where the hypothesis space H (which is an RKHS) is the same space in which the true PDE solution lives (Franke and Schaback, 1998a,b). Convergence to this solution is proved via fill-distance techniques (Wendland, 2004). Regularization with differential operators. When W = 0 (or equivalently q = 0 in Eq. (4)), the term of (2) can be interpreted as a regularizer. Indeed, regularizers of the form R second (Du(z))2 dz have been explored in the spline smoothing literature (Wahba, 1990), in inverse problems (Hanke, 1992; Engl et al., 1996; Arridge et al., 2019) as well as in machine learning (Poggio and Girosi, 1990; Smola et al., 1998). When D is simple, the estimator can be computed in closed form using Green’s functions. The non-homogeneous case (q ̸= 0) with more general PDEs is more complex and was explored more recently as spatial regression with PDE regularization (Azzimonti et al., 2015; Sangalli, 2021; Arnone et al., 2022). In this case the closed form is usually not available, but the estimator can be computed using finite element methods. Manifold Regularization. In semi-supervised learning, the differential regularizer is unknown and must be approximated from available data (Zhu et al., 2003; Zhou and Schölkopf, 2005). This leads to manifold regularization approaches. For instance, Belkin et al. (2006) use a gradient-based regularizer of the form ∥∇M u∥2 to penalize deviations of u from the data manifold M , which itself is estimated from data using the graph Laplacian. See also Slepcev and Thorpe (2019); Cabannes et al. (2021) for recent results in this direction. Misspecified kernel methods. The approximation capabilities of kernel methods depend on the size of the RKHS H, which itself depends on the kernel choice. Kernels for which H is dense in the space F in which the target lives (typically an Lp space) are called universal. Universality of many common kernels has been well studied (Micchelli et al., 2006; Sriperumbudur et al., 2011; Simon-Gabriel and Schölkopf, 2018), and for kernel ridge regression, universality implies asymptotic consistency in the misspecified setting (De Vito et al., 2005b). In the case of Sobolev RKHS, finer results exist in the misspecified setting, with convergence rates that depend on the exact smoothness of the target, characterized with a source condition (Steinwart et al., 2009; Lin et al., 2020; Fischer and Steinwart, 2020; Zhang et al., 2023). Convergence results have also been established for Gaussian processes in misspecified settings (Wynne et al., 2021; Wang and Jing, 2022). In the physics-informed machine learning context, the study of the misspecified setting is much more limited. Building on the work of Narcowich et al. (2006) who analyze the error of RBF interpolation in the misspecified setting, Schräder and Wendland (2012) apply these results to a PDE context, and establish error estimates when the solution is less smooth than the considered hypothesis space. Similarly to all the works that belong to the RBF interpolation literature, the results are established for a

deterministic and noiseless set of input points, and the estimates are expressed in terms of filldistance. Furthermore, the target function can be outside of the considered RKHS, but must still belong to a Sobolev RKHS, i.e. a space H s (Ω) with s > d/2. In contrast, the present work considers random design with noise, and the space F in which the target lives need not be an RKHS. More recently, Baptista et al. (2025) studied a kernel-based collocation method for PDEs with rough solutions, which could be outside the RKHS, but the method considers the PDE in a weak form, and relies on test functions, which is a different setting than the present one. Provided with the above discussion we next describe the approach we consider and analyze.

Physics-informed kernel methods (PIKS)

In this section, we introduce the PIKS estimator, specify some initial assumptions on D, and introduce the regularized problem and the closed-form expression for PIKS. We consider the hypothesis space H to be a reproducing kernel Hilbert space (RKHS) of functions from X to R. An RKHS is defined by a kernel function K : X × X → R such that for all x ∈ X , the function Kx : y 7→ K(x, y) belongs to H and satisfies the reproducing property ∀u ∈ H,

⟨u, Kx ⟩H = u(x).

(6)

For all the computations in this section, we need the following assumption. Assumption 1. For all x ∈ X , the map H → R, u 7→ Du(x) is a bounded functional. Under Assumption 1, the Riesz representation theorem guarantees that for all x ∈ X there exists a representer KxD ∈ H such that ∀u ∈ H,

⟨u, KxD ⟩H = Du(x).

(7)

Assumption 1 is formulated in a general setting where D can be any linear operator. We can however make it more concrete in the case of differential operators, for which the reproducing property is well studied (Zhou, 2008). In such a case, the smoothness of the kernel is enough to guarantee that Assumption 1 holds, as shown by the following result. Lemma 1. Assume that X is compact and satisfies Int(X ) = X , and consider a linear differential operator of the form X (8) D= cα ∂ α , |α|≤s

for some integer s ≥ 1 and where the cα : X → R are continuous coefficient functions. If K ∈ C 2s (X × X ), then Assumption 1 holds. Furthermore the representer KxD is available in closed form: KxD : y 7→ D1 K(x, y).

(9)

Proposition 1 is a simple consequence of a result from Zhou (2008, Theorem 1), and is proved in the appendix (see Proposition 8 in Section A.2). Then to verify Assumption 1, we simply ∥x−y∥2 need a smooth-enough kernel. For example the Gaussian kernel K(x, y) = exp(− 2σ2 2 ) is C ∞ over Rd × Rd , while the Matérn kernel Kν (x, y) of order ν > 0 is C 2s over Rd × Rd as soon as ν > s. Proposition 1 also provides the expression of the representer, which is key for computing the estimator in closed form.

If the kernel K satisfies Assumption 1, we define the PIKS estimator u bλ as the minimizer of the regularized physics-informed empirical risk: bλ (u), u bλ := arg min R u∈H

bλ (u) = 1 R n

n X i=1

1 X (u(xi ) − yi ) + (Du(zj ) − wj )2 + λ∥u∥2H . m 2

(10)

j=1

For any λ > 0, the regularized loss (10) is strongly convex, so the PIKS estimator is well-defined and unique, and one can show (see Proposition 9 in Section B) that it can be decomposed as u bλ =

n X

αi Kxi +

i=1

m X

βj KzDj ,

(11)

j=1

where (α, β) = (K + λJ)−1 Y ∈ Rn+m . In this formulation K is a (n + m) × (n + m) block kernel matrix, J a diagonal regularization matrix, and Y the regressed variables:   y1  ..     .  (   A C    yn  n if 1 ≤ i ≤ n    Ji,i = Y= K= (12)  w1  .  m if n + 1 ≤ i ≤ n + m   ⊤ C B  ..   .  wm The blocks of the kernel matrix are A ∈ Rn×n ,

B ∈ Rm×m , C ∈ Rn×m ,

Ai,i′ = ⟨Kxi , Kxi′ ⟩H ,

∀i, i′ ∈ J1, nK;

Ci,j ′ = ⟨Kxi , KzDj ′ ⟩H ,

∀i ∈ J1, nK, ∀j ′ ∈ J1, mK.

Bj,j ′ = ⟨KzDj , KzDj ′ ⟩H ,

∀j, j ′ ∈ J1, mK;

where we denoted J1, nK = {1, . . . , n}. Note that ⟨Kxi , Kxi′ ⟩H = K(xi , xi′ ), and in the case of differential operators discussed in Proposition 1, we also have ⟨KzDj , KzDj ′ ⟩H = D1 D2 K(zj , zj ′ ) and

⟨Kxi , KzDj ′ ⟩H = D2 K(xi , zj ′ ) where D1 denotes the differential operator with respect to the first variable of the kernel and D2 with respect to the second one. The structure of K arises from the presence of two types of data, (xi , yi )ni=1 and (zj , wj )m j=1 . Diagonal blocks A and B are the kernel matrices of the two separate datasets, and C represents the cross terms. Computationally the PIKS estimator (11) requires storing and later inverting an (n + m) × (n + m) matrix with a cost of O((n + m)2 ) space and O((n + m)3 ) time units. While this cost is high in general, there exists a rich literature on approximations which greatly reduce the computational complexity of kernel methods without compromising on accuracy (Rahimi and Recht, 2007; Rudi et al., 2015). In the particular case of structured matrices involving derivatives such as (12), we can notably cite the works of Eriksson et al. (2018); Padidar et al. (2021); De Roos et al. (2021) in the case of Hermite-Birkhoff interpolation and of Chen et al. (2025) in the case of PDEs. Remark 1. As we described in Section 2.2, hybrid regression settings in RKHS mixing different types of linear functionals (and in particular differential operators) were studied at least since Kimeldorf and Wahba (1971). Solving such problems leads to block matrices of the form (12), which were observed in Hermite-Birkhoff problems (Zongmin, 1992), in RBF collocation methods to approximate PDE solutions (Fasshauer, 1996; Franke and Schaback, 1998b; Wendland, 2004), as well as in recent uses of Gaussian processes for PDEs (Raissi et al., 2017a; Chen et al., 2021, 2025). 8

Theoretical analysis of the PIKS estimator

In this section we derive theoretical results about the asymptotic convergence of the PIKS estimator. Usual kernel analyses rely on the assumption that the model is well-specified, that is, that the target function u∗ belongs to the same RKHS from which the estimator is taken. However, in practice, one does not always know the exact regularity of the target function. Even worse, the target function might lack the smoothness required to belong to any RKHS. For instance, it is well known that a Sobolev space H s (Ω) is an RKHS if and only if s > d/2, a condition that is harder to satisfy in high dimensions. Such a situation arises naturally in PDE problems: for a second-order elliptic equation with sufficiently regular coefficients and L2 data, elliptic regularity theory guarantees that the solution belongs to H 2 (Ω) (Evans, 2010), but in general to no smoother Sobolev space. Since H 2 (Ω) is an RKHS only when d ≤ 3, the natural regularity class of the solution fails to be an RKHS as soon as d ≥ 4. This motivates studying the misspecified setting where u∗ ̸∈ H. Next, we first introduce the main technical assumptions in Section 4.1, then we state the main result in Section 4.2, which we illustrate in the case of differential operators and Sobolev spaces in Section 4.3, and in the particular case of PDE settings in Section 4.4. Finally, in Section 4.5, we derive convergence rates, which we illustrate in the example of the Laplacian on Sobolev spaces of periodic functions.

4.1

Functional assumptions

We collect here the functional assumptions used to prove our main result in Theorem 1. They are formulated for general Hilbert spaces of functions F and G. Later, in Section 4.3, we specialize these assumptions to Sobolev spaces and give concrete sufficient conditions under which they can be verified. We start by providing further details on the setting introduced in Section 2. Let us consider a compact domain X ⊂ Rd and a target u∗ ∈ F where F is a Hilbert space of functions embedded in L2 (X ) (endowed with the Lebesgue measure of Rd ). We consider a bounded linear operator D : F → G, where G is another Hilbert space of functions embedded in L2 (X ). Since the elements of F and G are only defined almost everywhere, we need the following compatibility condition: Assumption 2. The embeddings F ,→ L2 (ρX ) and G ,→ L2 (ρZ ) are well-defined and bounded. Assumption 2 ensures that the random variables u∗ (X) and Du∗ (Z) introduced in Section 2 are well defined. Indeed, suppose for example that u∗ ∈ L2 (X ). Such a function is only defined almost everywhere with respect to the Lebesgue measure, so for instance if ρX is equal to the Dirac distribution δx (that is, we only sample at a single, deterministic location x ∈ X ), then the variable u∗ (X) is not well defined, as the value of u∗ (x) is not uniquely defined a priori. A similar observation holds for Du∗ (Z). If instead, for this same u∗ , ρX and ρZ are absolutely continuous with respect to the Lebesgue measure on X , with bounded densities, we easily see that Assumption 2 holds. Note that Assumption 2 also covers less trivial settings. One important example is the one discussed in Sections 2.1 and 4.4 where ρX is only supported on the boundary of X . Then, the part of Assumption 2 concerning ρX is provided not by the embedding of L2 (X ) into L2 (ρX ) but rather by the trace theorems (Evans, 2010), as we discuss in Section 4.4. We now introduce the universality assumption, which is central to our study. Assumption 3 (Universality). The RKHS H is densely embedded in F . This assumption is standard in the kernel literature when F = L2 (X ), and is satisfied for a broad class of kernels which are referred to as universal (Micchelli et al., 2006). It guarantees sufficient 9

flexibility to learn the target function in misspecified settings, that is when u∗ ̸∈ H. However, in the physics-informed setting, we consider more structured spaces F, so the density of H in L2 is not enough, and Assumption 3 is thus a stronger requirement than classical universality. For Sobolev spaces F = H s (X ), a sufficient condition for Assumption 3 to hold is C0s -universality, as we discuss in Section 4.3, and it is satisfied for instance by Gaussian and Matérn kernels. We now make two boundedness assumptions on the data-generating process and on the features in H which are used to guarantee the validity of concentration inequalities. Assumption 4 (Bounded data). We have u∗ ∈ L∞ (ρX ) and Du∗ ∈ L∞ (ρZ ). Furthermore, both noise random variables ϵ and η are bounded almost surely. Assumption 4 restricts the class of admissible targets and noise distributions. Its role is to ensure that the empirical quantities appearing in our analysis are uniformly bounded, so that Hoeffdingtype concentration inequalities apply. Such boundedness assumptions are standard in convergence analyses of kernel methods. They can in principle be relaxed — for example, by imposing tail or moment conditions and using Bernstein-type inequalities, truncation, or other refined concentration tools — but pursuing these extensions is beyond the scope of the present work. Assumption 5 (Bounded features). There exist κ, κD > 0 such that for all u ∈ H, for all x ∈ X , |u(x)| ≤ κ∥u∥H ;

|Du(x)| ≤ κD ∥u∥H .

(13)

Under Assumption 1, we can define the kernel K D (x, y) := ⟨KxD , KyD ⟩H , and (13) can be reformulated as having, for all x ∈ X , K D (x, x) ≤ κ2D .

K(x, x) ≤ κ2 ;

(14)

The formulation (14) might be more familiar to the reader. Assumption 5 can typically be obtained as a consequence of the smoothness of the kernels and the compactness of X . It is in particular the case in the setting of Proposition 1 (for more details, see Proposition 32 in Section D).

4.2

Main result

Consider the physics-informed setting described in Section 2, as well as in the previous section. Consider the PIKS estimator u bλ ∈ H defined in (10). The following result characterizes the asymptotic behavior of the estimator as the dataset sizes increase. Theorem 1. Under Assumptions 1 to 5, the PIKS estimator is a universal learner: for any regularizing sequence (λn,m ) such that λn,m → 0,

log N →0 λ3n,m N

as n, m → ∞,

where N = min(n, m), almost surely, the estimator u bλn,m satisfies  ∥b uλn,m − u∗ ∥L2 (ρX ) −→ 0 n,m→∞

∥Db uλn,m − Du∗ ∥L2 (ρZ )

−→

n,m→∞

(15)

(16)

This theorem, which is proved as Proposition 29 in the appendix, shows that the PIKS estimator is able to learn the function’s values and satisfy the physical constraint at the same time. Here, the main technical difficulty lies in showing that it holds under weak assumptions on u∗ . As 10

discussed in Section 2.2, theoretical analyses usually assume u∗ belongs to the same RKHS H from which u bλ is taken. Indeed, in the well-specified setting, one can typically prove the stronger convergence ∥b uλ − u∗ ∥H → 0, from which we can derive both convergences ∥b uλ − u∗ ∥L2 (ρX ) → 0 and ∗ ∥Db uλ − Du ∥L2 (ρZ ) → 0 as direct consequences. In the misspecified setting instead, the standard analysis (see e.g. De Vito et al., 2005b) only guarantees the L2 convergence ∥b uλ − u∗ ∥L2 (ρX ) → 0 for the classical KRR estimator, which does not imply the convergence ∥Db uλ − Du∗ ∥L2 (ρZ ) → 0. Theorem 1 extends the analysis to the physics-informed setting, by showing that universality guarantees both convergences at the same time, without requiring u∗ ∈ H. Remark: Convergence with n → ∞ and m finite. An observation which could arise from looking at Equation (16) is that both dataset sizes (n and m) must tend to infinity to guarantee convergence, even for the function values alone; in other words, the two convergences are not decoupled. It is possible to decouple them by introducing a scaling parameter to the physics term. If this parameter were set to decay to 0 as n → ∞ while m remains finite, the PIKS estimator would become asymptotically equivalent to a classical KRR estimator, thereby recovering the standard L2 convergence on the values, but losing the physical consistency. However, this is not the focus of the present paper, and scaling parameters were left out to maintain a simplified analysis.

4.3

The case of differential operators on Sobolev spaces

The assumptions in Section 4.1 and the result in Section 4.2 are stated for a general linear operator D and function spaces F, G. Here, we specialize them to the case where D is a differential operator defined on a Sobolev space. In particular, we discuss the implications of the universality assumption, and analyse more concretely the theorem’s applications in this setting. Consider X = Ω with Ω a bounded Lipschitz domain, u∗ ∈ F = H s (Ω), for some integer s > 0, and G = L2 (Ω). Assume that D is a linear differential operator of order s as defined in Eq. (8). We first check that the operator D indeed defines a bounded operator H s (Ω) → L2 (Ω) (see Proposition 30 in Section D.1). For Theorem 1 to apply, we further need to check Assumptions 1 to 5. Assumption 4 is a standard boundedness assumption that we take independently from the setting. All other assumptions can be verified as consequences of the setting and the kernel choice, as we see below. Assume that K : X × X → R is a C 2s kernel. We already established in Section 3 that in such a case, Assumption 1 holds. It is also straightforward to check that Assumption 5 holds, as a consequence of the compactness of X (see Proposition 32 in Section D.1). Assumption 3 instead requires more care, as it requires the RKHS to be rich enough to approximate not only functions, but also their derivatives up to order s — a requirement that is stronger than the standard L2 universality. A sufficient condition, which holds for common kernels such as the Gaussian one, is C0s universality. Definition 1 (C0s universality). Let U be an open subset of Rd and let C0s (U ) be the set of continuous functions f : U → R such that f and all its derivatives up to order s tend to 0 at infinity (i.e. for any ϵ > 0 and α ∈ Nd , |α| ≤ s, there exists a compact Kϵ ⊂ U such that for all x ∈ U \Kϵ , |∂ α f (x)| ≤ ϵ). We say that a kernel K is C0s -universal on U if for any f ∈ C0s (U ), for any ϵ > 0, there exists fH in the RKHS H associated to K such that for any α ∈ Nd satisfying |α| ≤ s, sup |∂ α f (x) − ∂ α fH (x)| ≤ ϵ.

x∈U

The following lemma shows that if the kernel considered is C0s -universal on Rd , then Assumption 3 holds. 11

Lemma 2. Assume that X = Ω, with Ω a Lipschitz domain, and that F = H s (X ) = H s (Ω). If K is the restriction to X of a C0s -universal kernel on Rd , then Assumption 3 holds. We prove Proposition 2 in Section D.1, where it is restated as Proposition 31. The proof proceeds by extending a given function u ∈ H s (Ω) to Rd , approximating the extension by sufficiently regular compactly supported functions, and then invoking C0s universality. Restricting the resulting approximants to X yields the desired approximation of u. C0s universality was discussed by Simon-Gabriel and Schölkopf (2018), in a paper that characterizes it in particular for translation-invariant kernels on Rd (which are of the form K(x, y) = Φ(x−y)). Such kernels are C0s -universal if and only if they are C s,s (which is slightly weaker than the C 2s requirement of Proposition 1, that we consider in the present work — we refer the reader to SimonGabriel and Schölkopf (2018) for the precise definition) and the Fourier transform of Φ has full support. For instance, this is verified for the Gaussian kernel and the Matérn kernel of index ν > s. In-domain sampling. The only missing element to apply Theorem 1 is Assumption 2. Such assumption depends on ρX , ρZ and thus on the sampling setting. A natural choice is uniform sampling in Ω, as we study here; we will discuss another interesting choice in the next section. Assume now that ρX = ρZ = Unif(Ω). Then, we have L2 (ρX ) = L2 (ρZ ) ∼ L2 (Ω). The embedding of G into L2 (ρZ ) is then trivial, and the embedding of F into L2 (ρX ) is a consequence of the canonical embedding of H s (Ω) into L2 (Ω). We thus see that Assumption 2 is satisfied. Note that, more generally, it would be satisfied for any distributions with bounded density on Ω— and that Theorem 1 applies. Theorem 1 tells us that as the data grows to infinity, we get the convergences u b L2 (ρ

L2 (ρX )

−→

u∗

and Db u −→ Du∗ , which are (up to normalization) convergences in L2 (Ω) since ρX and ρZ are both uniform. The benefit of the PIKS estimator compared to standard kernel regression in this setting lies in the ability to learn with two kinds of data, and to guarantee that we learn both u∗ and Du∗ simultaneously, as the latter does not follow from the former in general. For instance, if we know that the target satisfies a PDE Du∗ = q, then, with enough data, with PIKS we will have Db u ≈ q (in the L2 sense), i.e. the estimator is physically consistent. In contrast, with kernel ridge regression one can only establish the L2 convergence of u b to u∗ , without a corresponding guarantee on physical consistency. The benefit of these two simultaneous convergences becomes even clearer in some cases where we have stability of the PDE, as we study in the next paragraph. This setting is related to the earlier fundamental work of Doumèche et al. (2024, 2025b) on physics-informed learning. The authors study a physics-informed constraint of the form ∥Du∥L2 in addition to the standard Tikhonov regularization, in a Sobolev RKHS. In our setting, this corresponds to the homogeneous case, Du∗ = 0, together with an infinite amount of physical information, namely m → ∞. They show that their estimator is equivalent to classical KRR with a modified, physics-informed kernel. Unlike the present work, they consider the well-specified setting u∗ ∈ H, in which the asymptotic convergence ∥b u − u∗ ∥H → 0 is straightforward to establish; the two convergence results in Theorem 1 then follow as consequences. Their focus is different, however: they study how incorporating the physics accelerates the convergence of u b to u∗ in the L2 sense, which is beyond the scope of the present work. Stronger convergence results in the elliptic case. For some differential operators, the convergence Db u → Du∗ provided by Theorem 1 can be used to deduce convergence of u b to u∗ in a stronger topology than that of L2 . This relies on PDE regularity estimates and depends strongly on the operator D. We illustrate this in the elliptic case. 12

We retain the Sobolev setting considered above, with uniform sampling in Ω. We further assume that s = 2 and that D is a uniformly elliptic operator satisfying the regularity assumptions of Section D.2.1. Ellipticity guarantees us (Evans, 2010, Section 6.3) that for any q ∈ L2 (Ω) and any weak solution u of the PDE Du = q on Ω, we have the estimate ∥u∥H 2 (V ) ≤ CV (∥u∥L2 (Ω) + ∥q∥L2 (Ω) ),

(17)

for any open V such that V ⊂ Ω. The regularity estimate (17) allows us to translate the two L2 convergences of Theorem 1 into a convergence of u b to the target u∗ in a stronger, Sobolev sense, as we see in the following result. Corollary 3. Under the preceding assumptions, for any sequence (λn,m ) satisfying (15), we have that almost surely, for any open V such that V ⊂ Ω, ∥b uλn,m − u∗ ∥H 2 (V )

−→

n,m→∞

Proposition 3 is restated as Proposition 37 and proved in Section D.2. Proposition 3 shows a benefit of physics-informed learning in the elliptic setting: while for a classical KRR estimator (using only the classical data (xi , yi )), theory typically only guarantees L2 convergence, here we obtain the asymptotic convergence of the PIKS estimator to the target in the H 2 sense on all compactly embedded domains.

4.4

Data on the boundary: solving PDEs

In the previous section, we focused on uniform sampling in Ω for both X and Z. Another important sampling setting is the one described in Section 2.1, where (xi )ni=1 are sampled on the domain’s boundary and (zj )m j=1 in its interior. This happens when instead of having data measurements of u∗ in the domain, we have a boundary condition at our disposal. This gets us closer to typical PDE settings, where we are trying to solve a boundary value problem ( Du∗ (x) = q(x) x ∈ Ω, (18) u∗ (x) = h(x) x ∈ ∂Ω. Let us indeed consider, same as Section 4.3, that Ω is a bounded Lipschitz domain, u∗ ∈ F = H s (Ω), G = L2 (Ω), and D : F → G is a differential operator defined as in (8). Consider this time that ρX = Unif(∂Ω) and ρZ = Unif(Ω). Let K ∈ C 2s (X × X ) be the restriction of a C0s -universal kernel on Rd . The only difference is the sampling setting, and we know from the previous section that Assumptions 1, 3 and 5 hold, and again we can assume that Assumption 4 holds independently. In order to apply Theorem 1, we only need to check Assumption 2. For that, we observe that while the first equality in (18) must be understood as holding almost everywhere, with respect to the Lebesgue measure, the second equality does not make sense a priori since ∂Ω has Lebesgue measure 0 and u∗ is only defined Lebesgue almost everywhere. Nevertheless it can be made meaningful with the trace operator (Evans, 2010, Sec. 5.5) which is a bounded operator T : H 1 (Ω) → H 1/2 (∂Ω) that coincides with the restriction to the boundary u 7→ u|∂Ω when u is a continuous function over Ω. The trace operator is the standard way of defining boundary conditions in PDE theory, as soon as we consider weak solutions. Composing T with the canonical Sobolev embeddings H s (Ω) ,→ H 1 (Ω),

H 1/2 (∂Ω) ,→ L2 (∂Ω), and with the bounded map L2 (∂Ω) → L2 (ρX ), we see that Assumption 2 holds. Hence, Theorem 1 applies. The conclusions of Theorem 1 have a different meaning than for in-domain sampling in SecL2 (ρX )

tion 4.3. Since ρX is supported on ∂Ω, the convergence u b −→ u∗ is merely a convergence on the boundary; by itself, it tells us nothing about what happens in Ω. The two convergences of Theorem 1 mean that the PIKS estimator approximately satisfies the PDE (18) when the data grow to infinity. For a strong convergence to u∗ on Ω, one further needs an appropriate stability estimate for the boundary-value problem, which is PDE-dependent. We give an example of stronger convergence results at the end of the current section. This setting has received considerably more attention in the literature, as it corresponds to a classical PDE formulation. In particular, the use of kernels to numerically solve linear PDEs by prescribing function values at given boundary points and the values of the differential operator at given domain points is not new and has been studied in the framework of interpolation with RBF functions (Wendland, 2004). In such literature, the points at which the differential operator is evaluated are deterministic, referred to as collocation points. The convergence analysis of these methods is then carried out using the notion of fill-distance, that is, the biggest distance that exists between a point and its closest neighbor. The analysis typically takes place in a noiseless setting, with the true solution of the PDE belonging to the RKHS defined by the chosen RBF. Our analysis addresses a different and complementary regime, building on KRR theory in misspecified settings. We assume that the sampling points xi and zj are drawn randomly from prescribed distributions, allow the observations to be noisy, and do not require the target u∗ to belong to the RKHS H. This makes it possible to derive statistical guarantees that account jointly for sampling variability, observation noise, and model misspecification. In turn, the resulting bounds are distribution-dependent and typically control an average error, rather than providing deterministic guarantees tied to the geometric coverage of a particular set of collocation points. The more recent work (Chen et al., 2021) studies nonlinear operators. They formulate the problem as a nested optimization and propose an iterative method (with a Gauss-Newton algorithm). On the theoretical side, they consider the well-specified setting as well, and prove convergence of the second member of the PDE. Similar to existing results on PINNs, they must assume the convergence of the optimization process, as it is difficult to obtain guarantees in the nonlinear case. Stronger convergence results. Assume now that s = 2 and D is uniformly elliptic with smooth coefficients over Ω (see the beginning of Section D.3.1 for the precise assumptions). Further assume that the boundary ∂Ω is a smooth (d − 1)-dimensional manifold, Ω being locally on one side of ∂Ω. We assume that 0 is not a Dirichlet eigenvalue for the operator D in Ω, so that using Fredholm alternative (Evans, 2010), we can guarantee that for any (q, h) ∈ L2 (Ω) × H 3/2 (∂Ω), there exists a unique solution u ∈ H 2 (Ω) to (18). This is the case for instance for the Laplacian D = −∆. Then, ellipticity and the regularity of the boundary guarantee us (see Lions and Magenes, 2012) that for any (q, h) ∈ L2 (Ω) × H 3/2 (∂Ω), for any u ∈ H 2 (Ω) solution of (18), we have the estimate ∥u∥H 1/2 (Ω) ≤ c ∥h∥L2 (∂Ω) + ∥q∥Ξ−3/2 (Ω) , (19) where c > 0 and ∥ · ∥Ξ−3/2 (Ω) is a norm weaker than ∥ · ∥L2 (Ω) (see Lions and Magenes, 2012, for more details). The estimate (19) allows us to translate the two L2 convergences of Theorem 1 into an H 1/2 convergence of u bλn,m to the target u∗ , as we see in the following result.

Corollary 4. Under the preceding assumptions, for any sequence (λn,m ) satisfying (15), almost surely, ∥b uλn,m − u∗ ∥H 1/2 (Ω) −→ 0, (20) n,m→∞

and, for every open V such that V ⊂ Ω, almost surely, ∥b uλn,m − u∗ ∥H 2 (V )

−→

n,m→∞

(21)

Proposition 4 is proved in Section D.3.1 under the form of two results, Proposition 42 and Proposition 43. Let us compare this result with Theorem 1: there, since ρX is supported on ∂Ω, the L2 convergence of u bλn,m to u∗ only happens on the boundary. Instead, (20) shows that, in the elliptic case with smooth boundary and uniqueness of the solution, we have in fact H 1/2 convergence on the whole domain Ω. Then, analogously to Proposition 3, (21) shows that we have H 2 convergence on all compactly embedded domains.

4.5

Convergence rates

Theorem 1 establishes asymptotic convergence, without any indication on the convergence speed. It is actually possible to get convergence rates if one further assumes a source condition, following other works in the KRR literature (see e.g. De Vito et al., 2005a; Steinwart et al., 2009; Blanchard and Mücke, 2018; Lin et al., 2020; Fischer and Steinwart, 2020; Zhang et al., 2023). To state the source condition, let us introduce some notation. Assumption 2 allows us to consider the bounded operator F −→ L2 (ρX ) × L2 (ρZ ) Aρ : (22) u 7−→ (u, Du), where we endow L2 (ρX ) × L2 (ρZ ) with the product norm defined by ∥(f, g)∥2ρ := ∥f ∥2L2 (ρX ) + ∥g∥2L2 (ρZ ) . Let us denote (h, q) := Aρ u∗ . In the boundary-value setting of Section 2.1 this is consistent with the notation, since h = u∗ |∂Ω and q = Du∗ ; in the in-domain setting h is simply u∗ viewed in L2 (ρX ). We can compose Aρ with the embedding i : H → F to define Aρ = Aρ ◦ i,

(23)

which allows us to write for any u ∈ H, ∥Aρ u − (h, q)∥2ρ = ∥u − h∥2L2 (ρX ) + ∥Du − q∥2L2 (ρZ ) . We can now define L := Aρ A∗ρ : L2 (ρX ) × L2 (ρZ ) → L2 (ρX ) × L2 (ρZ ).

(24)

The operator L is a composite integral operator: for (f, g) ∈ L2 (ρX ) × L2 (ρZ ) and (x, z) ∈ X 2 , it can be expressed as Z Z L(f, g)(x, z) = ⟨Ky , Kx ⟩H f (y)dρX (y) + ⟨KtD , Kx ⟩H g(t)dρZ (t) , Z X Z X D D D ⟨Ky , Kz ⟩H f (y)dρX (y) + ⟨Kt , Kz ⟩H g(t)dρZ (t) . X

We now have all the elements to introduce a source condition adapted to this physics informed setting. 15

Assumption 6 (Physics-informed source condition). There exists r ∈ (0, 1] such that (h, q) ∈ ran Lr , i.e. there exists (fe, ge) ∈ L2 (ρX ) × L2 (ρZ ) such that (h, q) = Lr (fe, ge). Given the above assumption we can derive explicit convergence rates. 1 Theorem 2. Let δ ∈ (0, 21 ). Under Assumptions 1 to 6, if λ ≥ C max n1 log nδ , m log m δ , then with probability at least 1 − δ, the PIKS estimator satisfies p p ln(1/δ) ln(1/δ) √ + 1−r̃ √ + λr , ∥b uλ − h∥L2 (ρX ) + ∥Db uλ − q∥L2 (ρZ ) ≲ (25) 1−r̃ λ n λ m where r̃ = min(r, 1/2). Theorem 2 is restated and proved as Proposition 24 in Section C.5. If we consider the minimum between n and m and select λ accordingly, the following corollary gives a simplified rate. 1

Corollary 5. Let us set λ = N −1/2 if r ≤ 1/2, and λ = N − 2r+1 if r > 1/2, where N := min(n, m). 1 Under Assumption 6, if δ ∈ (0, 12 ) is such that λ ≥ C max n1 log nδ , m log m δ , then with probability at least 1 − δ, we have p ∥b uλ − h∥L2 (ρX ) + ∥Db uλ − q∥L2 (ρZ ) ≲ ln(1/δ)N −r/2 if r ≤ 1/2, (26) p r − 2r+1 ln(1/δ)N if r > 1/2. (27) ∥b uλ − h∥L2 (ρX ) + ∥Db uλ − q∥L2 (ρZ ) ≲ Proposition 5 is restated and proved as Proposition 25 in Section C.5. Using integral operator techniques and source conditions like Assumption 6 is classical in KRR analysis; the distinctive feature here is that the rates depend on two indices n and m. For the convergence to be guaranteed, both need to tend to infinity. The practicality of the rates provided in (25) and (26) depends on the ability to interpret the source condition, which itself depends on the operator L. Such an operator is a composite integral operator which is more complex than the corresponding object in classical KRR. We provide below an example where the source condition is interpretable, while leaving a more general analysis for future work. Example 1 (Laplacian on periodic functions). Consider X = [0, 1]d , and denote, for any τ > 0, τ Hper (X ) ⊂ H τ (X )

the Sobolev space of order τ of periodic functions on X . Such spaces can be described with Fourier 2 (X ), G = L2 (X ) and D = ∆ = series, as we detail in Section D.4. Let us consider F = Hper Pd 2 ∂ τ i=1 (∂xi )2 the Laplacian, which defines a continuous operator from F to G. Let H = Hper (X ),

with τ > d/2 + 2. Let us consider ρX = ρZ = Unif([0, 1]d ) (which corresponds to the in-domain sampling setting of Section 4.3). In particular L2 (ρX ) = L2 (ρZ ) = L2 ([0, 1]d ). One can show that Assumptions 1 to 3 and 5 hold (see Proposition 44 in Section D.4). Consider furthermore a target u∗ ∈ F, and a noise model satisfying Assumption 4. Then Theorem 1 guarantees the asymptotic convergence. Furthermore, in such a setting, the source condition is interpretable in terms of smoothness of the target function u∗ , as shown by the following proposition. Proposition 6. Let u∗ ∈ F, and (h, q) = Aρ u∗ . Then, Assumption 6 is satisfied, with r ∈ (0, 1] σr (X ), with σ = 2 + 2r(τ − 2). In particular, for r = 1 , Assumption 6 is if and only if u∗ ∈ Hper r 2 equivalent to u∗ ∈ H, which corresponds to the well-specified setting. Proposition 6 is proved in Section D.4. It is classical in KRR to characterize source conditions in Sobolev RKHS as Sobolev smoothness of the target u∗ of a given order (see e.g. Fischer and Steinwart, 2020). We see with Proposition 6 that the same is possible in the case of the Laplacian on periodic functions. This gives a clear interpretation to Assumption 6. 16

Function value RMSE

100

only f f and ∇f

10−1

10−2

27 28 29 210 211 Num. function points

Figure 1: Improved learning accuracy with derivative data. A dataset of function values (in red) is augmented by gradient (in blue) information to improve accuracy.

Figure 2: Errors on the 1D wave equation. The colorscale is logarithmic and at this scale of errors only the standard PINN produces a solution which looks qualitatively incorrect.

Empirical verification

In this section we demonstrate the performance of PIKS in two scenarios: the first one is a supervised learning problem in which we show that gradient information helps improve the underlying model’s accuracy, decreasing the error on both the function values and derivatives. Such a setting is sometimes referred to as Sobolev training or Hermite learning. In the second scenario, we use PIKS as a PDE solver and compare it to FEM solvers, PINNs and alternative kernel-based PDE solvers on three linear PDEs. As discussed earlier, PIKS in this setting is close to classical meshless methods with kernels (Wendland, 2004), and our results emphasize that there are many cases in which kernel regression with a RBF kernel is competitive with more complex methods involving neural networks or kernel functions. We can furthermore show (see Fig. 5) that indeed it is possible to use kernel methods for learning in the misspecified setting.

5.1

Learning with gradient data

In traditional machine learning settings, where the goal is to learn a function f : X → Y given a dataset of pairs (x, y) ∈ X × Y, PIKS can be used whenever additional information about linear transformations of f is known. We demonstrate the effect of incorporating such information with a simple example of a smooth 2D function f (x) = sin(πx1 ) · sin(πx2 ) + 4 sin(4πx1 ) · sin(4πx2 ) with additive Gaussian noise. We compare the generalization error of KRR with that of PIKS which has additional access to derivative information. The number of training samples from f is fixed to n, and the additional samples from ∂x1 f and ∂x2 f are also fixed to n. The RMSE obtained as a function of n and of the available data is shown in Fig. 1. Gradient data is especially useful for this function: having access to n function plus n gradient points is better than having access to 2n function points.

5.2

PDE settings

We showcase PIKS for solving PDEs with three examples. The first two focus on a well-known failure mode of PINNs: they struggle to learn functions with a large range of frequency compo-

Table 1: Estimator comparison on the convection equation. Data for PINNs taken from Krishnapriyan et al. (2021), for PIKL from Doumèche et al. (2024). Method Vanilla PINN Curriculum PINN PIKL estimator PIKS estimator

Rel. RMSE (β = 30) 8.97 × 10−1 2.02 × 10−2 0.91 × 10−7 2.18 × 10−6

Table 2: Solver comparison on the 1D wave equation. Noisy initial conditions use σ = 0.1. Standard deviation is computed on 10 repetitions. Method

Rel. RMSE

Rel. RMSE (noisy IC)

Euler CN PINN PINN+NTK PIKL PIKS

6 × 10−4

(1.3 ± 0.1) × 10−1 (6.9 ± 0.5) × 10−2 (4.1 ± 0.2) × 10−1 (1.4 ± 0.5) × 10−2 (2.8 ± 1.2) × 10−2 (1.4 ± 0.2) × 10−2

6.4 × 10−2

(4.1 ± 0.3) × 10−1 (4.9 ± 3.6) × 10−3 (8.7 ± 0.1) × 10−4 (1.0 ± 0.4) × 10−6

nents (Rahaman et al., 2019), comparing also against other kernel estimators. In the last we focus on misspecification and compare against a FEM solver. Convection equation Following the setup of Krishnapriyan et al. (2021), we take the 1D convection equation with smoothness controlled by parameter β. This PDE, whose solution is f (t, x) = sin(x − βt), is defined as   ∀x ∈ [0, 2π], t ∈ [0, 1], ∂t f (t, x) + β∂x f (t, x) = 0 (28) ∀x ∈ [0, 2π], f (0, x) = sin(x)   ∀t ∈ [0, 1], f (t, 0) = f (t, 2π). Krishnapriyan et al. (2021) show that PINN performance degrades with increasing β and propose curriculum training to improve on this, while Doumèche et al. (2024) show that kernel-based PIKL is accurate to numerical precision. In Table 1, we show that, even with a standard RBF kernel and the same number of boundary points (100), PIKS achieves a relative RMSE nearly four orders of magnitude below that of curriculum-trained PINNs and approaches the accuracy of PIKL, which retains an edge in this experiment. Note that the periodic boundary condition needed for this problem can be implemented as an extra linear constraint on the estimator i.e., u(·) = . . . + Pn i=1 γi (K(xi , ·) − K(xi , ·)) where xi is simply xi on the other side of the boundary. 1D wave equation. Another case where PINNs struggle to have good accuracy is the high frequency 1D wave equation described in Wang et al. (2022):   ∀x ∈ [0, 1], t ∈ [0, 1], ∂tt f − c2 ∂xx f = 0    ∀x ∈ [0, 1], f (0, x) = sin(πx) + sin(4πx)/2 (29)  ∀x ∈ [0, 1], ∂t f (0, x) = 0    ∀t ∈ [0, 1], f (t, 0) = f (t, 1) = 0. Under both noiseless and noisy conditions we compare PIKS with standard PDE solvers (Euler and Crank–Nicolson), standard PINNs, PINNs augmented with the NTK kernel (Wang et al., 2022) and PIKL (Doumèche et al., 2024). The results presented in Table 2 show that PIKS performs better than all other methods on both clean and noisy data. Note that this is with a RBF kernel which only required small-scale tuning of regularization and length-scale. For this experiment PIKS used 22 000 data points in total, PIKL and the classic solvers used ∼ 100 000 while PINNs used ∼ 24 000 000 points. Figure 2 compares the error of the different methods on a logarithmic scale. 18

10−1

ν 0.5 1.5 2.5 ∞

10−7 10−10

FEM

RMSE

10−4 RMSE

FEM PIKS

100

PIKS

10−1

10−2

10−3

10−13 102

103 104 Number of training points

10−4

105

Figure 3: FEM vs. PIKS on noiseless data of different smoothness. All FEM results apart from ν = 0.5 overlap at the top of the plot. PIKS with ν = ∞ plateaus at numerical precision.

10−2 100 Noise standard deviation

Figure 4: FEM vs. PIKS with increasing noise on the boundary conditions. 500 data points were used; the true function comes from a Matérn 3/2 kernel.

Poisson equation. The final comparison is against a FEM solver (Baratta et al., 2023) with piecewise linear elements, with a focus on how accuracy scales with the training-set size and with the amount of noise. By taking multiple Poisson PDEs with decreasing levels of smoothness (more details about the definitions available in Section F) we can additionally evaluate the performance of PIKS in a misspecified setting. We take functions from the Matérn family which are parametrized by ν. A ν-Matérn function is ⌈ν⌉ − 1 times differentiable and since PIKS is used with a Gaussian kernel, only the ν = ∞ target matches the smoothness of the Gaussian kernel, with smaller ν signifying stronger misspecification. Meanwhile the FEM solver uses linear elements and as can be seen in Fig. 3 its performance does not depend on the function smoothness. PIKS performance instead gradually degrades based on the level of misspecification, from the well-specified setting where PIKS reaches perfect accuracy with as few as 1000 points (an RMSE difference with FEM of 10 orders of magnitude) down to ν = 0.5 where kernel learning performs essentially on par with FEM. We note that FEM scales much better with dataset size, but even with ten times as much data its performance is far from kernel regression. Similar reasoning follows with noisy data (see Fig. 4): when noise is small, PIKS greatly outperforms FEM (for a moderately misspecified problem, ν = 1.5); as noise increases the two algorithms attain equally accurate estimates.

Conclusions

We formulated and analyzed PIKS, a kernel-based framework for physics-informed learning that incorporates linear differential constraints into the standard regression objective. The resulting estimator can be computed in closed form by solving a linear system. Our principal contribution is to establish learning guarantees in the misspecified setting: PIKS is universally consistent, jointly recovering u∗ and Du∗ asymptotically even when the target function u∗ lies outside the native RKHS. We verify the required assumptions in Sobolev-space settings under two sampling schemes. Under suitable source conditions, we also derive convergence rates for estimating both function values and derivatives. Small-scale experiments illustrate these results, showing that PIKS can learn even under severe misspecification and can be competitive with existing methods.

Several questions remain open. The experiment of Section 5.1 suggests that derivative observations improve the estimation of the function values themselves (n value plus n gradient points outperform 2n value points) but our rate does not capture this effect. It would be interesting to establish theoretically when the two data sources influence each other, in either direction. On the algorithmic side, introducing separate weights for the value and differential loss terms could accommodate differences in noise levels, sample sizes, and physical scales. From a computational perspective, PIKS shares the principal limitation of standard kernel methods: solving the resulting linear system scales poorly with the dataset size. Developing efficient approximations is therefore an important direction for enabling large-scale applications. Extending our method to nonlinear PDEs would be an important and challenging research direction, as well as considering inverse problems for PDEs with kernel methods.

Acknowledgements The authors thank Rayan Autones for his contributions to the early experimental stages of the project, his exploration of the relevant literature, and many helpful discussions. This material is based upon work supported by the Air Force Office of Scientific Research under award number FA8655-23-1-7083. The research was supported in part by the MIUR Excellence Department Project awarded to Dipartimento di Matematica, Università di Genova, CUP D33C23001110001.

Notation and RKHS basics

A.1

Notation

We denote by N the set of nonnegative integers including 0. For functions Rd → R, where d, m ∈ N, and for α = (α1 , . . . , αd ) ∈ Nd0 we denote by ∂ α the multi-index derivative ∂α =

∂ α1 ∂ αd . · · · ∂xα1 1 ∂xαd d

We also define the modulus |α| = α1 + · · · + αd . For a function of two variables G(x, y), defined over Rd × Rd , we denote by ∂1α G(x, y) and ∂2α G(x, y) the derivative ∂ α applied to G(x, y) as a function of x and y, respectively. For α, β ∈ Nd0 , we denote by ∂ α,β = ∂1α ∂2β

(30)

the bilateral multi-index derivative. In particular, if we see a function of two d-dimensional variables as a function of one 2d-dimensional variable, we can write ∂ α,β =

∂ αd ∂ α1 ∂ β1 ∂ βd · · · . αd · α1 · · · β ∂x1 ∂xd ∂x 1 ∂xβd d+1

For a ≤ b ∈ N, we denote by Ja, bK the set {n ∈ N, a ≤ n ≤ b}.

A.2

RKHS basics

We consider a compact set X ⊂ Rd , for some integer d ≥ 1, such that X is the closure of its nonempty interior. We consider a continuous kernel K : X × X → R. We denote by H the RKHS associated to K, which is in particular a set of functions X → R. 20

For x ∈ X , we define Kx :

X −→ R y 7−→ K(x, y).

It is well-known that the elements Kx belong to the RKHS H (in fact, the span of all the Kx for x ∈ X is dense in H). They satisfy the reproducing property: ∀f ∈ H, ∀x ∈ X ,

⟨f, Kx ⟩H = f (x).

(31)

Definition 2. Assume that X ⊂ Rd is compact and is equal to the closure of its interior. For s ∈ N, we define C s (X ) as the space of functions u : X → R such that u ∈ C s (Int(X )) and for all multi-indices α ∈ Nd0 with |α| ≤ s, the derivative ∂ α u : Int(X ) → R can be extended to a continuous function of X . Proposition 7 (Reproducing property for the derivatives). Assume that X ⊂ Rd is compact and is equal to the closure of its interior, and assume that K ∈ C 2s (X × X ). Then, any function u in H is s times continuously differentiable on X , i.e. we have a natural embedding H ,→ C s (X ). Furthermore, for any index α ∈ Nd0 such that |α| ≤ s, for any x ∈ X , the function

belongs to H and we have

Kxα : y 7→ ∂1α K(x, y)

(32)

⟨u, Kxα ⟩H = ∂ α u(x).

(33)

Proposition 7 is proved in (Zhou, 2008, Theorem 1). We now restate and prove Proposition 1 from Section 3. Lemma 8. Assume that X is compact and satisfies Int(X ) = X , and consider a linear differential operator of the form X D= cα ∂ α , (34) |α|≤s

(35)

where D1 means that we applied the operator D with respect to the first variable of (x, y) 7→ K(x, y). P α Proof. For any x ∈ X , Du(x) = |α|≤s cα (x)∂ u(x) is a linear combination of the functionals α u 7→ ∂ u(x), which are bounded according to Proposition 7, and it is thus a bounded functional. We obtain the representer by linear combination of the representers Kxα : X KxD = cα (x)Kxα . |α|≤s

PIML with kernels

In this section, we work under Assumption 1. We define the sampling operator b: A

n m H −→ R ×R u 7−→ u(xi )ni=1 , Du(zj )m j=1 ,

(36)

where Rn × Rm is equipped with the normalized Euclidean norm ∥(a, b)∥2n,m = ∥a∥2n + ∥b∥2m =

i=1

j=1

1 X 2 1X 2 ai + bj . n m

b is a random bounded operator. Indeed, it depends on the random samples The sampling operator A xi , 1 ≤ i ≤ n and zj , 1 ≤ j ≤ m, and it is bounded thanks to (6) and (7). b∗ . Let u ∈ H and (a, b) ∈ Rn × Rm . Let us compute the adjoint A b∗ (a, b), u⟩H = ⟨(a, b), Au⟩ b n,m = 1 ⟨A n

n X

ai u(xi ) +

i=1 n X

1 X bj Du(zj ) m j=1

m 1 X 1 ⟨ai Kxi , u⟩H + ⟨bj KzDj , u⟩H n m i=1 j=1 * n + m 1X 1 X = ai Kxi + bj KzDj , u . n m

i=1

This proves that b∗ (a, b) = 1 A n

n X

j=1

ai Kxi +

i=1

1 X bj KzDj . m

(37)

j=1

Let us now denote n m Y = (yi )ni=1 , (wj )m j=1 ∈ R × R . We can reformulate the physics-informed empirical risk (10) defined in Section 3 as bλ (u) = ∥Au b − Y ∥2n,m + λ∥u∥2H . R

(38)

For any λ > 0, the empirical risk (38) admits a unique minimizer. Furthermore, this minimizer can be expressed in closed form as −1 b∗ A b + λI b∗ Y. u bλ := A A (39) Indeed, the function Lλ is strongly convex on the RKHS H. By solving ∇Lλ (u) = 0, we find u bλ as the unique solution, which shows that it is the unique minimizer. This proves that the PIKS estimator (10) is well-defined. Proposition 9. The estimator u bλ can be decomposed as u bλ =

n X

αi Kxi +

i=1

m X j=1

βj KzDj ,

(40)

where (α, β) = (K + λJ)−1 Y ∈ Rn+m , with J the diagonal (n + m) × (n + m) matrix satisfying ( n if 1 ≤ i ≤ n Ji,i = m if n + 1 ≤ i ≤ n + m and K the kernel matrix, which is a (n + m) × (n + m) block matrix:    A C , K=   ⊤ C B where the blocks are the following: A ∈ Rn×n ,

Ai,i′ = ⟨Kxi , Kxi′ ⟩H ,

∀i, i′ ∈ J1, nK;

C ∈ Rn×m ,

Ci,j ′ = ⟨Kxi , KzDj ′ ⟩H ,

∀i ∈ J1, nK, ∀j ′ ∈ J1, mK.

B ∈ Rm×m ,

Bj,j ′ = ⟨KzDj , KzDj ′ ⟩H ,

∀j, j ′ ∈ J1, mK;

bA b∗ = KJ−1 . Then, using (39) and the push-through Proof. First observe with (36) and (37) that A identity, we have −1 b∗ A b + λI b∗ Y u bλ = A A −1 bA b∗ + λI b∗ A Y =A =

n X

αi Kxi +

i=1

m X

βj KzDj ,

j=1

where we see with (37) that −1 −1 bA b∗ + λI (α, β) = J−1 A Y = J−1 KJ−1 + λI Y = (K + λJ)−1 Y.

Theoretical analysis

C.1

Assumptions

In this section, we make a slight change in the sampling setting with respect to Section 2, by temporarily forgetting about the target u∗ and replacing it with two decoupled constraints. This leads to a more general setting, which we will later specialize back to the setting of this paper to obtain the main results. Consider indeed the same sampling setting as in Section 2, with a compact set X ⊂ Rd and two variables X and Z over X with respective probability distributions ρX and ρZ . We then consider two functions h ∈ L2 (ρX ) ∩ L∞ (ρX ), q ∈ L2 (ρZ ) ∩ L∞ (ρZ ), and instead of (1), we now consider Y = h(X) + ϵ, W = q(Z) + η. (41)

The datasets (xi , yi )ni=1 and (zj , wj )m j=1 are then defined as i.i.d. copies of (X, Y ) and (Z, W ) respectively, and the PIKS estimator u bλ is defined — as in Section 2 — as: bλ (u) = 1 R n

bλ (u), u bλ := arg min R u∈H

n X i=1

1 X (u(xi ) − yi ) + (Du(zj ) − wj )2 + λ∥u∥2H . m 2

(42)

j=1

As we will see later, under our working assumptions, for any u ∈ H, we have u ∈ L2 (ρX ) and Du ∈ L2 (ρZ ) which allows us to study the error q E(u) := ∥u − h∥2L2 (ρX ) + ∥Du − q∥2L2 (ρZ ) . Throughout Section C, we work under Assumptions 1, 4 and 5. In particular, Assumption 4 guarantees that there exist My , Mw > 0 such that almost surely we have |yi | ≤ My and |wj | ≤ Mw .

C.2

Error decomposition

As is standard in theory of kernel methods, we introduce in this section an intermediate function uλ ∈ H, which allows us to decompose the error between an estimation term and an approximation term. Proposition 10. The two linear operators SρX :

H −→ L2 (ρX ) u 7−→ [u]ρX ,

Dρ Z :

H −→ L2 (ρZ ) u 7−→ [Du]ρZ .

and

are well-defined and bounded. Proof. For any u ∈ H, for any x ∈ supp(ρX ), we find using Assumption 5 that |u(x)| = |⟨u, Kx ⟩H | ≤ ∥u∥H ∥Kx ∥H ≤ κ∥u∥H , and thus we have ∥u∥2L2 (ρX ) =

Z X

|u(x)|2 dρX (x) ≤ κ2 ∥u∥2H .

This shows that [u]ρX belongs to L2 (ρX ), i.e. SρX is well defined, and SρX is bounded with ∥SρX ∥op ≤ κ. In a similar way, for any u ∈ H, for any z ∈ supp(ρZ ), we find using again Assumption 5 that |Du(z)| = |⟨u, KzD ⟩H | ≤ ∥u∥H ∥KzD ∥H ≤ κD ∥u∥H . We thus have ∥DρZ u∥2L2 (ρZ ) = which shows that [Du]ρZ ∥DρZ ∥op ≤ κD .

belongs to L2 (ρ

Z X

|Du(z)|2 dρZ (z) ≤ κ2D ∥u∥2H ,

Z ), i.e.

DρZ is well defined, and DρZ is bounded with

We can then define the operator Aρ :

H −→ L2 (ρX ) × L2 (ρZ ) u 7−→ (SρX u, DρZ u),

(43)

where L2 (ρX ) × L2 (ρZ ) is equipped with the standard (Hilbertian) product norm Z Z 2 2 2 2 |f2 (z)|2 dρZ (z), |f1 (x)| dρX (x) + ∥(f1 , f2 )∥ρ = ∥f1 ∥L2 (ρX ) + ∥f2 ∥L2 (ρZ ) =

(44)

where we also denote by ⟨·, ·⟩ρ the corresponding scalar product. Since SρX and DρZ are bounded, Aρ is bounded with norm q q ∥Aρ ∥op ≤ ∥SρX ∥2op + ∥DρZ ∥2op ≤ κ2 + κ2D . (45) Proposition 11. For λ > 0, we define the population estimator as uλ = arg min ∥Aρ u − (h, q)∥2ρ + λ∥u∥2H . u∈H

Its expression is uλ = A∗ρ Aρ + λI

−1

A∗ρ (h, q).

(46)

Proof. The regularized population risk Rλ (u) = ∥Aρ u − (h, q)∥2ρ + λ∥u∥2H is strongly convex on the RKHS H. By solving ∇Rλ (u) = 0, we find uλ as the unique solution, which shows that it is the unique minimizer. We can now decompose the error as E(u) = ∥Aρ u bλ − (h, q)∥ρ

≤

∥Aρ u bλ − Aρ uλ ∥ρ + ∥Aρ uλ − (h, q)∥ρ . {z } {z } | | estimation error

(47)

approximation error

The approximation error is studied in Section C.4 while the estimation error is studied in Section C.5.

C.3

Operator properties

In this section we define and study the integral and covariance operators that are key in the theoretical proofs. Let us define the integral operator as L = Aρ A∗ρ : L2 (ρX ) × L2 (ρZ ) → L2 (ρX ) × L2 (ρZ ).

(48)

Integral operators are common objects in kernel ridge regression analysis, where they can be used to obtain convergence results. We use the integral operator (48) for the same purposes. The following proposition establishes some of its basic properties. Proposition 12. For ρX -almost all x and ρZ -almost all z, we have Z Z L(f, g)(x, z) = ⟨Ky , Kx ⟩H f (y)dρX (y) + ⟨KtD , Kx ⟩H g(t)dρZ (t) , X Z Z X D D D ⟨Ky , Kz ⟩H f (y)dρX (y) + ⟨Kt , Kz ⟩H g(t)dρZ (t) . X

Furthermore, L is positive and trace class, and it admits a decomposition in an orthonormal system of eigenvectors (fi , gi )i∈I . X L(f, g) = µi ⟨(f, g), (fi , gi )⟩ρ (fi , gi ), (49) i∈I

where for all i ∈ I, µi > 0. Before proving Proposition 12, let us first prove the following Lemma which gives the expression of A∗ρ . Lemma 13. For any (f, g) ∈ L2 (ρX ) × L2 (ρZ ), we have Z Z ∗ Aρ (f, g) = f (y)Ky dρX (y) + g(t)KtD dρZ (t). X

(50)

Proof. Let (f, g) ∈ L2 (ρX ) × L2 (ρZ ). We have ⟨A∗ρ (f, g) , φ⟩H = ⟨(f, g) , Aρ φ⟩ρ Z Z = f (y)φ(y)dρX (y) + g(t)Dφ(t)dρZ (t) X X Z Z = f (y)⟨Ky , φ⟩H dρX (y) + g(t)⟨KtD , φ⟩H dρZ (t) X XZ Z D , = f (y)Ky dρX (y) + g(t)Kt dρZ (t) , φ X

where we used the reproducing properties (6) and (7), and the last equality is justified by the fact that y 7→ f (y)Ky is Bochner integrable with respect to ρX and t 7→ g(t)KtD is Bochner integrable with respect to ρZ . This proves formula (50). Proof of Proposition 12. Let us first establish the expression of L. Composing (50) with the expression (43) of Aρ , we obtain that for ρX -almost all x and ρZ -almost all z, we have Z Z L(f, g)(x, z) = f (y)Ky (x)dρX (y) + g(t)KtD (x)dρZ (t) , X Z X Z D D f (y)Ky dρX (y) (z) + D g(t)Kt dρZ (t) (z) X ZX Z = ⟨Ky , Kx ⟩H f (y)dρX (y) + ⟨KtD , Kx ⟩H g(t)dρZ (t) , X X Z Z D D D ⟨Ky , Kz ⟩H f (y)dρX (y) + ⟨Kt , Kz ⟩H g(t)dρZ (t) . X

The second equality is the desired expression. By definition, L = Aρ A∗ρ is self-adjoint. Let us show that Aρ is Hilbert-Schmidt. Let (ei )i∈N be an orthonormal basis of H, and let us prove that X ∥Aρ ∥2HS := ∥Aρ ei ∥2ρ < +∞. i∈N

Let us first observe that X i∈N

∥Aρ ei ∥2ρ = =

X Z X

i∈N

Z X X i∈N

Z X X i∈N 2

≤κ

|ei (x)| dρX (x) + 2

|ei (x)| dρX (x) +

Z X

|Dei (z)| dρZ (z)

Z X X i∈N

|⟨ei , Kx ⟩H | dρX (x) +

|Dei (z)|2 dρZ (z)

Z X X i∈N

|⟨ei , KzD ⟩H |2 dρZ (z)

+ κ2D ,

where the second equality is obtained by monotone convergence, and the final inequality comes from Assumption 5. This shows that Aρ is Hilbert-Schmidt, and as a consequence, L is trace class. To conclude the proof, trace class operators of separable Hilbert spaces (such as L2 (ρX )×L2 (ρZ )) are compact, and since L is also self-adjoint and positive, the spectral theorem guarantees that there exists an orthonormal system (fi , gi )i∈I , there exist µi > 0, i ∈ I such that for any (f, g) ∈ L2 (ρX ) × L2 (ρZ ), X L(f, g) = µi ⟨(f, g), (fi , gi )⟩ρ (fi , gi ). i∈I

Lemma 14. For any x ∈ X , we have ∥Kx ⊗ Kx ∥HS ≤ κ2 and ∥KxD ⊗ KxD ∥HS ≤ κ2D . Proof. Recall that we work under Assumption 5. Since H is a separable Hilbert space, we can consider a Hilbert basis (ei )i∈N of H. We have by definition Kx ⊗ Kx (ei ) = ⟨Kx , ei ⟩H Kx , so X i∈N

∥Kx ⊗ Kx (ei )∥2H = ∥Kx ∥2H

X ⟨Kx , ei ⟩2H = ∥Kx ∥4H ≤ κ4 . i∈N

This proves that Kx ⊗Kx is a Hilbert-Schmidt operator, with Hilbert-Schmidt norm ∥Kx ⊗Kx ∥HS ≤ κ2 . The proof of the second point is identical, replacing Kx by KxD and κ by κD . Let us now define the covariance operators b =A b∗ A, b Σ Lemma 15. We have b= 1 Σ n

n X

Σ = A∗ρ Aρ .

(51)

1 X D Kzj ⊗ KzDj , m

(52)

Σ = Ex∼ρX [Kx ⊗ Kx ] + Ez∼ρZ [KzD ⊗ KzD ],

(53)

i=1

Kxi ⊗ Kxi +

j=1

and where the expected values are defined as Bochner integrals in the space of Hilbert-Schmidt operators.

Proof. Combining (36) and (37), we get that for all u ∈ H, we have b = 1 Σu n

n X

u(xi )Kxi +

i=1

m n m 1 X 1X 1 X D Du(zj )KzDj = ⟨Kxi , u⟩H Kxi + ⟨Kzj , u⟩H KzDj m n m j=1 i=1 j=1   n m X X 1 1 = K xi ⊗ K xi + KzDj ⊗ KzDj  u, n m i=1

j=1

which proves (52). Combining (43) and (50), we get that for all u ∈ H, we have Z Z Z Z Σu = u(y)Ky dρX (y) + Du(t)KtD dρZ (t) = (Ky ⊗ Ky u)dρX (y) + (KtD ⊗ KtD u)dρZ (t) X X X X = Ex∼ρX [Kx ⊗ Kx ] + Ez∼ρZ [KzD ⊗ KzD ] u. This proves (53), provided we justify the last equality. Observe indeed that as proved in Proposition 14, for any x, z ∈ X , Kx ⊗ Kx is a Hilbert-Schmidt operator of norm ∥Kx ⊗ Kx ∥HS ≤ κ2 and KzD ⊗ KzD is a Hilbert-Schmidt operator of norm ∥KzD ⊗ KzD ∥HS ≤ κ2D . Since we have Z ∥Ky ⊗ Ky ∥HS dρX (y) ≤ κ2 X

and

Z X

∥KtD ⊗ KtD ∥HS dρZ (t) ≤ κ2D ,

the Bochner integrals Ex∼ρX [Kx ⊗ Kx ] and Ez∼ρZ [KzD ⊗ KzD ] are well defined. Finally, we use the fact that the Bochner integral commutes with bounded operators to factorize the expression.

C.4

Approximation error

In this section, we focus on the approximation part of the decomposition (47). In Section C.1, we defined (h, q) ∈ L2 (ρX ) × L2 (ρZ ), let us now define (hΠ , q Π ) ∈ L2 (ρX ) × L2 (ρZ ) as the orthogonal projection of (h, q) on ran Aρ , which is a closed subspace of L2 (ρX ) × L2 (ρZ ). For the error bounds, we introduce a slightly modified version of Assumption 6, which applies to the projection (hΠ , q Π ) of (h, q) onto ran Aρ . Assumption 7 (Source condition, variant). There exists r ∈ (0, 1] such that (hΠ , q Π ) ∈ ran Lr , i.e. there exists (fe, ge) ∈ L2 (ρX ) × L2 (ρZ ) such that (hΠ , q Π ) = Lr (fe, ge). Proposition 16. We have ∥Aρ uλ − (hΠ , q Π )∥ρ −→ 0. λ→0

Furthermore, under Assumption 7, we have the rate ∥Aρ uλ − (hΠ , q Π )∥ρ ≤ λr ∥(fe, ge)∥ρ . Proof. Let us denote Lλ = Aρ A∗ρ + λI. Recalling the definition (46) of uλ , we have uλ = (A∗ρ Aρ + λI)−1 A∗ρ (h, q) = (A∗ρ Aρ + λI)−1 A∗ρ (hΠ , q Π ) Π Π = A∗ρ L−1 λ (h , q ),

(54)

where the second equality holds because (hΠ , q Π ) is the orthogonal projection of (h, q) onto ran Aρ and the third equality is the so-called push-through identity. We thus have Π Π Π Π Aρ uλ − (hΠ , q Π ) = Aρ A∗ρ L−1 λ (h , q ) − (h , q ) −1 Π Π = (LL−1 λ − Lλ Lλ )(h , q ) Π Π = −λL−1 λ (h , q ).

(55)

Recall from (49) that there exists an orthonormal system of eigenvectors (fi , gi ) such that X L(f, g) = µi ⟨(f, g), (fi , gi )⟩ρ (fi , gi ), i∈I

with µi > 0. We thus see that (fi , gi )i∈I is a Hilbert basis of ran Aρ = (Ker L)⊥ . Since by definition (hΠ , q Π ) ∈ ran Aρ , we can decompose (hΠ , q Π ) in the basis (fi , gi )i∈I , hence we can write Π

∥Aρ uλ − (h , q

Π Π 2 )∥2ρ = ∥λL−1 λ (h , q )∥ρ =

X i∈I

λ λ + µi

⟨(hΠ , q Π ), (fi , gi )⟩2ρ .

For all i ∈ I, since µi > 0, we have

λ λ + µi

⟨(hΠ , q Π ), (fi , gi )⟩2ρ −→ 0. λ→0

By dominated convergence we thus have ∥Aρ uλ − (hΠ , q Π )∥2ρ −→ 0. λ→0

For the second point, using again the expression (55), we have

Π Π Aρ uλ − (hΠ , q Π ) = −λL−1 λ (h , q ) = −λL−1 Lr (fe, ge), λ

and thus Π

∥Aρ uλ − (h , q

)∥2ρ =

X i∈I

λ λ + µi

e e), (fi , gi )⟩2 µ2r i ⟨(f , g ρ

X λr λ1−r µr 2

i ⟨(fe, ge), (fi , gi )⟩2ρ λ + µi i∈I X λ1−r µr 2 2r i =λ ⟨(fe, ge), (fi , gi )⟩2ρ λ + µi i∈I X 2r ≤λ ⟨(fe, ge), (fi , gi )⟩2ρ

i∈I

≤ λ ∥(fe, ge)∥2ρ , where we observed that for any r ∈ [0, 1], we have either λ1−r µri ≤ λ or λ1−r µri ≤ µi , so λ1−r µri ≤ λ + µi .

C.5

Estimation error

In this section, we want to bound the estimation part of the decomposition (47), that is, ∥Aρ u bλ − Aρ uλ ∥ρ . To make proofs more compact, we will introduce the notation   if Assumption 7 does not hold 0 r̃ = r (56) if Assumption 7 holds with r ∈ (0, 1/2]   1/2 if Assumption 7 holds with r ∈ (1/2, 1]. bλ = A b∗ A+λI b When r̃ = 0, we take as a convention (fe, ge) := (hΠ , q Π ). Let us introduce the notation Σ ∗ and Σλ = Aρ Aρ +λI, which will be used at several moments and will help keeping the computations compact. We begin by proving two technical lemmas. b X = 1 Pn Kx ⊗ Kx , and define Lemma 17. Let us define the operators ΣX = E[Kx ⊗ Kx ] and Σ n

i=1

1 −1 b X )Σ− 2 . BX = Σλ 2 (ΣX − Σ λ

(57)

1 Pm

D D j=1 Kzj ⊗ Kzj , and define

bZ = Similarly, let ΣZ = E[KzD ⊗ KzD ] and Σ m

− b Z )Σ− 2 . BZ = Σλ 2 (ΣZ − Σ λ

Then, if ∥BX ∥op ≤ 41 and ∥BZ ∥op ≤ 41 , we have 1

−1

b 2 ∥op ≤ ∥Σλ2 Σ λ

√

(58)

−1

2 b Proof. Let us define the operator B = Σλ 2 (Σ − Σ)Σ λ . We observe that 1

−1

−1 b 2 ∥2 = ∥Σ 2 Σ b −1 2 ∥Σλ2 Σ op λ λ λ Σλ ∥op = ∥(I − B) ∥op .

If ∥B∥op ≤ 1/2, the Neumann series allows us to bound ∥(I − B)−1 ∥op ≤ (1 − ∥B∥op )−1 ≤ 2, which implies the statement of the lemma. We can write: −1

−1

b X )Σ 2 + Σ 2 (ΣZ − Σ b Z )Σ 2 B = Σλ 2 (ΣX − Σ λ λ λ = BX + BZ .

Thus: ∥B∥op ≤ ∥BX ∥op + ∥BZ ∥op ≤

1 1 1 + = . 4 4 2

Lemma 18. If we have 112κ2 λ≥ log 3n then, with probability at least 1 − δ, we have

n+4 δ

(59)

1 ∥BX ∥op ≤ . 4

Analogously, if we have 112κ2D λ≥ log 3m then, with probability at least 1 − δ, we have

m+4 δ

1 ∥BZ ∥op ≤ . 4 30

(60)

Proof. Let us first bound BX . Let us denote ΣX,λ := ΣX + λI. First note that Σλ = Σ + λI = ΣX + ΣZ + λI ⪰ ΣX + λI, so −1

−1

2 2 2 2 ∥op ∥ΣX,λ ΣX ΣX,λ ∥op ∥ΣX,λ Σλ 2 ∥op ∥Σλ 2 ΣX Σλ 2 ∥op ≤ ∥Σλ 2 ΣX,λ

−1

2 2 ΣX ΣX,λ ∥op . ≤ ∥ΣX,λ

The idea is now to apply Tropp’s concentration inequality (Proposition 49) to the random −1

−1

2 2 2 variables Zi := Ui ⊗ Ui with Ui := ΣX,λ Kxi . We have E[Zi ] = ΣX,λ ΣX ΣX,λ := T . We have

∥Kxi ∥2H κ2 ≤ , λ λ

∥Ui ⊗ Ui ∥op = ∥Ui ∥2H ≤ 2

so we can choose R = κλ . Now observe that E[(Ui ⊗ Ui − T )2 ] = E[∥Ui ∥2H Ui ⊗ Ui ] − T 2 ⪯ E[∥Ui ∥2H Ui ⊗ Ui ] ⪯ RT, 2

Let us define S := RT . We then have σ 2 := ∥S∥op ≤ κλ ∥T ∥op ≤ κλ (since T ⪯ I). We define (∥ΣX ∥op +λ)∥T ∥1 ∥S∥1 α = ∥S∥ = . We can thus use Proposition 49 to conclude that, with probability ∥ΣX ∥op op 1 − δ, r n βκ2 3βκ2 1X + , Zi − T ≤ ∥BX ∥op = n λn λn i=1

4(∥ΣX ∥op +λ)∥T ∥1 . By taking λ ≥ ∥ΣX ∥op δ check that we obtain the desired ∥BX ∥op ≤ 41 .

with β = 23 log

56βκ2 n

in the right member of the inequality, one can

Such a condition on λ is not explicit since β itself

depends on λ, but one can check that taking λ≥

112κ2 log 3n

n+4 δ

(61)

is enough to satisfy it, i.e. to have, with probability 1 − δ, 1 ∥BX ∥op ≤ . 4 We bound BZ by an identical argument for the second term. Utilizing the bound ∥KzD ∥H ≤ κD , we obtain that if 112κ2D m+4 λ≥ log (62) 3m δ then ∥BZ ∥op ≤ 1/4 with probability 1 − δ. Proposition 19. There exist four sequences of random variables (Un ), (Vn ), (Wm ), (Zm ) ∈ RN + such that (i) for all n, m ∈ N, we have ∥Aρ u bλ − Aρ uλ ∥ρ

≤

1 2

1 b − 2 ∥op ∥Σλ Σ λ

1 1

λ2

Un +

1 λ1−r̃

Vn +

1 1

λ2

Wm +

1 λ1−r̃

Zm ;

(ii) for all n ∈ N, for all τ > 0, with probability at least 1 − 2e−τ , we have √ τ Un ≤ CU √ , n

(63)

where CU = 4My κ; (iii) for all n ∈ N, for all τ > 0, with probability at least 1 − 2e−τ , we have √ τ Vn ≤ CV √ , n

(64)

where CV = 4κ2 (κ2 + κ2D )r−r̃ ∥(fe, ge)∥ρ ; (iv) for all m ∈ N, for all τ > 0, with probability at least 1 − 2e−τ , we have √ τ Wm ≤ CW √ , m

(65)

where CW = 4Mw κD ; (v) for all m ∈ N, for all τ > 0, with probability at least 1 − 2e−τ , we have √ τ Zm ≤ CZ √ , m

(66)

where CZ = 4κ2D (κ2 + κ2D )r−r̃ ∥(fe, ge)∥ρ . We prove Proposition 19 in this section. Let us denote b∗ Y ξb = A

ξ = A∗ρ (h, q).

(67)

b =A b∗ A b and Σ = A∗ Aρ defined in (51) and the operators Σ bλ = A b∗ A b + λI Recall the operators Σ ρ ∗ and Σλ = Aρ Aρ + λI introduced at the beginning of Section C.5. In order to prove Proposition 19, we first state and prove the following decomposition lemma. Lemma 20. We have

i h b −1 ξ . b −1 (ξb − ξ) + (Σ − Σ)Σ u bλ − uλ = Σ λ λ

(68)

b −1 ξb and uλ = Σ−1 ξ. We thus have Proof. We can rewrite the expressions (39) and (46) as u bλ = Σ λ λ b −1 ξb − Σ−1 ξ u bλ − uλ = Σ λ λ −1 b b b b −1 ξ − Σ−1 ξ = Σ ξ − Σ−1 ξ + Σ

λ λ λ λ −1 −1 −1 b (ξb − ξ) + (Σ b − Σ )ξ =Σ λ λ λ −1 b −1 b b b λ )Σ−1 ξ = Σλ (ξ − ξ) + Σλ (Σλ − Σ λ −1 b −1 −1 b b b = Σλ (ξ − ξ) + Σλ (Σ − Σ)Σλ ξ

h i b −1 (ξb − ξ) + (Σ − Σ)Σ b −1 ξ . =Σ λ λ

Lemma 21. We have ∥Σ∥r−r̃ 1 op b op ∥(fe, ge)∥ρ √ ∥ξb − ξ∥H + 1−r̃ ∥Σ − Σ∥ λ λ

1 2

1 b − 2 ∥op ∥Σλ Σ λ

∥Aρ (b uλ − uλ )∥ρ ≤

! .

(69)

Proof. First observe that 1

∥Aρ (b uλ − uλ )∥ρ = ∥Σ 2 (b uλ − uλ )∥H . Now, observe that ξ = A∗ρ (h, q) = A∗ρ (hΠ , q Π ) = A∗ρ Lr (fe, ge) = Σr A∗ρ (fe, ge). Thus, h i 1 1 b −1 (ξb − ξ) + (Σ − Σ)Σ b −1 ξ Σ 2 (b u λ − uλ ) = Σ 2 Σ λ λ i 1 1 h 1 b− 2 Σ b − 2 (ξb − ξ) + (Σ − Σ)Σ b −1 Σr A∗ρ (fe, ge) , = Σ2 Σ λ

thus 1

1 1 b − 2 ∥op b − 2 ∥op ∥Σ uλ − uλ )∥H ≤ ∥Σ 2 Σ ∥Σ 2 (b λ λ − 21 +r̃ − 21 −r̃ r̃ ∗ b e b · ∥ξ − ξ∥H + ∥Σ − Σ∥op ∥Σλ ∥op ∥Σλ Σ Aρ ∥ρ→H ∥(f , ge)∥ρ " # ∥Σ∥r−r̃ 1 1 − 21 − 12 −r̃ r̃ ∗ op b ∥op ∥ξb − ξ∥H + b op ∥Σ ≤ √ ∥Σ 2 Σ ∥Σ − Σ∥ Σ Aρ ∥ρ→H ∥(fe, ge)∥ρ , 1 λ λ λ λ 2 −r̃

where ∥ · ∥op denotes the standard operator norm on (H, ∥ · ∥H ), and ∥ · ∥ρ→H denotes the operator norm from (L2 (ρX ) × L2 (ρZ ), ∥ · ∥ρ ) to (H, ∥ · ∥H ). We used the fact that − 21 + r̃ ≤ 0 to bound − 1 +r̃

∥Σλ 2

∥op ≤ λ− 2 +r̃ .

1 b − 2 ∥op and ∥Σ− 2 −r̃ Σr̃ A∗ρ ∥ρ→H . First observe that In order to conclude, we need to bound ∥Σ 2 Σ λ λ s

− 1 −r

∥Σλ 2

− 1 −r̃

Σr̃ A∗ρ ∥ρ→H =

(Σλ 2

− 1 −r̃

Σr̃ A∗ρ )(Σλ 2

Σr̃ A∗ρ )∗ op

s =

− 1 −r̃

Σλ 2

− 1 −r̃

Σ1+2r̃ Σλ 2

. op

− 1 −r̃

µ1+2r̃

i The eigenvalues of the operator Σλ 2 Σ1+2r̃ Σλ 2 = Σ−1−2r̃ Σ1+2r̃ are (µi +λ) 1+2r̃ ≤ 1, where µi λ ∗ ∗ are the same eigenvalues from (49) (we use that Σ = Aρ Aρ and L = Aρ Aρ have the same non-zero eigenvalues.)). As a consequence, its operator norm is bounded by 1, and thus we have

− 1 −r̃

∥Σλ 2 1

Σr̃ A∗ρ ∥ρ→H ≤ 1.

−1

b 2 ∥op . We have Let us now bound ∥Σ 2 Σ λ 1

1 b − 2 ∥op ≤ ∥Σ 12 Σ− 2 ∥op ∥Σ 2 Σ b− 2 ∥Σ 2 Σ λ λ λ λ ∥op 1

−1

b 2 ∥op , ≤ ∥Σλ2 Σ λ 1

−1

√

i where we used ∥Σ 2 Σλ 2 ∥op ≤ 1 since the eigenvalues of Σ 2 Σλ 2 are √µ +λ . i

The following lemma allows us to express ξb as an empirical mean and ξ as an expectation. Lemma 22. We have

n X

1 ξb = n

yi Kxi +

i=1

1 X wj KzDj , m

(70)

j=1

and ξ = Ex∼ρX [h(x)Kx ] + Ez∼ρZ [q(z)KzD ].

(71)

Furthermore, we have the bound ∥ξ∥H ≤ ∥h∥L∞ (ρX ) κ + ∥q∥L∞ (ρZ ) κD .

(72)

Proof. Let u ∈ H. We have b u⟩H = ⟨A b∗ Y, u⟩H ⟨ξ, b n,m = ⟨Y, Au⟩ n

i=1 n X

j=1

1X 1 X = yi u(xi ) + wj Du(zj ) n m m 1 X 1 yi ⟨Kxi , u⟩H + wj ⟨KzDj , u⟩H n m i=1 j=1 * n + m X X 1 1 = yi Kxi + wj KzDj , u . n m

i=1

j=1

Since this is true for any u ∈ H, it proves (70). We prove the second formula the same way. Let u ∈ H, we have ⟨ξ, u⟩H = ⟨A∗ρ (h, q), u⟩H = ⟨(h, q), Aρ u⟩ρ

= Ex∼ρX [h(x)u(x)] + Ez∼ρZ [q(z)Du(z)] = Ex∼ρX [h(x)⟨Kx , u⟩H ] + Ez∼ρZ [q(z)⟨KzD , u⟩H ] = Ex∼ρX [h(x)Kx ] + Ez∼ρZ [q(z)KzD ] , u H . To prove (71), observe first that h(x) is ρX -almost surely bounded (as assumed at the beginning of Section C), and by Assumption 5, for all x ∈ X we have ∥Kx ∥H ≤ κ, so ∥h(x)Kx ∥H is bounded ρX -almost surely. As a consequence, the expected value Ex∼ρX [h(x)Kx ] is well defined as a Bochner integral with values in H. Analogously, ∥q(z)KzD ∥H is bounded ρZ -almost surely which allows us to define the expected value Ez∼ρZ [q(z)KzD ]. Finally, the Bochner integral commutes with bounded operators, so we can exchange the expectations and the scalar product ⟨·, u⟩H . Since this equality is true for any u ∈ H, it proves (71). Let us now prove the bound. We have ∥ Ex∼ρX [h(x)Kx ] ∥H ≤ Ex∼ρX [∥h(x)Kx ∥H ] ≤ ∥h∥L∞ (ρX ) κ. Analogously, we have ∥ Ez∼ρZ [q(z)KzD ] ∥H ≤ Ez∼ρZ [∥q(z)KzD ∥H ] ≤ ∥q∥L∞ (ρZ ) κD . Using the triangle inequality on (71) and the previous two inequalities yields (72). 34

Proof of Proposition 19. We know from Proposition 21 that ∥Aρ (b uλ − uλ )∥ρ ≤

1 2

1 b − 2 ∥op ∥Σλ Σ λ

∥Σ∥r−r̃ 1 op b op ∥(fe, ge)∥ρ √ ∥ξb − ξ∥H + 1−r̃ ∥Σ − Σ∥ λ λ

! .

(73)

b op . Let us prove (i) first by decomposing both ∥ξb − ξ∥H and ∥Σ − Σ∥ We begin with ∥ξb − ξ∥H . Using Proposition 22, we have 1 ξb = n

n X i=1

1 X yi Kxi + wj KzDj , m j=1

and ξ = Ex∼ρX [h(x)Kx ] + Ez∼ρZ [q(z)KzD ]. Let us define

1X Un = yi Kxi − Ex∼ρX [h(x)Kx ] n i=1

and m

1 X wj KzDj − Ez∼ρZ [q(z)KzD ] m

Wm =

j=1

. H

By the triangle inequality, we have ∥ξb − ξ∥H ≤ Un + Wm .

(74)

b op . Recall from Proposition 15 that we have Let us now proceed with ∥Σ − Σ∥ b= 1 Σ n

n X i=1

1 X D bX + Σ bZ, K xi ⊗ K xi + Kzj ⊗ KzDj = Σ m j=1

and Σ = Ex∼ρX [Kx ⊗ Kx ] + Ez∼ρZ [KzD ⊗ KzD ] = ΣX + ΣZ . We define b X − ΣX Vn = ∥(fe, ge)∥ρ Σ

b Z − ΣZ Zm = ∥(fe, ge)∥ρ Σ

and .

By the triangle inequality, we have b −Σ ∥(fe, ge)∥ρ Σ

≤ V n + Zm .

(75)

The inequalities (74) and (75) alongside with the decomposition (73) prove (i). The four remaining points are all proved the same way, by using Hoeffding’s inequality in Hilbert spaces. Let us prove (ii) and (iv) jointly, i.e. let us bound Un and Wm . We define the i.i.d. random variables ei = yi Kxi ,

and the i.i.d. random variables tj = wj KzDj . For i ∈ J1, nK, we have

E[ei ] = E [E[yi Kxi |xi ]] = E [h(xi )Kxi ] ,

and for j ∈ J1, mK, we have h i h i E[tj ] = E E[wj KzDj |zj ] = E q(zj )KzDj , so we can rewrite 1 Un = n and

n X i=1

! ei − E[ei ]

  m X 1  Wm = tj − E[tj ] m j=1

. H

We want to apply Hoeffding’s inequality in separable Hilbert spaces to bound Un and Wm . For completeness, we reproduced the inequality in this appendix as Proposition 48. For that, we need first to establish that the variables ei and tj are bounded. It is indeed the case since we have ∥ei ∥H ≤ |yi |∥Kxi ∥H ≤ My κ, where both My and κ are positive constants defined in Section C.1, and where the second inequality holds almost surely. As a consequence, the zero-mean variable ei − E[ei ] satisfies almost surely ∥ei − E[ei ]∥ ≤ 2My κ.

Similarly, we have ∥tj ∥H ≤ |wj |∥KzDj ∥H ≤ Mw κD , where again Mw and κD are defined in Section C.1 and the second inequality holds almost surely. As a consequence, the zero-mean variable tj − E[tj ] satisfies almost surely ∥tj − E[tj ]∥H ≤ 2Mw κD . We can thus apply the Hoeffding inequality in separable Hilbert spaces (see Proposition 48) and get that, with probability at least 1 − 2e−τ , we have r τ Un ≤ 4My κ , (76) n which proves (ii). Similarly, we find that with probability at least 1 − 2e−τ , we have r τ Wm ≤ 4Mw κD , m

(77)

which proves (iv). Let us now prove (iii) and (v). Bounding the operator norm by the Hilbert-Schmidt norm, we have n 1X (78) Vn ≤ ∥(fe, ge)∥ρ Kxi ⊗ Kxi − Ex [Kx ⊗ Kx ] n HS i=1

and

1 X D Zm ≤ ∥(fe, ge)∥ρ Kzj ⊗ KzDj − Ez [KzD ⊗ KzD ] . m HS j=1

(79)

To apply again Hoeffding inequality, we need the variables Kxi ⊗ Kxi and KzDj ⊗ KzDj to be bounded, which is the case as proved in Proposition 14. We can thus apply Hoeffding inequality in the separable Hilbert space of Hilbert-Schmidt operators (see Proposition 48) to find that with probability at least 1 − 2e−τ , we have r n 1X τ 2 Kxi ⊗ Kxi − Ex [Kx ⊗ Kx ] ≤ 4κ . (80) n n i=1

Together, (78), (80) and (72) prove (iii). Similarly, with probability at least 1 − 2e−τ , we have m

1 X D Kzj ⊗ KzDj − Ez [KzD ⊗ KzD ] m j=1

≤ 4κ2D

τ . m

(81)

Together, (79), (81) and (72) prove (v). Proposition 19 allows us to prove the following two results. 1 Corollary 23. Let δ ∈ (0, 12 ). There exists a constant C > 0 such that if λ ≥ C max n1 log nδ , m log m δ , and if λ < λ0 for an arbitrary upper bound λ0 > 0, then with probability at least 1 − δ, we have p p ln(1/δ) ln(1/δ) √ + 1−r̃ √ , ∥Aρ u bλ − Aρ uλ ∥ρ ≲ (82) 1−r̃ λ n λ m where the hidden constants in (82) depend on My , Mw , κ, κD and λ0 . Proof. We know from Proposition 19 that we can write 1 1 1 1 1 −1 2 b 2 Vn + 1 Wm + 1−r̃ Zm . ∥Aρ u bλ − Aρ uλ ∥ρ ≤ ∥Σλ Σλ ∥op 1 Un + λ1−r̃ λ λ2 λ2 1

b − 2 ∥op . We know that if we have Let us first bound ∥Σλ2 Σ λ 2 112 κ n + 4 κ2D m+4 λ≥ max log , log 3 n m δe δe e we have simultaneously then with probability 1 − 2δ, 1 ∥BX ∥op ≤ , 4

1 ∥BZ ∥op ≤ , 4

so Proposition 17 implies 1

−1

b 2 ∥op ≤ ∥Σλ2 Σ λ

√

Now using a union bound on this probability bound alongside with the bounds 2 from Proposition 19, 2 κ 112 κ n+4 m+4 D e we get that if λ ≥ where we pick τ = ln(1/δ) , then with probability 3 max n log δe , m log δe e we have at least 1 − 10δ, q q q  q  e e e e ln(1/δ) ln(1/δ) ln(1/δ) ln(1/δ) √ + CV 1−r̃ √ + CW + CZ 1−r̃ √  ∥Aρ u bλ − Aρ uλ ∥ρ ≤ 2  CU 1√ 1√ λ n λ m λ2 n λ2 m ! p p ln(1/δ) ln(1/δ) ′ √ + 1−r̃ √ ≤ C , 1−r̃ λ n λ m 37

e and the second inequality is true provided λ < λ0 for an arbitrary choice where we defined δ := 10δ, of λ0 > 0, and where the constant C ′ depends on My , Mw , κ, κD , λ0 and ∥(fe, ge)∥ρ . Using δ, the condition on λ can be written as 2 112 κ 10n + 40 κ2D 10m + 40 λ≥ , (83) max log , log 3 n δ m δ which we can simplify (using that δ < 1/2) as 1 n 1 m λ ≥ C max . log , log n δ m δ

Recall that in the current section, we work under Assumptions 1, 4 and 5. Consider furthermore the following assumption. Assumption 8. The pair (h, q) is in the closure of the range of Aρ taken in L2 (ρX ) × L2 (ρZ ): (h, q) ∈ ran Aρ . Remark 2. Since for any r > 0, ran Lr = ran Aρ, we see that Assumption 6 is equivalent to Assumptions 7 and 8 together. 1 log m Corollary 24. Let δ ∈ (0, 12 ). Under Assumption 7, if λ ≥ C max n1 log nδ , m δ , and if λ < λ0 for an arbitrary upper bound λ0 > 0, then with probability at least 1 − δ, we have p p ln(1/δ) ln(1/δ) Π Π √ + 1−r̃ √ + λr , (84) ∥Aρ u bλ − (h , q )∥ρ ≲ 1−r̃ λ n λ m where the hidden constants in (84) depend on My , Mw , κ, κD and λ0 . Furthermore, under Assumption 8, this simplifies to p p ln(1/δ) ln(1/δ) √ + 1−r̃ √ + λr . ∥Aρ u bλ − (h, q)∥ρ ≲ (85) 1−r̃ λ n λ m In particular, under Assumption 6, both Assumptions 7 and 8 hold, so (85) holds. Proof. Using together the approximation error from Equation (54) and Proposition 23 yields the first result. Then, Assumption 8 is equivalent to (hΠ , q Π ) = (h, q), which gives the second result. The final remark comes from the fact that Assumption 6 is equivalent to Assumptions 7 and 8 together, as mentioned in Remark 2. 1

Corollary 25. Consider Assumption 6, and let us pick λ = N −1/2 if r ≤ 1/2, and λ =N − 2r+1 1 log m if r > 1/2, where N := min(n, m). If δ ∈ (0, 12 ) is such that λ ≥ C max n1 log nδ , m δ , and if λ < λ0 for an arbitrary upper bound λ0 > 0, then with probability at least 1 − δ, we have p ∥Aρ u bλ − (h, q)∥ρ ≲ ln(1/δ)N −r/2 if r ≤ 1/2 (86) p r ∥Aρ u bλ − (h, q)∥ρ ≲ ln(1/δ)N − 2r+1 if r > 1/2, (87) where the hidden constants in (86) and (87) depend on My , Mw , κ, κD and λ0 .

Remark. Since λ ≥ N −1/2 and the condition on λ is of the form λ ≳ log(N/δ) , we see that for any N δ ∈ (0, 1), the condition on λ is always satisfied for N large enough. Proof. As mentioned in Remark 2, Assumption 6 is equivalent to Assumptions 7 and 8 together. We can thus use the rates in (85). Suppose first that r ≤ 1/2, i.e. r̃ = r, and let us set λ = N −1/2 . With probability 1 − δ, p p ln(1/δ) ln(1/δ) √ + 1−r √ + λr ∥Aρ u bλ − (h, q)∥ρ ≲ 1−r λ n λ m p p 1−r 1−r −1/2 ≲ + ln(1/δ)N 2 m−1/2 + N −r/2 ln(1/δ)N 2 n p ≲ ln(1/δ)N −r/2 + N −r/2 . 1

Suppose now that r > 1/2, i.e. r̃ = 1/2, and let us set λ = N − 2r+1 . With probability 1 − δ, p p ln(1/δ) ln(1/δ) √ √ + √ √ + λr ∥Aρ u bλ − (h, q)∥ρ ≲ λ n λ m p p 1 1 r ≲ ln(1/δ)N 4r+2 n−1/2 + ln(1/δ)N 4r+2 m−1/2 + N − 2r+1 p r r ln(1/δ)N − 2r+1 + N − 2r+1 . ≲

C.6

Asymptotic convergence

We take again the notation (hΠ , q Π ) from Section C.4 to denote the orthogonal projection (in L2 (ρX ) × L2 (ρZ )) of (h, q) on the closed subspace ran Aρ . Recall that for any fixed λ > 0, u bλ already depends on (n, m) (we dropped the index for lighter notation), which means that for a choice λn,m that depends on (n, m), the sequence u bλn,m depends on (n, m) in two different ways. Proposition 26. Let N = min(n, m). Let (λn,m ) be any regularization sequence such that λn,m → 0,

log N →0 λ3n,m N

as n, m → ∞.

Then, almost surely, ∥Aρ u bλn,m − (hΠ , q Π )∥ρ

−→

n,m→+∞

Equivalently, almost surely, for every ϵ > 0, there exists Nϵ ∈ N such that ∀(n, m) ∈ N2 , n, m ≥ Nϵ =⇒ ∥Aρ u bλn,m − (hΠ , q Π )∥ρ < ϵ . Proof. Using the triangle inequality and the decomposition of Proposition 19 with λ = λn,m , we have ∥Aρ u bλn,m − (hΠ , q Π )∥ρ ≤ ∥Aρ u bλn,m − Aρ uλn,m ∥ρ + ∥Aρ uλn,m − (hΠ , q Π )∥ρ . Proposition 16 shows that since λn,m → 0, the approximation error goes to 0: ∥Aρ uλn,m − (hΠ , q Π )∥ρ 39

−→

n,m→+∞

We thus only need to prove that, almost surely, ∥Aρ u bλn,m − Aρ uλn,m ∥ρ

−→

n,m→+∞

Since we do not assume any source condition here, we are in the case r̃ = 0 of Proposition 19, which gives, for any λ > 0, 1 1 1 1 1/2 b −1/2 U + V + W + Z ∥Aρ u bλ − Aρ uλ ∥ρ ≤ ∥Σλ Σ ∥ n n m m op λ λ λ λ1/2 λ1/2 1 1 1 1 1/2 ≤ ∥Σλ ∥op Un + 3/2 Vn + Wm + 3/2 Zm , λ λ λ λ where we used

1/2 ∥Σλ ∥op 1/2 b −1/2 √ ∥Σλ Σλ ∥op ≤ .

λ Let τk = 2 log k. By point (ii) of Proposition 19, with probability at least 1 − 2e−τk , √ τk Uk ≤ CU √ . k P −τ k Since k≥2 e < ∞, the Borel–Cantelli lemma gives ! r log k Uk = O almost surely. k The same argument applies to Vk , Wk , Zk . Hence, almost surely, for all large enough n, m, ! ! r r log n log m , Wm , Zm = O . Un , Vn = O n m Let N = min(n, m). Since k 7→ log k/k is decreasing for large k, we obtain, almost surely, ! s Un Wm log N + =O → 0, λn,m λn,m λ2n,m N because eventually λn,m ≤ 1 and

log N log N ≤ 3 → 0. 2 λn,m N λn,m N

Similarly, Vn 3/2

λn,m Finally, since λn,m → 0, we have

Zm 3/2

λn,m

s =O

log N λ3n,m N

! → 0.

1/2

∥Σλn,m ∥op → ∥Σ1/2 ∥op ,

and therefore, almost surely, ∥Aρ u bλn,m − Aρ uλn,m ∥ρ

−→

n,m→+∞

Both the approximation error and the estimation error tend to 0, hence ∥Aρ u bλn,m − (hΠ , q Π )∥ρ

−→

n,m→+∞

Corollary 27. If Assumptions 1, 4, 5 and 8 hold, for any choice of λn,m > 0, for all n, m ∈ N, such that log N λn,m → 0, →0 as n, m → ∞, λ3n,m N where N = min(n, m), then, almost surely, the following convergence holds: ∥Aρ u bλn,m − (h, q)∥ρ

−→

n,m→+∞

That is, almost surely, for any ϵ > 0, there exists Nϵ ∈ N such that ∀(n, m) ∈ N2 , n, m ≥ Nϵ =⇒ ∥Aρ u bλn,m − (h, q)∥ρ < ϵ . Proof. Assumptions 1, 4 and 5 are taken in the whole Section C and are usually omitted to ease the reading — we make them explicit again for the present result. In particular such assumptions are needed in order to apply Proposition 26. Under the additional Assumption 8 we have (h, q) ∈ ran Aρ , thus (h, q) is equal to its orthogonal projection (hΠ , q Π ) on the subspace ran Aρ . The result then follows from Proposition 26.

C.7

Proof of Theorem 1

Let us consider the setting of Section 2. Under Assumption 2, we can define the bounded operator Aρ :

F −→ L2 (ρX ) × L2 (ρZ ) u 7−→ (u, Du).

(88)

We can then define (h, q) := Aρ u∗ . Our goal is to show that the hypotheses of Proposition 27 hold. Proposition 28. If Assumptions 2 and 3 hold, then Assumption 8 holds, i.e. (h, q) ∈ ran Aρ ⊂ L2 (ρX ) × L2 (ρZ ), or in other words, for any ϵ > 0, there exists u ∈ H such that ∥Aρ u − (h, q)∥ρ ≤ ϵ. Proof. We have u∗ ∈ F. Assumption 3 implies that the inclusion i:

H −→ F f 7−→ f

is well-defined and bounded, and has a dense range. In particular, u∗ ∈ ran i. We thus have Aρ u∗ ∈ Aρ ran i which implies since Aρ is continuous Aρ u∗ ∈ ran(Aρ ◦ i). Now recalling the definition (43) of Aρ , we see that Aρ = Aρ ◦ i, thus we have (h, q) ∈ ran Aρ .

Corollary 29. Under Assumptions 1 to 5, Proposition 27 holds, i.e. for any choice of λn,m > 0, for all n, m ∈ N, such that λn,m → 0,

log N →0 λ3n,m N

as n, m → ∞,

where N = min(n, m), then, almost surely, we have the asymptotic convergence ∥Aρ u bλn,m − (h, q)∥ρ

−→

n,m→+∞

By definition of Aρ , (h, q) and the norm ∥ · ∥ρ , this is equivalent to  ∥b uλn,m − u∗ ∥L2 (ρX ) −→ 0 n,m→∞ ∗ ∥Db uλn,m − Du ∥L2 (ρZ ) −→ n,m→∞

(89)

i.e. Theorem 1 is satisfied. Proof. Proposition 27 holds under Assumptions 1, 4, 5 and 8, so the only missing part is Assumption 8. Proposition 28 shows that Assumptions 2 and 3 imply Assumption 8.

Sobolev setting

In this appendix, we focus on the setting where D is a linear differential operator and F is a Sobolev space. Such a mathematical setting is typical for PDEs, which are ubiquitous in scientific applications. We show that our assumptions — formulated with generality in mind — are satisfied in such a setting, which demonstrates their practicality.

D.1

Generalities

We now consider a bounded Lipschitz domain Ω ⊂ Rd and set X = Ω. Fix s ∈ N, s ≥ 1. Let F = H s (Ω) and G = L2 (Ω), and consider the differential operator D of (8), where cα ∈ C(Ω). The first thing to show is that D defines a bounded operator from F to G. Lemma 30. The operator D is well-defined and bounded from H s (Ω) to L2 (Ω). Proof. Let u ∈ H s (Ω). By definition of H s (Ω), for every multi-index α with |α| ≤ s, the weak derivative ∂ α u belongs to L2 (Ω). Since cα ∈ C(Ω) and Ω is compact, we have cα ∈ L∞ (Ω). Hence cα ∂ α u ∈ L2 (Ω),

|α| ≤ s,

and ∥cα ∂ α u∥L2 (Ω) ≤ ∥cα ∥L∞ (Ω) ∥∂ α u∥L2 (Ω) .

Since the sum defining Du is finite, it follows that Du ∈ L2 (Ω). Thus D is well-defined from H s (Ω) to L2 (Ω). Moreover, we have X X ∥Du∥L2 (Ω) ≤ ∥cα ∂ α u∥L2 (Ω) ≤ ∥cα ∥L∞ (Ω) ∥∂ α u∥L2 (Ω) . |α|≤s

|α|≤s

Using the Cauchy-Schwarz inequality, we obtain

∥Du∥L2 (Ω) ≤ 

1/2

1/2 

 X

|α|≤s

∥cα ∥2L∞ (Ω) 

X  |α|≤s

∥∂ α u∥2L2 (Ω) 

Therefore ∥Du∥L2 (Ω) ≤ CD ∥u∥H s (Ω) , where

1/2

 CD := 

X |α|≤s

∥cα ∥2L∞ (Ω) 

Hence D : H s (Ω) → L2 (Ω) is bounded. Definition 3 (C0s universality). Let U be an open subset of Rd , let K : U × U → R be a kernel, and let H denote its associated RKHS. We say that K is C0s -universal on U if for any f ∈ C0s (U ) and any ϵ > 0, there exists fH ∈ H such that for every α ∈ Nd with |α| ≤ s, sup |∂ α f (x) − ∂ α fH (x)| ≤ ϵ,

x∈U

i.e. H ∩ C0s (U ) is dense in C0s (U ) for the norm ∥f ∥C0s (U ) := max|α|≤s supx∈U |∂ α f (x)|. Lemma 31. Consider Ω a bounded Lipschitz domain and F := H s (Ω). If K ∈ C 2s (X × X ) is the restriction to X of a C0s -universal kernel over Rd , then Assumption 3 holds, i.e. H is continuously and densely embedded in H s (Ω). Proof. The embedding is well-defined and bounded: since K ∈ C 2s (X × X ), Proposition 7 gives 1/2 ∥u∥H for |α| ≤ s, and since Ω is H ,→ C s (X ) with supx∈X |∂ α u(x)| ≤ supx∈X ∂1α ∂2α K(x, x) bounded this yields ∥u∥H s (Ω) ≤ c ∥u∥H for a constant c depending on K, s, d and |Ω|. It remains to prove density. Let us denote by K both the kernel on X and its extension to Rd . Let H0 denote the RKHS associated to K over Rd . We can see H as a closed subspace of H0 , defined by: H = span{Kx , x ∈ X } ⊂ H0 . Consider a function f ∈ H s (Ω), and let us try to approximate it by an element of H. Since Ω is a Lipschitz domain, there exists f˜ ∈ H s (Rd ) such that f˜|Ω = f (Stein, 1970, Theorem 5). Let ϵ > 0, by density of Cc∞ (Rd ) in H s (Rd ), there exists ϕ ∈ Cc∞ (Rd ) such that ∥ϕ − f˜∥H s < ϵ. Since K is C0s universal over Rd , there now exists fH0 ∈ H0 such that ∀|α| ≤ s, ∀x ∈ Rd , |∂ α fH0 (x) − ∂ α ϕ(x)| < ϵ. We thus have XZ d+s α α 2 |∂ fH0 (x) − ∂ ϕ(x)| dx ≤ |Ω|ϵ2 . d Ω |α|≤s

Let us denote by fH the orthogonal projection of fH0 onto H. Since for x ∈ Ω we have Kx ∈ H, we have fH (x) = ⟨fH , Kx ⟩H = ⟨fH0 , Kx ⟩H0 = fH0 (x). If we denote by ϕ|Ω the restriction of ϕ to Ω, we thus have s d+s ∥fH − ϕ|Ω ∥H s ≤ |Ω|ϵ. d 43

We can finally bound ∥fH − f ∥H s ≤ ∥fH − ϕ|Ω ∥H s + ∥ϕ|Ω − f ∥H s ≤ ∥fH − ϕ|Ω ∥H s + ∥ϕ − f˜∥H s s ! d+s |Ω| + 1 ϵ. ≤ d Lemma 32. If K ∈ C 2s (X × X ), and D is a linear differential operator of the form (8) then Assumption 5 holds. Proof. Since by assumption X is compact, we can define q κ = max ∥Kx ∥H = max K(x, x). x∈X

(90)

x∈X

Now assume that the differential operator is of the form (8), i.e. it can be written as X D= cα ∂ α ,

(91)

|α|≤s

for some integer s ≥ 1 and where the cα : X → R are continuous coefficient functions. For x, y ∈ X , let us write X X ′ D1 D2 K(x, y) = cα (x)cα′ (y)∂1α ∂2α K(x, y). |α|≤s |α′ |≤s

′

Since K is C 2s , for any |α| ≤ s, |α′ | ≤ s, we observe that ∂1α ∂2α K(x, y) is continuous and thus bounded on the compact X × X . Since cα (x) and cα′ (y) are bounded by assumption, we can also define q (92) κD = max ∥KzD ∥H = max D1 D2 K(z, z). z∈X

z∈X

D.2

In-domain sampling

We now turn to Assumption 2. Let us first clarify the meaning of this assumption. The spaces F and G are Hilbert spaces of functions on X , possibly defined only up to almost-everywhere equivalence, as in L2 or Sobolev spaces. Hence, if ρX is singular with respect to the reference measure defining F, the expression [u]ρX ∈ L2 (ρX ) need not be well-defined for an arbitrary u ∈ F. For instance, if F = L2 (X ) and ρX = δx , then [u]ρX would be determined by the pointwise value u(x), which is not defined for a general L2 -equivalence class. We therefore define the map F → L2 (ρX ) by density, starting from continuous representatives. Let DF := F ∩ C 0 (X ), where this is understood as the subspace of elements of F admitting a continuous representative. For u ∈ DF , the class [u]ρX is well-defined, and we have ∥u∥L2 (ρX ) < +∞ because X is compact, so [u]ρX ∈ L2 (ρX ). If ∥u∥L2 (ρX ) ≤ C∥u∥F for all u ∈ DF , and DF is dense in F , then the map u ∈ DF 7−→ [u]ρX ∈ L2 (ρX ) 44

extends uniquely to a bounded linear operator from F to L2 (ρX ). The corresponding map G → L2 (ρZ ) is defined analogously. In Proposition 33 below, we prove that Assumption 2 holds for in-domain sampling, where both X and Z take values inside Ω. Lemma 33. Assume that the distributions ρX and ρZ are absolutely continuous with respect to the Lebesgue measure on Ω, with bounded densities. Then, Assumption 2 holds. Proof. Since ρX and ρZ have bounded densities, we can bound the L2 (ρX ) and L2 (ρZ ) norms by the standard L2 (Ω) norm. Combined with the canonical embedding H s (Ω) ,→ L2 (Ω), we see that Assumption 2 is satisfied. Proposition 34. Under the standing assumptions of Section D.1, assume that ρX and ρZ have bounded densities with respect to Lebesgue measure of Ω, that K ∈ C 2s (X × X ) is the restriction of a C0s -universal kernel on Rd , and that Assumption 4 holds. Then all the assumptions of Theorem 1 are satisfied. Proof. Proposition 8 shows that Assumption 1 holds. Proposition 33 shows that Assumption 2 holds. Proposition 31 shows that Assumption 3 holds. Proposition 32 shows that Assumption 5 holds. Finally, Assumption 4 is a standard assumption that we take independently from the rest. All the assumptions of Theorem 1 are thus satisfied, and the theorem applies. D.2.1

Stronger convergence in the elliptic case

In this section, we consider an operator D of order s = 2 in divergence form: X d d X ∂ ∂ ∂ ij u(x) + bi (x) u(x) + c(x)u(x), a (x) Du(x) = − ∂xj ∂xi ∂xi

(93)

i=1

i,j=1

where aij ∈ C 1 (Ω), bi ∈ C(Ω) and c ∈ C(Ω) are coefficient functions and aij = aji for all i, j. Note in particular that D is of the form (8). Definition 4 (Evans, 2010, Section 6.1). We say that D is uniformly elliptic if there exists a constant θ > 0 such that d X aij (x)ξi ξj ≥ θ|ξ|2 i,j=1

for almost every x ∈ Ω and all ξ ∈ Rd . In the remainder of this subsection, we assume that D is uniformly elliptic. The following classical result on elliptic PDEs allows us to control higher-order norms of u (here, H 2 (V ) for an open set V compactly embedded in Ω) thanks to the norm of Du (here, L2 ). Proposition 35 (Evans, 2010, Section 6.3, Theorem 1). Let f ∈ L2 (Ω) and assume that u ∈ H 1 (Ω) 2 (Ω) and for each open set V satisfying is a weak solution of the PDE Du = f on Ω. Then u ∈ Hloc V ⊂ Ω, there exists CV > 0 such that we have ∥u∥H 2 (V ) ≤ CV (∥u∥L2 (Ω) + ∥f ∥L2 (Ω) ).

Corollary 36. Assume that ρX and ρZ are both absolutely continuous with respect to the Lebesgue measure on Ω, with densities bounded away from 0, so that the norms ∥ · ∥L2 (ρX ) and ∥ · ∥L2 (ρZ ) are stronger than ∥ · ∥L2 (Ω) . Then for any open V such that V ⊂ Ω, there exists CV > 0 such that for all u ∈ H 2 (Ω), we have ∥u∥H 2 (V ) ≤ CV (∥u∥L2 (ρX ) + ∥Du∥L2 (ρZ ) ). (94) Proof. For any u ∈ H 2 (Ω), if we denote f := Du, we have f ∈ L2 (Ω) which allows us to apply Proposition 35 with u and f and then bound ∥u∥L2 (Ω) and ∥Du∥L2 (Ω) by ∥u∥L2 (ρX ) and ∥Du∥L2 (ρZ ) respectively to obtain (94). The preceding results allow us to state the following Sobolev convergence result in the elliptic case. Corollary 37. Assume the hypotheses of Proposition 34 with s = 2. Assume moreover that D is uniformly elliptic and the densities of ρX and ρZ with respect to the Lebesgue measure of Ω, which are bounded from above by the hypotheses of Proposition 34, are also bounded away from 0, so that the norms ∥ · ∥L2 (ρX ) and ∥ · ∥L2 (ρZ ) are equivalent to ∥ · ∥L2 (Ω) . Then, for any sequence (λn,m ) satisfying (15), for every open V such that V ⊂ Ω, we almost surely have ∥b uλn,m − u∗ ∥H 2 (V )

−→

n,m→∞

Proof. By Proposition 34, Theorem 1 applies: for any sequence (λn,m ) satisfying (15), we almost surely have  ∥b uλn,m − u∗ ∥L2 (ρX ) −→ 0 n,m→∞

∥Db uλn,m − Du∗ ∥L2 (ρZ )

−→

n,m→∞

For any open set V such that V ⊂ Ω, using Proposition 36 with vn,m = u bλn,m − u∗ ∈ H 2 (Ω) yields the result.

D.3

Boundary sampling

In Proposition 38 below, we prove that Assumption 2 holds for the boundary sampling scenario, where X takes values only on the boundary ∂Ω (while Z still takes values inside Ω). Such a scenario is slightly less straightforward as we cannot rely on absolute continuity (with respect to the Lebesgue measure on Ω), and we instead rely on results from trace theory. Lemma 38. Assume that ρX is absolutely continuous with respect to the Hausdorff measure of ∂Ω. Assume that ρZ is absolutely continuous with respect to the Lebesgue measure on Ω. Assume that both densities are bounded. Then, Assumption 2 holds. Proof. Since ρZ has a bounded density, we can bound the L2 (ρZ ) norm by the L2 (Ω) norm, so the embedding G → L2 (ρZ ) is indeed well-defined and bounded. Since ρX has a bounded density, we can bound the L2 (ρX ) norm by the standard L2 (∂Ω) norm. Since Ω is a Lipschitz domain, we can define the trace operator T : H 1 (Ω) → H 1/2 (∂Ω), which for any u ∈ H 1 (Ω) ∩ C 0 (Ω), coincides with the restriction of u to the boundary. By combining the trace operator with the embeddings H s (Ω) ,→ H 1 (Ω) (before the trace operator) and H 1/2 (∂Ω) ,→ L2 (∂Ω) (after the trace operator) and finally the bounded map L2 (∂Ω) → L2 (ρX ), we get that Assumption 2 holds. 46

Proposition 39. Under the standing assumptions of Section D.1, assume that ρX is absolutely continuous with respect to the Hausdorff measure of ∂Ω, with bounded density, and assume that ρZ is absolutely continuous with respect to the Lebesgue measure on Ω, with bounded density. Assume moreover that K ∈ C 2s (X × X ) is the restriction of a C0s -universal kernel on Rd , and that Assumption 4 holds. Then all the assumptions of Theorem 1 are satisfied. Proof. Proposition 8 shows that Assumption 1 holds. Proposition 38 shows that Assumption 2 holds. Proposition 31 shows that Assumption 3 holds. Proposition 32 shows that Assumption 5 holds. Finally, Assumption 4 is a standard assumption that we take independently from the rest. All the assumptions of Theorem 1 are thus satisfied, and the theorem applies. D.3.1

Strong convergence in the elliptic case

Throughout this subsubsection, we assume that s = 2 and that D is the operator defined in (93). We further assume that aij , bi , c ∈ C ∞ (Ω), that D is uniformly elliptic over Ω (cf. Definition 4), and that the boundary ∂Ω is a (d − 1)-dimensional smooth manifold, Ω being locally on one side of ∂Ω (such regularity is needed to apply the regularity estimates from Lions and Magenes (2012)). We consider the problem ( Du(x) = q(x) x ∈ Ω (95) u(x) = h(x) x ∈ ∂Ω, where (q, h) ∈ L2 (Ω) × H 3/2 (∂Ω). Let us assume that 0 is not a Dirichlet eigenvalue for the operator D in Ω, which, by the Fredholm alternative (see Evans, 2010), guarantees that for any (q, h) ∈ L2 (Ω) × H 3/2 (∂Ω), there exists a unique solution u ∈ H 2 (Ω) to (95). This is the case for instance for the Laplacian D = −∆. Proposition 40 (Lions and Magenes, 2012, Theorem 7.4). We have, for any (q, h) ∈ L2 (Ω) × H 3/2 (∂Ω), for any u ∈ H 2 (Ω) solution of (95), (96) ∥u∥H 1/2 (Ω) ≤ c ∥h∥L2 (∂Ω) + ∥q∥Ξ−3/2 (Ω) , where c > 0 and ∥ · ∥Ξ−3/2 (Ω) is a norm weaker than ∥ · ∥L2 (Ω) . Remark 3. Given the regularity of u, q and h, higher order estimates also hold, such as ∥u∥H 2 (Ω) ≤ c ∥h∥H 3/2 (∂Ω) + ∥q∥L2 (Ω) . However, although h is more regular than L2 , we only have the L2 convergence of T u to h, so the estimate (96) is the strongest we can use in our context. Corollary 41. Assume that ρX (resp. ρZ ) is absolutely continuous with respect to the (d − 1)dimensional Hausdorff measure on ∂Ω (resp. the Lebesgue measure on Ω), and assume that both densities are bounded away from 0. Then there exists C > 0 such that for any u ∈ H 2 (Ω), we have ∥u∥H 1/2 (Ω) ≤ C(∥ T u∥L2 (ρX ) + ∥Du∥L2 (ρZ ) ).

(97)

Proof. Let u ∈ H 2 (Ω). If we denote q = Du and h = T u, we have q ∈ L2 (Ω), h ∈ H 3/2 (∂Ω) and by definition, u is a solution of (95) for that specific choice of (q, h). Because the density of ρX is bounded away from 0, we can write ∥h∥L2 (∂Ω) = ∥ T u∥L2 (∂Ω) ≤ c1 ∥ T u∥L2 (ρX ) .

The L2 (Ω) norm being stronger than the Ξ−3/2 (Ω) norm and the density of ρZ being bounded away from 0, we can write ∥q∥Ξ−3/2 (Ω) = ∥Du∥Ξ−3/2 (Ω) ≤ c2 ∥Du∥L2 (Ω) ≤ c3 ∥Du∥L2 (ρZ ) . Combining these with Proposition 40 proves (97). Corollary 42. Assume the hypotheses of Proposition 39 with s = 2, together with the standing assumptions of Section D.3.1. Assume moreover that the densities of ρX and ρZ are bounded away from zero, so that the norms ∥ · ∥L2 (ρX ) and ∥ · ∥L2 (ρZ ) are respectively equivalent to ∥ · ∥L2 (∂Ω) and ∥ · ∥L2 (Ω) . Then, for any sequence (λn,m ) satisfying (15), we have that, almost surely ∥b uλn,m − u∗ ∥H 1/2 (Ω)

n,m→+∞

∥b uλn,m − u∗ ∥L2 (Ω)

n,m→+∞

−→

In particular, we have

Proof. As we have shown above, in this setting, Theorem 1 holds, i.e. for any regularizing sequence (λn,m ) satisfying (15), we have almost surely  ∥b uλn,m − h∥L2 (ρX ) −→ 0 n,m→∞

∥Db uλn,m − q∥L2 (ρZ )

−→

n,m→∞

Considering vn,m = u bλn,m − u∗ , we can rewrite this as  ∥ T vn,m ∥L2 (ρX )

n,m→∞

∥Dvn,m ∥L2 (ρZ )

n,m→∞

−→

Using Proposition 41 with vn,m , we immediately obtain ∥vn,m ∥H 1/2 (Ω)

−→

n,m→+∞

and in particular since the H 1/2 norm is stronger than the L2 norm, we obtain ∥vn,m ∥L2 (Ω)

−→

n,m→+∞

Proposition 42 establishes H 1/2 convergence on Ω. If we restrict to a smaller set V ⊂ Ω, analogously to Proposition 37 in Section D.2, we can prove a stronger (H 2 ) convergence on V . Corollary 43. Assume the hypotheses of Proposition 39 with s = 2, together with the standing assumptions of Section D.3.1. Assume moreover that the densities of ρX and ρZ are bounded away from zero, so that the norms ∥ · ∥L2 (ρX ) and ∥ · ∥L2 (ρZ ) are respectively equivalent to ∥ · ∥L2 (∂Ω) and ∥ · ∥L2 (Ω) . Then, for any regularizing sequence (λn,m ) satisfying (15), we have almost surely that for any open V such that V ⊂ Ω, we have ∥b uλn,m − u∗ ∥H 2 (V ) 48

−→

n,m→∞

Proof. Let us consider again vn,m = u bλn,m − u∗ , as in the proof of Proposition 42. We established that  ∥ T vn,m ∥L2 (ρX ) −→ 0 n,m→∞

∥Dvn,m ∥L2 (ρZ )

n,m→∞

−→

∥vn,m ∥L2 (Ω)

−→

and that n,m→+∞

which allows us to use Proposition 35, and the lower density bound for ρZ , to conclude that for any open V such that V ⊂ Ω, we have ∥b uλn,m − u∗ ∥H 2 (V )

D.4

−→

n,m→∞

Rates for the Laplacian on periodic functions

In this section, we consider the example of the Laplacian on the torus, i.e. applied to Sobolev spaces of periodic functions on [0, 1]d . We use this example to illustrate the source condition for the rates provided in Theorem 2 and Proposition 5. Even though it is a Sobolev setting, periodicity creates differences with the rest of Section D, and in particular, the results of Section D.1 do not apply. Let Ω = (0, 1)d , let X = Ω = [0, 1]d . Throughout this section, we identify a function u on X satisfying periodic boundary conditions with a function on the torus Td = Rd /Zd , and we write its Fourier expansion as Z X 2πik·x u(x) = u ek e , u ek := u(x) e−2πik·x dx. [0,1]d

k∈Zd

For τ ≥ 0, we define the Sobolev space of periodic functions of order τ as o n X 2 τ τ |e u | < ∞ , λk := 1 + 4π 2 |k|2 , λ (X ) := u ∈ L2 (X ) : ∥u∥2Hper := Hper τ k k k∈Zd

P τ endowed with the inner product ⟨u, v⟩Hper = k∈Zd λτk u ek vek . The weight λτk is equivalent to the more common (1 + |k|2 )τ , so this choice only rescales the norm; we fix it because λk is exactly the symbol of Id − ∆, i.e. λk u ek = (u^ − ∆u) , which simplifies the computations below. Note that all the k

τ statements of this section are invariant under replacing ∥ · ∥Hper by an equivalent norm, except for τ (X ) can be identified the exact values of the constants κ, κD in Proposition 44. In particular, Hper τ with the subspace of H (X ) whose elements satisfy periodic boundary conditions, and for τ > d/2 P −τ 1/2 τ it is an RKHS, since by Cauchy–Schwarz |u(x)| ≤ ∥u∥Hper < ∞ for every x ∈ X . k λk Pd ∂2 2 2 Let us consider F = Hper (X ), G = L (X ) and D = ∆ = i=1 (∂xi )2 the Laplacian, which τ (X ), with τ > d/2 + 2. In particular, defines a continuous operator from F to G. Let H = Hper since τ > d/2, H is an RKHS. Let us consider ρX = ρZ = Unif([0, 1]d ) (which corresponds to the in-domain sampling setting of Section 4.3, up to the periodicity difference). In particular L2 (ρX ) = L2 (ρZ ) = L2 ([0, 1]d ).

Lemma 44. Assumptions 1 to 3 and 5 hold.

Proof. We prove Assumptions 1 and 5 together. For u ∈ H, write the Fourier expansions X X |k|2 u ek e2πik·x . u ek e2πik·x , ∆u(x) = −4π 2 u(x) = k∈Zd

k∈Zd

By Cauchy–Schwarz, 1/2

 |u(x)| ≤ 

 λ−τ k

k∈Zd

τ , ∥u∥Hper

and similarly 1/2

 |∆u(x)| ≤ 4π 2 

X k∈Zd

 |k|4 λ−τ k

τ . ∥u∥Hper

The first series is finite since τ > d/2, and the second one is finite since −(τ −2)

|k|4 λ−τ k ≲ λk

and τ − 2 > d/2. Thus, since τ > d/2 + 2, both point evaluations u 7→ u(x),

u 7→ ∆u(x)

are bounded uniformly in x. Hence Assumption 1 holds, and Assumption 5 holds with  1/2  1/2 X X  ,  . κ= λ−τ κD = 4π 2  |k|4 λ−τ k k k∈Zd

k∈Zd

The remaining assumptions are standard. Since 2 Hper (X ) ,→ L2 (X ) and G = L2 (X ),

and since ρX = ρZ = Unif([0, 1]d ), Assumption 2 holds. Finally, s 2 Hper (X ) ,→ Hper (X )

is continuous and dense, by density of smooth periodic functions in Sobolev spaces. Hence Assumption 3 holds. In particular, Assumption 2 allows us to define, as in Section 4.5, the operators Aρ : F → L2 (ρX ) × L2 (ρZ )

and

u 7→ (u, ∆u)

Aρ = Aρ ◦ i,

where i : H → F, u 7→ u is the canonical Sobolev embedding. We can now consider (h, q) = Aρ u∗ , denote L = Aρ A∗ρ and C = A∗ρ Aρ . The following lemma allows us to interpret the source condition on (h, q) as an equivalent source condition on the rescaled target C 1/2 u∗ . Lemma 45. There exists an isometry U : F → L2 (ρX ) × L2 (ρZ ) with range ran Aρ and a bounded positive self-adjoint operator T : F → F such that L = U T U ∗ , and furthermore (h, q) ∈ L2 (ρX ) × L2 (ρZ ) satisfies a source condition with L if and only if C 1/2 u∗ ∈ F satisfies an equivalent source condition defined by T . 50

Proof. The map Aρ is a continuous injective linear map. Moreover, it is a topological isomorphism from F onto its image 2 ran(Aρ ) = {(u, ∆u) : u ∈ Hper (X )} ⊂ L2 (ρX ) × L2 (ρZ ).

Indeed, one can easily check using Fourier decomposition that ∥u∥2L2 (X ) + ∥∆u∥2L2 (X ) ≍ ∥u∥2Hper 2 . This shows that C = A∗ρ Aρ is an isomorphism. Now use the polar decomposition of Aρ . We may write Aρ = U C 1/2 , where U : F → L2 (ρX ) × L2 (ρZ ) is an isometry with range ran Aρ . Consequently, recalling Aρ = Aρ i, L = Aρ A∗ρ = U C 1/2 ii∗ C 1/2 U ∗ . Define T := C 1/2 ii∗ C 1/2 : F → F. Then T is bounded, self-adjoint, and positive, and by the continuous functional calculus, Lr = (U T U ∗ )r = U T r U ∗ , where we used the fact that U ∗ U = IF . Assume now that (h, q) = Lr (f, g), where we may assume (f, g) = U v, with v ∈ F (if that is not the case we may simply replace (f, g) by its orthogonal projection onto ran Aρ , which will leave its image by Lr unchanged, as U ∗ annihilates (ran Aρ )⊥ ). (h, q) = Lr (f, g) ⇐⇒ Aρ u∗ = U T r U ∗ (f, g) ⇐⇒ C 1/2 u∗ = T r v

(98)

with v := U ∗ (f, g). We thus see that (h, q) satisfies a source condition with L if and only if C 1/2 u∗ satisfies an equivalent source condition defined by T . Lemma 46. Let r > 0 and set σr := 2 + 2r(τ − 2). Then there exists v ∈ F such that C 1/2 u∗ = T r v if and only if σr u∗ ∈ Hper (X ).

Proof. The operator C = A∗ρ Aρ is the Fourier multiplier ] = ck u (Cu) ek , k

ck :=

1 + 16π 4 |k|4 . λ2k

t (X ). Since ck ≍ 1, the operators C 1/2 and C −1/2 are bounded isomorphisms on every Hper

Next, for the embedding τ 2 i : Hper (X ) → Hper (X ),

one checks from the Fourier definitions of the inner products that ∗ u) = λ2−τ u ^ (ii ek . k k

Hence T = C 1/2 ii∗ C 1/2 is the Fourier multiplier ] (T u)k = ck λ2−τ ek ≍ λ2−τ ek . k u k u Thus

r u) ≍ λ−r(τ −2) u ^ (T ek . k k

Set w := C 1/2 u∗ . 2 (X ) such that Using the Fourier coefficients of T , we see that there exists v ∈ F = Hper

w = T rv if and only if X

2+2r(τ −2)

λk

k∈Zd

|w ek |2 < ∞,

which is equivalent to 2+2r(τ −2) w ∈ Hper (X ).

Since C 1/2 is an isomorphism on Sobolev spaces, this is equivalent to 2+2r(τ −2) u∗ ∈ Hper (X ).

The result follows. We can now restate and prove Proposition 6 from Section 4.5. Proposition 47 (Proposition 6, restated). Let u∗ ∈ F , let (h, q) = Aρ u∗ , and let r ∈ (0, 1]. σr (X ), where Then Assumption 6 holds with exponent r, i.e. (h, q) ∈ ran Lr , if and only if u∗ ∈ Hper 1 σr := 2 + 2r(τ − 2). In particular, for r = 2 , Assumption 6 is equivalent to u∗ ∈ H. Proof. Fix r ∈ (0, 1]. By Proposition 45, and in particular by (98), since L = U T U ∗ with U an isometry with range ran Aρ , we have (h, q) ∈ ran Lr if and only if there exists v ∈ F such that σr (X ). C 1/2 u∗ = T r v. By Proposition 46, such a v exists if and only if u∗ ∈ Hper s (X ) = H, which is the Finally, for r = 12 we have σr = τ , so the condition reads u∗ ∈ Hper well-specified setting.

Useful results

Proposition 48 (Hoeffding inequality in separable Hilbert spaces). Take a family ξ1 , . . . , ξn : Ω → H of independent zero mean random variables such that ∥ξi ∥H ≤ c, then for all ϵ > 0 # " 2 n ϵ n 1X ξi > ϵ ≤ 2 exp − 2 P n 4c i=1

i.e. for all τ > 0, with probability at least 1 − 2e−τ n

1X ξi n i=1

√ 2c τ ≤ √ . n

For the proof, we refer the reader to Yurinsky (2006); Pinelis (1994, 1999). Proposition 49 (Tropp’s concentration inequality, Rudi et al., 2013, Theorem A.1). Let (Zi )1≤i≤n be independent copies of the random variable Z with values in the space of bounded self-adjoint operators S(H) over a separable Hilbert space H. Define T := E[Z], and let there be S ∈ S(H) such that E[(Z − T )2 ] ⪯ S, and a finite number R such that ∥Z∥op ≤ R almost surely. Define the quantities α := ∥S∥1 /∥S∥op and σ 2 := ∥S∥op . Then, for 0 < δ ≤ α, it holds   r n  1X 2 3βσ  βR Zi − T ≤ + ≥ 1−δ, (99) P  n n n  i=1

where β := 23 log 4α δ . Proposition 49 originally comes from Tropp (2012), but in its current form is a simple reproduction of Rudi et al. (2013, Theorem A.1).

Additional Experiment Information

F.1

FEM Comparisons

We describe here in more detail the Poisson PDE used for the FEM experiments. Define a sequence of support points {si }Si=1 belonging to the disk in R2 . Letting ν ∈ {0.5, 1.5, 2.5, ∞} denote the smoothness parameter, we define f as ( P (1 − ∥x∥2 ) Si=1 Kν (x, si ) if ∥x∥ ≤ 1 f (x) = 0 otherwise where Kν is the Matérn kernel with parameter ν (note that ν = 0.5 corresponds to the Laplacian kernel and ν = ∞ to the Gaussian). The Poisson PDE is then defined on a domain Ω = {x : R2 | ∥x∥ < 1} with boundary ∂Ω = {x : R2 | ∥x∥ = 1}: ( ∆u(x) = q(x) x ∈ Ω u(x) = h(x) x ∈ ∂Ω with q(x) = ∆f (x) and h(x) = f (x). The four Laplacians corresponding to ν = 0.5, 1.5, 2.5, ∞ are shown in Fig. 5. 53

Figure 5: Four Matérn Laplacians with increasing smoothness.

References Eleonora Arnone, Alois Kneip, Fabio Nobile, and Laura M Sangalli. Some first results on the consistency of spatial regression with partial differential equation regularization. Statistica Sinica, 32(1):209–238, 2022. Simon Arridge, Peter Maass, Ozan Öktem, and Carola-Bibiane Schönlieb. Solving inverse problems using data-driven models. Acta numerica, 28, 2019. Laura Azzimonti, Laura M Sangalli, Piercesare Secchi, Maurizio Domanin, and Fabio Nobile. Blood flow velocity field estimation via spatial regression with pde penalization. Journal of the American Statistical Association, 110(511):1057–1071, 2015. Ricardo Baptista, Edoardo Calvello, Matthieu Darcy, Houman Owhadi, Andrew M. Stuart, and Xianjin Yang. Solving roughly forced nonlinear PDEs via misspecified kernel methods and neural networks, 2025. Igor A. Baratta, Joseph P. Dean, Jørgen S. Dokken, Michal Habera, Jack S. Hale, Chris N. Richardson, Marie E. Rognes, Matthew W. Scroggs, Nathan Sime, and Garth N. Wells. DOLFINx: the next generation FEniCS problem solving environment. preprint, 2023. doi: 10.5281/zenodo.10447666. Pau Batlle, Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M Stuart. Error analysis of kernel/gp methods for nonlinear and parametric pdes. Journal of Computational Physics, 520: 113488, 2025. Mikhail Belkin, Partha Niyogi, and Vikas Sindhwani. Manifold regularization: A geometric framework for learning from labeled and unlabeled examples. Journal of machine learning research, 7 (11), 2006.

George David Birkhoff. General mean value and remainder theorems with applications to mechanical differentiation and quadrature. Transactions of the American Mathematical Society, 7(1):107–136, 1906. Gilles Blanchard and Nicole Mücke. Optimal rates for regularization of statistical inverse learning problems. Foundations of Computational Mathematics, 18(4):971–1013, 2018. Vivien Cabannes, Loucas Pillaud-Vivien, Francis Bach, and Alessandro Rudi. Overcoming the curse of dimensionality with laplacian regularization in semi-supervised learning. Advances in Neural Information Processing Systems, 34:30439–30451, 2021. Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis. Physicsinformed neural networks (pinns) for fluid mechanics: A review. Acta Mechanica Sinica, 37(12): 1727–1738, 2021. Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squares algorithm. Foundations of Computational Mathematics, 7:331–368, 2007. Yifan Chen, Bamdad Hosseini, Houman Owhadi, and Andrew M Stuart. Solving and learning nonlinear pdes with gaussian processes. Journal of Computational Physics, 447:110668, 2021. Yifan Chen, Houman Owhadi, and Florian Schäfer. Sparse cholesky factorization for solving nonlinear pdes via gaussian processes. Mathematics of Computation, 94(353):1235–1280, 2025. Salvatore Cuomo, Vincenzo Schiano Di Cola, Fabio Giampaolo, Gianluigi Rozza, Maziar Raissi, and Francesco Piccialli. Scientific machine learning through physics–informed neural networks: Where we are and what’s next. Journal of Scientific Computing, 92(3), 2022. Wojciech M Czarnecki, Simon Osindero, Max Jaderberg, Grzegorz Swirszcz, and Razvan Pascanu. Sobolev training for neural networks. Advances in neural information processing systems, 30, 2017. Filip De Roos, Alexandra Gessner, and Philipp Hennig. High-dimensional gaussian process inference with derivatives. In International Conference on Machine Learning, pages 2535–2545. PMLR, 2021. Tim De Ryck, Ameya D Jagtap, and Siddhartha Mishra. Error estimates for physics-informed neural networks approximating the navier–stokes equations. IMA Journal of Numerical Analysis, 44(1):83–119, 2024. Ernesto De Vito, Andrea Caponnetto, and Lorenzo Rosasco. Model selection for regularized leastsquares algorithm in learning theory. Foundations of Computational Mathematics, 5(1):59–85, 2005a. Ernesto De Vito, Lorenzo Rosasco, Andrea Caponnetto, Umberto De Giovannini, Francesca Odone, and Peter Bartlett. Learning from examples as an inverse problem. Journal of Machine Learning Research, 6(5), 2005b. Nathan Doumèche, Francis Bach, Claire Boyer, and Gérard Biau. Physics-informed machine learning as a kernel method. In Proceedings of the Thirty Seventh Annual Conference on Learning Theory, 2024.

Nathan Doumèche, Francis Bach, Gérard Biau, and Claire Boyer. Fast kernel methods: Sobolev, physics-informed, and additive models. arXiv preprint arXiv:2509.02649, 2025a. Nathan Doumèche, Francis Bach, Gérard Biau, and Claire Boyer. Physics-informed kernel learning. Journal of Machine Learning Research, 26(124):1–39, 2025b. Nathan Doumèche, Gérard Biau, and Claire Boyer. On the convergence of PINNs. Bernoulli, 31 (3):2127 – 2151, 2025c. doi: 10.3150/24-BEJ1799. Heinz Werner Engl, Martin Hanke, and Andreas Neubauer. Regularization of inverse problems, volume 375. Springer Science & Business Media, 1996. David Eriksson, Kun Dong, Eric Lee, David Bindel, and Andrew G Wilson. Scaling gaussian process regression with derivatives. Advances in neural information processing systems, 31, 2018. Lawrence C Evans. Partial differential equations, volume 19. American mathematical society, 2010. Gregory E Fasshauer. Solving partial differential equations by collocation with radial basis functions. In Proceedings of Chamonix, volume 1997, pages 1–8, 1996. Simon Fischer and Ingo Steinwart. Sobolev norm learning rates for regularized least-squares algorithms. Journal of Machine Learning Research, 21(205):1–38, 2020. Carsten Franke and Robert Schaback. Convergence order estimates of meshless collocation methods using radial basis functions. Advances in computational mathematics, 8:381–399, 1998a. Carsten Franke and Robert Schaback. Solving partial differential equations by collocation using radial basis functions. Applied Mathematics and Computation, 93(1):73–82, 1998b. Xiang Fu, Zhenghao Wu, Wujie Wang, Tian Xie, Sinan Keten, Rafael Gomez-Bombarelli, and Tommi Jaakkola. Forces are not enough: Benchmark and critical evaluation for machine learning force fields with molecular simulations. Transactions on Machine Learning Research, 2023. Martin Hanke. Regularization with differential operators: an iterative approach. Numerical functional analysis and optimization, 13(5-6):523–540, 1992. M Ch Hermite and M Borchardt. Sur la formule d’interpolation de lagrange. Journal für die reine und angewandte Mathematik (Crelles Journal), 1878(84):70–79, 1878. Edward J Kansa. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics—ii solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers & mathematics with applications, 19(8-9):147–161, 1990. George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6), 2021. Karthik Kashinath, Mustafa Mustafa, Adrian Albert, Jean-Luc Wu, C Jiang, Soheil Esmaeilzadeh, Kamyar Azizzadenesheli, R Wang, Ashesh Chattopadhyay, Aakanksha Singh, et al. Physicsinformed machine learning: case studies for weather and climate modelling. Philosophical Transactions of the Royal Society A, 379(2194):20200093, 2021. George Kimeldorf and Grace Wahba. Some results on tchebycheffian spline functions. Journal of mathematical analysis and applications, 33(1):82–95, 1971. 56

Georgios Kissas, Yibo Yang, Eileen Hwuang, Walter R Witschey, John A Detre, and Paris Perdikaris. Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from noninvasive 4d flow mri data using physics-informed neural networks. Computer methods in applied mechanics and engineering, 358:112623, 2020. Alexander Kovacs, Lukas Exl, Alexander Kornell, Johann Fischbacher, Markus Hovorka, Markus Gusenbauer, Leoni Breth, Harald Oezelt, Masao Yano, Noritsugu Sakuma, et al. Conditional physics informed neural networks. Communications in Nonlinear Science and Numerical Simulation, 104:106041, 2022. Aditi S. Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Characterizing possible failure modes in physics-informed neural networks. Advances in Neural Information Processing Systems, 34, 2021. Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021. URL https://openreview.net/forum?id=c8P9NQVtmnO. Junhong Lin, Alessandro Rudi, Lorenzo Rosasco, and Volkan Cevher. Optimal rates for spectral algorithms with least-squares regression over hilbert spaces. Applied and Computational Harmonic Analysis, 48(3):868–890, 2020. Jacques Louis Lions and Enrico Magenes. Non-homogeneous boundary value problems and applications: Vol. 1, volume 181. Springer Science & Business Media, 2012. Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nature machine intelligence, 3(3):218–229, 2021. Charles A Micchelli, Yuesheng Xu, and Haizhang Zhang. Universal kernels. Journal of Machine Learning Research, 7(12), 2006. Siddhartha Mishra and Roberto Molinaro. Estimates on the generalization error of physics-informed neural networks for approximating pdes. IMA Journal of Numerical Analysis, 43(1):1–43, 2023. Francis J. Narcowich, Joseph D. Ward, and Holger Wendland. Sobolev error estimates and a Bernstein inequality for scattered data interpolation via radial basis functions. Constructive Approximation, 24:175–186, 2006. Houman Owhadi. Bayesian numerical homogenization. Multiscale Modeling & Simulation, 13(3): 812–828, 2015. Misha Padidar, Xinran Zhu, Leo Huang, Jacob Gardner, and David Bindel. Scaling gaussian processes with derivative information using variational inference. Advances in Neural Information Processing Systems, 34:6442–6453, 2021. Iosif Pinelis. Optimum bounds for the distributions of martingales in banach spaces. The Annals of Probability, pages 1679–1706, 1994. Iosif Pinelis. Correction:“optimum bounds for the distributions of martingales in banach spaces”[ann. probab. 22 (1994), no. 4, 1679–1706; mr 96b: 60010]. The Annals of Probability, 27(4):2119–2119, 1999. 57

T. Poggio and F. Girosi. Regularization algorithms for learning that are equivalent to multilayer networks. Science, 247(4945), 1990. doi: 10.1126/science.247.4945.978. Alfio Quarteroni, Paola Gervasio, and Francesco Regazzoni. Combining physics-based and datadriven models: advancing the frontiers of research with scientific machine learning. Mathematical Models and Methods in Applied Sciences, 35(04):905–1071, 2025. doi: 10.1142/ S0218202525500125. Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal differential equations for scientific machine learning. arXiv preprint arXiv:2001.04385, 2020. Nasim Rahaman, Aristide Baratin, Devansh Arpit, Felix Draxler, Min Lin, Fred Hamprecht, Yoshua Bengio, and Aaron Courville. On the spectral bias of neural networks. In International conference on machine learning. PMLR, 2019. Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in neural information processing systems, 20, 2007. Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Machine learning of linear differential equations using gaussian processes. Journal of Computational Physics, 348:683–693, 2017a. Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations, 2017b. Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Numerical gaussian processes for timedependent and nonlinear partial differential equations. SIAM Journal on Scientific Computing, 40(1):A172–A198, 2018. Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019. Majid Rasht-Behesht, Christian Huber, Khemraj Shukla, and George Em Karniadakis. Physicsinformed neural networks (pinns) for wave propagation and full waveform inversions. Journal of Geophysical Research: Solid Earth, 127(5):e2021JB023120, 2022. Pratik Rathore, Weimu Lei, Zachary Frangella, Lu Lu, and Madeleine Udell. Challenges in training pinns: A loss landscape perspective. arXiv preprint arXiv:2402.01868, 2024. Alessandro Rudi, Guillermo D Canas, and Lorenzo Rosasco. On the sample complexity of subspace learning. Advances in Neural Information Processing Systems, 26, 2013. Alessandro Rudi, Raffaello Camoriano, and Lorenzo Rosasco. Less is more: Nyström computational regularization. Advances in neural information processing systems, 28, 2015. Francisco Sahli Costabal, Yibo Yang, Paris Perdikaris, Daniel E Hurtado, and Ellen Kuhl. Physicsinformed neural networks for cardiac activation mapping. Frontiers in Physics, 8:42, 2020. Laura M Sangalli. Spatial regression with partial differential equation regularisation. International Statistical Review, 89(3):505–531, 2021. Daniela Schräder and Holger Wendland. An extended error analysis for a meshfree discretization method of darcy’s problem. SIAM Journal on Numerical Analysis, 50(2):838–857, 2012. 58

Lei Shi, Xin Guo, and Ding-Xuan Zhou. Hermite learning with gradient data. Journal of computational and applied mathematics, 233(11):3046–3059, 2010. Yeonjong Shin, Jérôme Darbon, and George Em Karniadakis. On the convergence of physics informed neural networks for linear second-order elliptic and parabolic type pdes. Communications in Computational Physics, 28(5), 2020. doi: 10.4208/cicp.oa-2020-0193. Carl-Johann Simon-Gabriel and Bernhard Schölkopf. Kernel distribution embeddings: Universal kernels, characteristic kernels and kernel metrics on distributions. Journal of Machine Learning Research, 19(44):1–29, 2018. Dejan Slepcev and Matthew Thorpe. Analysis of p-laplacian regularization in semisupervised learning. SIAM Journal on Mathematical Analysis, 51(3):2085–2120, 2019. Steve Smale and Ding-Xuan Zhou. Learning theory estimates via integral operators and their approximations. Constructive approximation, 26(2):153–172, 2007. Alex J. Smola, Bernhard Schölkopf, and Klaus-Robert Müller. The connection between regularization operators and support vector kernels. Neural Networks, 11(4), 1998. doi: https: //doi.org/10.1016/S0893-6080(98)00032-X. Bharath K Sriperumbudur, Kenji Fukumizu, and Gert RG Lanckriet. Universality, characteristic kernels and rkhs embedding of measures. Journal of Machine Learning Research, 12(7), 2011. Elias M Stein. Singular integrals and differentiability properties of functions. Princeton university press, 1970. Ingo Steinwart and Andreas Christmann. Support Vector Machines. Springer Publishing Company, Incorporated, 1st edition, 2008. ISBN 0387772413. Ingo Steinwart, Don R Hush, and Clint Scovel. Optimal rates for regularized least squares regression. In COLT, pages 79–93, 2009. Juan Diego Toscano, Vivek Oommen, Alan John Varghese, Zongren Zou, Nazanin Ahmadi Daryakenari, Chenxi Wu, and George Em Karniadakis. From pinns to pikans: recent advances in physicsinformed machine learning. Machine Learning for Computational Science and Engineering, 1(1), 2025. doi: 10.1007/s44379-025-00015-1. Joel A Tropp. User-friendly tools for random matrices: An introduction, 2012. Grace Wahba. Spline models for observational data. SIAM, 1990. Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM Journal on Scientific Computing, 43(5):A3055– A3081, 2021. Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 449, 2022. ISSN 0021-9991. doi: https: //doi.org/10.1016/j.jcp.2021.110768. Wenjia Wang and Bing-Yi Jing. Gaussian process regression: Optimality, robustness, and relationship with kernel ridge regression. Journal of Machine Learning Research, 23(193):1–67, 2022. Holger Wendland. Scattered data approximation, volume 17. Cambridge university press, 2004. 59

George Wynne, François-Xavier Briol, and Mark Girolami. Convergence guarantees for gaussian process means with misspecified likelihoods and smoothness. Journal of Machine Learning Research, 22(123):1–40, 2021. Vadim Yurinsky. Sums and Gaussian vectors. Springer, 2006. Marius Zeinhofer, Rami Masri, and Kent-André Mardal. A unified framework for the error analysis of physics-informed neural networks. IMA Journal of Numerical Analysis, 45(5), 2025. Haobo Zhang, Yicheng Li, Weihao Lu, and Qian Lin. On the optimality of misspecified kernel ridge regression. In International Conference on Machine Learning, pages 41331–41353. PMLR, 2023. Zhiyuan Zhao, Xueying Ding, and B Aditya Prakash. Pinnsformer: A transformer-based framework for physics-informed neural networks. In International Conference on Learning Representations (ICLR), 2024. Dengyong Zhou and Bernhard Schölkopf. Regularization on discrete spaces. In Joint Pattern Recognition Symposium, pages 361–368. Springer, 2005. Ding-Xuan Zhou. Derivative reproducing properties for kernel methods in learning theory. Journal of computational and Applied Mathematics, 220(1-2):456–463, 2008. Xiaojin Zhu, Zoubin Ghahramani, and John D Lafferty. Semi-supervised learning using gaussian fields and harmonic functions. In Proceedings of the 20th International conference on Machine learning (ICML-03), pages 912–919, 2003. Wu Zongmin. Hermite-birkhoff interpolation of scattered data by radial basis functions. Approximation Theory and its Applications, 8(2):1–10, 1992.

Record · ID 411039 · SHA-256 225b078ac8200821
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.