Adaptive meshfree solution of linear partial differential equations with PDE-greedy kernel methods

We consider the meshless solution of PDEs via symmetric kernel collocation by using greedy kernel methods. In this way we avoid the need for mesh generation, which can be challenging for non-standard domains or manifolds. We introduce and discuss different kind of greedy selection criteria, such as the PDE-P -greedy and the PDE-f -greedy for collocation point selection. Subsequently we analyze the convergence rates of these algorithms and provide bounds on the approximation error in terms of the number of greedily selected points. Especially we prove that target-data dependent algorithms, i.e. those using knowledge of the right hand side functions of the PDE, exhibit faster convergence rates. The provided analysis is applicable to PDEs both on domains and manifolds. This fact and the advantages of target-data dependent algorithms are highlighted by numerical examples.


Introduction
Kernel methods [6,32,37] are a popular class of techniques in approximation theory, numerical mathematics and machine learning. They revolve around the notion of a kernel k. Given a non-empty set Ω, a real valued kernel is defined as a symmetric function k : Ω × Ω → R, and here we focus on Ω ⊂ R d .
For given data X n := {x 1 , . . . , x n } ⊂ Ω, the kernel matrix K Xn ∈ R n×n is defined as (K Xn ) ij := (k(x i , x j )) ij . A kernel is called strictly positive definite if the associated kernel matrix is positive definite for any choice of pairwise distinct points X n ⊂ Ω, n ∈ N. Those strictly positive kernels give rise to a unique Hilbert space, the so called Reproducing Kernel Hilbert Space H k (Ω) (RKHS). The RKHS is a space of real-valued functions from Ω to R with inner product ·, · H k (Ω) , and it is sometimes also called native space. Within the RKHS, the kernel k acts as a reproducing kernel, i.e. the following two properties hold: 1. k(·, x) ∈ H k (Ω), 2. f (x) = f, k(·, x) H k (Ω) ∀x ∈ Ω, ∀f ∈ H k (Ω).
Based on these properties, kernel interpolation can be characterized and analysed: For a given closed subspace V ⊂ H k (Ω) we denote the H k (Ω)-orthogonal projection onto V with Π V : H k (Ω) → V . Given n interpolation points X n , the unique minimum-norm interpolant s n to a function f ∈ H k (Ω) is equal to the H k (Ω)-orthogonal projection of f onto the subspace V Xn := span{k(·, x i ), x i ∈ X n } ⊂ H k (Ω), i.e. it holds s n (·) := Π VX n (f )(·) = n j=1 α j k(·, x j ).
The coefficients α j can be determined by the interpolation conditions s n (x i ) = f (x i ) ∀i = 1, ..., n ⇔ K Xn α = (f (x j )) n j=1 , where α = (α 1 , ..., α n ) ⊤ ∈ R n is the vector of the coefficients. A standard way to quantify the error between the function f and the interpolant s n in the · L ∞ (Ω) -norm is using the so called power function P Xn , which is given as P Xn (x) := k(·, x) − Π VX n (k(·, x)) H k (Ω) = sup and thus allows to estimate the residual r n := f − s n as In order to obtain a good approximation, a suitable choice of points X n ⊂ Ω is required, however usually unclear a priori. Therefore, greedy algorithms [34] are used which start with an empty set X 0 = ∅ and then add further points step by step as X n+1 := X n ∪ {x n+1 }. The next interpolation point x n+1 is usually chosen according to some optimality criterion. For greedy kernel interpolation, there are the so called P -greedy [2], f ·P -greedy (or psr -greedy) [4], f -greedy [31], and f /P -greedy [22] criteria, and they choose the next point according to the following strategies: i. P -greedy: x n+1 := arg max x∈Ω P n (x), ii. f · P -greedy: x n+1 := arg max x∈Ω |r n (x)| · P n (x), iii. f -greedy: x n+1 := arg max x∈Ω |r n (x)|, iv. f /P -greedy: x n+1 := arg max x∈Ω\Xn |r n (x)|/P n (x).
A unifying analysis of those standard greedy kernel interpolation algorithms was recently provided in [39] under the notion of β-greedy algorithms, that we will generalize in the current article.
The theory of kernel interpolation is not limited to the interpolation of function values. Indeed it is possible to consider arbitrary functionals Λ n := {λ 1 , ..., λ n } ∈ H k (Ω) ′ and obtain a function s n of minimum norm that satisfies where H k (Ω) ′ denotes the dual space of H k (Ω). This procedure is called generalized interpolation [37,Chapter 16]. It is possible to show that similarly to Eq. (1) it holds where we now have V n := span{v λi | λ i ∈ Λ n } and v λ1 , ..., v λn ∈ H k (Ω) are the unique Riesz representers of the functionals λ 1 , ..., λ n . Observe that for the special case of standard kernel interpolation, due to the reproducing property, the Riesz representer of the point evaluation functionals δ xj are simply given by k(·, x j ), such that the generalized interpolant from Eq. (5) again reduces to Eq.
This framework of generalized interpolation can be applied to solving linear PDEs [21,27]. For this we focus on linear elliptic PDEs with an elliptic differential operator L and a suitable boundary operator B, which will be discussed in more detail in Section 2. We are going to focus on Dirichlet boundary conditions, i.e. B = Id, and remarks will be provided that elaborate on different boundary conditions such as Neumann or mixed boundary conditions. The constraints imposed by the PDE can be collected in two sets of functionals, i.e. we introduce Due to Eq. (5), the Riesz representers for these functionals are required. We recall that under suitable conditions on the kernel k, the Riesz representer of the functional δ x • L is given by L (2) k(·, x), where L (2) means that L is applied to the second argument of k. The Riesz representer for the boundary functionals in the case of Dirichlet boundary conditions is simply given by k(·, x), x ∈ ∂Ω, due to the reproducing property.
In order to provide a gentle illustrative example, we focus in the remaining paragraph on the well known Poisson problem, i.e. we choose L as the negative Laplacian L = −∆. In this case, the required Riesz representers are given by −∆ (2) k(·, x) for x ∈ Ω and k(·, x) for x ∈ ∂Ω. Thus the overall equation system for given collocation points x 1 , ..., x n ∈ Ω ∪ ∂Ω is given by (in case of a suitable ordering of the points, rows and columns) with f | Xn = (f (x i )) n i=1 ∈ R n , similar for g| Xn and block matrices defined as The generalized interpolant s n finally has the form We remark that there are also other methods to solve PDEs via kernel methods and we elaborate on them in Subsection 2.4. The approach explained here so far is called symmetric collocation [5]. In the remainder of the presentation, we advance by extending the unifying framework and analysis of standard greedy kernel interpolation algorithms from [39] to generalized greedy kernel interpolation with focus on PDE approximation [1,[9][10][11].
The paper is structured as follows. Section 2 starts by giving required background information on PDEs, greedy kernel interpolation, generalized interpolation and pointing to the corresponding literature. Section 3 provides an analysis of Kolmogorov n-widths related to the considered PDE problems with focus on Dirichlet boundary conditions. Section 4 extends the standard target-data dependent greedy analysis of [39] to generalized kernel interpolation. Section 5 combines the results of the previous sections to derive convergence rates of the error u − s n L ∞ (Ω) between the approximant and the true solution of the PDE. Section 3 analyzes the optimality of the rates on the Kolmogorov widths presented before. Section 6 concludes with numerical experiments and Section 7 gives a conclusion and formulates an outlook.

Background results
In the following, we collect more required background information about partial differential equations (Subsection 2.1), interpolation with translational invariant kernels (Subsection 2.2) and greedy as well as generalized interpolation (Subsection 2.3) as well as other collocation methods to solve PDEs (Subsection 2.4).

Partial Differential Equations
Partial differential equations (PDEs) [21,27] are of interest since they describe various physical phenomena. We focus on elliptic boundary value problems, i.e.
where L is an elliptic operator and B j are boundary operators. We require multiple boundary operators B j for well-posedness of higher-order PDEs, while a single boundary operator suffices for the case of second-order PDEs. One specific example for such a boundary problem is the Poisson problem, i.e. m = 1, here L = −∆ is the Laplace operator and Id is the point evaluation on the boundary [20,Chapter 1.3]. The most important reason why we focus on elliptic problems is that they are well-posed (in the sense of Hadamard) under suitable conditions, i.e. 1) there exists a solution, 2) the solution is unique, 3) the solution depends continously on the data. The third property 3) especially implies a well-posedness ("WP") inequality, i.e.
Here · W P , · W PL and · W PB j denote some norms and C > 0. These inequalities will be important to prove convergence rates. The most important well-posedness inequality is a consequence of a theorem from PDE theory [20,Chapter 2]: If {L, B j } and Ω, ∂Ω are "nice" enough as in [20,Chapter 2] and deg(L) =: 2m, deg(B j ) =: m j ≤ 2m − 1 , 0 ≤ j ≤ m − 1, and if we assume for some r ∈ N 0 , then there exists a unique solution u ∈ H 2m+r (Ω) such that Eq. (10) holds. For details on what "nice" means and for other definitions like deg(L) of an operator L we refer the reader to [20]. Observe that if (11) holds for a certain r ∈ N 0 , then it holds also for any 0 ≤ r ′ ≤ r.
Then there exists a constant C r > 0 such that Inserting u − s n for u and f − Ls n for f (because u − s n is the solution of the linear PDE L(u − s n ) = f − Ls n in Ω, and B j u − B j s n = g j − B j s n on ∂Ω) gives This is a consequence of the theorem stated in [20, Chapter 2, Theorem 5.4], if one also assumes existence and uniqueness of the solution of the PDE. For parabolic problems a similiar statement can be made [21] and the analysis introduced later on works analogously.
Kernel methods are intrinsically related to the use of the · L ∞ (Ω) norm because of the reproducing properties of kernels. Therefore we also recall the maximum principle, which gives estimates for approximation in the · L ∞ (Ω) norm when using Dirichlet boundary conditions: Consider an elliptic operator L of second order and a function u ∈ C 2 (Ω) such that Eq. (10) with m = 1, B = Id holds. Then we obtain the maximum principle [27,Theorem 4.11] using the same argumentation as in (12). This will allow us to conclude convergence rates for the greedy methods that are analyzed in subsequent sections.

