arXiv CSopen access

Fast Learning Rates for Physics-Informed Kernel Methods

· arxiv_cs

arXiv CS · Papers · License: Open Access

machine learning, deep learning, neural networks

Fast Learning Rates for Physics-Informed Kernel Methods Luc Brogat-Motte1 , Joachim Bona-Pellissier2 , Giacomo Meanti2 , Lorenzo Rosasco1,2

arXiv:2609.18901v1 [stat.ML] 16 Sep 2026

1 2

SAIL Unit, Istituto Italiano di Tecnologia, Genoa, Italy

MaLGa Center, DIBRIS, Università degli Studi di Genova, Genoa, Italy [email protected], [email protected], [email protected], [email protected]

Abstract In physics-informed machine learning, a target function u∗ is learned from noisy value observations yi = u∗ (xi ) + εi , together with differential information, given either by noisy observations dj = (Du∗ )(zj ) + ξj or by a known physical constraint Du∗ = v. We consider the setting where D is a linear differential operator and analyze a physics-informed kernel estimator û combining n value observations and m differential observations. In this context, we ask how much can differential information improve predictions, and how does this improvement depend quantitatively on n, m, and D. We prove finite-sample bounds, supported by numerical simulations, revealing a two-regime structure for the prediction error. When m is limited, the rate depends jointly on n and m; when m exceeds a problem-dependent threshold, the rate saturates and matches the oracle rate obtained when the perfect constraint Dû = Du∗ is imposed. Examples are discussed for Sobolev spaces which are reproducing kernel Hilbert spaces and include partial Laplacian constraints on the torus and gradient observations on bounded domains. These examples illustrate the range of possible learning rate improvements – from the standard nonparametric n−1/4 to the parametric rate n−1/2 . Finally, we derive physically consistent rates in a stronger norm that jointly controls the errors in û and Dû.

Introduction

Supervised learning algorithms seek to estimate an unknown target function from noisy observations of its values. In many applications arising in physics, engineering, and scientific computing, however, one has access not only to value observations, but also to physical information involving derivatives or other differential quantities associated with the target function [1]. This additional information is often specified through functional constraints, such as gradients, Partial Differential Equation (PDE) residuals, or conservation laws [2–4]. Concrete examples include learning interatomic potentials, where energies and forces of atomic systems provide coupled value and gradient information through F = −∇E [5, 6]; 3D reconstruction problems, where surface normals give differential information about the underlying shape [7, 8]; or computational cardiology, where physical models involve differential constraints from fluid dynamics [9, 10]. In this paper we investigate physics-informed machine learning from the point of view of statistical learning theory. While there exists a comprehensive theory describing the statistical performance of 1

learning from value observations [11, 12], much less is known about how differential information affects learning rates and sample complexity. We are particularly interested in quantitative results of how differential information improves prediction: what is the improvement it can provide? how many differential observations are needed for this gain to appear? We study this question for target functions in a reproducing kernel Hilbert space (RKHS) H, which makes precise analyses feasible. We assume that there exists a target function u∗ ∈ H and that, in addition to noisy value observations yi = u∗ (xi ) + εi ,

i = 1, . . . , n,

also noisy differential quantities are available dj = (Du∗ )(zj ) + ξj ,

j = 1, . . . , m,

where D is a linear differential operator. Our contributions are as follows: • We derive finite-sample bounds that quantitatively describe how the prediction error depends on sample sizes n, m, differential operator D, and regularization parameters λ, γ. • We introduce a novel capacity assumption for physics-informed learning. This assumption decomposes the learning problem into a component that can be learned from differential data and a residual component not seen by D. It then quantifies how fast the first component can be learned from physics data and how much complexity remains in the residual one. We illustrate both quantities with two interpretable Sobolev RKHS examples, showing how they depend on the choice of D. • Under this assumption, we obtain learning rates by deriving the best choice of λ and γ. The rates exhibit two regimes: a differential-limited regime, where the error decreases jointly with n and m, and a saturation regime, where the rate matches the one obtained when the physical constraint is known exactly. In the Sobolev space examples, the saturated prediction rates range from n−1/4 to the parametric rate n−1/2 . Simulations further illustrate the predicted improvements and two-regime behavior. • Finally, we show that the same estimator is physically consistent beyond value prediction. We prove learning rates in a stronger norm that jointly controls the errors in û and Dû, ruling out accurate value prediction without learning the differential structure.

1.1

Related work

Settings. Learning from both function values and derivatives (or, more generally, differential operators) is classically known as Hermite-Birkhoff interpolation, whose study in a statistical setting dates back at least to the seminal work of Kimeldorf and Wahba [13]. This setting has recently resurfaced in machine learning under the names Sobolev training [14] and physics-informed learning [15]. Incorporating differential information, whether in discrete or continuous form, can serve various purposes. Differential operators can be used to regularize the learning procedure: penalties of R the form (Du(x))2 dx have been well studied for spline smoothing [16], inverse problems [17, 18], and machine learning [19]. Differential operators can also carry information about the data distribution.

This is the case in semi-supervised learning, where the operator D is learned from unlabeled data [20– 24]. Finally, differential constraints can come from the physical knowledge of the object of study itself, as is often the case in scientific applications [1]. Typically taking the form of a PDE, such constraints provide complementary data, either on the same domain as the value measurements [25] or on a different domain. A classical example of the latter is the one where Du∗ is known inside a domain Ω while u∗ is known only on the boundary ∂Ω, which is well studied in the PDE literature [26]. Methods. Deep learning-based algorithms to tackle problems with differential constraints include the Deep Ritz Method [27], neural operators [28, 29], Universal Differential Equations [4]. Closest to our approach are Physics-Informed Neural Networks (PINNs) [2], which minimize a loss function similar to ours. Kernel methods have also been considered to tackle the problems described above, from splines [13] to RBF collocation [30–33], semi-supervised learning [22] and, more recently, in scientific machine learning applications [34–38]. Our estimator belongs to this kernel-based family: it combines the usual kernel ridge loss on values with an empirical loss on differential observations Du(zj ). An orthogonal class of methods that is commonly used to estimate PDE solutions are finite elements/volumes [39, 25]. Theory. Many convergence results exist for learning PDE solutions with kernels [40, 33, 41, 42] or neural networks [43–46]. In the settings considered, data observations are only available on the boundary of the input domain (and not inside the domain). Such settings differ from ours and are closer to classical numerical methods for PDEs such as finite differences or finite element methods. Instead, we consider settings in which the value and differential observations lie in the same domain, and discuss the existing results for such settings, which are commonly of two types. A first type is concerned with physical consistency: showing that a physics-informed estimator û converges to the target u∗ in a physically consistent, stronger than L2 sense [47–49]. In Section 3.4, we provide a result of this nature with learning rates in a norm which controls both û and Dû. In particular, Doumèche et al. [47] shows, for a broad class of differential operators that PINNs are physically consistent. For kernel estimators, Shi et al. [48] and ul Abdeen et al. [49] show that learning with gradients allows to obtain convergence rates in the H 1 norm. A second type of result focuses on how the physics-informed penalty benefits the L2 convergence rates, which we show in Section 3.3. Previously, Fisher et al. [50] studied the impact of gradient information on L2 performance using a random feature model, and highlighted the existence of a regime where incorporating gradient information may hurt the prediction accuracy. Arnone et al. [51] investigated a regression estimator regularized with a 2-dimensional elliptic PDE, when u∗ belongs to the Sobolev space H 2 . Their estimator converges at least at the standard rate for H 2 functions (n−2/3 ), and can reach the faster n−4/5 rate when the PDE constraint is satisfied by the target. Doumèche et al. [37] reformulated empirical risk minimization with linear differential constraints as kernel regression in a physics-informed RKHS, yielding theoretical evidence that physical constraints can improve convergence rates, characterized through the effective dimension of the new kernel. Such improvements are illustrated for the setting of D = ∂x with homogeneous constraints (Du∗ = 0). Subsequent works [52, 53] focus on approximations and fast implementations, while preserving the underlying convergence rates. Overall, these theoretical results either study an ideal estimator for which continuous constraints are enforced (which typically is not computable in closed form) [51, 37], or tie the number of value and differential data points (m = n) [48–50]. In contrast, we provide results for any linear D which hold for arbitrary finite values of m and n, explicitly dealing with approximately enforced constraints. Finally, we must mention works which characterized the statistical performance of kernel regression [11, 54, 55], which we extend to handle physical constraints. They introduce the key 3

assumptions (source and capacity conditions) under which rates can be derived and the notion of effective dimension which controls the difficulty of the learning problem.

Physics-informed kernel regression

2.1

Problem setup

We now give a formal description of the learning problem. Let X be a measurable space and let u∗ : X → R denote the target function. We will consider linear differential operators D of order s ≥ 1, which take the form X D= cα ∂ α , |α|≤s

with coefficients cα ∈ R. We will assume (in (A2)) that Du∗ is well defined. Let ρ and ρD be distributions over X . We consider value and differential observations yi = u∗ (xi ) + εi , dj = (Du )(zj ) + ξj , ∗

xi ∼ ρ,

i.i.d.

i = 1, . . . , n,

i.i.d.

j = 1, . . . , m,

zj ∼ ρD ,

where (εi )i and (ξj )j are real-valued noise variables. In general, ρ and ρD can be different, reflecting the fact that value and differential information may be available on different regions of the domain.

2.2

Physics informed kernel estimator

Let (H, ⟨·, ·⟩H ) be an RKHS with associated kernel k : X × X → R. Given parameters λ > 0 and γ > 0, the regularized empirical risk is defined, for u ∈ H, as n m 2 2 1X 1 X b R(u) = u(xi ) − yi + γ Du(zj ) − dj + λ∥u∥2H , n i=1 m j=1

(1)

and the Physics Informed Kernel MethodS (PIKS) [38] estimator is the result of an optimization procedure in H b û := argmin R(u). u∈H

A small deviation from the estimator defined in [38] is the introduction of the scaling parameter γ. The parameter λ controls the regularity of the estimator, while γ controls the relative weight assigned to derivative observations. Setting γ = 0 recovers standard kernel ridge regression with value observations. b is a strictly convex quadratic functional on H, hence the minimizer û exists and Since λ > 0, R is unique. Moreover, the estimator admits a finite-dimensional representer form. In particular, û(·) =

n X

αi k(xi , ·) +

i=1

m X

βj D1 k(zj , ·),

j=1

where the coefficients α ∈ Rn and β ∈ Rm are the solutions to a linear system involving the value, mixed, and differential Gram matrices. The explicit system is given in Appendix A.1. 4

Main results

In this section, we present learning rates for PIKS and quantify the impact of derivative information. After listing the assumptions in Section 3.1, we state a general finite sample result in Section 3.2, from which, in Section 3.3, we derive learning rates, emphasizing the existence of two regimes. In Section 3.4 we study the rates in a stronger norm involving both u and Du.

3.1

Assumptions

Assumptions (A1) to (A5) are standard and ensure that derivative evaluations are well defined in the RKHS and that both value and derivative observations are generated by the same underlying target function. The main novel assumption is (A6), which quantifies how much of the function space is visible through derivative information. Assumption P (A1) (Linear differential operator). We consider linear differential operators D of the form D = |α|≤s cα ∂ α for some integer s ≥ 1 and coefficients cα ∈ R. Typical examples include directional derivatives, Laplacians, and more general linear partial differential operators. The order s determines the degree of smoothness required in the kernel such that derivative evaluations are well defined. Assumption (A2) (Kernel regularity). We assume that the kernel k is of class C s,s on X × X , meaning that all mixed derivatives ∂xα ∂xβ′ k(x, x′ ) with |α|, |β| ≤ s exist and are continuous. Under assumptions (A1) and (A2), the differential operator D can be applied to the kernel and gives rise to well-defined derivative representers in the RKHS. More precisely, let ϕ(x) := k(x, ·) be the RKHS feature map. Since derivative evaluations are continuous linear functionals on H [56, Corollary 4.36], and defining Dϕ(x) := Dk(x, ·), where D acts on the first argument of k, one has Dϕ(x) ∈ H and the reproducing property holds on both the kernel and its derivatives: u(x) = ⟨u, ϕ(x)⟩H ,

