{'fraction_non_alphanumeric': 0.060324942332765016, 'fraction_numerical': 0.02244007622104102, 'mean_word_length': 3.7521744310735134, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sparse Gaussian processes for solving nonlinear PDEs

In this article, we propose a numerical method based on sparse Gaussian processes (SGPs) to solve nonlinear partial differential equations (PDEs). The SGP algorithm is based on a Gaussian process (GP) method, which approximates the solution of a PDE with the maximum a posteriori probability estimator of a GP conditioned on the PDE evaluated at a finite number of sample points. The main bottleneck of the GP method lies in the inversion of a covariance matrix, whose cost grows cubically with respect to the size of samples. To improve the scalability of the GP method while retaining desirable accuracy, we draw inspiration from SGP approximations, where inducing points are introduced to summarize the information of samples. More precisely, our SGP method uses a Gaussian prior associated with a low-rank kernel generated by inducing points randomly selected from samples. In the SGP method, the size of the matrix to be inverted is proportional to the number of inducing points, which is much less than the size of the samples. The numerical experiments show that the SGP method using less than half of the uniform samples as inducing points achieves comparable accuracy to the GP method using the same number of uniform samples, which significantly reduces the computational cost. We give the existence proof for the approximation to the solution of a PDE and provide rigorous error analysis.


Introduction
This work develops a new numerical method based on sparse Gaussian processes (SGPs) [14,31] to solve nonlinear partial differential equations (PDEs). PDEs have been widely used to model applications in science, economics, and biology [20,30]. However, very few PDEs admit explicit solutions. Standard numerical approaches including finite difference [30] and finite element methods [9] for solving PDEs are prone to the curse of dimensions. Recently, to keep pace with increasing problem sizes, machine learning methods intrigue a lot of attention [1,15,16,24,36]. Probabilistic algorithms to solve linear problems are presented in [2,17,18]. The authors of [1] extend these ideas and propose a Gaussian process (GP) framework for solving general nonlinear PDEs. In particular, for time-dependent PDEs, GP methods based on time-discretization are considered in [11,23,32]. The computational costs of the above GP methods grow cubically with respect to the number of samples due to the inversion of covariance matrices. To improve the performance, based on [1], [16] proposes to approximate the solution of a PDE in the space generated by random Fourier features [22,35], where the corresponding covariance matrix is a low-rank approximation to that of the GP method in [1].
In this paper, we propose another low-rank approximation method based on SGPs, which leverages a few inducing points to summarize the information from observations. Our SGP method is based on [1], where a numerical approximation of the solution to a nonlinear PDE is viewed as the maximum a posteriori (MAP) estimator of a GP conditioned on solving the PDE at a finite set of sample points. In our SGP method, we replace the Gaussian prior by an alternative with a zero mean and a low-rank kernel induced from inducing points. The approximation to the solution of a PDE can be viewed as a MAP estimator of a SGP conditioned on a noisy observation about the values of linear operators of the true solution, where the observation satisfies the PDE at the samples. Meanwhile, the numerical approximation can also be viewed as a function with the minimum norm in the reproducing kernel Hilbert space (RKHS) associated with the low-rank kernel. The values of linear operators of the function are close to those satisfying the PDE at sample points. By a stable Woodbury matrix identity stated in Subsection 2.2, the dimension of the matrix to be inverted in the SGP method depends only on the number of inducing points selected.
This inducing-point approach is initially proposed in the deterministic training conditional (DTC) approximation [26], and then it has been fully studied in several works. [21] provides a unifying view for the approximation of GP regressions. This framework includes the DTC, the fully independent training conditional (FITC) approximation [29], and the partial independent training conditional (PITC) approximation. All above methods modify the joint prior using a conditional independent assumption between training and test data given inducing variables. Alternatively, [8,31] leverage the inducing-point approach into the computing of the evidence lower bound of the log marginal likelihood. They retain exact priors but approximate the posteriors via variational inference. On the other hand, the estimation of inducing inputs has been generalized into an augmented feature space in [12]. In particular, SGPs are extended in the spectral domain with and without variational inference [7,13]. Moreover, some works direct speed up the computation of GP regressions through fast matrix-vector multiplication and pre-conditional conjugate gradients [5], and through structured kernel interpolation [34].
Our contributions are summarized as follows: 1. We introduce a general framework based on SGPs to solve nonlinear PDEs. The numerical experiments in Section 4 show that the SGP method using less than half of the uniform samples as inducing points significantly reduces the computational cost without losing much accuracy compared to the GP method in [1] using the same number of uniform samples; 2. We give a rigorous existence proof for the numerical approximation to the solution of a general PDE and analyze error bounds in Section 3. Our analysis bases on the arguments in [1], learning theory [27,28], and Nyström approximations [10]. The error bound given in Theorem 3.7 provides qualitative hints to improve the accuracy of our method; 3. The numerical experiments show that the choice of the positions of inducing points and parameters of kernels has a profound impact on the accuracy of the numerical approximations. In Section 2, we give a probabilistic interpretation for the SGP method in the setting of [31], which motivates future work for hyperparameter learning.
The paper is organized as follows. Section 2 summarizes our SGP method by solving a nonlinear elliptic PDE. In Section 3, we provide a framework of the SGP method to solve general PDEs. Numerical experiments are presented in Section 4. Further discussions and future work appear in Section 5.
Notations. Vectors are in bold font. For a real-valued vector v, we represent by |v| the Euclidean norm of v and by v T its transpose. For a matrix A, we denote by A the spectral norm of A. Let u, v or u T v be the inner product of vectors u and v. For a normed vector space V , let · V be the norm of V . Let U be a Banach space associated with a quadratic norm · U and let U * be the dual of U. Denote by [·, ·] the duality pairing between U * and U. We assume that there exists a linear, bijective, symmetric Let {φ i } P i=1 be a set of P ∈ N elements in U * and let φ := (φ 1 , . . . , φ P ) be in the product space (U * ) P . Then, for u ∈ U, we denote the pairing [φ, u] by Furthermore, for u := (u 1 , . . . , u S ) ∈ U S , S ∈ N, we represent by [φ, u] ∈ R P ×S the matrix with entries [φ i , u j ]. We write · as the induced norm for linear operators on U, i.e., for any F : U → U, we define F = max u U =1 F(u) U . Given any R ∈ N, we denote by [R] the set of nonnegative integers less equal to R. Finally, we represent by C a positive real number whose value may change line by line.

A Recapitulation of the SGP Method
In Subsection 2.1, we demonstrate the SGP method by dealing with a nonlinear elliptic PDE given in [1]. The SGP method approximates the solution of the PDE by a minimizer of an optimal recovery problem, which searching a solution with the minimum norm in the RKHS associated with a low-rank kernel generated by inducing points and the values of linear operators of the solution are close to those satisfying the PDE at finite samples. From a probabilistic perspective, the optimal recovery problem can also be interpreted as the MAP estimation for a SGP constrained by a noisy observation of values of linear operators of the true solution. The main bottleneck of the SGP method is the inversion of a covariance matrix. In Subsection 2.2, we describe a robust way to compute the inverse of the covariance matrix, which is the cornerstone of the SGP method to reduce computational complexity. Finally, in Subsection 2.3, we give a probabilistic explanation for the SGP method, which motivates a way to choose hyperparameters in future work. Given continuous functions f : Ω → R, τ : Ω → R, and g : ∂Ω → R, we find the strong solution u * solving where τ is chosen such that (2.1) admits a unique strong solution. We propose to approximate the solution to (2.1) by functions in a RKHS generated by inducing points. First, we show the construction of such RKHS. We assume that the solution of (2.1) is in the space U, where U is the RKHS associated with a nondegenerate, symmetric, and positive definite kernel K : Ω × Ω → R. We view U as the ambient space. Let U * be the dual space of U and let [·, ·] : U * × U → R be the duality pairing. Next, we sample M inducing Then, we define the functionals φ (1) j = δ xj for 1 j M and φ (2) k = δ x k •∆ for 1 k M Ω . Denote by φ (1) and φ (2) the vectors with entries φ (1) j and φ (2) k , respectively. Let φ be the vector concatenating φ (1) and φ (2) . For simplicity, we denote by φ j the j th component of φ. Define K(·, φ) as the vector with entries K(·, x )φ j (x ) dx and define K(φ, φ) as the matrix with elements Suppose that the inducing points are chosen such that K(φ, φ) are invertible. Then, we define a new kernel Q : Ω × Ω → R as Using the positivity of the kernel K, we see that Q is semi-definite (see Lemma 3.1 in Section 3). Thus, by the Moore-Aronszajn theorem, Q induces a unique RKHS, denoted by U Q . We call U Q the RKHS generated by inducing points.
Next, we approximate the solution u in U Q . We take N samples Then, our SGP method approximates the solution u * of (2.1) with a minimizer of the following regularized optimal recovery problem where γ > 0 is a given small regularization parameter and z := (z NΩ ). Let (u † , z † ) be a minimizer of (2.2). A minimizer u † of (2.2) can be viewed as a MAP estimator of a SGP ξ ∼ N (0, Q) conditioned on a noisy observation z † about values of the linear operators of u † , where z † satisfies the PDE at the sample points. A probabilistic interpretation of (2.2) is given in Subsection 2.3.
Then, we define the functionals ψ (1) j = δ xj for 1 j N and ψ ψ as the vector consisting of ψ (1) j and ψ (2) k . We represent by ψ j the j th component of ψ. Let Q(·, ψ) be the vector containing entries Q(·, x )ψ j (x ) dx and let Q(ψ, ψ) be the matrix consisting of elements Q(x, x )ψ i (x)ψ j (x ) dx dx . Theorem 3.4 in Section 3 shows that (2.2) admits a minimizer u † such that where I is the identity matrix and z † is a minimizer of Next, we use the technique of eliminating variables in [1] to remove the constraints of (2.4). More precisely, observing z . Then, we obtain z † by solving (2.5) using the Gauss-Newton method and get u † by (2.3).
Remark 2.1. In the first attempt, we consider a similar formulation to that of the GP method in [1] and solve the optimal recovery problem u(x j ) = g(x j ), ∀j = N Ω + 1, . . . , N.
(2.6) However, since Q(ψ, ψ) = K(ψ, φ)(K(φ, φ)) −1 K(φ, ψ), the matrix Q(ψ, ψ) is not invertible when the number of samples is larger than the number of inducing points. Hence, the representer theorem [19, Section 17.8] does not apply. Furthermore, since U Q is a subspace of U, the solution u * to (2.1) may not lie in U Q . Hence, we may not find a function in U Q satisfying the constraints in (2.6).

2.2.
Computing (γI + Q(ψ, ψ)) −1 . The bottleneck of solving (2.5) is to calculate the inverse of γI + Q(ψ, ψ), whose size is proportional to the number of sample points. The computational complexity of standard Cholesky decomposition algorithms is O(N 3 ) in general. The algorithm in [25] computes approximate Cholesky factors of (γI + Q(ψ, ψ)) −1 in O(N log 2d (N/ )). Using the structure of Q(ψ, ψ), we show that the computational complexity can be reduced and depends only on the number of inducing points. The technique is from Section 3.4.3 of [33] and is well explained in the documents of GPflows * . We provide the full details below for the sake of completeness.
Using the Woodbury identity, we get The equation (2.7) is not numerically stable when γ is small. To get a better conditioned matrix, we perform the Cholesky decomposition on K(φ, φ) and get K(φ, φ) = LL T (see Remark 2.2 below). Then, we consider The identity (2.9) is the cornerstone for the SGP method to reduce the computational complexity. We note that the dimension of the matrix needed to be inverted at the right-hand side of (2.9) depends only on the length of φ, which is proportional to the number of inducing points. Hence, if M << N , using (2.9), we are able to save computational time and memory storage. The numerical results in Section 4 show that a small number of inducing points is sufficient for the SGP approach to achieve comparable accuracy to the GP method in [1].
Remark 2.2. In general, K(φ, φ) is ill-conditioned. Thus, to deal with the inverse of K(φ, φ), we perform the Cholesky decomposition on K(φ, φ) + ηR, where R is a block diagonal nugget built using the approach in [1] and η > 0 is a small regularization parameter.

A Probabilistic Perspective of the SGP Method.
Here, we give a probabilistic interpretation for the SGP method. First, we put a Gaussian prior on the solution u * to (2.1), i.e., let u * ∼ N (0, K), where K is kernel associated with the RKHS U in Subsection 2.1. Let {x i } N i=1 be the set of sample points in Subsection 2.1. Define u = (u * (x 1 ), . . . , u * (x N ), ∆u * (x 1 ), . . . , ∆u * (x NΩ )).
Then, u is a multivariate Gaussian random variable following the distribution N (0, K(ψ, ψ)). Let z be a noisy observation of u, i.e., where is a noise following the distribution N (0, γI), which is independent with u * . Thus, z ∼ N (0, K(ψ, ψ)+ γI) and the probability density of z is Then, solving (2.1) by the GP method in [1] is equivalent to finding z with the MAP, where z satisfies the PDE system. More precisely, the GP method seeks to find z solving (2.10) In our SGP method, we replace the log-likelihood in (2.10) with its lower bound. The derivation follows [31]. (1) , v (2) ). Then, v and z are joint Gaussian random variables. Let q(v) be the variational distribution of v and let q(u, v) = p(u|v)q(v) be the probability density function of the variational distribution of (u, v). Then, by Jensen's inequality, we have (2.11) Thus, F(z, q) provides a lower bound for ln p(z). To get an optimal lower bound, we take the derivative of F w.r.t. q and get Letting δF(q(v))/δq(v) = 0, we get the optimal variational distribution where Z is a real number guaranteeing the unit mass of q * . Then, plugging (2.12) into (2.11), we get Hence, if we replace ln p(z) by F(z, q * ) in (2.10), we note that (2.4) is equivalent to (2.13) Remark 2.3. The probabilistic interpretation of the SGP method above suggests a way for hyperparameter learning. Let θ be the parameter vector consisting of the positions of the inducing points and the kernel's parameters. Write q * θ to emphasize the dependence of q * in (2.12) on θ. Then, the SGP method with hyperparameter learning for solving (2.1) considers the maximization problem (2.14) However, the objective function in (2.14) is highly nonlinear. Thus, the computation of (2.14) is very costly. In the numerical experiments of Section 4, we fix the positions of inducing points. Then, we use the Gauss-Newton method to solve (2.14) for different values of the parameters of the kernel on the uniform grid in the domain of parameters and accept the result that achieves the largest value of the objective function in (2.14). If we choose Gaussian kernels for solving PDEs, which contain only one or two lengthscale parameters, the above grid search method for choosing hyperparameters is efficient. We leave the full optimization over z and θ in (2.14) to future work.