Interpolation with translational invariant kernels
Kernels are widely used because they appear naturally in different areas of mathematics. Particularly worth mentioning are the translation invariant kernels, i.e. when there exists a function Φ : holds. The class of translation invariant kernels has nice numerical as well as theoretical properties, including especially the statement of Bochner's theorem: If k is a translation invariant kernel then there exists a measure µ : where B(R d ) denotes the Borel σ-algebra of R d . The measure µ is called the spectral density of k. If the measure µ has a density function with respect to the Lebesgue measure λ, i.e. there exists a function m : R d → R such that dµ(ξ) = m(ξ)dλ(ξ) holds, we call m the spectral density function of k. Since we demand that k is real valued, the spectral density of k needs to be symmetric. Translation invariant kernels are also closely related to Sobolev spaces: Assume that a translation invariant kernel has a spectral density function m : R d → R such that there exists c, C > 0, τ > d/2 for all ξ ∈ R d . Then the native space H k (R d ) is norm equivalent to the Sobolev space H τ (R d ) = W τ 2 (R d ) of possibly fractional order τ > d/2 [19]. This norm equivalence can be extended to RKHS H k (Ω) on bounded domains Ω ⊂ R d under suitable conditions on the boundary ∂Ω, such as a Lipschitz boundary. In this case there is also a norm equivalence of the RKHS H k (Ω) to the Sobolev space H τ (Ω). We call those kind of kernels Sobolev kernels of smoothness τ . A special subclass of the translation invariant kernels are the radial-basis-function (RBF) kernels, i.e. there exists a function ϕ : [0, ∞) → R such that holds. Some examples of radial basis functions ϕ which give rise to RBF kernels and which are also used in the numerical experiments in Section 6 are given by In order to estimate the error f − s n L ∞ (Ω) between the true function f ∈ H k (Ω) and the kernel interpolant s n , one usually makes use of sampling inequalities and corresponding zero Lemmas, which exist for Sobolev spaces or spaces of analytic functions [23,38]. They bound the error in terms of the fill distance h = h X,Ω , which describes the largest hole within Ω where no point from X exists, i.e.
h := h X := h X,Ω := sup The same definition of the fill distance can also be used for the boundary ∂Ω of Ω. Here we recall [13, Theorem 2.2]: Suppose Ω ⊂ R d is a bounded domain satisfying an interior cone condition and having a Lipschitz boundary. Let X ⊂ Ω be a discrete set with Here, | · | W m q (Ω) denotes the Sobolev semi-norm, which uses only the highest derivatives. By using asymptotically quasi uniform points (i.e. c 1 n −1/d ≤ h Xn ≤ c 2 n −1/d for some c 1 , c 2 > 0, which we will abbreviate as h n ≍ n −1/d ), this statement can be transformed into a statement based on the number of interpolation points. For Sobolev kernels in the sense of Eq. (14), the sampling inequality from Theorem 1 can be applied to the residual u − s n , whereby s n interpolates u pointwise in X = X n . This yields . When using m = 0, q = ∞, i.e. bounds in the · L ∞ (Ω) norm, it was recently shown in [40] that this convergence rate for interpolation in RKHS (in contrast to interpolation in Sobolev spaces) does essentially not depend on the shape of the domain. In particular a Lipschitz boundary is not required. We point to [40,Theorem 7] for a detailed formulation.
Similar statements as in Theorem 1 can also be derived for interpolation on manifolds under some assumptions on the shape of the manifold. In the following we are particularly interested in statements for the boundary of the domain Ω, i.e. the manifold ∂Ω. For this, we recall the result [12, Lemma 10]: Here h X,M is defined as the fill distance h X , but using the intristic distance on the manifold instead of the Euclidean one, and W t p (M) is a Sobolev space on the manifold M.
We remark that similar statements as in Theorem 1 and Theorem 2 are also available for interpolation with analytic kernels such as the Gaussian kernel. However, as we focus on Sobolev kernels (see also Remark 11), we do not recall those estimates here.