Du(z) = ⟨u, Dϕ(z)⟩H .

Assumption (A3) (Boundedness of features and outputs). There exist constants κV , κD , LV , LD > 0 such that, almost surely, ∥ϕ(x)∥H ≤ κV ,

∥Dϕ(z)∥H ≤ κD ,

|y| ≤ LV ,

|d| ≤ LD .

Such boundedness assumptions are standard in statistical analyses of kernel methods and are satisfied under assumptions (A1) and (A2) when X is compact. Assumption (A4) (Attainability). We assume that the regression function u∗ (x) := E[y | x] belongs to H. Equivalently, there exists w∗ ∈ H such that u∗ (x) = ⟨w∗ , ϕ(x)⟩H . This is a standard assumption for well-specified problems: it states that the target function belongs to the RKHS induced by the kernel, so it can be represented by the model class under consideration.

Assumption (A5) (Consistency of derivative observations). We assume that the derivative observations are unbiased measurements of the derivative of the regression function, in the sense that E[d | z] = Du∗ (z). This assumption states that, when Du∗ is known, it can be interpreted as having a physical constraint which is well-specified. It is also referred to as using a perfect measure ρD [48]. To quantify both the effective dimension of the hypothesis space and its alignment with value and derivative observations, we introduce the following covariance operators Σ := E[ϕ(x) ⊗ ϕ(x)],

ΣD := E[Dϕ(z) ⊗ Dϕ(z)],

where, for a, b ∈ H, the rank-one operator a ⊗ b : H → H is defined as (a ⊗ b)w := ⟨w, b⟩H a. For self-adjoint operators A, B, we write A ⪯ B when B − A is positive semi-definite. Assumption (A6) (Value-derivative capacity decomposition). Assume that there exist positive semi-definite operators Σ1 , Σ2 ⪰ 0 such that Σ = Σ1 + Σ2 , and that there exist constants c1 , c2 > 0 and α, r ∈ [0, 1] such that, for all λ ∈ (0, 1], Tr Σ1 (Σ1 + λI)−1 ≤ c1 λ−α , Tr (ΣD + λI)−1 Σ2 ≤ c2 λ−r .

(2)

This assumption decomposes the geometry of function values into directions that are invisible, or only weakly visible, to the operator D (captured by Σ1 ) and directions that are detectable through differential observations (captured by Σ2 ). The exponent α governs the effective dimension of the invisible component, while r measures how well the detectable component is separated from the null-space of D. In the limiting case Σ1 = Σ and Σ2 = 0, the assumption reduces to the usual effective-dimension condition of standard kernel regression. Together, these conditions formalize how differential information can reduce the statistical complexity of learning function values via capacity reduction, depending on interactions between kernel, differential operator, and data distribution.

3.2

Finite-sample bounds

Theorem 3.1 (High-probability finite-sample bounds for PIKS). Under assumptions (A1) to (A5), let λ > 0, γ > 0 and let δ ∈ (0, 1/2). Define Σγ := Σ + γΣD ,

B(λ, γ) := (Σγ + λI)−1/2 w∗ H ,

and the (value and derivative) capacity terms 1/2

NV (λ, γ) := (Σγ + λI)−1/2 Σ1/2 HS ,

ND (λ, γ) := (Σγ + λI)−1/2 ΣD

36 2n 36 2m log ≤ λ ≤ ∥Σ∥∞ and log ≤ λγ −1 ≤ ∥ΣD ∥∞ , (3) n δ m δ then there exists a constant c > 0, depending only on κV , κD , LV , LD and ∥u∗ ∥H , but not on λ, γ, n, m, δ, such that, with probability at least 1 − 2δ, h i1 4 1 1 ND (λ, γ) γ NV (λ, γ) 2 2 ∗ √ √ E (û(x) − u (x)) ≤ c log +γ √ + λB(λ, γ) . (4) + + δ n m λ n m 6

1.0

m < mcrit 0.8

m < mcrit(n)

Error

0.6

n 3 = n 2 m = n m = m

0.4

m ≥ mcrit

0.2

saturation regime

0.0 0.0

0.2

0.4

0.6

0.8

α = 0.2 α = 0.4 α = 0.6 α = 0.8

m ≥ mcrit(n)

1.0

Figure 1: (left) Rate-regime changes in the (α, r) plane. On the lower-right of each contour line the problem parameters are such that m is larger than the critical threshold. (right) Log-linear plot of the expected error as a function of n and α for fixed m, r. For smaller α a larger m/n ratio is required to be in the fast saturated regime. Sketch of the proof. The proof stems from decomposing the error into separate approximation and estimation terms, and controlling the stochastic parts through concentration inequalities for covariance operators in Hilbert spaces. The main novelty compared to the classical kernel ridge regression analysis lies in handling of the derivative observations. We introduce the covariance operator Σγ = Σ + γΣD , and decompose all deviation terms into value and derivative contributions while carefully tracking their dependence on γ, n, and m. Once these additional decompositions are established, the remainder of the argument closely parallels the standard KRR proof [11]. The bound in Theorem 3.1 separates the contributions of value observations, derivative observations, and regularization bias. The quantities NV (λ, γ) and ND (λ, γ) play the role of effective dimensions associated with the value and derivative components of the problem. They control the corresponding variance contributions, while the bias term λB(λ, γ) captures the approximation error induced by regularization. The contributions of NV (λ, γ) and ND (λ, γ) interact through the shared parameters λ and γ. Increasing γ places more weight on the derivative constraints, which reduces the effective complexity of the value component by shrinking the set of admissible functions in directions that are observable through D, but also amplifies the variance contribution coming from noisy derivative observations. Setting γ = 0 recovers the standard KRR bound [11]. Setting m = ∞ instead recovers the rates of Doumèche et al. [37].

3.3

Learning rate acceleration from differential information

We now instantiate the finite-sample bound under assumption (A6) and optimize over the regularization parameters λ and γ to precisely characterize the learning rates of the PIKS estimator. Corollary 3.1 (Learning rates for PIKS). Under assumptions (A1) to (A6), there exist choices of regularization parameters λ = λ(n, m) and γ = γ(n, m), given explicitly in the proof and which

depend polynomially on n, m, log n, and log m, such that, with high probability,  1 1−r n− 4 m− 8 , m ≲ mcrit , i1/2 h 2 ∗ E (û(x) − u (x)) ≲ 1  − 2(1+α) n , m ≳ mcrit , 2(1−α)

where mcrit ≍ n (1+α)(1−r) . Here and throughout, a(n, m) ≲ b(n, m) means that a(n, m) ≤ c b(n, m) for some c > 0 independent of n and m, up to logarithmic factors in n, m, and δ −1 (the high-probability parameter). The error bound reveals two regimes depending on the amount of differential samples m, relative to the number of function-value samples n as depicted also in Fig. 1: • Differential-limited regime (m ≲ mcrit ). With low m, the error decreases as both n and m increase. In particular as the number of differential observations m increases, the estimation error along directions captured by Σ2 decreases. The rate of decrease depends on how well the Σ2 directions are aligned with the differential covariance ΣD . The exponent r quantifies this effect, with smaller values of r leading to faster decay of the error with respect to m. • Saturation threshold (m ≍ mcrit ). Increasing m, a threshold is reached where differential and value contributions to the error are of the same order: 1

1−r

n− 4 m− 8

≍ n− 2(1+α) ,

2(1−α)

which yields mcrit ≍ n (1+α)(1−r) . The dependence of mcrit on α and r follows from this relation: smaller values of r lead to a slower decay of the differential term in m, and smaller values of α lead to a faster decay of the value term in n, so that a larger m is required to balance the two contributions. • Saturation regime (m ≳ mcrit ). Once m exceeds the threshold, the rate saturates: increasing the number of differential observations only reduces lower-order terms, but it no longer improves the leading error term. At this point, the error on the differential-accessible component is below the error driven by Σ1 , which captures directions that are either in the null-space of D, and therefore invisible to differential observations, or only weakly visible through them. These directions can be estimated only with value samples. In this regime, the leading rate matches the physical oracle rate obtained when the differential information is known exactly (see Appendix A.8).

3.4

Learning rates in physically consistent norms

In many physics-informed settings, respecting the differential structure is as important as minimizing prediction error. For instance, in learning interatomic potentials, predicting differential quantities such as atomic forces is essential for molecular dynamics simulations [57]. However, improved learning rates for function values do not by themselves ensure that û captures this structure. Indeed, convergence in value alone does not control derivatives: on any bounded domain, uε (x) = ε sin(ε−2 x) satisfies ∥uε ∥L2 → 0 as ε → 0, whereas for D = ∂x , Duε (x) = ε−1 cos(ε−2 x) and ∥Duε ∥L2 → ∞. We therefore study convergence in a stronger L2 -type norm that jointly controls function-value and differential errors. 8

Corollary 3.2 (Learning rates in the physically consistent norm). Under assumptions (A1) to (A6) there exist choices of regularization parameters λ = λ(n, m) and γ = γ(n, m), given explicitly in the proof with γ ≥ 1 such that, with high probability,  1 1−r n− 4 m− 8 + m−1/4 , m ≲ mcrit h i1/2 2 2 ∗ ∗ E (û(x) − u (x)) + (Dû(z) − Du (z)) ≲ 1 n− 2(1+α) + m−1/4 , m ≳ mcrit , 2(1−α)

where mcrit ≍ n (1+α)(1−r) . The stronger norm bounds rule out the possibility that the PIKS estimator has a good prediction accuracy on u∗ while failing to accurately learn Du∗ . This can be interpreted as a form of physical consistency: the learned function respects both the values and the differential structure of the target. Compared with the value-only rates in Corollary 3.1, this stronger guarantee incurs an additional term m−1/4 , which depends only on the number of derivative observations and therefore cannot be reduced by increasing n alone. This is a reasonable price to pay in settings where differential predictions are themselves of interest or are used in downstream tasks.

Examples

Here we illustrate the improved rates of Corollary 3.1 in concrete learning settings. We focus in Section 4.1 on Sobolev spaces on the torus where Laplacian information can provide large benefits, and in Section 4.2 on bounded domains, which are particularly relevant for their connection to PDEs.

4.1

Example 1: Sobolev spaces on the torus

We define the periodic Matérn kernel on Td = [0, 1]d with Fourier expansion X ′ k(x − x′ ) = µk e2πik·(x−x ) , µk ≍ (1 + |k|2 )−s , k∈Zd

for some s > d2 + 2, where ≍ denotes the equality P up to positive multiplicative constants. The associated RKHS H consists of functions f (x) = k∈Zd ck e2πik·x with norm ∥f ∥2H =

X |ck |2 k∈Zd

µk

≍

(1 + |k|2 )s |ck |2 ,

k∈Zd

where ck ∈ C are the Fourier coefficients. H is norm-equivalent to the Sobolev space H s (Td ). Proposition 4.1 (Capacity decomposition for SobolevPspaces). Let X = Td , let k be a periodic Matérn kernel of smoothness s > d2 + 2, and let D = i∈S ∂z2i be a partial Laplacian with S ⊆ {1, . . . , d}. Assume that value and differential samples are drawn uniformly on Td . Then the value covariance Σ admits a decomposition Σ = Σ1 + Σ2 satisfying assumption (A6), with exponents α=

d − |S| , 2s 9

d . 2s