A General Framework for the SGP method
In this section, we provide a framework of the SGP method for solving a general PDE. Then, we give the existence argument for a minimizer of the optimal recovery problem used to approximate the solution of the PDE. Finally, we analyze the approximation errors.

A RKHS Generated by Linear Operators.
In this subsection, we show the construction of a RKHS generated by linear operators based on the abstract theory of RKHSs (see [19]). The RKHS forms a space where we seek approximations for solutions to general PDEs.
Let U be a Banach space, let U * be its dual, and let [·, ·] be their duality paring. We assume further that there exists a covariance operator K : U * → U that is linear, bijective, symmetric ([Kφ, ϕ] = [Kϕ, φ]), and positive ([Kφ, φ] > 0 for φ = 0), and suppose that the norm of U is given by Meanwhile, the inner products of U and U * are given by Then, U coincides with the RKHS space of the kernel K defined by where δ x is the Dirac delta function centered at The GP method proposed in [1] approximates a solution of a nonlinear PDE by a MAP point for the GP ξ ∼ N (0, K) conditioned on PDE constraints at the collocation points. To reduce the computation time, we consider a Gaussian prior with a different covariance operator generated by selected linear operators.
Suppose that Θ is invertible. We define the induced covariance operator Q as The next lemma guarantees that Q is a covariance operator and associates with a RKHS.
be the collection of R linearly independent non-trivial elements of U * . Then, the operator Q defined in (3.2) is linear, symmetric, and nonnegative. Hence, there exists a unique RKHS U Q associated with the kernel Q defined by Proof. It is easy to see that Q is linear and symmetric. Since K is positive, Θ is positive definite. Then, The equality holds if and only if ϕ is perpendicular to φ w.r.t. K. Thus, Q is nonnegative. Hence, the kernel Q in (3.3) is symmetric and positive semi-definite. According to the Moore-Aronszajn theorem, there exists a unique RKHS U Q associated with the kernel Q. Next, we show the relation between U Q and U.
Since Q(x, ·) ∈ U for any x ∈ Ω, the linear span of {Q(x, ·), x ∈ Ω}, denoted by U Q0 , is a subset of U. By the Moore-Aronszajn theorem, U Q is the completion of U Q0 . Thus, where the last two equities follows by (3.2) and (3.3). Thus, we conclude (3.4).
3.2. The SGP method. In this subsection, we present a general framework of the SGP method to solve nonlinear PDEs. Let Ω ⊂ R d be a bounded open domain with the boundary ∂Ω for d 1. Let L 1 , . . . , L D ∈ L(U; C(Ω)) and L 1 , . . . , L D b ∈ L(U; C(∂Ω)) be bounded linear operators for 1 D b D. We seek to find a function u * solving the nonlinear PDE where P, B represent nonlinear operators, and f , g are given data. Throughout this section, we assume that (3.5) admits a unique strong solution in a quadratic Banach space U associated with the covariance operator K, which means that u * has enough regularity for the linear operators to be well defined pointwisely in (3.5).
We propose to approximate the solution to (3.5) in a RKHS generated by linear operators associated with inducing points. To do that, we take M inducing points Let φ (i) be the vector concatenating the functionals φ (i) j for fixed i and define The next corollary gives a RKHS generated by φ.
Corollary 3.2. Let φ be as in (3.6). Assume that the elements of φ are linearly independent in U * . Define Then, there exists a RKHS U Q associated with the kernel Q defined by such that U Q ⊂ U. Denote by · U Q the norm of U Q . Then, Proof. The claim is a direct result of Lemma 3.1.
Next, we propose to approximate the solution u * of (3.5) in the space U Q . More precisely, we take N samples Given a regularization parameter γ > 0, we approximate u * by the minimizer of the following optimization problem For ease of presentation, we rewrite (3.7) into a compact form. To do that, we define the functionals ψ Meanwhile, let ψ (i) be the vector consisting of ψ (i) j for fixed i and define Next, we define the data vector y ∈ R N by Furthermore, we define the nonlinear map Thus, we reformulate (3.7) as (3.10) Next, we prove the existence of minimizers of (3.10) and analyze the error between a minimizer u † of (3.10) and the strong solution u * to (3.5). The arguments are based on the framework in [27]. Let 2 (x) be the set of sequences a = (a x ) x∈x with a, b = x∈x a x b x defining the inner product. Define the sampling operator S x : (3.11) One can view S x as a generalization of the sampling operator in [27], which only samples values of a function at the sample points. Denote by ·, · U Q the inner product of U Q and by S * x the adjoint of S x . Then, for each c ∈ 2 (x) and u ∈ U Q , we have The next lemma gives some basic properties of S x and S * x . For a little abuse of notations, in the rest of this section, we denote by I the identity map or the identity matrix, whose meaning is easy to recognize from the context.  .11) and (3.12). Then, for any γ > 0, S * x S x + γI and S x S * x + γI are bijective. Meanwhile, For any c ∈ 2 (x), we have Furthermore, Proof. We first show that S * x S x + γI is invertible for any γ > 0. Let u, v ∈ U Q and S * x S x u + γu = v. Then, according to (3.11) and (3.12), we have (3.16) Acting ψ on both sides of (3.16), we get Since Q(ψ, ψ) is positive semi-definite, Q(ψ, ψ) + γI is invertible. Thus, We get from (3.16) and (3.17) that Similarly, one can show that S x S * x + γI is bijective. To prove (3.13), we see that For any c ∈ 2 (x), using (3.18) and (3.12), we get which concludes (3.14). For any c ∈ 2 (x), we use (3.14) and obtain Since (3.19) holds for any c ∈ 2 (x), we conclude (3.15).
The following theorem gives the existence of minimizers to (3.10).
Theorem 3.4. Let S x and S * x be defined in (3.11) and (3.12). Given γ > 0, the system (3.10) admits a solution
Proof. In (3.10), given z ∈ R R , we first consider the minimization problem over u By the definition of S x in (3.11), we get (3.23) Taking the function derivative w.r.t. u, we obtain from (3.23) that Thus, using (3.23) and (3.24), we have where the last equality results from (3.15). Thus, (3.25) implies that (3.10) is equivalent to (3.21). Using a similar compactness argument as in the proof of Theorem 1.1 in [1], one can show that (3.21) admits a minimizer z † . Therefore, we conclude (3.20) by (3.24).
Remark 3.5. Theorem 3.4 implies that in order to find the solution u * of the PDE system (3.5), we need to solve (3.21). To handle the constraints of (3.21), we use the methods of eliminating variables or relaxation described in Subsection 3.3 of [1] when necessary.
Remark 3.6. To compute the inversion of Q(ψ, ψ) + γI, we use the techniques introduced in Subsection 2.2 and Remark 2.2.
Next, we estimate the error between the approximation solution u † given in (3.20) and the strong solution u * of (3.5).
Theorem 3.7. Let u * be the strong solution of (3.5). Given γ > 0, let (u † , z † ) be a solution of (3.10), where u † is given in (3.20), and z † solves (3.21). Define the sampling operator S x : U → 2 (x) by Denote S * x as the adjoint of S x . Then, S * x c = c T Kψ, ∀c ∈ 2 (x). Define Then, (3.28) Proof. Since S * x is the adjoint of S x , for each c ∈ 2 (x), we have Thus, we conclude (3.26). Meanwhile, for γ > 0, by the same arguments as the proof of Lemma 3.3, S x admits similar properties to those of S x , i.e., and We define the operator Given γ > 0, let u † be as in (3.20) and z † be a solution to (3.21). We first bound u † − Lz † in U.
Remark 3.8. The term λ x in (3.27) is a straightforward extension of Definition 1 in [27] and takes into account the information about values of linear operators evaluated at sample points. Similar to [27], we say that the set x consisting of samples provides rich data if λ x > 0 and it provides poor data if λ x = 0. Indeed, when we use the uniform sampling and the number of samples N is large enough, we will have λ x > 0. To see this, we suppose that S x v contains the values of v evaluated at the sample points for v ∈ U. Then, Thus, as N → +∞ and γ → 0, the first term of (3.28) tends to 0.
3.3. The Nyström Approximation Error. As discussed in Remark 3.9, the upper bound in (3.28) is small when Q(ψ, ψ) approximates K(ψ, ψ) well. In this subsection, we bound the spectral norm of K(ψ, ψ) − Q(ψ, ψ), denoted by K(ψ, ψ) − Q(ψ, ψ) . We recall that Q(ψ, ψ) = K(ψ, φ)(K(φ, φ)) −1 K(φ, ψ). Thus, we can view Q(ψ, ψ) as the Nyström approximation to K(ψ, ψ). Hence, we call K(ψ, ψ) − Q(ψ, ψ) the Nyström approximation error in the SGP method. A typical Nyström approximation approximates a symmetric positive semi-definite matrix G ∈ R n×n , n ∈ N, by sampling m columns from G and does not assume that columns of G have close relations with each other. The authors in [3] show that for any m uniformly sampled columns, with a high probability, the approximation error is O(n/ √ m). Later, [10] improves the bound from O(n/ √ m) to O(n/m 1−ρ ) for ρ ∈ (0, 1/2) under the assumption that the eigenvalues of G have a big eigengap. In our setting, since ψ has linear operators corresponding to the same sample points, the columns of K(ψ, ψ) are highly correlated. Thus, we approximate K(ψ, ψ) using a small number of inducing points instead of using columns from K(ψ, ψ). In general, instead of being taken from the samples, inducing points can be selected anywhere in the input space [4]. The error analysis has been thoroughly studied in [37], which indicates that the approximation error in the matrix Frobenious norm has an upper bound that is influenced by the quantization error defined in [37]. Since that quantization error is the objective function in the k-means clustering [6], k-means algorithm is typically used for the initialization of inducing inputs. In our settings, we assume that the inducing points are sampled from the set of samples. We adapt the framework in [10] to our settings and show that the Nyström approximation error has an upper bound that depends only on the numbers of samples and inducing points. In other words, though the dimension of K(ψ, ψ) is proportional to the product of the size of samples and the number of linear operators, the error bound is not influenced by the size of linear operators.
Proof. Let h ∈ H b and v ∈ H a . Then, there exists α and β such that v = α T Kφ and h = β T Kψ. We have Hence, we obtain which concludes (3.44).
Next, we estimate the upper bound for E(H a ). For simplification, we write ψ := (ψ Ω , ψ ∂Ω ), where ψ Ω represents the collection of linear operators for interior points and ψ ∂Ω consists of linear operators for boundary samples. We write x := (x 1 , x 2 ), where we recall that x contains all the sample points, x 1 is the set of interior points in Ω, and x 2 stands for the collection of boundary points. We define the operators L x1 and L x2 such that for any function v ∈ U, (3.45) where N Ω is the size of x 1 and N ∂Ω is the length of x 2 . The next proposition connects the eigenvalues and the eigenfunctions of L x1 and L x2 to the eigenvalues and the eigenvectors of K(ψ Ω , ψ Ω ) and K(ψ ∂Ω , ψ ∂Ω ). Let {λ j } R1 j=1 be the set of eigenvalues of K(ψ Ω , ψ Ω ), which are ranked in the descending order. Define the operator L x1 as in (3.45). Then, the eigenvalues of L x1 are {λ j /N Ω } R1 j=1 . Moreover, let ϕ 1 , . . . , ϕ R1 be the corresponding eigenfunctions of L x1 that have normalized functional norms, i.e. ϕ i , ϕ j = δ ij , 1 i, j R 1 . Let V = [v 1 , . . . , v R1 ] be the orthonormal eigenvector matrix of K(ψ Ω , ψ Ω ). Then, we have and where ψ Ω,j is the j th component of ψ Ω . Similarly, let {τ j } R2 j=1 be the set of eigenvalues of K(ψ ∂Ω , ψ ∂Ω ) sorted in the descending order and let L x2 be as in (3.45). Then, the eigenvalues of L x2 are {τ j /N ∂Ω } R2 j=1 . Meanwhile, let ζ 1 , . . . , ζ R2 be the eigenfunctions of L x2 that have normalized functional norms. Let W = [w 1 , . . . , w R2 ] be the orthonormal eigenvector matrix of K(ψ ∂Ω , ψ ∂Ω ). Then, we obtain and where ψ ∂Ω,j is the j th element of ψ ∂Ω .
Proof. We only give proofs related to L x1 since the arguments for L x2 are similar. Let λ j be j th eigenvalue of L x1 and let ϕ j be the corresponding eigenfunction with a normalized norm. Then, we have Acting ψ Ω on both sides of (3.54), we get (3.55) which implies that N Ω λ j is the j th eigenvalue of K(ψ Ω , ψ Ω ) and [ψ Ω , ϕ j ] is the corresponding eigenvector.
Hence, we confirm that if {λ j } R1 j=1 contains the eigenvalues of K(ψ Ω , ψ Ω ), {λ j /N Ω } R1 j=1 is the set of the eigenvalues of L x1 . Furthermore, let V = [v 1 , . . . , v R1 ] be the orthonormal eigenvector matrix of K(ψ Ω , ψ Ω ). Then, (3.55 which yields Since ϕ j is normalized, we obtain from (3.57) and (3.58) that (3.60) Thus, we get from (3.60) that where the last equality follows (3.46). Thus, we conclude (3.49). Using , For any r := (r 1 , which is the worst error in approximating any h ∈ H r b by functions in H a . The next proposition bounds E(H a ) by E(H a , r).
The above proposition implies that it is sufficient to bound E(H a , r) in order to estimate the Nyström approximation error. To proceed with the argument, we make the following assumptions about the inducing points and the kernel associated with the RKHS U.
Assumption 1. Assume that the inducing points are uniformly sampled from the collection of sample points. More precisely, let x := (x 1 , x 2 ), where x 1 is the set of interior sample points in Ω and x 2 is the collection of boundary samples, and let x := ( x 1 , x 2 ) be the set of inducing points. Assume that x 1 and x 2 are uniformly sampled from x 1 and x 2 , separately.
Assumption 2. Let ψ be as in (3.9). Denote by ψ x := (δ x • L 1 , . . . , δ x • L D ) the vector of linear operators in ψ that is associated with the sample point at x ∈ x. Assume that there exists a constant C > 0 such that Remark 3.13. We note that if the kernel K in (3.1) generated by K is smooth enough and the domain Ω is compact, Assumption 2 holds.
Let φ be as in (3.6). For clarification, we denote by φ := (φ Ω , φ ∂Ω ), where φ Ω represents the set of linear operators for interior points and φ ∂Ω includes boundary linear operators. Let x := ( x 1 , x 2 ) be as in Assumption 1, let M Ω be the size of x 1 , and let M ∂Ω be the length of x 2 . We define the operators L x1 and L x2 such that for any function v ∈ U, The following proposition provides upper bounds for the Hilbert-Schmidt norms of L x1 − L x1 and L x2 − L x2 . We recall that for any linear operator L : H → H, where H is a Hilbert space, we denote by L HS and L 2 the Hilbert-Schmidt norm and the spectral norm, respectively, i.e.
where {e i , i = 1, . . . } is a complete orthogonal basis of H. We note that for any linear operator L : H → H, L 2 L HS .
Proposition 3.14. Suppose that Assumptions 1 and 2 hold. Let L x1 , L x2 be as in (3.45) and let L x1 , L x2 be given in (3.65). Then, with a probability 1 − δ, 0 < δ < 1, there exists a constant C > 0 such that Proof. Here, we adapt the proof of Corollary 8 in [10]. Let ψ x be given in Assumption 2 for x ∈ x. For any v ∈ U, we define Then, , which are endowed with the norm of U. Under Assumption 1, we observe that L x1 and L x1 map H 1 to H 1 , and that L x2 and L x2 map H 2 to H 2 . Thus, by the definition of the Hilbert-Schmidt norm in (3.66), there exists a constant C such that where the last equality follows by (3.49) and we use Assumption 2 in the last inequality. Similarly, we get , with a probability at least 1 − δ, 0 < δ < 1, there exists a constant C such that Proof. The proof is an adaption of the arguments of Theorem 7 in [10].
Finally, we are ready to give a bound for E(H a ) in the following theorem.
Remark 3.17. Proposition 3.10 and Theorem 3.16 imply that K(ψ, ψ) − K(ψ, φ)(K(φ, φ)) −1 K(φ, ψ) admits an upper bound that depends only on the size of sample points and is not influenced by the number of linear operators in the PDE.

Numerical Results
In this section, we demonstrate the efficacy of the SGP method by solving a nonlinear elliptic equation in Besides, we randomly choose M points from the samples and treat them as inducing points. To show the performance of the SGP method, we first plot samples, inducing points, loss histories of the Gauss-Newton iterations, and contours of pointwise solution errors for fixed values of N and M . Though we solve an equation at a collection of uniform random samples, we compute the solution error by reevaluating the numerical approximation function at a 100 × 100 equally spaced grid using (3.20) and by comparing the results with the values of the true solution. Then, we record the L ∞ errors of the GP method in [1] and the SGP algorithm for different values of N and M . Results are averaged over 10 realizations of random sample points.
Our implementation is based on the code † of [1], which leverages Python with the JAX package for automatic differentiation. Our experiments are conducted only on CPUs. Additional performance speedups can be obtained by considering accelerated hardware such as Graphics Processing Units (GPU).

A Nonlinear Elliptic Equation.
In this example, we reconsider the nonlinear elliptic equation in Section 3.5 of [1] on a larger domain. More precisely, we consider (2.1) with d = 2, Ω = (0, 3) 2 , τ (u) = u 3 , and g(x) = 0. We prescribe the solution u to be sin(πx 1 ) sin(πx 2 ) + 4 sin(4πx 1 ) sin(4πx 2 ) for x := (x 1 , x 2 ) ∈ Ω and compute f accordingly. We use the Gaussian kernel K(x, y) = exp(− |x−y| 2 2σ 2 ) with the lengthscale parameter σ = 0.2. Given N and M , we set N Ω = 0.75 × N and M Ω = 0.75 × M . In this example, we choose γ = η = 10 −12 (see Subsection 2.2 and Remark 2.2). The algorithm stops once the error between two successive steps is less than 10 −5 . Figure 1 shows the numerical results for N = 4800 and M = 1200. We choose γ = η = 10 −12 . The uniform sample points used are plotted in Figure 1a, from which we randomly choose M including points shown in Figure 1b. We solve the unconstrained minimization problem in (2.5). The convergence history of the Gauss-Newton iteration in 1c shows that the SGP method converges in 5 iterations. In Figure 1d, we plot the contour of pointwise errors between the numerical approximation and the true solution on a 100 × 100 equally spaced grid.
To further explore the efficacy of our algorithm. We compare the L ∞ errors of the SGP method and the GP algorithm in [1] for different numbers of sample points and inducing points in Table 1. We see that when N = M , both algorithms achieve the accuracy of the same magnitude. Meanwhile, when N >> M , the SGP method with N sample points and M inducing points is more accurate than the GP method with M sample points. Furthermore, the SGP method with N sample points, half of which are used as inducing points, achieves comparable accuracy to the GP method with the same number of sample points. Hence, by using (2.9), the size of the matrix needed to be inverted in the SGP method can be significantly reduced. In other words, the SGP approach consumes much less computational time while retaining desirable accuracy compared to the GP method.

Burgers' equation.
In this example, we consider the same Burger's equation as in [1]. More precisely, given ν = 0.02, we seek to find u solving We uniformly sample N points in the space-time domain, from which M samples are randomly chosen as inducing points. All experiments in this example use a fixed ratio of interior points, N Ω /N = M Ω /M = 5/6. To take into consideration the space and time variability of the solution to Burger's equation, as in [1], we choose the anisotropic kernel with (σ 1 , σ 2 ) = (0.3, 0.05). For comparison, as in [1], the true solution is calculated from the Cole-Hopf transformation and the numerical quadrature. We use the technique of eliminating variables to deal with the constraints in (3.21) (see Subsection 3.3 in [1]) and use the Gauss-Newton method to solve the resulting unconstrained minimization problem. All the experiments in this example stop once the errors between two successive steps are less than 10 −5 . Figure 2 shows the numerical results for N = 2400 and M = 600 with parameters γ = η = 10 −6 . The samples and the inducing points are given in Figures 2a and 2b. The history of Gauss-Newton iterations is plotted in Figure 2c, which implies the fast convergence of the SGP method. The absolute values of the pointwise errors between the numerical solution and the true solution are shown in Figure 2d. Similar to the results of the GP method given in [1], the maximum errors occurred close to the shock. In Figures 2e-2g, we compare the numerical and true solutions at times t = 0.25, 0.5, 0.75 to highlight the accuracy of our SGP method.
In Table 2, we present the L ∞ errors of the SGP method and the GP algorithm in [1] for different choices of N and M . Similar to the nonlinear elliptic example, the SGP algorithm yields comparable accurate numerical solutions to the GP method by using less than half of the samples as inducing points.    Table 2. Burger's equation: L ∞ errors of the SGP method and the GP algorithm for different numbers of sample points and inducing points. The N/A means not available. We set γ = η = 10 −6 for N 1200 and choose γ = η = 10 −8 for N = 2400.
We prescribe the true solution u(t, x) = (sin(πx) + 2 cos(2πx))e −t and compute f , g, h 1 , and h 2 accordingly. We fix the ratios N Ω /N = M Ω /M = 6/7. Meanwhile, we use the anisotropic kernel in (4.1) with the same lengthscales, i.e., (σ 1 , σ 2 ) = (0.3, 0.05). We use the technique of eliminating variables to handle the constraints in (3.21) (see Subsection 3.3 of [1]). All the experiments in this example stop once the errors between two successive steps are less than 10 −5 or the iteration numbers exceed 20.  Figure 3c, which shows that the algorithm converges within 6 iterations. The absolute values of the pointwise errors between the numerical solution and the true solution are given in Figure 3d. Table 3 records the L ∞ errors of the SGP method and the GP algorithm in [1] for different values of N and M , which justifies the efficacy of the SGP method.  Table 3. The parabolic PDE: L ∞ errors of the SGP method and the GP algorithm for different numbers of sample points and inducing points. The N/A means not available. We set γ = η = 10 −10 .

Conclusion and Future Work
In this paper, we present a SGP framework to solve nonlinear PDEs. The SGP method finds a numerical solution to a PDE in the RKHS associated with a low-rank kernel generated by inducing points. Meanwhile, the approximation solution can be viewed as a maximum a posteriori probability estimator of a SGP conditioned a noisy observation of values of linear operators satisfying the PDE at a finite set of samples.
Our numerical experiments imply that the SGP algorithm consumes less computational time than the GP method in [1] without losing much accuracy when both methods use uniform samples. We notice that the positions of inducing points greatly influence the performance of our algorithm. Hence, a potential future work is to investigate a better way to put inducing points. Meanwhile, the choice of hyperparameters has a  profound impact on the performance of our method. The probabilistic interpretation of the SGP approach in Subsection 2.3 provides a way for hyperparameter learning in future work.