Greedy kernel interpolation
As presented in the introduction, the literature introduces the four greedy kernel algorithms P -greedy, f ·P -greedy, f -greedy and f /P -greedy, which were recently unified in the framework of β-greedy algorithms [39]. Those β-greedy algorithms for a fixed β ∈ [0, ∞) select the next interpolation point according to such that one obtains the standard algorithms for β = 0, 1/2, 1 and the f /Pgreedy as a corresponding limit case β → ∞.
In [28] the convergence rate of the P -greedy algorithm was analyzed. In the case of a Sobolev kernel of smoothness τ , one obtains a decay of the Power function as P n L ∞ (Ω) ≤ Cn 1/2−τ /d , which directly translates into an error bound via the standard power function bound from Eq. (3). For the target-data dependent greedy kernel algorithms, i.e. those which rely on the data of the target function values f , i.e. β > 0, we recall the case of finitely smooth kernels from [39,Corollary 12] for corresponding error bounds: be a strictly positive definite translational invariant kernel of finite smoothness τ > d/2 with associated reproducing kernel Hilbert space H k (Ω), where the domain Ω ⊂ R d is assumed to be bounded with a Lipschitz boundary and to satisfy an interior cone condition. Then any β-greedy algorithm with β ∈ {0, 1/2, 1} applied to some function f ∈ H k (Ω) satisfies the following error bound for n ∈ N 0 and r i ≡ f − s i : We refer the reader to [39,Corollary 12] for the exact values of the constants, and the result in its full generality can be found in [39,Corollary 11]. Theorem 3 especially shows that the use of a target-data dependent selection criterion yields a theoretically faster convergence rate than the target-data independent P -greedy criterion, as is often observed in practice.