Referring back to Corollary 3.1, this Sobolev example provides a functional-analytic interpretation of assumption (A6) and of the learning rates. For small m, the learning rate coincides with that of estimating a function in H s (Td ) from value observations, reflecting the full d-dimensional complexity of the data. When m exceeds the saturation threshold, all n function-value data points can be used to learn the nullspace of the partial Laplacian D, which depends only on d0 = d − |S| variables. The saturated rate therefore matches the minimax rate for Sobolev regression on Td0 , making explicit how the differential constraint reduces the problem dimensionality.

4.2

Example 2: Gradients on bounded domain

Consider now the case when D is the gradient operator ∇, and the RKHS H is a Sobolev space of smoothness s > d2 + 1 on domain X ⊂ Rd . We consider the reproducing kernel k associated with H to be a Matérn kernel of smoothness s [58]. Note that, while in Sections 2 and 3 the analysis is limited to the case of scalar-valued D, it still holds for vector-valued operators. Denoting Pk by Di each scalar component of D, we decompose the differential-data covariance as ΣD := ℓ=1 ΣDℓ , where ΣDℓ := E[Dℓ ϕ(z) ⊗ Dℓ ϕ(z)]. Provided this new ΣD satisfies assumption (A6) (which is unchanged), Corollary 3.1 still holds. We give more details on the extension to the vector-valued case in Appendix A.7. hP i d When D = ∇, we have ΣD = E i=1 ∂i ϕ(z) ⊗ ∂i ϕ(z) . The effect of learning with differential constraints is shown in the following proposition. Proposition 4.2. Let X ⊂ Rd be a bounded, connected and Lipschitz domain. Fix s > d2 + 1. Let H = H s (Ω) and let D = ∇ be the gradient operator. Then the value covariance operator Σ admits a decomposition Σ = Σ1 + Σ2 satisfying assumption (A6), with exponents α = 0,

d . 2(s − 1) 2(s−1)

The saturated rate from Corollary 3.1 is obtained as long as m ≥ mcrit = n (s−1)−d/2 ≥ n2 . Then, h i1/2 1 2 E (û(x) − u∗ (x)) ≲ n− 2(1+α) , 1

which proves that the gradient information allows to obtain the parametric rate n− 2 .

Numerical experiments

We show how the different rate-regimes of Corollary 3.1 look like in practice through simulations on two learning problems. We will show i) the error saturation effect when m increases and ii) the effect of decreasing the capacity of Σ1 (decreasing α) via the partial Laplacian. More details on the implementation are available in Appendix C. Saturation effect on Matérn data. We sample synthetic data on the 2d disk according to the following function, letting xsupp be a support point in the domain and p := ∥x − xsupp ∥ √ i √ h u∗ (x) = 1 − ∥x∥2 1 + 5p + (5/3)p2 exp − 5p .

RMSE kb u − u∗ k

0.1

10−3

n=4 n=8 n = 32 n = 128

0.2

0.0

D0 (|S| = 0) D1 (|S| = 1) D2 (|S| = 2) D3 (|S| = 3) D4 (|S| = 4)

10−4

101

102

101

103

m Figure 2: Theoretical and experimental rates for fixed n, increasing m showing the saturation effect with gradient information on Matérn data.

|S|

Figure 3: Test error rates increasing n and partial Laplacian dimensions |S|. Dashed lines show a best fit for the rates of Proposition 4.1. On the right are the inferred α exponents.

u∗ is in H s with s = 3, and we will use the corresponding Matérn (ν = 5/2) kernel. In Fig. 2 we show that, for fixed n and increasing m, the error rates switch between the two regimes of Corollary 3.1: an initial phase in which m < mcrit and the error decreases as m grows, followed by a phase in which m has grown beyond the threshold and the error has reached a saturation with respect to m in which additional data does not improve accuracy. Here we used D = ∇ and plotted empirical rates using PIKS next to the theoretical results from Section 4.2. Despite some differences – notably the experimental curves saturate at similar values of m – due to the unknown problem-dependent constants affecting the rates, the saturation effect is clear as well as the dependence of the initial rate on both n and m. Partial Laplacian. We now consider Sobolev spaces on the torus with partial Laplacian information. Unlike the example in Section 4.1, we approximate the infinite-dimensional Fourier construction by truncating to a finite number of frequencies. With X = T4 , we randomly sample 4 frequency √ vectors kℓ ∈ Z such that a single dimension is active (non-zero) at a time. Writing ϕℓ (x) = 2 cos(2π⟨kℓ , x⟩) for the feature-map, target function and kernel are defined as u∗ (x) = u0 +

F X

cℓ ϕℓ (x),

k(x, x′ ) = 1 +

ℓ=1

F X

µℓ ϕℓ (x)ϕℓ (x′ ),

ℓ=1

with F frequency coefficients cℓ and µℓ which depend on a smoothness parameter of the problem. kℓ are sampled such that u∗ decomposes into frequency blocks u∗ (x) = u0 +u1 (x)+u2 (x)+u3 (x)+u4 (x), each of which only Psdepends on a single dimension of x and cannot be observed unless the partial Laplacian Ds = i=1 ∂z2i includes that specific dimension. This creates a simple setting in which the incrementing the partial Laplacian order s, increases the amount of information which can be learned with the differential data. Moreover the Sobolev capacity decomposition from Proposition 4.1 remains valid after truncation, with constants independent of F . In Fig. 3 we observe the learning rate in n, as the type of differential data changes: from having access to function data only (D0 ) to having access to the full Laplacian (D4 ). The number of 11

differential points m is set as a constant factor of n. Best fits are obtained for the theoretical rates of Theorem 3.1 to obtain the exponent term of the error rate in n. Despite the high variance of the experiments, due to the noise ε, ξ and the random sampling of points xi , zj , it is evident that with each increase in the amount of differential information |S| the learning rate increases. From the fitted curves we can infer the values of α which – as predicted by Proposition 4.1 – have a linear relationship with the number of partial Laplacian dimensions (denoted by |S|), despite the mismatch caused by the finite-dimensional setting.

Conclusion and research directions

In this paper, we provide a precise characterization of how kernel regression benefits from learning with both function-value and differential observations. Our analysis quantifies how the prediction error depends jointly on the number of value samples, the number of differential samples, and the structure of the underlying differential operator. The resulting rates reveal two regimes: a differential-limited regime, in which the error decreases with both types of samples, and a saturation regime, in which the leading rate matches the physical oracle rate attainable with exact differential information. We also establish convergence in a stronger norm that controls both function-value and differential errors. Several questions remain open: are these rates minimax optimal? Can similar gains be obtained for nonlinear differential operators or misspecified physical constraints? Is it possible to preserve the statistical gains with approximate algorithms which are more efficient?

Acknowledgements This material is based upon work supported by the European Commission (Horizon Europe grant ELIAS 101120237), and the Ministry of Education, University and Research (FARE grant ML4IP R205T7J2KP).

References [1] 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. [2] 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. [3] 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. [4] 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. [5] Jörg Behler and Michele Parrinello. Generalized neural-network representation of high-dimensional potential-energy surfaces. Physical review letters, 98(14):146401, 2007.

[6] Albert P. Bartók, Mike C. Payne, Risi Kondor, and Gábor Csányi. Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons. Phys. Rev. Lett., 104, 2010. doi: 10.1103/ PhysRevLett.104.136403. [7] Amos Gropp, Lior Yariv, Niv Haim, Matan Atzmon, and Yaron Lipman. Implicit geometric regularization for learning shapes. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 3789–3799. PMLR, 13–18 Jul 2020. URL https://proceedings.mlr.press/v119/gropp20a.html. [8] Aarya Patel, Hamid Laga, and Ojaswa Sharma. Normal-guided detail-preserving neural implicit function for high-fidelity 3D surface reconstruction. Proceedings of the ACM on Computer Graphics and Interactive Techniques, 8(1), 2025. doi: https://doi.org/10.1145/3728293. [9] Alfio Quarteroni, Paola Gervasio, and Francesco Regazzoni. Combining physics-based and data-driven models: advancing the frontiers of research with scientific machine learning. Mathematical Models and Methods in Applied Sciences, 35(4):905–1071, 2025. doi: 10.1142/S0218202525500125. [10] Tommaso Botarelli, Marco Fanfani, Paolo Nesi, and Lorenzo Pinelli. Using physics-informed neural networks for solving navier-stokes equations in fluid dynamic complex scenarios. Engineering Applications of Artificial Intelligence, 148, 2025. doi: https://doi.org/10.1016/j.engappai.2025.110347. [11] Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squares algorithm. Foundations of Computational mathematics, 7(3):331–368, 2007. [12] László Györfi, Michael Kohler, Adam Krzyżak, and Harro Walk. A distribution-free theory of nonparametric regression. Springer, 2002. [13] George Kimeldorf and Grace Wahba. Some results on Tchebycheffian spline functions. Journal of mathematical analysis and applications, 33(1):82–95, 1971. [14] 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. [15] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations, 2017. [16] Grace Wahba. Spline models for observational data. SIAM, 1990. [17] Martin Hanke. Regularization with differential operators: an iterative approach. Numerical functional analysis and optimization, 13(5-6):523–540, 1992. [18] Heinz Werner Engl, Martin Hanke, and Andreas Neubauer. Regularization of inverse problems, volume 375. Springer Science & Business Media, 1996. [19] 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. [20] 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. [21] Dengyong Zhou and Bernhard Schölkopf. Regularization on discrete spaces. In Joint Pattern Recognition Symposium, pages 361–368. Springer, 2005. [22] 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.

[23] Dejan Slepcev and Matthew Thorpe. Analysis of p-Laplacian regularization in semisupervised learning. SIAM Journal on Mathematical Analysis, 51(3):2085–2120, 2019. [24] 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. [25] Laura M Sangalli. Spatial regression with partial differential equation regularisation. International Statistical Review, 89(3):505–531, 2021. [26] Lawrence C Evans. Partial differential equations, volume 19. American mathematical society, 2010. [27] Bing Yu et al. The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems. Communications in Mathematics and Statistics, 6(1), 2018. [28] 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), 2021. [29] 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. [30] 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. [31] Wu Zongmin. Hermite-Birkhoff interpolation of scattered data by radial basis functions. Approximation Theory and its Applications, 8(2):1–10, 1992. [32] Gregory E Fasshauer. Solving partial differential equations by collocation with radial basis functions. In Proceedings of Chamonix, volume 1997, pages 1–8, 1996. [33] Holger Wendland. Scattered data approximation, volume 17. Cambridge university press, 2004. [34] Houman Owhadi. Bayesian numerical homogenization. Multiscale Modeling & Simulation, 13(3): 812–828, 2015. [35] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Machine learning of linear differential equations using Gaussian processes. Journal of Computational Physics, 348:683–693, 2017. [36] 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. [37] Nathan Doumèche, Francis Bach, Gérard Biau, and Claire Boyer. Physics-informed machine learning as a kernel method. In The Thirty Seventh Annual Conference on Learning Theory, pages 1399–1450. PMLR, 2024. [38] Joachim Bona-Pellissier, Giacomo Meanti, Matteo Santacesaria, and Lorenzo Rosasco. Piks: Universal physics-informed kernel methods. arXiv preprint arXiv:2607.27062, 2026. [39] 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.