Generalized kernel interpolation
As already mentioned in the introduction in Section 1, kernels can also be used for interpolation of functionals and we point the reader to [37,Chapter 16] for more detailed background information. Let Λ n = {λ 1 , ..., λ n } ∈ H k (Ω) ′ be n linearly independent functionals with associated Riesz representers v λ1 , ..., v λn ∈ H k (Ω) and V n = span{v λi |λ i ∈ Λ n }. Then the unique minimal-norm generalized interpolant can be obtained as in Eq. (5), which is the projection of the target function f onto the subspace V n spanned by those Riesz representers.
In order to bound the approximation error, a generalized power function similar to the standard power function in Eq. (2) can be defined as which immediately gives for all u ∈ H k (Ω) We like to point out again that for λ 1 := δ x1 , ..., λ n := δ xn the generalized approximant as well as the generalized power function boil down to their standard counterparts by identifying x and δ x . Similar to the above task of standard kernel interpolation, the search for suitable functionals is a challenging question. Again greedy methods are a suitable approach and first steps to the analysis of the corresponding version of the P -greedy algorithm for generalized interpolation were done in [30]. However the target data dependent greedy algorithms for generalized interpolation were not investigated so far in the literature and we aim at advancing in this direction in this paper.
Independently of the way the next interpolation points / functionals are chosen, it makes sense not to use the standard basis given by the Riesz representers v λ1 , ..., v λn for greedy approximation, because otherwise the coefficients α j within the standard kernel expansion from Eq. (1) respective Eq. (5) need to be recomputed after each basis extension step. Instead one usually makes use of a Newton basis [25,30], which is derived iteratively by applying a Gram-Schmidt orthonormalization process to the selected Riesz representers of the corresponding functionals. This allows for an incremental update procedure without the need to recompute the coefficients.
Let v 1 , ..., v n be the resulting Newton basis functions derived from the Riesz representers v λ1 , ..., v λn . Then the (standard or generalized) interpolant s n to u ∈ H k (Ω) can be written as and it can be shown that the Newton basis coefficients satisfy Therefore, it also holds that where the inequality is in fact an equality iff s n n→∞ −→ u in H k (Ω).

Solving PDEs by collocation
As motivated in the introduction in Section 1, the choice of Λ according to Eq. (7) naturally yields an approximation of PDEs in the setting of generalized interpolation. This approach is usually called symmetric kernel collocation since it is obtained by a symmetric application of the differential operators to the two arguments of the kernel or, in other words, since the differential operators are applied both to construct the ansatz in Eq. (5) (or (9)), and to apply the generalized interpolation conditions of Eq. (4). Since this method is based on the optimality principles described in Section 1, it can be analyzed by extending in a systematic manner the ideas used for plain interpolation.
If the approximant s n in Eq. (5) is instead assumed to be represented as s n := n j=1 α j k(·, x j ) independently of the specific problem to be solved, the resulting method is usually referred as unsymmetric kernel collocation, or method of Kansa (see [17,18]). This approach has the advantage of requiring less regularity on the kernel and of being usually computationally more flexible, but the corresponding theoretical analysis is more subtle [16,29].
These ideas have been significantly extended by introducing a local approximation scheme known as RBF Finite Difference (RBF-FD) [35], which is becoming the state of the art for the solution of PDEs by kernel collocation (see e.g. [8]). The convergence of this method has been extensively studied only very recently [36].

Analysis of the Kolmogorov widths for PDEfunctionals
For the analysis of the greedy kernel approximation algorithms in Section 4 we will extend an analysis introduced in [28,39]. For this, knowledge about the Kolmogorov n-widths [26] that are associated to our approximation problems are required, whereby we focus on Dirichlet boundary conditions for simplicity.
The shorthand notation d n (Λ) will be used whenever the Hilbert space H respectively its dual H ′ are clear. We remark that we can rewrite this Kolmogorov width: It is possible to move this definition from the dual space H ′ to the primal space H, since the Riesz representer Theorem gives an isometry. Furthermore we calculate the distance in the Hilbert space by using orthogonal projections and obtain that For the problem of PDE approximation the Hilbert space is chosen as H := H k (Ω), i.e. as the RKHS of some kernel k. The set Λ is defined via the PDE equations in Eq. (7), i.e.
With x λ ∈ Ω ∪ ∂Ω we will denote in the following the point related to a functional λ ∈ Λ ⊂ H k (Ω) ′ , i.e. where the evaluation takes place within Ω (for evaluations of λ ∈ Λ L ) or on its boundary ∂Ω (for evaluations of the boundary functionals λ ∈ Λ B ). Since we are focusing on Sobolev kernels, we have the norm equivalence H k (Ω) ≍ H τ (Ω) for some τ > d/2. Therefore we immediately obtain i.e. instead of analyzing the Kolmogorov n-widths related to the RKHS H k (Ω) we can analyze the n-widths related to the Sobolev space H τ (Ω), because they provide the same asymptotics. As remarked in [30], a thorough analysis of those Kolmogorov n-widths for PDEs does not seem to be available so far in the literature. Nevertheless [30] conjectured that there should be a τ − (smoothness of the Hilbert space) and Indeed it will be shown in the following that the Kolmogorov n-widths related to our PDE problems exactly behave according to this relation. We stress that the difficulty in obtaining bounds for d n (Λ) is due to the fact that Λ is the union of two different sets of functionals.
Remark 4 (Relation to Kolmogorov n widths for parametric PDEs). We want to briefly relate and discriminate the present notion of Kolmogorov widths for PDEs to other usage in literature. Especially in the context of parametric PDEs and projection-based model order reduction, solution manifolds of the parametric solutions appear and need to be approximated, e.g. [15]. A crucial property for successful approximation in these scenarios is a rapid decay of the Kolmogorov widths of those solution manifolds. Especially in transport-dominated cases, the decay can be very slow for the wave equation [14] or the convection of discontinuities [24]. In contrast, we here are not aiming at approximating solution manifolds of parametric PDE solutions, but addressing the Kolmogorov width of the set of the PDE collocation functionals of a single non-parametric PDE.
The next Subsection 3.1 provides results to obtain lower and upper bounds on d n ( M j=1 Λ j ) with help of d n (Λ j ), j = 1, . . . , M . Subsequently we derive in Subsection 3.2 precise decay rates for our special case of Sobolev spaces.

Abstract setting: Bounds on d n (Λ)
We call this subsection the abstract setting, because the following Theorem 5 and Theorem 6 provide bounds on the Kolmogorov n-widths of general sets that are given as the union of two or more (not necessarily disjoint) sets. Moreover, although we use the same notation used before for orthogonality, it should be noted that Theorem 5 and Theorem 6 work also for Kolmogorov n-widths in Banach spaces because we only use distance-based arguments and no dot products.
The following Theorem 5 shows that it is possible to lower bound the Kolmogorov n-width d n (Λ) of Λ = M j=1 Λ j with help of the Kolmogorov n-widths d n (Λ j ), j = 1, .., M .
Proof. For any j = 1, . . . , M we can calculate: Therefore the statement directly follows.
In a similar way it is possible to derive an upper bound on d n (Λ) via the quantities d n (Λ j ), j = 1, . . . , M .
Proof. We calculate: Theorem 5 and 6 enable us to analyze the asymptotics of d n (Λ) with help of the knowledge of d n (Λ j ), 1 ≤ j ≤ n. This will be applied in the following subsection for the particular example of Sobolev kernels.

Kernel setting: Bounds on d n (Λ)
We start by analyzing the Kolmogorov n-widths d n (Λ L ) and d n (Λ B ) as defined in Eq. (7). To the best of our knowledge there are no similar estimates on d n (Λ L ) and d n (Λ B ) available in the literature.

Upper bound on
For d n (Λ L ) we obtain the following theorem.
where C L > 0 is a constant depending on L but not on n.
Proof. Based on Eq. (21) and Eq. (17) we have the following representation of the Kolmogorov n-width d n (Λ L ) (using the Riesz representer v µ of a functional µ): For the third equality Λ L = {δ x • L | x ∈ Ω} was used and in the final step the infimum was upper bounded by the infimum computed over spaces spanned by any n elements of Λ L . Since u, Π Hn (u) ∈ H k (Ω) ≍ H τ (Ω) and the differential operator L is of degree deg(L) = m, it follows that Lu, LΠ Hn (u) ∈ H τ −m (Ω). Recall that where the evaluation takes place within Ω (for evaluations of λ ∈ Λ L ) or on its boundary ∂Ω (for evaluations of the boundary functionals λ ∈ Λ B ). Then To bound the infimum over all H n ⊂ H k (Ω) that are induced by Λ L -functionals, one can now consider functionals related to asymptotically uniformly distributed points. For such a choice of points the fill distance h Xn decays according to h Xn ≤ cn −1/d , such that we obtain This result can be plugged into the (in)equality chain above for the Kolmogorov n-width d n (Λ L ), such that we obtain Hereby Lemma 21 from the Appendix B was leveraged for the second inequality.
where C B > 0 is a constant depending on B = Id but not on n.
Proof. The proof follows the same strategy as the one of Theorem 7: Based on Eq. (21) and Eq. (17) we have analogously to the calculation in Eq. (22) the following representation of the Kolmogorov n-width d n (Λ B ) (using the Riesz representer v µ of a functional µ): Since the functionals λ 1 , . . . , λ n are associated to the interpolation points x λ1 , ..., x λn , u − Π Hn (u) is a function with zeros in these interpolation points. Furthermore note that u − Π Hn (u) ∈ H k (Ω) ≍ H τ (Ω) ֒→ H τ −1/2 (∂Ω) by Lemma 23. Therefore the sampling inequalities from Theorem 2 can be applied to the smooth manifold M = ∂Ω, provided h Xn,M ≤ C M (using k = d − 1, µ = 0, q = ∞, p = 2, t = τ −1/2 and u−Π Hn (u) instead of u, which suffice the assumptions) where . This result can be plugged into the (in-)equality chain for the Kolmogorov nwidth d n (Λ B ), such that we obtain We remark that 1 for m ≥ 0, i.e. the decay in Theorem 8 is faster than the decay in Theorem 7.

Bounds on d n (Λ)
This subsection leverages the abstract results from Subsection 3.1 and combines them with the concrete estimates for the kernel setting from Subsection 3. Proof. For this proof we introduce α : For n ≥ 2 we have 0 < 1−2 α/β−1 ≤ 1−n α/β−1 ≤ 1, therefore we have n 1 (n) ≍ n. In more details, there are constants c 1 , c 2 > 0 such that c 1 n ≤ n 1 (n) ≤ c 2 n for n ≥ 2, c 1 , c 2 are independent of n, and c 1 < 1. Thus, using Theorem 6, 7 and 8 gives A tiny extra calculation helps to resolve the expressions n 1 (n) −α and n 2 (n) −β : Therefore we finally obtain Recalling the abbreviation −α = 1 2 − τ −m d we obtain the statement.
Applied to the PDE setting of [30] this means that the Kolmogorov width related to the solution of the PDE Lu = f (here given by d n (Λ L ) ≍ n −α ) decides about the Kolmogorov widths for the whole PDE (here d n (Λ)), as it was already conjectured in [30]. The additional boundary conditions do not influence the decay of the Kolmogorov n-width, because the Kolmogorov n-widths related to the boundary conditions decay faster.
This observation can also be verified in the numerical experiments related to the PDE-P -greedy algorithm in Subsection 6.2.
The proposed analysis in this section is in principle not restricted to only two sets of functionals Λ L and Λ B . This is explained in the following remark.
Remark 10. We remark that the results of this section can be extended to more sets of functionals. This was already done for the abstract setting in Subsection 3.1 by considering Λ = M j=1 Λ j . In the kernel setting this refers e.g. to more types of boundary conditions, i.e. we can incorporate both Dirichlet boundary conditions and Neumann boundary conditions by using Λ B1 (for the Dirichletrelated functionals) and Λ B2 (for the Neumann-related functionals) instead of only using Λ B .
Furthermore we have another remark on the applicability of the proposed analysis for analytic kernels such as the Gaussian.
Remark 11. The same analysis and used techniques are also applicable to analytic kernels. In particular the results of Subsection 3.1 still hold true, since they were completely independent of kernels. Also the following analysis of Section 4 works, because it holds for general kernels without any smoothness asumption. However the analysis of Subsection 3.2 needs to be modified, as sampling inequalities for analytic functions need to be leveraged. We leave the details to future work.