[40] Carsten Franke and Robert Schaback. Solving partial differential equations by collocation using radial basis functions. Applied Mathematics and computation, 93(1):73–82, 1998. [41] Robert Schaback and Holger Wendland. Kernel techniques: from machine learning to meshless methods. Acta numerica, 15, 2006. [42] 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, 2025. [43] Yeonjong Shin, Jerome 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. [44] 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. [45] Yeonjong Shin, Zhongqiang Zhang, and George Em Karniadakis. Error estimates of residual minimization using neural networks for linear PDEs. Journal of Machine Learning for Modeling and Computing, 4 (4), 2023. [46] 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. [47] Nathan Doumèche, Gérard Biau, and Claire Boyer. On the convergence of PINNs. Bernoulli, 31(3): 2127 – 2151, 2025. doi: 10.3150/24-BEJ1799. [48] Lei Shi, Xin Guo, and Ding-Xuan Zhou. Hermite learning with gradient data. Journal of computational and applied mathematics, 233(11):3046–3059, 2010. [49] Zain ul Abdeen, Ruoxi Jia, Vassilis Kekatos, and Ming Jin. A theoretical analysis of using gradient data for Sobolev training in RKHS. IFAC-PapersOnLine, 56(2), 2023. doi: https://doi.org/10.1016/j. ifacol.2023.10.1491. 22nd IFAC World Congress. [50] Katharine E Fisher, Matthew TC Li, Youssef Marzouk, and Timo Schorlepp. Precise asymptotic analysis of Sobolev training for random feature models. arXiv preprint arXiv:2511.03050, 2025. [51] 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. [52] Nathan Doumèche, Francis Bach, Gérard Biau, and Claire Boyer. Physics-informed kernel learning. Journal of Machine Learning Research, 26(124):1–39, 2025. [53] 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, 2025. [54] Ingo Steinwart, Don R Hush, Clint Scovel, et al. Optimal rates for regularized least squares regression. In COLT, pages 79–93, 2009. [55] Simon Fischer and Ingo Steinwart. Sobolev norm learning rates for regularized least-squares algorithms. Journal of Machine Learning Research, 21(205):1–38, 2020. [56] Ingo Steinwart and Andreas Christmann. Support vector machines. Springer Science & Business Media, 2008.

[57] Frank Noé, Alexandre Tkatchenko, Klaus-Robert Müller, and Cecilia Clementi. Machine learning for molecular simulation. Annual review of physical chemistry, 71(1), 2020. [58] Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006. [59] Alessandro Rudi, Guillermo D Canas, and Lorenzo Rosasco. On the sample complexity of subspace learning. Advances in Neural Information Processing Systems, 26, 2013. [60] Luc Brogat-Motte, Riccardo Bonalli, and Alessandro Rudi. Learning controlled stochastic differential equations. arXiv preprint arXiv:2411.01982, 2024.

Proofs

This section provides the proofs of the main results. We first derive the closed-form expressions of the empirical and population minimizers in Section A.1, then establish an error decomposition in Section A.2, derive high-probability bounds for its constituent terms in Section A.3, and finally combine these results to obtain finite-sample guarantees and optimized learning rates in Sections A.4 and A.5.

A.1

Closed forms

We derive here the closed-form expressions of the empirical minimizer of the regularized empirical b introduced in Section 2, risk R n m 2 2 1X 1 X b R(u) = u(xi ) − yi + γ Du(zj ) − dj + λ∥u∥2H , n i=1 m j=1

and of the minimizer of its population counterpart Rλ,γ (u) = E (u(x) − y)2 + γ E (Du(z) − d)2 + λ∥u∥2H , where the expectations are taken over x ∼ ρ, z ∼ ρD , and the corresponding observation noises. The empirical minimizer is given in operator form in Lemma A.1 and in finite-dimensional form in Lemma A.2. The population minimizer is given in Lemma A.3. Lemma A.1 (Empirical solution). Let λ > 0 and γ > 0. Define n X b := 1 Σ ϕ(xi ) ⊗ ϕ(xi ), n i=1 n 1X yi ϕ(xi ), Vb := n i=1

m X b D := 1 Σ Dϕ(zj ) ⊗ Dϕ(zj ), m j=1 m 1 X VbD := dj Dϕ(zj ). m j=1

b has a unique minimizer û ∈ H of the form Then R û(x) = ⟨w, b ϕ(x)⟩H ,

b + γΣ b D + λI −1 Vb + γ VbD . w b= Σ

b b Proof. Using the representation u(x) = ⟨w, ϕ(x)⟩H , Du(z) = ⟨w, Dϕ(z)⟩H , we write R(u) = R(w) as n m 2 2 1X 1 X b R(w) = ⟨w, ϕ(xi )⟩H − yi + γ ⟨w, Dϕ(zj )⟩H − dj + λ⟨w, w⟩H . n i=1 m j=1 Expanding the least-squares terms in H, we obtain for the value part n n 2 1X 1 X ⟨w, ϕ(xi )⟩ − yi = ⟨w, ϕ(xi )⟩2 − 2yi ⟨w, ϕ(xi )⟩ + yi2 n i=1 n i=1 * + * n + n X X 1 1 = w, n ϕ(xi ) ⊗ ϕ(xi ) w −2 n yi ϕ(xi ), w + const i=1

b H − 2⟨Vb , w⟩H + const, = ⟨w, Σw⟩ 17

i=1

and for the derivative part m

2 1 X b D w⟩H − 2⟨VbD , w⟩H + const . ⟨w, Dϕ(zj )⟩ − dj = γ ⟨w, Σ m j=1