Target-data dependent analysis of PDE-greedy
First we define what a PDE-greedy kernel algorithm is:

Decay of power function quantities
In [39, Section 3] the abstract analysis of greedy algorithms in Hilbert spaces of [3] was extended by considering a broader class of greedy algorithms. Here we will make use of those results and apply them to obtain bounds on see Corollary 14, with Λ i = {λ 1 , ..., λ i } as introduced in Subsection 2.3. This will later enable us to derive convergence rates for a range of target data dependent PDE-greedy kernel algorithms. For this we recall the notation of [39, Section 3]: Let H be a Hilbert space with norm · = · H and F ⊂ H a compact subset. Without loss of generality we assume that f H ≤ 1 for all f ∈ F . We consider algorithms that select elements f 0 , f 1 , · · · ∈ F , without yet specifying any particular selection criterion. We define V n := span{f 0 , . . . , f n−1 } and the following quantities, where H n is any n-dimensional subspace of H: Now we recall [39, Corollary 2] including its assumptions which were stated before in [39, Theorem 1]): Corollary 13 (Corollary 2 from [39]). Consider a compact set F in a Hilbert space H, and a greedy algorithm that selects elements from F according to any arbitrary selection rule.
These choices allow us to formulate the abstract definitions of Eq. (24) in the kernel setting for PDE approximation. For this we point again to the definition of the generalized power function from Eq. (17). The quantity σ n is then simply the maximal value sup λ∈Λ P Λn (λ) of the generalized power function P Λn , while ν n is the value of the generalized power function P Λn at the next selected functional λ n+1 corresponding to the point x λn+1 .
Furthermore, d n can be similarly bounded as Due to Eq. (26), any bound on sup λ∈ΛL∪ΛB P Λn (λ) for a given set of functionals Λ n ⊂ F (which corresponds to a distribution of points in Ω and on the boundary ∂Ω) gives a bound on d n . Additionally, observe that the assumption f H ≤ 1 for f ∈ F from the abstract setting is satisfied using the kernel setting as soon as which can always be ensured by rescaling the kernel k to ǫk, 0 < ǫ < 1 (which does not change the RKHS, but only rescales its norm). Applying this established link to Corollary 13 immediately gives the following result: , and a greedy algorithm that selects elements from F according to any arbitrary selection rule.
0 e α log(n) α n −α , n ≥ 3, Proof. This is a direct consequence of Corollary 13 using the choices of H, F and f 0 , f 1 , ..., f i as introduced before and summarized in Eq. (25).
The assumptions of Corollary 14 about the decay of the Kolmogorov widths d n (F ) can be ensured via Theorem 9 because of the choice F ≡ {v λ , λ ∈ Λ} and the Riesz representer isometry: Again we note that the results of this subsection simplify to the results of [39,Section 5] if we choose L = Id, i.e. {v λ | λ ∈ Λ 1 ∪ Λ 2 } = {k(·, x), x ∈ Ω}.

Definition of PDE-β-greedy algorithms
We generalize the notion of β-greedy algorithms as introduced in [39] for standard interpolation to PDE-β-greedy algorithms via the following definition.

Analysis of β-greedy algorithms
We remark that the following lemma and theorems are generalizations of the corresponding statements within [39], when we pick Λ = {δ x | x ∈ Ω} or equivalently F = {k(·, x) | x ∈ Ω}. Consequently also the proofs for the theorems within this section are generalizations of the corresponding proofs in [39,Section 4], made possible by using properties of generalized interpolation instead of standard interpolation. Therefore they are not included in the main text, but can be found in the Appendix A.

Convergence rates for PDE-greedy algorithms
This section serves as an amalgamation of the results from the previous sections : We apply the analysis of the PDE-β-greedy algorithms from Section 4 to concrete kernels and PDEs. By using estimates on the corresponding Kolmogorov widths from Section 3 we obtain convergence rates for PDE-greedy algorithms. Using the maximum principle for elliptic operators (see Eq. (13)) we can transfer the obtained convergence rates directly to convergence rates of the approximation of the true solution.
We make use of the following assumptions, which especially give rise to a well-posed PDE: Let Ω ⊂ R d be a bounded domain with smooth boundary ∂Ω of dimension d − 1. Let L be an elliptic operator of second order with smooth coefficients (i.e., deg(L) = 2) and B = Id, i.e. the Dirichlet boundary operator. Let k be a Sobolev kernel of smoothness τ > d/2 + 2, i.e. H k (Ω) ≍ H τ (Ω). Also let f ∈ H τ −2 (Ω) and g ∈ H τ −1/2 (∂Ω), and assume there exists a solution u ∈ H τ (Ω).
Theorem 19. We consider the PDE as in Eq. (10) under the Assumption 1. Then any β-greedy algorithm with β ∈ [0, ∞] applied to this PDE satisfies the following error bound for n = 1, 2, . . . (using 1/∞ = 0): with C := 2 α+1/2 max{1, C 0 }e α 1 max{1,β} and α = 1 Proof. The proof is structured into several steps: • Since k is a kernel of finite smoothness such that H k (Ω) ≍ H τ (Ω), Theorem 9 gives a decay of the corresponding Kolmogorov width d n (Λ) of This bound satisfies the assumptions of Corollary 14, such that it can be applied to derive a decay of the geometric mean of subsequent power function values as for any choice of functionals, thus especially for the β-greedy choices.
• Now we make use of the maximum principle from Eq. (13) which is available due to our assumptions on the PDE. Therefore it holds Since r j = u−s j we can directly combine this last equation with the result of the previous step to obtain min j=n+1,..,2n For β ∈ (0, 1] it can be observed that we obtain a faster rate of convergence in the number of selected functionals λ i ∈ Λ L ∪ Λ B than for β = 0 and therefore in the number of selected interpolation points x i ∈ Ω ∪ ∂Ω. This additional decay of n −β/2 increases with increasing β ∈ (0, 1], which corresponds to a more target data dependent choice of the next interpolation point according to the definition of PDE-β-greedy algorithms, see Definition 15. In particular the proven decay for PDE-f -greedy (β = 1) is better than the proven decay for the PDE-Pgreedy. This makes intuitively sense, since the PDE-f -greedy choice is adapted to the right hand side of the PDE, while the PDE-P -greedy is independent of the right hand side of the PDE and only worst-case optimal [30]. It is furthermore optimal in the sense of realizing the same decay rate as the Kolmogorov n-width (which is based on the abstract result of [3]). A visualization of the different point distributions can be found in Subsection 6.2.
Moreover, the additional n −β/2 convergence breaks the curse of dimensionality because it does not depend on the dimension d of the space, in contrast to (log(n) · n −1 ) α/ max{1,β} because α = 1 2 − τ −2 d deteriorates in the dimension d. This breaking of the curse of dimensionality is theoretically and practically appealing for high dimensional PDE problems.
As already remarked in [39], we conjecture that the analysis so far is suboptimal. In particular we do not expect that the additional log(n) α factor is really required, which occurs for β > 0 -however we did not find a way to circumvent it or show an example where it is necessary. Similar statements as in Theorem 19 can be derived for the case of analytic kernels such as the Gaussian kernel. We did not include such results because of the following reasons: 1. For analytic kernels such as the Gaussian kernel, spectral convergence can be expected. Therefore the additional convergence rate of up to n −1/2 which can be obtained for target data dependent algorithms using our analysis does not seem that appealing as it is way inferior to exp(−n 1/d ).
Nevertheless also in this case, the additional convergence rate of up to n −1/2 is interesting because it still breaks the curse of dimensionality in contrast to exp(−n 1/d ). A comparison and connection with the dimension independent convergence rates within [7] for the Gaussian kernel will be part of future work.
2. The RKHSs of analytic kernels such as the Gaussian are difficult to describe, as they cannot simply be characterized with help of an equivalence to a Sobolev space [33]. In particular these RKHS are usually subsets of Sobolev spaces of arbitrary orders. In order to derive statements for a suitable large class of functions, we need to apply results for functions outside the RKHS, which is usually treated under the notion of escaping the native space. This will be part of future work.
Finally we remark that the first step of the proof of Theorem 19 does neither require the uniqueness of a solution nor the existence of a solution. In fact the PDE-β-greedy algorithms can be applied to any kind of linear PDE, i.e. also to non well-posed PDEs: Also in this case we obtain the convergence of the geometric mean of subsequent power function values, since it depends mostly on the kernel k and the degree deg(L) of the PDE, but not on the exact form of the differential operator L or the right hand side functions. However this convergence cannot be related to a convergence of the corresponding interpolant s n . This is elaborated in a numerical example in Subsection 6.4.
Remark 20. We want to remark that the whole analysis also works for PDEs on manifolds, e.g. Theorem 18 holds for both Ω a domain or a manifold. The analysis of target-data dependent greedy kernel algorithms works in very general terms, because it simply leveraged the abstract analysis of [3] and [39,Section 3] for Hilbert spaces and connects it to kernel approximation in RKHS. In order to derive practically usable convergence rates, the Kolmogorov n-widths need to be analyzed (Section 3). In order to derive those for manifolds, sampling inequalities and corresponding zero lemmas for manifolds can be leveraged, as e.g. in Theorem 2. In order to transfer the resulting convergence rates on sup λ∈Λ |λ(r i )| to convergence rates on r i → 0 in some norm, PDE theory is required, e.g. some well-posedness estimate or a maximum principle.

.1 Convergence rates of the power function
This subsection considers the PDE-P -greedy algorithm using four different kernels on the domain Ω = (0, 1) 2 for L = −∆ in conjunction with Dirichlet boundary values. The domain Ω was discretized using 10 4 uniformly randomly sampled points, while the boundary ∂Ω was discretized using 10 3 uniformly randomly sampled points. The power function values on the interior and the boundary were adjusted (importance reweighting): Instead of P Λn (λ) for λ ∈ Λ 1 we considered P Λn (λ)/(L (1) L (2) k)(0, 0), which was necessary because it holds (L (1) L (2) k)(x, x) ≫ k(x, x) for the kernels we considered. The results are visualized in Figure 2. The number of selected interior and boundary points are noted in Table 1. The black lines within these plots visualize the theoretically expected rate of decay according to (for finitely smooth kernels) The corresponding values for τ can be found e.g. in [37].
In Figure 2 some decay rates seem to be slightly too fast (i.e. faster than the black dashed lines). The explanation for this might be that only 10 4 points were used for the discretization of the interior, which is not significantly more than the 10 3 points which were selected.