Adding the regularization term λ⟨w, w⟩H , we get b + γΣ b D + λI)w − 2 Vb + γ VbD , w + const. b R(w) = w, (Σ H H b and Σ b D are positive self-adjoint and λ > 0, the operator Σ b + γΣ b D + λI is strictly positive. Since Σ b Thus R(w) is a strictly convex quadratic functional of w, whose unique minimizer w b is characterized by the normal equation b + γΣ b D + λI) w (Σ b = Vb + γ VbD , which yields

b + γΣ b D + λI −1 Vb + γ VbD . w b= Σ

The corresponding minimizer in function space is û(x) = ⟨w, b ϕ(x)⟩H . Lemma A.2 (Finite-dimensional form of the PIKS estimator). Let λ > 0 and γ > 0. Define the Gram matrices (KXX )ii′ = k(xi , xi′ ), (KXZ )ij = D1 k(zj , xi ), (KZX )ji = D2 k(xi , zj ),

D (KZZ )jj ′ = D1 D2 k(zj , zj ′ ).

Then the PIKS estimator can be computed as r n m 1 X γ X û(·) = √ αi k(xi , ·) + βj D1 k(zj , ·), m j=1 n i=1

θ=

α , β

where θ ∈ Rn+m solves "

1 KXX pnγ nm KZX

p γ

nm KXZ

γ D m KZZ

! + λI

√ y/ √ n θ= √ . γ d/ m

n m Proof. Let (ei )ni=1 and (fj )m j=1 denote the canonical bases of R and R , and define

Φei = k(xi , ·),

DΦfj = D1 k(zj , ·).

Set Γγ = Then The optimality condition gives

√1 Φ n

γ m DΦ ,

√ y/ √ n bγ = √ . γ d/ m

b R(u) = ∥Γ∗γ u − bγ ∥22 + λ∥u∥2H . (Γγ Γ∗γ + λI)û = Γγ bγ . 18

Using

(Γγ Γ∗γ + λI)−1 Γγ = Γγ (Γ∗γ Γγ + λI)−1 ,

we get

û = Γγ θ,

Finally,

θ = (Γ∗γ Γγ + λI)−1 bγ .

" Γ∗γ Γγ =

p γ

1 KXX pn γ nm KZX

nm KXZ γ D m KZZ

# ,

which gives the claimed system. Moreover, since 1 û = Γγ θ = √ Φα + n we obtain

1 X αi k(xi , ·) + û(·) = √ n i=1

γ DΦβ, m m

γ X βj D1 k(zj , ·). m j=1

Lemma A.3 (Population solution). Let x ∼ ρ and z ∼ ρD , and let the expectations involving y and d be taken with respect to the joint laws induced by the observation models. Define Σ := Ex∼ρ ϕ(x) ⊗ ϕ(x) , ΣD := Ez∼ρD Dϕ(z) ⊗ Dϕ(z) , V := E y ϕ(x) , VD := E d Dϕ(z) . Then Rλ,γ has a unique minimizer uλ,γ ∈ H of the form uλ,γ (x) = ⟨wλ,γ , ϕ(x)⟩H ,

wλ,γ = Σ + γΣD + λI

−1

(V + γVD ).

Proof. The derivations are the same as in the empirical case, (Lemma A.1) with empirical averages b Σ b D , Vb , VbD replaced by Σ, ΣD , V, VD , and with the same quadratic replaced by expectations and Σ, expansion of Rλ,γ (w).

A.2

Error decomposition

We next establish a standard error decomposition adapted to our setting. The idea is to separate the error into a bias term and an estimation term, controlled by the deviations of the empirical covariance operators from their population counterparts. Lemma A.4 (Error decomposition). Let Σγ := Σ + γΣD ,

b γ := Σ b + γΣ bD, Σ

Vγ := V + γVD ,

Vbγ := Vb + γ VbD ,

and define wλ,γ := (Σγ + λI)−1 Vγ ,

b γ + λI)−1 Vbγ , w b := (Σ

u∗ (x) = ⟨w∗ , ϕ(x)⟩H .

Assume Assumption (A4) holds. Define the quantities b γ + λI)−1/2 ∆Σ,1 := ∥(Σγ + λI)1/2 (Σ bγ) ∆Σ,2 := (Σγ + λI)−1/2 (Σγ − Σ ∆V := (Σγ + λI)−1/2 (Vbγ − Vγ ) B(λ, γ) := (Σγ + λI)−1/2 w∗ . Then ∥Σ1/2 (w b − w∗ )∥ ≤ ∆2Σ,1 (∆V + w∗ ∆Σ,2 ) + λB(λ, γ). Proof. We first split the error into an estimation and a bias term: b − w∗ )∥ ≤ ∥Σ1/2 (w b − wλ,γ )∥ + ∥Σ1/2 (wλ,γ − w∗ )∥. ∥Σ1/2 (w We now bound each term. Bias term. Under the Assumption (A4) we have V = Σw∗ , hence

VD = ΣD w∗ ,

Vγ = V + γVD = (Σ + γΣD )w∗ = Σγ w∗ .

Therefore

wλ,γ − w∗ = (Σγ + λI)−1 Vγ − w∗ = (Σγ + λI)−1 Σγ w∗ − w∗ = (Σγ + λI)−1 Σγ − I w∗ = −λ(Σγ + λI)−1 w∗ .

As a consequence, ∥Σ1/2 (wλ,γ − w∗ )∥ = λ Σ1/2 (Σγ + λI)−1 w∗ ≤ λ (Σγ + λI)−1/2 w∗ , where we used ∥Σ1/2 (Σγ + λI)−1/2 ∥ ≤ 1 since Σ ≼ Σγ + λI. Estimation term. We have h i b γ + λI)−1 (Vbγ − Vγ ) + (Σ b γ + λI)−1 − (Σγ + λI)−1 Vγ , w b − wλ,γ = (Σ b γ + λI, B = Σγ + λI, and, using A−1 − B −1 = A−1 (B − A)B −1 with A = Σ b γ + λI)−1 − (Σγ + λI)−1 = (Σ b γ + λI)−1 (Σγ − Σ b γ )(Σγ + λI)−1 . (Σ Multiplying by Σ1/2 and inserting (Σγ + λI)±1/2 , we obtain b γ + λI)−1/2 (Σ b γ + λI)−1/2 (Σγ + λI)1/2 ∥Σ1/2 (w b − wλ,γ )∥ ≤ Σ1/2 (Σ h × (Σγ + λI)−1/2 (Vbγ − Vγ ) bγ) + (Σγ + λI)−1/2 (Σγ − Σ 20

(Σγ + λI)−1 Vγ

i .

Then, b γ + λI)−1/2 = Σ1/2 (Σγ + λI)−1/2 (Σγ + λI)1/2 (Σ b γ + λI)−1/2 Σ1/2 (Σ b γ + λI)−1/2 ≤ Σ1/2 (Σγ + λI)−1/2 (Σγ + λI)1/2 (Σ b γ + λI)−1/2 , ≤ (Σγ + λI)1/2 (Σ where we used again that ∥Σ1/2 (Σγ + λI)−1/2 ∥ ≤ 1. Moreover, using Vγ = Σγ w∗ , we have (Σγ + λI)−1 Vγ = (Σγ + λI)−1 Σγ w∗ ≤ w∗ . Conclusion. We define b γ + λI)−1/2 ∆Σ,1 := (Σγ + λI)1/2 (Σ bγ) ∆Σ,2 := (Σγ + λI)−1/2 (Σγ − Σ ∆V := (Σγ + λI)−1/2 (Vbγ − Vγ ) B(λ, γ) := (Σγ + λI)−1/2 w∗ . Collecting the terms, we obtain ∥Σ1/2 (w b − wλ,γ )∥ ≤ ∆2Σ,1 (∆V + w∗ ∆Σ,2 ), and

∥Σ1/2 (wλ,γ − w∗ )∥ ≤ λB(λ, γ),

which together yields ∥Σ1/2 (w b − w∗ )∥ ≤ ∆2Σ,1 (∆V + w∗ ∆Σ,2 ) + λB(λ, γ).

A.3

High-probability bounds for the error decomposition terms ∆Σ,1 , ∆Σ,2 , and ∆V

We next control the random quantities appearing in the error decomposition of Lemma A.4. Compared with the standard proof for kernel ridge regression, the only additional difficulty is the presence of the derivative term and the need to keep track of its dependence on γ. This leads us to decompose both covariance and moment deviations into value and derivative contributions. The remainder of the argument then follows the standard kernel ridge regression proof strategy: we derive high-probability bounds for ∆V and ∆Σ,2 , and control the multiplicative factor ∆Σ,1 through the auxiliary normalized covariance deviation b γ )(Σγ + λI)−1/2 . ∆Σ,3 := (Σγ + λI)−1/2 (Σγ − Σ The key point is that when ∆Σ,3 < 1, the empirical covariance operator remains well conditioned relative to Σγ + λI, which yields a bound on ∆Σ,1 .

Lemma A.5 (High-probability bound on ∆Σ,3 ). We define b γ )(Σγ + λI)−1/2 . ∆Σ,3 := (Σγ + λI)−1/2 (Σγ − Σ Let δ ∈ (0, 1) and λ > 0 satisfy 36 n log ≤ λ ≤ ∥Σ∥∞ , n δ Then, with probability at least 1 − 2δ,

36 m log ≤ λγ −1 ≤ ∥ΣD ∥∞ . m δ

∆Σ,3 ≤ 21 .

Proof. We first relate ∆Σ,3 to the value and derivative parts separately. Using the decomposition b γ = (Σ − Σ) b + γ(ΣD − Σ b D ), Σγ − Σ and the triangle inequality, we obtain b γ )(Σγ + λI)−1/2 ∆Σ,3 = (Σγ + λI)−1/2 (Σγ − Σ −1/2 b b D )(Σγ + λI)−1/2 . ≤ (Σγ + λI)−1/2 (Σ − Σ)(Σ + γ (Σγ + λI)−1/2 (ΣD − Σ γ + λI)

For the first term, write −1/2 b b (Σγ + λI)−1/2 (Σ − Σ)(Σ = A (Σ + λI)−1/2 (Σ − Σ)(Σ + λI)−1/2 A∗ , γ + λI)

where

A := (Σγ + λI)−1/2 (Σ + λI)1/2 .

From Σ ≼ Σγ we have (Σγ + λI)−1 ≼ (Σ + λI)−1 , hence ∥A∥ ≤ 1. Therefore −1/2 b b (Σγ + λI)−1/2 (Σ − Σ)(Σ ≤ (Σ + λI)−1/2 (Σ − Σ)(Σ + λI)−1/2 . γ + λI)

For the second term we argue analogously, but comparing with γΣD . Using γΣD ≼ Σγ we obtain b D )(Σγ + λI)−1/2 ≤ (γΣD + λI)−1/2 (ΣD − Σ b D )(γΣD + λI)−1/2 , (Σγ + λI)−1/2 (ΣD − Σ and hence ∆Σ,3 ≤

b (Σ + λI)−1/2 (Σ − Σ)(Σ + λI)−1/2

b D )(γΣD + λI)−1/2 . + γ (γΣD + λI)−1/2 (ΣD − Σ

We now control each term probabilistically using Lemma 3.6 in [59]. Although the lemma yields a bound of 1/2, this constant can be reduced by rescaling the regularization parameter and applying the same proof. In particular, applying the proof of Lemma 3.6 with target bound 1/4, and using the admissibility conditions 36 n log ≤ λ ≤ ∥Σ∥∞ , n δ

36 m log ≤ λ ≤ γ∥ΣD ∥∞ , m δ

we obtain that, with probability at least 1 − δ, b (Σ + λI)−1/2 (Σ − Σ)(Σ + λI)−1/2 22

< 41 ,

and, again with probability at least 1 − δ, b D )(γΣD + λI)−1/2 γ (γΣD + λI)−1/2 (ΣD − Σ

< 41 .

Here the improvement from 12 to 14 comes from running the same concentration proof with the smaller target threshold 14 , at the cost of replacing the constant 9 in the lower admissibility condition by 36. By a union bound, both events hold simultaneously with probability at least 1 − 2δ, and on this event we have ∆Σ,3 ≤ 41 + 14 = 12 . This concludes the proof.

Lemma A.6 (Bound on ∆Σ,1 ). Recall b γ + λI)−1/2 . ∆Σ,1 := (Σγ + λI)1/2 (Σ We have ∆Σ,1 ≤ (1 − ∆Σ,3 )−1/2 . Proof. We have Define

b γ + λI)−1 (Σγ + λI)1/2 . ∆2Σ,1 = (Σγ + λI)1/2 (Σ b γ )(Σγ + λI)−1/2 . b := (Σγ + λI)−1/2 (Σγ − Σ B

By definition of ∆Σ,3 , b = ∆Σ,3 . ∥B∥ We can rewrite 1/2 b γ + λI = Σγ + λI − (Σγ − Σ b γ ) = (Σγ + λI)1/2 (I − B)(Σ b Σ , γ + λI)

b γ + λI)−1 = (Σγ + λI)−1/2 (I − B) b −1 (Σγ + λI)−1/2 . (Σ

Plugging this into the expression for ∆2Σ,1 gives b −1 ∥. ∆2Σ,1 = ∥(I − B) b = ∆Σ,3 < 1, then the spectrum of B b lies in [−∆Σ,3 , ∆Σ,3 ], hence the spectrum of (I − B) b −1 If ∥B∥ −1 is contained in {(1 − µ) : |µ| ≤ ∆Σ,3 }, and therefore b −1 ∥ ≤ (1 − ∥B∥) b −1 = (1 − ∆Σ,3 )−1 . ∆2Σ,1 = ∥(I − B)

Lemma A.7 (High-probability bound on ∆V ). Recall ∆V := (Σγ + λI)−1/2 (Vbγ − Vγ ) . Then ∆V ≤

+ γ (Σγ + λI)−1/2 (VbD − VD ) .

(Σγ + λI)−1/2 (Vb − V )

For any δ ∈ (0, 1), with probability at least 1 − 2δ, √ log(2/δ) ≤ 2 βV + 2 σV n

log(2/δ) n ! r log(2/δ) log(2/δ) √ + γ 2 βD + 2 σD , m m

∆V

where the explicit constants are 2

βV := 2 λ−1/2 LV κV ,

σV2 := L2V κ2V (Σγ + λI)−1/2 Σ1/2 HS ,

βD := 2 λ−1/2 LD κD ,

2 σD := L2D κ2D (Σγ + λI)−1/2 ΣD

Proof. The decomposition

1/2 2 . HS

Vbγ − Vγ = (Vb − V ) + γ(VbD − VD )

immediately gives ∆V ≤ (Σγ + λI)−1/2 (Vb − V ) + γ (Σγ + λI)−1/2 (VbD − VD ) . Value part.

Write

(Σγ + λI)

−1/2

n 1X b M (xi , yi ), (V − V ) = n i=1

M (x, y) := (Σγ + λI)−1/2 (yϕ(x) − V ).

Since |y| ≤ LV and ∥ϕ(x)∥ ≤ κV , ∥yϕ(x) − V ∥ ≤ LV κV + ∥V ∥ ≤ 2LV κV , hence

∥M (x, y)∥ ≤ ∥(Σγ + λI)−1/2 ∥ 2LV κV ≤ 2λ−1/2 LV κV =: βV .

The second moment satisfies E∥M (x, y)∥2 ≤ L2V κ2V ∥(Σγ + λI)−1/2 Σ1/2 ∥2HS =: σV2 . Applying Hilbert-space Bernstein (Prop. 7.15 of [60]) gives, with probability at least 1 − δ, r log(2/δ) √ log(2/δ) −1/2 b (Σγ + λI) (V − V ) ≤ 2βV + 2 σV . n n

Derivative part.

Exactly the same argument with MD (z, d) := (Σγ + λI)−1/2 (d Dϕ(z) − VD )

yields

∥MD (z, d)∥ ≤ 2λ−1/2 LD κD =: βD ,

and

1/2

2 σD = L2D κ2D ∥(Σγ + λI)−1/2 ΣD ∥2HS .

Bernstein again gives, with probability at least 1 − δ, (Σγ + λI)

−1/2

Conclusion.

log(2/δ) (VbD − VD ) ≤ 2 βD + m

√

r 2 σD

log(2/δ) . m

A union bound yields the stated bound with probability at least 1 − 2δ.

Lemma A.8 (High-probability bound on ∆Σ,2 ). Recall bγ) , ∆Σ,2 := (Σγ + λI)−1/2 (Σγ − Σ

Σγ = Σ + γΣD ,

bγ = Σ b + γΣ bD. Σ

Then ∆Σ,2 ≤

b (Σγ + λI)−1/2 (Σ − Σ)

bD) . + γ (Σγ + λI)−1/2 (ΣD − Σ

Moreover, for any δ ∈ (0, 1), with probability at least 1 − 2δ, ! r r √ √ log(2/δ) log(2/δ) log(2/δ) log(2/δ) + 2 σΣ + γ 2 βΣD + 2 σΣD , ∆Σ,2 ≤ 2 βΣ n n m m where βΣ := λ−1/2 2κ2V , βΣD := λ−1/2 2κ2D , Proof. The decomposition

2 σΣ := 4κ2V (Σγ + λI)−1/2 Σ1/2 HS , 1/2 2 . HS

2 σΣ := 4κ2D (Σγ + λI)−1/2 ΣD D

b γ = (Σ − Σ) b + γ(ΣD − Σ bD) Σγ − Σ

and the triangle inequality give bγ) ∆Σ,2 = (Σγ + λI)−1/2 (Σγ − Σ b + γ (Σγ + λI)−1/2 (ΣD − Σ bD) . ≤ (Σγ + λI)−1/2 (Σ − Σ) We now bound these two terms separately, with the same Hilbert-space Bernstein argument as in Lemma A.7, but now in the Hilbert space of Hilbert-Schmidt operators and with yϕ(x) replaced by ϕ(x) ⊗ ϕ(x).

Value part.

Recall Σ = E[ϕ(x) ⊗ ϕ(x)],

Define

U (x) := Σ − ϕ(x) ⊗ ϕ(x),

n X b= 1 Σ ϕ(xi ) ⊗ ϕ(xi ). n i=1

MΣ (x) := (Σγ + λI)−1/2 U (x).

Then E[U (x)] = 0, and

b = (Σγ + λI)−1/2 (Σ − Σ)

1X MΣ (xi ). n i=1

We view MΣ (x) as a random element of the Hilbert space of Hilbert-Schmidt operators. We first bound its norm. Using the boundedness ∥ϕ(x)∥ ≤ κV we have ∥ϕ(x) ⊗ ϕ(x)∥HS = ∥ϕ(x)∥2 ≤ κ2V , Thus

∥Σ∥HS ≤ E∥ϕ(x) ⊗ ϕ(x)∥HS ≤ κ2V .

∥U (x)∥HS ≤ ∥ϕ(x) ⊗ ϕ(x)∥HS + ∥Σ∥HS ≤ 2κ2V ,

and since ∥(Σγ + λI)−1/2 ∥ ≤ λ−1/2 , ∥MΣ (x)∥HS ≤ λ−1/2 2κ2V = βΣ . Next, we bound the second moment. Using that Σγ + λI is self-adjoint and positive, we have 2 ∥MΣ (x)∥2HS = (Σγ + λI)−1/2 U (x) HS = Tr U (x)(Σγ + λI)−1 U (x) . Take expectation and use linearity of the trace: E∥MΣ (x)∥2HS = Tr (Σγ + λI)−1 E U (x)2 . We now bound E[U (x)2 ] in terms of Σ. First note that 2 2 U (x)2 = Σ − ϕ(x) ⊗ ϕ(x) ≼ 2Σ2 + 2 ϕ(x) ⊗ ϕ(x) , by the elementary operator inequality (A − B)2 = A2 + B 2 − AB − BA ≼ 2A2 + 2B 2 for self-adjoint A, B. Moreover, 2 ϕ(x) ⊗ ϕ(x) = ∥ϕ(x)∥2 ϕ(x) ⊗ ϕ(x) ≼ κ2V ϕ(x) ⊗ ϕ(x), and Σ2 ≼ ∥Σ∥ Σ ≼ κ2V Σ, since ∥Σ∥ ≤ κ2V for a bounded feature map. Combining these, U (x)2 ≼ 2κ2V ϕ(x) ⊗ ϕ(x) + 2κ2V Σ. Taking expectations yields E[U (x)2 ] ≼ 2κ2V E[ϕ(x) ⊗ ϕ(x)] + 2κ2V Σ = 4κ2V Σ.

Plugging this into the expression for the second moment, we get E∥MΣ (x)∥2HS ≤ Tr (Σγ + λI)−1 4κ2V Σ = 4κ2V Tr (Σγ + λI)−1/2 Σ(Σγ + λI)−1/2 2

= 4κ2V (Σγ + λI)−1/2 Σ1/2 HS . Thus we set

2 σΣ := 4κ2V (Σγ + λI)−1/2 Σ1/2 HS

and obtain

2 E∥MΣ (x)∥2HS ≤ σΣ .

We can now apply the Hilbert-space Bernstein inequality (Proposition 7.15 in [60]) in the Hilbert Pn space of Hilbert-Schmidt operators to the empirical average n1 i=1 MΣ (xi ). This gives: for any δ ∈ (0, 1), with probability at least 1 − δ, r log(2/δ) √ log(2/δ) −1/2 b (Σ − Σ) HS ≤ 2 βΣ (Σγ + λI) + 2 σΣ . n n Since ∥T ∥ ≤ ∥T ∥HS for every Hilbert-Schmidt operator T , the same bound holds for b . (Σγ + λI)−1/2 (Σ − Σ) Derivative part.

The argument is identical with m X bD = 1 Dϕ(zj ) ⊗ Dϕ(zj ). Σ m j=1

ΣD = E[Dϕ(z) ⊗ Dϕ(z)], Define

UD (z) := ΣD − Dϕ(z) ⊗ Dϕ(z),

MΣD (z) := (Σγ + λI)−1/2 UD (z).

Then

bD) = (Σγ + λI)−1/2 (ΣD − Σ

1 X MΣD (zj ), m j=1

with E[MΣD (z)] = 0. Using ∥Dϕ(z)∥ ≤ κD we obtain ∥UD (z)∥HS ≤ 2κ2D ,

∥MΣD (z)∥HS ≤ λ−1/2 2κ2D = βΣD ,

and, by the same operator inequalities as above, 1/2 2 2 =: σΣ . D HS

E∥MΣD (z)∥2HS ≤ 4κ2D (Σγ + λI)−1/2 ΣD

Another application of the Hilbert-space Bernstein inequality in the Hilbert space of Hilbert-Schmidt operators yields, with probability at least 1 − δ, r log(2/δ) √ log(2/δ) −1/2 b (Σγ + λI) (ΣD − ΣD ) HS ≤ 2 βΣD + 2 σΣD . m m Again the same bound holds for bD) . (Σγ + λI)−1/2 (ΣD − Σ 27

Conclusion.

We have shown that, with probability at least 1 − δ, (Σγ + λI)

−1/2

√ b ≤ 2 βΣ log(2/δ) + 2 σΣ (Σ − Σ) n

log(2/δ) , n

and, with probability at least 1 − δ, (Σγ + λI)

−1/2

b D ) ≤ 2 βΣ log(2/δ) + (ΣD − Σ D m

√

r 2 σΣD

log(2/δ) . m

By a union bound, both events hold simultaneously with probability at least 1 − 2δ, and combining these with the initial decomposition of ∆Σ,2 yields the claimed high-probability bound.

A.4

Finite-sample bounds

Combining Lemma A.4 with the concentration results of the previous subsection yields the following finite-sample bound. Theorem A.9 (High-probability finite-sample bounds for PIKS). Let λ > 0, γ > 0 and let δ ∈ (0, 1/2). Assume that 36 2n log ≤ λ ≤ ∥Σ∥∞ , n δ

36 2m log ≤ λγ −1 ≤ ∥ΣD ∥∞ . m δ

(5)

Define Σγ := Σ + γΣD ,

B(λ, γ) := (Σγ + λI)−1/2 w∗ ,

and the (value and derivative) capacity terms 1/2

NV (λ, γ) := (Σγ + λI)−1/2 Σ1/2 HS ,

ND (λ, γ) := (Σγ + λI)−1/2 ΣD

Then there exists a constant c > 0, depending only on κV , κD , LV , LD and ∥w∗ ∥, but not on λ, γ, n, m, δ, such that, with probability at least 1 − 2δ, # " r log(4/δ) NV (λ, γ) 4 ND (λ, γ) 1/2 ∗ −1/2 log(4/δ) √ +γ + log +λ B(λ, γ) . Σ (w−w b ) ≤ c λ +γ √ n m δ n m (6) Proof. From the error decomposition Lemma A.4, we have ∥Σ1/2 (w b − w∗ )∥ ≤ ∆2Σ,1 ∆V + ∥w∗ ∥∆Σ,2 where

b γ + λI)−1/2 , ∆Σ,1 := (Σγ + λI)1/2 (Σ

+ λ B(λ, γ),

bγ) , ∆Σ,2 := (Σγ + λI)−1/2 (Σγ − Σ

∆V := (Σγ + λI)−1/2 (Vbγ − Vγ ) . Step 1: bound on ∆Σ,1 . Recall b γ )(Σγ + λI)−1/2 . ∆Σ,3 := (Σγ + λI)−1/2 (Σγ − Σ 28

Applying Lemma A.5 with confidence parameter η, and the condition (3), we obtain an event of probability at least 1 − 2η on which ∆Σ,3 ≤ 1/2. On this event, Lemma A.6 yields ∆2Σ,1 ≤ (1 − ∆Σ,3 )−1 ≤ 2. Step 2: bounds on ∆V and ∆Σ,2 . We use Lemmas A.7 and A.8. From Lemma A.7, for any η ∈ (0, 1), with probability at least 1 − 2η, r r log(2/η) √ log(2/η) log(2/η) log(2/η) √ +γ 2βD , ∆V ≤ 2βV + 2σV + 2σD n n m m | | {z } {z } value part

with

derivative part

βV = λ−1/2 2LV κV ,

σV2 = L2V κ2V (Σγ + λI)−1/2 Σ1/2 HS ,

βD = λ−1/2 2LD κD ,

2 σD = L2D κ2D (Σγ + λI)−1/2 ΣD

1/2 2 . HS

Hence there exist constants aV , bV , aD , bD > 0 depending only on (κV , LV ) and (κD , LD ) such that r r log(2/η) log(2/η) log(2/η) log(2/η) −1/2 ∆V ≤ λ aV +γ aD +γ bD ND (λ, γ) + bV NV (λ, γ) . n m n m Similarly, from Lemma A.8, for the same confidence parameter, with probability at least 1 − 2η, r r log(2/η) log(2/η) log(2/η) log(2/η) −1/2 ∆Σ,2 ≤ λ ãV +γ ãD +γ b̃D ND (λ, γ) + b̃V NV (λ, γ) , n m n m for some constants ãV , ãD , b̃V , b̃D > 0 depending only on κV and κD . Combining these two bounds, we obtain constants A1 , A2 , B1 , B2 > 0, depending only on (κV , κD , LV , LD ), such that r r log(2/η) log(2/η) log(2/η) log(2/η) +A2 γ +B2 γ ND (λ, γ) ∆V ≤ λ−1/2 A1 + B1 NV (λ, γ) , n m n m (7) and r r log(2/η) log(2/η) log(2/η) log(2/η) −1/2 ∆Σ,2 ≤ λ A1 +A2 γ + B1 NV (λ, γ) +B2 γ ND (λ, γ) , n m n m (8) where we have simply taken maxima of the corresponding constants from ∆V and ∆Σ,2 . Step 3: combine all bounds. Set η := δ/3 and apply Lemmas A.5, A.7, and A.8 with confidence parameter η. Each holds with probability at least 1 − 2η, so by a union bound, with probability at least 1 − 3 · 2η = 1 − 2δ, we simultaneously have ∆2Σ,1 ≤ 2,

and the bounds (7)-(8) for ∆V , ∆Σ,2 . 29

On this event,

∥Σ1/2 (w b − w∗ )∥ ≤ ∆2Σ,1 ∆V + ∥w∗ ∥∆Σ,2 + λ B(λ, γ) ≤ 2 ∆V + ∥w∗ ∥∆Σ,2 + λ B(λ, γ).

Substituting the bounds (7)-(8), we obtain " log(2/η) log(2/η) 1/2 ∗ ∗ + A2 γ ∥Σ (w b − w )∥ ≤ 2(1 + ∥w ∥) λ−1/2 A1 n m # r r log(2/η) log(2/η) + B1 NV (λ, γ) + B2 γ ND (λ, γ) + λ B(λ, γ). n m Since η = δ/3, we have log(2/η) = log(6/δ) ≤ c1 log(4/δ) for some constant c1 > 0. Absorbing √ the constants 2(1 + ∥w∗ ∥)A1 , 2(1 + ∥w∗ ∥)A2 , 2(1 + ∥w∗ ∥)B1 , 2(1 + ∥w∗ ∥)B2 and c1 into a single ∗ constant c > 0 depending only on κV , κD , LV , LD and ∥w ∥, we rewrite the bound as # " r log(4/δ) NV (λ, γ) ND (λ, γ) 4 1/2 ∗ −1/2 log(4/δ) √ +γ + +γ √ log +λ B(λ, γ) , ∥Σ (w−w b )∥ ≤ c λ n m δ n m which is exactly (4).

A.5

Learning rate acceleration from differential information

We now instantiate the finite-sample bound under the Assumption (A6) and optimize over λ, γ. Corollary A.1 (Learning rates for PIKS). Assume Assumptions (A1)-(A6) hold. Then there exist choices of regularization parameters λ = λ(n, m) and γ = γ(n, m), given explicitly in the proof and depending polynomially on n, m, log n, and log m, such that, with high probability:  1 2(1−α) 1−r n− 4 m− 8 , m ≲ n (1+α)(1−r) , 1/2 ∗ ∥Σ (w b − w )∥ ≲ 2(1−α) 1  − 2(1+α) n , m ≳ n (1+α)(1−r) . Proof. The proof proceeds in four steps. Step 1: Capacity bounds. As a consequence of Assumption (A6), using Σγ + λI ⪰ Σ1 + λI and Σγ + λI ⪰ γΣD + λI, we obtain NV (λ, γ)2 = Tr Σ1 (Σγ + λI)−1 + Tr Σ2 (Σγ + λI)−1 ≲ λ−α + γ r−1 λ−r . (9) Moreover, since Σγ + λI ⪰ λI, ND (λ, γ)2 ≤ λ−1 Tr(ΣD ) ≲ λ−1 . Step 2: Error bound. From Assumption (A4), the model is well specified, i.e. ∥w∗ ∥H < +∞. Then √ B(λ, γ)2 = ⟨w∗ , (Σγ + λI)−1 w∗ ⟩ ≤ λ−1 ∥w∗ ∥2H , λ B(λ, γ) ≲ λ.

Throughout the remainder of the proof, we restrict to regularization parameters satisfying log n log m max , γ ≲ λ ≲ γ. n m √ √ √ This ensures the admissibility condition of Theorem 3.1. Under this regime, using a + b ≤ a + b, the finite-sample bound simplifies to r−1

∥Σ

1/2

√ λ−α/2 γ 2 λ−r/2 γ λ−1/2 √ (w b − w )∥ ≲ √ + + √ + λ. n n m ∗

(10)

This bound makes explicit how derivative information can reduce the effective value complexity through the aligned component Σ2 , while simultaneously introducing a variance cost that grows with γ through the derivative observations. Step 3: Optimization over γ. For fixed λ, only the middle two terms in (10) depend on γ. Increasing γ reduces the variance associated with estimating the function values in directions constrained by the derivative operator (by shrinking the contribution of the corresponding eigenspaces), but simultaneously amplifies the variance coming from the noisy estimation of derivative information. The optimal choice of γ is obtained by balancing these two effects by solving r−1

γ 2 λ−r/2 γ λ−1/2 √ = √ . n m This is equivalent to r−3

γ 2

which yields γ ⋆ (λ) ≍

n 1/2 m

r−1

λ 2 ,

1 m 3−r

1−r

λ 3−r . n Substituting (11) into (10), the two γ–dependent terms become equal and we obtain √ ∥Σ1/2 (w b − w∗ )∥ ≲ A(λ) + B(λ) + λ, where

λ−α/2 A(λ) := √ , n 1

1−r

(11)

(12)

1+r

B(λ) := m− 2(3−r) n− 3−r λ− 2(3−r) .

1−r

Substituting γ ⋆ (λ) ≍ (m/n) 3−r λ 3−r into the working regime max{log n/n, γ log m/m} ≲ λ ≲ γ yields the equivalent λ-only conditions r 3−r m log n m 21 log m 2 ≲ λ ≲ max , . n n m n √ Step 4: Optimization over λ and regime analysis. We choose λ by balancing the bias λ with one of the variance terms in (12). (i) Balance with A(λ). We solve

√

λ =

λ−α/2 √ , n

which gives

λ⋆A ≍ n− 1+α ,

λ⋆A ≍ A(λ⋆A ) ≍ n− 2(1+α) .

(13)

(ii) Balance with B(λ). We solve √ 1+r 1−r 1 λ = m− 2(3−r) n− 3−r λ− 2(3−r) , which gives

1−r

λ⋆B ≍ m− 4 n− 2 ,

1−r

λ⋆B ≍ B(λ⋆B ) ≍ n− 4 m− 8 .

(14)

Validity of the regimes. The bound (12) is controlled by the dominant variance term at the chosen λ. Saturation regime. We choose λ = λ⋆A when B(λ⋆A ) ≲ A(λ⋆A ). Substituting λ⋆A ≍ n−1/(1+α) into B(λ) gives 1+2α−r

1−r

B(λ⋆A ) ≍ m− 2(3−r) n− 2(1+α)(3−r) . Hence

B(λ⋆A ) ≲ A(λ⋆A )

is equivalent to

2(1−α)

m ≳ n (1+α)(1−r) . Derivative-limited regime. We choose λ = λ⋆B when A(λ⋆B ) ≲ B(λ⋆B ). which yields the same transition value for m. Moreover, to ensure our working regime condition, we take the regularization parameter to be clipped as ⋆ log n ⋆ ⋆ log m λ := max λ , , γ (λ ) , n m where λ⋆ ∈ {λ⋆A , λ⋆B } is the optimizer in the corresponding regime. One can check that this choice of λ satisfies the admissibility conditions, including the upper condition λ ≲ γ, in both regimes. Since the variance terms are nonincreasing in only increases the bias term, yielding an pλ, clippingp additional lower-order contribution of order log n/n + γ log m/m. Resulting rate.

Let

2(1−α)

mcrit ≍ n (1+α)(1−r) . The optimized learning rate is ∥Σ1/2 (w b − w∗ )∥ ≲

 1 1−r n− 4 m− 8 , n

1 − 2(1+α)

m ≲ mcrit , m ≳ mcrit .

A.6

Learning rates in physically consistent norms

We finally show that the estimator also converges in the stronger norm induced by Σ + ΣD , which jointly controls the prediction error and the error on the differential quantities. Corollary A.2 (Learning rates in the physically consistent norm). Assume Assumptions (A1)-(A6) hold. Then there exist choices of regularization parameters λ = λ(n, m) and γ = γ(n, m) such that, with high probability,  1 1−r n− 4 m− 8 + m−1/4 , m ≲ mcrit , ∥(Σ + ΣD )1/2 (w b − w∗ )∥ ≲ 1 n− 2(1+α) + m−1/4 , m ≳ mcrit , where

2(1−α)

mcrit ≍ n (1+α)(1−r) . Proof. Since γ ≥ 1, we have and therefore

Σγ = Σ + γΣD ⪰ Σ + ΣD ,

∥(Σ + ΣD )1/2 (w b − w∗ )∥ ≤ ∥Σ1/2 b − w∗ )∥. γ (w 1/2

Finite-sample bound. Repeating the proof of Theorem 3.1 with Σ1/2 replaced by Σγ yields the same finite-sample bound, since both the estimation and bias terms are handled exactly as before, 1/2 but using ∥Σγ (Σγ + λI)−1/2 ∥ ≤ 1 instead of ∥Σ1/2 (Σγ + λI)−1/2 ∥ ≤ 1. Then, the same argument as in the proof of Corollary 3.1 gives r−1

∥Σ1/2 b − w∗ )∥ ≲ γ (w

γ 2 λ−r/2 γ λ−1/2 √ λ−α/2 √ √ + + √ + λ, n n m

up to logarithmic factors and lower-order terms. The optimization over γ is identical to that of Corollary 3.1, except that we now impose the constraint γ ≥ 1. Let 1 m 3−r 1−r γ ⋆ (λ) ≍ λ 3−r n denote the unconstrained optimizer from Corollary 3.1. Under the constraint γ ≥ 1, we therefore choose γ(λ) ≍ max{1, γ ⋆ (λ)}. Optimizing γ and λ. When γ ⋆ (λ) ≥ 1, we recover exactly the same bound as in Corollary 3.1. When γ ⋆ (λ) < 1, the constraint is active and we set γ = 1, which yields the additional term λ−1/2 √ . m Balancing this term with the bias term

√

λ gives λ ≍ m−1/2 , 33

and therefore an additional contribution of order m−1/4 . Therefore, the optimal regularization parameter is λ ≍ max{λ⋆A , λ⋆B , m−1/2 }, where λ⋆A and λ⋆B are the optimizers obtained in the proof of Corollary 3.1. This yields  1 1−r n− 4 m− 8 + m−1/4 , m ≲ mcrit , ∥(Σ + ΣD )1/2 (w b − w∗ )∥ ≲ 1 n− 2(1+α) + m−1/4 , m ≳ mcrit , as claimed.

A.7

Adaptation to the vector-valued case

The setting of this paper considers scalar-valued differential operators, for simplicity of the argument. The adaptation to vector-valued operator is straightforward as we detail in this section. Consider a vector-valued differential operator D = (D1 , . . . , Dk ), where each Di is of the form given in Assumption (A1). Then, for any u ∈ H and any x ∈ X , Du(x) is a vector in Rk . The differential data points dj are also in Rk and we replace the risk (1) by the following one: n m 2 1X 1 X b R(u) = u(xi ) − yi + γ ∥Du(zj ) − dj ∥22 + λ∥u∥2H . n i=1 m j=1

(15)

We can define component-wise feature maps Dℓ ϕ(z), and define the differential covariance operator as # " k X ΣD := Ez∼ρD Dℓ ϕ(z) ⊗ Dℓ ϕ(z) . (16) ℓ=1

The key point is that the analysis is carried out using covariance operators Σ and ΣD , and the key assumption, Assumption (A6), is an abstract condition on Σ and ΣD , regardless of their expressions. As such, even if going from scalar-valued to vector-valued changes the expression of the covariance operator ΣD as we just saw, as long as one guarantees that Assumption (A6) is satisfied, the rest of the analysis holds. The other assumptions must be adapted in a straightforward way: the smoothness s required in Assumption (A2) is the maximum of the orders of each Dℓ , and for Assumption (A3) we need the boundedness of each feature ∥Dℓ ϕ(z)∥H . b D , VbD by summing over ℓ: For the proofs, we also adapt VD , Σ " k # X VD := E dℓ Dℓ ϕ(z) ℓ=1 m X k X b D := 1 Σ Dℓ ϕ(zj ) ⊗ Dℓ ϕ(zj ) m j=1 ℓ=1

m k 1 XX VbD := dj,ℓ Dℓ ϕ(zj ) m j=1 ℓ=1

Then, the concentration inequalities are still valid and the proofs are the same. 34

A.8

Learning rates with exact physical constraint

We derive the PIKS rate in the idealized setting where the differential information is known exactly. In this physical oracle setting, the value observations remain empirical, while the differential part of the risk is replaced by its population counterpart: n 2 1X b Ror (u) = u(xi ) − yi + γ Ez∼ρD (Du(z) − Du∗ (z))2 + λ∥u∥2H . n i=1

b D and VbD by their population counterparts ΣD and Equivalently, in operator form, this replaces Σ VD . Corollary A.3 (Oracle learning rate). Assume Assumptions (A1)-(A6) hold, with r < 1. Then there exist choices of λ = λ(n) and γ = γ(n) such that, with high probability, 1

∥Σ1/2 (w bor − w∗ )∥ ≲ n− 2(1+α) . Proof. The oracle estimator has the closed form b + γΣD + λI)−1 (Vb + γVD ). w bor = (Σ Let

b or := Σ b + γΣD , Σ γ

Vbγor := Vb + γVD ,

Σγ := Σ + γΣD .

Under Assumptions (A4) and (A5), V = Σw∗ ,

VD = ΣD w∗ ,

V + γVD = Σγ w∗ .

The proof of Theorem 3.1 applies verbatim, except that the differential part is no longer empirical. Hence the only stochastic deviations are b − Σ, Σ

Vb − V,

b D − ΣD and VbD − VD disappear. Up to logarithmic factors, and all terms involving Σ ∥Σ1/2 (w bor − w∗ )∥ ≲ where

N (λ, γ) 1 √ + V√ + λB(λ, γ), n n λ

NV (λ, γ) = (Σγ + λI)−1/2 Σ1/2 HS ,

B(λ, γ) = (Σγ + λI)−1/2 w∗ .

By Assumption (A6), as in the proof of Corollary 3.1, NV (λ, γ)2 ≲ λ−α + γ r−1 λ−r . Moreover, by attainability,

√ λB(λ, γ) ≲

Therefore,

λ. r−1

1/2

∥Σ

1 λ−α/2 γ 2 λ−r/2 √ √ (w bor − w )∥ ≲ √ + √ + + λ. n n n λ ∗

Since the differential information is known at the population level, there is no derivative-sample variance term increasing with γ. We choose γ large enough so that γ r−1 λ−r ≲ λ−α . Equivalently, when r > α, it suffices to take r−α

γ ≳ λ− 1−r , while when r ≤ α, any γ ≳ 1 is enough. With this choice, ∥Σ1/2 (w bor − w∗ )∥ ≲

λ−α/2 √ 1 √ + √ + λ. n n λ

Balancing the dominant variance term with the bias, √ λ−α/2 √ = λ, n gives

λ ≍ n− 1+α . For this choice,

while

√ 1 λ−α/2 √ ≍ λ ≍ n− 2(1+α) , n 1 1 √ = n−1+ 2(1+α) n λ

is lower order since α > 0. Hence 1

∥Σ1/2 (w bor − w∗ )∥ ≲ n− 2(1+α) .

Examples

We now illustrate Assumption (A6) on concrete examples where the decomposition can be computed explicitly.

B.1

Partial Laplacian and periodic Sobolev RKHS

Lemma B.1 (Capacity decomposition spaces). Let X = Td , let k be a Matérn kernel P for Sobolev 2 of smoothness s > d/2, and let D = i∈S ∂xi be a partial Laplacian with S ⊆ {1, . . . , d}. Assume that value and differential samples are drawn uniformly on Td . Then the value covariance operator Σ admits a decomposition Σ = Σ1 + Σ2 satisfying Assumption (A6), with exponents d − |S| d α= , r= . 2s 2s 36

Proof. Step 1: diagonalization in the Fourier basis. Let D be a constant-coefficient differential operator of order q, of the form X D= cα ∂ α , ∂ α = ∂xα11 · · · ∂xαdd , |α|≤q

with real coefficients cα . Its Fourier symbol is given by X P (k) = cα (2πik)α ,

k ∈ Zd ,

|α|≤q

so that De2πik·x = P (k) e2πik·x . Since the kernel k is translation invariant and both x and z are sampled uniformly on Td , the covariance operators Σ = E[ϕ(x) ⊗ ϕ(x)],

ΣD = E[Dϕ(z) ⊗ Dϕ(z)]

are convolution operators and are therefore diagonal in the Fourier basis {ek (x) := e2πik·x }k∈Zd . The corresponding eigenvalues are σk = µk ,

τk = |P (k)|2 µk .

Step 2: partial Laplacian and visible/invisible frequencies. For D = symbol is P (k) = −(2π)2 ∥kS ∥2 , so τk = (2π)4 ∥kS ∥4 µk ,

2 i∈S ∂xi , the

where kS denotes the restriction of k to coordinates in S. Define the index sets Z := {k ∈ Zd : kS = 0},

Z c := Zd \ Z.

Thus Z consists of frequencies invisible to D (since P (k) = 0), whereas Z c corresponds to frequencies detectable by D. Step 3: the decomposition Σ = Σ1 + Σ2 . Define X X Σ1 := µk ek ⊗ ek , Σ2 := µk ek ⊗ ek . k∈Z c

k∈Z

Then Σ = Σ1 + Σ2 and Σ1 , Σ2 ⪰ 0. Step 4: capacity of the invisible component Σ1 . Let d0 = d − |S| and write k−S ∈ Zd−S for the subvector of coordinates outside S. The map k−S 7→ (k−S , 0S ) is a bijection Zd−S → Z and |(k−S , 0S )| = |k−S |, hence X k∈Z

µk = µk + λ

µ(k−S ,0S ) , µ (k−S ,0S ) + λ d−S

X k−S ∈Z

µ(k−S ,0S ) ≍ (1 + |k−S |2 )−s .

Therefore, by a standard comparison between lattice sums and integrals, for all λ ∈ (0, 1], X µk d0 Tr Σ1 (Σ1 + λI)−1 = ≲ λ− 2s . µk + λ k∈Z

This verifies the first trace condition with exponent α = d0 /(2s). Step 5: alignment of Σ2 with ΣD . On Z c we have ∥kS ∥ ≥ 1 and hence τk = (2π)4 ∥kS ∥4 µk ≥ c µk for some c > 0. Therefore, for all t ∈ (0, 1], X X µk µk ≤ . Tr (ΣD + tI)−1 Σ2 = τ + t cµ k k +t c c k∈Z

k∈Z

Using µk ≍ (1 + |k|2 )−s and standard comparison between lattice sums and integrals gives X µk d ≲ t− 2s , t ∈ (0, 1], µ + t k d k∈Z

and hence

Tr (ΣD + tI)−1 Σ2

≲ t− 2s .

This verifies the second trace condition with exponent r = d/(2s). Combining Steps 3-5 yields Assumption (A6) with the claimed exponents.

B.2

Gradients

Let Ω ⊂ Rd be a bounded, connected, Lipschitz domain. Fix s > d2 + 1. Let H = H s (Ω), with norm ∥ · ∥H equivalent to the standard H s (Ω) norm. We denote by K the reproducing kernel associated to H, and by ϕ : Ω → H the feature map defined, for all x ∈ Ω, by ϕ(x) = K(x, ·) ∈ H. Choose ρ = ρD the uniform distribution on Ω. B.2.1

Covariance operator decomposition

Definition of Σ

Define the covariance operator Σ := Ex∼ρ [ϕ(x) ⊗ ϕ(x)].

We have, for all u ∈ H, Definition of ΣD

⟨Σu, u⟩H = ∥u∥2L2 (ρ) .

(17)

Define the gradients covariance operator # " d X ΣD = Ez∼ρD ∂i ϕ(z) ⊗ ∂i ϕ(z) . i=1

We have, for all u ∈ H, Definition of Σ1 that Denote Then for all u ∈ H,

⟨ΣD u, u⟩H = ∥∇u∥2L2 (ρ) .

(18)

Denote by ϕΩ = Ex∼ρ [ϕ(x)] ∈ H the representer of the averaging functional, so ⟨u, ϕΩ ⟩H = Ex∼ρ [u(x)] =: uΩ . Σ1 = ϕΩ ⊗ ϕΩ . ⟨Σ1 u, u⟩H = u2Ω . 38

(19)

Definition of Σ2

Next, define the bounded operator A : H → L2 (ρ),

Set

Au := u − uΩ ,

Σ2 := A∗ A,

so that for all u ∈ H,

⟨Σ2 u, u⟩H = ∥u − uΩ ∥2L2 (ρ) = Varx∼ρ (u(x)).

(20)

Proposition B.2 (Decomposition of Σ). One has Σ = Σ1 + Σ2 . Proof. Let u ∈ H. By the variance decomposition formula, 2 Ex∼ρ u(x)2 = Ex∼ρ [u(x)] + Varx∼ρ (u(x)). Equivalently,

∥u∥2L2 (ρ) = u2Ω + ∥u − uΩ ∥2L2 (ρ) = u2Ω + Varx∼ρ (u(x)).

Using (17), (19), and (20), we obtain ⟨Σu, u⟩H = ⟨Σ1 u, u⟩H + ⟨Σ2 u, u⟩H for every u ∈ H. Since Σ, Σ1 , and Σ2 are bounded self-adjoint operators on H, equality of their quadratic forms implies, by polarization, that Σ = Σ1 + Σ2 .

B.2.2

Computation of the coefficients

Let us first bound the effective dimension N1 (λ) := Tr Σ1 (Σ1 + λI)−1 . Since Σ1 has rank one, the effective dimension N1 (λ) is uniformly bounded. In particular, it satisfies N1 (λ) ≤ Cα λ−α , λ ∈ (0, 1], with Let us now bound

α=0. Tr (ΣD + λI)−1 Σ2 .

Since Ω is bounded, connected, and Lipschitz, the Poincaré–Wirtinger inequality [26] holds: there exists C > 0 such that, for every u ∈ H 1 (Ω), ∥u − uΩ ∥L2 (ρ) ≤ C∥∇u∥L2 (ρ) . 39

Therefore, for every u ∈ H,

⟨Σ2 u, u⟩H = ∥u − uΩ ∥2L2 (ρ) ≤ C 2 ∥∇u∥2L2 (ρ) = C 2 ⟨ΣD u, u⟩H .

Equivalently, in the Loewner order, It follows that

0 ⪯ Σ2 ⪯ C 2 ΣD .

Tr (ΣD + λI)−1 Σ2 ≤ C 2 Tr ΣD (ΣD + λI)−1 .

We now bound the last trace. Let (ηj )j≥1 be the nonzero eigenvalues of ΣD , arranged in nonincreasing order. Since ΣD is the covariance operator associated with the map H s (Ω) → L2 (ρ; Rd ),

u 7−→ ∇u,

the standard eigenvalue estimate for Sobolev embeddings on bounded Lipschitz domains gives ηj ≤ CD j −2(s−1)/d . Set β :=

2(s − 1) . d

Because s > d2 + 1, we have β > 1. Hence, for 0 < λ ≤ 1, X Tr ΣD (ΣD + λI)−1 = j≥1

≤

X j≥1

ηj ηj + λ CD j −β CD j −β + λ

≲ λ−1/β . Since 1/β = d/(2(s − 1)), we obtain Tr (ΣD + λI)−1 Σ2 ≲ λ−d/(2(s−1)) . Thus one may take r=

d . 2(s − 1)

In particular, since s > d2 + 1, one has r < 1.

Detailed experimental setup

In this section we provide more details about the setup used for the experiments of Section 5.

Implementation To maximize the flexibility of our system, we implemented the PIKS estimator using the Jax framework, which allows to efficiently compute arbitrary derivatives. While very flexible, any automatic differentiation framework does introduce some computational overhead and requires particular care when implementing kernel functions. For example, Matérn kernels require differentiating through a square-root which is numerically unstable at 0. We introduce a small additive nugget term to ensure stability everywhere. The kernel solver is a straightforward implementation of the equations in Lemma A.2 which yield a space complexity of O((n + m)2 ) and time complexity of O((n + m)3 ). Bounded domain gradient example We repeated the experiment 10 times to obtain standard deviations reported in Fig. 2. We set the variance on the differential data to be low in order to avoid having to scale m too much before seeing a saturation effect. All data points were uniformly sampled on a unit disk with 1000 samples used for validation (to select hyperparameters) and 10000 for computing the test error. Partial Laplacian example We give a few more details on the kernel and target functions for the example on Sobolev spaces on the torus with the partial Laplacian operator. Here we used X = T4 . We sampled uniformly at random F = 64 frequency vectors kℓ ∈ Z4 such that for ℓ ∈ [0, 15], kℓ ∼ [U([−8, 8]\{0}), 0, 0, 0]; for ℓ ∈ [16, 31], kℓ ∼ [0, √ U([−8, 8]\{0})] and so on (only one dimension active for each vector). For feature map ϕℓ (x) = 2 cos(2π⟨kℓ , x⟩), target function and kernel are defined as F F X X u∗ (x) = u0 + cℓ ϕℓ (x), k(x, x′ ) = 1 + µℓ ϕℓ (x)ϕℓ (x′ ), ℓ=1

ℓ=1

with the smoothness coefficients µℓ ≍ (1 + ∥kℓ ∥2 )−β ,

cℓ ≍ (1 + ∥kℓ ∥2 )−δ ,

β > 1,

δ>

1 . 2

In particular, we set β = δ = 4.1 to have a well-specified problem in a smooth enough space. We use additive Gaussian noise, with standard-deviation equal to 10% to the range of the data. This was done because the different data components have vastly different numerical ranges, and setting a fixed variance would have introduced artifacts in the results. The best-fit lines of Fig. 3 are obtained by least-squares regression of the experimental data against the function anb to find the exponential coefficient b. Then, using Corollary 3.1, we get 1 α = − 2b − 1. Hyperparameters The two hyperparameters of PIKS (λ and γ), as well as the hyperparameters of the kernel (notably the length-scale of the Matérn kernel) were determined by a coarse grid-search using small, noisy validation sets to estimate their performance at generalization time.

Record · ID 965410 · SHA-256 5bb4cd57966310a1
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.