Comparison of PDE-f -greedy and P -greedy
In this subsection we consider the PDE −∆u = 0, on Ω, on a pacman shaped domain as given in Figure 3. The unique solution u to this PDE is given by u(x) = x 1 x 2 + 1. The cubic Matern kernel was used for approximation and 200 centers were chosen in each algorithm from an overall set of 52000 points (50000 points in Ω and 2000 points in ∂Ω, both uniformly randomly sampled). Figure 3 visualizes the point distribution of the selected PDE-f -greedy and PDE-P -greedy points within the domain and on its boundary. It is clearly visible that 1. the PDE-P -greedy selected points are quasi uniformly distributed. Way more points were selected in the interior than on the boundary.
2. the PDE-f -greedy selected points cluster towards some parts of the boundary. Also way more points are selected in the interior than on the boundary.
A surprising observation can be made: The points selected by the target-data dependent PDE-f -greedy algorithms do not cluster around the cone at zero, but next to the remaining parts of the boundary. This seems in contrast to typical results of e.g. adaptive finite element discretization results, where grid adaptivity, hence mesh-point-clustering is expected towards the corner singularity at the origin. The reason for this behaviour is likely that the right hand side function and solution x 1 x 2 + 1 is very benign around zero. The errors reached by PDE-f -greedy and PDE-P -greedy are given by which clearly shows that the target-data adaptive choice of PDE-f -greedy gives a more accurate approximation of the true solution u. We remark that it holds u L ∞ (Ω) ≈ 1, i.e. the residual values reached by PDE-P -greedy and PDE-fgreedy are several orders of magnitude smaller. Furthermore we like to remark that the set Ω used here does not have a smooth boundary ∂Ω, but still the numerical results match the analytical results in the case of a smooth boundary. A crucial point -the existence of a solution -is satisfied nevertheless because it is explicitly given. We apply PDE-β-greedy algorithms using the quadratic Matérn kernel for β ∈ {0, 0.2, 0.4, 0.6, 0.8, 1} to this problem. Since there are only two boundary points ∂Ω = {0, 1}, this example mostly focusses on the interior points.

Scale of PDE-β-greedy algorithms
The results are depicted in Figure 4, which shows both the decay of sup λ∈Λ |λ(u − s β n )| and u − s β n L ∞ (Ω) in the number n of greedily chosen functionals. It is clearly visible that the decay of sup λ∈Λ |λ(u − s β n )| has the fastest decay for PDE-f -greedy, while PDE-P -greedy provides the slowest decay. For intermediate values of β ∈ (0, 1), intermediate convergence rates can be observed.
Surprisingly these different convergence rates do not directly transfer to the decay of u − s β n L ∞ (Ω) : Here we can observe the same (maximal) rate of convergence n −2 independent of the β value. The advantages of target-data dependent PDE-greedy algorithms, i.e. β > 0 is that they provide faster convergence rates in stronger norms, e.g. in the W 2 ∞ (Ω) (and thus W 2 2 (Ω)) norm: • Because of the exact fulfillment of the Dirichlet boundary conditions (after picking the two Λ B functionals) the maximums principle gives • Applying the Poincaré-Wirtinger inequality to ∂ x (u − s β n ) (which satisfies This together implies and thus a faster decay of sup λ∈Λ |λ(u − s β n )| translates into a faster convergence in the W 2 ∞ (Ω) norm. This can also be seen in the top of Figure 4. However we remark that this phenomenon described here seems to be unique for dimension d = 1 and we did not observe it in d ≥ 2.
Nevertheless we can apply the PDE-f -greedy to this problem by choosing L = −∆−µ 1 ·Id. For this we used the Wendland kernel Φ 3,3 and a discretization of Ω into 1001 equidistant points.
The result is depicted in Figure 5. We can observe that we obtain an interpolation of the right-hand side f (x) = sin(πx) in the interpolation points, however apparently no convergence on the full domain. By computing the dot product f, Ls n L 2 (Ω) numerically we obtain values very close to zero, which is in accordance with R(−∆ − µ 1 Id) ∋ Ls n ⊥ f ≡ sin(πx). The reason for this non-convergence is the non-existence of some solution u ∈ H k (Ω): This was the basic assumption on which the whole theory was built (see beginning of Section 2 or Assumption 1).

Conclusion and outlook
This paper considered the approximation of solutions of linear PDEs via symmetric kernel collocation with help of greedy kernel methods. After deriving estimates on the Kolmogorov n-widths related to those PDEs, an abstract analysis of greedy kernel methods was applied to derive worst-case optimal approximation rates. Extending an analysis of target-data dependent greedy kernel algorithms for standard interpolation to symmetric kernel collocation, it was possible to derive better convergence results depending on the degree of targetdata dependency. These algorithms yield point sets which are adapted to both the domain Ω and the right hand side of the considered PDE. Open questions include the optimality of the proven decay rates as well as approximation in different norms and work on analytic kernels such as the Gaussian kernel. Future work will also address the analysis of greedy algorithms for other methods of solving PDEs via kernel methods, such as unsymmetric kernel collocation or RBF-FD methods.  As expected there is obviously no convergence. Right: Visualization of the right hand side f (black) and well as the prediciton Ls 20 (pink) after 20 steps over the domain Ω = (0, 1). Although we interpolate at several points, we do not have an overall approximation.
For β ∈ (0, 1], we start by reorganizing the estimate (31) of Lemma 17 to get Proof. This is a consequence from [20, Chapter 1: Lemma 7.2+Remark 8.1] and the triangle inequality. Proof. This is a consequence of the previous Theorem 22.